time_to_fc_v2

i_have_no_life.txt

Manual work debugging the damn thing /(>_<)\


'''
@ 0.5 prob & max_misses = 0
hit 1: 100 ms     Time to fc: E[  0 + 100, 0.5] = 200       actual: 200 ms
hit 2: 100 ms     Time to fc: E[200 + 100, 0.5] = 600       actual: 600 ms
hit 3: 100 ms     Time to fc: E[600 + 100, 0.5] = 1400      actual: 1400 ms


@ max_misses = 1
prob of pass note 1 = not (P(miss note 1) and P(miss == 1))
    P(hit note 1) = 0.5
    P(miss note 1) = 1 - P(hit note 1) = 0.5
    P(misses == 1) = nCr(note, max_misses)/sum(m: 0 -> max_misses, nCr(note, m)) = nCr(0, 1)/(nCr(0, 0) + nCr(0, 1)) = 0.0
        Total:  [] -> 0
            0 miss = [] -> 0
            1 miss = [] -> 0

    not (P(miss note 1) and P(miss == 1)) = 1 - (0.5*0.0) = 1.0 - 0.0 = 1.0

    E[100, P(pass note 1)] = E[100, 1.0]

prob of pass note 2 = not (P(hit note 2) and P(miss == 1))
    P(hit note 2) = 0.5
    P(miss note 2) = 1 - P(hit note 2) = 0.5
    P(misses == 1) = nCr(note, max_misses)/sum(m: 0 -> max_misses, nCr(note, m)) = nCr(1, 1)/(nCr(1, 0) + nCr(1, 1)) = 0.5
        Total:  [0], [1] -> 2
            0 miss = [0] -> 1/2
            1 miss = [1] -> 1/2

    not (P(miss note 1) and P(miss == 1)) = 1 - (0.5*0.5) = 1 - 0.25 = 0.75

    E[100, P(pass note 2)] = E[100, 0.75]

prob of pass note 3 = not (P(hit note 3) and P(miss == 1))
    P(hit note 3) = 0.5
    P(miss note 3) = 1 - P(hit note 3) = 0.5
    P(misses == 1) = nCr(note, max_misses)/sum(m: 0 -> max_misses, nCr(note, m)) = nCr(2, 1)/(nCr(2, 0) + nCr(2, 1)) = 0.6666
        Total:  [0, 0], [1, 0], [0, 1]  -> 3
        0 miss = [0, 0] -> 1/3
        1 miss = [1, 0], [0, 1] -> 2/3
        2 miss = [] -> 0/3

    not (P(miss note 3) and P(miss == 1)) = 1 - (0.5*0.6666) = 1 - 0.3333 = 0.6666

    E[100, P(pass note 3)] = E[100, 0.666‬]

@ 0.5 prob & max_misses = 1
hit 1: 100 ms     Time to fc: E[100, 1.0]             = 100           actual: 100    ms
hit 2: 100 ms     Time to fc: E[100 + 100, 0.75]      = 266.666       actual: 266.66 ms
hit 2: 100 ms     Time to fc: E[266.666 + 100, 0.666] = 550           actual: 550    ms

total FC time:
    E[ E[ E[ 100, 1.0 ] + 100, 0.75 ] + 100, 0.666‬ ]


@ max_misses = 2
prob of pass note 1 = not (P(miss note 1) and P(miss == 2))
    P(hit note 1) = 0.5
    P(miss note 1) = 1 - P(hit note 1) = 0.5
    P(misses == 2) = nCr(note, max_misses)/sum(m: 0 -> max_misses, nCr(note, m)) = nCr(0, 2)/(nCr(0, 0) + nCr(0, 1) + nCr(0, 2)) = 0.0
        Total:  [] -> 0
            0 miss = [] -> 0
            1 miss = [] -> 0

    not (P(miss note 1) and P(miss == 2)) = 1 - (0.5*0.0) = 1.0 - 0.0 = 1.0

    E[100, P(pass note 1)] = E[100, 1.0]

prob of pass note 2 = not (P(hit note 2) and P(miss == 2))
    P(hit note 2) = 0.5
    P(miss note 2) = 1 - P(hit note 2) = 0.5
    P(misses == 2) = nCr(note, max_misses)/sum(m: 0 -> max_misses, nCr(note, m)) = nCr(1, 2)/(nCr(1, 0) + nCr(1, 1) + nCr(1, 2)) = 0.0
        Total:  [0], [1] -> 2
            0 miss = [0] -> 1/2
              0.5
            1 miss = [1] -> 1/2
              0.5
            2 miss = []  -> 0/2
              0.0

        0.0/(0.5 + 0.5 + 0.0) = 0.0

    not (P(miss note 1) and P(miss == 1)) = 1 - (0.5*0.0) = 1.0 - 0.0 = 1.0

    E[100, P(pass note 2)] = E[100, 1.0]

prob of pass note 3 = not (P(hit note 3) and P(miss == 2))
    P(hit note 3) = 0.5
    P(miss note 3) = 1 - P(hit note 3) = 0.5
    P(misses == 2) = nCr(note, max_misses)/sum(m: 0 -> max_misses, nCr(note, m)) = nCr(2, 2)/(nCr(2, 0) + nCr(2, 1) + nCr(2, 2)) = 0.25
        0 miss = [0, 0]
            0.5*0.5 = 0.25
        1 miss = [1, 0], [0, 1]
            0.5*0.5 = 0.25
            0.5*0.5 = 0.25
        2 miss = [1, 1]
            0.5*0.5 = 0.25

        0.25/(0.25 + 0.25 + 0.25 + 0.25) = 0.25/1.0 = 0.25

    not (P(miss note 3) and P(miss == 2)) = 1 - (0.5*0.25) = 1 - 0.125 = 0.875

    E[100, P(pass note 3)] = E[100, 0.875‬]

prob of pass note 4 = not (P(hit note 4) and P(miss == 2))
    P(hit note 4) = 0.5
    P(miss note 4) = 1 - P(hit note 4) = 0.5
    P(misses == 2) = nCr(note, max_misses)/sum(m: 0 -> max_misses, nCr(note, m)) = nCr(3, 2)/(nCr(3, 0) + nCr(3, 1) + nCr(3, 2)) = 0.42857142857142855
        Total:  [0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1], [1, 1, 0], [1, 0, 1], [0, 1, 1] -> 7
        0 miss = [0, 0, 0] -> 1/7
            0.5*0.5*0.5 = 0.125
        1 miss = [1, 0, 0], [0, 1, 0], [0, 0, 1] -> 3/7
            0.5*0.5*0.5 = 0.125
            0.5*0.5*0.5 = 0.125
            0.5*0.5*0.5 = 0.125
        2 miss = [1, 1, 0], [1, 0, 1], [0, 1, 1] -> 3/7
            0.5*0.5*0.5 = 0.125
            0.5*0.5*0.5 = 0.125 => 0.375
            0.5*0.5*0.5 = 0.125

        0.375/(0.125 + 0.125 + 0.125 + 0.125 + 0.125 + 0.125 + 0.125) = 0.375/0.875 = 0.429

    not (P(miss note 3) and P(miss == 2)) = 1 - (0.5*0.429) = 1 - 0.214 = 0.785

    E[100, P(pass note 3)] = E[100, 0.785‬]

@ 0.5 prob & max_misses = 2
hit 1: 100 ms     Time to fc: E[  0 + 100, 1.0]       = 100       actual: 100 ms
hit 2: 100 ms     Time to fc: E[100 + 100, 1.0]       = 200       actual: 200 ms
hit 3: 100 ms     Time to fc: E[200 + 100, 0.875‬]     = 342.857   actual: 342.857 ms
hit 4: 100 ms     Time to fc: E[342.857 + 100, 0.785‬] = 563.636   actual: 563.636 ms

total FC time:
    E[ E[ E[ E[ 100, 1.0 ] + 100, 1.0 ] + 100, 0.875‬ ] + 100, 0.785 ]



@ max_misses = 2
prob of pass note 1 = not (P(miss note 1) and P(miss == 2))
    P(hit note 1) = 0.0
    P(miss note 1) = 1 - P(hit note 1) = 1.0
    P(misses == 2) = nCr(note, max_misses)/sum(m: 0 -> max_misses, nCr(note, m)) = nCr(0, 2)/(nCr(0, 0) + nCr(0, 1) + nCr(0, 2)) = 0.0
        Total:  [] -> 0
            0 miss = [] -> 0
            1 miss = [] -> 0

        0.0

    not (P(miss note 1) and P(miss == 2)) = 1 - (1.0*0.0) = 1 - 0.0 = 1.0

    E[100, P(pass note 1)] = E[100, 1.0]

prob of pass note 2 = not (P(hit note 2) and P(miss == 2))
    P(hit note 2) = 0.0
    P(miss note 2) = 1 - P(hit note 2) = 1.0
    P(misses == 2) = nCr(note, max_misses)/sum(m: 0 -> max_misses, nCr(note, m)) = nCr(1, 2)/(nCr(1, 0) + nCr(1, 1) + nCr(1, 2)) = 0.0
        Total:  [0], [1] -> 2
            0 miss = [0] -> 1/2
              0.0
            1 miss = [1] -> 1/2
              1.0
            2 miss = []  -> 0/2
              0.0

        0.0/(0.0 + 1.0 + 0.0) = 0.0

    not (P(miss note 1) and P(miss == 1)) = 1 - (1.0*0.0) = 1.0 - 0.0 = 1.0

    E[100, P(pass note 2)] = E[100, 1.0]

prob of pass note 3 = not (P(hit note 3) and P(miss == 2))
    P(hit note 3) = 0.5
    P(miss note 3) = 1 - P(hit note 3) = 0.5
    P(misses == 2) = nCr(note, max_misses)/sum(m: 0 -> max_misses, nCr(note, m)) = nCr(2, 2)/(nCr(2, 0) + nCr(2, 1) + nCr(2, 2)) = 0.25
        Total:  [0, 0], [1, 0], [0, 1], [1, 1]  -> 4
            0 miss = [0, 0] -> 1/4
                0.0*0.0 * 1 = 0.0
            1 miss = [1, 0], [0, 1] -> 2/4 = 1/2
                0.0*0.0 = 0.0
                0.0*0.0 = 0.0
            2 miss = [1, 1] -> 1/4
                1.0*1.0 * 1 = 1.0

            1.0/(0.0 + 0.0 + 0.0 + 1.0) = 1.0

    not (P(miss note 3) and P(miss == 2)) = 1 - (0.5*1.0) = 1 - 0.5 = 0.5

    E[100, P(pass note 3)] = E[100, 0.5‬]


@ max_misses = 2
hit 1: 100 ms     Time to fc: E[  0 + 100, 0.0] = 100   actual: 100 ms
hit 2: 100 ms     Time to fc: E[100 + 100, 0.0] = 200   actual: 200 ms
hit 3: 100 ms     Time to fc: E[200 + 100, 0.5‬] = 600   actual: 600 ms

total FC time:
    E[ E[ E[ 100, 1.0 ] + 100, 1.0 ] + 100, 0.5‬ ]

    ((100/1.0 + 100)/1.0 + 100)/0.5
    (100 + 100*1.0 + 100*1.0*1.0)/(1.0*1.0*0.5)
    100/(1.0*1.0*0.5) + 100/(1.0*0.5) + 100/(0.5)
    100/0.5 + 100/0.5 + 100/0.5
'''

