import numpy as np
import matplotlib.pyplot as plt
from math import *
import random
import seaborn as sns


"""
Setting the parameters for the process and the algorithm.
aplha and epsilon are the parameters of the OU process.
Ta is the simulation time of the trajectories. T is the averaged time of the observable.
Num_traj is the number of trajectories in the ensemble. dt is the timestep of the simulation.
tau is the cloning time of the GKTL algorithm. k is sampling parameter used.
run_count is the number of independent algorithm runs.
"""
alpha=1
epsilon= 0.5
Ta=100
Num_traj=600
T_obs=50
dt=1/1000
N=int(Ta/dt)
tau=0.5
k=0.3
run_count=100


"""
This function calculates the statistics of the ensemble generated by the algorithm.
Returns the score of the trajectories in the ensemble and the probability associated.
"""
def Statistics(data,Z_value):
    def running_mean(x, N):
        """ x == an array of data. N == number of samples per average """
        cumsum = np.cumsum(np.insert(x, 0, 0))
        return (cumsum[N:] - cumsum[:-N]) / float(N)

    Sc=[]
    for i in range(Num_traj):
        X_mean=running_mean(data[i],int(T_obs/dt))
        prob= (1/Num_traj)*exp(-k*np.mean(data[i])*Ta)*Z_value
        Sc.append([max(X_mean),prob/run_count])

    return(Sc)


"""
Function used to plot the trajectories if needed.
"""
def F_plot(T,data):
    i=0
    while i<Num_traj:
        plt.plot(T,data[i])
        i=i+1


"""
This function scores the trajectories in the ensemble.
Called by the algorithm iteration function to calculate the weights of the trajectories in the ensemble.
"""
def Data_scorer(data,iter):

    norm=0
    score=[]
    integ=0
    start=int((iter-1)*tau/dt)
    stop=int(iter*tau/dt)-1
    print(start,stop)
    for i in range(Num_traj):
        integ = np.mean(data[i][start:stop])*tau
        score.append(exp(k*integ))
        norm=norm+exp(k*integ)

    return(np.array(score),norm)


"""
This function is the algorithm iteration.
Calls the Data_scorer function to calculate the weights. Clones the existing trajectories, proportional to their weight.
Returns the new ensemble.
"""
def GKTL_iter(data,Z,iter):

    new_data=[]
    weight,norm=Data_scorer(data,iter)
    mutation_count=Num_traj*weight/norm
    Z=Z*norm/Num_traj

    if iter==Ta/tau:
        return(Z)

    for i in range(len(mutation_count)):
        mutation_count[i]=mutation_count[i]+random.random()

    mutation_count=np.floor(mutation_count)

    for i in range(Num_traj):
        for j in range(int(mutation_count[i])):
            X_new=data[i].copy()
            m=int((iter)*tau/dt)
            stop=int((iter+1)*tau/dt)
            while m<stop:
                X_new[m] = X_new[m-1] - alpha*X_new[m-1]*dt + sqrt(dt*2*epsilon)*np.random.normal(0,1)
                m=m+1
            new_data.append(X_new)

    S=int(sum(mutation_count))
    if S>Num_traj:
        for c in range(S-Num_traj):
            traj_random=np.random.randint(0,S-c-1)
            new_data.pop(traj_random)

    if S<Num_traj:
        for c in range(Num_traj-S):
            traj_random=np.random.randint(0,S-1)
            traj=new_data[traj_random].copy()
            new_data.append(traj)

    return(new_data,Z)


"""
This function does one run of the entire algorithm and returns the run statistics.
Creates the initial trajectories and iterates the algorithm till time Ta.
"""
def Return_time():

    data=[]
    for i in range(Num_traj):
        X = np.zeros(N)
        X[0]=np.random.normal(0,sqrt(epsilon/alpha))

        end=int(tau/dt)
        for j in range(end):
            X[j] = X[j-1] - alpha*X[j-1]*dt + sqrt(dt*2*epsilon)*np.random.normal(0,1)

        data.append(X)

    iter=1
    Z=1
    while iter<=Ta/tau:
        if iter==Ta/tau:
            Z=GKTL_iter(data,Z,iter)
            break
        data,Z = GKTL_iter(data,Z,iter)
        iter=iter+1

    S=Statistics(data,Z)
    return(S)


