abstract.py

from statistics import NormalDist
from itertools import accumulate
from functools import cache
from operator import or_, and_
from math import comb, log2, floor, ceil, e
from typing import Iterator

from .shared import *


class StoreBHV:
    def __deepcopy__(self, memodict=None):
        return self.from_bytes(self.to_bytes())

    def to_bytes(self):
        raise NotImplementedError()

    @classmethod
    def from_bytes(cls, bs: bytes):
        raise NotImplementedError()

    def bits(self):
        for b in self.to_bytes():
            yield b & 0x1
            yield (b >> 1) & 0x1
            yield (b >> 2) & 0x1
            yield (b >> 3) & 0x1
            yield (b >> 4) & 0x1
            yield (b >> 5) & 0x1
            yield (b >> 6) & 0x1
            yield (b >> 7) & 0x1

    @classmethod
    def from_bitstream(cls, bits_s: 'Iterator[bool]') -> Self:
        it = iter(bits_s)
        bs = bytes(
            next(it) | (next(it) << 1) | (next(it) << 2) | (next(it) << 3) | (next(it) << 4) | (next(it) << 5) | (next(it) << 6) | next(it)  << 7
            for _ in range(DIMENSION//8)
        )
        return cls.from_bytes(bs)

    def bitstring(self) -> str:
        return "".join("1" if b else "0" for b in self.bits())

    @classmethod
    def from_bitstring(cls, s: str) -> Self:
        return cls.from_bitstream(c == "1" for c in s)


class BaseBHV(StoreBHV):
    @classmethod
    def rand(cls) -> Self:
        raise NotImplementedError()

    @classmethod
    def nrand(cls, n: int) -> list[Self]:
        return [cls.rand() for _ in range(n)]

    @classmethod
    def majority(cls, vs: list[Self]) -> Self:
        raise NotImplementedError()

    @classmethod
    def majority3(cls, x, y, z) -> Self:
        return cls.majority([x, y, z])

    def __eq__(self, other: Self) -> bool:
        raise NotImplementedError()

    def __xor__(self, other: Self) -> Self:
        raise NotImplementedError()

    @classmethod
    def parity(cls, vs: list[Self]) -> Self:
        a = cls.ZERO
        for v in vs:
            a = a ^ v
        return a

    def active(self) -> int:
        raise NotImplementedError()

    def permute(self, permutation_id: int) -> Self:
        raise NotImplementedError()

    def active_fraction(self) -> float:
        return self.active()/DIMENSION

    def hamming(self, other: Self) -> int:
        return (self ^ other).active()

    def closest(self, vs: list[Self]) -> int:
        return min(range(len(vs)), key=lambda i: self.hamming(vs[i]))

    def top(self, vs: list[Self], k: int) -> list[int]:
        return list(sorted(range(len(vs)), key=lambda i: self.hamming(vs[i])))[:k]

    def within(self, vs: list[Self], d: int) -> list[int]:
        return list(filter(lambda i: self.hamming(vs[i]) <= d, range(len(vs))))

    def within_std(self, vs: list[Self], d: float, relative: bool = False) -> list[int]:
        return self.within(vs, int(DIMENSION*(self.std_to_frac(d + self.EXPECTED_RAND_APART*relative))))

    def distribution(self, vs: list[Self], metric=lambda x, y: x.std_apart(y), softmax: bool = False, base: int = e):
        ds = [metric(self, v) for v in vs]
        if softmax:
            eds = [base**d for d in ds]
            sed = sum(eds)
            return [ed/sed for ed in eds]
        return ds

    def bit_error_rate(self, other: Self) -> float:
        return self.hamming(other)/DIMENSION

    EXPECTED_RAND_APART: float = NormalDist(0, (DIMENSION/4)**.5).zscore(DIMENSION/2)
    SPACE_WIDTH: float = NormalDist(0, (DIMENSION/4)**.5).zscore(DIMENSION)
    def std_apart(self, other: Self, relative=False) -> float:
        return self.frac_to_std(self.bit_error_rate(other) - .5*relative, relative)

    def std_relation(self, other: Self, compact=False, stdev=6, rounding=1) -> str:
        #    d=0   d < stdev   d < EXPECTED - stdev   EXPECTED - stdev < d    d = 0    d < EXPECTED + stdev  EXPECTED + stdev < d   WIDTH - stdev < d   d=WIDTH
        #   equal    noisy     closer than average       about average       halfway       about avergage       further opposite      noisy opposite    opposite
        #    NP       dNP                dN                   dNE               E               dSE                    dS                  dSP            SP
        r = int if rounding == 0 else lambda x: round(x, rounding)
        d = self.std_apart(other)
        dr = d - self.EXPECTED_RAND_APART
        if d == 0.:
            return "NP" if compact else "equal"
        elif d == self.SPACE_WIDTH:
            return "SP" if compact else "opposite"
        elif d < stdev:
            return f"{r(d)}NP" if compact else f"{r(d)}σ away from equal"
        elif self.SPACE_WIDTH - d < stdev:
            return f"{r(self.SPACE_WIDTH - d)}SP" if compact else f"{r(self.SPACE_WIDTH - d)}σ away from opposite"
        elif dr < 0. and abs(dr) < stdev:
            return f"{r(abs(dr))}NE" if compact else f"{r(abs(dr))}σ random closer to equal"
        elif dr < 0.:
            return f"{r(abs(dr))}N" if compact else f"{r(abs(dr))}σ closer to equal"
        elif dr > 0. and abs(dr) < stdev:
            return f"{r(dr)}SE" if compact else f"{r(dr)}σ random closer to opposite"
        elif dr > 0.:
            return f"{r(dr)}S" if compact else f"{r(dr)}σ closer to opposite"
        else:
            return "halfway"

    @staticmethod
    def normal(mean=0., af=.5):
        return NormalDist(mean, (DIMENSION*af*(1 - af))**.5)

    @staticmethod
    def maj_frac(n):
        return comb(n - 1, (n - 1)//2)/2**n

    @staticmethod
    def frac_to_std(frac, relative=False):
        n = AbstractBHV.normal(.5*DIMENSION if relative else 0)
        stdvs = n.zscore((frac*DIMENSION + .5*DIMENSION) if relative else frac*DIMENSION)
        return stdvs

    @staticmethod
    def std_to_frac(std, relative=False):
        n = AbstractBHV.normal(.5*DIMENSION if relative else 0)
        frac = (std*n.stdev + n.mean)/DIMENSION
        return frac - .5 if relative else frac

    def zscore(self) -> float:
        n = AbstractBHV.normal(0.5*DIMENSION)
        return n.zscore(self.active())

    def close(self, other: Self, stdvs=6) -> bool:
        return abs(self.std_apart(other)) <= stdvs

    def related(self, other: Self, stdvs=6) -> bool:
        return abs(self.std_apart(other, relative=True)) > stdvs

    @classmethod
    def _majority5_via_3(cls, a: Self, b: Self, c: Self, d: Self, e: Self) -> Self:
        return cls.majority3(a, cls.majority3(b, c, d), cls.majority3(e, d, cls.majority3(c, b, a)))

    @classmethod
    def _majority7_via_3(cls, a: Self, b: Self, c: Self, d: Self, e: Self, f: Self, g: Self) -> Self:
        mdef = cls.majority3(d, e, f)
        mabc = cls.majority3(a, b, c)
        return cls.majority3(g, cls.majority3(c, mdef, cls.majority3(a, b, mdef)), cls.majority3(f, mabc, cls.majority3(d, e, mabc)))

    @classmethod
    def _majority9_via_3(cls, a: Self, b: Self, c: Self, d: Self, e: Self, f: Self, g: Self, h: Self, i: Self) -> Self:
        mabc = cls.majority3(a, b, c)
        mdef = cls.majority3(d, e, f)
        mghi = cls.majority3(g, h, i)

        l_ = cls.majority3(cls.majority3(c, mdef, b), mdef, a)
        l = cls.majority3(cls.majority3(i, l_, h), l_, g)

        r_ = cls.majority3(cls.majority3(c, mghi, b), mghi, a)
        r = cls.majority3(cls.majority3(f, r_, e), r_, d)

        m = cls.majority3(mabc, mdef, mghi)
        return cls.majority3(l, m, r)

    ZERO: Self


class BooleanAlgBHV(BaseBHV):
    def __and__(self, other: Self) -> Self:
        raise NotImplementedError()

    def __or__(self, other: Self) -> Self:
        raise NotImplementedError()

    def __invert__(self) -> Self:
        return self ^ self.ONE

    @classmethod
    def rand2(cls, power: int) -> Self:
        assert power >= 0
        r = cls.rand()
        return r if power == 0 else r & cls.rand2(power - 1)

    @classmethod
    def nrand2(cls, n: int, power: int) -> list[Self]:
        assert power >= 0
        return [cls.rand2(power) for _ in range(n)]

    def select(self, when1: Self, when0: Self) -> Self:
        return when0 ^ (self & (when0 ^ when1))

    def select_rand(self, other: Self) -> Self:
        return self.rand().select(self, other)

    def select_rand2(self, other: Self, power_left: int) -> Self:
        return self.rand2(power_left).select(self, other)

    def randomize_half(self) -> Self:
        return self.rand().select_rand(self)

    def randomize_pow(self, pow_randomize: int) -> Self:
        return self.rand().select_rand2(self, pow_randomize)

    def flip_half(self) -> Self:
        return self ^ self.rand()

    def flip_pow(self, pow_flip: int) -> Self:
        return self ^ self.rand2(pow_flip)

    def flip_pow_on(self, pow_flip_on: int) -> Self:
        return self | ~self.rand2(pow_flip_on)

    def flip_pow_off(self, pow_flip_off: int) -> Self:
        return self & self.rand2(pow_flip_off)

    def overlap(self, other: Self, distance=False) -> float:
        active = (self & other).active_fraction()
        return 1. - active if distance else active

    def jaccard(self, other: Self, distance=False) -> float:
        union_active = (self | other).active()
        if union_active == 0:
            raise RuntimeError("Jaccard index where both vectors are zero")
        res = (self & other).active() / union_active
        return 1. - res if distance else res

    def cosine(self, other: Self, distance=False) -> float:
        s_active = self.active()
        o_active = other.active()
        if not s_active or not o_active:
            raise RuntimeError("Cosine similarity where one of the two vectors is zero")
        res = (self & other).active() / (s_active*o_active)**.5
        return 1. - res if distance else res

    def mutual_information(self, other: Self, distance=False) -> float:
        nself = ~self
        nother = ~other

        # Probabilities
        P00 = nself.overlap(nother)
        P01 = nself.overlap(other)
        P10 = self.overlap(nother)
        P11 = self.overlap(other)

        # Marginal probabilities
        PX0 = (P00 + P01)
        PX1 = (P10 + P11)
        PY0 = (P00 + P10)
        PY1 = (P01 + P11)

        # Mutual Information
        MI = 0
        if P00 and PX0 and PY0:
            MI += P00 * log2(P00 / (PX0 * PY0))
        if P01 and PX0 and PY1:
            MI += P01 * log2(P01 / (PX0 * PY1))
        if P10 and PX1 and PY0:
            MI += P10 * log2(P10 / (PX1 * PY0))
        if P11 and PX1 and PY1:
            MI += P11 * log2(P11 / (PX1 * PY1))

        return 1 - MI if distance else MI

    def tversky(self, other: Self, l: float, r: float) -> float:
        o = self.overlap(other)
        lr = self.overlap(~other)
        rl = other.overlap(~self)
        return o/(l*lr + r*rl + o)

    @classmethod
    def majority(cls, vs: list[Self]) -> Self:
        n = len(vs)

        if n == 0:
            return cls.rand()
        elif n == 1:
            return vs[0]
        elif n % 2 == 0:
            return cls.majority([cls.rand()] + vs)
        else:
            return cls._logic_majority(vs)

    @classmethod
    def representative(cls, vs: list[Self]) -> Self:
        n = len(vs)

        if n == 0:
            return cls.rand()
        elif n == 1:
            return vs[0]
        elif n == 2:
            return cls.select(cls.rand(), vs[0], vs[1])
        else:
            return cls._logic_representative(vs)

    @classmethod
    def _logic_representative(cls, vs: list[Self]) -> Self:
        n = len(vs)

        def smallest_power2_larger(x):
            # returns log2(x) on x power of 2, else its adds the smallest d such that log2(x + d) is integer
            k = 0
            while 2 ** k < x:
                k += 1
            return k

        circumscribed = smallest_power2_larger(n)
        inscribed = circumscribed - 1

        rs = cls.nrand(circumscribed)

        def interval(p):
            lower = (int(p, 2) if p else 0) * 2 ** (circumscribed - len(p))
            upper = min(n, lower + 2 ** (circumscribed - len(p)))
            return lower, upper

        def length(lu):
            d = lu[1] - lu[0]
            if d >= 1:
                return d
            else:
                return 0

        def rec(p):
            L0, L1 = length(interval(p + '0')), length(interval(p + '1'))

            if L0 and L1:
                # TODO could be replaced by fraction (created with AND, OR, and RAND)
                d = rs[len(p)] if L0 == L1 else cls.random(L1 / (L0 + L1))
                return d.select(rec(p + '1'), rec(p + '0'))
            elif L0:
                return rec(p + '0')
            elif L1:
                return rec(p + '1')
            else:
                i = int(p, 2)
                return vs[i]

        return rec('')

    @classmethod
    def _logic_majority(cls, vs: list[Self]) -> Self:
        assert len(vs) % 2 == 1,  "The true majority function which is only defined on uneven length inputs, see `majority`"
        half = len(vs)//2
        disjunction = list(accumulate(vs[:half-1:-1], or_))[1:]
        conjunction = list(accumulate(vs[:half-1:-1], and_))[1:]

        @cache
        def rec(ons, offs):
            if ons + 1 > half:
                return disjunction[-offs - 1]
            if offs + 1 > half:
                return conjunction[-ons - 1]
            return vs[ons + offs].select(
                rec(ons + 1, offs),
                rec(ons, offs + 1)
            )

        return rec(0, 0)

    # Alternative implementations

    @classmethod
    def _majority3_select(cls, a: Self, b: Self, c: Self) -> Self:
        # C:  1 0 0 1 0 1 1

        # |:  1 1 1 1 0 1 0
        # &:  1 0 0 0 0 0 0

        # M:  1 0 0 1 0 1 0

        return a.select(b | c, b & c)

    @classmethod
    def _majority5_select(cls, x4: Self, x3: Self, x2: Self, x1: Self, x0: Self) -> Self:
        t1 = x1 | x0
        b1 = x1 & x0
        m2 = x2.select(t1, b1)
        t2 = x2 | t1
        b2 = x2 & b1
        m3h = x3.select(t2, m2)
        m3l = x3.select(m2, b2)
        m4 = x4.select(m3h, m3l)
        return m4

    @classmethod
    def _majority7_select(cls, x6: Self, x5: Self, x4: Self, x3: Self, x2: Self, x1: Self, x0: Self) -> Self:
        t1 = x1 | x0
        b1 = x1 & x0
        m2 = x2.select(t1, b1)
        t2 = x2 | t1
        b2 = x2 & b1
        m3h = x3.select(t2, m2)
        m3l = x3.select(m2, b2)
        m4 = x4.select(m3h, m3l)
        t3 = x3 | t2
        b3 = x3 & b2
        m4h = x4.select(t3, m3h)
        m4l = x4.select(m3l, b3)
        m5h = x5.select(m4h, m4)
        m5l = x5.select(m4, m4l)
        m6 = x6.select(m5h, m5l)
        return m6

    @classmethod
    def _majority9_select(cls, x8: Self, x7: Self, x6: Self, x5: Self, x4: Self, x3: Self, x2: Self, x1: Self, x0: Self) -> Self:
        t1 = x0 | x1
        b1 = x0 & x1
        m2 = x2.select(t1, b1)
        t2 = x2 | t1
        b2 = x2 & b1
        m3h = x3.select(t2, m2)
        m3l = x3.select(m2, b2)
        t3 = x3 | t2
        b3 = x3 & b2
        m4h = x4.select(t3, m3h)
        m4 = x4.select(m3h, m3l)
        m4l = x4.select(m3l, b3)
        t4 = x4 | t3
        b4 = x4 & b3
        m5hh = x5.select(t4, m4h)
        m5h = x5.select(m4h, m4)
        m5l = x5.select(m4, m4l)
        m5ll = x5.select(m4l, b4)
        m6h = x6.select(m5hh, m5h)
        m6 = x6.select(m5h, m5l)
        m6l = x6.select(m5l, m5ll)
        m7h = x7.select(m6h, m6)
        m7l = x7.select(m6, m6l)
        m8 = x8.select(m7h, m7l)
        return m8

    @classmethod
    def _majority3(cls, a: Self, b: Self, c: Self) -> Self:
        # a:  1 0 0 1 0 1 1
        # b:  1 1 0 1 0 1 0
        # c:  1 0 1 0 0 0 0
        # M:  1 0 0 1 0 1 0

        # at least 2/3 agreeing on TRUE
        abh = a & b
        bch = b & c
        cah = c & a
        h = abh | bch | cah
        return h

    @classmethod
    def _majority5(cls, a: Self, b: Self, c: Self, d: Self, e: Self) -> Self:
        # at least 3/5 agreeing on TRUE

        # 2*10 AND
        # 9*1 OR
        # (b ∧ d ∧ e) ∨ (a ∧ d ∧ e) ∨ (b ∧ c ∧ e) ∨ (a ∧ c ∧ e) ∨ (b ∧ c ∧ d) ∨ (a ∧ c ∧ d) ∨ (c ∧ d ∧ e) ∨ (a ∧ b ∧ e) ∨ (a ∧ b ∧ d) ∨ (a ∧ b ∧ c)

        ab = a & b
        cd = c & d
        de = d & e
        ce = c & e

        bde = b & de
        ade = a & de

        bce = b & ce
        ace = a & ce

        bcd = b & cd
        acd = a & cd
        cde = e & cd

        abe = ab & e
        abd = ab & d
        abc = ab & c

        h = bde | ade | bce | ace | bcd | acd | cde | abe | abd | abc
        return h

    @classmethod
    def _majority_via_custom(cls, vs: list[Self]):
        n = len(vs)
        if n == 0:
            return cls.rand()
        elif n == 1:
            return vs[0]
        elif n == 2:
            return vs[0].select_rand(vs[1])
        elif n == 3:
            return cls._majority3(vs[0], vs[1], vs[2])
        elif n == 4:
            return cls._majority5_via_3(vs[0], vs[1], vs[2], vs[3], cls.rand())
        elif n == 5:
            return cls._majority5_via_3(vs[0], vs[1], vs[2], vs[3], vs[4])
        elif n == 6:
            return cls._majority7_via_3(vs[0], vs[1], vs[2], vs[3], vs[4], vs[5], cls.rand())
        elif n == 7:
            return cls._majority7_via_3(vs[0], vs[1], vs[2], vs[3], vs[4], vs[5], vs[6])
        elif n == 8:
            return cls._majority9_via_3(vs[0], vs[1], vs[2], vs[3], vs[4], vs[5], vs[6], vs[7], cls.rand())
        elif n == 9:
            return cls._majority9_via_3(vs[0], vs[1], vs[2], vs[3], vs[4], vs[5], vs[6], vs[7], vs[8])
        else:
            raise RuntimeError(f"No optimal implementations beyond MAJORITY-9 (tried to take MAJORITY-{n})")

    ONE: Self


class FractionalBHV(BooleanAlgBHV):
    @classmethod
    def random(cls, active: float) -> Self:
        raise NotImplementedError()

    @classmethod
    def nrandom(cls, n: int, active: float) -> list[Self]:
        assert 0. <= active <= 1.
        return [cls.random(active) for _ in range(n)]

    def select_random(self, other: Self, frac_left: float) -> Self:
        return self.random(frac_left).select(self, other)

    def randomize_frac(self, frac_randomize: float) -> Self:
        return self.rand().select_random(self, frac_randomize)

    def flip_frac(self, frac_flip: float) -> Self:
        return self ^ self.random(frac_flip)

    def flip_frac_on(self, frac_flip_on: float) -> Self:
        return self | self.random(frac_flip_on)

    def flip_frac_off(self, frac_flip_off: float) -> Self:
        return self & self.random(1. - frac_flip_off)


class CryptoBHV(FractionalBHV):
    def swap_even_odd(self) -> Self:
        raise NotImplementedError()

    def rehash(self) -> Self:
        raise NotImplementedError()

    _FEISTAL_ROUNDS = 8
    _FEISTAL_SUBKEYS: list[Self]

    @classmethod
    def _feistal_round(cls, block: Self, round_key: Self) -> Self:
        L = block & cls.EVEN
        R = (block ^ (L ^ round_key).rehash().swap_even_odd()) & cls.ODD

        return (L | R).swap_even_odd()

    def feistal(self, k: Self, inv=False) -> Self:
        block = self
        rounds = range(self._FEISTAL_ROUNDS)

        for r in reversed(rounds) if inv else rounds:
            block = self._feistal_round(block, self._FEISTAL_SUBKEYS[r] & k)

        return block.swap_even_odd()

    EVEN: Self
    ODD: Self


class LevelBHV(FractionalBHV):
    @classmethod
    def level(cls, active_fraction: float) -> Self:
        raise NotImplementedError()

    @classmethod
    def threshold(cls, vs: list[Self], t: int) -> Self:
        # threshold(vs, t) = t < counts(vs)
        return cls._logic_threshold(vs, t)

    @classmethod
    def _logic_threshold(cls, vs: list[Self], t: int) -> Self:
        if t >= len(vs): return cls.ZERO
        if t < 0: return cls.ONE

        size = len(vs)
        t_ = size - t - 1

        grid: list[list[Self]] = [[None]*(t_+1) for _ in range(t + 1)]

        grid[t][t_] = vs[size - 1]

        for i in range(t):
            grid[t - i - 1][t_] = grid[t - i][t_] & vs[size - i - 2]

        for i in range(t_):
            grid[t][t_ - i - 1] = grid[t][t_ - i] | vs[size - i - 2]

        for i in range(t - 1, -1, -1):
            for j in range(t_ - 1, -1, -1):
                grid[i][j] = vs[i + j].select(grid[i + 1][j], grid[i][j + 1])

        return grid[0][0]

    @classmethod
    def window(cls, vs: list[Self], b: int, t: int) -> Self:
        # assert 0 <= b <= t <= len(vs)
        # window(vs, b, t) = threshold(vs, b) and not threshold(not vs, t)
        # window(vs, b, t) = b <= counts(vs) <= t
        #                  = b <= counts(vs) and counts(vs) <= t
        #                  = b <= counts(vs) and not (t < counts(vs))
        #                  = (b - 1 < counts(vs)) and not (t < counts(vs))
        return cls.threshold(vs, b - 1) & ~cls.threshold(vs, t)

    @classmethod
    def _ite_window(cls, vs: list[Self], b: int, t: int) -> Self:
        #   b   t
        #   |   |
        # --+---+----
        #   +---+      Sure that you'll land in the inclusive interval
        #              ons >= b and offs >= len(vs) - t
        # --     ----  Sure that you'll land in one of the two exclusive intervals
        #              ons > t or offs > len(vs) - b
        l = len(vs)

        @cache
        def rec(ons, offs):
            on1 = ons + 1 >= b and offs >= l - t
            on0 = ons + 1 > t or offs > l - b
            off1 = ons >= b and offs + 1 >= l - t
            off0 = ons > t or offs + 1 > l - b

            if on1 and off0:
                return vs[ons + offs]
            elif on0 and off1:
                return ~vs[ons + offs]
            elif on1:
                return vs[ons + offs] | rec(ons, offs + 1)
            elif on0:
                return ~vs[ons + offs] & rec(ons, offs + 1)
            elif off1:
                return ~vs[ons + offs] | rec(ons + 1, offs)
            elif off0:
                return vs[ons + offs] & rec(ons + 1, offs)
            else:
                return vs[ons + offs].select(
                    rec(ons + 1, offs),
                    rec(ons, offs + 1)
                )

        return rec(0, 0)

    @classmethod
    def agreement(cls, vs: list[Self], p: float) -> Self:
        assert .5 <= p <= 1.
        # agreement(vs, 100%) = all(v_i) or not any(v_i) = counts(vs) >= len(vs) or counts(vs) <= 0
        # agreement(vs, 3/5) = 5*counts(vs) >= 3*len(vs) or 5*counts(vs) <= 2*len(vs)
        # agreement(vs, 1/2) = true
        # agreement(vs, p) = not window(vs, p*len(vs), (1 - p)*len(vs))
        if p == .5: return cls.ONE
        t = ceil(p*len(vs))
        b = len(vs) - t
        return ~cls.window(vs, b + 1, t - 1)

    @staticmethod
    def best_division_window_bounds(n: int, f: float = .5, af=.5):
        # Note, uses Normal approximation and may deviate from the actual best p
        d = NormalDist(n*af, (n*af*(1. - af))**.5)
        b = round(d.inv_cdf((1. - f)/2))
        t = round(d.inv_cdf((1. + f)/2))
        return b, t

    @classmethod
    def variant(cls, vs: list[Self], f: float = .5, af=.5):
        # best_division(vs, f) = min_p E[(f - avg(agreement(vs, p)))]
        b, t = cls.best_division_window_bounds(len(vs), f, af)
        return cls.window(vs, b + 1, t - 1)


class AbstractBHV(CryptoBHV, LevelBHV):
    pass


class Permutation:
    @classmethod
    def nrandom(cls, n) -> list[Self]:
        return [cls.random() for _ in range(n)]

    @classmethod
    def random(cls) -> Self:
        raise NotImplementedError()

    def __mul__(self, other: Self) -> Self:
        raise NotImplementedError()

    def __invert__(self) -> Self:
        raise NotImplementedError()

    def __pow__(self, power):
        permutation = self
        for _ in range(power - 1):
            permutation = permutation*permutation
        return permutation

    def __call__(self, hv: 'AbstractBHV') -> 'AbstractBHV':
        raise NotImplementedError()

    IDENTITY: Self


class MemoizedPermutation(Permutation):
    _permutations: 'dict[int | tuple[int, ...], Self]'

    @classmethod
    def _get_singular(cls, permutation_id: int) -> Self:
        if permutation_id == 0:
            return cls.IDENTITY
        elif permutation_id in cls._permutations:
            return cls._permutations[permutation_id]
        elif -permutation_id in cls._permutations:
            inv_permutation = cls._permutations[-permutation_id]
            permutation = ~inv_permutation
            cls._permutations[permutation_id] = permutation
            return permutation
        else:
            permutation = cls.random()
            cls._permutations[permutation_id] = permutation
            return permutation

    @classmethod
    def _get_composite(cls, permutation_id: 'tuple[int, ...]') -> Self:
        # this can be optimized a lot by looking for partial compositions
        composite_permutation = cls._get_singular(permutation_id[0])
        for component_permutation_id in permutation_id[1:]:
            component_permutation = cls.get(component_permutation_id)
            composite_permutation = composite_permutation * component_permutation

        cls._permutations[permutation_id] = composite_permutation
        return composite_permutation

    @classmethod
    def get(cls, permutation_id: 'int | tuple[int, ...]') -> Self:
        permutation = cls._get_singular(permutation_id) if isinstance(permutation_id, int) else \
            cls._get_composite(permutation_id)
        return permutation