main.py

"""
Authors: abraker, Joseph Ireland, Tr3sLeches
"""

import random
import itertools
from scipy.special import comb  # nCr
from scipy.special import factorial
from scipy.stats import norm
import numpy as np
import time
import datetime
from poibin import PoiBin


class PoissonBinomialApprox:

    def __init__(self):
        self.mu    = 0
        self.var   = 0
        self.gamma = 0


    def set_p(self, pList):
        self.mu    = np.sum(pList)
        fail       = 1 - pList
        self.var   = np.sum(pList * fail)
        self.gamma = np.sum(pList * fail * (1 - 2*pList))


    def add(self, p):
        self.mu    += p
        self.var   += p * (1 - p)
        self.gamma += p * (1 - p) * (1 - 2*p)


    def cdf(self):
        sigma = np.sqrt(self.var)
        v     = self.gamma/(6 * sigma**3)
        mu    = self.mu

        def func(n):
            k = (n + 0.5 - mu)/sigma
            result = norm.cdf(k) + v*(1 - k*k) * norm.pdf(k)
            return max(0, min(1, result))

        return func


class PoissonBinomialCumulative:

    def __init__(self, pList):
        self.mu    = np.add.accumulate(pList)                           # mean
        fail       = 1 - pList
        self.var   = np.add.accumulate(pList * fail)                    # variance
        self.gamma = np.add.accumulate(pList * fail * (1 - 2*pList))
        self.sigma = np.sqrt(self.var)
        self.v     = self.gamma/(6 * self.sigma**3)             # skewness


    def cdf(self, n):
        k = (n + 0.5 - self.mu)/self.sigma
        result = norm.cdf(k) + self.v*(1 - k*k) * norm.pdf(k)
        return np.maximum(0, np.minimum(1, result))