"""
Obtaining stats over the different runs of the algorithm.
"""
Stats=[]
for i in range(run_count):
    S=Return_time()
    for j in S:
        Stats.append(j)

Stats.sort()
Stats=np.array(Stats)

"""
Calculating return time for different values of a from the Stats.
"""
A=[]
T=[]
a=0
while a<=max(Stats[:,0]):
    p=0
    for i in range(len(Stats)):
        if Stats[i][0]>a:
            p=p+Stats[i][1]

    A.append(a)
    T.append(-(Ta-T_obs)/log(1-p))
    a=a+0.01


#np.save('a_alg.npy', A)
#np.save('t_alg.npy', T)

plt.plot(np.log10(T),A)
#plt.xlim([0,14])
#plt.ylim([0,1])
plt.grid(True)
plt.show()

64f616d4ecf88ef5d68a4744a8f7eaf769baaf26
<mask token>


def Statistics(data, Z_value):

    def running_mean(x, N):
        """ x == an array of data. N == number of samples per average """
        cumsum = np.cumsum(np.insert(x, 0, 0))
        return (cumsum[N:] - cumsum[:-N]) / float(N)
    Sc = []
    for i in range(Num_traj):
        X_mean = running_mean(data[i], int(T_obs / dt))
        prob = 1 / Num_traj * exp(-k * np.mean(data[i]) * Ta) * Z_value
        Sc.append([max(X_mean), prob / run_count])
    return Sc


<mask token>


def F_plot(T, data):
    i = 0
    while i < Num_traj:
        plt.plot(T, data[i])
        i = i + 1


<mask token>


def Data_scorer(data, iter):
    norm = 0
    score = []
    integ = 0
    start = int((iter - 1) * tau / dt)
    stop = int(iter * tau / dt) - 1
    print(start, stop)
    for i in range(Num_traj):
        integ = np.mean(data[i][start:stop]) * tau
        score.append(exp(k * integ))
        norm = norm + exp(k * integ)
    return np.array(score), norm


<mask token>


def GKTL_iter(data, Z, iter):
    new_data = []
    weight, norm = Data_scorer(data, iter)
    mutation_count = Num_traj * weight / norm
    Z = Z * norm / Num_traj
    if iter == Ta / tau:
        return Z
    for i in range(len(mutation_count)):
        mutation_count[i] = mutation_count[i] + random.random()
    mutation_count = np.floor(mutation_count)
    for i in range(Num_traj):
        for j in range(int(mutation_count[i])):
            X_new = data[i].copy()
            m = int(iter * tau / dt)
            stop = int((iter + 1) * tau / dt)
            while m < stop:
                X_new[m] = X_new[m - 1] - alpha * X_new[m - 1] * dt + sqrt(
                    dt * 2 * epsilon) * np.random.normal(0, 1)
                m = m + 1
            new_data.append(X_new)
    S = int(sum(mutation_count))
    if S > Num_traj:
        for c in range(S - Num_traj):
            traj_random = np.random.randint(0, S - c - 1)
            new_data.pop(traj_random)
    if S < Num_traj:
        for c in range(Num_traj - S):
            traj_random = np.random.randint(0, S - 1)
            traj = new_data[traj_random].copy()
            new_data.append(traj)
    return new_data, Z


<mask token>


def Return_time():
    data = []
    for i in range(Num_traj):
        X = np.zeros(N)
        X[0] = np.random.normal(0, sqrt(epsilon / alpha))
        end = int(tau / dt)
        for j in range(end):
            X[j] = X[j - 1] - alpha * X[j - 1] * dt + sqrt(dt * 2 * epsilon
                ) * np.random.normal(0, 1)
        data.append(X)
    iter = 1
    Z = 1
    while iter <= Ta / tau:
        if iter == Ta / tau:
            Z = GKTL_iter(data, Z, iter)
            break
        data, Z = GKTL_iter(data, Z, iter)
        iter = iter + 1
    S = Statistics(data, Z)
    return S


