<jupyter_start><jupyter_text>Setup Load packages<jupyter_code>import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from scipy.special import binom
import itertools
import helpers as hp
import simulations as sims
import importlib
importlib.reload(sims)

plt.style.use('ggplot')

FIG_PATH = "./../../figures/"
# copied and pasted from mathematica
def M(t,p,q):
    return -((2**t*(1 - p)**t*(-1 + 2*p - q))/(1 - 2*p + q))
def B(t,p,q):
    return (2*p*(1 - q)**t)/(-1 + 2*p - q) + (2**(1 + t)*(1 - p)**t*p)/(1 - 2*p + q)
def G(t,p,q):
    return (2**t*q)/(1 + q) - (2*p*(1 - q)**t*q)/((-1 + 2*p - q)*(1 + q)) - (2**t*(1 - p)**t*q)/(1 - 2*p + q)
def Sigma(t,p,q):
    return np.array([[2*p*(1-p)*M(t-1,p,q),2*p*(1-p)*M(t-1,p,q),0],\
                     [2*p*(1-p)*M(t-1,p,q),q*(1-q)*B(t-1,p,q)+2*p*(1-p)*M(t-1,p,q),q*(1-q)*B(t-1,p,q)],\
                    [0,q*(1-q)*B(t-1,p,q),q*(1-q)*B(t-1,p,q)]])
def d(p,q):
    return np.array([[2*(1-p),0,0],[2*p,1-q,0],[0,q,2]])
def L(p,q):
    return np.array([[1,1,1],[0,0,1]])

def Omega(t,p,q):
    Om = Sigma(t,p,q)
    for j in range(1,t-1):
        dp = np.linalg.matrix_power(d(p,q),t-j)
        Om = Om + dp @ Sigma(j,p,q) @ np.transpose(dp)
    return Om

def Omega_bar(t,p,q):
    return L(p,q) @ Omega(t,p,q) @ np.transpose(L(p,q))

# this is the complicated analytical formula for the variance conditioned on 
# the total number of cells
def varGN_theory(t,p,q):
    M = Omega_bar(t,p,q)
    return M[1,1]-M[0,1]**2/M[0,0]


def varGN_theory_nobroken(t,p,q):
    if p>0.5:
        return p/(2*p-1)*(2**t)**(1+np.log(1-p)/np.log(2))
    else:
        return p/(1-2*p)*(2**t)**(2+2*np.log(1-p)/np.log(2))
Omega_bar(10,0.2,0.8)
fig,ax = plt.subplots(figsize=(10,5))


ax.semilogy([varGN_theory_nobroken(t,0.2,0.999) for t in range(30)],'C0--',label="$p = 0.2$ (no broken state)")
ax.semilogy([varGN_theory(t,0.2,0.1) for t in range(30)],'C0-',label="$p = 0.2$,$q=0.1$")

ax.semilogy([varGN_theory_nobroken(t,0.4,0.999) for t in range(30)],'C1--',label="$p = 0.4$ (no broken state)")
ax.semilogy([varGN_theory(t,0.4,0.1) for t in range(30)],'C1-',label="$p = 0.4$,$q=0.1$")

ax.semilogy([varGN_theory_nobroken(t,0.6,0.999) for t in range(30)],'C2--',label="$p = 0.6$ (no broken state)")
ax.semilogy([varGN_theory(t,0.6,0.1) for t in range(30)],'C2-',label="$p = 0.6$,$q=0.1$")
ax.set_xlabel(r"$t$")
ax.set_ylabel(r"${\rm var}(G_t|N_t)$")
ax.grid(False)
ax.legend()
plt.savefig(FIG_PATH+"var_GN_broken_state_theory.pdf",bbox_inches = "tight")
[varGN_theory(t,0.7,0.1) for t in range(50)]
qs = np.linspace(0,1,100)
ps = np.linspace(0,1,100)
gr = np.zeros((100,100))
for k in range(100):
    for j in range(100):
        v = varGN_theory(40,qs[k],ps[j])
        if v>1:
            gr[k,j] = 1
        else:
            gr[k,j] = 0
plt.matshow(gr.T)
plt.grid(False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.2588480222068008, 'fraction_numerical': 0.06766134628730049, 'mean_word_length': 3.7339901477832513, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '611077', 'license': 'Unlicense', 'n_tokens_mistral': 1512, 'n_tokens_neox': 1407, 'n_words': 231, 'path': 'code/analysis/.ipynb_checkpoints/broken_state_model_fluctuations-checkpoint.ipynb', 'repo_name': 'elevien/yeast_dna_repair', 'size': '2882'}