kt.pyx

from numpy cimport *
from numpy import empty, zeros, ones, corrcoef, tile, concatenate, fromiter, uint32, uint8, double, array, arange, dot, finfo, seterr, ndindex, r_, repeat, copy, loadtxt, outer, random, vstack, hstack, zeros_like, ones_like, empty_like, vectorize, transpose, pi, linspace, e, ceil, nan
import numpy as np
seterr(over='ignore', invalid='ignore', divide='ignore')
from scipy.special import gammainc, gamma, erf
import functools

cdef extern from "c/exponential.h":
    double exponential()
    void mt_init()

mt_init()               # see final.h

cdef extern from "math.h":
    double log( double x)
    double sqrt( double x)
    double exp( double x)
    double pow(double x, double y)

def getUdUp(u=1e-8, **otherParams):
    Tdf = 2*700
    Tpf = 2*5e6
    TdTp = otherParams.pop('TdTp', Tdf/Tpf)
    Td = sqrt(TdTp*Tpf*Tdf)
    Tp = sqrt(Tdf*Tpf/TdTp)
    return dict(Ud=u*Td, Up=u*Tp, **otherParams)

class first:
    Ud = 1e-8*700
    Up = 1e-8*5e6
    sd = 0.1
    sp = 1e-3
    N_0 = 1e3
    t_max = 10000
    verbose = False

    def p(self, txt):
        if self.verbose: print txt

    def __init__(self, **params):
        self.__dict__.update(params)

    @property
    def f(self): return self.Ud*self.sd/(1 + self.sd)

    @property
    def dN(self): return self.sd

    @property
    def vp(self): return self.sp*self.Up

    @property
    def N_eq(self): 
        assert self.vp > 0, 'vp: {x.vp:g}, Up: {x.Up:g}, sp: {x.sp:g}, N: {x.N_0}'.format(x = self )
        assert self.f > 0,  ' f: {x.f:g} , Up: {x.Up:g}, sp: {x.sp:g}, N: {x.N_0}'.format(x = self )
        assert self.dN > 0, 'dN: {x.dN:g}, Up: {x.Up:g}, sp: {x.sp:g}, N: {x.N_0}'.format(x = self )
        return self.vp/(self.f*log(1 + self.dN))
    
    def __repr__(self):
        return """Ud =    {x.Ud:g}
Up =    {x.Up:g}
sd =    {x.sd:g}
sp =    {x.sp:g}
N_0 =   {x.N_0:g}
t_max = {x.t_max:g}
N_eq =  {x.N_eq:g}
""".format(x=self)

    def velocity(model, x): return x*model.N_eq*model.vp*(x - 1)

    def potential(model, x): return model.N_eq*model.vp*x*x*(0.5 - x/3)

    def pC(model, x=None):
        D = 2/log(1+model.dN)
        if x is None:
            x = model.N_0/model.N_eq
        return 1 - gammainc(D, D/x) 

def clonalInterference(model):
    """ CI -> d<d>/dt : sp = 0, Ud ~ sd"""
    from scipy.optimize import fsolve
    cdef double sd = model.sd
    def f(V):
        A = V/model.Ud; B = log(A);
        return model.N - (exp(V/2/sd*((B-1)**2+1)-0.5*(log(sd**3/(A*B))- B if V>=sd else 0) ) if V>sd/B else 0.5*A/sd)
    return fsolve(f, sd)

@functools.lru_cache(maxsize=None)
def fAnddN(double Up, double sd, double sp, double N):
    from scipy.stats import poisson
    cdef:
        double L = Up/sp
        int i
        int i_max = int(ceil(log(1+sd)/log(1+sp)))
        int k = (poisson.cdf(arange(L), L) < 0.5*pow(N,-0.5)).sum()
        ndarray[int64_t,  ndim=1] X = arange(i_max)
        ndarray[double_t, ndim=1] W = (1+sd)*(1+sp)**-X
        ndarray[double_t, ndim=1] f = poisson.pmf(X+k, L)/(1-poisson.cdf(k-1, L))
        ndarray[double_t, ndim=1] fP = f*(W-1)/W
        ndarray[double_t, ndim=1] dNp = empty(i_max)
        
    Lambda = fP.sum()
    assert Lambda < sd/(1+sd), 'wtf!'
    
    dNp[i_max - 1] = i_max - 1
    for i in range(i_max-2, -1, -1): dNp[i] = (i + sp*dNp[i+1])/(1 + sp)
    V = (1 + sd)*(1 + sp)**-dNp
    assert (V <= 1+sd).all()
    return Lambda, fP.dot(V)/Lambda - 1
    