<mask token>

<mask token>


def Statistics(data, Z_value):

    def running_mean(x, N):
        """ x == an array of data. N == number of samples per average """
        cumsum = np.cumsum(np.insert(x, 0, 0))
        return (cumsum[N:] - cumsum[:-N]) / float(N)
    Sc = []
    for i in range(Num_traj):
        X_mean = running_mean(data[i], int(T_obs / dt))
        prob = 1 / Num_traj * exp(-k * np.mean(data[i]) * Ta) * Z_value
        Sc.append([max(X_mean), prob / run_count])
    return Sc


<mask token>


def F_plot(T, data):
    i = 0
    while i < Num_traj:
        plt.plot(T, data[i])
        i = i + 1


<mask token>


def Data_scorer(data, iter):
    norm = 0
    score = []
    integ = 0
    start = int((iter - 1) * tau / dt)
    stop = int(iter * tau / dt) - 1
    print(start, stop)
    for i in range(Num_traj):
        integ = np.mean(data[i][start:stop]) * tau
        score.append(exp(k * integ))
        norm = norm + exp(k * integ)
    return np.array(score), norm


<mask token>


def GKTL_iter(data, Z, iter):
    new_data = []
    weight, norm = Data_scorer(data, iter)
    mutation_count = Num_traj * weight / norm
    Z = Z * norm / Num_traj
    if iter == Ta / tau:
        return Z
    for i in range(len(mutation_count)):
        mutation_count[i] = mutation_count[i] + random.random()
    mutation_count = np.floor(mutation_count)
    for i in range(Num_traj):
        for j in range(int(mutation_count[i])):
            X_new = data[i].copy()
            m = int(iter * tau / dt)
            stop = int((iter + 1) * tau / dt)
            while m < stop:
                X_new[m] = X_new[m - 1] - alpha * X_new[m - 1] * dt + sqrt(
                    dt * 2 * epsilon) * np.random.normal(0, 1)
                m = m + 1
            new_data.append(X_new)
    S = int(sum(mutation_count))
    if S > Num_traj:
        for c in range(S - Num_traj):
            traj_random = np.random.randint(0, S - c - 1)
            new_data.pop(traj_random)
    if S < Num_traj:
        for c in range(Num_traj - S):
            traj_random = np.random.randint(0, S - 1)
            traj = new_data[traj_random].copy()
            new_data.append(traj)
    return new_data, Z


<mask token>


def Return_time():
    data = []
    for i in range(Num_traj):
        X = np.zeros(N)
        X[0] = np.random.normal(0, sqrt(epsilon / alpha))
        end = int(tau / dt)
        for j in range(end):
            X[j] = X[j - 1] - alpha * X[j - 1] * dt + sqrt(dt * 2 * epsilon
                ) * np.random.normal(0, 1)
        data.append(X)
    iter = 1
    Z = 1
    while iter <= Ta / tau:
        if iter == Ta / tau:
            Z = GKTL_iter(data, Z, iter)
            break
        data, Z = GKTL_iter(data, Z, iter)
        iter = iter + 1
    S = Statistics(data, Z)
    return S


<mask token>
for i in range(run_count):
    S = Return_time()
    for j in S:
        Stats.append(j)
Stats.sort()
<mask token>
while a <= max(Stats[:, 0]):
    p = 0
    for i in range(len(Stats)):
        if Stats[i][0] > a:
            p = p + Stats[i][1]
    A.append(a)
    T.append(-(Ta - T_obs) / log(1 - p))
    a = a + 0.01
plt.plot(np.log10(T), A)
plt.grid(True)
plt.show()

<mask token>
alpha = 1
epsilon = 0.5
Ta = 100
Num_traj = 600
T_obs = 50
dt = 1 / 1000
N = int(Ta / dt)
tau = 0.5
k = 0.3
run_count = 100
<mask token>


