text:
# Simulations parameters
include("/home/kevin/Documents/1 - PhD Projects/1 - Sequential Decision making/3-SequentialEffects/1-Inputs/src/decision_making.jl")
include("/home/kevin/Documents/1 - PhD Projects/1 - Sequential Decision making/3-SequentialEffects/1-Inputs/src/reaction_time_analysis_functions.jl")
include("/home/kevin/Documents/1 - PhD Projects/1 - Sequential Decision making/3-SequentialEffects/1-Inputs/src/distributions.jl")

param = Dict("T_rest"=>10000,
"mu0"=>30,
"noise_amp"=>0.02,
"I0E1"=>0.3255,
"I0E2"=>0.3255,
"JN11"=>0.2609,
"JN22"=>0.2609,
"JN12"=>0.0497,
"JN21"=>0.0497,
"JAext" => 0.00052,
"Tnmda" => 100,
"Tampa" => 2,
"gamma" => 0.641,
"Tstim"=>5000,
"coh"=>5,
"false_time"=>300,
"Threshold"=>20
)

fixed_param=Dict("a"=>270,"b"=>108,"d"=>0.1540)
#Parameters of the simulation for EUler methods
simu_param=Dict("dt"=>0.5,"time_wind"=>4,"slide_wind"=>4) #4=2/dt

using DataFrames
using CSV

#a=2
Tr=500
println(Tr)
# de 0.02 à 0.025 déjàa sauvegardé (voir 0.02,0.3 50)

Irlist=[0.035,0.047]
Cohlist=[2,10,20]

using Plots
gr()

############################
## Fig neuronale PES
############################
a=300
PEsameL=zeros(251+a)#
PEdiffL=zeros(251+a)
PEsameW=zeros(251+a)
PEdiffW=zeros(251+a)
PCsameL=zeros(251+a)#
PCdiffL=zeros(251+a)
PCsameW=zeros(251+a)
PCdiffW=zeros(251+a)
numPEsame=0
numPEdiff=0
numPEsame=0
numPEdiff=0
numPCsame=0
numPCdiff=0
numPCsame=0
numPCdiff=0
Ir=Irlist[2]
c=Cohlist[3]

numr=0
numr2=1

for i=1:25
    r=CSV.read("./datatempAll/Coh$c-Final200Exp$i-T$Tr-Ir$Ir.csv")
    r2=CSV.read("./datatempAll/NondynCoh$c-Final200Exp$i-T$Tr-Ir$Ir.csv")

temp = 250
for i=1:2500
temp= temp + Int(r2[:RTs][i]/2)
if r2[:NotError][i] == 1 && r2[:Distri][i] >0 && r2[:Distri][i+1] >0
PCsameL=PCsameL + r[:R2][temp:(temp+250+a)]
PCsameW=PCsameW + r[:R1][temp:(temp+250+a)]

numPCsame=numPCsame+1
#rlist2=rlist2 + r[:R2][temp:(temp+250)]

elseif r2[:NotError][i] == -1 && r2[:Distri][i] <0 && r2[:Distri][i+1] >0
    #rlist=rlist + r[:R2][temp:(temp+250)]
    PEsameL=PEsameL + r[:R2][temp:(temp+250+a)]
    PEsameW=PEsameW + r[:R1][temp:(temp+250+a)]

    numPEsame=numPEsame+1

elseif r2[:NotError][i] == -1 && r2[:Distri][i] >0 && r2[:Distri][i+1] <0
    #rlist=rlist + r[:R2][temp:(temp+250)]
    PEsameL=PEsameL + r[:R1][temp:(temp+250+a)]
    PEsameW=PEsameW + r[:R2][temp:(temp+250+a)]

    numPEsame=numPEsame+1
elseif r2[:NotError][i] == 1 && r2[:Distri][i] <0 && r2[:Distri][i+1] <0
    #rlist=rlist + r[:R2][temp:(temp+250)]
    PCsameL=PCsameL + r[:R1][temp:(temp+250+a)]
    PCsameW=PCsameW + r[:R2][temp:(temp+250+a)]

    numPCsame=numPCsame+1

elseif r2[:NotError][i] == 1 && r2[:Distri][i] >0 && r2[:Distri][i+1] <0 && r2[:NotError][i+1] == 1
    PCdiffL=PCdiffL + r[:R2][temp:(temp+250+a)]
    PCdiffW=PCdiffW + r[:R1][temp:(temp+250+a)]

    numPCdiff=numPCdiff+1
    #rlist2=rlist2 + r[:R2][temp:(temp+250)]

elseif r2[:NotError][i] == -1 && r2[:Distri][i] <0 && r2[:Distri][i+1] <0 && r2[:NotError][i+1] == 1
        #rlist=rlist + r[:R2][temp:(temp+250)]
        PEdiffL=PEdiffL + r[:R2][temp:(temp+250+a)]
        PEdiffW=PEdiffW + r[:R1][temp:(temp+250+a)]

        numPEdiff=numPEdiff+1

    elseif r2[:NotError][i] == -1 && r2[:Distri][i] >0 && r2[:Distri][i+1] >0 && r2[:NotError][i+1] == 1
        #rlist=rlist + r[:R2][temp:(temp+250)]
        PEdiffL=PEdiffL + r[:R1][temp:(temp+250+a)]
        PEdiffW=PEdiffW + r[:R2][temp:(temp+250+a)]

        numPEdiff=numPEdiff+1
    elseif r2[:NotError][i] == 1 && r2[:Distri][i] <0 && r2[:Distri][i+1] >0 && r2[:NotError][i+1] == 1
        #rlist=rlist + r[:R2][temp:(temp+250)]
        PCdiffL=PCdiffL + r[:R1][temp:(temp+250+a)]
        PCdiffW=PCdiffW + r[:R2][temp:(temp+250+a)]

        numPCdiff=numPCdiff+1

end
temp = temp+250
end

end

PEsameL=PEsameL./numPEsame
PEdiffL=PEdiffL./numPEdiff
PEsameW=PEsameW./numPEsame
PEdiffW=PEdiffW./numPEdiff
PCsameL=PCsameL./numPCsame
PCdiffL=PCdiffL./numPCdiff
PCsameW=PCsameW./numPCsame
PCdiffW=PCdiffW./numPCdiff

plot(PEsameL,linewidth=3)
plot!(PCsameL,linewidth=3)

plot(PEdiffL[5:250],linewidth=4,color=:darkorange,line=:dash)
plot!(PCdiffL[5:250],linewidth=4,color=:blue,line=:dash)
plot!(PEdiffW[5:250],linewidth=4,color=:darkorange)
plot!(PCdiffW[5:250],linewidth=4,color=:blue)
plot!(legend=:bottomleft)

savefig("./DiffPEQTrue.svg")

mean(r2[:RTs])

plot!(PEsameW,linewidth=3)
plot!(PCsameW,linewidth=3)

plot(PEsameL[5:250],linewidth=4,color=:darkorange,line=:dash)
plot!(PCsameL[5:250],linewidth=4,color=:blue,line=:dash)
plot!(PEsameW[5:250],linewidth=4,color=:darkorange)
plot!(PCsameW[5:250],linewidth=4,color=:blue)
plot!(legend=:topleft)

savefig("./SamePEQ.svg")

plot(PEdiffL[30:320]-PCdiffL[30:320],linewidth=3,color=:lightgray)
plot!(-PEsameW[30:320]+PCsameW[30:320],linewidth=3,color=:black)
vline!([220],linewidth=4,color=:black,line=:dash)

savefig("./DifferencePES.svg")

############################
## Fig neuronale PES: time from sacccade
############################
a=80
PEsameL=zeros(a+1)#
PEdiffL=zeros(a+1)
PEsameW=zeros(a+1)
PEdiffW=zeros(a+1)
PCsameL=zeros(a+1)#
PCdiffL=zeros(a+1)
PCsameW=zeros(a+1)
PCdiffW=zeros(a+1)
numPEsame=0
numPEdiff=0
numPEsame=0
numPEdiff=0
numPCsame=0
numPCdiff=0
numPCsame=0
numPCdiff=0
Ir=Irlist[2]
c=Cohlist[1]

numr=0
numr2=1

for i=1:25
    r=CSV.read("./datatempAll/Coh$c-Final200Exp$i-T$Tr-Ir$Ir.csv")
    r2=CSV.read("./datatempAll/NondynCoh$c-Final200Exp$i-T$Tr-Ir$Ir.csv")

temp = 250
for i=1:2500
temp= temp + Int(r2[:RTs][i]/2) + 250
if r2[:NotError][i] == 1 && r2[:Distri][i] >0 && r2[:Distri][i+1] >0
PCsameL=PCsameL + r[:R2][temp:temp+a]
PCsameW=PCsameW + r[:R1][temp:temp+a]

numPCsame=numPCsame+1
#rlist2=rlist2 + r[:R2][temp:(temp+250)]

elseif r2[:NotError][i] == -1 && r2[:Distri][i] <0 && r2[:Distri][i+1] >0
    #rlist=rlist + r[:R2][temp:(temp+250)]
    PEsameL=PEsameL + r[:R2][temp:temp+a]
    PEsameW=PEsameW + r[:R1][temp:temp+a]

    numPEsame=numPEsame+1

elseif r2[:NotError][i] == -1 && r2[:Distri][i] >0 && r2[:Distri][i+1] <0
    #rlist=rlist + r[:R2][temp:(temp+250)]
    PEsameL=PEsameL + r[:R1][temp:temp+a]
    PEsameW=PEsameW + r[:R2][temp:temp+a]

    numPEsame=numPEsame+1
elseif r2[:NotError][i] == 1 && r2[:Distri][i] <0 && r2[:Distri][i+1] <0
    #rlist=rlist + r[:R2][temp:(temp+250)]
    PCsameL=PCsameL + r[:R1][temp:temp+a]
    PCsameW=PCsameW + r[:R2][temp:temp+a]

    numPCsame=numPCsame+1

elseif r2[:NotError][i] == 1 && r2[:Distri][i] >0 && r2[:Distri][i+1] <0 && r2[:NotError][i+1] == 1
    PCdiffL=PCdiffL + r[:R2][temp:temp+a]
    PCdiffW=PCdiffW + r[:R1][temp:temp+a]

    numPCdiff=numPCdiff+1
    #rlist2=rlist2 + r[:R2][temp:(temp+250)]

elseif r2[:NotError][i] == -1 && r2[:Distri][i] <0 && r2[:Distri][i+1] <0 && r2[:NotError][i+1] == 1
        #rlist=rlist + r[:R2][temp:(temp+250)]
        PEdiffL=PEdiffL + r[:R2][temp:temp+a]
        PEdiffW=PEdiffW + r[:R1][temp:temp+a]

        numPEdiff=numPEdiff+1

    elseif r2[:NotError][i] == -1 && r2[:Distri][i] >0 && r2[:Distri][i+1] >0 && r2[:NotError][i+1] == 1
        #rlist=rlist + r[:R2][temp:(temp+250)]
        PEdiffL=PEdiffL + r[:R1][temp:temp+a]
        PEdiffW=PEdiffW + r[:R2][temp:temp+a]

        numPEdiff=numPEdiff+1
    elseif r2[:NotError][i] == 1 && r2[:Distri][i] <0 && r2[:Distri][i+1] >0 && r2[:NotError][i+1] == 1
        #rlist=rlist + r[:R2][temp:(temp+250)]
        PCdiffL=PCdiffL + r[:R1][temp:temp+a]
        PCdiffW=PCdiffW + r[:R2][temp:temp+a]

        numPCdiff=numPCdiff+1

end

end

end

PEsameL=PEsameL./numPEsame
PEdiffL=PEdiffL./numPEdiff
PEsameW=PEsameW./numPEsame
PEdiffW=PEdiffW./numPEdiff
PCsameL=PCsameL./numPCsame
PCdiffL=PCdiffL./numPCdiff
PCsameW=PCsameW./numPCsame
PCdiffW=PCdiffW./numPCdiff

using Plots
gr()
plot(PEsameL,linewidth=3,color=:darkorange,line=:dash)
plot!(PCsameL,linewidth=3,color=:blue,line=:dash)

plot!(PEsameW,color=:darkorange,linewidth=3)
plot!(PCsameW,color=:blue,linewidth=3)
plot!(legend=:topleft)

#savefig("./PEQSameMotiononset.svg")
savefig("./NoSameMotiononset.svg")

plot(PEdiffL,linewidth=3,color=:darkorange,line=:dash)
plot!(PCdiffL,linewidth=3,color=:blue,line=:dash)
plot!(PEdiffW,color=:darkorange,linewidth=3)
plot!(PCdiffW,color=:blue,linewidth=3)
plot!(legend=:topleft)
savefig("./NoDiffTrueMotiononset.svg")

c_min=-30
 c_max = 30
 number = 40

   distri = make_distributions(c_min,c_max,number,20000)

 param["T_rest"] = 1500
 Ir=0.035
 If=0

  result,r2= sim_mf_decisionmaking_pattern(distri,Ir,If,
    param,fixed_param,simu_param)

Accrept=[]
Accaltt=[]

for a in linspace(-30,30,40)
    Accrep=Float64[]
    Accalt=Float64[]
for i=2:1000
    if r2[:Distri][i]==a
        if r2[:NotError][i-1]*sign(r2[:Distri][i-1])==1

            push!(Accrep,r2[:NotError][i]*sign(r2[:Distri][i]))
        else
            push!(Accalt,r2[:NotError][i]*sign(r2[:Distri][i]))
        end
    end
end

push!(Accrept,mean(Accrep))
push!(Accaltt,mean(Accalt))

end
r2

scatter(linspace(-30,30,40),0.5+0.5.*Accrept,color=:blue)
scatter!(linspace(-30,30,40),0.5+0.5.*Accaltt,color=:darkorange)

using LsqFit

model(x, p) = 1./(1+exp.(-(p[1]+p[2].*x)))

xdata = linspace(-30,30,40)
ydata = 0.5+0.5.*Accrept
p0 = [0.5, 0.5]
fit = curve_fit(model, xdata, ydata, p0)

fit.param

plot!(linspace(-30,30,40),1./(1+exp.(-(fit.param[1]+fit.param[2].*linspace(-30,30,40)))),color=:blue,label="Repetition",linewidth=3)

fit = curve_fit(model, xdata, 0.5+0.5.*Accaltt, p0)

fit.param

plot!(linspace(-30,30,40),1./(1+exp.(-(fit.param[1]+fit.param[2].*linspace(-30,30,40)))),color=:darkorange,label="Alternation",linewidth=3)
ylabel!("Percentage Rightward choice")
plot!(legend=:topleft)
savefig("./accuracy.pdf")