from scipy.integrate import simps
def moran(double Up, double sp, double N, int r1=150, int r2=600):
    cdef double x0 = exp(-Up/sp), spp =sp/(e-1)
    a, x1, x2 = N*spp/(2*x0), linspace(x0/r1, x0, r1), linspace(x0, 1, r2)
    erfi = r_[erf(1j*sqrt(a)*(x1-x0-x1[0])).imag,0]
    I, G1, G2 = sqrt(pi/(4*a))*exp(-a*x0**2)*(erfi[1:]-erfi[0]), np.exp(a*x1*(x1-2*x0)), np.exp(a*x2*(x2-2*x0))
    Tclick = max(1/sp*(1-e*spp/Up), 0)
    Tclick += N/2*(simps(r_[1,I/(x1*G1)],dx=x1[0])+I[-1]*simps(1./(x2*G2),dx=x2[1]-x2[0]))
    return Up*1/Tclick

def moran(double Up, double sp, double N, int r1=150, int r2=600):
    cdef double x0 = exp(-Up/sp), spp =sp/(e-1)
    Tclick = exp(0.5*spp*N*x0)/spp
    return Up*1/Tclick

#def moran(double Up, double sp, double N, double alpha=1/(e-1)):
#    cdef:
#        double No = N*exp(-Up/sp)
#        double Nosa = No*sp*alpha
#    return sqrt(Nosa*alpha*alpha/pi)*exp(-Nosa)

def wave(double Up, double sp, double N, double v_0=1e-12):
    from scipy.optimize import brentq
    x2 = sp*log(N)
    if Up < x2:
        No = N*exp(-Up/sp)
        pi_p = sp/((1 + sp)**No - 1) #if sp > 5e-5 else 1/(No*(1+2*sp))
        #y1 = Up*N*sp/((1 + sp)**N - 1)    
        #y2 = wave(x2, sp, N)*Up/x2
        assert No*Up*pi_p > 0
        return  No*Up*pi_p #y2
    
    cdef:
        double iL = sp/Up
        double Nd2spSqrtiL = N/2*sp*sqrt(iL)
        double fiveSixthsiL = 5*iL/6
    def f(double v):
        cdef double A = log(e*Up/v), C = v/Up 
        return 1 - 0.5*C*(A*A + 1) - iL*log(Nd2spSqrtiL*sqrt(v*C*C/(Up - v))*A/(1 - C*A + fiveSixthsiL))
    cdef double vp = brentq(f, v_0, Up - v_0)
    assert vp > 0
    return min(vp, Up)

@functools.lru_cache(maxsize=None)
def Vp(double Up, double sp, double N):
    cdef double L = Up/sp
    if L > 1: 
        return wave(Up, sp, N)
    N_0 = N*exp(-L)
    return moran(Up, sp, N) if sp*N_0 > 1 else (1-L)*moran(Up, sp, N) + L*wave(Up, sp, N)

class second(first):
    @property
    def f(self): 
        print('New f')
        return fAnddN(self.Up, self.sd, self.sp, self.N_0)[0]*self.Ud

    @property
    def dN(self): return fAnddN(self.Up, self.sd, self.sp, self.N_0)[1]

    @property
    def vp(self): 
        print('New vp')
        return self.sp*Vp(self.Up, self.sp, self.N_0)

@functools.lru_cache(maxsize=None)
def simplified_fanddN(double sd, double sp, double Up, double N_0):
    from scipy.stats import poisson
    cdef:
        double Lambda = Up/sp
        ndarray[int64_t,  ndim=1] I = arange( int(sd/sp) + 1 )
        int k = (poisson.cdf(arange(int(Lambda)), Lambda) < 0.5*pow(N_0, -0.5)  ).sum()
        ndarray[double_t, ndim=1] f = poisson.pmf(I + k, Lambda)
        ndarray[double_t, ndim=1] W = sd - sp*I
    return dot(f, W/(1+W)), dot(f, W)/f.sum() 

class third(first):    
    @property
    def f(self): return self.Ud*simplified_fanddN(self.sd, self.sp, self.Up, self.N_0)[0]

    @property
    def dN(self): return simplified_fanddN(self.sd, self.sp, self.Up, self.N_0)[1]

    @property
    def vp(self):
        from scipy.stats import poisson
        cdef:
            double Up = self.Up
            double sp = self.sp
            double sd = self.sd
            double Lambda = Up/sp
            ndarray[int64_t,  ndim=1] I = arange( int(sd/sp) + 1 )
            int k = (poisson.cdf(arange(int(Lambda)), Lambda) < 0.5*pow(self.N_0, -0.5)  ).sum()
            double No = self.N_0*poisson.pmf(k, Lambda)            #exp(-Lambda) 
            double pi_p = sp/((1 + sp)**No - 1) if sp > 5e-5 else 1/(No*(1+2*sp))
        return Up*sp*No*pi_p if No > 1e-1 else Up*sp