class ApproxPoiBin():
    
    def __init__(self, probs):
        self.probs = probs

        self.mean     = np.sum(self.probs)
        self.var      = np.sum(self.probs * (1 - self.probs))
        self.skew     = np.sum(self.probs * (1 - self.probs) * (1 - 2*self.probs))/(np.sqrt(self.var)**3)
        self.kurtosis = np.sum(self.probs * (1 - 6*(1 - self.probs)*self.probs)/(np.sqrt(self.var)**4))

        self.offset = ApproxPoiBin.mid_cdf(-self.skew/2)
        self.model_skew = -self.skew/2


    def cdf(self, x):
        x = (1/self.var**0.5)*(x - self.mean + 0.5) + self.offset
        a = self.model_skew

        asq = (a*a)**(1/2)
        ax  = (a*a*x*x)
        a1  = a*a - 1

        term1 = 1/(x*x + 1)**0.5
        term2 = (a*(ax + 1)**0.5)/a1 + x
        term3 = (a*asq)/a1

        return 0.5*(term1*term2 - term3 + 1)


class PoiBin2():
    
    def __init__(self, probs):
        omega = 2 * np.pi / (len(probs) + 1)

        half_number_trials = int(len(probs) / 2 + len(probs) % 2)
        num_values = min(250, half_number_trials)

        idx_array = np.arange(1, num_values + 1)
        exp_value = np.exp(omega * idx_array * 1j)

        # Multiplying each probability by each complex number forming an arc of a circle
        xy = 1 - probs + probs*exp_value.reshape(num_values, 1)
        
        argz_sum = np.arctan2(xy.imag, xy.real).sum(axis=1)

        d_value = np.zeros(num_values)
        d_value = np.abs(xy).prod(axis=1)

        chi = np.zeros(len(probs) + 1, dtype=complex)
        chi[0] = 1
        chi[1:num_values + 1] = d_value*np.exp(argz_sum * 1j)
        chi[half_number_trials + 1:len(probs) + 1] = np.conjugate(chi[1:len(probs) - half_number_trials + 1] [::-1])
        chi /= len(probs) + 1

        self.pmf = np.fft.fft(chi).real
        self.cmf = np.add.accumulate(self.pmf)


    def pdf(self, x):
        return self.pmf[x]


    def cdf(self, x):
        return self.cmf[x]


    @staticmethod
    def mid_cdf(a):
        term1 = a*a - a**6
        term2 = 2*abs(a**3)*(a**6 - a**4 - a*a + 1)**0.5
        term3 = 3*a**6 - a**4 - 3*a*a + 1

        if  1 < a:      return  ((term1 + term2)/term3)**0.5
        if  0 <= a < 1: return  ((term1 - term2)/term3)**0.5
        if -1 < a < 0:  return -((term1 - term2)/term3)**0.5
        if  a < -1:     return -((term1 + term2)/term3)**0.5

        return np.nan