def Statistics(data, Z_value):

    def running_mean(x, N):
        """ x == an array of data. N == number of samples per average """
        cumsum = np.cumsum(np.insert(x, 0, 0))
        return (cumsum[N:] - cumsum[:-N]) / float(N)
    Sc = []
    for i in range(Num_traj):
        X_mean = running_mean(data[i], int(T_obs / dt))
        prob = 1 / Num_traj * exp(-k * np.mean(data[i]) * Ta) * Z_value
        Sc.append([max(X_mean), prob / run_count])
    return Sc


<mask token>


def F_plot(T, data):
    i = 0
    while i < Num_traj:
        plt.plot(T, data[i])
        i = i + 1


<mask token>


def Data_scorer(data, iter):
    norm = 0
    score = []
    integ = 0
    start = int((iter - 1) * tau / dt)
    stop = int(iter * tau / dt) - 1
    print(start, stop)
    for i in range(Num_traj):
        integ = np.mean(data[i][start:stop]) * tau
        score.append(exp(k * integ))
        norm = norm + exp(k * integ)
    return np.array(score), norm


<mask token>


def GKTL_iter(data, Z, iter):
    new_data = []
    weight, norm = Data_scorer(data, iter)
    mutation_count = Num_traj * weight / norm
    Z = Z * norm / Num_traj
    if iter == Ta / tau:
        return Z
    for i in range(len(mutation_count)):
        mutation_count[i] = mutation_count[i] + random.random()
    mutation_count = np.floor(mutation_count)
    for i in range(Num_traj):
        for j in range(int(mutation_count[i])):
            X_new = data[i].copy()
            m = int(iter * tau / dt)
            stop = int((iter + 1) * tau / dt)
            while m < stop:
                X_new[m] = X_new[m - 1] - alpha * X_new[m - 1] * dt + sqrt(
                    dt * 2 * epsilon) * np.random.normal(0, 1)
                m = m + 1
            new_data.append(X_new)
    S = int(sum(mutation_count))
    if S > Num_traj:
        for c in range(S - Num_traj):
            traj_random = np.random.randint(0, S - c - 1)
            new_data.pop(traj_random)
    if S < Num_traj:
        for c in range(Num_traj - S):
            traj_random = np.random.randint(0, S - 1)
            traj = new_data[traj_random].copy()
            new_data.append(traj)
    return new_data, Z


<mask token>


def Return_time():
    data = []
    for i in range(Num_traj):
        X = np.zeros(N)
        X[0] = np.random.normal(0, sqrt(epsilon / alpha))
        end = int(tau / dt)
        for j in range(end):
            X[j] = X[j - 1] - alpha * X[j - 1] * dt + sqrt(dt * 2 * epsilon
                ) * np.random.normal(0, 1)
        data.append(X)
    iter = 1
    Z = 1
    while iter <= Ta / tau:
        if iter == Ta / tau:
            Z = GKTL_iter(data, Z, iter)
            break
        data, Z = GKTL_iter(data, Z, iter)
        iter = iter + 1
    S = Statistics(data, Z)
    return S


<mask token>
Stats = []
for i in range(run_count):
    S = Return_time()
    for j in S:
        Stats.append(j)
Stats.sort()
Stats = np.array(Stats)
<mask token>
A = []
T = []
a = 0
while a <= max(Stats[:, 0]):
    p = 0
    for i in range(len(Stats)):
        if Stats[i][0] > a:
            p = p + Stats[i][1]
    A.append(a)
    T.append(-(Ta - T_obs) / log(1 - p))
    a = a + 0.01
plt.plot(np.log10(T), A)
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from math import *
import random
import seaborn as sns
<mask token>
alpha = 1
epsilon = 0.5
Ta = 100
Num_traj = 600
T_obs = 50
dt = 1 / 1000
N = int(Ta / dt)
tau = 0.5
k = 0.3
run_count = 100
<mask token>


def Statistics(data, Z_value):

    def running_mean(x, N):
        """ x == an array of data. N == number of samples per average """
        cumsum = np.cumsum(np.insert(x, 0, 0))
        return (cumsum[N:] - cumsum[:-N]) / float(N)
    Sc = []
    for i in range(Num_traj):
        X_mean = running_mean(data[i], int(T_obs / dt))
        prob = 1 / Num_traj * exp(-k * np.mean(data[i]) * Ta) * Z_value
        Sc.append([max(X_mean), prob / run_count])
    return Sc