c_min=-30
 c_max = 30
 number = 20

   distri = make_distributions(c_min,c_max,number,10000)

 param["T_rest"] = 500
 Ir=0.04
 If=0

  result= sim_mf_decisionmaking_pattern_fast(distri,Ir,If,
    param,fixed_param,simu_param)

using LsqFit
r2=0

function list_elements(l)
      a = Float64[]
        i=1
        while (length(a) <20)
           if any(a.==l[i])

            else
                push!(a,l[i])
            end
            i=i+1
        end
        return sort(a)
    end

colorsT=[:red,:blue,:darkorange,:violet,:green]
o=1

param0=zeros(10)
paramR0=zeros(10)

paramR1=zeros(10)
param1=zeros(10)
fig=plot(legend=false)
for Tr in [200,500,1500,2500,5000]
    fig=plot()
    o=2
    Accrept=zeros(20,10)
    Accaltt=zeros(20,10)
for u=1:10
    r2 = CSV.read("/home/kevin/Documents/1 - PhD Projects/1 - Sequential Decision making/3-SequentialEffects/2-Data/AccuracySequential/Accuracyc51Tau200Exp$u-T$Tr-Ir0.04.csv")
j=0
for a in list_elements(r2[:Distri])
    j=j+1
    Accrep=Float64[]
    Accalt=Float64[]
for i=2:10000
    if r2[:Distri][i]==a
        if r2[:NotError][i-1]*sign(r2[:Distri][i-1])==1

            push!(Accrep,r2[:NotError][i]*sign(r2[:Distri][i]))
        else
            push!(Accalt,r2[:NotError][i]*sign(r2[:Distri][i]))
        end
    end
end

Accrept[j,u]=mean(Accrep)
Accaltt[j,u]=mean(Accalt)

end
end

scatter!(list_elements(r2[:Distri]),0.5+0.5.*mean(Accrept,2),yerr=sqrt.(var(Accrept,2)),color=colorsT[o])
scatter!(list_elements(r2[:Distri]),0.5+0.5.*mean(Accaltt,2),yerr=sqrt.(var(Accaltt,2)),color=colorsT[o])

model(x, p) = 1./(1+exp.(-(p[1]+p[2].*x)))
for u=1:10
xdata = linspace(-51,51,20)
ydata = 0.5+0.5.*Accrept[:,u]
p0 = [0.5, 0.5]
fit = curve_fit(model, xdata, ydata, p0)

paramR0[u]=fit.param[1]
paramR1[u]=fit.param[2]
plot!(linspace(-51,51,40),1./(1+exp.(-(fit.param[1]+fit.param[2].*linspace(-51,51,40)))),alpha=0.3,label="Repetition",linewidth=3,color=colorsT[o])
ydata = 0.5+0.5.*Accaltt[:,u]
p0 = [0.5, 0.5]
fit = curve_fit(model, xdata, ydata, p0)
plot!(linspace(-51,51,40),1./(1+exp.(-(fit.param[1]+fit.param[2].*linspace(-51,51,40)))),alpha=0.3,label="Repetition",linewidth=3,line=:dash,color=colorsT[o])
param0[u]=fit.param[1]
param1[u]=fit.param[2]
end
plot!(legend=false)
fig
savefig(fig,"./Acc$Tr.svg")
end

#### Left right indecision stuff
# compute param0/param1

histogram(+paramR0./paramR1-param0./param1,bins=10)
-param0./param1

scatter(linspace(-30,30,20),0.5+0.5.*Accrept,color=:blue)
scatter!(linspace(-30,30,20),0.5+0.5.*Accaltt,color=:darkorange)

using LsqFit

model(x, p) = 1./(1+exp.(-(p[1]+p[2].*x)))
for u=1:10
xdata = linspace(-20,20,12)
ydata = 0.5+0.5.*Accrept[:,u]
p0 = [0.5, 0.5]
fit = curve_fit(model, xdata, ydata, p0)

fit.param
plot!(linspace(-20,20,40),1./(1+exp.(-(fit.param[1]+fit.param[2].*linspace(-20,20,40)))),alpha=0.1,color=:blue,label="Repetition",linewidth=3)
end

fig

fit = curve_fit(model, xdata, 0.5+0.5.*Accaltt, p0)

fit.param

plot!(linspace(-30,30,40),1./(1+exp.(-(fit.param[1]+fit.param[2].*linspace(-30,30,40)))),alpha=0.5,color=:darkorange,label="Alternation",linewidth=3)
ylabel!("Percentage Rightward choice")
plot!(legend=:topleft)
savefig("./accuracy.pdf")

#######################################################
##########################################
#Mean stuff for neuronal activity
################################################
#################################################

############################
## Fig neuronale PES
############################
a=300
PEsameL=zeros(251+a,25)#
PEdiffL=zeros(251+a,25)
PEsameW=zeros(251+a,25)
PEdiffW=zeros(251+a,25)
PCsameL=zeros(251+a,25)#
PCdiffL=zeros(251+a,25)
PCsameW=zeros(251+a,25)
PCdiffW=zeros(251+a,25)

j=2
for u=2:25
numPEsame=0
numPEdiff=0
numPEsame=0
numPEdiff=0
numPCsame=0
numPCdiff=0
numPCsame=0
numPCdiff=0
Ir=Irlist[1]
c=Cohlist[2]

numr=0
numr2=1

    r=CSV.read("/home/kevin/Documents/1 - PhD Projects/1 - Sequential Decision making/3-SequentialEffects/2-Data/FiringRates/datatempAll/V$j-Coh$c-Final200Exp$u-T$Tr-Ir$Ir.csv")
    r2=CSV.read("/home/kevin/Documents/1 - PhD Projects/1 - Sequential Decision making/3-SequentialEffects/2-Data/FiringRates/datatempAll/V$j-NondynCoh$c-Final200Exp$u-T$Tr-Ir$Ir.csv")

temp = 250
for i=1:2990
temp= temp + Int(r2[:RTs][i]/2)
if r2[:NotError][i] == 1 && r2[:Distri][i] >0 && r2[:Distri][i+1] >0
PCsameL[:,u]=PCsameL[:,u] + r[:R2][temp:(temp+250+a)]
PCsameW[:,u]=PCsameW[:,u] + r[:R1][temp:(temp+250+a)]

numPCsame=numPCsame+1
#rlist2=rlist2 + r[:R2][temp:(temp+250)]

elseif r2[:NotError][i] == -1 && r2[:Distri][i] <0 && r2[:Distri][i+1] >0
    #rlist=rlist + r[:R2][temp:(temp+250)]
    PEsameL[:,u]=PEsameL[:,u] + r[:R2][temp:(temp+250+a)]
    PEsameW[:,u]=PEsameW[:,u] + r[:R1][temp:(temp+250+a)]

    numPEsame=numPEsame+1

elseif r2[:NotError][i] == -1 && r2[:Distri][i] >0 && r2[:Distri][i+1] <0
    #rlist=rlist + r[:R2][temp:(temp+250)]
    PEsameL[:,u]=PEsameL[:,u] + r[:R1][temp:(temp+250+a)]
    PEsameW[:,u]=PEsameW[:,u] + r[:R2][temp:(temp+250+a)]

    numPEsame=numPEsame+1
elseif r2[:NotError][i] == 1 && r2[:Distri][i] <0 && r2[:Distri][i+1] <0
    #rlist=rlist + r[:R2][temp:(temp+250)]
    PCsameL[:,u]=PCsameL[:,u] + r[:R1][temp:(temp+250+a)]
    PCsameW[:,u]=PCsameW[:,u] + r[:R2][temp:(temp+250+a)]

    numPCsame=numPCsame+1

elseif r2[:NotError][i] == 1 && r2[:Distri][i] >0 && r2[:Distri][i+1] <0 && r2[:NotError][i+1] == 1
    PCdiffL[:,u]=PCdiffL[:,u] + r[:R2][temp:(temp+250+a)]
    PCdiffW[:,u]=PCdiffW[:,u] + r[:R1][temp:(temp+250+a)]

    numPCdiff=numPCdiff+1
    #rlist2=rlist2 + r[:R2][temp:(temp+250)]

elseif r2[:NotError][i] == -1 && r2[:Distri][i] <0 && r2[:Distri][i+1] <0 && r2[:NotError][i+1] == 1
        #rlist=rlist + r[:R2][temp:(temp+250)]
        PEdiffL[:,u]=PEdiffL[:,u] + r[:R2][temp:(temp+250+a)]
        PEdiffW[:,u]=PEdiffW[:,u] + r[:R1][temp:(temp+250+a)]

        numPEdiff=numPEdiff+1

    elseif r2[:NotError][i] == -1 && r2[:Distri][i] >0 && r2[:Distri][i+1] >0 && r2[:NotError][i+1] == 1
        #rlist=rlist + r[:R2][temp:(temp+250)]
        PEdiffL[:,u]=PEdiffL[:,u] + r[:R1][temp:(temp+250+a)]
        PEdiffW[:,u]=PEdiffW[:,u] + r[:R2][temp:(temp+250+a)]

        numPEdiff=numPEdiff+1
    elseif r2[:NotError][i] == 1 && r2[:Distri][i] <0 && r2[:Distri][i+1] >0 && r2[:NotError][i+1] == 1
        #rlist=rlist + r[:R2][temp:(temp+250)]
        PCdiffL[:,u]=PCdiffL[:,u] + r[:R1][temp:(temp+250+a)]
        PCdiffW[:,u]=PCdiffW[:,u] + r[:R2][temp:(temp+250+a)]

        numPCdiff=numPCdiff+1

end
temp = temp+250
end

PEsameL[:,u]=PEsameL[:,u]./numPEsame
PEdiffL[:,u]=PEdiffL[:,u]./numPEdiff
PEsameW[:,u]=PEsameW[:,u]./numPEsame
PEdiffW[:,u]=PEdiffW[:,u]./numPEdiff
PCsameL[:,u]=PCsameL[:,u]./numPCsame
PCdiffL[:,u]=PCdiffL[:,u]./numPCdiff
PCsameW[:,u]=PCsameW[:,u]./numPCsame
PCdiffW[:,u]=PCdiffW[:,u]./numPCdiff

end

#############################"
##########" Same with S
#################

a=300
SPEsameL=zeros(251+a,25)#
SPEdiffL=zeros(251+a,25)
SPEsameW=zeros(251+a,25)
SPEdiffW=zeros(251+a,25)
SPCsameL=zeros(251+a,25)#
SPCdiffL=zeros(251+a,25)
SPCsameW=zeros(251+a,25)
SPCdiffW=zeros(251+a,25)

j=2
for u=2:25
numPEsame=0
numPEdiff=0
numPEsame=0
numPEdiff=0
numPCsame=0
numPCdiff=0
numPCsame=0
numPCdiff=0
Ir=Irlist[1]
c=Cohlist[2]

numr=0
numr2=1

    r=CSV.read("/home/kevin/Documents/1 - PhD Projects/1 - Sequential Decision making/3-SequentialEffects/2-Data/FiringRates/datatempAll/V$j-Coh$c-Final200Exp$u-T$Tr-Ir$Ir.csv")
    r2=CSV.read("/home/kevin/Documents/1 - PhD Projects/1 - Sequential Decision making/3-SequentialEffects/2-Data/FiringRates/datatempAll/V$j-NondynCoh$c-Final200Exp$u-T$Tr-Ir$Ir.csv")

temp = 250
for i=1:2990
temp= temp + Int(r2[:RTs][i]/2)
if r2[:NotError][i] == 1 && r2[:Distri][i] >0 && r2[:Distri][i+1] >0
SPCsameL[:,u]=SPCsameL[:,u] + r[:S2][temp:(temp+250+a)]
SPCsameW[:,u]=SPCsameW[:,u] + r[:S1][temp:(temp+250+a)]

numPCsame=numPCsame+1
#rlist2=rlist2 + r[:R2][temp:(temp+250)]

elseif r2[:NotError][i] == -1 && r2[:Distri][i] <0 && r2[:Distri][i+1] >0
    #rlist=rlist + r[:R2][temp:(temp+250)]
    SPEsameL[:,u]=SPEsameL[:,u] + r[:S2][temp:(temp+250+a)]
    SPEsameW[:,u]=SPEsameW[:,u] + r[:S1][temp:(temp+250+a)]

    numPEsame=numPEsame+1

elseif r2[:NotError][i] == -1 && r2[:Distri][i] >0 && r2[:Distri][i+1] <0
    #rlist=rlist + r[:R2][temp:(temp+250)]
    SPEsameL[:,u]=SPEsameL[:,u] + r[:S1][temp:(temp+250+a)]
    SPEsameW[:,u]=SPEsameW[:,u] + r[:S2][temp:(temp+250+a)]

    numPEsame=numPEsame+1
elseif r2[:NotError][i] == 1 && r2[:Distri][i] <0 && r2[:Distri][i+1] <0
    #rlist=rlist + r[:R2][temp:(temp+250)]
    SPCsameL[:,u]=SPCsameL[:,u] + r[:S1][temp:(temp+250+a)]
    SPCsameW[:,u]=SPCsameW[:,u] + r[:S2][temp:(temp+250+a)]

    numPCsame=numPCsame+1

elseif r2[:NotError][i] == 1 && r2[:Distri][i] >0 && r2[:Distri][i+1] <0 && r2[:NotError][i+1] == 1
    SPCdiffL[:,u]=SPCdiffL[:,u] + r[:S2][temp:(temp+250+a)]
    SPCdiffW[:,u]=SPCdiffW[:,u] + r[:S1][temp:(temp+250+a)]

    numPCdiff=numPCdiff+1
    #rlist2=rlist2 + r[:R2][temp:(temp+250)]

elseif r2[:NotError][i] == -1 && r2[:Distri][i] <0 && r2[:Distri][i+1] <0 && r2[:NotError][i+1] == 1
        #rlist=rlist + r[:R2][temp:(temp+250)]
        SPEdiffL[:,u]=SPEdiffL[:,u] + r[:S2][temp:(temp+250+a)]
        SPEdiffW[:,u]=SPEdiffW[:,u] + r[:S1][temp:(temp+250+a)]

        numPEdiff=numPEdiff+1

    elseif r2[:NotError][i] == -1 && r2[:Distri][i] >0 && r2[:Distri][i+1] >0 && r2[:NotError][i+1] == 1
        #rlist=rlist + r[:R2][temp:(temp+250)]
        SPEdiffL[:,u]=SPEdiffL[:,u] + r[:S1][temp:(temp+250+a)]
        SPEdiffW[:,u]=SPEdiffW[:,u] + r[:S2][temp:(temp+250+a)]

        numPEdiff=numPEdiff+1
    elseif r2[:NotError][i] == 1 && r2[:Distri][i] <0 && r2[:Distri][i+1] >0 && r2[:NotError][i+1] == 1
        #rlist=rlist + r[:R2][temp:(temp+250)]
        SPCdiffL[:,u]=SPCdiffL[:,u] + r[:S1][temp:(temp+250+a)]
        SPCdiffW[:,u]=SPCdiffW[:,u] + r[:S2][temp:(temp+250+a)]

        numPCdiff=numPCdiff+1

