text:
using LinearAlgebra
using Polynomials
using Optim
using Distributions
using LsqFit
using Plots
using CSV
using DifferentialEquations
using DelimitedFiles
using StatsPlots
using Measures
using LaTeXStrings
using PyPlot
using GR
using DataFrames
using Dates
using StatsBase
using Pkg
#Pkg.update()

dataFreq_t=readdlm("../data/freq/Data_cleansed/TransnetBW/2018.csv")
dataFreq=dataFreq_t[:,2]
for i in 1:length(dataFreq)
    dataFreq[i]=parse(Float64,dataFreq_t[i,2][10:length(dataFreq_t[i,2])])
end
fs=(dataFreq[:].-50).*2*pi
turnOn=20
posOfdbcrossing=20
solu=Array{Float64,2}(undef,900,2)
sol4=Array{Float64,2}(undef,900,5)
sol5=Array{Float64,2}(undef,900,2)
sol5t=Array{Float64,2}(undef,900,3)

############# read results of all 3 models in #############

dataf3=readdir("./model_atPeak/RES_model_atPeak/")[28:36]
mleres3=Array{Float64,2}(undef,parse(Int64,dataf3[length(dataf3)][15:18]),16)  #8760
for i in 1:length(dataf3)
    le=parse(Int64,dataf3[i][10:13])
    ue=parse(Int64,dataf3[i][15:18])
    mleres3[le:ue,:]=readdlm("./model_atPeak/RES_model_atPeak/"*dataf3[i])[le:ue,:]
end
mleres3

dataf_t=readdir("./model_atDB/RES_model_atDB/")
dataf2=dataf_t[2:length(dataf_t)]
mleres2=Array{Float64,2}(undef,length(dataf2),16)
i=1
for i in 1:length(dataf2)
    mleres2[parse(Int64,dataf2[i][1:length(dataf2[i])-4]),:]=readdlm("./model_atDB/RES_model_atDB/"*dataf2[i])
end
mleres2
Plots.plot(1:8760,mleres2[1:8760,2])

dataf_t1=readdir("./model_always/RES_model_always/")
dataf1=dataf_t1[2:length(dataf_t1)-1]
mleres1=Array{Float64,2}(undef,length(dataf1),16)
i=1
for i in 1:length(dataf1)
    mleres1[parse(Int64,dataf1[i][1:length(dataf1[i])-4]),:]=readdlm("./model_always/RES_model_always/"*dataf1[i])
end
mleres1
Plots.plot(1:4500,mleres1[1:4500,2])

function dyn_moments12!(dy,y,par,t)  #par will look like: [gamma1,gamma2,q,r,D]
    dy[1] = par[3]+(par[4]*t)-par[1]*y[1]-par[2]*y[2]
    dy[2] = y[1]
    dy[3] = par[5]^2-(2*par[1]*y[3])-(2*par[2]*y[5])
    dy[4] = 2*y[5]
    dy[5] = -par[1]*y[5]-par[2]*y[4]+y[3]
end
function dyn_moments1_22!(dy,y,par,t)  #par will look like: [gamma1,gamma2,q,r,D]
    dy[1] = par[3]+(par[4]*t)-par[1]*((t-posOfdbcrossing+1)/turnOn)*y[1]
    dy[2] = y[1]
    dy[3] = par[5]^2-(2*par[1]*((t-posOfdbcrossing+1)/turnOn)*y[3])
    dy[4] = 2*y[5]
    dy[5] = -par[1]*((t-posOfdbcrossing+1)/turnOn)*y[5]+y[3]
end
function dyn_moments02!(dy,y,par,t)  #par will look like: [gamma1,gamma2,q,r,D]
    dy[1] = par[3]+(par[4]*t)
    dy[2] = y[1]
    dy[3] = par[5]^2
    dy[4] = 2*y[5]
    dy[5] = y[3]
end
function timeseries_fpes(posOfdbcrossing,turnOn,par) #with smooth turn on but no theta-values
    if posOfdbcrossing==0       #always in db-->control off
        so=solve(ODEProblem(dyn_moments02!,[par[6],par[7],par[8],par[9],par[10]],(0.0,899.0),par[1:5]))
        for j in 1:900
            global sol4[j,:]=so(j-1.0)
        end
    elseif posOfdbcrossing<0       #always in db-->control off
        so=solve(ODEProblem(dyn_moments12!,[par[6],par[7],par[8],par[9],par[10]],(0.0,899.0),par[1:5]))
        for j in 1:900
            global sol4[j,:]=so(j-1.0)
        end
    else                             # posOfdbcrossing>0 (first inside db, then outside)
        so=solve(ODEProblem(dyn_moments02!,[par[6],par[7],par[8],par[9],par[10]],(0.0,posOfdbcrossing-1),par[1:5]))
        for j in 1:posOfdbcrossing
            global sol4[j,:]=so(j-1.0)
        end
        so=solve(ODEProblem(dyn_moments1_22!,[sol4[posOfdbcrossing,1],sol4[posOfdbcrossing,2],sol4[posOfdbcrossing,3],sol4[posOfdbcrossing,4],sol4[posOfdbcrossing,5]],(posOfdbcrossing-1.0,posOfdbcrossing+(turnOn-1.0)),par[1:5]))
        for j in posOfdbcrossing:(posOfdbcrossing+turnOn)
            global sol4[j,:]=so(j-1.0)
        end
        so=solve(ODEProblem(dyn_moments12!,[sol4[(posOfdbcrossing+turnOn),1],sol4[(posOfdbcrossing+turnOn),2],sol4[(posOfdbcrossing+turnOn),3],sol4[(posOfdbcrossing+turnOn),4],sol4[(posOfdbcrossing+turnOn),5]],(posOfdbcrossing+(turnOn-1.0),899.0),par[1:5]))
        for j in (posOfdbcrossing+turnOn):900
            global sol4[j,:]=so(j-1.0)
        end
    end
    global sol5[:,1]=sol4[:,1]
    global sol5[:,2]=sol4[:,3]
    return sol5
end
function timeseries_fpest(posOfdbcrossing,turnOn,par) #with smooth turn on and theta-values
    if posOfdbcrossing==0       #always in db-->control off
        so=solve(ODEProblem(dyn_moments02!,[par[6],par[7],par[8],par[9],par[10]],(0.0,899.0),par[1:5]))
        for j in 1:900
            global sol4[j,:]=so(j-1.0)
        end
    elseif posOfdbcrossing<0       #always in db-->control off
        so=solve(ODEProblem(dyn_moments12!,[par[6],par[7],par[8],par[9],par[10]],(0.0,899.0),par[1:5]))
        for j in 1:900
            global sol4[j,:]=so(j-1.0)
        end
    else                             # posOfdbcrossing>0 (first inside db, then outside)
        so=solve(ODEProblem(dyn_moments02!,[par[6],par[7],par[8],par[9],par[10]],(0.0,posOfdbcrossing-1),par[1:5]))
        for j in 1:posOfdbcrossing
            global sol4[j,:]=so(j-1.0)
        end
        so=solve(ODEProblem(dyn_moments1_22!,[sol4[posOfdbcrossing,1],sol4[posOfdbcrossing,2],sol4[posOfdbcrossing,3],sol4[posOfdbcrossing,4],sol4[posOfdbcrossing,5]],(posOfdbcrossing-1.0,posOfdbcrossing+(turnOn-1.0)),par[1:5]))
        for j in posOfdbcrossing:(posOfdbcrossing+turnOn)
            global sol4[j,:]=so(j-1.0)
        end
        so=solve(ODEProblem(dyn_moments12!,[sol4[(posOfdbcrossing+turnOn),1],sol4[(posOfdbcrossing+turnOn),2],sol4[(posOfdbcrossing+turnOn),3],sol4[(posOfdbcrossing+turnOn),4],sol4[(posOfdbcrossing+turnOn),5]],(posOfdbcrossing+(turnOn-1.0),899.0),par[1:5]))
        for j in (posOfdbcrossing+turnOn):900
            global sol4[j,:]=so(j-1.0)
        end
    end
    global sol5t[:,1]=sol4[:,1]
    global sol5t[:,2]=sol4[:,3]
    global sol5t[:,3]=sol4[:,2]
    return sol5t
end
function mu_sigma2_t(x,para::AbstractArray)
    solu[:,1]=(para[5]-para[3]/(para[1]-para[2]))*exp.(-x.*para[1]).+para[3]/(para[1]-para[2])*exp.(-x.*para[2])
    solu[:,2]=(para[6]^2-para[4]/(2*para[1]))*exp.(x.*-2*para[1]).+para[4]/(2*para[1])
    return solu
end
function sort_out1(mleres,fs) #gives a list of all intervals where mle worked (no NaNs) and where mle makes sense
    #      1) L small ---> var negative
    j=1
    Lsmall=zeros(length(mleres[:,1]))
    for i in 1:length(mleres[:,1])
        if mleres[i,16]<-1e13 # since all intervaks w/o control have L=-1e-25
            global Lsmall[j]=i
            j+=1
        end
    end
    Lsmall2=Lsmall[1:j-1] # intervals with negative vav
    #      2) comparison of MLE to smoothing model
    #            use AIC Value: 2*k-2ln(L),  with k=number of estimated parameters and ln(L): Log-Likelihood
    #                k_mle= 10 from mle + turn-on-procedure, which makes 1 more??
    #                k_movAv= var
    smoothfs=Array{Float64,1}(undef,900)
    smoothfs2=Array{Float64,1}(undef,900)
    vari=Array{Float64,1}(undef,length(mleres[:,1]))
    likeli=Array{Float64,1}(undef,length(mleres[:,1]))
    likelidiff=Array{Float64,1}(undef,length(mleres[:,1]))
    #AIC_diff=Array{Float64,1}(undef,length(mleres[:,1]))
    for h in 1:length(mleres[:,1])
        for i in 1:70
            global smoothfs[i]=sum(fs[(((h-1)*3600)+1):(((h-1)*3600)+(2*i))])/(2*i)
        end
        for i in 71:830
            global smoothfs[i]=sum(fs[(((h-1)*3600)+(i-70)):(((h-1)*3600)+(i+70))])/140
        end
        for i in 831:900
            global smoothfs[i]=sum(fs[(((h-1)*3600)+900-(2*(900-i))):(((h-1)*3600)+900)])/(2*(900-i)+1)
        end
        for i in 880:900
            global smoothfs[i]=smoothfs[879]
        end
        for i in 1:40
            global smoothfs2[i]=sum(smoothfs[1:(2*i)])/(2*i)
        end
        for i in 41:860
            global smoothfs2[i]=sum(smoothfs[(i-40):(i+40)])/80
        end
        for i in 861:880
            global smoothfs2[i]=sum(smoothfs[(900-(2*(900-i))):900])/(2*(900-i))
        end
        for i in 880:900
            global smoothfs2[i]=smoothfs2[879]
        end
        global vari[h]=var(fs[((h-1)*3600+1):((h-1)*3600+900)].-smoothfs2)
        global likeli[h]= -sum((1/2)*log(2*pi*vari[h]).-(((fs[(3600*(h-1)+1):(3600*(h-1)+900)].-smoothfs2).^2)/(2*vari[h])))
        likelidiff[h]=abs(likeli[h]-mleres[h,16])
        #AIC_diff[h]=abs(2*10-2*mleres[h,16]-(2*1-2*likeli[h]))
        if likelidiff[h]>1e4
            likelidiff[h]=NaN
        end
        #if AIC_diff[h]>2e4
        #    AIC_diff[h]=NaN
        #end
    end
    #plot(1:length(mleres[:,1]),likelidiff[1:length(mleres[:,1])],seriestype = :scatter,ylims=(700,2000))
    #boxplot(filter(!isnan,likelidiff[1:length(mleres[:,1])]),ylims=(700,2000))
    #plot(1:length(mleres[:,1]),AIC_diff[1:length(mleres[:,1])],seriestype = :scatter,ylims=(700,6000))
    #boxplot(filter(!isnan,AIC_diff[1:length(mleres[:,1])]),ylims=(700,6000))
    q10_L=quantile(filter(!isnan,likelidiff[1:length(mleres[:,1])]),0.9)
    #q10_AIC=quantile(filter(!isnan,AIC_diff[1:length(mleres[:,1])]),0.9)
    j=1
    #k=1
    likelidifflarge=Array{Float64,1}(undef,length(mleres[:,1]))
    #AICdifflarge=Array{Float64,1}(undef,length(mleres[:,1]))
    for i in 1:length(mleres[:,1])
        if abs(likelidiff[i])>q10_L || (isnan(likelidiff[i])==true)
            likelidifflarge[j]=i
            j+=1
        end
        #if abs(AIC_diff[i])>q10_AIC
        #    AIC_diff[k]=i
        #    k+=1
        #end
    end
    likelidifflarge2=likelidifflarge[1:j-1]
    #AIC_diff2=AIC_diff[1:j-1]
    #AIC_diff2==likelidifflarge2
    #    bring likelilarge and neg var together:
    goodones2=Array{Float64,1}(undef,length(mleres[:,1]))
    throuwout2=Array{Float64,1}(undef,length(Lsmall2)+length(likelidifflarge2))
    throuwout2[1:length(Lsmall2)]=Lsmall2
    throuwout2[length(Lsmall2)+1:(length(Lsmall2)+length(likelidifflarge2))]=likelidifflarge2
    j=1
    for i in 1:length(mleres[:,1])
        if i in throuwout2
        else
            global goodones2[j]=i
            j+=1
        end
    end
    goodones2
    return goods1=convert.(Int64,goodones2[1:j-1]),likelidiff