<mask token>


def F_plot(T, data):
    i = 0
    while i < Num_traj:
        plt.plot(T, data[i])
        i = i + 1


<mask token>


def Data_scorer(data, iter):
    norm = 0
    score = []
    integ = 0
    start = int((iter - 1) * tau / dt)
    stop = int(iter * tau / dt) - 1
    print(start, stop)
    for i in range(Num_traj):
        integ = np.mean(data[i][start:stop]) * tau
        score.append(exp(k * integ))
        norm = norm + exp(k * integ)
    return np.array(score), norm


<mask token>


def GKTL_iter(data, Z, iter):
    new_data = []
    weight, norm = Data_scorer(data, iter)
    mutation_count = Num_traj * weight / norm
    Z = Z * norm / Num_traj
    if iter == Ta / tau:
        return Z
    for i in range(len(mutation_count)):
        mutation_count[i] = mutation_count[i] + random.random()
    mutation_count = np.floor(mutation_count)
    for i in range(Num_traj):
        for j in range(int(mutation_count[i])):
            X_new = data[i].copy()
            m = int(iter * tau / dt)
            stop = int((iter + 1) * tau / dt)
            while m < stop:
                X_new[m] = X_new[m - 1] - alpha * X_new[m - 1] * dt + sqrt(
                    dt * 2 * epsilon) * np.random.normal(0, 1)
                m = m + 1
            new_data.append(X_new)
    S = int(sum(mutation_count))
    if S > Num_traj:
        for c in range(S - Num_traj):
            traj_random = np.random.randint(0, S - c - 1)
            new_data.pop(traj_random)
    if S < Num_traj:
        for c in range(Num_traj - S):
            traj_random = np.random.randint(0, S - 1)
            traj = new_data[traj_random].copy()
            new_data.append(traj)
    return new_data, Z


<mask token>


def Return_time():
    data = []
    for i in range(Num_traj):
        X = np.zeros(N)
        X[0] = np.random.normal(0, sqrt(epsilon / alpha))
        end = int(tau / dt)
        for j in range(end):
            X[j] = X[j - 1] - alpha * X[j - 1] * dt + sqrt(dt * 2 * epsilon
                ) * np.random.normal(0, 1)
        data.append(X)
    iter = 1
    Z = 1
    while iter <= Ta / tau:
        if iter == Ta / tau:
            Z = GKTL_iter(data, Z, iter)
            break
        data, Z = GKTL_iter(data, Z, iter)
        iter = iter + 1
    S = Statistics(data, Z)
    return S


<mask token>
Stats = []
for i in range(run_count):
    S = Return_time()
    for j in S:
        Stats.append(j)
Stats.sort()
Stats = np.array(Stats)
<mask token>
A = []
T = []
a = 0
while a <= max(Stats[:, 0]):
    p = 0
    for i in range(len(Stats)):
        if Stats[i][0] > a:
            p = p + Stats[i][1]
    A.append(a)
    T.append(-(Ta - T_obs) / log(1 - p))
    a = a + 0.01
plt.plot(np.log10(T), A)
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from math import *
import random
import seaborn as sns


"""
Setting the parameters for the process and the algorithm.
aplha and epsilon are the parameters of the OU process.
Ta is the simulation time of the trajectories. T is the averaged time of the observable.
Num_traj is the number of trajectories in the ensemble. dt is the timestep of the simulation.
tau is the cloning time of the GKTL algorithm. k is sampling parameter used.
run_count is the number of independent algorithm runs.
"""
alpha=1
epsilon= 0.5
Ta=100
Num_traj=600
T_obs=50
dt=1/1000
N=int(Ta/dt)
tau=0.5
k=0.3
run_count=100