class PoissonBinomialTr3sleches:
    def __init__(self, pList):
        self.mu    = np.add.accumulate(pList)                           # mean
        self.var   = np.add.accumulate(pList * (1 - pList))             # variance
        self.sigma = np.sqrt(self.var)                          # standard deviation
        self.theta = self.mu/self.sigma*np.sqrt(self.mu)
        self.lam   = 1 - np.sqrt(self.mu/self.var)


    def cdf(self, x):
        def pdf(n):
            result = self.theta*np.power(self.theta + self.lam*n, n - 1)*np.exp(-self.theta - self.lam*n)/factorial(n)
            return np.maximum(0, np.minimum(1, result))

        result = 0
        for n in range(x + 1): 
            result += pdf(n)

        return np.maximum(0, np.minimum(1, result))
        


# Used for debugging
def calc_one_hit(time, hit_prob):
    return time/hit_prob


# Unused
def calc_prob_misses(note, misses):
    # nCr(note, misses)/sum(m: 0 -> misses, nCr(note, m))
    return comb(note, misses)/np.sum(comb(note, np.arange(misses + 1)))


def calc_prob_misses_orig(hit_probs, misses):
    total_mask  = np.asarray(list(map(list, itertools.product([0, 1], repeat=len(hit_probs)))))
    total_mask  = total_mask[np.count_nonzero(total_mask, axis=1) <= misses]

    misses_mask = total_mask[np.count_nonzero(total_mask, axis=1) == misses]

    prob_total  = np.sum(np.prod(np.abs(total_mask - hit_probs), axis=1))
    prob_misses = np.sum(np.prod(np.abs(misses_mask - hit_probs), axis=1))

    try: return prob_misses/prob_total
    except ZeroDivisionError: return 1.0


def calc_prob_misses_exact(hit_probs, misses):
    if len(hit_probs) < misses: return 0
    elif misses == 0:           return 1

    pb = PoiBin(1 - hit_probs)
    prob_reach_note    = pb.cdf([misses])[0]
    prob_at_max_misses = pb.pmf([misses])[0]

    try: return prob_at_max_misses/prob_reach_note
    except ZeroDivisionError: return 1.0


def calc_prob_misses_approx(hit_probs, misses):
    if len(hit_probs) < misses: return 0
    elif misses == 0:           return 1
    
    pb = PoissonBinomialApprox()
    pb.set_p(1 - hit_probs)
    cdf = pb.cdf()
    
    prob_reach_note    = cdf(misses)
    prob_at_max_misses = prob_reach_note - cdf(misses - 1)
   
    try: return prob_at_max_misses/prob_reach_note
    except ZeroDivisionError: return 1.0