end
temp = temp+250
end

SPEsameL[:,u]=SPEsameL[:,u]./numPEsame
SPEdiffL[:,u]=SPEdiffL[:,u]./numPEdiff
SPEsameW[:,u]=SPEsameW[:,u]./numPEsame
SPEdiffW[:,u]=SPEdiffW[:,u]./numPEdiff
SPCsameL[:,u]=SPCsameL[:,u]./numPCsame
SPCdiffL[:,u]=SPCdiffL[:,u]./numPCdiff
SPCsameW[:,u]=SPCsameW[:,u]./numPCsame
SPCdiffW[:,u]=SPCdiffW[:,u]./numPCdiff

end

using Bootstrap
PEL = Float64[]
PELB=Float64[]
PELU = Float64[]

PCL = Float64[]
PCLB=Float64[]
PCLU = Float64[]

PEW = Float64[]
PEWB=Float64[]
PEWU = Float64[]

PCW = Float64[]
PCWB=Float64[]
PCWU = Float64[]

for i=5:300
    n_boot = 2000

      ## basic bootstrap
      #bs1 = bootstrap(PEsameL[i,2:end], mean, BasicSampling(n_boot))
      bs1 = bootstrap(PEdiffL[i,2:end], mean, BasicSampling(n_boot))

      cil = 0.95;

      ## basic CI
      bci1 = ci(bs1, BasicConfInt(cil));
      push!(PEL,bci1[1][1])
      push!(PELB,-bci1[1][3]+bci1[1][1])
      push!(PELU,bci1[1][2]-bci1[1][1])

     # bs1 = bootstrap(PCsameL[i,2:end], mean, BasicSampling(n_boot))
      bs1 = bootstrap(PCdiffL[i,2:end], mean, BasicSampling(n_boot))
      cil = 0.95;

      ## basic CI
      bci1 = ci(bs1, BasicConfInt(cil));
      push!(PCL,bci1[1][1])
      push!(PCLB,-bci1[1][3]+bci1[1][1])
      push!(PCLU,bci1[1][2]-bci1[1][1])

     # bs1 = bootstrap(PEsameW[i,2:end], mean, BasicSampling(n_boot))
     bs1 = bootstrap(PEdiffW[i,2:end], mean, BasicSampling(n_boot))

      cil = 0.95;

      ## basic CI
      bci1 = ci(bs1, BasicConfInt(cil));
      push!(PEW,bci1[1][1])
      push!(PEWB,-bci1[1][3]+bci1[1][1])
      push!(PEWU,bci1[1][2]-bci1[1][1])

    #  bs1 = bootstrap(PCsameW[i,2:end], mean, BasicSampling(n_boot))

      bs1 = bootstrap(PCdiffW[i,2:end], mean, BasicSampling(n_boot))
      cil = 0.95;

      ## basic CI
      bci1 = ci(bs1, BasicConfInt(cil));
      push!(PCW,bci1[1][1])
      push!(PCWB,-bci1[1][3]+bci1[1][1])
      push!(PCWU,bci1[1][2]-bci1[1][1])
end

fig=plot()
plot!(PCW,ribbon=(PCWB,PCWU),linewidth=1,color=:blue)
plot!(PEW,ribbon=(PEWB,PEWU),linewidth=1,color=:darkorange)

plot!(PEL,ribbon=(PELB,PELU),linewidth=1,line=:dash,color=:darkorange)
plot!(PCL,ribbon=(PCLB,PCLU),linewidth=1,line=:dash,color=:blue)

plot!(legend=false)

plot!(mean(PEsameL[5:300,2:end],2),line=:dash,ribbon=sqrt.(var(PEsameL[5:300,2:end],2)),color=:darkorange)
plot!(mean(PCsameL[5:300,2:end],2),line=:dash,ribbon=sqrt.(var(PCsameL[5:300,2:end],2)),color=:blue)

plot(linspace(5,300,296),mean(PEsameW[5:300,2:end],2),ribbon=sqrt.(var(PEsameW[5:300,2:end],2)),color=:darkorange,linewidth=3)
plot!(linspace(5,300,296),mean(PCsameW[5:300,2:end],2),ribbon=sqrt.(var(PCsameW[5:300,2:end],2)),color=:blue,linewidth=3)

plot!(legend=false)
savefig("./PESsame.svg")

plot(mean(log.(PEsameW[200:300,2:end]),2),ribbon=sqrt.(var(log.(PEsameW[200:300,2:end]),2)),color=:darkorange,linewidth=3)
plot!(mean(log.(PCsameW[200:300,2:end]),2),ribbon=sqrt.(var(log.(PCsameW[200:300,2:end]),2)),color=:blue,linewidth=3)

plot!(legend=false)

n_boot=1000
MSPE=Float64[]
MSPC=Float64[]
VUSPE=Float64[]
VDSPE=Float64[]
VUSPC=Float64[]
VDSPC=Float64[]
DMSPE=Float64[]
DVUSPE=Float64[]
DVDSPE=Float64[]

using Bootstrap
for i=10:350
    bs1 = bootstrap(SPEsameW[i,2:end], mean, BasicSampling(n_boot))

    cil = 0.95;

    ## basic CI
    bci1 = ci(bs1, BasicConfInt(cil));
    push!(MSPE,bci1[1][1])
    push!(VUSPE,-bci1[1][3]+bci1[1][1])
    push!(VDSPE,bci1[1][2]-bci1[1][1])

    bs1 = bootstrap(SPCsameW[i,2:end], mean, BasicSampling(n_boot))

    cil = 0.95;

    ## basic CI
    bci1 = ci(bs1, BasicConfInt(cil));
    push!(MSPC,bci1[1][1])
    push!(VUSPC,-bci1[1][3]+bci1[1][1])
    push!(VDSPC,bci1[1][2]-bci1[1][1])

    bs1 = bootstrap(SPEsameW[i,2:end]-SPCsameW[i,2:end], mean, BasicSampling(n_boot))

    cil = 0.95;

    ## basic CI
    bci1 = ci(bs1, BasicConfInt(cil));
    push!(DMSPE,bci1[1][1])
    push!(DVUSPE,-bci1[1][3]+bci1[1][1])
    push!(DVDSPE,bci1[1][2]-bci1[1][1])

end
MSPE
gr()

p1=plot(linspace(10,300,141),MSPE[1:141],ribbon=(VDSPE[1:141],VUSPE[1:141]),color=:darkorange,linewidth=2)
plot!(p1,linspace(10,300,141),MSPC[1:141],ribbon = (VDSPC[1:141],VUSPC[1:141]),color=:blue,linewidth=2)
ps=twinx(p1)
plot!(ps,linspace(600,700,51),MSPE[291:end],ribbon=(VDSPE[291:end],VUSPE[291:end]),color=:darkorange,linewidth=2)
plot!(ps,linspace(600,700,51),MSPC[291:end],ribbon=(VDSPC[291:end],VUSPC[291:end]),color=:blue,linewidth=2)
plot!(legend=false)
vline!([500],line=:dash,color=:black)
plot!(legend=false)

savefig("./PESsameS.svg")

plot(linspace(10,700,341),DMSPE./MSPE.*100,ribbon=(DVDSPE./MSPE.*100,DVUSPE./MSPE.*100),linewidth=3)

plot!(twinx(),linspace(10,700,341),DMSPE./MSPE+sqrt.(linspace(10,700,341)),ribbon=(DVDSPE./MSPE,DVUSPE./MSPE),linewidth=3)

plot!(legend=false)
savefig("./PESsameDiff.svg")

plot(mean(log.(SPEsameW[200:300,2:end]),2),ribbon=sqrt.(var(log.(SPEsameW[200:300,2:end]),2)),color=:darkorange,linewidth=3)
plot!(mean(log.(SPCsameW[200:300,2:end]),2),ribbon=sqrt.(var(log.(SPCsameW[200:300,2:end]),2)),color=:blue,linewidth=3)

plot!(legend=false)

plot(MSPE,ribbon=(VDSPE,VUSPC))
plot!(MSPC,ribbon = (VDSPC,VUSPC))

n_boot=1000
MSPE=Float64[]
MSPC=Float64[]
DMSPE=Float64[]
DVUSPE=Float64[]
DVDSPE=Float64[]

VUSPE=Float64[]
VDSPE=Float64[]
VUSPC=Float64[]
VDSPC=Float64[]
using Bootstrap
for i=10:350
    bs1 = bootstrap(PEsameW[i,2:end], mean, BasicSampling(n_boot))

    cil = 0.95;

    ## basic CI
    bci1 = ci(bs1, BasicConfInt(cil));
    push!(MSPE,bci1[1][1])
    push!(VUSPE,-bci1[1][3]+bci1[1][1])
    push!(VDSPE,bci1[1][2]-bci1[1][1])

    bs1 = bootstrap(PCsameW[i,2:end], mean, BasicSampling(n_boot))

    cil = 0.95;

    ## basic CI
    bci1 = ci(bs1, BasicConfInt(cil));
    push!(MSPC,bci1[1][1])
    push!(VUSPC,-bci1[1][3]+bci1[1][1])
    push!(VDSPC,bci1[1][2]-bci1[1][1])

    bs1 = bootstrap(PEsameW[i,2:end]-PCsameW[i,2:end], mean, BasicSampling(n_boot))

    cil = 0.95;

    ## basic CI
    bci1 = ci(bs1, BasicConfInt(cil));
    push!(DMSPE,bci1[1][1])
    push!(DVUSPE,-bci1[1][3]+bci1[1][1])
    push!(DVDSPE,bci1[1][2]-bci1[1][1])

end

p1=plot(linspace(10,300,141),MSPE[1:141],ribbon=(VDSPE[1:141],VUSPE[1:141]),color=:darkorange,linewidth=2)
plot!(p1,linspace(10,300,141),MSPC[1:141],ribbon = (VDSPC[1:141],VUSPC[1:141]),color=:blue,linewidth=2)
ps=twinx(p1)
plot!(ps,linspace(420,700,141),MSPE[201:end],ribbon=(VDSPE[201:end],VUSPE[201:end]),color=:darkorange,linewidth=2)
plot!(ps,linspace(420,700,141),MSPC[201:end],ribbon=(VDSPC[201:end],VUSPC[201:end]),color=:blue,linewidth=2)
vline!([500],line=:dash,color=:black)
plot!(legend=false)
savefig("./PESBootstrapsame.svg")

plot(linspace(10,700,341),MSPE,ribbon=(VDSPE,VUSPC),linewidth=3,color=:darkorange)
plot!(linspace(10,700,341),MSPC,ribbon = (VDSPC,VUSPC),linewidth=3,color=:blue)
plot!(legend=false)
savefig("./PESBootstrapsame.svg")

plot(linspace(10,700,341),DMSPE./MSPE.*100,ribbon=(DVDSPE./MSPE.*100,DVUSPE./MSPE.*100),linewidth=3)

plot(linspace(10,700,341),DMSPE,ribbon=(DVDSPE,DVUSPE),linewidth=3)
plot!(legend=false)
savefig("./DiffPESsame.svg")

savefig(fig,"./PEQdiffdec.png")

plot(mean(PEdiffL[5:300,2:end],2),ribbon=sqrt.(var(PEdiffL[5:300,2:end],2)))
plot!(mean(PCdiffL[5:300,2:end],2),ribbon=sqrt.(var(PCdiffL[5:300,2:end],2)))

plot!(legend=false)

plot!(mean(PEdiffW[5:300,2:end],2),ribbon=sqrt.(var(PEdiffW[5:300,2:end],2)))
plot!(mean(PCdiffW[5:300,2:end],2),ribbon=sqrt.(var(PCdiffW[5:300,2:end],2)))

plot!(legend=false)

mean(PEsameL[:,2:3],2)
var(PEsameL[:,2:3],2)

#### Indiv stuffs
i=1
r=CSV.read("/home/kevin/Documents/1 - PhD Projects/1 - Sequential Decision making/3-SequentialEffects/2-Data/FiringRates/datatempAll/Coh$c-Final200Exp$i-T$Tr-Ir$Ir.csv")
r2=CSV.read("/home/kevin/Documents/1 - PhD Projects/1 - Sequential Decision making/3-SequentialEffects/2-Data/FiringRates/datatempAll/NondynCoh$c-Final200Exp$i-T$Tr-Ir$Ir.csv")

for i=1:2300
if r2[:NotError][i]==-1
    print(i)
end
end

r2[:NotError][35]
r2[:Distri][35]
r2[:Distri][36]

r2[:NotError][270]
r2[:Distri][144]
r2[:Distri][145]
tinit = sum(r2[:RTs][1:357]./2)+250*357
fig=plot()
#plot(r[:R1][Int(tinit)-100:Int(tinit)+450],linewidth=2,color=:darkorange)
plot!(linspace(0,1000,501),r[:R2][Int(tinit)-50:Int(tinit)+450],linewidth=2,color=:darkorange)
r2[:NotError][103]
r2[:Distri][108]
r2[:Distri][119]
tinit = sum(r2[:RTs][1:118]./2)+250*118
plot!(linspace(0,1000,501),r[:R1][Int(tinit)-50:Int(tinit)+450],linewidth=2,color=:blue)
#plot!(r[:R2][Int(tinit)-100:Int(tinit)+450],linewidth=2,color=:blue)
plot!(legend=false)
vline!([100],line=:dash,linewidth=2,color=:black)
vline!([600],line=:dash,linewidth=2,color=:black)

fig
savefig(fig,"./PEQdiffFR.svg")

r2[:NotError][521]
r2[:Distri][521]
r2[:Distri][522]
tinit = sum(r2[:RTs][1:521]./2)+250*521
fig=plot()
#plot(r[:R1][Int(tinit)-100:Int(tinit)+450],linewidth=2,color=:darkorange)
plot!(linspace(0,900,451),r[:R2][Int(tinit)-50:Int(tinit)+400],linewidth=2,color=:darkorange)
r2[:NotError][103]
r2[:Distri][108]
r2[:Distri][111]
tinit = sum(r2[:RTs][1:111]./2)+250*111
plot!(linspace(0,900,451),r[:R2][Int(tinit)-50:Int(tinit)+400],linewidth=2,color=:blue)
#plot!(r[:R2][Int(tinit)-100:Int(tinit)+450],linewidth=2,color=:blue)
plot!(legend=false)
vline!([100],line=:dash,linewidth=2,color=:black)
vline!([600],line=:dash,linewidth=2,color=:black)

fig
savefig(fig,"./neeffFR.svg")