end
function sort_out2(mleres,fs) #gives a list of all intervals where mle worked (no NaNs) and where mle makes sense
    #      1) L small ---> var negative
    j=1
    Lsmall=zeros(length(mleres[:,1]))
    for i in 1:length(mleres[:,1])
        if mleres[i,16]<-1e13 # since all intervaks w/o control have L=-1e-25
            global Lsmall[j]=i
            j+=1
        end
    end
    Lsmall2=Lsmall[1:j-1] # intervals with negative vav
    #      2) comparison of MLE to smoothing model
    #            use AIC Value: 2*k-2ln(L),  with k=number of estimated parameters and ln(L): Log-Likelihood
    #                k_mle= 10 from mle + turn-on-procedure, which makes 1 more??
    #                k_movAv= var
    smoothfs=Array{Float64,1}(undef,900)
    smoothfs2=Array{Float64,1}(undef,900)
    vari=Array{Float64,1}(undef,length(mleres[:,1]))
    likeli=Array{Float64,1}(undef,length(mleres[:,1]))
    likelidiff=Array{Float64,1}(undef,length(mleres[:,1]))
    #AIC_diff=Array{Float64,1}(undef,length(mleres[:,1]))
    for h in 1:length(mleres[:,1])
        for i in 1:70
            global smoothfs[i]=sum(fs[(((h-1)*3600)+1):(((h-1)*3600)+(2*i))])/(2*i)
        end
        for i in 71:830
            global smoothfs[i]=sum(fs[(((h-1)*3600)+(i-70)):(((h-1)*3600)+(i+70))])/140
        end
        for i in 831:900
            global smoothfs[i]=sum(fs[(((h-1)*3600)+900-(2*(900-i))):(((h-1)*3600)+900)])/(2*(900-i)+1)
        end
        for i in 880:900
            global smoothfs[i]=smoothfs[879]
        end
        for i in 1:40
            global smoothfs2[i]=sum(smoothfs[1:(2*i)])/(2*i)
        end
        for i in 41:860
            global smoothfs2[i]=sum(smoothfs[(i-40):(i+40)])/80
        end
        for i in 861:880
            global smoothfs2[i]=sum(smoothfs[(900-(2*(900-i))):900])/(2*(900-i))
        end
        for i in 880:900
            global smoothfs2[i]=smoothfs2[879]
        end
        global vari[h]=var(fs[((h-1)*3600+1):((h-1)*3600+900)].-smoothfs2)
        global likeli[h]= -sum((1/2)*log(2*pi*vari[h]).-(((fs[(3600*(h-1)+1):(3600*(h-1)+900)].-smoothfs2).^2)/(2*vari[h])))
        likelidiff[h]=abs(likeli[h]-mleres[h,16])
        #AIC_diff[h]=abs(2*10-2*mleres[h,16]-(2*1-2*likeli[h]))
        if likelidiff[h]>1e4
            likelidiff[h]=NaN
        end
        #if AIC_diff[h]>2e4
        #    AIC_diff[h]=NaN
        #end
    end
    #plot(1:length(mleres[:,1]),likelidiff[1:length(mleres[:,1])],seriestype = :scatter,ylims=(700,2000))
    #boxplot(filter(!isnan,likelidiff[1:length(mleres[:,1])]),ylims=(700,2000))
    #plot(1:length(mleres[:,1]),AIC_diff[1:length(mleres[:,1])],seriestype = :scatter,ylims=(700,6000))
    #boxplot(filter(!isnan,AIC_diff[1:length(mleres[:,1])]),ylims=(700,6000))
    q10_L=quantile(filter(!isnan,likelidiff[1:length(mleres[:,1])]),0.9)
    #q10_AIC=quantile(filter(!isnan,AIC_diff[1:length(mleres[:,1])]),0.9)
    j=1
    #k=1
    likelidifflarge=Array{Float64,1}(undef,length(mleres[:,1]))
    #AICdifflarge=Array{Float64,1}(undef,length(mleres[:,1]))
    for i in 1:length(mleres[:,1])
        if abs(likelidiff[i])>q10_L || (isnan(likelidiff[i])==true)
            likelidifflarge[j]=i
            j+=1
        end
        #if abs(AIC_diff[i])>q10_AIC
        #    AIC_diff[k]=i
        #    k+=1
        #end
    end
    likelidifflarge2=likelidifflarge[1:j-1]
    #AIC_diff2=AIC_diff[1:j-1]
    #AIC_diff2==likelidifflarge2
    #    bring likelilarge and neg var together:
    goodones2=Array{Float64,1}(undef,length(mleres[:,1]))
    throuwout2=Array{Float64,1}(undef,length(Lsmall2)+length(likelidifflarge2))
    throuwout2[1:length(Lsmall2)]=Lsmall2
    throuwout2[length(Lsmall2)+1:(length(Lsmall2)+length(likelidifflarge2))]=likelidifflarge2
    j=1
    for i in 1:length(mleres[:,1])
        if i in throuwout2
        else
            global goodones2[j]=i
            j+=1
        end
    end
    goodones2
    return goods2=convert.(Int64,goodones2[1:j-1]),likelidiff
end
function sort_out3(mleres,fs) #gives a list of all intervals where mle worked (no NaNs) and where mle makes sense
    #      1) L small ---> var negative
    j=1
    Lsmall=zeros(length(mleres[:,1]))
    for i in 1:length(mleres[:,1])
        if mleres[i,16]<-1e13 # since all intervaks w/o control have L=-1e-25
            global Lsmall[j]=i
            j+=1
        end
    end
    Lsmall2=Lsmall[1:j-1] # intervals with negative vav
    #      2) comparison of MLE to smoothing model
    #            use AIC Value: 2*k-2ln(L),  with k=number of estimated parameters and ln(L): Log-Likelihood
    #                k_mle= 10 from mle + turn-on-procedure, which makes 1 more??
    #                k_movAv= var
    smoothfs=Array{Float64,1}(undef,900)
    smoothfs2=Array{Float64,1}(undef,900)
    vari=Array{Float64,1}(undef,length(mleres[:,1]))
    likeli=Array{Float64,1}(undef,length(mleres[:,1]))
    likelidiff=Array{Float64,1}(undef,length(mleres[:,1]))
    #AIC_diff=Array{Float64,1}(undef,length(mleres[:,1]))
    for h in 1:length(mleres[:,1])
        for i in 1:70
            global smoothfs[i]=sum(fs[(((h-1)*3600)+1):(((h-1)*3600)+(2*i))])/(2*i)
        end
        for i in 71:830
            global smoothfs[i]=sum(fs[(((h-1)*3600)+(i-70)):(((h-1)*3600)+(i+70))])/140
        end
        for i in 831:900
            global smoothfs[i]=sum(fs[(((h-1)*3600)+900-(2*(900-i))):(((h-1)*3600)+900)])/(2*(900-i)+1)
        end
        for i in 880:900
            global smoothfs[i]=smoothfs[879]
        end
        for i in 1:40
            global smoothfs2[i]=sum(smoothfs[1:(2*i)])/(2*i)
        end
        for i in 41:860
            global smoothfs2[i]=sum(smoothfs[(i-40):(i+40)])/80
        end
        for i in 861:880
            global smoothfs2[i]=sum(smoothfs[(900-(2*(900-i))):900])/(2*(900-i))
        end
        for i in 880:900
            global smoothfs2[i]=smoothfs2[879]
        end
        global vari[h]=var(fs[((h-1)*3600+1):((h-1)*3600+900)].-smoothfs2)
        global likeli[h]= -sum((1/2)*log(2*pi*vari[h]).-(((fs[(3600*(h-1)+1):(3600*(h-1)+900)].-smoothfs2).^2)/(2*vari[h])))
        likelidiff[h]=abs(likeli[h]-mleres[h,16])
        #AIC_diff[h]=abs(2*10-2*mleres[h,16]-(2*1-2*likeli[h]))
        if likelidiff[h]>1e4
            likelidiff[h]=NaN
        end
        #if AIC_diff[h]>2e4
        #    AIC_diff[h]=NaN
        #end
    end
    #plot(1:length(mleres[:,1]),likelidiff[1:length(mleres[:,1])],seriestype = :scatter,ylims=(700,2000))
    #boxplot(filter(!isnan,likelidiff[1:length(mleres[:,1])]),ylims=(700,2000))
    #plot(1:length(mleres[:,1]),AIC_diff[1:length(mleres[:,1])],seriestype = :scatter,ylims=(700,6000))
    #boxplot(filter(!isnan,AIC_diff[1:length(mleres[:,1])]),ylims=(700,6000))
    q10_L=quantile(filter(!isnan,likelidiff[1:length(mleres[:,1])]),0.9)
    #q10_AIC=quantile(filter(!isnan,AIC_diff[1:length(mleres[:,1])]),0.9)
    j=1
    #k=1
    likelidifflarge=Array{Float64,1}(undef,length(mleres[:,1]))
    #AICdifflarge=Array{Float64,1}(undef,length(mleres[:,1]))
    for i in 1:length(mleres[:,1])
        if abs(likelidiff[i])>q10_L || (isnan(likelidiff[i])==true)
            likelidifflarge[j]=i
            j+=1
        end
        #if abs(AIC_diff[i])>q10_AIC
        #    AIC_diff[k]=i
        #    k+=1
        #end
    end
    likelidifflarge2=likelidifflarge[1:j-1]
    #AIC_diff2=AIC_diff[1:j-1]
    #AIC_diff2==likelidifflarge2
    #    bring likelilarge and neg var together:
    goodones2=Array{Float64,1}(undef,length(mleres[:,1]))
    throuwout2=Array{Float64,1}(undef,length(Lsmall2)+length(likelidifflarge2))
    throuwout2[1:length(Lsmall2)]=Lsmall2
    throuwout2[length(Lsmall2)+1:(length(Lsmall2)+length(likelidifflarge2))]=likelidifflarge2
    j=1
    for i in 1:length(mleres[:,1])
        if i in throuwout2
        else
            global goodones2[j]=i
            j+=1
        end
    end
    goodones2
    return goods3=convert.(Int64,goodones2[1:j-1]),likelidiff
end
function density_plot(xli,yli,xdat,ydat,zdat,gridx,gridy)#(xl,xu)=xlims(plo)
    xd=(xli[2]-xli[1])/gridx
    yd=(yli[2]-yli[1])/gridy
    xval=zeros(length(xdat))
    yval=zeros(length(xdat))
    zval=zeros(length(xdat))
    c=1
    for i in 1:gridx
        for j in 1:gridy
            counte=0
            ztemp=0
            for k in 1:length(xdat)
                if (xli[1]+xd*(i-1))<xdat[k]<(xli[1]+xd*i)
                    if (yli[1]+yd*(j-1))<ydat[k]<(yli[1]+yd*j)
                        ztemp+=zdat[k]
                        counte+=1
                    end
                end
            end
            ztemp=ztemp/counte
            if isnan(ztemp)==false
                xval[c]=(xli[1]+xd*(i-0.5))
                yval[c]=(yli[1]+yd*(j-0.5))
                zval[c]=ztemp
                c+=1
            end
        end
    end
    return (xval[1:c-1],yval[1:c-1],zval[1:c-1])