def calc_all_prob_misses_approx(hit_probs, misses):
    pb = PoissonBinomialApprox()
    yield 0

    for i, prob in enumerate(hit_probs[:-1]):
        pb.add(1 - prob)
        if (i + 1 < misses):
          yield 0
          continue

        cdf = pb.cdf()
        prob_reach_note    = cdf(misses)
        prob_at_max_misses = prob_reach_note - cdf(misses - 1)
        
        try: yield prob_at_max_misses/prob_reach_note
        except ZeroDivisionError: yield 1.0


def calc_all_prob_misses_numpy(hit_probs, misses):
    pb = PoissonBinomialCumulative(1 - hit_probs[:-1])
    
    prob_reach_note    = pb.cdf(misses)
    #print('prob_reach_note', prob_reach_note)
    prob_at_max_misses = prob_reach_note - pb.cdf(misses - 1)
    #print('prob_at_max_misses', prob_at_max_misses)
    result = prob_at_max_misses / prob_reach_note
    result[np.isnan(result)] = 1
    #print('result', result)
    
    result[:misses - 1] = 0
    result = np.insert(result, 0, 0)
    return result


def calc_all_prob_misses_numpy_tr3sLeches(hit_probs, misses):
    pb = PoissonBinomialTr3sleches(1 - hit_probs[:-1])
    
    prob_reach_note    = pb.cdf(misses)
    #print('prob_reach_note', prob_reach_note)
    prob_at_max_misses = prob_reach_note - pb.cdf(misses - 1)
    #print('prob_at_max_misses', prob_at_max_misses)
    result = prob_at_max_misses / prob_reach_note
    result[np.isnan(result)] = 1
    #print('result', result)
    
    result[:misses - 1] = 0
    result = np.insert(result, 0, 0)
    return result


def calc_prob_reach_note_exact(miss_probs, misses):
    if len(miss_probs) < misses: return 1
    return PoiBin(miss_probs).cdf([misses])[0]  # prob_reach_note


def calc_prob_reach_note_poibin_approx(miss_probs, misses):
    if len(miss_probs) < misses: return 1
    return PoiBin2(miss_probs).cdf(misses)  # prob_reach_note


def calc_prob_reach_note_exact_numpy(miss_probs, misses):
    if miss_probs.shape[0] < misses: return 1
    return PoiBin(miss_probs).cdf([misses])[0]  # prob_reach_note


def calc_fc_time(times, hit_probs, max_misses, prob_misses_at_max):
    prob_fail    = (1 - hit_probs)*prob_misses_at_max
    prob_success = 1 - prob_fail

    # [::-1] reverses the array
    prob_success_til_end = np.cumprod(prob_success[::-1])[::-1]
    return np.sum(times/prob_success_til_end)


def calc_fc_time_orig(times, hit_probs, max_misses=0):
    prob_misses_at_max = [ calc_prob_misses_orig(hit_probs[:k], max_misses) for k in range(len(hit_probs)) ] if max_misses else 1
    return calc_fc_time(times, hit_probs, max_misses, prob_misses_at_max)


def calc_fc_time_exact(times, hit_probs, max_misses=0):
    prob_misses_at_max = [ calc_prob_misses_exact(hit_probs[:k], max_misses) for k in range(len(hit_probs)) ] if max_misses else 1
    return calc_fc_time(times, hit_probs, max_misses, prob_misses_at_max)


def calc_fc_time_from_playtime_exact(times, miss_probs, max_misses):
    prob_reach_note = [ calc_prob_reach_note_exact(miss_probs[:k], max_misses) for k in range(len(miss_probs)) ]
    prob_success    = PoiBin(miss_probs).cdf([max_misses])[0]

    return np.sum(times*prob_reach_note)/prob_success


def calc_fc_time_from_playtime_poibin_approx(times, miss_probs, max_misses):
    #prob_reach_note = [ calc_prob_reach_note_poibin_approx(miss_probs[:k], max_misses) for k in range(len(miss_probs)) ]
    prob_reach_note = [ (PoiBin2(miss_probs[:k]).cdf(max_misses) if k >= misses else 1) for k in range(len(hit_probs)) ]
    prob_success    = PoiBin2(miss_probs).cdf(max_misses)
    
    return np.sum(times*prob_reach_note)/prob_success


def calc_fc_time_from_playtime_poibin_approx_numba(times, miss_probs, max_misses):
    #prob_reach_note = [ calc_prob_reach_note_poibin_approx(miss_probs[:k], max_misses) for k in range(len(miss_probs)) ]
    prob_reach_note = [ (PoiBin2(miss_probs[:k]).cdf(max_misses) if k >= misses else 1) for k in range(len(hit_probs)) ]
    prob_success    = PoiBin2(miss_probs).cdf(max_misses)
    
    return np.sum(times*prob_reach_note)/prob_success


def calc_fc_time_approx(times, hit_probs, max_misses=0):
    prob_misses_at_max = [ calc_prob_misses_approx(hit_probs[:k], max_misses) for k in range(len(hit_probs)) ] if max_misses else 1
    return calc_fc_time(times, hit_probs, max_misses, prob_misses_at_max)