#### Indiv stuffs
i=1
c=Cohlist[2]
Ir=Irlist[1]
r=CSV.read("/home/kevin/Documents/1 - PhD Projects/1 - Sequential Decision making/3-SequentialEffects/2-Data/FiringRates/datatempAll/Coh$c-Final200Exp$i-T$Tr-Ir$Ir.csv")
r2=CSV.read("/home/kevin/Documents/1 - PhD Projects/1 - Sequential Decision making/3-SequentialEffects/2-Data/FiringRates/datatempAll/NondynCoh$c-Final200Exp$i-T$Tr-Ir$Ir.csv")

for i=1:2300
if r2[:NotError][i]==-1
    print(i)
end
end

r2[:NotError][53]
r2[:Distri][52]
r2[:Distri][53]

r2[:NotError][270]
r2[:Distri][144]
r2[:Distri][145]
fig=plot()
tinit = sum(r2[:RTs][1:499]./2)+250*499
#plot(r[:R1][Int(tinit)-100:Int(tinit)+450],linewidth=2,color=:darkorange)
#linspace(0,1000,501),
plot!(linspace(0,1000,501),r[:R2][Int(tinit)-50:Int(tinit)+450],linewidth=2,color=:darkorange)
r2[:NotError][161]
r2[:Distri][160]
r2[:Distri][161]
tinit = sum(r2[:RTs][1:365]./2)+250*365
plot!(linspace(0,1000,501),r[:R2][Int(tinit)-50:Int(tinit)+450],linewidth=2,color=:blue)
#plot!(r[:R2][Int(tinit)-100:Int(tinit)+450],linewidth=2,color=:blue)
plot!(legend=false)
vline!([100],line=:dash,linewidth=2,color=:black)
vline!([600],line=:dash,linewidth=2,color=:black)
savefig(fig,"./PESsameFR.svg")

for i=1:length(r2[:NotError])-1
    if r2[:NotError][i]==1 && r2[:NotError][i+1]==1 &&
         sign(r2[:Distri][i])==sign(r2[:Distri][i+1])
print(i)
print("\n")
end
end

for i=1:length(r2[:NotError])-1
    if r2[:NotError][i]==-1 && r2[:NotError][i+1]==1 &&
         sign(r2[:Distri][i])!=sign(r2[:Distri][i+1])
print(i)
print("\n")
end
end

fig

savefig(fig,"./PESdiffFR.svg")


text:
using TradeModels, MarketData, TradeBlotter

signal = (sma(AAPL["Close"],50) .> sma(AAPL["Close"],200));
#signal = (sma(Cl,50) .> sma(Cl,200));
#signal = Cl .> sma(Cl,200);
aapl    = FinancialTimeSeries(AAPL, "AAPL")
golden = Blotter(signal, aapl)

trades = tradearray(golden, aapl["Close"])

pnls = convert(Array{Float64}, ([(t.close-t.open) for t in trades].*100))
pos  = pnls[pnls.>0]
neg  = pnls[pnls.<0]
profitfactor = sum(pos)/abs(sum(neg))
mae = convert(Array{Float64}, ([t.open - minimum(t.timeseries["Low"].values) for t in trades]))
trades

# read from csv, parse into blotter and analyze

# qux = readcsv(Pkg.dir("FinancialBlotter/src/dev/account.csv"));
#
# fee = qux[:,[1,4,5,6,10]];
# fum = fee[2:11, :]
# eurtstamp = DateTime{ISOCalendar, UTC}[parsedatetime_from_TOS(t) for t in fum[:,1]]
# eurvals = convert(Array{Float64},fum[:,[2,5]])
# eurtrades = Blotter(eurtstamp, eurvals, blottercolnames, Ticker("EURUSD"))
#
# # get euro data
#
# fname = Pkg.dir("FinancialBlotter/src/dev/data.csv")
# blob    = readcsv(fname)
# tstamps  = Date{ISOCalendar}[date(i) for i in blob[:, 1]]
# vals  = convert(Array{Float64}, blob[:,2])
# eurta = TimeArray(tstamps, vals, ["Close"])

ext: jl


text:
# Benchmark suite for Molly
# Run with something like:
#   using Molly, PkgBenchmark
#   results = benchmarkpkg(Molly, BenchmarkConfig(env=Dict("JULIA_NUM_THREADS" => 16)))
#   export_markdown(out_file, results)

using Molly
using BenchmarkTools
using CUDA

using Base.Threads
using DelimitedFiles

run_parallel_tests = nthreads() > 1
if run_parallel_tests
    @info "The parallel benchmarks will be run as Julia is running on $(nthreads()) threads"
else
    @warn "The parallel benchmarks will not be run as Julia is running on 1 thread"
end

run_gpu_tests = CUDA.functional()
if run_gpu_tests
    @info "The GPU benchmarks will be run as a CUDA-enabled device is available"
else
    @warn "The GPU benchmarks will not be run as a CUDA-enabled device is not available"
end

CUDA.allowscalar(false) # Check that we never do scalar indexing on the GPU

const SUITE = BenchmarkGroup(
    [],
    "interactions" => BenchmarkGroup(),
    "spatial"      => BenchmarkGroup(),
    "simulation"   => BenchmarkGroup(),
    "protein"      => BenchmarkGroup(),
)

c1 = SVector(1.0, 1.0, 1.0)u"nm"
c2 = SVector(1.4, 1.0, 1.0)u"nm"
a1 = Atom(charge=1.0, σ=0.3u"nm", ϵ=0.2u"kJ * mol^-1")
box_size = SVector(2.0, 2.0, 2.0)u"nm"
coords = [c1, c2]
dr = vector(c1, c2, box_size)
b1 = HarmonicBond(b0=0.6u"nm", kb=100_000.0u"kJ * mol^-1 * nm^-2")

SUITE["interactions"]["LennardJones force" ] = @benchmarkable force($(LennardJones()), $(dr), $(c1), $(c2), $(a1), $(a1), $(box_size))
SUITE["interactions"]["LennardJones energy"] = @benchmarkable potential_energy($(LennardJones()), $(dr), $(c1), $(c2), $(a1), $(a1), $(box_size))
SUITE["interactions"]["Coulomb force"      ] = @benchmarkable force($(Coulomb()), $(dr), $(c1), $(c2), $(a1), $(a1), $(box_size))
SUITE["interactions"]["Coulomb energy"     ] = @benchmarkable potential_energy($(Coulomb()), $(dr), $(c1), $(c2), $(a1), $(a1), $(box_size))
SUITE["interactions"]["HarmonicBond force" ] = @benchmarkable force($(b1), $(c1), $(c2), $(box_size))
SUITE["interactions"]["HarmonicBond energy"] = @benchmarkable potential_energy($(b1), $(c1), $(c2), $(box_size))

SUITE["spatial"]["vector1D"] = @benchmarkable vector1D($(4.0u"nm"), $(6.0u"nm"), $(10.0u"nm"))
SUITE["spatial"]["vector"  ] = @benchmarkable vector($(SVector(4.0, 1.0, 1.0)u"nm"), $(SVector(6.0, 4.0, 3.0)u"nm"), $(SVector(10.0, 5.0, 3.5)u"nm"))

n_atoms = 400
atom_mass = 10.0u"u"
box_size = SVector(6.0, 6.0, 6.0)u"nm"
starting_coords = place_diatomics(n_atoms ÷ 2, box_size, 0.2u"nm", 0.2u"nm")
starting_velocities = [velocity(atom_mass, 1.0u"K") for i in 1:n_atoms]
starting_coords_f32 = [Float32.(c) for c in starting_coords]
starting_velocities_f32 = [Float32.(c) for c in starting_velocities]

function test_sim(nl::Bool, parallel::Bool, gpu_diff_safe::Bool, f32::Bool, gpu::Bool)
    n_atoms = 400
    n_steps = 200
    atom_mass = f32 ? 10.0f0u"u" : 10.0u"u"
    box_size = f32 ? SVector(6.0f0, 6.0f0, 6.0f0)u"nm" : SVector(6.0, 6.0, 6.0)u"nm"
    simulator = VelocityVerlet(dt=f32 ? 0.02f0u"ps" : 0.02u"ps")
    b0 = f32 ? 0.2f0u"nm" : 0.2u"nm"
    kb = f32 ? 10_000.0f0u"kJ * mol^-1 * nm^-2" : 10_000.0u"kJ * mol^-1 * nm^-2"
    bonds = [HarmonicBond(b0=b0, kb=kb) for i in 1:(n_atoms ÷ 2)]
    specific_inter_lists = (InteractionList2Atoms(collect(1:2:n_atoms), collect(2:2:n_atoms),
                            repeat([""], length(bonds)), gpu ? cu(bonds) : bonds),)

    neighbor_finder = NoNeighborFinder()
    cutoff = DistanceCutoff(f32 ? 1.0f0u"nm" : 1.0u"nm")
    pairwise_inters = (LennardJones(nl_only=false, cutoff=cutoff),)
    if nl
        if gpu_diff_safe
            neighbor_finder = DistanceVecNeighborFinder(nb_matrix=gpu ? cu(trues(n_atoms, n_atoms)) : trues(n_atoms, n_atoms),
                                                        n_steps=10, dist_cutoff=f32 ? 1.5f0u"nm" : 1.5u"nm")
        else
            neighbor_finder = DistanceNeighborFinder(nb_matrix=trues(n_atoms, n_atoms), n_steps=10,
                                                        dist_cutoff=f32 ? 1.5f0u"nm" : 1.5u"nm")
        end
        pairwise_inters = (LennardJones(nl_only=true, cutoff=cutoff),)
    end

    if gpu
        coords = cu(deepcopy(f32 ? starting_coords_f32 : starting_coords))
        velocities = cu(deepcopy(f32 ? starting_velocities_f32 : starting_velocities))
        atoms = cu([Atom(charge=f32 ? 0.0f0 : 0.0, mass=atom_mass, σ=f32 ? 0.2f0u"nm" : 0.2u"nm",
                            ϵ=f32 ? 0.2f0u"kJ * mol^-1" : 0.2u"kJ * mol^-1") for i in 1:n_atoms])
    else
        coords = deepcopy(f32 ? starting_coords_f32 : starting_coords)
        velocities = deepcopy(f32 ? starting_velocities_f32 : starting_velocities)
        atoms = [Atom(charge=f32 ? 0.0f0 : 0.0, mass=atom_mass, σ=f32 ? 0.2f0u"nm" : 0.2u"nm",
                        ϵ=f32 ? 0.2f0u"kJ * mol^-1" : 0.2u"kJ * mol^-1") for i in 1:n_atoms]
    end

    s = System(
        atoms=atoms,
        pairwise_inters=pairwise_inters,
        specific_inter_lists=specific_inter_lists,
        coords=coords,
        velocities=velocities,
        box_size=box_size,
        neighbor_finder=neighbor_finder,
        gpu_diff_safe=gpu_diff_safe,
    )

    simulate!(s, simulator, n_steps; parallel=parallel)
    return s.coords
end

runs = [
    ("in-place"        , [false, false, false, false, false]),
    ("in-place NL"     , [true , false, false, false, false]),
    ("in-place f32"    , [false, false, false, true , false]),
    ("out-of-place"    , [false, false, true , false, false]),
    ("out-of-place NL" , [true , false, true , false, false]),
    ("out-of-place f32", [false, false, true , true , false]),
]
if run_parallel_tests
    push!(runs, ("in-place parallel"   , [false, true , false, false, false]))
    push!(runs, ("in-place NL parallel", [true , true , false, false, false]))
end
if run_gpu_tests
    push!(runs, ("out-of-place gpu"       , [false, false, true , false, true ]))
    push!(runs, ("out-of-place gpu f32"   , [false, false, true , true , true ]))
    push!(runs, ("out-of-place gpu NL"    , [true , false, true , false, true ]))
    push!(runs, ("out-of-place gpu f32 NL", [true , false, true , true , true ]))
end

for (name, args) in runs
    test_sim(args...) # Run once for setup
    SUITE["simulation"][name] = @benchmarkable test_sim($(args[1]), $(args[2]), $(args[3]), $(args[4]), $(args[5]))
end

data_dir = normpath(@__DIR__, "..", "data")
ff_dir = joinpath(data_dir, "force_fields")
openmm_dir = joinpath(data_dir, "openmm_6mrr")

ff = OpenMMForceField(joinpath.(ff_dir, ["ff99SBildn.xml", "tip3p_standard.xml", "his.xml"])...)
velocities = SVector{3}.(eachrow(readdlm(joinpath(openmm_dir, "velocities_300K.txt"))))u"nm * ps^-1"
s = System(joinpath(data_dir, "6mrr_equil.pdb"), ff; velocities=velocities)
simulator = VelocityVerlet(dt=0.0005u"ps")
n_steps = 25

simulate!(s, simulator, n_steps; parallel=true)
SUITE["protein"]["in-place NL parallel"] = @benchmarkable simulate!($(s), $(simulator), $(n_steps); parallel=true)


text:
using WCS
using Unitful

export AbstractSpectrum
#--------------------------------------------------------------------------------------

"""
    AbstractSpectrum{W<:Number, F<:Number}

An abstract holder for astronomical spectra. All types inheriting from this must have the following fields

- wave::Array{W, N}
- flux::Array{F, N}
- meta::Dict{Symbol, Any}
"""
abstract type AbstractSpectrum{W,F} end

function Base.getproperty(spec::AbstractSpectrum, nm::Symbol)
    if nm in keys(getfield(spec, :meta))
        return getfield(spec, :meta)[nm]
    else
        return getfield(spec, nm)
    end
end

function Base.propertynames(spec::AbstractSpectrum)
    natural = (:wave, :flux, :meta)
    meta = keys(spec.meta)
    return (natural..., meta...)
end

"""
    wave(::AbstractSpectrum)

Return the wavelengths of the spectrum.
"""
wave(spec::AbstractSpectrum) = spec.wave

"""
    flux(::AbstractSpectrum)

Return the flux of the spectrum.
"""
flux(spec::AbstractSpectrum) = spec.flux

# Collection
Base.eltype(spec::AbstractSpectrum) = eltype(spec.flux)
Base.size(spec::AbstractSpectrum) = size(spec.flux)
Base.size(spec::AbstractSpectrum, i) = size(spec.flux, i)
Base.length(spec::AbstractSpectrum) = length(spec.flux)
Base.maximum(spec::AbstractSpectrum) = maximum(spec.flux)
Base.minimum(spec::AbstractSpectrum) = minimum(spec.flux)
Base.argmax(spec::AbstractSpectrum) = argmax(spec.flux)
Base.argmin(spec::AbstractSpectrum) = argmin(spec.flux)
Base.findmax(spec::AbstractSpectrum) = findmax(spec.flux)
Base.findmin(spec::AbstractSpectrum) = findmin(spec.flux)