end
function plotInt1(hour)
    global turnOn=convert(Int64,mleres1[hour,1])
    global posOfdbcrossing=convert(Int64,mleres1[hour,2])
    mleres1[hour,8]
    TSdp=Array{Float64,1}(undef,900)
    TSdpt=Array{Float64,1}(undef,900)
    #TSdp[1:100].=NaN
    #TSdpt[1:100].=NaN
    TS=timeseries_fpest(convert(Int64,mleres1[hour,2]),convert(Int64,mleres1[hour,1]),mleres1[hour,3:12])
    for i in 1:900
         TSdp[i]=mleres1[hour,5]+mleres1[hour,6]*(i-1)
         if i<posOfdbcrossing+turnOn
             TSdpt[i]=TSdp[i]
         else
             TSdpt[i]=TSdp[i]-mleres1[hour,4]*TS[i,3]
         end
    end
    TSsigma = sqrt.(abs.(TS[:,2]))
    TSmuPlus = TS[:,1]+TSsigma[:]
    TSmuMinus = TS[:,1]-TSsigma[:]
    yupper=maximum([maximum(TSdp[1:900]),maximum(TSdpt[1:900]),0])
    ylower=minimum([minimum(TSdp[1:900]),minimum(TSdpt[1:900]),0])
    inform="g1= "*string(round(mleres1[hour,3];digits=5))*"   g2= "*string(round(mleres1[hour,4];digits=5))*"   D= "*string(round(mleres1[hour,7];digits=3))*"   r= "*string(round(mleres1[hour,6];digits=7))*"   L= "*string(round(mleres1[hour,16];digits=0))*"   AIC= "*string(round((2*10-2*mleres1[hour,16]);digits=0))
    ueberschrift="\$\\textrm{model \\ 1, \\ hour \\ }"*string(hour)*"\$"
    #path=".\plots\first70\combi"*string(hour)
    p1=Plots.plot(0:899,[TS[:,1] TSmuPlus TSmuMinus],color=[:blue :lightblue :lightblue],style=[:solid :dash :dash],label=false,ylabel=L"\omega",xaxis=false,title=ueberschrift)
    Plots.plot!(0:899,fs[((hour-1)*3600+1):((hour-1)*3600+900)],label=false,color=:black)
    p2=Plots.plot(0:899,[TSdp TSdpt],color=[:blue :lightblue],label=[L"\Delta \textrm{P}" L"\Delta \textrm{P}+2\textrm{nd}"],ylims=(ylower,yupper),xlabel=L"\textrm{time \ in \ s}") #,annotations = (0,yupper,Plots.text(inform, :left,10))
    Plots.plot(p1,p2,layout = (2, 1))
end
function plotInt2(hour)
    global turnOn=convert(Int64,mleres2[hour,1])
    global posOfdbcrossing=convert(Int64,mleres2[hour,2])
    mleres2[hour,8]
    TSdp=Array{Float64,1}(undef,900)
    TSdpt=Array{Float64,1}(undef,900)
    #TSdp[1:100].=NaN
    #TSdpt[1:100].=NaN
    TS=timeseries_fpest(convert(Int64,mleres2[hour,2]),convert(Int64,mleres2[hour,1]),mleres2[hour,3:12])
    for i in 1:900
         TSdp[i]=mleres2[hour,5]+mleres2[hour,6]*(i-1)
         if i<posOfdbcrossing+turnOn
             TSdpt[i]=TSdp[i]
         else
             TSdpt[i]=TSdp[i]-mleres2[hour,4]*TS[i,3]
         end
    end
    TSsigma = sqrt.(abs.(TS[:,2]))
    TSmuPlus = TS[:,1]+TSsigma[:]
    TSmuMinus = TS[:,1]-TSsigma[:]
    yupper=maximum([maximum(TSdp[1:900]),maximum(TSdpt[1:900]),0])
    ylower=minimum([minimum(TSdp[1:900]),minimum(TSdpt[1:900]),0])
    inform="g1= "*string(round(mleres2[hour,3];digits=5))*"   g2= "*string(round(mleres2[hour,4];digits=5))*"   D= "*string(round(mleres2[hour,7];digits=3))*"   r= "*string(round(mleres2[hour,6];digits=7))*"   L= "*string(round(mleres2[hour,16];digits=0))*"   AIC= "*string(round((2*10-2*mleres2[hour,16]);digits=0))
    ueberschrift="\$\\textrm{model \\ 2, \\ hour \\ }"*string(hour)*"\$"
    #path=".\plots\first70\combi"*string(hour)
    p1=Plots.plot(0:899,[TS[:,1] TSmuPlus TSmuMinus],color=[:blue :lightblue :lightblue],style=[:solid :dash :dash],label=false,ylabel=L"\omega",xaxis=false,title=ueberschrift)
    Plots.plot!(0:899,fs[((hour-1)*3600+1):((hour-1)*3600+900)],label=false,color=:black)
    p2=Plots.plot(0:899,[TSdp TSdpt],color=[:blue :lightblue],label=[L"\Delta \textrm{P}" L"\Delta \textrm{P}+2\textrm{nd}"],ylims=(ylower,yupper),xlabel=L"\textrm{time \ in \ s}") #,annotations = (0,yupper,Plots.text(inform, :left,10))
    Plots.plot(p1,p2,layout = (2, 1))
end
function plotInt3(hour)
    global turnOn=convert(Int64,mleres3[hour,1])
    global posOfdbcrossing=convert(Int64,mleres3[hour,2])
    mleres3[hour,8]
    TSdp=Array{Float64,1}(undef,1000)
    TSdpt=Array{Float64,1}(undef,1000)
    TSdp[1:100].=NaN
    TSdpt[1:100].=NaN
    TS=timeseries_fpest(convert(Int64,mleres3[hour,2]),convert(Int64,mleres3[hour,1]),mleres3[hour,3:12])
    for i in 1:900
         TSdp[i+100]=mleres3[hour,5]+mleres3[hour,6]*(i-1)
         if i<posOfdbcrossing+turnOn
             TSdpt[i+100]=TSdp[i+100]
         else
             TSdpt[i+100]=TSdp[i+100]-mleres3[hour,4]*TS[i,3]
         end
    end
    TSsigma = sqrt.(abs.(TS[:,2]))
    TSmuPlus = TS[:,1]+TSsigma[:]
    TSmuMinus = TS[:,1]-TSsigma[:]
    yupper=maximum([maximum(TSdp[101:1000]),maximum(TSdpt[101:1000]),0])
    ylower=minimum([minimum(TSdp[101:1000]),minimum(TSdpt[101:1000]),0])
    inform="g1= "*string(round(mleres3[hour,3];digits=5))*"   g2= "*string(round(mleres3[hour,4];digits=5))*"   D= "*string(round(mleres3[hour,7];digits=3))*"   r= "*string(round(mleres3[hour,6];digits=7))*"   L= "*string(round(mleres3[hour,16];digits=0))*"   AIC= "*string(round((2*10-2*mleres3[hour,16]);digits=0))
    ueberschrift="\$\\textrm{model \\ 3, \\ hour \\ }"*string(hour)*"\$"
    #path=".\plots\first70\combi"*string(hour)
    p1=Plots.plot(0:899,[TS[:,1] TSmuPlus TSmuMinus],color=[:blue :lightblue :lightblue],style=[:solid :dash :dash],label=false,ylabel=L"\omega",xaxis=false,title=ueberschrift)
    Plots.plot!(0:899,fs[((hour-1)*3600+1):((hour-1)*3600+900)],label=false,color=:black)
    p2=Plots.plot(-100:899,[TSdp TSdpt],color=[:blue :lightblue],label=[L"\Delta \textrm{P}" L"\Delta \textrm{P}+2\textrm{nd}"],ylims=(ylower,yupper),xlabel=L"\textrm{time \ in \ s}") #,annotations = (0,yupper,Plots.text(inform, :left,10))
    Plots.plot(p1,p2,layout = (2, 1))
end

function plotInt1w(hour)
    global turnOn=convert(Int64,mleres1[hour,1])
    global posOfdbcrossing=convert(Int64,mleres1[hour,2])
    mleres1[hour,8]
    TSdp=Array{Float64,1}(undef,900)
    TSdpt=Array{Float64,1}(undef,900)
    #TSdp[1:100].=NaN
    #TSdpt[1:100].=NaN
    TS=timeseries_fpest(convert(Int64,mleres1[hour,2]),convert(Int64,mleres1[hour,1]),mleres1[hour,3:12])
    for i in 1:900
         TSdp[i]=mleres1[hour,5]+mleres1[hour,6]*(i-1)
         if i<posOfdbcrossing+turnOn
             TSdpt[i]=TSdp[i]
         else
             TSdpt[i]=TSdp[i]-mleres1[hour,4]*TS[i,3]
         end
    end
    TSsigma = sqrt.(abs.(TS[:,2]))
    TSmuPlus = TS[:,1]+TSsigma[:]
    TSmuMinus = TS[:,1]-TSsigma[:]
    yupper=maximum([maximum(TSdp[1:900]),maximum(TSdpt[1:900]),0])
    ylower=minimum([minimum(TSdp[1:900]),minimum(TSdpt[1:900]),0])
    inform="g1= "*string(round(mleres1[hour,3];digits=5))*"   g2= "*string(round(mleres1[hour,4];digits=5))*"   D= "*string(round(mleres1[hour,7];digits=3))*"   r= "*string(round(mleres1[hour,6];digits=7))*"   L= "*string(round(mleres1[hour,16];digits=0))*"   AIC= "*string(round((2*10-2*mleres1[hour,16]);digits=0))
    ueberschrift="\$\\textrm{model \\ 1, \\ hour \\ }"*string(hour)*"\$"
    #path=".\plots\first70\combi"*string(hour)
    p1=Plots.plot(0:899,[TS[:,1] TSmuPlus TSmuMinus],color=[:blue :lightblue :lightblue],style=[:solid :dash :dash],label=false,xlabel=L"\textrm{time \ in \ s}",ylabel=L"\omega",title=ueberschrift)
    Plots.plot!(0:899,fs[((hour-1)*3600+1):((hour-1)*3600+900)],label=false,color=:black)
    #p2=Plots.plot(0:899,[TSdp TSdpt],color=[:blue :lightblue],label=[L"\Delta \textrm{P}" L"\Delta \textrm{P}+2\textrm{nd}"],ylims=(ylower,yupper),xlabel=L"\textrm{time \ in \ s}") #,annotations = (0,yupper,Plots.text(inform, :left,10))
    #Plots.plot(p1,p2,layout = (2, 1))
end
function plotInt2w(hour)
    global turnOn=convert(Int64,mleres2[hour,1])
    global posOfdbcrossing=convert(Int64,mleres2[hour,2])
    mleres2[hour,8]
    TSdp=Array{Float64,1}(undef,900)
    TSdpt=Array{Float64,1}(undef,900)
    #TSdp[1:100].=NaN
    #TSdpt[1:100].=NaN
    TS=timeseries_fpest(convert(Int64,mleres2[hour,2]),convert(Int64,mleres2[hour,1]),mleres2[hour,3:12])
    for i in 1:900
         TSdp[i]=mleres2[hour,5]+mleres2[hour,6]*(i-1)
         if i<posOfdbcrossing+turnOn
             TSdpt[i]=TSdp[i]
         else
             TSdpt[i]=TSdp[i]-mleres2[hour,4]*TS[i,3]
         end
    end
    TSsigma = sqrt.(abs.(TS[:,2]))
    TSmuPlus = TS[:,1]+TSsigma[:]
    TSmuMinus = TS[:,1]-TSsigma[:]
    yupper=maximum([maximum(TSdp[1:900]),maximum(TSdpt[1:900]),0])
    ylower=minimum([minimum(TSdp[1:900]),minimum(TSdpt[1:900]),0])
    inform="g1= "*string(round(mleres2[hour,3];digits=5))*"   g2= "*string(round(mleres2[hour,4];digits=5))*"   D= "*string(round(mleres2[hour,7];digits=3))*"   r= "*string(round(mleres2[hour,6];digits=7))*"   L= "*string(round(mleres2[hour,16];digits=0))*"   AIC= "*string(round((2*10-2*mleres2[hour,16]);digits=0))
    ueberschrift="\$\\textrm{model \\ 2, \\ hour \\ }"*string(hour)*"\$"
    #path=".\plots\first70\combi"*string(hour)
    p1=Plots.plot(0:899,[TS[:,1] TSmuPlus TSmuMinus],color=[:blue :lightblue :lightblue],style=[:solid :dash :dash],label=false,ylabel=L"\omega",xlabel=L"\textrm{time \ in \ s}",title=ueberschrift)
    Plots.plot!(0:899,fs[((hour-1)*3600+1):((hour-1)*3600+900)],label=false,color=:black)
    #p2=Plots.plot(0:899,[TSdp TSdpt],color=[:blue :lightblue],label=[L"\Delta \textrm{P}" L"\Delta \textrm{P}+2\textrm{nd}"],ylims=(ylower,yupper),xlabel=L"\textrm{time \ in \ s}") #,annotations = (0,yupper,Plots.text(inform, :left,10))
    #Plots.plot(p1,p2,layout = (2, 1))