def calc_fc_time_from_playtime_numpy_exact(times, miss_probs, max_misses):
    def get_missprobs(data):
        data = np.vstack([data]*(len(data) + 1))
        out  = np.zeros((len(data) - 1, len(data) - 1))
        inds = np.mask_indices(data.shape[1], np.tril, -1)

        out[inds] = data[inds]
        return out

    prob_reach_note = np.apply_along_axis(calc_prob_reach_note_exact_numpy, 1, get_missprobs(miss_probs), max_misses)
    prob_success    = PoiBin(miss_probs).cdf([max_misses])[0]

    return np.sum(times*prob_reach_note)/prob_success


def calc_fc_time_approx_fast(times, hit_probs, max_misses=0):
    prob_misses_at_max = [ p for p in calc_all_prob_misses_approx(hit_probs, max_misses) ] if max_misses else 1
    return calc_fc_time(times, hit_probs, max_misses, prob_misses_at_max)


def calc_fc_time_approx_numpy(times, hit_probs, max_misses=0):
    prob_misses_at_max = calc_all_prob_misses_numpy(hit_probs, max_misses) if max_misses else 1
    return calc_fc_time(times, hit_probs, max_misses, prob_misses_at_max)


def calc_fc_time_approx_numpy_tr3sleches(times, hit_probs, max_misses=0):
    prob_misses_at_max = calc_all_prob_misses_numpy_tr3sLeches(hit_probs, max_misses) if max_misses else 1
    return calc_fc_time(times, hit_probs, max_misses, prob_misses_at_max)


def calc_fc_time_from_playtime(times, miss_probs, max_misses):
    pb  = PoissonBinomialCumulative(miss_probs)
    cdf = pb.cdf(misses)

    prob_success    = cdf[-1]
    prob_reach_note = cdf[:-1]

    prob_reach_note[:misses-1] = 1
    prob_reach_note = np.insert(prob_reach_note, 0, 1)

    return np.sum(times*prob_reach_note)/prob_success


def sim(times, hit_probs, max_misses=0):
    # Simulation
    random.seed()
    trials = 1000

    time_wait_avg = 0
    trial  = 0
    misses = 0

    while trial < trials:
        time_wait = 0
        note      = 0
        misses    = 0

        while note < len(hit_probs):
            time_wait += times[note]

            if random.random() > hit_probs[note]:
                misses += 1
                if misses > max_misses:
                    misses = 0
                    note = 0
                    continue
                
            note += 1
          
        time_wait_avg += time_wait
        trial += 1

    return time_wait_avg/trials


def time_it(f, map_length, *args, **kwargs):
    start      = time.time()
    time_to_fc = f(*args,**kwargs)
    retries    = 'N/A'
    
    if not np.isnan(time_to_fc):
      try: 
          time_to_fc = datetime.timedelta(milliseconds=time_to_fc) 
          retries = round(time_to_fc/map_length, 1)
      except OverflowError: 
          time_to_fc = 'FOREVER'
          retries = 'Infinite'

    print(f'{f.__name__}: {time_to_fc} ({retries}x),   elapsed: {time.time() - start}')


times     = np.repeat(np.asarray([ 100, 100, 100 ]), 300)
hit_probs = np.repeat(np.asarray([ .995, .99, .98 ]), 300)

misses = 3


try: map_length = datetime.timedelta(milliseconds=float(np.sum(times)))
except OverflowError: map_length = 'Too damn long'
print(f'Map length: {map_length},   Notes: {len(times)}')


#time_it(sim, times, map_length, hit_probs, max_misses=misses)                                 # Simulation, used as ground truth
#time_it(calc_fc_time_orig, map_length, times, hit_probs, max_misses=misses)                   # Exact calculation, 100% trusted to be correct

#time_it(calc_fc_time_from_playtime_numpy_exact, map_length, times, 1 - hit_probs, max_misses=misses) # Exact calculation, 100% trusted to be correct
###time_it(calc_fc_time_exact, map_length, times, hit_probs, max_misses=misses)                   # Exact calculation, 100% trusted to be correct
#time_it(calc_fc_time_from_playtime_exact, map_length, times, 1 - hit_probs, max_misses=misses) # Exact calculation, 100% trusted to be correct
time_it(calc_fc_time_from_playtime_poibin_approx, map_length, times, 1 - hit_probs, max_misses=misses) # Approximate estimation, trusted to be acceptable for most cases