# Arithmetic
Base.:+(s::T, A) where {T <: AbstractSpectrum} = T(s.wave, s.flux .+ A, s.meta)
Base.:*(s::T, A) where {T <: AbstractSpectrum} = T(s.wave, s.flux .* A, s.meta)
Base.:/(s::T, A) where {T <: AbstractSpectrum} = T(s.wave, s.flux ./ A, s.meta)
Base.:-(s::T) where {T <: AbstractSpectrum} = T(s.wave, -s.flux, s.meta)
Base.:-(s::AbstractSpectrum, A) = s + -A
Base.:-(A, s::AbstractSpectrum) = s - A

# Multi-Spectrum
Base.:+(s::T, o::T) where {T <: AbstractSpectrum} = T(s.wave, s.flux .+ o.flux, s.meta)
Base.:*(s::T, o::T) where {T <: AbstractSpectrum} = T(s.wave, s.flux .* o.flux, s.meta)
Base.:/(s::T, o::T) where {T <: AbstractSpectrum} = T(s.wave, s.flux ./ o.flux * unit(s)[2], s.meta)
Base.:-(s::T, o::T) where {T <: AbstractSpectrum} = T(s.wave, s.flux .- o.flux, s.meta)

"""
    Unitful.ustrip(::AbstractSpectrum)

Remove the units from a spectrum. Useful for processing spectra in tools that don't play nicely with `Unitful.jl`

# Examples
```jldoctest
julia> using Unitful, UnitfulAstro

julia> wave = range(1e4, 3e4, length=1000);

julia> flux = wave .* 10 .+ randn(1000);

julia> spec = spectrum(wave*u"angstrom", flux*u"W/m^2/angstrom")
Spectrum(Quantity{Float64,𝐋,Unitful.FreeUnits{(Å,),𝐋,nothing}}, Quantity{Float64,𝐌*𝐋^-1*𝐓^-3,Unitful.FreeUnits{(Å^-1, m^-2, W),𝐌*𝐋^-1*𝐓^-3,nothing}})

julia> ustrip(spec)
Spectrum(Float64, Float64)

```
"""
Unitful.ustrip(spec::AbstractSpectrum) = spectrum(ustrip.(spec.wave), ustrip.(spec.flux); spec.meta...)

"""
    Unitful.unit(::AbstractSpectrum)

Get the units of a spectrum. Returns a tuple of the wavelength units and flux/sigma units

# Examples
```jldoctest
julia> using Unitful, UnitfulAstro

julia> wave = range(1e4, 3e4, length=1000);

julia> flux = wave .* 10 .+ randn(1000);

julia> spec = spectrum(wave * u"angstrom", flux * u"W/m^2/angstrom");

julia> w_unit, f_unit = unit(spec)
(Å, W Å^-1 m^-2)

```
"""
Unitful.unit(spec::AbstractSpectrum) = unit(eltype(spec.wave)), unit(eltype(spec.flux))


text:
import Base: *

export adjacency_matrix,
    laplacian_matrix,
    incidence_matrix,
    coo_sparse,
    spectral_distance

"""
    adjacency_matrix(g, dir=:out, T::DataType=Int)

Returns a sparse boolean adjacency matrix for a graph, indexed by `[u, v]`
vertices. `true` values indicate an edge between `u` and `v`. Users may
specify a direction (`:in`, `:out`, or `:all` are currently supported; `:out`
is default for both directed and undirected graphs) and a data type for the
matrix (defaults to `Int`).
"""
function adjacency_matrix(g::AGraphOrDiGraph, dir::Symbol=:out, T::DataType=Int)
    n_v = nv(g)
    nz = ne(g) * (is_directed(g) ? 1 : 2)
    colpt = ones(Int, n_v + 1)

    # see below - we iterate over columns. That's why we take the
    # "opposite" neighbor function. It's faster than taking the transpose
    # at the end.
    if dir == :out
        neighborfn = in_neighbors
    elseif dir == :in
        neighborfn = out_neighbors
    elseif dir == :all
        if is_directed(g)
            #TODO use flatten and check performance
            neighborfn = (h,j)->[collect(in_neighbors(h, j)); collect(out_neighbors(h, j))]
            nz *= 2
        else
            neighborfn = out_neighbors
        end
    else
        error("Wrong dir ($dir)")
    end
    rowval = sizehint!(Vector{Int}(), nz)
    selfloops = Vector{Int}()
    for j in 1:n_v  # this is by column, not by row.
        if has_edge(g,j,j)
            push!(selfloops, j)
        end
        dsts = collect(neighborfn(g, j))
        colpt[j+1] = colpt[j] + length(dsts)
        append!(rowval, sort!(dsts))
    end
    spmx = SparseMatrixCSC(n_v, n_v, colpt, rowval, ones(T,nz))

    # this is inefficient. There should be a better way of doing this.
    # the issue is that adjacency matrix entries for self-loops are 2,
    # not one(T).
    for i in selfloops
        if !(T <: Bool)
            spmx[i,i] += one(T)
        end
    end
    return spmx
end

adjacency_matrix(g::AGraph, T::DataType=Int) = adjacency_matrix(g, :out, T)

"""
    laplacian_matrix(g, dir::Symbol=:out, T::DataType=Int)

Returns a sparse [Laplacian matrix](https://en.wikipedia.org/wiki/Laplacian_matrix)
for a graph `g`, indexed by `[u, v]` vertices. `dir` has to be `:in, :out` or `:all`.
"""
function laplacian_matrix(g::AGraphOrDiGraph, dir::Symbol=:out, T::DataType=Int)
    A = adjacency_matrix(g, dir, T)
    D = convert(SparseMatrixCSC{T}, Diagonal(sparse(sum(A, dims=2)[:])))
    return D - A
end

"""
    incidence_matrix(g::AGraphOrDiGraph, T::DataType=Int; oriented=false)

Returns a sparse node-arc incidence matrix for a graph, indexed by
`[v, i]`, where `i` is in `1:ne(g)`, indexing an edge `e`. For
directed graphs, a value of `-1` indicates that `src(e) == v`, while a
value of `1` indicates that `dst(e) == v`. Otherwise, the value is
`0`.

For undirected graphs, both entries are `1` if `oriented=false`,
otherwise `[v, i] -> -1` and `[u, i] -> 1` if `v < u`.
"""
function incidence_matrix(g::AGraphOrDiGraph, T::DataType=Int; oriented::Bool = false)
    isdir = is_directed(g)
    n_v = nv(g)
    n_e = ne(g)
    nz = 2 * n_e

    # every col has the same 2 entries
    colpt = collect(1:2:(nz + 1))
    nzval = repeat([(isdir || oriented) ? -one(T) : one(T), one(T)], n_e)

    # iterate over edges for row indices
    rowval = zeros(Int, nz)
    i = 1
    for u in vertices(g)
        for v in out_neighbors(g, u)
            if isdir || u < v # add every edge only once
                rowval[2*i - 1] = u
                rowval[2*i] = v
                i += 1
            end
        end
    end

    spmx = SparseMatrixCSC(n_v, n_e, colpt, rowval, nzval)
    return spmx
end

"""
    spectral_distance(G₁, G₂ [, k])

Compute the spectral distance between undirected n-vertex
graphs G₁ and G₂ using the top k ≤ n greatest eigenvalues.
If k is ommitted, uses full spectrum.

For further details, please refer to:

JOVANOVIC, I.; STANIC, Z., 2014. Spectral Distances of
Graphs Based on their Different Matrix Representations
"""
function spectral_distance(G₁::AGraph, G₂::AGraph, k::Integer)
  A₁ = adjacency_matrix(G₁)
  A₂ = adjacency_matrix(G₂)

  λ₁ = k < nv(G₁)-1 ? Arpack.eigs(A₁, nev=k, which=:LR)[1] : eigvals(Matrix(A₁))[end:-1:end-(k-1)]
  λ₂ = k < nv(G₂)-1 ? Arpack.eigs(A₂, nev=k, which=:LR)[1] : eigvals(Matrix(A₂))[end:-1:end-(k-1)]

  sum(abs, λ₁ - λ₂)
end

function spectral_distance(G₁::AGraph, G₂::AGraph)
  @assert nv(G₁) == nv(G₂) "spectral distance not defined for |G₁| != |G₂|"
  spectral_distance(G₁, G₂, nv(G₁))
end


text:
module GalacticEvolutionary

using Reexport, GalacticOptim, GalacticOptim.SciMLBase
@reexport using Evolutionary

decompose_trace(trace::Evolutionary.OptimizationTrace) = last(trace)
decompose_trace(trace::Evolutionary.OptimizationTraceRecord) = trace

function Evolutionary.trace!(record::Dict{String,Any}, objfun, state, population, method::Evolutionary.AbstractOptimizer, options)
    record["x"] = population
end

function __map_optimizer_args(prob::OptimizationProblem, opt::Evolutionary.AbstractOptimizer;
    callback=nothing,
    maxiters::Union{Number,Nothing}=nothing,
    maxtime::Union{Number,Nothing}=nothing,
    abstol::Union{Number,Nothing}=nothing,
    reltol::Union{Number,Nothing}=nothing,
    kwargs...)

    mapped_args = (;)

    mapped_args = (; mapped_args..., kwargs...)

    if !isnothing(callback)
        mapped_args = (; mapped_args..., callback=callback)
    end

    if !isnothing(maxiters)
        mapped_args = (; mapped_args..., iterations=maxiters)
    end

    if !isnothing(maxtime)
        mapped_args = (; mapped_args..., time_limit=maxtime)
    end

    if !isnothing(abstol)
        mapped_args = (; mapped_args..., abstol=abstol)
    end

    if !isnothing(reltol)
        mapped_args = (; mapped_args..., reltol=reltol)
    end

    return Evolutionary.Options(; mapped_args...)
end

function SciMLBase.__solve(prob::OptimizationProblem, opt::Evolutionary.AbstractOptimizer, data=GalacticOptim.DEFAULT_DATA;
    callback=(args...) -> (false),
    maxiters::Union{Number,Nothing}=nothing,
    maxtime::Union{Number,Nothing}=nothing,
    abstol::Union{Number,Nothing}=nothing,
    reltol::Union{Number,Nothing}=nothing,
    progress=false, kwargs...)
    local x, cur, state

    if data != GalacticOptim.DEFAULT_DATA
        maxiters = length(data)
    end

    cur, state = iterate(data)

    function _cb(trace)
        cb_call = callback(decompose_trace(trace).metadata["x"], trace.value...)
        if !(typeof(cb_call) <: Bool)
            error("The callback should return a boolean `halt` for whether to stop the optimization process.")
        end
        cur, state = iterate(data, state)
        cb_call
    end

    maxiters = GalacticOptim._check_and_convert_maxiters(maxiters)
    maxtime = GalacticOptim._check_and_convert_maxtime(maxtime)

    _loss = function (θ)
        x = prob.f(θ, prob.p, cur...)
        return first(x)
    end

    opt_args = __map_optimizer_args(prob, opt, callback=_cb, maxiters=maxiters, maxtime=maxtime, abstol=abstol, reltol=reltol; kwargs...)

    t0 = time()
    if isnothing(prob.ub) | isnothing(prob.ub)
        opt_res = Evolutionary.optimize(_loss, prob.u0, opt, opt_args)
    else
        cons = Evolutionary.BoxConstraints(prob.lb, prob.ub)
        opt_res = Evolutionary.optimize(_loss, cons, prob.u0, opt, opt_args)
    end
    t1 = time()
    opt_ret = Symbol(Evolutionary.converged(opt_res))

    SciMLBase.build_solution(prob, opt, Evolutionary.minimizer(opt_res), Evolutionary.minimum(opt_res); original=opt_res, retcode=opt_ret)
end

end


text:
logic = [
    Axiom("Definition of And", P ∧ FALSE, FALSE),
    Axiom("Definition of And", FALSE ∧ P, FALSE),
    Axiom("Definition of And", P ∧ TRUE, P),
    Axiom("Definition of And", TRUE ∧ P, P),
    Axiom("And is associative", @left_associative(∧, P, Q, R)...),
    Axiom("And is associative", @right_associative(∧, P, Q, R)...),
    Axiom("And is commutative", @commutative(∧, P, Q)...),
    Axiom("And is idempotent", @idempotent(∧, P)...),
    Axiom("De Morgan's laws", P ∧ Q, ~(~P ∨ ~Q)),
    Axiom("De Morgan's laws", ~(~P ∨ ~Q), P ∧ Q),
    Axiom("And is absorbing over Or", P ∧ (P ∨ Q), Q),
    Axiom("And is distributive over Or", @distributive(∧, ∨, P, Q, R)...),
    Axiom("And is factorisable over Or", @factorisable(∧, ∨, P, Q, R)...),

    Axiom("Law of noncontradiction", P ∧ ~P, FALSE),
    # now prove Q ∨ ~Q is true which is the law of excluded middle

    Axiom("Definition of Or", P ∨ TRUE, TRUE),
    Axiom("Definition of Or", TRUE ∨ P, TRUE),
    Axiom("Definition of Or", P ∨ FALSE, P),
    Axiom("Definition of Or", FALSE ∨ P, P),
    Axiom("Or is associative", @left_associative(∨, P, Q, R)...),
    Axiom("Or is associative", @right_associative(∨, P, Q, R)...),
    Axiom("Or is commutative", @commutative(∨, P, Q)...),
    Axiom("Or is idempotent", @idempotent(∨, P)...),
    Axiom("De Morgan's laws", P ∨ Q, ~(~P ∧ ~Q)),
    Axiom("De Morgan's laws", ~(~P ∧ ~Q), P ∨ Q),
    Axiom("Or is absorbing over And", P ∨ (P ∧ Q), Q),
    Axiom("Or is distributive over And", @distributive(∨, ∧, P, Q, R)...),
    Axiom("Or is factorisable over And", @factorisable(∨, ∧, P, Q, R)...),
    Axiom("Or is the weaker And", P ∧ Q, P ∨ Q),

    Axiom("Definition of Not", ~FALSE, TRUE),
    Axiom("Definition of Not", ~TRUE, FALSE),
    Axiom("De Morgan's laws", ~(P ∨ Q), ~P ∧ ~Q),
    Axiom("De Morgan's laws", ~(P ∨ Q), ~P ∧ ~Q),
    Axiom("Not is involute", @involute(~, P)...),
    Axiom("Not is involute", @re_involute(~, P)...),

    Axiom("Definition of Implies", P ⟹ Q, ~P ∨ Q),
    Axiom("Definition of Implies", ~P ∨ Q, P ⟹ Q),
    Theorem("Contrapositive", P ⟹ Q, ~Q ⟹ ~P, Proof([])),

    Axiom("Definition of Equivalent", P ≡ Q, (P ⟹ Q) ∧ (Q ⟹ P)),
    Axiom("Definition of Equivalent", (P ⟹ Q) ∧ (Q ⟹ P), P ≡ Q),
]