end
function plotInt3w(hour)
    global turnOn=convert(Int64,mleres3[hour,1])
    global posOfdbcrossing=convert(Int64,mleres3[hour,2])
    mleres3[hour,8]
    TSdp=Array{Float64,1}(undef,1000)
    TSdpt=Array{Float64,1}(undef,1000)
    TSdp[1:100].=NaN
    TSdpt[1:100].=NaN
    TS=timeseries_fpest(convert(Int64,mleres3[hour,2]),convert(Int64,mleres3[hour,1]),mleres3[hour,3:12])
    for i in 1:900
         TSdp[i+100]=mleres3[hour,5]+mleres3[hour,6]*(i-1)
         if i<posOfdbcrossing+turnOn
             TSdpt[i+100]=TSdp[i+100]
         else
             TSdpt[i+100]=TSdp[i+100]-mleres3[hour,4]*TS[i,3]
         end
    end
    TSsigma = sqrt.(abs.(TS[:,2]))
    TSmuPlus = TS[:,1]+TSsigma[:]
    TSmuMinus = TS[:,1]-TSsigma[:]
    yupper=maximum([maximum(TSdp[101:1000]),maximum(TSdpt[101:1000]),0])
    ylower=minimum([minimum(TSdp[101:1000]),minimum(TSdpt[101:1000]),0])
    inform="g1= "*string(round(mleres3[hour,3];digits=5))*"   g2= "*string(round(mleres3[hour,4];digits=5))*"   D= "*string(round(mleres3[hour,7];digits=3))*"   r= "*string(round(mleres3[hour,6];digits=7))*"   L= "*string(round(mleres3[hour,16];digits=0))*"   AIC= "*string(round((2*10-2*mleres3[hour,16]);digits=0))
    ueberschrift="\$\\textrm{model \\ 3, \\ hour \\ }"*string(hour)*"\$"
    #path=".\plots\first70\combi"*string(hour)
    p1=Plots.plot(0:899,[TS[:,1] TSmuPlus TSmuMinus],color=[:blue :lightblue :lightblue],style=[:solid :dash :dash],label=false,ylabel=L"\omega",xlabel=L"\textrm{time \ in \ s}",title=ueberschrift)
    Plots.plot!(0:899,fs[((hour-1)*3600+1):((hour-1)*3600+900)],label=false,color=:black)
    #p2=Plots.plot(-100:899,[TSdp TSdpt],color=[:blue :lightblue],label=[L"\Delta \textrm{P}" L"\Delta \textrm{P}+2\textrm{nd}"],ylims=(ylower,yupper),xlabel=L"\textrm{time \ in \ s}") #,annotations = (0,yupper,Plots.text(inform, :left,10))
    #Plots.plot(p1,p2,layout = (2, 1))
end
################## sort out the bad fits ##################

(goods1,likelidiff)=sort_out1(mleres1,fs)
Plots.plot(goods1,likelidiff[goods1],seriestype=:scatter)
goods1

(goods2,likelidiff)=sort_out2(mleres2,fs)
Plots.plot(goods2,likelidiff[goods2],seriestype=:scatter)
goods2

(goods3,likelidiff)=sort_out3(mleres3,fs)
Plots.plot(goods3,likelidiff[goods3],seriestype=:scatter)
goods3

################## make the the histograms ################

goods1
goods2
goods3
goods12=intersect(goods1,goods2)
goods23=intersect(goods2,goods3)
goods13=intersect(goods1,goods3)

likelidiff12=mleres1[goods12,16]-mleres2[goods12,16]
Plots.histogram(-likelidiff12,title="\$\\textrm{Difference \\ of \\ the \\ fit \\ quality}\$",label=false,xlabel="\$\\textrm{AIC}_1-\\textrm{AIC}_2\$")
png("/Users/raphaelbiertz/Documents/masterarbeit/coding/codes/model_comparison/plots/hist12.png")
mean(-likelidiff12)

likelidiff23=mleres2[goods23,16]-mleres3[goods23,16]
Plots.histogram(-likelidiff23,title="\$\\textrm{Difference \\ of \\ the \\ fit \\ quality}\$",label=false,xlabel="\$\\textrm{AIC}_2-\\textrm{AIC}_3\$")
png("/Users/raphaelbiertz/Documents/masterarbeit/coding/codes/model_comparison/plots/hist23.png")
mean(-likelidiff23)

likelidiff13=mleres1[goods13,16]-mleres3[goods13,16]
Plots.histogram(-likelidiff13,title="\$\\textrm{Difference \\ of \\ the \\ fit \\ quality}\$",label=false,xlabel="\$\\textrm{AIC}_1-\\textrm{AIC}_3\$")
png("/Users/raphaelbiertz/Documents/masterarbeit/coding/codes/model_comparison/plots/hist13.png")
mean(-likelidiff13)

################ check, why model 1 and 2 give nearly always the same result ###############
mleres1[:,2]
whereStart_t=zeros(length(mleres2[:,2]))
j=1
for i in 1:8760
    if mleres2[i,2] > -1
        whereStart_t[j]=convert(Int64,i)
        j+=1
    end
end
controlOn=convert.(Int64,whereStart_t[1:j-1])
goods12Con=intersect(goods1,goods2,controlOn)

i=21
plotInt1(goods12Con[i])
plotInt2(goods12Con[i])
mleres2[i,:]

likelidiff12Con=mleres1[goods12Con,16]-mleres2[goods12Con,16]
Plots.histogram(-likelidiff12Con,title="\$\\textrm{Difference \\ of \\ fit \\ quality}\$",label=false,xlabel="\$\\textrm{AIC}_1-\\textrm{AIC}_2\$")
png("/Users/raphaelbiertz/Documents/masterarbeit/coding/codes/model_comparison/plots/hist12Con.png")
mean(-likelidiff12Con)

########### comparison of 2,3 #############

better2=Array{Int64,1}(undef,length(filter!(x -> x > 300, likelidiff23)))
likelidiff23=mleres2[goods23,16]-mleres3[goods23,16]
j=1
for i in 1:length(goods23)
    if likelidiff23[i]>300
        better2[j]=goods23[i]
        j+=1
    end
end
better2

better3=Array{Int64,1}(undef,length(filter!(x -> x < -800, likelidiff23)))
likelidiff23=mleres2[goods23,16]-mleres3[goods23,16]
j=1
for i in 1:length(goods23)
    if likelidiff23[i]<-800
        better3[j]=goods23[i]
        j+=1
    end
end
better3

i=17
plotInt2(better2[i])
plotInt3(better2[i])
beautiful: better2[2,3,4,5,6,7,10,13]
i=20
plotInt2(better3[i])
plotInt3(better3[i])
beautiful: better3[2,4,5,13,17,20   ,109]

i=17
plotInt2w(better2[i])
plotInt3w(better2[i])
beautiful: better2[2,3,4,5,6,7,10,13]
i=17
plotInt2w(better3[i])
plotInt3w(better3[i])
beautiful: better3[2,4,5,13,17,20   ,109]
png("/Users/raphaelbiertz/Documents/masterarbeit/coding/codes/model_comparison/plots/23_m3better2.png")
png("/Users/raphaelbiertz/Documents/masterarbeit/coding/codes/model_comparison/plots/23_m2better2.png")

plotInt1w(226)
plotInt3w(226)
mleres3[2]
########### comparison of 1,2 #############

better1=Array{Int64,1}(undef,length(filter!(x -> x > 300, likelidiff12)))
likelidiff12=mleres1[goods12,16]-mleres2[goods12,16]
j=1
for i in 1:length(goods12)
    if likelidiff12[i]>300
        better1[j]=goods12[i]
        j+=1
    end
end
better2

better3=Array{Int64,1}(undef,length(filter!(x -> x < -800, likelidiff23)))
likelidiff23=mleres2[goods23,16]-mleres3[goods23,16]
j=1
for i in 1:length(goods23)
    if likelidiff23[i]<-800
        better3[j]=goods23[i]
        j+=1
    end
end
better3

i=17
plotInt2(better2[i])
plotInt3(better2[i])
beautiful: better2[2,3,4,5,6,7,10,13]
i=17
plotInt2(better3[i])
plotInt3(better3[i])
beautiful: better3[2,4,5,13,17,20   ,109]


text:
using Revise
using OptiMode
using LinearAlgebra, Statistics, StaticArrays, HybridArrays, GeometryPrimitives, BenchmarkTools
using ChainRules, Zygote, ForwardDiff, FiniteDifferences

p = [
       1.7,                #   top ridge width         `w_top`         [μm]
       0.7,                #   ridge thickness         `t_core`        [μm]
       π / 14.0,           #   ridge sidewall angle    `θ`             [radian]
       # 0.5,                #   vacuum gap at boundaries `edge_gap`     [μm]
               ];
Δx,Δy,Δz,Nx,Ny,Nz = 6.0, 4.0, 1.0, 128, 128, 1;
Δx,Δy,Δz,Nx,Ny,Nz = 6.0, 4.0, 1.0, 256, 256, 1;
gr = Grid(Δx,Δy,Nx,Ny)
rwg(x) = ridge_wg(x[1],x[2],x[3],0.5,MgO_LiNbO₃,SiO₂,Δx,Δy) # fourth param is gap between modeled object and sidewalls
geom = rwg(p)
ms = ModeSolver(1.45, geom, gr)

## prototype backprop through NLsolve of PDE
# from https://github.com/JuliaNLSolvers/NLsolve.jl/issues/205

using NLsolve
using Zygote
using Zygote: @adjoint
using IterativeSolvers
using LinearMaps
using SparseArrays
using LinearAlgebra
using BenchmarkTools

# nlsolve maps f to the solution x of f(x) = 0
# We have ∂x = -(df/dx)^-1 ∂f, and so the adjoint is df = -(df/dx)^-T dx
@adjoint nlsolve(f, x0; kwargs...) =
    let result = nlsolve(f, x0; kwargs...)
        result, function(vresult)
            dx = vresult[].zero
            x = result.zero
            _, back_x = Zygote.pullback(f, x)

            JT(df) = back_x(df)[1]
            # solve JT*df = -dx
            L = LinearMap(JT, length(x0))
            df = gmres(L,-dx)

            _, back_f = Zygote.pullback(f -> f(x), f)
            return (back_f(df)[1], nothing, nothing)
        end
    end

const NN = 10000
const nonlin = 0.1
const A = spdiagm(0 => fill(10.0, NN), 1 => fill(-1.0, NN-1), -1 => fill(-1.0, NN-1))
const p0 = randn(NN)
f(x, p) = A*x + nonlin*x.^2 - p
solve_x(p) = nlsolve(x -> f(x, p), zeros(NN), method=:anderson, m=10).zero
obj(p) = sum(solve_x(p))