"""
This function calculates the statistics of the ensemble generated by the algorithm.
Returns the score of the trajectories in the ensemble and the probability associated.
"""
def Statistics(data,Z_value):
    def running_mean(x, N):
        """ x == an array of data. N == number of samples per average """
        cumsum = np.cumsum(np.insert(x, 0, 0))
        return (cumsum[N:] - cumsum[:-N]) / float(N)

    Sc=[]
    for i in range(Num_traj):
        X_mean=running_mean(data[i],int(T_obs/dt))
        prob= (1/Num_traj)*exp(-k*np.mean(data[i])*Ta)*Z_value
        Sc.append([max(X_mean),prob/run_count])

    return(Sc)


"""
Function used to plot the trajectories if needed.
"""
def F_plot(T,data):
    i=0
    while i<Num_traj:
        plt.plot(T,data[i])
        i=i+1


"""
This function scores the trajectories in the ensemble.
Called by the algorithm iteration function to calculate the weights of the trajectories in the ensemble.
"""
def Data_scorer(data,iter):

    norm=0
    score=[]
    integ=0
    start=int((iter-1)*tau/dt)
    stop=int(iter*tau/dt)-1
    print(start,stop)
    for i in range(Num_traj):
        integ = np.mean(data[i][start:stop])*tau
        score.append(exp(k*integ))
        norm=norm+exp(k*integ)

    return(np.array(score),norm)


"""
This function is the algorithm iteration.
Calls the Data_scorer function to calculate the weights. Clones the existing trajectories, proportional to their weight.
Returns the new ensemble.
"""
def GKTL_iter(data,Z,iter):

    new_data=[]
    weight,norm=Data_scorer(data,iter)
    mutation_count=Num_traj*weight/norm
    Z=Z*norm/Num_traj

    if iter==Ta/tau:
        return(Z)

    for i in range(len(mutation_count)):
        mutation_count[i]=mutation_count[i]+random.random()

    mutation_count=np.floor(mutation_count)

    for i in range(Num_traj):
        for j in range(int(mutation_count[i])):
            X_new=data[i].copy()
            m=int((iter)*tau/dt)
            stop=int((iter+1)*tau/dt)
            while m<stop:
                X_new[m] = X_new[m-1] - alpha*X_new[m-1]*dt + sqrt(dt*2*epsilon)*np.random.normal(0,1)
                m=m+1
            new_data.append(X_new)

    S=int(sum(mutation_count))
    if S>Num_traj:
        for c in range(S-Num_traj):
            traj_random=np.random.randint(0,S-c-1)
            new_data.pop(traj_random)

    if S<Num_traj:
        for c in range(Num_traj-S):
            traj_random=np.random.randint(0,S-1)
            traj=new_data[traj_random].copy()
            new_data.append(traj)

    return(new_data,Z)


"""
This function does one run of the entire algorithm and returns the run statistics.
Creates the initial trajectories and iterates the algorithm till time Ta.
"""
def Return_time():

    data=[]
    for i in range(Num_traj):
        X = np.zeros(N)
        X[0]=np.random.normal(0,sqrt(epsilon/alpha))

        end=int(tau/dt)
        for j in range(end):
            X[j] = X[j-1] - alpha*X[j-1]*dt + sqrt(dt*2*epsilon)*np.random.normal(0,1)

        data.append(X)

    iter=1
    Z=1
    while iter<=Ta/tau:
        if iter==Ta/tau:
            Z=GKTL_iter(data,Z,iter)
            break
        data,Z = GKTL_iter(data,Z,iter)
        iter=iter+1

    S=Statistics(data,Z)
    return(S)


"""
Obtaining stats over the different runs of the algorithm.
"""
Stats=[]
for i in range(run_count):
    S=Return_time()
    for j in S:
        Stats.append(j)

Stats.sort()
Stats=np.array(Stats)

"""
Calculating return time for different values of a from the Stats.
"""
A=[]
T=[]
a=0
while a<=max(Stats[:,0]):
    p=0
    for i in range(len(Stats)):
        if Stats[i][0]>a:
            p=p+Stats[i][1]

    A.append(a)
    T.append(-(Ta-T_obs)/log(1-p))
    a=a+0.01


#np.save('a_alg.npy', A)
#np.save('t_alg.npy', T)

plt.plot(np.log10(T),A)
#plt.xlim([0,14])
#plt.ylim([0,1])
plt.grid(True)
plt.show()