text:
using Dash
using DashHtmlComponents
using DashCoreComponents

app = dash()

app.layout = html_div() do
    dcc_input(id="input", value = "initial value"),
    html_div(id="output1"),
    html_div(id="output2")
end
callback!(app, callid"input.value => output1.children, output2.children") do input
    return ("$input first", "$input second")
end

run_server(app)


text:
module kMeansClustering

function squared_euclidean_distance(v::Array{Float64},w::Array{Float64})
    dist = 0
    for d in 1:length(v)
        dist+= (v[d]-w[d])^2
    end
    return dist
end

function agroup(matrix::Matrix{Float64})
    l = [Float64[]]
    pop!(l)
    dim,n_points = size(matrix)
    for p in 1:n_points
        v = Float64[]
        for d in 1:dim
            push!(v,matrix[dim*(p-1)+d])
        end
        push!(l,v)
    end
    return l,dim
end

function compare_distance(p::Array{Float64,1},old_ctr::Array{Float64,1},new_ctr::Array{Float64,1})
    old_dist = squared_euclidean_distance(p,old_ctr)
    new_dist = squared_euclidean_distance(p,new_ctr)
    if new_dist<=old_dist
        return new_ctr
    end
    return old_ctr
end

function dic(l::Array{Array{Float64,1},1},centers::Array{Array{Float64,1},1})
    dic_centers = Dict()
    for ctr in centers
        dic_centers[ctr] = [Float64[]]
        pop!(dic_centers[ctr])
    end
    for p in l
        center = centers[1]
        for ctr in centers
            center = compare_distance(p,center,ctr)
        end
        c = dic_centers[center]
        push!(c,p)
        dic_centers[center] = c
    end
    return dic_centers
end

function reorganize(clusters::Array{Array{Array{Float64,1},1},1},dim::Int64)
    centers = [Float64[]]
    pop!(centers)
    for cluster in clusters
        s = zeros(dim)
        for p in cluster
            for n in 1:dim
                s[n]+= p[n]
            end
        end
        means = Float64[]
        for d in s
            push!(means,d/length(cluster))
        end
        push!(centers,means)
    end
    return centers
end

function cost(dic_centers::Dict{Any,Any})
    cst = 0
    for (centr,clus) in dic_centers
        for p in clus
            cst+=squared_euclidean_distance(p,centr)
        end
    end
    return cst
end

function assign(clusters::Array{Array{Array{Float64,1},1},1},l::Array{Array{Float64,1},1})
    dic_assign = Dict()
    for n in 1:length(clusters)
        for el in clusters[n]
            dic_assign[el] = n
        end
    end
    l_assign = Int64[]
    for p in l
        push!(l_assign,dic_assign[p])
    end
    return l_assign
end

function new_format(centers::Array{Array{Float64,1},1})
    new = hcat(centers[1],centers[2])
    for p in centers[3:length(centers)]
        new = hcat(new,p)
    end
    return new
end

function build_clusters(dic_centers::Dict{Any,Any})
    clusters = [[Float64[]]]
    pop!(clusters)
    for (center,cluster) in dic_centers
        push!(clusters,cluster)
    end
    return clusters
end

function mykmeansclustering(matrix::Matrix{Float64},k::Int64)
    l,dim = agroup(matrix)
    centers = [Float64[]]
    pop!(centers)
    while length(unique(centers))<k
        centers = rand(l,k) # Forgy Method
    end
    initials = centers
    cop = [Float64[]]
    pop!(cop)
    while centers!=cop
        dic_centers = dic(l,centers)
        cop = centers[1:k]
        clusters = build_clusters(dic_centers)
        centers = reorganize(clusters,dim)
    end
    dic_centers = dic(l,centers)
    clusters = build_clusters(dic_centers)
    cst = cost(dic_centers)
    assignment = assign(clusters,l)
    centers = new_format(centers)
    initials = new_format(initials)
    return centers,cst,assignment,initials
end

export mykmeansclustering
end


text:
function FixedRate(arg0::Attitude)
    return FixedRate((Attitude,), arg0)
end

function get_attitude(obj::FixedRate, arg0::FieldPVCoordinatesProvider, arg1::FieldAbsoluteDate, arg2::Frame)
    return jcall(obj, "getAttitude", FieldAttitude, (FieldPVCoordinatesProvider, FieldAbsoluteDate, Frame), arg0, arg1, arg2)
end

function get_attitude(obj::FixedRate, arg0::PVCoordinatesProvider, arg1::AbsoluteDate, arg2::Frame)
    return jcall(obj, "getAttitude", Attitude, (PVCoordinatesProvider, AbsoluteDate, Frame), arg0, arg1, arg2)
end

function get_reference_attitude(obj::FixedRate)
    return jcall(obj, "getReferenceAttitude", Attitude, ())
end


text:
# Refined Löwner-John ellipsoid example with Mosek, MosekTools, and JuMP

# Computes the inner and outer Löwner-John ellipsoids around the feasible
# space of the JuMP model.

# Uses Polyhedra and CDDLib packages to express the polyhedral set of the JuMP
# model by inequalities and a convex hull of a finite set of points

using JuMP, Mosek, MosekTools
using LinearAlgebra, Polyhedra, CDDLib
using Plots