Zygote.refresh()
g_auto, = gradient(obj, p0)
g_analytic = gmres((A + Diagonal(2*nonlin*solve_x(p0)))', ones(NN))
display(g_auto)
display(g_analytic)
@show sum(abs.(g_auto - g_analytic))

@btime gradient(obj, p0);
@btime gmres((A + Diagonal(2*nonlin*solve_x(p0)))', ones(NN));

##
f(x, p) = A*x + nonlin*x.^2 - p
g(x) = f(x,p0)
x0 = solve_x(p0)
_, back_x = Zygote.pullback(g, x0)
_, back_f = Zygote.pullback(ff -> ff(x0), g)
_, back_f2 = Zygote.pullback(pp -> f(x0,pp), p0)
dx = ones(length(x0))
JT(df) = back_x(df)[1]
# solve JT*df = -dx
L = LinearMap(JT, length(x0))
df = gmres(L,-dx)
back_f(df)[1]
back_f2(df)[1]

## Try to rewrite _solve_Δω² & solve_k for use with nlsolve instead of roots
# and check if this adjoint still works
using Zygote: dropgrad

# f(k,ωₜ,ε⁻¹) = solve_ω²(dropgrad(ms),k,ε⁻¹)[1] - ωₜ^2

function findk(ms::ModeSolver{ND},om::Real,epsi::AbstractArray{<:SMatrix{3,3},ND},k₀::Real) where ND
    f(k,ωₜ,ε⁻¹) = solve_ω²(dropgrad(ms),k,ε⁻¹)[1] - ωₜ^2
    res = nlsolve(x->f(x[1],om,epsi),[k₀,])
    k = res.zero[1]
end

function findk(ms::ModeSolver{ND},om::Real,geom::AbstractVector{S},k₀::Real) where {ND,S<:Shape}
    epsi = εₛ⁻¹(om,geom;ms)
    f(k,ωₜ,ε⁻¹) = real(solve_ω²(dropgrad(ms),real(k),real.(ε⁻¹))[1]) - real(ωₜ)^2
    res = nlsolve(x->real(f(real(x[1]),real(om),real.(epsi))),[k₀,])
    k = real(res.zero[1])
end

ei = ms.M̂.ε⁻¹

findk(ms,1/1.55,ei,1.46)
findk(ms,1/1.58,ei,1.26)

findk(ms,1/1.61,rwg(p),1.23)
findk(ms,1/1.64,rwg(p),1.20)

Zygote.gradient((om,pp)->findk(ms,om,rwg(pp),1.20), 1/1.64,p)

##
om = 1/ 1.55
k₀ = 1.2
geom= rwg(p)
epsi = εₛ⁻¹(om,geom;ms)
f1(k,ωₜ,ε⁻¹) = real(solve_ω²(dropgrad(ms),real(k),real.(ε⁻¹))[1]) - real(ωₜ)^2
f2(x) = real(f1(real(x[1]),real(om),real.(epsi)))
res = nlsolve(f2,[k₀,])

dk = SMatrix{1,1}(1.0) #SVector{1}(1.0) #1.0
k = res.zero #[1]
f2_val, back_k = Zygote.pullback(f2, k)
JT(df) = back_k(df)[1]
L = SMatrix{1,1}(ms.∂ω²∂k[1]) # LinearMap(JT, 1)
df = gmres(L,-dk) # -1/ms.∂ω²∂k[1]
_, back_f = Zygote.pullback(ff->ff(k),f2)
back_f(df)[1]

f2(k)
Zygote.gradient(ff->ff(k),f2)

f1(1.2,1/1.55,epsi)
Zygote.gradient(f1,1.2,1/1.55,epsi)
f3(k,ωₜ,p) = real(solve_ω²(dropgrad(ms),real(k),εₛ⁻¹(ωₜ,rwg(p);ms))[1]) - real(ωₜ)^2
f3(1.2,1/1.55,p)
Zygote.gradient(f3,1.2,1/1.55,p)
f3(1.26,1/1.55,p)
Zygote.gradient(f3,1.26,1/1.55,p)
f4(x) = f3(x[1],1/1.55,p)
f4([1.26,])

res2 = nlsolve(f4,[k₀,])
k2 = res2.zero
f4_val, back_k2 = Zygote.pullback(f4, k2)
JT2(df) = back_k2(df)[1]
JT2([1.0])

JT([1.0+0.1im,])

Zygote.gradient(f4,k2)

ms.∂ω²∂k[1]
ms.∂ω²∂k
domsq,domsq_domsqp = _solve_Δω²(ms,k2[1],1/1.55)

Zygote.gradient((om,pp)->f3(1.26,om,pp),1/1.55,p)
_, solve_domsq_pb = Zygote.pullback((om,pp)->real(solve_ω²(dropgrad(ms),real(k[1]),εₛ⁻¹(om,rwg(pp);ms))[1]) - real(om)^2, 1/1.55, p)
solve_domsq_pb(-1/ms.∂ω²∂k[1])

FiniteDifferences.grad(central_fdm(3,1),x->solve_k(ms,1/1.55,rwg(x))[1],p)[1]

FiniteDifferences.grad(central_fdm(3,1),x->solve_k(ms,x,rwg(p))[1],1/1.55)[1]

_, solve_domsq_pb = Zygote.pullback(1/1.55, p) do om, pp
    omsq,HH = solve_ω²(dropgrad(ms),real(k[1]),εₛ⁻¹(om,rwg(pp);ms))
    domsq = real(omsq) - real(om)^2
end

Zygote.gradient(ff->ff(1.26),f4)

omsqH, omsqH_pb = Zygote.pullback(k[1], epsi) do kk, einv
    omsq,HH = solve_ω²(dropgrad(ms),real(kk),einv)
end

H̄₁ = randn(ComplexF64,length(omsqH[2]))

k̄₂,eī₁ = omsqH_pb( (0.0, H̄₁) )

ω²errH, ω²errH_pb = Zygote.pullback(k[1], 1/1.55, epsi) do kk, om_in, einv
    omsq_out,HH = solve_ω²(dropgrad(ms),real(kk),einv)
    omsq_err = real(omsq_out) - real(om_in)^2
    return (omsq_err, HH)
end

k̄₂2,ω̄₁ , eī₁2 = ω²errH_pb( (0.0, H̄₁) )
k̄₂3, ω̄₁2 , eī₁3 = ω²errH_pb( (-1 / ms.∂ω²∂k[1], nothing) )

k̄₂2 ≈ k̄₂
eī₁2 ≈ eī₁

k,H = solve_k(ms,0.8,epsi)
ω²err, ω²err_pb = Zygote.pullback(0.8, epsi) do om, einv
    omsq,HH = solve_ω²(dropgrad(ms),real(k),einv)
    omsq_err = real(omsq) - real(om)^2
end

ω²err_pb(1.0)

ω̄ , eīₖ = ω²err_pb( -1 / ms.∂ω²∂k[1]  )

ω̄  ≈ ω̄₁2
eīₖ ≈ eī₁3

solve_k

ext: jl


text:
using Test, OpenStreetMapX
import LightGraphs

@testset "maps" begin

m = get_map_data("data/reno_east3.osm",use_cache=false);

@test length(m.nodes) == 9032

using Random
Random.seed!(0);
pA = generate_point_in_bounds(m)
@test all(isapprox.(pA,(39.53584630184622, -119.71506095062803)))
pB = generate_point_in_bounds(m)
@test all(isapprox.(pB,(39.507242155639005, -119.78506509516248)))
pointA = point_to_nodes(pA, m)
pointB = point_to_nodes(pB, m)

@test pointA == 3052967037
@test pointB == 140393352

sr1, shortest_distance1, shortest_time1 = shortest_route(m, pointA, pointB)
@test (sr1[1], sr1[end]) == (pointA, pointB)

@test shortest_route(m, pointA, pointB; routing = :astar) == shortest_route(m, pointA, pointB; routing = :dijkstra)
@test fastest_route(m, pointA, pointB; routing = :astar) == fastest_route(m, pointA, pointB; routing = :dijkstra)

#######################################################
### nodes.jl ###
#nodes.jl/nearest_node
@test nearest_node(m.nodes, ENU(3300,1500,-1)) == 3052967130

#nodes.jl/centroid
node_list = [pointA,pointB]
@test centroid(m.nodes,node_list) == ENU(309.1010361600104, -35.35610349960848, -0.9191984126487114)

#nodes.jl/nodes_within_range
@test nodes_within_range(m.nodes, ENU(3300,1500,-1), 100.0) == [3052967037, 3052967180, 3052967130, 6050217400, 3052966904, 3052967140, 6050217409, 3052967199]

#nodes.jl/add_new_node!
@test haskey(m.nodes,-1) == false
OpenStreetMapX.add_new_node!(m.nodes,ENU(3300,1500,-1),-1)
@test length(m.nodes) == 9033
@test nearest_node(m.nodes,ENU(3300,1500,-1)) == -1

### intersections.jl ###
#intersections.jl/oneway
@test OpenStreetMapX.oneway(m.roadways[1]) == true

#intersections.jl/distance
@test distance(m.nodes,sr1) == 9019.07204040599

#intersections.jl/find_intersections
@test OpenStreetMapX.find_intersections(m.roadways[1:2]) == Dict(139988738=>Set([14370413]),2975020216=>Set([14370407]),2441017888=>Set([14370407]),385046328=>Set([14370413]))

### routing.jl ###
#routing.jl/get_edges
@test OpenStreetMapX.get_edges(m.nodes,m.roadways[1:2]) == (Tuple{Int64,Int64}[(139988738, 385046327), (2441017870, 2975020216), (385046327, 385046328), (2441017888, 2441017878), (2441017878, 2441017870)], [1, 4, 1, 4, 4])

#routing.jl/get_vertices
@test OpenStreetMapX.get_vertices(OpenStreetMapX.get_edges(m.nodes,m.roadways[1:2])[1]) == Dict(2441017878=>7,139988738=>1,2975020216=>4,2441017870=>3,385046328=>5,2441017888=>6,385046327=>2)

#routing.jl/distance
#Returns seem to be equal yet returning false (?)
@test distance(m.nodes,OpenStreetMapX.get_edges(m.nodes,m.roadways[1:2])[1]) == [30.2013937293296, 7.243941886194111, 35.492758006997796, 12.29992029473937, 11.290063259013777]

#intersections.jl/features_to_graph

#######################################################

function perftest()
  sr_len=length(shortest_route(m, pointA, pointB; routing = :astar)[1])
  fr_len=length(fastest_route(m, pointA, pointB; routing = :astar)[1])
  shortest_route(m, pointA, pointB; routing = :dijkstra)
  fastest_route(m, pointA, pointB; routing = :dijkstra)

  print("shortest_route(...; routing = :astar) of $sr_len nodes")
  @time shortest_route(m, pointA, pointB; routing = :astar)
  print("shortest_route(...; routing = :dijkstra)  of $sr_len nodes")
  @time shortest_route(m, pointA, pointB; routing = :dijkstra)

  print("fastest_route(...; routing = :astar) of $fr_len nodes")
  @time fastest_route(m, pointA, pointB; routing = :astar)
  print("fastest_route(...; routing = :dijkstra) of $fr_len nodes")
  @time fastest_route(m, pointA, pointB; routing = :dijkstra)

end

perftest()

end;


text:
module EMLHeatExchangerDetailed_cost
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...heat_exchangers.EMLHeatExchangerDetailed
	import EMLtypes.length
	include("HeatExchangerDetailed_cost/ShellandTubes_NTU_cost.jl")
	include("HeatExchangerDetailed_cost/ShellandTubes_LMTD_cost.jl")
end


text:
"""
    This package has been inspired and even sometimes directly copied from two projects:
    - https://github.com/gridap/GridapPETSc.jl
    - https://github.com/JuliaParallel/PETSc.jl
"""
module PetscWrap

using Libdl
using MPI

include("const_arch_ind.jl")
export  PetscErrorCode, PETSC_DECIDE
# export all items of some enums
for item in Iterators.flatten((instances(InsertMode), instances(MatAssemblyType), instances(PetscViewerFormat)))
    @eval export $(Symbol(item))
end

include("load.jl")

include("const_arch_dep.jl")
export PetscReal, PetscScalar, PetscInt, PetscIntOne

include("init.jl")
export PetscInitialize, PetscFinalize

include("viewer.jl")
export  PetscViewer, CViewer,
        PetscViewerASCIIOpen,
        PetscViewerCreate,
        PetscViewerDestroy,
        PetscViewerPopFormat,
        PetscViewerPushFormat,
        PetscViewerStdWorld

include("vec.jl")
export  PetscVec, CVec,
        VecAssemble,
        VecAssemblyBegin,
        VecAssemblyEnd,
        VecCreate,
        VecDestroy,
        VecDuplicate,
        VecGetArray,
        VecGetLocalSize,
        VecGetOwnershipRange,
        VecGetSize,
        VecRestoreArray,
        VecSetFromOptions,
        VecSetSizes,
        VecSetUp,
        VecSetValue,
        VecSetValues,
        VecView

include("mat.jl")
export  PetscMat, CMat,
        MatAssemble,
        MatAssemblyBegin,
        MatAssemblyEnd,
        MatCreate,
        MatCreate,
        MatCreateVecs,
        MatDestroy,
        MatGetOwnershipRange,
        MatSetFromOptions,
        MatSetSizes,
        MatSetUp,
        MatSetValue,
        MatSetValues,
        MatView,
        MatCreateComposite

include("ksp.jl")
export  PetscKSP, CKSP,
        KSPCreate,
        KSPDestroy,
        KSPSetFromOptions,
        KSPSetOperators,
        KSPSetUp,
        KSPSolve

# fancy
include("fancy/vec.jl")
export  assemble!,
        create_vector,
        destroy!,
        duplicate,
        get_range,
        set_local_size!, set_global_size!,
        set_from_options!,
        set_up!,
        vec2array

include("fancy/mat.jl")
export  create_matrix,
        set_values!

include("fancy/ksp.jl")
export  create_ksp,
        solve, solve!,
        set_operators!

end


text:
"""
    runhandlers(handlers, eff)
    runhandlers((Vector, Option), eff)::Vector{Option{...}}

run all handlers such that the first handler will define the most outer container
"""
runhandlers(single_handler, eff::Eff) = runlast_ifpossible(runhandler(single_handler, eff))
runhandlers(all_handlers::Vector, eff::Eff) = runhandlers(tuple(all_handlers...), eff)
runhandlers(all_handlers::Tuple, eff::Eff) = runlast_ifpossible(_runhandlers(all_handlers, eff))
_runhandlers(all_handlers::Tuple{}, eff::Eff) = eff
function _runhandlers(all_handlers::Tuple, eff::Eff)
  subresult = _runhandlers(Base.tail(all_handlers), eff)
  runhandler(first(all_handlers), subresult)
end

runhandlers(any, not_eff) = error("can only apply runhandlers onto an `Eff`, got a `$(typeof(not_eff))`")

"""
    runlast(eff::Eff)

extract final value from Eff with all effects (but Identity) already run
"""
function runlast(final::Eff)
  @assert isempty(final.cont) "expected eff without continuation, but found cont=$(final.cont)"
  @assert final.effectful isa NoEffect "not all effects have been handled, found $(final.value)"
  final.effectful.value
end

"""
    runlast_ifpossible(eff::Eff)

like `ExtensibleEffects.runlast`, however if the Eff is not yet completely handled, it just returns it.
"""
function runlast_ifpossible(final::Eff)
  if isempty(final.cont) && final.effectful isa NoEffect
    final.effectful.value
  else
    final
  end
end

"""
    runhandler(handler, eff::Eff)
    runhandler(handler, eff::Eff, context)

key method to run an effect on some effecthandler Eff

note that we represent effectrunners as plain types in order to associate
standard effect runners with types like Vector, Option, ...
"""
function runhandler(handler, eff::Eff)
  eff_applies(handler, eff.effectful) || return runhandler_not_applies(handler, eff)

  interpreted_continuation = if isempty(eff.cont)
    # you may think we could simplify this, however for eff `eff_flatmap(handler, x -> eff_pure(handler,  x), eff) != eff`
    # because there is the handler which may have extra information
    Continuation(x -> _eff_pure(handler, x))
  else
    Continuation(x -> runhandler(handler, eff.cont(x)))
  end
  _eff_flatmap(handler, interpreted_continuation, eff.effectful)
end

"""
    runhandler_not_applies(handler, eff)

if your handler does not apply, use this as the fallback to handle the unknown effect.
"""
function runhandler_not_applies(handler, eff::Eff)
  interpreted_continuation = if isempty(eff.cont)
    Continuation(x -> _eff_pure(handler, x))
  else
    Continuation(x -> runhandler(handler, eff.cont(x)))
  end
  # if we don't know how to handle the current eff, we return it with the new continuation
  # this ensures the handler is applied recursively
  Eff(eff.effectful, interpreted_continuation)
end

"""
    @runhandlers handlers eff

For convenience we provide `runhandlers` function also as a macro.

With this you can easier run left-over handlers from an `@syntax_eff` autorun.

Example
-------
```
@runhandlers WithCall(args, kwargs) @syntax_eff begin
  a = Callable(x -> 2x)
  @pure a
end
```
"""
macro runhandlers(handlers, eff)
  esc(:(ExtensibleEffects.runhandlers($handlers, $eff)))
end

# effecthandler interface
# =======================
# to support an effect, a type needs to implement three functions
# 1. eff_flatmap
# 2. eff_pure
# 3. eff_applies

# eff_flatmap
# -----------

function _eff_flatmap(handler, interpreted_continuation::Continuation, value)
  # provide convenience wrapper if someone forgets to return an Eff
  result = eff_flatmap(handler, interpreted_continuation, value)
  noeffect(result)
end

"""
    ExtensibleEffects.eff_flatmap(handler, interpreted_continuation, value)
    ExtensibleEffects.eff_flatmap(interpreted_continuation, value)

Overwrite this for your custom effect handler to handle your effect.
This function is only called if `eff_applies(handler, value)==true`.

While for custom effects it is handy to dispatch on the handler itself, in simple cases
`handler == typeof(value)` and hence, we allow to ommit it.

Parameters
----------
The arg `interpreted_continuation` is guaranteed to return an Eff of the handled type.
E.g. if you might handle the type `Vector`, you are guaranteed that `interpreted_continuation(x)::ExtensibleEffects.Eff{Vector}`

Return
------
If you do not return an `ExtensibleEffects.Eff`, the result will be wrapped into `noeffect` automatically,
i.e. assuming the effect is handled afterwards.
"""
eff_flatmap(handler, interpreted_continuation, value) = eff_flatmap(interpreted_continuation, value)

# eff_pure
# --------

function _eff_pure(handler, value)
  result = eff_pure(handler, value)
  noeffect(result)
end

"""
    ExtensibleEffects.eff_pure(handler, value)

Overwrite this for your custom effect handler, return either `ExtensibleEffects.Eff` type, or a plain value.
Plain values will be wrapped with `noeffect` automatically.
"""
function eff_pure end

# eff_applies
# -----------

"""
    ExtensibleEffects.eff_applies(handler::YourHandlerType, value::ToBeHandledEffectType) = true

Overwrite this function like above to indicate that a concrete effect is handled by a handler.
In most cases you will have `YourHandlerType = Type{ToBeHandledEffectType}`, like for `Vector` or similar.

Sometimes you need extra information without which you cannot run a specific effect. Then you need to link
the specific handler containing the required information. E.g. `Callable` needs `args` and `kwargs` to be run,
which are captured in the handler type `CallableHandler(args, kwargs)`.
Hence above you would choose YourHandlerType = `CallableHandler`, and ToBeHandledEffectType = `Callable`.
"""
eff_applies(handler, value) = false

max_stars_repo_name: schlichtanders/ExtensibleEffects

max_issues_repo_name: schlichtanders/ExtensibleEffects

max_forks_repo_name: schlichtanders/ExtensibleEffects


text:
using Flux, DiffEqFlux, DifferentialEquations, Plots

## Setup DDE to optimize
function delay_lotka_volterra(du,u,h,p,t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = (α - β*y)*h(p,t-0.1)[1]
  du[2] = dy = (δ*x - γ)*y
end

h(p,t) = ones(eltype(p),2)
prob = DDEProblem(delay_lotka_volterra,[1.0,1.0],h,(0.0,10.0),constant_lags=[0.1])
p = param([2.2, 1.0, 2.0, 0.4])
function predict_rd_dde()
  diffeq_rd(p,prob,MethodOfSteps(Tsit5()),saveat=0.1)[1,:]
end
loss_rd_dde() = sum(abs2,x-1 for x in predict_rd_dde())

# Optimize the parameters so the DDE's solution stays near 1

data = Iterators.repeated((), 100)
opt = ADAM(0.1)
cb = function () #callback function to observe training
  display(loss_rd_dde())
  # using `remake` to re-create our `prob` with current parameters `p`
  display(plot(solve(remake(prob,p=Flux.data(p)),MethodOfSteps(Tsit5()),saveat=0.1),ylim=(0,6)))
end
# Display the DDE with the initial parameter values.
cb()
Flux.train!(loss_rd_dde, [p], data, opt, cb = cb)

max_stars_repo_name: c-p-murphy/model-zoo

max_issues_repo_name: c-p-murphy/model-zoo

max_forks_repo_name: c-p-murphy/model-zoo


text:
# example 4: read with iterators,
#            custom return type

module FastaIOExample4
using FastaIO

function read_fasta_file(filename::String)
    # We set the reader so that it returns a Vector{Char};
    # If you have BioSeq package, you can use more
    # specialized uints, e.g. Vector{AminoAcid}
    fr = FastaReader{Vector{Char}}(filename)

    for (desc, seq) in fr
        println("num=$(fr.num_parsed) desc=$desc seq=$seq")
    end
    println("read $(fr.num_parsed) entries")
end

examples_dir = dirname(Base.source_path())
read_fasta_file(joinpath(examples_dir, "example.fasta.gz"))
end


text:
using Distributions, Distances
export SIS, SisPosterior, DimensionRange

struct DimensionRange
    l::Real
    u::Real
end
# NOTE - both these functions perform worse if called repeatedly.
# If we have r::DimensionRange better to assign values such as..
# l, u = (r.l, r.u)
# and then use these in the code. This is what I have done in MCMC samplers.
function isin(
    r::DimensionRange,
    val::Real
    )
    return (val >= r.l) & (val <= r.u)
end

function notin(
    r::DimensionRange,
    val::Real
    )
    return (val < r.l) | (val > r.u)
end

struct SIS{T<:Metric}
    mode::Vector{Path{Int}} # Mode
    γ::Float64 # Precision
    dist::T # Distance metric
    V::UnitRange{Int} # Vertex Set
    K_inner::DimensionRange # Maximum interaction sequence size
    K_outer::DimensionRange # Maximum path (interaction) length
end

SIS(
    mode::InteractionSequence{Int},
    γ::Float64,
    dist::Metric,
    V::UnitRange{Int}
) = SIS(
        mode, γ, dist, V,
        DimensionRange(1,Inf),
        DimensionRange(1,Inf)
)

SIS(
    mode::InteractionSequence{Int},
    γ::Float64,
    dist::Metric,
    V::UnitRange{Int},
    K_inner::Real, K_outer::Real
) = SIS(
        mode, γ, dist, V,
        DimensionRange(1,K_inner),
        DimensionRange(1,K_outer)
)

function Base.show(
    io::IO, model::SIS
    )

    title = "$(typeof(model))"
    n = length(title)
    println(io, title)
    println(io, "-"^n)
    println(io, "Parameters:")
    for par in fieldnames(typeof(model))
        println(io, par, " = $(getfield(model, par))")
    end
end

function Base.similar(
    model::SIS{T},
    mode::Vector{Vector{Int}},
    γ::Float64
    ) where {T}
    return SIS(mode, γ, model.dist, model.V, model.K_inner, model.K_outer)
end

struct SisPosterior
    data::InteractionSequenceSample{Int}
    S_prior::SIS
    γ_prior::ContinuousUnivariateDistribution
    dist::Metric
    V::UnitRange{Int}
    K_inner::DimensionRange
    K_outer::DimensionRange
    sample_size::Int
    function SisPosterior(
        data::InteractionSequenceSample{Int},
        S_prior::SIS,
        γ_prior::ContinuousUnivariateDistribution
        )

        dist = S_prior.dist
        V = S_prior.V
        K_inner = S_prior.K_inner
        K_outer = S_prior.K_outer
        sample_size = length(data)
        new(data, S_prior, γ_prior, dist, V, K_inner, K_outer, sample_size)
    end
end

max_stars_repo_path: src/Metric Models/SIS.jl

max_issues_repo_path: src/Metric Models/SIS.jl

max_forks_repo_path: src/Metric Models/SIS.jl


text:
#Define the baseline ct at beginning
baseline_ct = ones(24)

#Define a penalty for having too many turning points in ct fit
function ct_penalty(λ,cts)
        λ*sum([(cts[t-1] + cts[t+1] - 2*cts[t])^2 for t = 2:(length(cts)-1)])
end

"""
        function ct_neg_log_density_grad!(grad,log_cts,θ,model::KenyaCoVSD.CoVAreaModel,baseline_contact,λ,p_waning)

In-place calculation of gradient of contact rate negative log-posterior density (first `length(baseline_contact)` days are fixed as Google estimate).
With sero-waning `p_waning`. This is based on the log-ct value (to maintain positivity)
"""
function ct_neg_log_density_grad!(grad,log_cts,θ,model::KenyaCoVSD.CoVAreaModel,baseline_contact,λ,p_waning)
        l = log_cts -> -ll_onegroup_newvariant_infboost_ct(θ,model,p_waning,vcat(baseline_contact,exp.(log_cts))) + ct_penalty(λ,vcat(baseline_contact,exp.(log_cts)))
        ForwardDiff.gradient!(grad,l,log_cts)
        return nothing
end

"""
	function random_grad_log_ct!(grad,x₀,model_ct,baseline_contact,λ,trans_name,p_serowaning)

In-place calculation for gradient of the log-contact rate with sero-waning for a randomly chosen draw from the MCMC
"""
function random_grad_log_ct!(grad,log_cts,model::KenyaCoVSD.CoVAreaModel,baseline_contact,λ,trans_name,p_serowaning)
	k = rand(1:size(model.MCMC_results.chain,1))
	θ = NamedTuple{Tuple(keys(model.MCMC_results.chain))}([model.MCMC_results.chain[k,n,1] for n = 1:size(model.MCMC_results.chain,2)])
	ct_neg_log_density_grad!(grad,log_cts,θ,model,baseline_contact,λ,p_serowaning)
	return nothing
end

"""
	ADAM_step!(X_n::Vector{Float64},
						∇V::Vector{Float64},
						∇V_n::Vector{Float64},
						∇V_sq_n::Vector{Float64},
						n::Int64,
						ADAM_params::NamedTuple)

In-place calculation of parameter, gradient and squared gradient estimates using ADAM SGD method.
"""
function ADAM_step!(X_n::Vector{Float64},
					∇V::Vector{Float64},
					∇V_n::Vector{Float64},
					∇V_sq_n::Vector{Float64},
					hat_∇V_n::Vector{Float64},
					hat_∇V_sq_n::Vector{Float64},
					n::Int64,
					ADAM_params::NamedTuple)
	@unpack β₁,β₂,η,ϵ = ADAM_params
	∇V_n .= β₁.*∇V_n .+ (1-β₁).*∇V
	∇V_sq_n .= β₂.*∇V_sq_n .+ (1-β₂).*(∇V.^2)
	hat_∇V_n .= ∇V_n./(1 - β₁^(n+1))
	hat_∇V_sq_n .= ∇V_sq_n./(1 - β₂^(n+1))
	X_n .= X_n .- η.*hat_∇V_n./(sqrt.(hat_∇V_sq_n) .+ ϵ)
	return nothing
end

"""
	function ADAM_optim(X_opt::Vector{Float64},
                        grad_func!::Function,
                        ADAM_params::NamedTuple,
                        model::KenyaSerology.CoVAreaModel,
                        baseline_contact::Vector{Float64},
                        trans_name,
                        p_waning)

Optimize the vector initialized at `X_opt` using ADAM method with `grad_func!` gradient with sero-waning.
"""
function ADAM_optim(X_opt::Vector{Float64},
			grad_func!::Function,
			ADAM_params::NamedTuple,
			model::KenyaCoVSD.CoVAreaModel,
			baseline_contact::Vector{Float64},
			trans_name,
			p_waning)
	println("Starting ADAM SGD optimisation for effective contact rates.")
	@unpack total_num_steps,averaging_num_steps,λ = ADAM_params

	∇V_n = zeros(length(X_opt))
	∇V_sq_n = zeros(length(X_opt))
	hat_∇V_n = zeros(length(X_opt))
	hat_∇V_sq_n = zeros(length(X_opt))
	grad = zeros(length(X_opt))
	X̂ = zeros(length(X_opt))

	for n = 0:(total_num_steps-1)
		if n%10 == 0
			println("On ADAM step $(n)")
		end
		grad_func!(grad,X_opt,model,baseline_contact,λ,trans_name,p_waning)
		ADAM_step!(X_opt,
					grad,
					∇V_n,∇V_sq_n,
					hat_∇V_n,hat_∇V_sq_n,
					n,ADAM_params)
		if n > total_num_steps - averaging_num_steps #
			m = n - (total_num_steps - averaging_num_steps)
			X̂ .+= (X_opt - X̂)/m
		end

	end
	final_LL = 0
	for k = 1:size(model.MCMC_results.chain,1)
		θ = NamedTuple{Tuple(keys(model.MCMC_results.chain))}([model.MCMC_results.chain[k,n,1] for n = 1:size(model.MCMC_results.chain,2)])
		final_LL += ll_onegroup_newvariant_infboost_ct(θ,model,p_waning,vcat(baseline_contact,exp.(X̂)))
	end
	final_LL = (final_LL/size(model.MCMC_results.chain,1)) - ct_penalty(λ,vcat(baseline_contact,exp.(X̂)))

	return X̂,final_LL
end


text:
module ClusteredGa

using StatsBase
using Clustering

import Base.Threads

const GA_TYPE_CLUSTER = 0
const GA_TYPE_CLUSTER_SIM = 1
const GA_TYPE_CLASSIC = 2

mutable struct Chromosome
    genes::Array{Float64, 1}
    cost::Float64
    clusterid::Int
end

struct GAResult
    population::Array{Chromosome, 1}
    best::Chromosome
    generations::Int
end

# Create chromosome with given genes, cost of Infinity,
# and cluster id of -1.
function Chromosome(genes::Array{Float64, 1})
    return Chromosome(
        genes,
        Inf64,
        -1
    )
end

function Chromosome(lower::Array{Float64, 1}, upper::Array{Float64, 1})
    L = length(lower)
    genes = zeros(Float64, L)
    for i in 1:L
        genes[i] = lower[i] + rand() * (upper[i] - lower[i])
    end
    return Chromosome(genes)
end

function uniformcrossover(ch1::Chromosome, ch2::Chromosome)::Chromosome
    L = length(ch1.genes)
    offspring = Chromosome(zeros(Float64, L))
    for i in 1:L
        if rand() < 0.5
            offspring.genes[i] = ch1.genes[i]
        else
            offspring.genes[i] = ch2.genes[i]
        end
    end
    return offspring
end

function makeuniformcrossover()::Function
    function tempfn(ch1::Chromosome, ch2::Chromosome)
        return uniformcrossover(ch1, ch2)
    end
    return tempfn
end

function weightedcrossover(ch1::Chromosome, ch2::Chromosome)::Chromosome
    alpha = rand()
    genes = alpha .* ch1.genes .+ (1.0 - alpha) .* ch2.genes
    return Chromosome(genes)
end

function makeweightedcrossover()::Function
    function tempfn(ch1::Chromosome, ch2::Chromosome)
        return weightedcrossover(ch1, ch2)
    end
    return tempfn
end

function linearcrossover(costfn::Function, ch1::Chromosome, ch2::Chromosome)::Chromosome
    genes1 = 0.5 .* ch1.genes .+ 0.5 .* ch2.genes
    genes2 = 1.5 .* ch1.genes .- 0.5 .* ch2.genes
    genes3 = 1.5 .* ch2.genes .- 0.5 .* ch1.genes
    allgenes = [genes1, genes2, genes3]
    costs = map(costfn, allgenes)
    bestindice = costs |> sortperm |> first
    return Chromosome(
        allgenes[bestindice],
        costs[bestindice],
        -1
    )
end

function makelinearcrossover(costfn::Function)::Function
    function tmpfn(ch1::Chromosome, ch2::Chromosome)::Chromosome
        return linearcrossover(costfn, ch1, ch2)
    end
    return tmpfn
end

function blxalphacrossover(ch1::Chromosome, ch2::Chromosome)::Chromosome
    L = length(ch1.genes)
    offspring = Chromosome(zeros(Float64, L))
    for i in 1:L
        d = abs(ch1.genes[i] - ch2.genes[i])
        alpha = rand()
        ad = d * alpha
        umin = min(ch1.genes[i], ch2.genes[i]) - ad
        umax = max(ch1.genes[i], ch2.genes[i]) + ad
        u = umin + rand() * (umax - umin)
        offspring.genes[i] = u
    end
    return offspring
end

function makeblxalphacrossover()::Function
    function tmpfn(ch1::Chromosome, ch2::Chromosome)::Chromosome
        return blxalphacrossover(ch1, ch2)
    end
    return tmpfn
end

function normalmutation(stddev::Float64, mutationprob::Float64, ch::Chromosome)::Chromosome
    newgenes = copy(ch.genes)
    for i in 1:length(length(newgenes))
        if rand() < mutationprob
            newgenes[i] += randn() * stddev
        end
    end
    return Chromosome(newgenes)
end

function makenormalmutation(stddev::Float64, mutationprob::Float64)::Function
    function tmpfn(ch::Chromosome)::Chromosome
        return normalmutation(stddev, mutationprob, ch)
    end
    return tmpfn
end

function randommutation(lower::Array{Float64, 1}, upper::Array{Float64, 1}, mutationprob::Float64, ch::Chromosome)::Chromosome
    newgenes = zeros(Float64, length(ch.genes))
    for i in 1:length(newgenes)
        if rand() < mutationprob
            newgenes[i] = lower[i] + rand() * (upper[i] - lower[i])
        else
            newgenes[i] = ch.genes[i]
        end
    end
    return Chromosome(newgenes)
end

function makerandommutation(lower::Array{Float64, 1}, upper::Array{Float64, 1}, mutationprob::Float64)::Function
    function tmpfn(ch::Chromosome)::Chromosome
        return randommutation(lower, upper, mutationprob, ch)
    end
end

function tournamentselection(pop::Array{Chromosome, 1}, k::Int)::Array{Chromosome, 1}
    fathers = sort(sample(pop, k, replace = false), by = ch -> ch.cost)
    mothers = sort(sample(pop, k, replace = false), by = ch -> ch.cost)
    return [fathers[1], mothers[1]]
end

function maketournamentselection(k::Int)::Function
    function tmpfn(pop::Array{Chromosome, 1})::Array{Chromosome, 1}
        return tournamentselection(pop, k)
    end
end

function distance(ch1::Chromosome, ch2::Chromosome)::Float64
    (ch1.genes .- ch2.genes) .|> (x -> x * x) |> sum
end

function simulatedkmeanstournamentselection(pop::Array{Chromosome, 1}, k::Int)::Array{Chromosome, 1}
    fathers = sort(sample(pop, k, replace = false), by = ch -> ch.cost)
    mothers = sort(sample(pop, k, replace = false), by = ch -> ch.cost)

    thefather = fathers[1]
    distances = map(ch -> distance(thefather, ch), mothers)
    furthestindice = distances |> sortperm |> last
    themother = mothers[furthestindice]

    return [thefather, themother]
end

function makesimulatedkmeanstournamentselection(k::Int)::Function
    function tmpfn(pop::Array{Chromosome, 1})::Array{Chromosome, 1}
        return simulatedkmeanstournamentselection(pop, k)
    end
end

function kmeanstournamentselection(pop::Array{Chromosome, 1}, k::Int)::Array{Chromosome, 1}
    father = sort(sample(pop, k, replace = false), by = ch -> ch.cost)[1]
    mothers = filter(ch -> ch.clusterid != father.clusterid, pop)
    if length(mothers) == 0
        mothers = pop
    end
    mother = sort(sample(mothers, k, replace = true), by = ch -> ch.cost)[1]
    return [father, mother]
end

function makekmeanstournamentselection(k::Int)::Function
    function tmpfn(pop::Array{Chromosome, 1})::Array{Chromosome, 1}
        return kmeanstournamentselection(pop, k)
    end
end

function randompopulation(popsize::Int, lower::Array{Float64, 1}, upper::Array{Float64, 1})::Array{Chromosome, 1}
    population = Array{Chromosome, 1}(undef, popsize)
    for i in 1:popsize
        population[i] = Chromosome(lower, upper)
    end
    return population
end

function calculatefitness(population::Array{Chromosome, 1}, costfn::Function)
    Threads.@threads for ch in population
        ch.cost = costfn(ch.genes)
    end
end

function population2matrix(population::Array{Chromosome, 1})::Array{Float64, 2}
    popsize = length(population)
    chsize = length(population[1].genes)
    mat = Array{Float64, 2}(undef, popsize, chsize)
    for i in 1:popsize
        mat[i,:] = population[i].genes
    end
    return mat
end

function assignclusterids(population::Array{Chromosome, 1})
    mat = population2matrix(population)
    n, _ = size(mat)
    k = Int(ceil(sqrt(n)))
    clust::ClusteringResult = kmeans(transpose(mat), k)
    clusterids = clust.assignments
    for i in 1:n
        population[i].clusterid = clusterids[i]
    end
end

function generation(
    population::Array{Chromosome, 1}, costfn::Function, crossfn::Function,
    mutatefn::Function, gatype::Int; elitism::Int = 0)::Array{Chromosome, 1}

    popsize = length(population)
    calculatefitness(population, costfn)
    if gatype == GA_TYPE_CLUSTER
        assignclusterids(population)
        selectfn = makekmeanstournamentselection(3)
    elseif  gatype == GA_TYPE_CLASSIC
        selectfn = maketournamentselection(3)
    elseif gatype == GA_TYPE_CLUSTER_SIM
        selectfn = makesimulatedkmeanstournamentselection(3)
    end
    newpop = Array{Chromosome, 1}(undef, 0)

    if elitism > 0
        sort!(population, by = ch -> ch.cost)
        for i in 1:elitism
            push!(newpop, population[i])
        end
    end

    for _ in (elitism + 1):popsize
        father, mother = selectfn(population)
        offspring = crossfn(father, mother) |> mutatefn
        push!(newpop, offspring)
    end

    return newpop
end

function ga(
    popsize::Int, generations::Int, lower::Array{Float64, 1}, upper::Array{Float64, 1},
    costfn::Function, crossfn::Function, mutatefn::Function, gatype::Int;
    elitism::Int = 0, initialpopulation = nothing, verbose = true, targetcost = nothing)

    generations_done = 0

    if isnothing(targetcost)
        targetcost = generations
    end

    if isnothing(initialpopulation)
        population = randompopulation(popsize, lower, upper)
    else
        population = initialpopulation
    end

    for g in 1:generations
        population = generation(population, costfn, crossfn, mutatefn, gatype, elitism = elitism)
        if verbose
            @info g population[1].cost
        end
        bestcost = costfn(population[1].genes)
        generations_done = generations_done + 1
        if bestcost < targetcost
            break
        end
    end

    calculatefitness(population, costfn)
    sort!(population, by = ch -> ch.cost)
    return GAResult(
        population,
        population[1],
        generations_done
    )
end

function hybridga(popsize::Int, generations::Array{Int, 1}, lower::Array{Float64, 1}, upper::Array{Float64, 1},
    costfn::Function, crossfn::Function, mutatefn::Function; elitism::Int = 0)

    population = randompopulation(popsize, lower, upper)

    for _ in 1:generations[1]
        population = generation(population, costfn, crossfn, mutatefn, GA_TYPE_CLUSTER, elitism = elitism)
    end

    for _ in 1:generations[2]
        population = generation(population, costfn, crossfn, mutatefn, GA_TYPE_CLASSIC, elitism = elitism)
    end

    calculatefitness(population, costfn)
    sort!(population, by = ch -> ch.cost)
    return population
end

end # end of module ClusteredGa

max_stars_repo_name: jbytecode/divpopga

max_issues_repo_name: jbytecode/divpopga

max_forks_repo_name: jbytecode/divpopga


text:
"""
    UniformBinning(nbins::Int)

Binning scheme of the probability simplex with `nbins` bins of uniform width for
each component.
"""
struct UniformBinning <: AbstractBinningAlgorithm
    nbins::Int

    function UniformBinning(nbins::Int)
        nbins > 0 || error("number of bins must be positive")

        return new(nbins)
    end
end

function perform(
    binning::UniformBinning,
    predictions::AbstractVector{<:Real},
    targets::AbstractVector{Bool},
)
    @unpack nbins = binning

    # create dictionary of bins
    T = eltype(float(zero(eltype(predictions))))
    bins = Dict{Int,Bin{T}}()

    # reserve some memory (very rough guess)
    nsamples = length(predictions)
    sizehint!(bins, min(nbins, nsamples))

    # for all other samples
    @inbounds for (prediction, target) in zip(predictions, targets)
        # compute index of bin
        index = binindex(prediction, nbins)

        # create new bin or update existing one
        bin = get(bins, index, nothing)
        if bin === nothing
            bins[index] = Bin(prediction, target)
        else
            adddata!(bin, prediction, target)
        end
    end

    return values(bins)
end

function perform(
    binning::UniformBinning,
    predictions::AbstractVector{<:AbstractVector{<:Real}},
    targets::AbstractVector{<:Integer},
)
    return _perform(binning, predictions, targets, Val(length(predictions[1])))
end

function _perform(
    binning::UniformBinning,
    predictions::AbstractVector{<:AbstractVector{T}},
    targets::AbstractVector{<:Integer},
    nclasses::Val{N},
) where {T<:Real,N}
    @unpack nbins = binning

    # create bin for the initial sample
    binindices = NTuple{N,Int}[binindex(predictions[1], nbins, nclasses)]
    bins = [Bin(predictions[1], targets[1])]

    # reserve some memory (very rough guess)
    nsamples = length(predictions)
    guess = min(nbins, nsamples)
    sizehint!(bins, guess)
    sizehint!(binindices, guess)

    # for all other samples
    @inbounds for i in 2:nsamples
        # obtain prediction and corresponding target
        prediction = predictions[i]
        target = targets[i]

        # compute index of bin
        index = binindex(prediction, nbins, nclasses)

        # create new bin or update existing one
        j = searchsortedfirst(binindices, index)
        if j > length(binindices) || (binindices[j] !== index)
            insert!(binindices, j, index)
            insert!(bins, j, Bin(prediction, target))
        else
            bin = bins[j]
            adddata!(bin, prediction, target)
        end
    end

    return bins
end

function binindex(
    probs::AbstractVector{<:Real}, nbins::Int, ::Val{N}
)::NTuple{N,Int} where {N}
    ntuple(N) do i
        binindex(probs[i], nbins)
    end
end

function binindex(p::Real, nbins::Int)
    # check argument
    zero(p) ≤ p ≤ one(p) || throw(ArgumentError("predictions must be between 0 and 1"))

    # handle special case p = 0
    iszero(p) && return 1

    return ceil(Int, nbins * p)
end


text:
module VLConstraintBasedModelGenerationUtilities

# Include my files -
include("Include.jl")

# export methods and types -
export check
export build_gene_table
export build_protein_table
export build_metabolic_reaction_table
export build_transcription_reaction_table
export build_translation_reaction_table
export build_transport_reaction_table
export build_txtl_transport_reaction_table
export transcribe
export translate
export count

# metabolic methods -
export build_reaction_id_array
export build_flux_bounds_array
export build_species_symbol_array
export build_species_bounds_array
export build_species_table
export build_species_symbol_table
export build_stoichiometric_matrix
export build_atom_matrix
export build_objective_coefficient_array
export parse_molecular_formula_string

# file io methods -
export write_system_model_file

# private: but has a unit test
# export extract_stochiometric_coefficient_from_phrase

# export complement function -
export !

end # module


text:
### TODO : zseries linescans-pointscans support. Voltage stim and recording import.

### Reads a Prairie xml file and extracts the useful metadata.
function importPrairie(xmlFile)
    scanConfigDoc = root(parse_file(xmlFile))
    baseDir = dirname(xmlFile)

    prairieRelease = attribute(scanConfigDoc,"version")

    expDate = attribute(scanConfigDoc,"date")

   ### Reading the file per repeat
    runElements = get_elements_by_tagname(scanConfigDoc,"Sequence")
    seqParameters = [dictSequenceParameters(runE) for runE=runElements]

    ### Getting frame related information
    frames= Array{Dict}(length(seqParameters))
    for i in 1:length(seqParameters)
        frame = get_elements_by_tagname(runElements[i],"Frame")
        files = [get_elements_by_tagname(fr,"File") for fr=frame]
        files = [[attributes_dict(f) for f=fi] for fi=files]
        timing = [Meta.parse(attribute(fr,"absoluteTime")) for fr=frame]
        frames[i] = Dict("timing" => timing,"File" => files)
    end

    laserParameters = parseLaserConfiguration(scanConfigDoc)
    serieParameters = parseSerieParameters(scanConfigDoc,prairieRelease)

    ### Getting DAC related information
    dacParameters = [getDAC(runE,prairieRelease,baseDir) for runE=runElements]

    Dict("laser" => laserParameters,"globalConfig" => serieParameters,"sequences" => seqParameters,"version" => prairieRelease, "frames" => frames, "dataDir" => baseDir, "DAC" => dacParameters)
end

### Functions to access given elements of the xml doc ###########################
    ### Extract the laser settings (only useful if there's more than one laser being used)
    function parseLaserConfiguration(xmlPrairieDoc)
        systemConfigItem = find_element(xmlPrairieDoc,"SystemConfiguration")
        laserItem =  find_element(systemConfigItem,"Lasers")
        laserItems = get_elements_by_tagname(laserItem,"Laser")
        rigConfig = Dict("laser.$i" => attributes_dict(laserItems[i]) for i=1:length(laserItems))
    end

    ### Global PVStateShard configuration of the experiment
    function parseSerieParameters(xmlPrairieDoc,version,dictOut=true)
        version = Meta.parse(version[1:3])
          ### In Prairie 4, parameters are given for each frame, so we take the very first frame parameters for the  global ones (and just ignore if there's anything frame specific, like in ZSeries)
        if 4 <= version < 5
            params = find_element(find_element(find_element(xmlPrairieDoc,"Sequence"),"Frame"),"PVStateShard")
        elseif version >= 5
            params = find_element(xmlPrairieDoc,"PVStateShard")
        end
        if dictOut
            params = pvParameters2dict(params)
        end
        params
    end

### Functions parsing sequence elements
    ### Sequence configuration and attributes
    function dictSequenceParameters(seqElement)
        ### In some cases (?) there are StateShard elements per sequence
        Dict("attributes" => attributes_dict(seqElement), "state" => pvParameters2dict(find_element(seqElement,"PVStateShard")))
    end

### Other utilities
### Creates a dict from an xml "PVStateShard" node
function pvParameters2dict(prairieConfNode)
    if  prairieConfNode == nothing
        return(false)
    elseif !has_children(prairieConfNode)
        return(false)
    else
    keyElement = get_elements_by_tagname(prairieConfNode,"Key")
    ### false is written False in the xml file (hence the lowercase), everything else is numeric except the magnification that has an "x" written next to it (so we remove it to get the objective mag as a number).
    params = Dict(attribute(keyElement[i],"key") => Meta.parse(replace(lowercase(attribute(keyElement[i],"value")),"x","")) for i=1:length(keyElement))
    return params
    end
end

### Functions reading the output of an import
### Return image frames of a given sequence

function getPrairieFrames(prairieImport;seqN=1,channel=2,frameN="All")

    if frameN=="All"
        frames = prairieImport["frames"][seqN]["File"]
    else
        frames = prairieImport["frames"][seqN]["File"][[frameN]]
    end

    channels = [Meta.parse(fr["channel"]) for fr=frames[1]]
    whichFile = findfirst(channels,channel)
    filenames = [joinpath(prairieImport["dataDir"],fX[whichFile]["filename"]) for fX=frames]
    seqParams = prairieImport["globalConfig"]

    width = seqParams["pixelsPerLine"]
    height = seqParams["linesPerFrame"]

    im = pmap(load,filenames)
    #im = reinterpret(Normed{UInt16,16},cat(3,im...))
    im = reinterpret(N4f12,cat(3,im...))  ## Output from Prairie is 12 bits, in a 16 bits file.
    #protocolType = prairieImport["sequences"][seqN]["attributes"]["type"]

    xStep = seqParams["micronsPerPixel_XAxis"]
    yStep = seqParams["micronsPerPixel_YAxis"]
    timeStep = seqParams["framePeriod"]*seqParams["rastersPerFrame"]

    im = AxisArray(im,Axis{:x}(xStep*(1/2:size(im,1))),Axis{:y}(yStep*(1/2:size(im,2))),Axis{:time}(timeStep*(0:(size(im,3)-1))))
    #im["timedim"] = 3
    #im["spatialorder"]=["x","y"]
    #im["pixelspacing"] = [seqParams["micronsPerPixel_XAxis"],seqParams["micronsPerPixel_YAxis"],seqParams["framePeriod"]*seqParams["rastersPerFrame"]]
    #im["period"]=seqParams["framePeriod"]*seqParams["rastersPerFrame"]
    im
end

function getDAC(seqItem,version,baseDir)
    version = Meta.parse(version[1:3])
    if version>=5
        if find_element(seqItem,"VoltageOutput")!=nothing
            vOut = get_elements_by_tagname(seqItem,"VoltageOutput")
            voltOutPath = attribute(vOut[1],"filename")
            voltOutPath = "$(baseDir)/$(voltOutPath)"
            voutXml = root(parse_file(voltOutPath))
            voutCh = get_elements_by_tagname(voutXml,"Waveform")
            output = findall(map((x) -> content(get_elements_by_tagname(x,"Enabled")[1])=="true",voutCh))

            outList = map(voutCh[output]) do dac

            #dacList <- xmlToList(dac)
        #    dacList <- dacList[-which(grepl("text",names(dacList)) | grepl("PlotColor",names(dacList)))]
                pulseTrain = get_elements_by_tagname(dac,"WaveformComponent_PulseTrain")
                stimDict = Dict{String,Any}()
                for e in child_elements(pulseTrain[1])
                    val = Meta.parse(replace(content(e)," ",""))
                    if typeof(val) == Symbol
                        val = string(val)
                    end
                    stimDict[name(e)] = val
                end

                dacDict = Dict{String,Any}()
                for d in child_elements(dac)
                    if "PlotColor" != name(d) != "WaveformComponent_PulseTrain"
                        val = Meta.parse(replace(content(d)," ",""))
                        if typeof(val) == Symbol
                            val = string(val)
                        end
                        dacDict[name(d)] = val
                    end
                end

                dacDict["triggerMode"] = attribute(vOut[1],"triggerMode")
                dacDict["relativeTime"] = Meta.parse(attribute(vOut[1],"relativeTime"))
                dacDict["absoluteTime"] = Meta.parse(attribute(vOut[1],"absoluteTime"))
                Dict("DACDeviceParameters"=>dacDict, "DACStimulusParameters"=>stimDict)
            end
        else
            nothing
        end
    elseif 4<=version<5
       ##TO DO, for now just return the prm file
        [joinpath(baseDir,filter(x->contains(x,".prm"), readdir(baseDir))[1])]
    end

end


text:
# certain.jl
# ==========
#
# Define Certain sites for the site-counting framework.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

@inline function nibble_mask(::Type{Certain}, x::UInt64)
    return nibble_mask(enumerate_nibbles(x), 0x1111111111111111)
end

@inline function nibble_mask(::Type{Certain}, a::UInt64, b::UInt64)
    return nibble_mask(Certain, a) & nibble_mask(Certain, b)
end


text:
###ODE for a predator-prey model
### t = timepoints to evaluate the ODE at
### y = the model's state variables (prey population size - predator population size)
### ydot = the derivate values
### paras = the equation parameters (r,K,s,a,u,v)
function predatorpreyode(t,y,ydot,paras)

    r = paras[1]
    k = paras[2]
    s = paras[3]
    a = paras[4]
    u = paras[5]
    v = paras[6]

    tmp = y[1]*y[2]/(a + y[1])

    ydot[1] = r*y[1]*(1-y[1]/k) - s*tmp
    ydot[2] = u*tmp - v*y[2]

    nothing
end

###Convenience function to create a Predator-Prey ODEModel
### y0: initial conditions of the dynamic system
### variance: the variance of the noise model
### paraminit: the default and boundary values of the parameters
function predatorpreymodel(y0::AbstractVector,variance::AbstractVector,paraminit...)
    modelparameters = parameters([:r,:k,:s,:a,:u,:v],paraminit...)
    noisemodel = noise(:gaussian,variance)
    measurementdata = data(:array,predatorpreydataset1()...)
    model(:ode,modelparameters,measurementdata,noisemodel,predatorpreyode,y0,2,[1,2];
          name="Predator-Prey",abstol=10e-9,reltol=10e-9)
end