def tC(model, x, n_low=0.35, n_high=7):
    from scipy.integrate import quadrature
    D = 2/model.dN
    N_eq = model.N_eq
    vp = model.vp
    p = model.pC(x)
    gD = gamma(D)

    def I1(x):
        p = model.pC(x)
        de_mN = (D/x)**D*exp(-D/x)/gD
        return (1 - p)*p/(de_mN*N_eq*N_eq)

    def I2(x):
        p = model.pC(x)
        de_mN = (D/x)**D*exp(-D/x)/gD
        return p*p/(de_mN*N_eq*N_eq)

    if x < n_high:
        C = D*N_eq/vp
        if x < n_low:
            return C*quadrature(I1, n_low, n_high)[0] + log(n_low/x)/vp + x/(D*vp)
        else:
            return C*(quadrature(I1, x, n_high)[0] + (1-p)/p*(gamma(D-1)*(1 - gammainc(D - 1, D/n_low))/(D*N_eq*gD) + quadrature(I2, n_low, x)[0])) + 1/(n_high*vp)
    else:
        return 1/(vp*x)

def hGillespie(model, long trials=10000000, double x_max=20):
    """Hybrid gillespie model"""
    cdef:
        double theta = 1 + model.dN
        double tau_max = model.t_max*model.vp        # tau = t*vp
        double dimmensionless_if = model.vp/model.f/model.N_0    # In units of 
        double x, Q
        long n_cancers = 0
        long i
        ndarray[double_t, ndim=1] Taus = empty(trials, dtype=double)
    for i in range(trials):
        x = 1
        Taus[n_cancers] = 0
        while x < x_max:
            Q = 1 - exponential()*dimmensionless_if/x
            x *= Q*theta
            Taus[n_cancers] -= log(Q)
            if Taus[n_cancers] > tau_max or x < 0: 
                break
        else:    
            n_cancers += 1
    return {'P_cancer':double(n_cancers)/trials, 'generations':Taus[0:n_cancers]/model.vp}

#def hGompGillespie(model, long trials=10000000, double y_max=20):
# Solution to waiting-time jump is transcendental equation, with Exponential Integral function in it!
#    """hybrid gillespie using gomp-ex death rate"""
#
#    cdef:
#        double vp = model.vp
#        double theta = 1 + model.dN
#        double jumpRate = model.f*model.N_0/vp
#        double em1dJump = (e - 1)/jumpRate
#        double omega_d
#        double sigma_d
#        double eSigma_d
#        double t
#    
#    for i in range(trials):
#        y = 1
#        sigma_d = 0
#        omega_d = -sigma_d*vp
#        t = 0
#        eSigma_d = exp(sigma_d)

B2 = array([5.72640947,  -4.84536097])
B3 = array([5.77629529,  -4.94655988,    0.04817712])
B4 = array([10.82164987, -16.63749275,   8.31660758,  -1.75101626])
B5 = array([10.78285714, -16.57518638,   8.32469893,  -1.79888812,   0.01717027])
B6 = array([12.38219463, -22.35712125,  15.97711017,  -6.31622209,   1.13227449,  -0.07230064])

def empiricalPcancer10(y, f=lambda y, B: 1/(1+exp((B*pow(y, arange(len(B)))).sum() )), B=B6):
    """assumes sd=0.1 """
    return f(y, B)




#Pcancer = {0.1: 0.415*N_eq,
#0.2:67, 
#0.05: 0.94*N_eq =
#0.96*N_eq = 
#0.98*N_eq =
#}


#        while y < y_max:
#            dt = log(em1dJump*omega_d*exponential() + 1)/omega_d
#            y *= eTheta 

#            t = log(omega_d*em1/(jumpRate*eSigma_d)*exponential() + exp(omega_d*t))/omega_d
#            sigma_d += theta
#            eSigma_d = nan
#            y *= eSigma 







#def ndnp(double Ud=1e-8, double sd=0.1, double sp=0.001, double K=1e3, double Td=700, double Tp=5000000, double Genome_size = 149948690/4, long trials=50000000, double Nmax=2):
#
#random.poisson(T*u*(Genome_size-2*Tp))
#    cdef double Ud = 2*u*Td, Up = 2*u*Tp
#    out = lambdaG(Ud, Up, sd, sp, No)
#    cdef double vp = sp*Vp(sp, Up, No), L = out[0], dN = out[1]+1, N, Q, T, m_hitchhikers = (sd - out[1])/sp
#    cdef ndarray[uint32_t, ndim=2] Trial = empty((Iter,3),uint32)
#    cdef uint32_t i, j = 0, d
#    for i in range(Iter):
#        N, T, d = No, 0, 0
#        while N > 0 and T < Time:
#            Q = 1 - exponential()*vp/(N*L)
#            N *= Q*dN
#            d += 1
#            T -= log(Q)/vp
#            if N > No*Top:
#                j+=1
#                break
#    print j
#    return Trial[:j,:]