time_it(calc_fc_time_approx, map_length, times, hit_probs, max_misses=misses)                  # Approximate estimation, trusted to be acceptable for most cases
time_it(calc_fc_time_approx_fast, map_length, times, hit_probs, max_misses=misses)             # Approximate estimation, trusted to be acceptable for most cases
time_it(calc_fc_time_approx_numpy, map_length, times, hit_probs, max_misses=misses)            
time_it(calc_fc_time_approx_numpy_tr3sleches, map_length, times, hit_probs, max_misses=misses)  
#time_it(calc_fc_time_from_playtime, map_length, times, 1 - hit_probs, max_misses=misses)       # Rough estimation, I wouldn't trust this

poibin.py

# -*- coding: utf-8 -*-
"""
Created on Tue Mar 29, 2016
Module:
    poibin - Poisson Binomial Distribution
Author:
    Mika Straka
Description:
    Implementation of the Poisson Binomial distribution for the sum of
    independent and not identically distributed random variables as described
    in the reference [Hong2013]_.
    Implemented method:
        * ``pmf``: probability mass function
        * ``cdf``: cumulative distribution function
        * ``pval``: p-value (1 - cdf)
Usage:
    Be ``p`` a list or  NumPy array of success probabilities for ``n``
    non-identically distributed Bernoulli random variables.
    Import the module and create an instance of the distribution with::
        >>> from poibin import PoiBin
        >>> pb = PoiBin(p)
    Be ``x`` a list or NumPy array of different number of successes.
    To obtain the:
    * probability mass function of x, use::
        >>> pb.pmf(x)
    * cumulative distribution function of x, use::
        >>> pb.cdf(x)
    * p-values of x, use::
        >>> pb.pval(x)
    The functions are applied component-wise and a NumPy array of the same
    length as ``x`` is returned.
References:
.. [Hong2013] Yili Hong, On computing the distribution function for the Poisson
    binomial distribution,
    Computational Statistics & Data Analysis, Volume 59, March 2013,
    Pages 41-51, ISSN 0167-9473,
    http://dx.doi.org/10.1016/j.csda.2012.10.006.

MIT License

Copyright (c) 2016 Mika J. Straka

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

"""

import collections
import numpy as np


class PoiBin(object):
    """Poisson Binomial distribution for random variables.
    This class implements the Poisson Binomial distribution for Bernoulli
    trials with different success probabilities. The distribution describes
    thus a random variable that is the sum of independent and not identically
    distributed single Bernoulli random variables.
    The class offers methods for calculating the probability mass function, the
    cumulative distribution function, and p-values for right-sided testing.
    """

    def __init__(self, probabilities):
        """Initialize the class and calculate the ``pmf`` and ``cdf``.
        :param probabilities: sequence of success probabilities :math:`p_i \\in
            [0, 1] \\forall i \\in [0, N]` for :math:`N` independent but not
            identically distributed Bernoulli random variables
        :type probabilities: numpy.array
        """
        self.success_probabilities = np.array(probabilities)
        self.number_trials = self.success_probabilities.size
        self.check_input_prob()
        self.omega = 2 * np.pi / (self.number_trials + 1)
        self.pmf_list = self.get_pmf_xi()
        self.cdf_list = self.get_cdf(self.pmf_list)

# ------------------------------------------------------------------------------
# Methods for the Poisson Binomial Distribution
# ------------------------------------------------------------------------------

    def pmf(self, number_successes):
        """Calculate the probability mass function ``pmf`` for the input values.
        The ``pmf`` is defined as
        .. math::
            pmf(k) = Pr(X = k), k = 0, 1, ..., n.
        :param number_successes: number of successful trials for which the
            probability mass function is calculated
        :type number_successes: int or list of integers
        """
        self.check_rv_input(number_successes)
        return self.pmf_list[number_successes]

    def cdf(self, number_successes):
        """Calculate the cumulative distribution function for the input values.
        The cumulative distribution function ``cdf`` for a number ``k`` of
        successes is defined as
        .. math::
            cdf(k) = Pr(X \\leq k), k = 0, 1, ..., n.
        :param number_successes: number of successful trials for which the
            cumulative distribution function is calculated
        :type number_successes: int or list of integers
        """
        self.check_rv_input(number_successes)
        return self.cdf_list[number_successes]

    def pval(self, number_successes):
        """Return the p-values corresponding to the input numbers of successes.
        The p-values for right-sided testing are defined as
        .. math::
            pval(k) = Pr(X \\geq k ),  k = 0, 1, ..., n.
        .. note::
            Since :math:`cdf(k) = Pr(X <= k)`, the function returns
            .. math::
                1 - cdf(X < k) & = 1 - cdf(X <= k - 1)
                               & = 1 - cdf(X <= k) + pmf(X = k),
                               k = 0, 1, .., n.
        :param number_successes: number of successful trials for which the
            p-value is calculated
        :type number_successes: int, numpy.array, or list of integers
        """
        self.check_rv_input(number_successes)
        i = 0
        try:
            isinstance(number_successes, collections.Iterable)
            pvalues = np.array(number_successes, dtype='float')
            # if input is iterable (list, numpy.array):
            for k in number_successes:
                pvalues[i] = 1. - self.cdf(k) + self.pmf(k)
                i += 1
            return pvalues
        except TypeError:
            # if input is an integer:
            if number_successes == 0:
                return 1
            else:
                return 1 - self.cdf(number_successes - 1)