function lowner_john_inner_ellipsoid(A,b)
    """
    Löwner-John inner ellipsoid:
        the function takes H-representation (A,b) of a polytope and
        outputs the parameters of the maximum volume inscribed ellipsoid.

        polyhedral set S:
        S = { x ∈ Rn | ai' x ⩽ bi, i= 1,...,m }

        Proposition 3.7.1 [BenTalN01]
            ... Then the largest volume ellipsoid contained in S is
            E = { x = Z∗u + z∗ | uTu ≤1 }, where (Z∗,z∗) are from

            maximize    t                                           (1.0)
                s.t.    t ⩽ (DetZ)^1/n                              (1.1)
                        Z - PSD                                     (1.2)
                        || Z * ai ||_2 ≤ bi −ai' * z, i = 1,...,m   (1.3)
                        Z ∈ S^n, z ∈ Rn, t ∈ R                      (1.4)

        To model (1.1), refer to [Mosek, 6.2.3 Log-determinant]

            t ⩽ (DetZ)^1/n <==> [Z B; B' diag(B)] - PSD             (2.1)
                                B is lower triangular               (2.2)
                                t ⩽ (Π_i={1,...,n} Bii)^1/n         (2.3)

            where (2.3) maximizes the geometric mean of the Bii variables. It
            can be modelled using power cone [Mosek, 4.2.4 Geometric mean] and
            MOI.GeometricMeanCone available with the MathOptInterface.jl

    [BenTalN01]	    A. Ben-Tal and A. Nemirovski. Lectures on Modern Convex
                    Optimization: Analysis, Algorithms, and Engineering Applications.
                    MPS/SIAM Series on Optimization. SIAM, 2001.
    [Mosek]         https://docs.mosek.com/modeling-cookbook/index.html
    """
    # model size
    m,n = size(A)
    # model declaraton
    model = Model(optimizer_with_attributes(Mosek.Optimizer, "LOG" => 0))
    # variables
    @variable(model, t)
    @variable(model, Z[1:n,1:n], PSD)
    @variable(model, z[1:n])
    @variable(model, B[1:n,1:n])
    # maximize the volume
    @objective(model, Max, t)
    # SOC constraint (1.3)
    @constraint(model, SOC[i=1:m], [ b[i] - (A[i,:]' * z) ; (A[i,:]'*Z)'] in SecondOrderCone())
    # PSD constraint (2.1)
    @constraint(model, [Z LowerTriangular(B); LowerTriangular(B)' diagm(diag(B))] in PSDCone())
    # Power cone costraint (2.3)
    @constraint(model, [t; vec([B[i,i] for i in 1:n])] in MOI.GeometricMeanCone(n+1))
    # optimize
    optimize!(model)
    status=termination_status(model)
    CPU_time = MOI.get(model, MOI.SolveTime())
    # return solution
    return JuMP.value.(Z), JuMP.value.(z), CPU_time
end

function lownerjohn_outer(x)
    """
    Löwner-John outer ellipsoid:
        the function takes V-representation (x) of a polytope and
        outputs the parameters of the maximum volume inscribed ellipsoid.

        Polyhedral set S = Conv{x_1,...,x_n}

        Uses the same Mosek features as the inner ellipsoid function,
        refer to Proposition 3.7.2 [BenTalN01] for further details.
    """

    m,n = size(x)
    model = Model(optimizer_with_attributes(Mosek.Optimizer, "LOG" => 0))

    @variable(model,t >= 0)
    @variable(model,P[1:n,1:n],PSD)
    @variable(model,B[1:n,1:n])
    @variable(model,c[1:n])

    # (1, Px-c) in cone
    for i in 1:m
        @constraint(model,[1 ; (x[i,:]' * P)' - c] in MOI.SecondOrderCone(n+1))
    end

    @constraint(model, [P LowerTriangular(B); LowerTriangular(B)' diagm(diag(B))] in PSDCone())
    @constraint(model, [t; vec([B[i,i] for i in 1:n])] in MOI.GeometricMeanCone(n+1))

    @objective(model,Max,t)

    optimize!(model)
    status=termination_status(model)
    CPU_time = MOI.get(model, MOI.SolveTime())

    JuMP.value.(P), JuMP.value.(c), CPU_time
end

# create a JuMP model
function JuMP_model(n)
    model = Model(optimizer_with_attributes(Mosek.Optimizer, "LOG" => 0))
    @variable(model, x[1:n])
    @objective(model, Min, x'diagm(ones(n))'x)
    @constraint(model, con1[i=1:n], 0 <= x[i] <= 1)
    @constraint(model, ones(n)'x >= n/3)
    optimize!(model)
    return model
end

# model dimension
n = 2
# get JuMP model
model = JuMP_model(n)
@info("done exporting JuMP model")
# create a polyhedron from the feasible set of a JuMP model
poly = polyhedron(model, CDDLib.Library(:exact))
@info("done exporting the polyhedron")

# obtain H-representation of the polyhedron
hr = MixedMatHRep(hrep(poly))
A = hr.A
b = hr.b
@info("done obtaining H-representation")
# solve ellipsoid problem
Ci, di, CPU_time = lowner_john_inner_ellipsoid(A,b)
@info("done solving inner ellipsoid problem in $(CPU_time) seconds")

# obtain V-representation of the polyhedron
vr = vrep(poly)
vr = MixedMatVRep(vr)
p = vr.V
@info("done obtaining V-representation")
Po, co, CPU_time = lownerjohn_outer(p)
@info("done solving outter ellipsoid problem in $(CPU_time) seconds")
Poinv = Po^-1

# make a plot
if n == 2
    println("Plotting solution...")
    Plots.gr()
    Plots.plot([p[:,1] ; p[1,1]], [p[:,2] ; p[1,2]])
    Plots.plot!(t -> Ci[1,:]' * [ cos(t),sin(t) ] + di[1],
                t -> Ci[2,:]' * [ cos(t),sin(t) ] + di[2],
                0,2*pi)
    Plots.plot!(t -> Poinv[1,:]' * ( [cos(t),sin(t)] .+ co[1] ),
                t -> Poinv[2,:]' * ( [cos(t),sin(t)] .+ co[2] ),
                0,2*pi)
end


text:
# F16 constants
const MASS = 20500.0 * LB2KG # kgf

const AXX = 9496.0 * SLUG2KG * FT2M^2 # slug·ft² -> kg·m²
const AYY = 55814.0 * SLUG2KG * FT2M^2  # slug·ft² -> kg·m²
const AZZ = 63100.0 * SLUG2KG * FT2M^2 # slug·ft² -> kg·m²
const AXZ = 982.0 * SLUG2KG * FT2M^2 # slug·ft² -> kg·m²

const INERTIA = @SMatrix [
    AXX 0.0 AXZ;
    0.0 AYY 0.0;
    AXZ 0.0 AZZ
]  # Kg·m²

const XCGR = 0.35  # units MAC

const S = 300 * FT2M^2  # m²
const CBAR = 11.32 * FT2M # m
const B = 30 * FT2M  # m

const HX = 160.0 * SLUG2KG * FT2M^2 # slug·ft² -> kg·m²

const DE_MAX = 25.0  # deg
const DA_MAX = 20.0  # deg  #XXX: In Stevens' book says 21.5 deg (Appendix A Section A.4)
const DR_MAX = 30.0  # deg

# Get the GD value used in Stevens to pass tests
# TODO: Note that this is different from gD * M2FT!!
const GD = floor(gD * M2FT, digits=2)

# TODO: doc
struct F16{T}<:Aircraft
    mass::T
    inertia::SMatrix{3, 3, T}
    cg_mac::T
end

F16() = F16(MASS, INERTIA, XCGR)

get_mass(ac::F16) = ac.mass
get_inertia_tensor(ac::F16) = ac.inertia
get_xcg_mac(ac::F16) = ac.cg_mac

get_surface(ac::F16) = S
get_chord(ac::F16) = CBAR
get_wing_span(ac::F16) = B


text:
using LegibleLambdas: @λ
using SimpleWorkflows

j = AtomicJob(`ls`)
k = AtomicJob(() -> sleep(5); desc = "Sleep for 5 seconds")
l = AtomicJob(`sleep 3`; desc = "Sleep for 3 seconds")
m = AtomicJob(@λ(() -> sin(1)))
n = AtomicJob(`pwd` & `sleep 3`)
j ▷ l ▷ k ▷ n ▷ m
j ▷ k
k ▷ m
w = Workflow(l, k, j, m, n)
@test w.jobs == Workflow(k, j, l, n, m).jobs == Workflow(k, l, m, n, j).jobs
run!(w)
@test something(getresult(j)) isa Base.Process
@test something(getresult(k)) === nothing
@test something(getresult(l)) isa Base.Process
@test something(getresult(m)) == 0.8414709848078965
@test something(getresult(n)) isa Base.ProcessChain


text:
Text provided under a Creative Commons Attribution license, CC-BY.  All code is made available under the FSF-approved BSD-3 license.  (c) Lorena A. Barba, Gilbert F. Forsyth 2017. Thanks to NSF for support via CAREER award #1149784.
[@LorenaABarba](https://twitter.com/LorenaABarba)

12 steps to Navier–Stokes
=====
***

We continue our journey to solve the Navier–Stokes equation with Step 4. But don't continue unless you have completed the previous steps! In fact, this next step will be a combination of the two previous ones. The wonders of *code reuse*!

Step 4: Burgers' Equation
----
***

You can read about Burgers' Equation on its [wikipedia page](http://en.wikipedia.org/wiki/Burgers'_equation).

Burgers' equation in one spatial dimension looks like this:

$$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = \nu \frac{\partial ^2u}{\partial x^2}$$

As you can see, it is a combination of non-linear convection and diffusion. It is surprising how much you learn from this neat little equation!

We can discretize it using the methods we've already detailed in Steps [1](./01_Step_1.ipynb) to [3](./04_Step_3.ipynb).  Using forward difference for time, backward difference for space and our 2nd-order method for the second derivatives yields:

$$\frac{u_i^{n+1}-u_i^n}{\Delta t} + u_i^n \frac{u_i^n - u_{i-1}^n}{\Delta x} = \nu \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\Delta x^2}$$

As before, once we have an initial condition, the only unknown is $u_i^{n+1}$. We will step in time as follows:

$$u_i^{n+1} = u_i^n - u_i^n \frac{\Delta t}{\Delta x} (u_i^n - u_{i-1}^n) + \nu \frac{\Delta t}{\Delta x^2}(u_{i+1}^n - 2u_i^n + u_{i-1}^n)$$

### Initial and Boundary Conditions

To examine some interesting properties of Burgers' equation, it is helpful to use different initial and boundary conditions than we've been using for previous steps.

Our initial condition for this problem is going to be:

\begin{eqnarray}
u &=& -\frac{2 \nu}{\phi} \frac{\partial \phi}{\partial x} + 4 \\\
\phi &=& \exp \bigg(\frac{-x^2}{4 \nu} \bigg) + \exp \bigg(\frac{-(x-2 \pi)^2}{4 \nu} \bigg)
\end{eqnarray}

This has an analytical solution, given by:

\begin{eqnarray}
u &=& -\frac{2 \nu}{\phi} \frac{\partial \phi}{\partial x} + 4 \\\
\phi &=& \exp \bigg(\frac{-(x-4t)^2}{4 \nu (t+1)} \bigg) + \exp \bigg(\frac{-(x-4t -2 \pi)^2}{4 \nu(t+1)} \bigg)
\end{eqnarray}

Our boundary condition will be:

$$u(0) = u(2\pi)$$

This is called a *periodic* boundary condition. Pay attention! This will cause you a bit of headache if you don't tread carefully.

### Saving Time with SymPy

The initial condition we're using for Burgers' Equation can be a bit of a pain to evaluate by hand.  The derivative $\frac{\partial \phi}{\partial x}$ isn't too terribly difficult, but it would be easy to drop a sign or forget a factor of $x$ somewhere, so we're going to use SymPy to help us out.

[SymPy](http://sympy.org/en/) is the symbolic math library for Python.  It has a lot of the same symbolic math functionality as Mathematica with the added benefit that we can easily translate its results back into our Python calculations (it is also free and open source).

Start by loading the SymPy library, together with our favorite library, NumPy.

```python
import numpy
import sympy
```

We're also going to tell SymPy that we want all of its output to be rendered using $\LaTeX$. This will make our Notebook beautiful!

```python
from sympy import init_printing
init_printing(use_latex=True)
```

Start by setting up symbolic variables for the three variables in our initial condition and then type out the full equation for $\phi$.  We should get a nicely rendered version of our $\phi$ equation.

```python
x, nu, t = sympy.symbols('x nu t')
phi = (sympy.exp(-(x - 4 * t)**2 / (4 * nu * (t + 1))) +
       sympy.exp(-(x - 4 * t - 2 * sympy.pi)**2 / (4 * nu * (t + 1))))
phi
```

$$e^{- \frac{\left(- 4 t + x - 2 \pi\right)^{2}}{4 \nu \left(t + 1\right)}} + e^{- \frac{\left(- 4 t + x\right)^{2}}{4 \nu \left(t + 1\right)}}$$

It's maybe a little small, but that looks right.  Now to evaluate our partial derivative $\frac{\partial \phi}{\partial x}$ is a trivial task.

```python
phiprime = phi.diff(x)
phiprime
```

$$- \frac{e^{- \frac{\left(- 4 t + x\right)^{2}}{4 \nu \left(t + 1\right)}}}{4 \nu \left(t + 1\right)} \left(- 8 t + 2 x\right) - \frac{1}{4 \nu \left(t + 1\right)} \left(- 8 t + 2 x - 4 \pi\right) e^{- \frac{\left(- 4 t + x - 2 \pi\right)^{2}}{4 \nu \left(t + 1\right)}}$$

If you want to see the unrendered version, just use the Python print command.

```python
print(phiprime)
```

    -(-8*t + 2*x)*exp(-(-4*t + x)**2/(4*nu*(t + 1)))/(4*nu*(t + 1)) - (-8*t + 2*x - 4*pi)*exp(-(-4*t + x - 2*pi)**2/(4*nu*(t + 1)))/(4*nu*(t + 1))

### Now what?

Now that we have the Pythonic version of our derivative, we can finish writing out the full initial condition equation and then translate it into a usable Python expression.  For this, we'll use the *lambdify* function, which takes a SymPy symbolic equation and turns it into a callable function.

```python
from sympy.utilities.lambdify import lambdify

u = -2 * nu * (phiprime / phi) + 4
print(u)
```

    -2*nu*(-(-8*t + 2*x)*exp(-(-4*t + x)**2/(4*nu*(t + 1)))/(4*nu*(t + 1)) - (-8*t + 2*x - 4*pi)*exp(-(-4*t + x - 2*pi)**2/(4*nu*(t + 1)))/(4*nu*(t + 1)))/(exp(-(-4*t + x - 2*pi)**2/(4*nu*(t + 1))) + exp(-(-4*t + x)**2/(4*nu*(t + 1)))) + 4

### Lambdify

To lambdify this expression into a useable function, we tell lambdify which variables to request and the function we want to plug them in to.

```python
ufunc = lambdify((t, x, nu), u)
print(ufunc(1, 4, 3))
```

    3.49170664206

### Back to Burgers' Equation

Now that we have the initial conditions set up, we can proceed and finish setting up the problem.  We can generate the plot of the initial condition using our lambdify-ed function.

```python
from matplotlib import pyplot
%matplotlib inline

###variable declarations
nx = 101
nt = 100
dx = 2 * numpy.pi / (nx - 1)
nu = .07
dt = dx * nu

x = numpy.linspace(0, 2 * numpy.pi, nx)
un = numpy.empty(nx)
t = 0

u = numpy.asarray([ufunc(t, x0, nu) for x0 in x])
u
```

    array([ 4.        ,  4.06283185,  4.12566371,  4.18849556,  4.25132741,
            4.31415927,  4.37699112,  4.43982297,  4.50265482,  4.56548668,
            4.62831853,  4.69115038,  4.75398224,  4.81681409,  4.87964594,
            4.9424778 ,  5.00530965,  5.0681415 ,  5.13097336,  5.19380521,
            5.25663706,  5.31946891,  5.38230077,  5.44513262,  5.50796447,
            5.57079633,  5.63362818,  5.69646003,  5.75929189,  5.82212374,
            5.88495559,  5.94778745,  6.0106193 ,  6.07345115,  6.136283  ,
            6.19911486,  6.26194671,  6.32477856,  6.38761042,  6.45044227,
            6.51327412,  6.57610598,  6.63893783,  6.70176967,  6.76460125,
            6.82742866,  6.89018589,  6.95176632,  6.99367964,  6.72527549,
            4.        ,  1.27472451,  1.00632036,  1.04823368,  1.10981411,
            1.17257134,  1.23539875,  1.29823033,  1.36106217,  1.42389402,
            1.48672588,  1.54955773,  1.61238958,  1.67522144,  1.73805329,
            1.80088514,  1.863717  ,  1.92654885,  1.9893807 ,  2.05221255,
            2.11504441,  2.17787626,  2.24070811,  2.30353997,  2.36637182,
            2.42920367,  2.49203553,  2.55486738,  2.61769923,  2.68053109,
            2.74336294,  2.80619479,  2.86902664,  2.9318585 ,  2.99469035,
            3.0575222 ,  3.12035406,  3.18318591,  3.24601776,  3.30884962,
            3.37168147,  3.43451332,  3.49734518,  3.56017703,  3.62300888,
            3.68584073,  3.74867259,  3.81150444,  3.87433629,  3.93716815,  4.        ])

```python
pyplot.figure(figsize=(11, 7), dpi=100)
pyplot.plot(x, u, marker='o', lw=2)
pyplot.xlim([0, 2 * numpy.pi])
pyplot.ylim([0, 10]);
```

![png](output_23_0.png)

This is definitely not the hat function we've been dealing with until now. We call it a "saw-tooth function".  Let's proceed forward and see what happens.

### Periodic Boundary Conditions

One of the big differences between Step 4 and the previous lessons is the use of *periodic* boundary conditions.  If you experiment with Steps 1 and 2 and make the simulation run longer (by increasing `nt`) you will notice that the wave will keep moving to the right until it no longer even shows up in the plot.

With periodic boundary conditions, when a point gets to the right-hand side of the frame, it *wraps around* back to the front of the frame.

Recall the discretization that we worked out at the beginning of this notebook:

$$u_i^{n+1} = u_i^n - u_i^n \frac{\Delta t}{\Delta x} (u_i^n - u_{i-1}^n) + \nu \frac{\Delta t}{\Delta x^2}(u_{i+1}^n - 2u_i^n + u_{i-1}^n)$$

What does $u_{i+1}^n$ *mean* when $i$ is already at the end of the frame?

Think about this for a minute before proceeding.

```python
for n in range(nt):
    un = u.copy()
    for i in range(1, nx-1):
        u[i] = un[i] - un[i] * dt / dx *(un[i] - un[i-1]) + nu * dt / dx**2 *\
                (un[i+1] - 2 * un[i] + un[i-1])
    u[0] = un[0] - un[0] * dt / dx * (un[0] - un[-2]) + nu * dt / dx**2 *\
                (un[1] - 2 * un[0] + un[-2])
    u[-1] = u[0]

u_analytical = numpy.asarray([ufunc(nt * dt, xi, nu) for xi in x])
```

```python
pyplot.figure(figsize=(11, 7), dpi=100)
pyplot.plot(x,u, marker='o', lw=2, label='Computational')
pyplot.plot(x, u_analytical, label='Analytical')
pyplot.xlim([0, 2 * numpy.pi])
pyplot.ylim([0, 10])
pyplot.legend();
```

![png](output_27_0.png)

***

What next?
----

The subsequent steps, from 5 to 12, will be in two dimensions. But it is easy to extend the 1D finite-difference formulas to the partial derivatives in 2D or 3D. Just apply the definition — a partial derivative with respect to $x$ is the variation in the $x$ direction *while keeping $y$ constant*.

Before moving on to [Step 5](./07_Step_5.ipynb), make sure you have completed your own code for steps 1 through 4 and you have experimented with the parameters and thought about what is happening. Also, we recommend that you take a slight break to learn about [array operations with NumPy](./06_Array_Operations_with_NumPy.ipynb).

```python
from IPython.core.display import HTML
def css_styling():
    styles = open("../styles/custom.css", "r").read()
    return HTML(styles)
css_styling()
```

<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>
<style>
    @font-face {
        font-family: "Computer Modern";
        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
    }
    div.cell{
        width:800px;
        margin-left:16% !important;
        margin-right:auto;
    }
    h1 {
        font-family: 'Alegreya Sans', sans-serif;
    }
    h2 {
        font-family: 'Fenix', serif;
    }
    h3{
		font-family: 'Fenix', serif;
        margin-top:12px;
        margin-bottom: 3px;
       }
	h4{
		font-family: 'Fenix', serif;
       }
    h5 {
        font-family: 'Alegreya Sans', sans-serif;
    }
    div.text_cell_render{
        font-family: 'Alegreya Sans',Computer Modern, "Helvetica Neue", Arial, Helvetica, Geneva, sans-serif;
        line-height: 135%;
        font-size: 120%;
        width:600px;
        margin-left:auto;
        margin-right:auto;
    }
    .CodeMirror{
            font-family: "Source Code Pro";
			font-size: 90%;
    }
/*    .prompt{
        display: None;
    }*/
    .text_cell_render h1 {
        font-weight: 200;
        font-size: 50pt;
		line-height: 100%;
        color:#CD2305;
        margin-bottom: 0.5em;
        margin-top: 0.5em;
        display: block;
    }
    .text_cell_render h5 {
        font-weight: 300;
        font-size: 16pt;
        color: #CD2305;
        font-style: italic;
        margin-bottom: .5em;
        margin-top: 0.5em;
        display: block;
    }

    .warning{
        color: rgb( 240, 20, 20 )
        }
</style>
<script>
    MathJax.Hub.Config({
                        TeX: {
                           extensions: ["AMSmath.js"]
                           },
                tex2jax: {
                    inlineMath: [ ['$','$'], ["\\(","\\)"] ],
                    displayMath: [ ['$$','$$'], ["\\[","\\]"] ]
                },
                displayAlign: 'center', // Change this to 'center' to center equations.
                "HTML-CSS": {
                    styles: {'.MathJax_Display': {"margin": 4}}
                }
        });
</script>

max_stars_repo_name: mkitti/CFDJulia

max_issues_repo_name: mkitti/CFDJulia

max_forks_repo_name: mkitti/CFDJulia


text:
module FunctionWrapper

using JavaCall

const Abs = @jimport org.hipparchus.analysis.function.Abs
const Acos = @jimport org.hipparchus.analysis.function.Acos
const Acosh = @jimport org.hipparchus.analysis.function.Acosh
const Add = @jimport org.hipparchus.analysis.function.Add
const Asin = @jimport org.hipparchus.analysis.function.Asin
const Asinh = @jimport org.hipparchus.analysis.function.Asinh
const Atan = @jimport org.hipparchus.analysis.function.Atan
const Atan2 = @jimport org.hipparchus.analysis.function.Atan2
const Atanh = @jimport org.hipparchus.analysis.function.Atanh
const Cbrt = @jimport org.hipparchus.analysis.function.Cbrt
const Ceil = @jimport org.hipparchus.analysis.function.Ceil
const Class = @jimport java.lang.Class
const Constant = @jimport org.hipparchus.analysis.function.Constant
const Cos = @jimport org.hipparchus.analysis.function.Cos
const Cosh = @jimport org.hipparchus.analysis.function.Cosh
const Derivative = @jimport org.hipparchus.analysis.differentiation.Derivative
const DerivativeStructure = @jimport org.hipparchus.analysis.differentiation.DerivativeStructure
const Divide = @jimport org.hipparchus.analysis.function.Divide
const Exp = @jimport org.hipparchus.analysis.function.Exp
const Expm1 = @jimport org.hipparchus.analysis.function.Expm1
const Floor = @jimport org.hipparchus.analysis.function.Floor
const Gaussian = @jimport org.hipparchus.analysis.function.Gaussian
const HarmonicOscillator = @jimport org.hipparchus.analysis.function.HarmonicOscillator
const Identity = @jimport org.hipparchus.analysis.function.Identity
const Inverse = @jimport org.hipparchus.analysis.function.Inverse
const JString = @jimport java.lang.String
const Log = @jimport org.hipparchus.analysis.function.Log
const Log10 = @jimport org.hipparchus.analysis.function.Log10
const Log1p = @jimport org.hipparchus.analysis.function.Log1p
const Logistic = @jimport org.hipparchus.analysis.function.Logistic
const Logit = @jimport org.hipparchus.analysis.function.Logit
const Max = @jimport org.hipparchus.analysis.function.Max
const Min = @jimport org.hipparchus.analysis.function.Min
const Minus = @jimport org.hipparchus.analysis.function.Minus
const Multiply = @jimport org.hipparchus.analysis.function.Multiply
const Object = @jimport java.lang.Object
const Pow = @jimport org.hipparchus.analysis.function.Pow
const Power = @jimport org.hipparchus.analysis.function.Power
const Rint = @jimport org.hipparchus.analysis.function.Rint
const Sigmoid = @jimport org.hipparchus.analysis.function.Sigmoid
const Signum = @jimport org.hipparchus.analysis.function.Signum
const Sin = @jimport org.hipparchus.analysis.function.Sin
const Sinc = @jimport org.hipparchus.analysis.function.Sinc
const Sinh = @jimport org.hipparchus.analysis.function.Sinh
const Sqrt = @jimport org.hipparchus.analysis.function.Sqrt
const StepFunction = @jimport org.hipparchus.analysis.function.StepFunction
const Subtract = @jimport org.hipparchus.analysis.function.Subtract
const Tan = @jimport org.hipparchus.analysis.function.Tan
const Tanh = @jimport org.hipparchus.analysis.function.Tanh
const Ulp = @jimport org.hipparchus.analysis.function.Ulp

export Abs
export Acos
export Acosh
export Add
export Asin
export Asinh
export Atan
export Atan2
export Atanh
export Cbrt
export Ceil
export Constant
export Cos
export Cosh
export Divide
export Exp
export Expm1
export Floor
export Gaussian
export HarmonicOscillator
export Identity
export Inverse
export Log
export Log10
export Log1p
export Logistic
export Logit
export Max
export Min
export Minus
export Multiply
export Pow
export Power
export Rint
export Sigmoid
export Signum
export Sin
export Sinc
export Sinh
export Sqrt
export StepFunction
export Subtract
export Tan
export Tanh
export Ulp
export equals
export get_class
export hash_code
export notify
export notify_all
export to_string
export value
export wait

include("abs.jl")
include("acos.jl")
include("acosh.jl")
include("add.jl")
include("asin.jl")
include("asinh.jl")
include("atan.jl")
include("atan2.jl")
include("atanh.jl")
include("cbrt.jl")
include("ceil.jl")
include("constant.jl")
include("cos.jl")
include("cosh.jl")
include("divide.jl")
include("exp.jl")
include("expm1.jl")
include("floor.jl")
include("gaussian.jl")
include("harmonic_oscillator.jl")
include("identity.jl")
include("inverse.jl")
include("log.jl")
include("log10.jl")
include("log1p.jl")
include("logistic.jl")
include("logit.jl")
include("max.jl")
include("min.jl")
include("minus.jl")
include("multiply.jl")
include("pow.jl")
include("power.jl")
include("rint.jl")
include("sigmoid.jl")
include("signum.jl")
include("sin.jl")
include("sinc.jl")
include("sinh.jl")
include("sqrt.jl")
include("step_function.jl")
include("subtract.jl")
include("tan.jl")
include("tanh.jl")
include("ulp.jl")

end


text:
import HoTRG_lemon.SpinModule:
	initialize!

"""
	Quantum2dFractalLattice

# arguments

* `spinmodel::QuantumSpinModel`
* `legextension::Int`
  default value is one. set zero for the normal square lattice.
"""
type Quantum2dFractalLattice{T} <: Quantum2dLattice{T}
	spinmodel::QuantumSpinModel{T}
	legextension::Int
	tensorT::Array{T,6} # ; x,xp,y,yp,z,zp
	tensorTtilde::Array{T,6}
	tensorP::Array{Array{T,5},1}# leg tensor; x,xp,s,z,zp
	tensorQ::Array{T,4} # corner tensor; x,y,z,zp
	function Quantum2dFractalLattice{T}(spinmodel::QuantumSpinModel{T}, legextension::Int)
		this = new{T}()
 		this.spinmodel = spinmodel
		this.legextension = legextension
 		initialize!(this)
 		return this
 	end
end

Quantum2dFractalLattice{T}(spinmodel::QuantumSpinModel{T}, legextension::Int) =
Quantum2dFractalLattice{T}(spinmodel, legextension)
Quantum2dFractalLattice{T}(spinmodel::QuantumSpinModel{T}) =
Quantum2dFractalLattice{T}(spinmodel, 1)

#---
# initialize functions for constructor
"""
	initialize(lattice::Quantum2dFractalLattice)
initialize and set tensorT, tensorTtilde, tensorP and tensorQ.
"""
function initialize!{T}(lattice::Quantum2dFractalLattice{T})
	lattice.tensorT = initializeTensorT2dQ(lattice.spinmodel)
	lattice.tensorTtilde = initializeTensorT2dQ(lattice.spinmodel; tilde = true)
	tensorP = initializeTensorP_q2d(lattice)
	lattice.tensorP = fill(tensorP, getSpaceDimension(lattice))
	tensorQ = initializeTensorQ_q2d(lattice)
	setTensorQ!(lattice, tensorQ)
end

function initializeTensorP_q2d{T}(lattice::Quantum2dFractalLattice{T})
	factorW = getFactorW(lattice)
	factorWp = getFactorWp(lattice)
	leng = size(factorW,1)==size(factorW,2) ? size(factorW,1) : throw(DomainError())
	tensorPx = zeros(T,leng,leng,leng,leng,leng) #; y yp s z zp
	for a = 1:leng, l = 1:leng, r = 1:leng, f=1:leng, u=1:leng, v=1:leng
		tensorPx[l,r,f,u,v] +=
		factorW[a,l] * factorW[a,r] * factorW[a,f] *
		factorWp[a,u] * factorWp[a,v]
	end
	return tensorPx
end

function initializeTensorQ_q2d{T}(lattice::Quantum2dFractalLattice{T})
	factorW = getFactorW(lattice)
	factorWp = getFactorWp(lattice)
	leng = size(factorW,1)==size(factorW,2) ? size(factorW,1) : throw(DomainError())
	tensorQ = zeros(T,leng,leng,leng,leng) #; x y z zp
	for a = 1:leng, x = 1:leng, y = 1:leng, z = 1:leng, zp = 1:leng
		tensorQ[x,y,z,zp] +=
		factorW[a,x] * factorW[a,y] *
		factorWp[a,z] * factorWp[a,zp]
	end
	return tensorQ
end

#---
# functions for Quantum2dFractalLattice
function getCoarserate(lattice::Quantum2dFractalLattice)
	##### maybe not right... or even I don't need this
	L = getLegextension(lattice)
	return 8*(L+1)
end

#---
# functions for test
function testQuantum2dFractalLattice()
	spinmodel = testQuantumIsingModel()
	latticecode = "quantum_2d_fractal_1"
	lattice = buildLattice(latticecode, spinmodel)
	return lattice
end


text:
include("../constants.jl")
include("../CESMReader.jl")
include("../LinearRegression.jl")
include("../PCA.jl")

using .PCA
using NCDatasets
using Formatting
using ArgParse
using Statistics
using JSON
using .CESMReader

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--data-file-timestamp-form"
            help = "Data filename timestamp form. Either `YEAR` or `YEAR_MONTH`."
            arg_type = String
            required = true

        "--dims"
            help = "The form of dimensions. Now can be `XYZT`, `XYT`, `YZT`"
            arg_type = String
            required = true

        "--data-file-prefix"
            help = "Input data filename prefix including folder and path until the timestamp. File extension `nc` is assumed."
            arg_type = String
            required = true

        "--output-file"
            help = "Output data filename"
            arg_type = String
            required = true

        "--domain-file"
            help = "Domain file."
            arg_type = String
            required = true

        "--SLP-varname"
            help = "Variable name of sea-level pressure"
            arg_type = String
            default = "PSL"

        "--beg-year"
            help = "Year of begin."
            arg_type = Int64
            required = true

        "--end-year"
            help = "Year of end."
            arg_type = Int64
            required = true

        "--overwrite"
            action = :store_true

        "--sparsity"
            help = "Because solving PCA matrix is time consuming. This parameter let user to use coarse data to derive PCA. Sparseity `n` means to skip every `n-1` other grid point in lon / lat. So sparsity 1 (default) means do not skip any data point. Sparsity 1 means means data density will drop to 1/4 of the original (1/2 * 1/2)."
            arg_type = Int64
            default = 1

    end

    return parse_args(ARGS, s)
end

parsed = parse_commandline()
print(json(parsed, 4))

Dataset(parsed["domain-file"], "r") do ds

    global lon = replace(ds["xc"][:], missing=>NaN)
    global lat = replace(ds["yc"][:], missing=>NaN)
    global mask = replace(ds["mask"][:], missing=>NaN)
    global area = replace(ds["area"][:], missing=>NaN)

    global Nx, Ny = size(lat)

    # sparsity
    global sparsity_mask = copy(mask)
    sparsity_mask .= 0.0

    skip = parsed["sparsity"]

    for i=1:Nx, j=1:Ny
        if i % skip == 0 && j % skip == 0
            sparsity_mask[i, j] = 1.0
        end
    end

    # https://climatedataguide.ucar.edu/climate-data/hurrell-north-atlantic-oscillation-nao-index-pc-based
    # "... leading Empirical Orthogonal Function (EOF) of SLP
    #  anomalies over the Atlantic sector, 20°-80°N, 90°W-40°E."

    global EOF_valid_pts  = (sparsity_mask .== 1) .& ( 20 .<= lat .<= 80 ) .& (( 270 .<= lon .<= 360) .| (0 .<= lon .<= 40.0))

end

let
    if parsed["data-file-timestamp-form"] == "YEAR"
        filename_format = format("{:s}{{:04d}}.nc", joinpath(parsed["data-file-prefix"]))
        form = :YEAR
    elseif parsed["data-file-timestamp-form"] == "YEAR_MONTH"
        filename_format = format("{:s}{{:04d}}-{{:02d}}.nc", joinpath(parsed["data-file-prefix"]))
        form = :YEAR_MONTH
    end

    fh = FileHandler(filename_format=filename_format, form=form)

    beg_t = (parsed["beg-year"] - 1) * 12 + 1
    end_t = (parsed["end-year"] - 1) * 12 + 12

    global SLP = getData(fh, parsed["SLP-varname"], (parsed["beg-year"], parsed["end-year"]), (:, :))

    global Nt = end_t - beg_t + 1
end

modes = 4
PCAs, PCAs_ts = PCA.findPCAs(SLP, EOF_valid_pts; modes=modes)

Dataset(parsed["output-file"], "c") do ds

    defDim(ds, "Nx", Nx)
    defDim(ds, "Ny", Ny)
    defDim(ds, "time", Inf)
    defDim(ds, "modes", modes)

    for (varname, vardata, vardim, attrib) in [
        ("PCAs",  reshape(PCAs, Nx, Ny, modes), ("Nx", "Ny", "modes",), Dict()),
        ("PCAs_ts",  PCAs_ts, ("modes", "time",), Dict()),
    ]

        if ! haskey(ds, varname)
            var = defVar(ds, varname, Float64, vardim)
            var.attrib["_FillValue"] = 1e20
        end

        var = ds[varname]

        for (k, v) in attrib
            var.attrib[k] = v
        end

        rng = []
        for i in 1:length(vardim)-1
            push!(rng, Colon())
        end
        push!(rng, 1:size(vardata)[end])
        var[rng...] = vardata

    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
using Oceananigans.Models.HydrostaticFreeSurfaceModels: ColumnEnsembleSize

const OneDimensionalEnsembleGrid = RegularRectilinearGrid{<:Any, Flat, Flat, Bounded}

"""
    OneDimensionalEnsembleGrid(datapath; ensemble_size = (1,1), Nz = nothing)

Construct a 3D grid of type `Oceananigans.RegularRectilinearGrid`
with independent columns, each with `Nz` grid points and
arranged horizontally according to `(Nx, Ny)`.
The grid is identical in extent to that of the LES simulation
stored in `datapath`. `Nz` defaults to the simulation resolution.
"""
function OneDimensionalEnsembleGrid(datapath; size = (1,1,nothing))

    _, _, Nz_, _, _, Lz_ = get_grid_params(datapath)

    Nz = isnothing(size[3]) ? Nz_ : size[3]

    @assert Nz <=  Nz_ "Desired grid resolution exceeds the simulation resolution!"
    sz = ColumnEnsembleSize(Nz=Nz, ensemble=size[1:2])
    halo = ColumnEnsembleSize(Nz=Nz)
    grid = RegularRectilinearGrid(size=sz, halo=halo, z=(-Lz_, 0), topology=(Flat, Flat, Bounded))

    return grid
end


text:
using Gtk.ShortNames
using GtkUtilities, GtkUtilities.Graphics

c = @Canvas()
win = @Window(c, "PanZoomCanvas")
showall(win)
panzoom(c, (0,1), (0,1))

draw(c) do widget
    xview, yview = guidata[widget, :xview], guidata[widget, :yview]
    ctx = getgc(widget)
    set_coords(ctx, xview, yview)
    # Paint red rectangle
    rectangle(ctx, 0, 0, 0.5, 0.5)
    set_source_rgb(ctx, 1, 0, 0)
    fill(ctx)
    # Paint blue rectangle
    rectangle(ctx, 0.5, 0, 0.5, 0.5)
    set_source_rgb(ctx, 0, 0, 1)
    fill(ctx)
    # Paint white rectangle
    rectangle(ctx, 0, 0.5, 1, 0.75)
    set_source_rgb(ctx, 1, 1, 1)
    fill(ctx)
    # Paint green rectangle
    rectangle(ctx, 0, 0.75, 1, 0.25)
    set_source_rgb(ctx, 0, 1, 0)
    fill(ctx)
end

panzoom_key(c)
panzoom_mouse(c)

nothing


text:
using .NestedSamplers

samplertype(::Type{<:Nested}) = IsNested()

function AbstractMCMC.sample(post::TransformedPosterior, sampler::Nested, args...; kwargs...)
    ℓ(x) = logdensity(post, x)
    model = NestedModel(ℓ, identity)

    samples, stats = sample(model, sampler, args...; chain_type=Array, kwargs...)
    weights = samples[:, end]
    chain = transform.(Ref(post), eachrow(samples[:,1:end-1]))
    return TupleVector(chain), merge((;weights,), stats)
end