# ------------------------------------------------------------------------------
# Methods to obtain pmf and cdf
# ------------------------------------------------------------------------------

    def get_cdf(self, event_probabilities):
        """Return the values of the cumulative density function.
        Return a list which contains all the values of the cumulative
        density function for :math:`i = 0, 1, ..., n`.
        :param event_probabilities: array of single event probabilities
        :type event_probabilities: numpy.array
        """
        cdf = np.empty(self.number_trials + 1)
        cdf[0] = event_probabilities[0]
        for i in range(1, self.number_trials + 1):
            cdf[i] = cdf[i - 1] + event_probabilities[i]
        return cdf

    def get_pmf_xi(self):
        """Return the values of the variable ``xi``.
        The components ``xi`` make up the probability mass function, i.e.
        :math:`\\xi(k) = pmf(k) = Pr(X = k)`.
        """
        chi = np.empty(self.number_trials + 1, dtype=complex)
        chi[0] = 1
        half_number_trials = int(
            self.number_trials / 2 + self.number_trials % 2)
        # set first half of chis:
        chi[1:half_number_trials + 1] = self.get_chi(
            np.arange(1, half_number_trials + 1))
        # set second half of chis:
        chi[half_number_trials + 1:self.number_trials + 1] = np.conjugate(
            chi[1:self.number_trials - half_number_trials + 1] [::-1])
        chi /= self.number_trials + 1
        xi = np.fft.fft(chi)
        if self.check_xi_are_real(xi):
            xi = xi.real
        else:
            raise TypeError("pmf / xi values have to be real.")
        xi += np.finfo(type(xi[0])).eps
        return xi

    def get_chi(self, idx_array):
        """Return the values of ``chi`` for the specified indices.
        :param idx_array: array of indices for which the ``chi`` values should
            be calculated
        :type idx_array: numpy.array
        """
        # get_z:
        exp_value = np.exp(self.omega * idx_array * 1j)
        xy = 1 - self.success_probabilities + \
            self.success_probabilities * exp_value[:, np.newaxis]
        # sum over the principal values of the arguments of z:
        argz_sum = np.arctan2(xy.imag, xy.real).sum(axis=1)
        # get d value:
        exparg = np.log(np.abs(xy)).sum(axis=1)
        d_value = np.exp(exparg)
        # get chi values:
        chi = d_value * np.exp(argz_sum * 1j)
        return chi

# ------------------------------------------------------------------------------
# Auxiliary functions
# ------------------------------------------------------------------------------

    def check_rv_input(self, number_successes):
        """Assert that the input values ``number_successes`` are OK.
        The input values ``number_successes`` for the random variable have to be
        integers, greater or equal to 0, and smaller or equal to the total
        number of trials ``self.number_trials``.
        :param number_successes: number of successful trials
        :type number_successes: int or list of integers """
        try:
            for k in number_successes:
                assert (type(k) == int or type(k) == np.int64), \
                        "Values in input list must be integers"
                assert k >= 0, 'Values in input list cannot be negative.'
                assert k <= self.number_trials, \
                    'Values in input list must be smaller or equal to the ' \
                    'number of input probabilities "n"'
        except TypeError:
            assert (type(number_successes) == int or \
                type(number_successes) == np.int64), \
                'Input value must be an integer.'
            assert number_successes >= 0, "Input value cannot be negative."
            assert number_successes <= self.number_trials, \
                'Input value cannot be greater than ' + str(self.number_trials)
        return True

    @staticmethod
    def check_xi_are_real(xi_values):
        """Check whether all the ``xi``s have imaginary part equal to 0.
        The probabilities :math:`\\xi(k) = pmf(k) = Pr(X = k)` have to be
        positive and must have imaginary part equal to zero.
        :param xi_values: single event probabilities
        :type xi_values: complex
        """
        return np.all(xi_values.imag <= np.finfo(float).eps)

    def check_input_prob(self):
        """Check that all the input probabilities are in the interval [0, 1]."""
        if self.success_probabilities.shape != (self.number_trials,):
            raise ValueError(
                "Input must be an one-dimensional array or a list.")
        if not np.all(self.success_probabilities >= 0):
            raise ValueError("Input probabilities have to be non negative.")
        if not np.all(self.success_probabilities <= 1):
            raise ValueError("Input probabilities have to be smaller than 1.")

################################################################################
# Main
################################################################################

if __name__ == "__main__":
    pass

requirements.txt

numpy < 1.24
scipy