text:
# 6D Vlasov - Poisson Single Species

include("pif3d.jl");

# importall "./pif3d.jl
using Sobol
using ProgressMeter

Np=Int(1e4) #Number of particles
Ns=1 #Number of species
dt=0.05
tmax=5
gamma=1.0 # Gravity <0
q=-1.0 # Electrons
m=1.0

Nx1=1
Nx2=0
Nx3=0

# Default parameters
v01=0;v02=0;eps=0;delta=0;epsi=0; betai=0; betar=0;
sigma1=1;sigma2=1;sigma3=1;
c=1.0; k=0.5; eps=0; B0=zeros(3);

testcase="weibel"

#
if testcase=="weibel"
  eps=1e-3; # Amplitude of perturbation, 0.05 for linear, 0.5 for nonlinear
  betar=sign(q)*1e-3; betai=0;
  k=1.25;    # Wave vector
  sigma1=0.02/sqrt(2);
  sigma2=sqrt(12)*sigma1;
  sigma3=sigma2;
  v01=0;
  v02=0;
  delta=0;
  tmax=150;
  dt=0.2;
elseif testcase=="weibels"
  #Streaming Weibel instability
  sigma1=0.1/sqrt(2);
  sigma2=sigma1;
  k=0.2;
  betai=sign(q)*1e-3; betar=0;
  v01=0.5;
  v02=-0.1;
  delta=1/6.0;
  eps=0;
  tmax=150;
elseif testcase=="landau"
  c=1.0
  k=0.5; eps=0.5; B0=zeros(3);
  tmax=20

end

L1=2*pi/k
L2=2*pi/k
L3=2*pi/k

# Particles
x1n=zeros(Float64,Np)
x2n=zeros(Float64,Np)
x3n=zeros(Float64,Np)

v1n=zeros(Float64,Np)
v2n=zeros(Float64,Np)
v3n=zeros(Float64,Np)
#qn=q*ones(Float64,Np) # Charge
#mn=ones(Float64,Np) # Mass

fn=zeros(Float64,Np)
gn=zeros(Float64,Np)

# Initial condition for one species
f0(x1,x2,x3,v1,v2,v3)=((1+ eps*sin.(k*x1)).*
    exp.(- (v1.^2./sigma1^2 + v3.^2/sigma3^2)/2 ).*
    (delta*exp.(-(v2-v01).^2/2/sigma2.^2) +
    (1-delta)*exp.( - (v2-v02).^2/2/sigma2^2))*
    (2*pi)^(-1.5))/sigma3/sigma1/sigma2

g0(x1,x2,x3,v1,v2,v3)=exp.(-v1.^2./sigma1^2/2 -v3.^2./sigma3^2/2).*
    (delta*exp.(-(v2-v01).^2/2/sigma2.^2) +
    (1-delta)*exp.( - (v2-v02).^2/2/sigma2^2))*
    (2*pi)^(-1.5)/sigma3/sigma1/sigma2/L1/L2/L3;
# Fill with uniform numbers
using Sobol
sob = SobolSeq(6)
skip(sob, 4) # Skip some entries
for i=1:Np
  x1n[i],x2n[i],x3n[i],v1n[i],v2n[i],v3n[i]=Sobol.next(sob)
end

using SpecialFunctions
# Sampling
norminv=p->sqrt(2.).*erfinv.(2.*p - 1.0)
x1n*=L1;x2n*=L2;x3n*=L3;
v1n.=norminv(v1n)*sigma1
v2n.=norminv(v2n)*sigma2
v3n.=norminv(v3n)*sigma3
#Shifts
v2n[1:floor(Int,delta*Np)].+=v01
v2n[floor(Int,delta*Np)+1:end].+=v02
qn=ones(Float64,Np)*q
mn=ones(Float64,Np)*m

# Likelihoods
fn=f0(x1n,x2n,x3n,v1n,v2n,v3n)
gn=g0(x1n,x2n,x3n,v1n,v2n,v3n)

# Coefficients for Runge Kutta
comp_gamma=[1/(2-2^(1/3)), -2^(1/3)/(2-2^(1/3)),  1/(2-2^(1/3))]
comp_gamma=1.0

#Particle mesh coupling
pm=pif3d{Float64}(Nx1,Nx2,Nx3,L1,L2,L3)

print( pm.N1*pm.N2*pm.N3 )
# Fields
rho=Array{Complex{Float64}}(pm.N1,pm.N2,pm.N3)
Phi=Array{Complex{Float64}}(pm.N1,pm.N2,pm.N3)
E1=similar(Phi);E2=similar(Phi);E3=similar(Phi)
B1=similar(Phi);B2=similar(Phi);B3=similar(Phi)

E1n=zeros(Complex{Float64},Np)
E2n=zeros(Complex{Float64},Np)
E3n=zeros(Complex{Float64},Np)
wn=fn./gn

Nt=Int(ceil(tmax/dt))
Epot=zeros(Float64,Nt,3)
Ekin=zeros(Float64,Nt,3)
Momentum=zeros(Float64,Nt,3)
kineticenergy=zeros(Float64,Nt)
ttime=collect((0:Nt-1)*dt)

#Initial poisson solve
rho=accum_osde(pm,x1n,x2n,x3n, wn.*qn)./Np
solve_poisson!(pm,rho,Phi)
E1,E2,E3=gradient(pm,-Phi)
#Initial electric field
B1[:]=0;B2[:]=0;B3[:]=0
# Initial magnetic field
B3[pm.K1.==k,pm.K2.==0.,pm.K3.==0.]=betar + im*betai
B3[pm.K1.==-k,pm.K2.==0.,pm.K3.==0.]=betar - im*betai
B1[pm.N1n,pm.N2n,pm.N3n]=B0[1]
B2[pm.N1n,pm.N2n,pm.N3n]=B0[2]
B3[pm.N1n,pm.N2n,pm.N3n]=B0[3]

# rho[1,pm.Nx2+1,pm.Nx3+1]
print(gauss_error(pm,rho,E1,E2,E3),"\n" )

# Allocate temporary arrays
J1=zeros(Complex{Float64},pm.N1,pm.N2,pm.N3 )
J2=zeros(Complex{Float64},pm.N1,pm.N2,pm.N3 )
J3=zeros(Complex{Float64},pm.N1,pm.N2,pm.N3 )
dt2=dt

DTYPE=Float64
#Diagnostics
Nt=Int(ceil(tmax/dt))
Epot=zeros(DTYPE,Nt,3)
Bpot=zeros(DTYPE,Nt,3)
Ekin=zeros(DTYPE,Nt,3)
Momentum=zeros(DTYPE,Nt,3)
kineticenergy=zeros(DTYPE,Nt)
ttime=collect((0:Nt-1)*dt)
tic()
using ProgressMeter
# pg = Progress(10, 1,"Vlasov")
# for tdx=1:10
# sleep(0.1)
# ProgressMeter.next!(pg)
# end
# stop
#
# integrate_Hp_midpoint(pm,x1n,x2n,x3n,v1n,v2n,v3n,qn,mn,wn,
#                 J1,J2,J3,B1,B2,B3,dt2)

# @time integrate_Hp_midpoint(pm,x1n,x2n,x3n,v1n,v2n,v3n,qn,mn,wn,
#                  J1,J2,J3,B1,B2,B3,dt2)

# @time integrate_Hp_split_sym(pm,x1n,x2n,x3n,v1n,v2n,v3n,qn,mn,wn,
#                                    J1,J2,J3,B1,B2,B3,dt2/2,dt2/2)
tic()
@showprogress 1 "Vlasov-Maxwell 6D" for tdx=1:Nt
  Bpot[tdx,1]=L2norm(pm,B1)/2.
  Bpot[tdx,2]=L2norm(pm,B2)/2.
  Bpot[tdx,3]=L2norm(pm,B3)/2.
  Epot[tdx,1]=L2norm(pm,E1)/2.
  Epot[tdx,2]=L2norm(pm,E2)/2.
  Epot[tdx,3]=L2norm(pm,E3)/2.
  Ekin[tdx,1]=sum(v1n.^2.*wn)/Np/2.
  Ekin[tdx,2]=sum(v2n.^2.*wn)/Np/2.
  Ekin[tdx,3]=sum(v3n.^2.*wn)/Np/2.
  Momentum[tdx,1]=sum(v1n.*wn)/Np + (dot(pm,E2,B3)-dot(pm,E3,B2))
  Momentum[tdx,2]=sum(v2n.*wn)/Np + (dot(pm,E3,B1)-dot(pm,E1,B3))
  Momentum[tdx,3]=sum(v3n.*wn)/Np + (dot(pm,E1,B2)-dot(pm,E2,B1))

 # integrate_H_E(pm,E1,E2,E3,B1,B2,B3,x1n,x2n,x3n,v1n,v2n,v3n,qn,mn,dt)
 # integrate_H_B(pm, E1,E2,E3,B1,B2,B3,c,dt)
 #  J1[:]=0.0;J2[:]=0.0;J3[:]=0.0 #mandatory
 # # Hamiltonian splitting
 # integrate_Hp(pm,x1n,x2n,x3n,v1n,v2n,v3n,qn,mn,wn,
 #                   J1,J2,J3,B1,B2,B3,dt)
 #    J1[pm.N1n,pm.N2n,pm.N3n].=0
 #   J2[pm.N1n,pm.N2n,pm.N3n].=0
 #   J3[pm.N1n,pm.N2n,pm.N3n].=0
 #  E1.-=J1/Np;E2.-=J2/Np; E3.-=J3/Np

  for gdx=1:length(comp_gamma)
     dt2=comp_gamma[gdx]*dt
  #
  #
    integrate_H_E(pm,E1,E2,E3,B1,B2,B3,x1n,x2n,x3n,v1n,v2n,v3n,qn,mn,dt2/2.)
    integrate_H_B(pm, E1,E2,E3,B1,B2,B3,c,dt2/2.)

    J1[:]=0.0;J2[:]=0.0;J3[:]=0.0 #mandatory
   # Hamiltonian splitting
   #integrate_Hp_split_sym(pm,x1n,x2n,x3n,v1n,v2n,v3n,qn,mn,wn,
  #                   J1,J2,J3,B1,B2,B3,[dt2/2.],[dt2/2.])
   # Exponential boris - not symplectic, but good for strong field
    integrate_Hp_boris_exp(pm,x1n,x2n,x3n,v1n,v2n,v3n,qn,mn,wn,
                           J1,J2,J3,B1,B2,B3,[dt2/2],[dt2/2])

     # integrate_Hp_midpoint(pm,x1n,x2n,x3n,v1n,v2n,v3n,qn,mn,wn,
     #                  J1,J2,J3,B1,B2,B3,dt2)

    # Subtract steady background for momentum conservation
    # J1[pm.N1n,:,:].=0
    # J2[:,pm.N2n,:].=0
    # J3[:,:,pm.N3n].=0

     #
     J1[pm.N1n,pm.N2n,pm.N3n].=0
     J2[pm.N1n,pm.N2n,pm.N3n].=0
     J3[pm.N1n,pm.N2n,pm.N3n].=0
    E1.-=J1/Np;E2.-=J2/Np; E3.-=J3/Np

    integrate_H_B(pm, E1,E2,E3,B1,B2,B3,c,dt2/2.)
    integrate_H_E(pm,E1,E2,E3,B1,B2,B3,x1n,x2n,x3n,v1n,v2n,v3n,qn,mn,dt2/2.)

  end

  #x1n=mod.(x1n,L1);  x2n=mod.(x2n,L2);  x3n=mod.(x3n,L3)

end
toc()

# T=Float64;pdx=1;
#
# rho=accum_osde(pm,x1n,x2n,x3n, wn.*qn)./Np
# solve_poisson!(pm,rho,Phi)
# E1,E2,E3=gradient(pm,-Phi)
#
# J1[:]=0.0;J2[:]=0.0;J3[:]=0.0
#
 # integrate_Hp_midpoint(pm,x1n,x2n,x3n,v1n,v2n,v3n,qn,mn,wn,
                #  J1,J2,J3,B1,B2,B3,dt)
# #J1[pm.N1n,:,:].=0; J2[:,pm.N2n,:].=0; J3[:,:,pm.N3n].=0
# J1[pm.N1n,pm.N2n,pm.N3n].=0
# J2[pm.N1n,pm.N2n,pm.N3n].=0
# J3[pm.N1n,pm.N2n,pm.N3n].=0
# E1.-=J1/Np;E2.-=J2/Np; E3.-=J3/Np
# #
# # rho=accum_osde(pm,x1n,x2n,x3n, wn.*qn)./Np
# # solve_poisson!(pm,rho*gamma,Phi)
# # E1_,E2_,E3_=gradient(pm,-Phi)
# #
# # print(E1-E1_, "\n" )
# # print(J1/Np,"\n")
rho=accum_osde(pm,x1n,x2n,x3n, wn.*qn)./Np
print("Gauss Error: ",gauss_error(pm,rho,E1,E2,E3),"\n")

# Variance
# rho_var=accum_osde_var(pm,rho,x1n,x2n,x3n, wn.*qn)./Np
#
# rdx=sortperm(abs.(rho)[:])
#
# semilogy(abs.(rho[rdx]))
# plot(rho_var[rdx]/sqrt(Np))

energy=sum(Epot,2)+sum(Ekin,2)+sum(Bpot,2)

using PyPlot
figure()
semilogy(ttime,abs.((energy-energy[1])./energy[1]))
grid()
ylabel("rel. energy error")
xlabel("time ")

# plot(ttime, Ekin[:,1]+Epot[:,1])
#
figure()
 semilogy(ttime, Epot)
 figure()
 semilogy(ttime, Bpot)

# figure()
# semilogy(ttime, Bpot)

moment_error=abs.(Momentum[:,:].-reshape(Momentum[1,:],1,3));

figure()
semilogy(ttime,moment_error)
xlabel("time "); grid()
ylabel("absolute momentum error")

# print(test_eval_basis(pm,10000))
#
# function push_Hp(pm,x1n,x2n,x3n,v1n,v2n,v3n,qn,mn,wn,dt,J1,J2,J3,B1,B2,B3)
#   Np=length(x1n)
#   @inbounds begin
#   for pdx=1:Np
#
#    eval_basis(pm,x1n[pdx],x2n[pdx],x3n[pdx],psin0 )
#   #  (x3n[pdx],v1n[pdx],v2n[pdx])=push_Hp3(pm, x1n[pdx],x2n[pdx],x3n[pdx],
#   #                         v1n[pdx],v2n[pdx],v3n[pdx],
#   #                          qn[pdx],mn[pdx],wn[pdx],
#   #                          psin0,dt/2,J3, B1,B2 )
#   #   (x2n[pdx],v1n[pdx],v3n[pdx])=push_Hp2(pm, x1n[pdx],x2n[pdx],x3n[pdx],
#   #                          v1n[pdx],v2n[pdx],v3n[pdx],
#   #                           qn[pdx],mn[pdx],wn[pdx],
#   #                           psin0,dt/2,J2, B1,B3 )
#      (x1n[pdx],v2n[pdx],v3n[pdx])=push_Hp1(pm, x1n[pdx],x2n[pdx],x3n[pdx],
#                 v1n[pdx],v2n[pdx],v3n[pdx],
#                  qn[pdx],mn[pdx],wn[pdx],
#                  psin0,dt,J1, B2,B3 )
#     # (x2n[pdx],v1n[pdx],v3n[pdx])=push_Hp2(pm, x1n[pdx],x2n[pdx],x3n[pdx],
#     #                        v1n[pdx],v2n[pdx],v3n[pdx],
#     #                         qn[pdx],mn[pdx],wn[pdx],
#     #                         psin0,dt/2,J2, B1,B3 )
#     # (x2n[pdx],v1n[pdx],v3n[pdx])=push_Hp2(pm, x1n[pdx],x2n[pdx],x3n[pdx],
#     #                        v1n[pdx],v2n[pdx],v3n[pdx],
#     #                         qn[pdx],mn[pdx],wn[pdx],
#     #                         psin0,dt/2,J2, B1,B3 )
#
#   end
#   end
#
#
#
# end
# using BenchmarkTools

#tic()
# push_Hp(pm,x1n,x2n,x3n,v1n,v2n,v3n,qn,mn,wn,dt,J1,J2,J3,B1,B2,B3)

# @benchmark integrate_Hp(pm,x1n,x2n,x3n,v1n,v2n,v3n,qn,mn,wn,J1,J2,J3,B1,B2,B3,dt)
# @time integrate_Hp(pm,x1n,x2n,x3n,v1n,v2n,v3n,qn,mn,wn,J1,J2,J3,B1,B2,B3,dt)

#
#
# maximum(abs(rho-rho2))
#
#  solve_poisson!(pm,rho*gamma,Phi)
#  E1_,E2_,E3_=gradient(pm,-Phi)
# print(gauss_error(pm,rho,E1_,E2_,E3_),"\n")

# print(E1)
# print(E1_)
#J1[1,:,:]=0.0
# J2[:,p,:]=0.0
# dE1[1,pm.Nx2+1,pm.Nx3+1]=0.0

 # (E1_-E1)[:,pm.Nx2+1,pm.Nx3+1]

 # DE=J1[:,pm.Nx2+1,pm.Nx3+1]/Np
# J1

# (dE1/Np/pm.L1/pm.L2/pm.L3)[:,pm.Nx2+1,pm.Nx3+1]

# end
# tic()
# for tdx=1:Nt
#
#   @inbounds for rkdx=1:length(rksd)
#
#
#     rho=accum_osde(pm,x1n,x2n,x3n, wn.*qn)./Np
#
#     solve_poisson!(pm,rho*gamma,Phi)
#
#     E1,E2,E3=gradient(pm,-Phi)
#
#     if rkdx==1
#
#     #fieldenergy[tdx]=H1seminorm(pm,Phi)/2.
#     #kineticenergy[tdx]=sum( (v1n.^2 .+v2n.^2.+v3n.^2).*wn)/Np/2.
#     Epot[tdx,1]=L2norm(pm,E1)/2.
#     Epot[tdx,2]=L2norm(pm,E2)/2.
#     Epot[tdx,3]=L2norm(pm,E3)/2.
#     Ekin[tdx,1]=sum(v1n.^2.*wn)/Np/2.
#     Ekin[tdx,2]=sum(v2n.^2.*wn)/Np/2.
#     Ekin[tdx,3]=sum(v3n.^2.*wn)/Np/2.
#     Momentum[tdx,1]=sum(v1n.*wn)/Np
#     Momentum[tdx,2]=sum(v2n.*wn)/Np
#     Momentum[tdx,3]=sum(v3n.*wn)/Np
#     end
#
#
#     eval_vectorfield(pm,x1n,x2n,x3n,E1,E2,E3,E1n,E2n,E3n)
#     v1n.+=dt*rksc[rkdx]*real(E1n).*(qn./mn)
#     v2n.+=dt*rksc[rkdx]*real(E2n).*(qn./mn)
#     v3n.+=dt*rksc[rkdx]*real(E3n).*(qn./mn)
#
#     # eval_scalar(pm,x1n,x2n,x3n,E1,En)
#     # v1n.+=dt*rksc[rkdx]*real(En)
#     # eval_scalar(pm,x1n,x2n,x3n,E2,En)
#     # v2n.+=dt*rksc[rkdx]*real(En)
#     # eval_scalar(pm,x1n,x2n,x3n,E3,En)
#     # v3n.+=dt*rksc[rkdx]*real(En)
#
#     # v1n+=dt*rksc[rkdx]*eval_scalar(pm,x1n,x2n,x3n,E1*q/m)
#     # v2n+=dt*rksc[rkdx]*eval_scalar(pm,x1n,x2n,x3n,E2*q/m)
#     # v3n+=dt*rksc[rkdx]*eval_scalar(pm,x1n,x2n,x3n,E3*q/m)
#
#     x1n.+=dt*rksd[rkdx].*v1n
#     x2n.+=dt*rksd[rkdx].*v2n
#     x3n.+=dt*rksd[rkdx].*v3n
#
#
#   end
#   x1n=mod(x1n,L1)
#   x2n=mod(x2n,L2)
#   x3n=mod(x3n,L3)
#
# end
# toc()
#
# using PyPlot
# # semilogy(ttime,Epot)
#
#
# energy=gamma*sum(Epot,2)+sum(Ekin,2)
# moment_error=sum(abs(Momentum[:,:].-reshape(Momentum[1,:],1,3)),2);
#
# # E1n_=zeros(Float64,Np)
# # using BenchmarkTools
# #
# # @benchmark eval_scalar(pm,x1n,x2n,x3n,E1*q/m,E1n)
# # @benchmark eval_scalar2(pm,x1n,x2n,x3n,E1*q/m,E1n_)
# # print((E1n_[1]-real(E1n[1]))/real(E1n[1]) )
# figure()
# semilogy(ttime,abs((energy-energy[1])./energy[1]))
# grid()
# ylabel("rel. energy error")
# xlabel("time ")
# # savefig("$prefix\_energy_error.png")
# figure()
# semilogy(ttime,moment_error)
# xlabel("time "); grid()
# ylabel("absolute momentum error")
# # savefig("$prefix\_momentum_error.png")
#
#
# # figure()
# # plot(energy)
# #
# #
# # figure()
# # semilogy(moment_error)
# #
# # grid()
# #
# #
# figure()
# plot(abs(Phi[:]))
# # using PyPlot
# (L2norm(pm,E1)+L2norm(pm,E2)+L2norm(pm,E3))/2
#
# fieldenergy[end]
# Nviz=Np
# figure()
# sc=scatter3D(x1n[1:Nviz],x2n[1:Nviz],x3n[1:Nviz],s=2,c=fn[1:Nviz],depthshade="True",lw=0)
# colorbar(sc)
# tight()
# plot(fieldenergy)
#
# scatter(x1n,x2n,s=3,lw=0)
#
# figure()
# semilogy(ttime,fieldenergy)
# fig=figure()
# ax = gca(projection="3d")
#
# ax[:scatter](x1n,x2n,x3n)
# @time eval_scalar(pm,x1n,x2n,x3n,E1*q/m,En)

# figure()
# plot(x1n,x2n,x3n)

#@time rho=rhs_particle(pm,x1n,x2n,x3n, wn)./Np*q

#Damping rate for linear Landau
#omega=1.415661888604536 - 0.153359466909605im;
#semilogy(ttime[ttime.<10.],0.5*fieldenergy[1]*abs(exp(-im*omega*(ttime[ttime.<10.]-0.4))).^2)

#plot(abs(rho)[:])

#rho

# plot(abs(Phi[:]))

# type vlasov6d
#
#   Np::Int #Number of particles
#   Ns::Int #Number of species
#
#   x1n::Array{Float64,1}
#   x2n::Array{Float64,1}
#   x3n::Array{Float64,1}
#
#   v1n::Array{Float64,1}
#   v2n::Array{Float64,1}
#   v3n::Array{Float64,1}
#   qn::Array{Float64,1} # Charge
#   mn::Array{Float64,1} # Mass
#
#   fn::Array{Float64,1}
#   gn::Array{Float64,1}
#
#   function vlasov6d(Np,Ns)
#     x=zeros(Float64,Int(Np)*Int(Ns))
#     new(Int(Np),Int(Ns),x,x,x, x,x,x, x,x,x,x)
#
#   end
# end
#
#
# function std_Maxwellian(x1,x2,x3,v1,v2,v3)
#   return exp(-(v1.^2+v2.^2+v3.^2)/2 )/(sqrt(2*pi))^3
# end
#
#
#
# function sample_std_Maxwellian(sim::vlasov6d, L1,L2,L3   )
#   sim.x1n=rand(Float64,sim.Np)*L1
#   sim.x2n=rand(Float64,sim.Np)*L2
#   sim.x3n=rand(Float64,sim.Np)*L3
#
#   sim.v1n=randn(Float64,sim.Np)
#   sim.v2n=randn(Float64,sim.Np)
#   sim.v3n=randn(Float64,sim.Np)
#
#   gn=std_Maxwellian(sim.x1n,sim.x2n,sim.x3n,sim.v1n,sim.v2n,sim.v3n)/L1/L2/L3
# end
#
#
# f0=(x1,x2,x3,v1,v2,v3)->(exp(-(v1.^2+v2.^2+v3.^2)/2 )/(sqrt(2*pi))^3
#
#
# #Electron only
# sim=vlasov6d(1e5,1)
# k=0.5;
# epse=0.1
# #Particle mesh coupling
# pm=pif3d(2,2,2,2*pi/k)
#
# sim.qn[:]=1.;sim.mn[:]=1.
# sample_std_Maxwellian(sim,pm.L1,pm.L2,pm.L3)
# sim.fn=(1+epse*cos(k*x)  std_Maxwellian(sim.x1n,sim.x2n,sim.x3n,sim.v1n,sim.v2n,sim.v3n)
#
# rhs=rhs.particle(sim.x1n,sim.x2n,sim.x3n, sim.fn/sim.gn) /Np*Ns
# sim.x1n

# using BenchmarkTools
# @benchmark rho=accum_osde(pm,x1n,x2n,x3n, wn.*qn)./Np
# @benchmark rho=accum_osde2(pm,x1n,x2n,x3n, wn.*qn)./Np

max_stars_repo_name: ameresj/VlasovPIF

max_issues_repo_name: ameresj/VlasovPIF

max_forks_repo_name: ameresj/VlasovPIF


text:
function convert_image_files_to_weights(bus1, bus2)

    # create dictionary with rbg rgb_values
    rgb_values = Dict{String, Any}()
    boundaries = Dict{String, Any}()
    nodes_lp = Dict{String, Any}()
    plot_dictionary = Dict{String, Any}()

    #To Do add rnage based on input location
    img_ag = Images.load("../src/spatial_image_files/clc_agricultural.tif")
    img_overlay = zeros(size(img_ag))
    rgb_values, img_overlay = convert2integer(img_ag, rgb_values, "agricultural", 1, img_overlay)
    img_sea = Images.load("../src/spatial_image_files/sea.tif")
    rgb_values, img_overlay = convert2integer(img_sea, rgb_values, "sea", 2, img_overlay)
    img_grid = Images.load("../src/spatial_image_files/grid_225kv_400kv.tif")
    rgb_values, img_overlay = convert2integer(img_grid, rgb_values, "grid", 3, img_overlay)
    img_roads = Images.load("../src/spatial_image_files/roads.tif")
    rgb_values, img_overlay = convert2integer(img_roads, rgb_values, "roads", 3, img_overlay)
    img_railroads = Images.load("../src/spatial_image_files/railroads.tif")
    rgb_values, img_overlay = convert2integer(img_railroads, rgb_values, "railroads", 3, img_overlay)
    img_nature = Images.load("../src/spatial_image_files/clc_natural.tif")
    rgb_values, img_overlay = convert2integer(img_nature, rgb_values, "nature", 4, img_overlay)
    img_mountains = Images.load("../src/spatial_image_files/clc_mountains.tif")
    rgb_values, img_overlay = convert2integer(img_mountains, rgb_values, "mountain", 5, img_overlay)
    img_urban = Images.load("../src/spatial_image_files/clc_urban.tif")
    rgb_values, img_overlay = convert2integer(img_urban, rgb_values, "urban", 6, img_overlay)
    img_natura2000 = Images.load("../src/spatial_image_files/Natura2000.tif")
    rgb_values, img_overlay = convert2integer(img_natura2000, rgb_values, "Natura2000", 7, img_overlay)

    X1,Y1 = convert_latlon_to_etrs89(bus1["longitude"], bus1["latitude"])
    X2,Y2 = convert_latlon_to_etrs89(bus2["longitude"], bus2["latitude"])

    # http://www.eea.europa.eu/data-and-maps/data/corine-land-cover-2006-raster-2
    x0=1500000;
    xmax=7400000;
    y0=900000;
    ymax=5500000;
    resolution=2500;

    # Calculate the positions pof the nodes
    xpos=round.([(X1-x0),(X2-x0)] / resolution)
    ypos=round.([(ymax-Y1),(ymax-Y2)] / resolution)

    # Determine range: Euclidian distance +/- 30%
    d = round(0.3 * sqrt((xpos[1] - xpos[1])^2 + (ypos[1] - ypos[2])^2))
    x_min = min(xpos[1], xpos[2]) - d
    x_max = max(xpos[1], xpos[2]) + d
    y_min = min(ypos[1], ypos[2]) - d
    y_max = max(ypos[1], ypos[2]) + d

    boundaries["xmin"] = x_min
    boundaries["xmax"] = x_max
    boundaries["ymin"] = y_min
    boundaries["ymax"] = y_max

    xpos_plot=round.([(X1-x0); (X2-x0)]/resolution .- x_max) #range(1,1);
    ypos_plot=round.([(ymax-Y1); (ymax-Y2)]/resolution .- y_min) #- range(2,1);

    nodes_lp["x1"] = xpos[1]
    nodes_lp["x2"] = xpos[2]
    nodes_lp["y1"] = ypos[1]
    nodes_lp["y2"] = ypos[2]

    plot_dictionary["overlay_image"] = img_overlay
    plot_dictionary["x_position"] = xpos
    plot_dictionary["y_position"] = ypos

    return rgb_values, nodes_lp, boundaries, plot_dictionary
end

function convert2integer(image, rgb_values, imagename, factor, img_overlay)
    rw = Images.rawview(Images.channelview(image))
    bitarray = (rw  .!==  0x00000000)
    rgb_values[imagename] = ones(size(bitarray)) .* bitarray
    img = findall(x->x==1, rgb_values[imagename])

    for idx = 1:length(img)
        img_overlay[img[idx]] = max(factor,  img_overlay[img[idx]])
    end

    return rgb_values, img_overlay
end

function convert_latlon_to_etrs89(longitude, latitude)

    # Coordinate transformation based on IOGP Geomatics
    # Guidance Note Number 7, part 2
    # Coordinate Conversions and Transformations including Formulas
    # Section 3.4.2 Lambert Azimuthal Equal Area (EPSG Dataset coordinate operation method code 9820)
    # https://www.iogp.org/bookstore/product/coordinate-conversions-and-transformation-including-formulas/

    # fixed parameters
    a = 6378137.0
    e = 0.081819191
    phi0 = 0.907571211 # rad = 52 degree North
    lambda0 = 0.174532925 # rad = 10 degree East
    FE = 4321000 # False easting in m
    FN = 3210000 # False nortinh in m

    phi = latitude * pi / 180
    lambda = longitude * pi / 180

    q = (1 - e^2) * ((sin(phi) / (1 - e^2 * sin(phi)^2)) - ((1 / (2 * e)) * log((1 - e * sin(phi)) / (1 + e * sin(phi)))))
    q0 = (1 - e^2) * ((sin(phi0) / (1 - e^2 * sin(phi0)^2)) - ((1 / (2 * e)) * log((1 - e * sin(phi0)) / (1 + e * sin(phi0)))))
    qp = (1 - e^2) * ((1 / (1 - e^2)) - ((1 / (2 * e) * log((1 - e) / (1 + e)))))

    Rq = a * (qp / 2)^0.5
    beta = asin(q / qp)
    beta0 = asin(q0 / qp)

    B = Rq * (2 / (1+ sin(beta0) * sin(beta) + (cos(beta0) * cos(beta) * cos(lambda - lambda0))))^0.5
    D = a * (cos(phi0) / (1 - e^2 * sin(phi0)^2)^0.5) / (Rq * cos(beta0))

    E = FE + (B*D) * (cos(beta) * sin(lambda - lambda0))
    N = FN + (B/D) * (cos(beta0) * sin(beta) - (sin(beta0) * cos(beta) * cos(lambda - lambda0)))
    return E, N
end

ext: jl


text:
n,m = 3,2

# Custom Equality Constraint
p1 = 3
c(v,x,u) = begin v[1]  = x[1]^2 + x[2]^2 - 5; v[2:3] =  u - ones(2,1) end
c(v,x) = c(v,x,zeros(2))
jacob_c(x,u) = [2x[1] 2x[2] 0 0 0;
                0     0     0 1 0;
                0     0     0 0 1];
jacob_c(x) = [2x[1] 2x[2] 0;
                0     0     0;
                0     0     0];

v = zeros(p1)
x = [1.,2,3]
u = [-5.,5]
c(v,x,u)
@test v == [0,-6,4]

c(v,x)
@test v == [0, -1, -1]

J = zeros(p1,n+m)
J_term = zeros(p1,n)
ForwardDiff.jacobian!(J_term,c,v,x)
@test J_term == jacob_c(x)

import TrajectoryOptimization: generate_jacobian
∇c, = generate_jacobian(c, n, m, p1)
methods(∇c)
∇c(J,x,u)
@test J == jacob_c(x,u)
@test ∇c(J_term,x) == jacob_c(x)

# Test constraint function
con = Constraint{Equality}(c,n,m,p1,:custom)
evaluate!(v, con, x, u)
@test v == [0,-6,4]
evaluate!(v, con, x)
@test v == [0, -1, -1]
@test length(con) == p1

# Test constraint jacobian
con.∇c(J,x,u);
jacobian!(J, con, x, u)
@test J == jacob_c(x,u)
jacobian!(J_term, con, x)
@test J_term == jacob_c(x)

# Custom inequality constraint
p2 = 2
c2(v,x,u) = begin v[1] = sin(x[1]); v[2] = sin(x[3]) end
∇c2(Z,x,u) = begin Z[1,1] = cos(x[1]); Z[2,3] = cos(x[3]); end
con2 = Constraint{Inequality}(c2,∇c2,n,m,p2,:ineq)

# Bound constraint
x_max = [5,5,Inf]
x_min = [-10,-5,0.]
u_max = 0.
u_min = -10.
#
# p3 = 2(n+m)
# bnd = BoundConstraint(n,m,x_max=x_max,x_min=x_min,u_min=u_min,u_max=u_max, trim=false)
# v = zeros(p3)
# evaluate!(v, bnd, x, u)
# @test v == [-4,-3,-Inf,-5,5,-11,-7,-3,-5,-15]
# C = PartedArray(zeros(p3,n+m),create_partition2((p3,),(n,m),Val((:xx,:xu))))
# jacobian!(C, bnd, x, u)
# @test C.xx == [Diagonal(1I,n); zeros(m,n); -Diagonal(1I,n); zeros(m,n)]
# @test C.xu == [zeros(n,m); Diagonal(1I,m); zeros(n,m); -Diagonal(1I,m)]
# @test length(bnd) == p3
#
# v = zeros(2n)
# evaluate!(v, bnd, x)
# @test v == [-4,-3,-Inf,-11,-7,-3]
# C = zeros(2n, n)
#
# C = zeros(2n-1,n)
# bnd.active.x_all
# bnd.jac[bnd.active.x_all,1:n]
# bnd.active.x_min
# jacobian!(C, bnd, x)
# @test jacobian!(C, bnd, x) == [Diagonal(1.0I,n)[1:(n-1),:]; -Diagonal(1.0I,n)]
# @test length(bnd,:terminal) == 2n

# Trimmed bound constraint
bnd = BoundConstraint(n,m,x_max=x_max,x_min=x_min,u_min=u_min,u_max=u_max)
p3 = 2(n+m)-1
p3_N = 2n-1
v = zeros(p3)
evaluate!(v, bnd, x, u)
@test v == [-4,-3,-5,5,-11,-7,-3,-5,-15]
C = zeros(p3, n+m)
jac_bnd = [Diagonal(1I,n+m)[[true;true;false;true;true],:]; -Diagonal(1I,n+m)]
@test jacobian!(C, bnd, x, u) == jac_bnd
@test length(bnd,:stage) == p3

v = zeros(p3_N)
evaluate!(v, bnd, x)
@test v == [-4,-3,-11,-7,-3]
C = zeros(p3_N, n)
jacobian!(C, bnd, x)
jac_bnd_term = [Diagonal(1I,n)[[true;true;false],:]; -Diagonal(1I,n)]
@test C == jac_bnd_term

# Create Constraint Set
C = [con,con2,bnd]
@test C isa ConstraintSet
C = con + con2 + bnd
@test C isa ConstraintSet

@test findall(C,Inequality) == [false,true,true]
@test split(C) == ([con2,bnd],[con,])
@test count_constraints(C) == (p2+p3,p1)
@test inequalities(C) == [con2,bnd]
@test equalities(C) == [con,]
@test bounds(C) == [bnd,]
@test labels(C) == [:custom,:ineq,:bound]
@test TrajectoryOptimization.is_terminal.(C) == [true,false,true]

# Terminal Constraint
cterm(v,x) = begin v[1] = x[1] - 5; v[2] = x[1]*x[2] end
∇cterm(x) = [1 0 0; x[2] x[1] 0]
∇cterm(A,x) = copyto!(A,∇cterm(x))
p_term = 2
v = zeros(p_term)
cterm(v,x)
con_term = Constraint{Equality}(cterm,∇cterm,n,p_term,:terminal)
v2 = zeros(p_term)
con_term.c(v2,x)
@test v == v2
A = zeros(p_term,n)
@test con_term.∇c(A,x) == ∇cterm(x)

C_term = [con,bnd,con_term]
C2 = [con,con2,bnd,con_term]
@test C2 isa ConstraintSet
@test terminal(C2) == C_term
@test terminal(C_term) == C_term
@test stage(C2) == C
count_constraints(C2)
@test split(C2) == ([con2,bnd],[con,con_term])
@test split(C2) == (inequalities(C2),equalities(C2))
@test bounds(C2) == [bnd]
@test labels(C2) == [:custom,:ineq,:bound,:terminal]
@test TrajectoryOptimization.num_stage_constraints(C2) == p1+p2+p3
@test TrajectoryOptimization.num_terminal_constraints(C2) == p3_N + p_term + p1

import TrajectoryOptimization: num_stage_constraints, num_terminal_constraints
num_stage_constraints(C2)
num_constraints(C2,:stage)

v_stage = PartedVector(C2)
v_term = PartedVector(C2,:terminal)
v_stage2 = PartedVector(C2,:stage)
v_term2 = PartedVector(C2,:terminal)
evaluate!(v_stage,C,x,u)
evaluate!(v_term,C_term,x)
evaluate!(v_stage2,C2,x,u)
evaluate!(v_term2,C2,x)
@test v_stage == v_stage2
@test v_term == v_term2

c_jac = PartedMatrix(C,n,m)
@test size(c_jac) == (p1+p2+p3,n+m)
@test size(c_jac.custom) == (p1,n+m)
@test size(c_jac.x) == (p1+p2+p3,n)
@test size(c_jac.u) == (p1+p2+p3,m)
jacobian!(c_jac, C, x, u)
@test c_jac.custom == jacob_c(x,u)
@test c_jac.bound == jac_bnd

c_jac = PartedMatrix(C2,n,m,:terminal)
@test size(c_jac) == (p1+p3_N+p_term,n)
@test size(c_jac.custom) == (p1,n)
@test size(c_jac.ineq) == (0,n)
@test size(c_jac.u) == (p1+p3_N+p_term,0)
@test size(c_jac.terminal) == (p_term,n)
jacobian!(c_jac, C2, x)
jacobian!(c_jac.custom, con, x)
@test c_jac.custom == jacob_c(x)
@test c_jac.bound == jac_bnd_term

# Augment State
m_inf = m+n
con_inf = TrajectoryOptimization.infeasible_constraints(n,m)
u_inf = [u; 5; -5; 10]
v = zeros(n)
con_inf.c(v,x,u_inf[con_inf.inds[2]])
@test v == [5,-5,10]

C_inf = [con,con2,bnd,con_inf]
v_stage = PartedVector(stage(C))
v_inf = PartedVector(stage(C_inf))
TrajectoryOptimization.evaluate!(v_stage,C,x,u_inf)
TrajectoryOptimization.evaluate!(v_inf,C_inf,x,u_inf)
@test v_inf == [v_stage;5;-5;10]

C2
C3 = copy(C2)
@test TrajectoryOptimization.remove_bounds!(C3) == [bnd]
@test C3 == [con,con2,con_term]
@test pop!(C3,:ineq) == con2
@test C3 == [con,con_term]

# Constraints
CS = [[con], con+con2+bnd, con+con2+bnd, con+bnd]
N = length(CS)
@test CS isa Vector{<:ConstraintSet}
PC = Constraints(CS)
PC[1] += bnd
@test PC[1] == [con,bnd]
@test num_constraints(PC) == [p1+p3, p1+p2+p3, p1+p2+p3, p1+p3_N]
PC[N] += con_term
@test num_constraints(PC) == [p1+p3, p1+p2+p3, p1+p2+p3, p1+p3_N+p_term]
PC2 = copy(PC)
pop!(PC2[2],:ineq)
@test num_constraints(PC2) == [p1+p3, p1+p3, p1+p2+p3, p1+p3_N+p_term]
@test num_constraints(PC) == [p1+p3, p1+p2+p3, p1+p2+p3, p1+p3_N+p_term]

Cval, ∇Cval = TrajectoryOptimization.init_constraint_trajectories(PC,n,m,N)
@test length.(Cval) == num_constraints(PC)
X = [float.(x) for k = 1:N]
U = [float.(u) for k = 1:N-1]
TrajectoryOptimization.update_constraints!(Cval, PC, X, U)
v = zeros(p1)
c(v,x,u)
@test Cval[1].custom == v
c(v,x)
@test Cval[N].custom == v

jacobian!(∇Cval, PC, X, U)
@test ∇Cval[N] == c_jac

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
#=
Types and basic constructors for phylogenetic trees
=#

abstract type AbstractTree end
abstract type AbstractNode end
abstract type AbstractBranch end

const NodeVector = Vector{N} where N <: AbstractNode
const TreeVector = Vector{T} where T <: AbstractTree

struct NodeLink
    to::AbstractNode
    branch::AbstractBranch
end # struct NodeLink

mutable struct Branch <: AbstractBranch
    length::Union{Float64,Nothing}
    label::String
    annotations::Dict{String,Any}
    datablocks::Vector{AbstractBranchDataBlock}
    bipart::Bipartition

    Branch() = new(nothing, "", Dict{String,Any}(), AbstractBranchDataBlock[], Bipartition(falses(1)))
end # struct Branch

function Branch(brlength::Union{Nothing, Float64})
    br = Branch()
    br.length = brlength

    return br
end # function Branch

mutable struct Node <: AbstractNode
    label::String
    links::Vector{NodeLink}
    species::Int
    annotations::Dict{String,Any}
    datablocks::Vector{AbstractNodeDataBlock}

    function Node()
        links = NodeLink[]
        sizehint!(links, 3)
        new("", links, 0, Dict{String,Any}(), AbstractNodeDataBlock[])
    end
end # struct Node

function Node(l::String)
    p = Node()
    p.label = l

    return p
end # function Node

mutable struct Tree <: AbstractTree
    origin::AbstractNode
    label::String
    rooted::Bool
    dir::SpeciesDirectory
    annotations::Dict{String,Any}
    datablocks::Vector{AbstractTreeDataBlock}

    Tree(p::AbstractNode) =
        new(p, "", false, SpeciesDirectory(), Dict{String,Any}(), AbstractTreeDataBlock[])
end # struct Tree

# mutable struct FixedTree
#     start::Node
#     label::String
#     rooted::Bool
#     preorder::Matrix{Node}
# end


text:
#=
Author:= Shreeda Bhat
License:= MIT
UTF-8 =#
println("*********************** 1st Example *************************")
println(methods(+))

println("*********************** 3rd  Example *************************")
import Base:+

+(x::String,y::String)=string(x,y)
"hello"+"world!"


text:
"""
```
struct MMALagModel <: AbstractModel
    parent::VecModel
    quadweight::Base.RefValue{<:Real}
    linweights::AbstractVector
    origobjval::Base.RefValue{<:Real}
    origconstrval::AbstractVector
end
```

A model with aggregated block constraints. Any block constraint in `parent` is aggregated using the corresponding weights in `linweights`. Non-block constraints are treated as block constraints of dimension 1. Let `violations` be the sum of constraint violations at a point `x`. The objective value of an instance of `MMALagModel` is the original objective value of its `parent` plus `quadweight[] * sum(max.(violations, 0)^2)`. The Lagrangian function of the new problem is the augmented Lagrangian function of the original problem. `origobjval` and `origconstrval` store the original objective and constraint values respectively when evaluating the aggregated objective and constraints.
"""
@params struct MMALagModel <: AbstractModel
    parent::VecModel
    quadweight::Base.RefValue{<:Real}
    linweights::AbstractVector
    origobjval::Base.RefValue{<:Real}
    origconstrval::AbstractVector
end
function MMALagModel(
    model::VecModel;
    quadweight = 1e-5,
    linweights = ones(2*getdim(getineqconstraints(model))),
    kwargs...,
)
    return MMALagModel(
        model,
        Ref(quadweight),
        linweights,
        Ref(zero(quadweight)),
        zeros(getdim(getineqconstraints(model))),
    )
end
MMALagModel(args...; kwargs...) = MMALagModel(Model(args...); kwargs...)

getparent(m::MMALagModel) = m.parent

function getobjectiveconstraints(m::MMALagModel)
    parent = getparent(m)
    quadpenalty = NonNegSumOfSquares(m.quadweight[])

    dims = getdim.(getineqconstraints(getparent(m)).fs)
    M = sum(dims)
    nblocks = length(dims)
    func = x -> begin
        origobjval = getobjective(parent)(x)
        origconstrval = getineqconstraints(parent)(x)
        setvals!(m, origobjval / get_objective_multiple(parent), origconstrval)

        obj = origobjval + quadpenalty(origconstrval)
        start = 1
        constr1 = map(1:nblocks) do i
            inds = start:(start+dims[i]-1)
            start += dims[i]
            return @views WeightedSum(m.linweights[inds])(origconstrval[inds])
        end
        constr2 = map(1:nblocks) do i
            inds = start:(start+dims[i]-1)
            start += dims[i]
            return @views WeightedSum(m.linweights[inds])(origconstrval[inds .- M])
        end
        @assert start == 2*length(origconstrval) + 1
        return [obj; constr1; constr2]
    end
    return FunctionWrapper(func, 2nblocks + 1)
end

function getobjective(m::MMALagModel)
    objf(x) = getobjectiveconstraints(m)(x)[1]
    return Objective(objf)
end

function getineqconstraint(m::MMALagModel, i::Integer)
    throw("`getineqconstraint` is not defined for `MMALagModel`.")
end

function geteqconstraint(m::MMALagModel, i::Integer)
    throw("`geteqconstraint` is not defined for `MMALagModel`.")
end

function getineqconstraints(m::MMALagModel)
    dim = getdim(getobjectiveconstraints(m)) - 1
    constrf(x) = getobjectiveconstraints(m)(x)[2:end]
    return FunctionWrapper(constrf, dim)
end

function geteqconstraints(m::MMALagModel)
    throw("`geteqconstraint` is not defined for `MMALagModel`.")
end

get_objective_multiple(model::MMALagModel) = get_objective_multiple(getparent(model))

function set_objective_multiple!(model::MMALagModel, m::Real)
    set_objective_multiple!(getparent(model), m)
    return model
end

getmin(m::MMALagModel)= getmin(m.parent)

getmax(m::MMALagModel) = getmax(m.parent)

function addvar!(m::MMALagModel, lb, ub; kwargs...)
    addvar!(getparent(m), lb, ub; kwargs...)
    return m
end

function add_ineq_constraint!(m::MMALagModel, f::IneqConstraint)
    add_ineq_constraint!(getparent(m), f)
    append!(m.linweights, ones(2*getdim(f)))
    append!(m.origconstrval, zeros(getdim(f)))
    return m
end

function add_ineq_constraint!(m::MMALagModel, fs::Vector{<:IneqConstraint})
    add_ineq_constraint!(getparent(m), fs)
    append!(m.linweights, ones(2*getdim(VectorOfFunctions(fs))))
    append!(m.origconstrval, zeros(getdim(VectorOfFunctions(fs))))
    return m
end

"""
    getquadweight(m::MMALagModel)

Returns the quadratic penalty of `m`.
"""
getquadweight(m::MMALagModel) = m.quadweight[]

"""
    getlinweights(m::MMALagModel)

Returns the linear weights of `m`. These are the weights used to aggregate the block constraints in the `MMALag` algorithm.
"""
getlinweights(m::MMALagModel) = m.linweights

"""
    getorigobjval(m::MMALagModel)

Returns the last cached objective value.
"""
getorigobjval(m::MMALagModel) = m.origobjval[]

"""
    getorigconstrval(m::MMALagModel)

Returns the last cached constraint values.
"""
getorigconstrval(m::MMALagModel) = m.origconstrval

"""
    setvals!(m::MMALagModel, obj, constr)

Caches the objective and constraint values `obj` and `constr` in `m`.
"""
function setvals!(m::MMALagModel, obj, constr)
    m.origobjval[] = obj
    m.origconstrval .= constr
    return m
end
ChainRulesCore.@non_differentiable setvals!(m::MMALagModel, obj, constr)

"""
    setquadweight!(m::MMALagModel, w)

Sets the quadratic penalty in `m` to `w`.
"""
function setquadweight!(m::MMALagModel, w)
    m.quadweight[] = w
    return m
end

"""
    setlinweights!(m::MMALagModel, λ)

Sets linear weights in `m` to `λ`.
"""
function setlinweights!(m::MMALagModel, λ)
    m.linweights .= λ
    return m
end


text:
# This file was generated, do not modify it. # hide
knn = machine(knn_model, X, y)


text:
@Base.kwdef mutable struct Supplier
    name::String
    cost_per_cut::MoneyType
    kerf::LengthType
    available_stock::Vector{AvailablePanel}
end

export Supplier


text:
__precompile__()

module Pages

using HTTP, JSON

import HTTP.WebSockets.WebSocket

export Endpoint, Callback, Plotly

mutable struct Endpoint
    handler::Function
    route::String
    sessions::Dict{String,WebSocket}

    function Endpoint(handler,route)
        p = new(handler,route,Dict{String,WebSocket}())
        !haskey(pages,route) || warn("Page $route already exists.")
        pages[route] = p
        finalizer(p, p -> delete!(pages, p.route))
        p
    end
end
function Base.show(io::Base.IO,endpoint::Endpoint)
    print(io,"Endpoint created at $(endpoint.route).")
end
const pages = Dict{String,Endpoint}() # url => page

include("callbacks.jl")
include("server.jl")
include("api.jl")
include("displays/plotly.jl")
# include("ijulia.jl")

include("../examples/examples.jl")

end


text:
using FactCheck

include("unknowns.jl")

exitstatus()


text:
input_data_path = "/home/user/Downloads/raw_data_old"
output_data_path = "/home/user/Downloads/raw_data"

input_dir = joinpath(input_data_path, "csv");
output_dir = joinpath(output_data_path, "csv");

input_files = joinpath.(input_dir, readdir(input_dir));
output_files = joinpath.(output_dir, readdir(input_dir));

# %%

for i in 1:length(input_files)
	open(input_files[i], "r") do f
		data = read(f, String)
		new_data = replace(data,
			"graph_type,evolution_method,n,t,seed,r_0,r_1,m_0,m_1" =>
			"graph_type,evolution_method,N,t,seed,r_0,r_1,m_0,m_1")
		open(output_files[i], "w") do g
			write(g, new_data)
		end
	end
end

# %%

for t in ["t_0", "t_1"]
	input_dir = joinpath(input_data_path, "cluster_size_distribution", t);
	output_dir = joinpath(output_data_path, "cluster_size_distribution", t);

	input_files = joinpath.(input_dir, readdir(input_dir));
	output_files = joinpath.(output_dir, readdir(input_dir));

	for i in 1:length(input_files)
		open(output_files[i], "w") do file
			for (j, line) in enumerate(eachline(input_files[i]))
				if j == 1
					write(file, string(replace(line, "# n = " => "# N = "), "\n"))
					a    = maximum(findfirst("_seed_", input_files[i])) + 1
					b    = minimum(findfirst(".csv", input_files[i])) - 1
					seed = parse(Int, input_files[i][a:b])
					write(file, "# seed = $(seed)\n")
					continue
				end
				write(file, string(line, "\n"))
			end
		end
	end
end

max_stars_repo_name: cameronperot/explosive-percolation

max_issues_repo_name: cameronperot/explosive-percolation

max_forks_repo_name: cameronperot/explosive-percolation


text:
#=
mask.jl
Methods related to a image support mask:
	embed, maskit, mask_or, mask_outline
2019-06-22 Jeff Fessler, University of Michigan
=#

export embed, mask_or, mask_outline, maskit, mask_test

using Test: @test, @inferred
using ImageFiltering: imfilter, centered
#using SparseArrays: sparse, findnz, AbstractSparseVector

"""
    mask_or(mask)
compress 3D mask to 2D by logical `or` along `z` direction
"""
function mask_or(mask::AbstractArray{Bool})
	ndim = ndims(mask)
	return	ndim == 2 ? mask :
			ndim == 3 ? dropdims(sum(mask, dims=3) .> 0, dims=3) :
			throw("ndims(mask) = $ndim")
end

"""
    mask_or(:test)
self test
"""
function mask_or(test::Symbol)
	test != :test && throw(ArgumentError("test $test"))
	mask2 = trues(3,4)
	@test (@inferred mask_or(mask2)) == mask2
	mask3 = trues(3,4,5)
	@test (@inferred mask_or(mask3)) == trues(3,4)
	true
end

"""
    mask_outline(mask)
return outer boundary of 2D mask (or mask_or for 3D)
"""
function mask_outline(mask::AbstractArray{Bool})
	mask2 = mask_or(mask)
	tmp = imfilter(mask2, centered(ones(Int32,3,3))) # dilate
#	tmp = tmp[2:end-1,2:end-1] # 'same'
	return (tmp .> 1) .& (.! mask2)
end

"""
    mask_outline(:test)
self test
"""
function mask_outline(test::Symbol)
	test != :test && throw(ArgumentError("test $test"))
	mask2 = trues(3,4)
	@test (@inferred mask_outline(mask2)) == falses(3,4)
	true
end

"""
    array = embed(v, mask)

embed vector `v` of length `sum(mask)`
into elements of an array where `mask` is `true`
"""
function embed(v::AbstractVector{<:Number}, mask::AbstractArray{Bool})
	array = zeros(eltype(v), size(mask))
	array[mask] .= v
	return array
end

"""
    array = embed(matrix::AbstractMatrix{<:Number}, mask::AbstractArray{Bool})

Embed each column of `matrix` into `mask` then `cat` along next dimension
In:
* `matrix [sum(mask) L]`
* `mask [(N)]`

Out:
* `array [(N) L]`
"""
function embed(x::AbstractMatrix{<:Number}, mask::AbstractArray{Bool})
	L = size(x,2)
	out = zeros(eltype(x), prod(size(mask)), L)
	for l=1:L
		out[:,l] = vec(embed(x[:,l], mask))
	end
	reshape(out, size(mask)..., L)
end

#=
this is needed only for "unpacking" sparse system matrices so ignore for now
# image_geom_embed_sparse
# called by image_geom_embed
# function _embed_sparse(x::Array{T}) where {T <: Number}
function embed(x::AbstractSparseVector{<:Number},
				mask::AbstractArray{Bool,N} where N)
	i, v = findnz(x)
	ind = findall(mask)
	j = ind(j)
	return sparsevec(i, j, a, size(x,1), length(mask))
end
=#

"""
    embed(:test)
"""
function embed(test::Symbol)
	test != :test && throw(ArgumentError("test $test"))
	mask = [false true true; true false false]
	a = [0 2 3; 1 0 0]; v = a[vec(mask)]
	@test (@inferred embed(v,mask)) == a
	@test (@inferred embed([v 2v], mask)) == cat(dims=3, a, 2a)
#	@test embed(sparse(1:3),mask) == sparse([0 2 3; 1 0 0]) # later
	true
end

"""
    maskit(x::AbstractArray{<:Number})
opposite of embed
"""
function maskit(x::AbstractArray{<:Number}, mask::Array{Bool})
	dim = size(x)
	x = reshape(x, length(mask), :)
	x = x[vec(mask),:] # reshape(mask, prod(_dim()))]
	if length(dim) == ndims(mask)
		x = dropdims(x, dims=2) # squeeze
	elseif length(dim) > ndims(mask)
		x = reshape(x, :, dim[(1+ndims(mask)):end])
	else
		throw(DimensionMismatch("size(x) = $(size(x))"))
	end
	return x
end

"""
    maskit(:test)
"""
function maskit(test::Symbol)
	test != :test && throw(ArgumentError("test $test"))
	mask = [false true true; true false false]
	@inferred maskit([7 2 3; 1 7 7], mask)
	@test maskit([7 2 3; 1 7 7], mask) == [1,2,3]
	true
end

"""
    mask_test()
self tests
"""
function mask_test()
	@test mask_or(:test)
	@test mask_outline(:test)
	@test embed(:test)
	@test maskit(:test)
	true
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct TensorDescription <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function TensorDescription(; kwargs...)
        obj = new(meta(TensorDescription), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct TensorDescription
const __meta_TensorDescription = Ref{ProtoMeta}()
function meta(::Type{TensorDescription})
    ProtoBuf.metalock() do
        if !isassigned(__meta_TensorDescription)
            __meta_TensorDescription[] = target = ProtoMeta(TensorDescription)
            fnum = Int[1,2,4]
            allflds = Pair{Symbol,Union{Type,String}}[:dtype => Int32, :shape => TensorShapeProto, :allocation_description => AllocationDescription]
            meta(target, TensorDescription, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_TensorDescription[]
    end
end
function Base.getproperty(obj::TensorDescription, name::Symbol)
    if name === :dtype
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :shape
        return (obj.__protobuf_jl_internal_values[name])::TensorShapeProto
    elseif name === :allocation_description
        return (obj.__protobuf_jl_internal_values[name])::AllocationDescription
    else
        getfield(obj, name)
    end
end

export TensorDescription


text:
module BostonHousing_Tests
using Test
using DataDeps
using MLDatasets

@testset "Boston Housing" begin
    X  = BostonHousing.features()
    Y  = BostonHousing.targets()
    names = BostonHousing.feature_names()
    @test X isa Matrix{Float64}
    @test Y isa Matrix{Float64}
    @test names == ["crim","zn","indus","chas","nox","rm","age","dis","rad","tax","ptratio","b","lstat"]
    @test size(X) == (13, 506)
    @test size(Y) == (1, 506)
end

end #module


text:
export EulerSymp
export step!

"""
    EulerSymp(arguments;Niter,implicit,realdata)

Symplectic Euler solver.

Constructs an object of type `TimeSolver` to be used in `Problem(model, initial, param; solver::TimeSolver)`

Arguments can be either
0. an object of type `AbstractModel`;
1. an `Array` of size `(N,m)` where `N` is the number of collocation points, and `m` is irrelevant;
2. an integer `N` being the number of collocation points;
3. a `NamedTuple` containing a key `N`.

The keyword argument `Niter` (optional, defaut value = 10) determines the number of steps in the Neumann iteration solver of the implicit step.
The keyword argument `implicit` (optional, defaut value = 1) determines which equation is implicit (must be `1` or `2`).
The keyword argument `realdata` is optional, and determines whether pre-allocated vectors are real- or complex-valued.
By default, they are either determined by the model or the type of the array in case `0.` and `1.`, complex-valued otherwise.

"""
struct EulerSymp <: TimeSolver

    U1 :: Array
    U2   :: Array
    Niter:: Int
    implicit:: Int

    function EulerSymp( U :: Array; Niter = 10, implicit = 1, realdata=nothing )
        U1 = copy(U[:,1])
        U2 = copy(U[:,1])
        if realdata==true
            U1 = real.(U1);U2 = real.(U2)
        end
        if realdata==false
            U1 = complex.(U1);U2 = complex.(U2)
        end
        if implicit!=1 && implicit!=2
            @warn "the keyword `implicit` must be 1 or 2. solve! will not work."
        end
        new( U1, U2, Niter, implicit)
    end

    function EulerSymp( model :: AbstractModel; Niter = 10, implicit = 1, realdata=nothing )
        U=model.mapto(Init(x->0*x,x->0*x))
        EulerSymp(U; Niter = Niter, realdata=realdata, implicit=implicit)
    end
    function EulerSymp( N::Int; Niter=10, implicit = 1, realdata=false )
        U = zeros(Float64, (N,2))
        EulerSymp(U; Niter = Niter, realdata=realdata,implicit=implicit)
    end
    function EulerSymp( param::NamedTuple, datasize=2::Int; Niter = 10, implicit = 1, realdata=false )
        EulerSymp(param.N; Niter = Niter, realdata=realdata,implicit=implicit)
    end
end

function step!(solver :: EulerSymp,
                model :: AbstractModel,
                U  ,
                dt )

    solver.U1 .= copy(U[:,1])
    solver.U2 .= copy(U[:,2])

    if solver.implicit == 1
        for i=1:solver.Niter
            model.f1!( solver.U1, solver.U2 )
            solver.U1 .= U[:,1] + dt * solver.U1
        end
        U[:,1] .= solver.U1
        model.f2!( solver.U1, solver.U2 )
        U[:,2] .+= dt * solver.U2
    elseif solver.implicit == 2
        for i=1:solver.Niter
            model.f2!( solver.U1, solver.U2 )
            solver.U2 .= U[:,2] + dt * solver.U2
        end
        U[:,2] .= solver.U2
        model.f1!( solver.U1, solver.U2 )
        U[:,1] .+= dt * solver.U1
    else
        error("when defining `EulerSymp`, the keyword `implicit` must be either 1 or 2.")
    end
end


text:
n = parse(Int64, readline())

for i in 1:n
    S = readline()

    len = length(S)

    odd = ""
    even = ""

    for j in 1:2:len
        odd = odd*S[j]
    end
    for j in 2:2:len
        even = even*S[j]
    end
    println(odd," ",even)
    # println(even)
end

max_stars_repo_path: Hackerrank/30 Days of Code/Julia/day 06.jl

max_stars_repo_name: Next-Gen-UI/Code-Dynamics

max_issues_repo_name: Next-Gen-UI/Code-Dynamics

max_forks_repo_name: Next-Gen-UI/Code-Dynamics


text:
# run this script to update the pages for Documenter.jl

# to generate doc pages
using Literate
using Dates

function update_date(content)
    content = replace(content, "DATEOFTODAY" => Date(now()))
    return content
end

output_folder = joinpath(@__DIR__, "..", "docs", "src", "generated")
output_test_folder = joinpath(@__DIR__, "generated")
foreach(x -> rm(x, recursive = true), readdir(output_test_folder, join = true))
foreach(x -> rm(x, recursive = true), readdir(output_folder, join = true))

# clearing files in the output folders

function run_literate(names, folder)
    for t in names
        fp = joinpath(@__DIR__, folder, string(t, ".jl"))
        of = joinpath(output_folder, folder)
        otf = joinpath(output_test_folder, folder)
        mkpath(of)
        mkpath(otf)
        println("writing files to: \n", of, "\n", otf)
        Literate.markdown(
            fp,
            of,
            documenter = true,
            credit = false,
            preprocess = update_date,
        )
        Literate.script(
            fp,
            otf,
            documenter = true,
            credit = false,
            preprocess = update_date,
        )
        #		Literate.notebook(fp, output_folder, documenter=true,credit=false)
    end
end

run_literate(["doc"], "media")

run_literate(["doc"], "ageom")

run_literate(["doc"], "srcwav")

run_literate(["doc"], "data")

# run_literate(["doc","forw", "test_born"], "Poisson")

#run_literate(["gradient_accuracy","born_map"], "fwi")

# run_literate(["doc", "pizza", "born_tutorial"], "fwi")
run_literate(["doc", "create_snaps", "reuse_expt"], "fdtd")


text:
export SolutionInv

# Solution wrappers
@compat struct SolutionInvWrapper{QR<:ApproxFunBase.QROperator,F<:Fun,T} <: Operator{T}
  op::QR
  u::F
end
SolutionInvWrapper(op::ApproxFunBase.QROperator,f::Fun) = SolutionInvWrapper{typeof(op),typeof(f),eltype(op)}(op,f)
ApproxFun.qr!(op::SolutionInvWrapper) = op
ApproxFun.qr(op::SolutionInvWrapper) = op
#SolutionInvWrapper(op::Operator) = SolutionInvWrapper(op,uniform(domainspace(op)))
ApproxFun.@wrapper SolutionInvWrapper

#ApproxFun.\{S,T,DD,dim}(A::SolutionInvWrapper,b::Fun{MatrixSpace{S,T,DD,dim}};kwds...) = \(L.op,b;kwds...) # avoid method ambiguity
(\)(L::SolutionInvWrapper,b::ApproxFun.Fun;kwds...) = \(L.op,b;kwds...)

function uniform(S::Space)
  u = Fun(one,S)
  @compat rmul!(u.coefficients,1/sum(u))
  u
end
uniform(D::Domain) = uniform(Space(D))

function SolutionInv(L::Operator,u::Fun=uniform(domainspace(L)))
  @assert domain(L) == rangedomain(L)
  if isa(domainspace(L),ApproxFun.TensorSpace) #TODO: put this into ApproxFun
    di = DefiniteIntegral(domainspace(L).spaces[1])⊗DefiniteIntegral(domainspace(L).spaces[2])
    for i = 3:length(domainspace(L).spaces)
      di = di⊗DefiniteIntegral(domainspace(L).spaces[2])
    end
  else
    di = DefiniteIntegral(domainspace(L))
  end

  SolutionInvWrapper(ApproxFun.qr(I-L + cache(u/sum(u) * di)),u)
end
SolutionInv(M::AbstractMarkovMap,u::Fun=uniform(Space(domain(M)))) = SolutionInv(Transfer(M,space(u)),u)

Transfer(K::SolutionInvWrapper) = K.op.R_cache.op.ops[2].op # eek??

# function DefiniteIntegral(sp::ProductSpace) # likely a hack
#   di = DefiniteIntegral(sp.spaces[1])
#   for i = 2:length(sp.spaces)
#     di = di⊗DefiniteIntegral(sp.spaces[i])
#   end
#   di
# end


text:
using Sctp
using Test
import Sockets: IPv4, @ip_str

@testset "Sctp.jl" begin
    # Write your tests here.
    sctp_init()
    sock = SctpSocket()

    bind(sock, ip"0.0.0.0", 5050)

    # TODO listen()
end

@testset "Bind socket" begin
end


text:
# https://genieframework.github.io/Genie.jl/dev/tutorials/7--Using_JSON_Payloads.html
using Genie, Genie.Router, Genie.Requests, Genie.Renderer.Json

route("/hello") do
  "Hello World"
end

route("/jsonpayload", method = POST) do
  @show jsonpayload()
  @show rawpayload()

  #json("Hello $(jsonpayload()["name"])")
  json("Hello $(jsonpayload()["name"])")
end

up()

max_stars_repo_name: RoyCCWang/WebAppRKHS

max_issues_repo_name: RoyCCWang/WebAppRKHS

max_forks_repo_name: RoyCCWang/WebAppRKHS


text:
"""
    fit_quasi!(gcm::GaussianCopulaVCModel, solver=Ipopt.IpoptSolver)

Fit an `GaussianCopulaVCModel` object by MLE using a nonlinear programming solver. Start point
should be provided in `gcm.β`, `gcm.θ`, `gcm.τ` this is for Normal base.

# Arguments
- `gcm`: A `GaussianCopulaVCModel` model object.
- `solver`: Specified solver to use. By default we use IPOPT with 100 quas-newton iterations with convergence tolerance 10^-6.
    (default `solver = Ipopt.IpoptSolver(print_level=3, max_iter = 100, tol = 10^-6, limited_memory_max_history = 20, warm_start_init_point="yes", hessian_approximation = "limited-memory")`)
"""
function fit!(
        gcm::GaussianCopulaVCModel,
        solver=Ipopt.IpoptSolver(print_level = 3, tol = 10^-6, max_iter = 100,
        limited_memory_max_history = 20, warm_start_init_point="yes", hessian_approximation = "limited-memory")
    )
    initialize_model!(gcm)
    npar = gcm.p + gcm.m + 1
    optm = MathProgBase.NonlinearModel(solver)
    # set lower bounds and upper bounds of parameters
    lb   = fill(-Inf, npar)
    ub   = fill( Inf, npar)
    offset = gcm.p + 1
    for k in 1:gcm.m + 1
        lb[offset] = 0
        offset += 1
    end
    MathProgBase.loadproblem!(optm, npar, 0, lb, ub, Float64[], Float64[], :Max, gcm)
    # starting point
    par0 = zeros(npar)
    modelpar_to_optimpar!(par0, gcm)
    MathProgBase.setwarmstart!(optm, par0)
    # optimize
    MathProgBase.optimize!(optm)
    optstat = MathProgBase.status(optm)
    optstat == :Optimal || @warn("Optimization unsuccesful; got $optstat")
    # update parameters and refresh gradient
    optimpar_to_modelpar!(gcm, MathProgBase.getsolution(optm))
    loglikelihood!(gcm, true, false)
    # gcm
end

"""
    modelpar_to_optimpar!(par, gcm)

Translate model parameters in `gcm` to optimization variables in `par` for Normal base.
"""
function modelpar_to_optimpar!(
        par :: Vector,
        gcm :: GaussianCopulaVCModel
    )
    # β
    copyto!(par, gcm.β)
    # L
    offset = gcm.p + 1
    @inbounds for k in 1:gcm.m
        par[offset] = gcm.θ[k]
        offset += 1
    end
    par[offset] = gcm.τ[1]
    par
end

"""
    optimpar_to_modelpar_quasi!(gcm, par)

Translate optimization variables in `par` to the model parameters in `gcm`.
"""
function optimpar_to_modelpar!(
        gcm :: GaussianCopulaVCModel,
        par :: Vector
    )
    # β
    copyto!(gcm.β, 1, par, 1, gcm.p)
    offset = gcm.p + 1
    @inbounds for k in 1:gcm.m
        gcm.θ[k] = par[offset]
        offset   += 1
    end
    gcm.τ[1] = par[offset]
    gcm
end

function MathProgBase.initialize(
    gcm::GaussianCopulaVCModel,
    requested_features::Vector{Symbol})
    for feat in requested_features
        if !(feat in [:Grad, :Hess])
            error("Unsupported feature $feat")
        end
    end
end

MathProgBase.features_available(gcm::GaussianCopulaVCModel) = [:Grad, :Hess]

function MathProgBase.eval_f(
        gcm :: GaussianCopulaVCModel,
        par :: Vector
    )
    optimpar_to_modelpar!(gcm, par)
    loglikelihood!(gcm, false, false) # don't need gradient here
end

function MathProgBase.eval_grad_f(
    gcm  :: GaussianCopulaVCModel,
    grad :: Vector,
    par  :: Vector
    )
    optimpar_to_modelpar!(gcm, par)
    obj = loglikelihood!(gcm, true, false)
    # gradient wrt β
    copyto!(grad, gcm.∇β)
    # gradient wrt variance comps
    offset = gcm.p + 1
    @inbounds for k in 1:gcm.m
        grad[offset] = gcm.∇θ[k]
        offset += 1
    end
    grad[offset] = gcm.∇τ[1]
obj
end

MathProgBase.eval_g(gcm::GaussianCopulaVCModel, g, par) = nothing
MathProgBase.jac_structure(gcm::GaussianCopulaVCModel) = Int[], Int[]
MathProgBase.eval_jac_g(gcm::GaussianCopulaVCModel, J, par) = nothing

function MathProgBase.hesslag_structure(gcm::GaussianCopulaVCModel)
    m◺ = ◺(gcm.m)
    # we work on the upper triangular part of the Hessian
    arr1 = Vector{Int}(undef, ◺(gcm.p) + m◺ + 1)
    arr2 = Vector{Int}(undef, ◺(gcm.p) + m◺ + 1)
    # Hββ block
    idx  = 1
    for j in 1:gcm.p
        for i in j:gcm.p
            arr1[idx] = i
            arr2[idx] = j
            idx += 1
        end
    end
    # variance components
    for j in 1:gcm.m
        for i in 1:j
            arr1[idx] = gcm.p + i
            arr2[idx] = gcm.p + j
            idx += 1
        end
    end
    arr1[idx] = gcm.p + gcm.m + 1
    arr2[idx] = gcm.p + gcm.m + 1
    return (arr1, arr2)
end

function MathProgBase.eval_hesslag(
        gcm   :: GaussianCopulaVCModel,
        H   :: Vector{T},
        par :: Vector{T},
        σ   :: T,
        μ   :: Vector{T}
    )where {T <: BlasReal}
    optimpar_to_modelpar!(gcm, par)
    loglikelihood!(gcm, true, true)
    # Hβ block
    idx = 1
    @inbounds for j in 1:gcm.p, i in 1:j
        H[idx] = gcm.Hβ[i, j]
        idx   += 1
    end
    # Haa block
    @inbounds for j in 1:gcm.m, i in 1:j
        H[idx] = gcm.Hθ[i, j]
        idx   += 1
    end
    H[idx] = gcm.Hτ[1, 1]
    # lmul!(σ, H)
    H .*= σ
end


text:
using RvSpectMLBase
using Test

@testset "Utility functions that don't use package's types." begin
    include("physics.jl")
    include("cs_algs.jl")
    include("files.jl")
end

@testset "Utility functions for RvSpecMLBase's types" begin
    include("chunk_of_spectrum.jl")
    include("theoretical_instrument.jl")
    #include("find_line.jl")                                # TODO: Add tests
    include("chunk_list.jl")
    include("spectra.jl")
    include("spectra_timeseries.jl")
    #include("spectra_timeseries_common_wavelengths.jl")    # TODO: Add tests
end

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
# module LoRDeckCodes

"""
```julia
struct CardCodeAndCount
    code::String
    count::Int
end
```
"""
struct CardCodeAndCount
    code::String
    count::Int

    function CardCodeAndCount(code::String, count::Int)
        new(code, count)
    end

    function CardCodeAndCount(set::Int, faction::String, number::Int, count::Int)
        new(string(lpad(set, 2, '0'), faction, lpad(number, 3, '0')), count)
    end

    function CardCodeAndCount(card_string::AbstractString)
        count, code = split(card_string, ':')
        new(code, parse(Int, count))
    end
end

function Base.getproperty(card::CardCodeAndCount, prop::Symbol)
    if prop === :set
        parse(Int, card.code[1:2])
    elseif prop === :faction
        card.code[3:4]
    elseif prop === :number
        parse(Int, card.code[5:7])
    else
        getfield(card, prop)
    end
end

"""
```julia
struct Deck
    cards::Vector{CardCodeAndCount}
    version::UInt8
end
```
"""
struct Deck
    cards::Vector{CardCodeAndCount}
    version::UInt8
end

function Base.:(==)(a::Deck, b::Deck)
    by = x -> (x.code, x.count)
    length(a.cards) == length(b.cards) && sort(a.cards, by = by) == sort(b.cards, by = by) && a.version == b.version
end

struct ArgumentException <: Exception
    msg::String
end

# module LoRDeckCodes


text:
using Test
using JLD
mutable struct problem
    data::Array{Float64}
    Sigma::Array{Float64}
end
###############################################################################
### We have hard-coded the test data here, in order to avoid Julia dependency issues
###############################################################################
include(normpath(joinpath(@__FILE__,"..",".."))*"Algorithm/branchAndBound.jl")
include(normpath(joinpath(@__FILE__,"..",".."))*"Algorithm/utilities.jl")
include(normpath(joinpath(@__FILE__,"..",".."))*"Algorithm/multiComponents.jl")

################################################################################
######## Test #1: Pitprops data set
################################################################################
@testset "Berk and Bertsimas SPCA Tables 6-7" begin

pitprops=[[1,0.954,0.364,0.342,-0.129,0.313,0.496,0.424,0.592,0.545,0.084,-0.019,0.134];
       [0.954,1,0.297,0.284,-0.118,0.291,0.503,0.419,0.648,0.569,0.076,-0.036,0.144];
       [0.364,0.297,1,0.882,-0.148,0.153,-0.029,-0.054,0.125,-0.081,0.162,0.22,0.126];
       [0.342,0.284,0.882,1,0.22,0.381,0.174,-0.059,0.137,-0.014,0.097,0.169,0.015];
       [-0.129,-0.118,-0.148,0.22,1,0.364,0.296,0.004,-0.039,0.037,-0.091,-0.145,-0.208];
       [0.313,0.291,0.153,0.381,0.364,1,0.813,0.09,0.211,0.274,-0.036,0.024,-0.329];
       [0.496,0.503,-0.029,0.174,0.296,0.813,1,0.372,0.465,0.679,-0.113,-0.232,-0.424];
       [0.424,0.419,-0.054,-0.059,0.004,0.09,0.372,1,0.482,0.557,0.061,-0.357,-0.202];
       [0.592,0.648,0.125,0.137,-0.039,0.211,0.465,0.482,1,0.526,0.085,-0.127,-0.076];
       [0.545,0.569,-0.081,-0.014,0.037,0.274,0.679,0.557,0.526,1,-0.319,-0.368,-0.291];
       [0.084,0.076,0.162,0.097,-0.091,-0.036,-0.113,0.061,0.085,-0.319,1,0.029,0.007];
       [-0.019,-0.036,0.22,0.169,-0.145,0.024,-0.232,-0.357,-0.127,-0.368,0.029,1,0.184];
       [0.134,0.144,0.126,0.015,-0.208,-0.329,-0.424,-0.202,-0.076,-0.291,0.007,0.184,1];]
pitprops=reshape(pitprops, (13,13));
B=sqrt(pitprops);
theProb=problem(B, pitprops)
println("Running pitprops data set k=5")
obj, xVal, timetoBound, timetoConverge, timeOut, explored, toPrint=branchAndBound(theProb, 5)
@show explored, timetoConverge
println("Testing pitprops data set k=5")
@show @test(abs(obj-3.40615)<1e-4)
println("Running pitprops data set k=10")
obj, xVal, timetoBound, timetoConverge, timeOut, explored, toPrint=branchAndBound(theProb, 10)
@show explored, timetoConverge
println("Testing pitprops data set k=10")
@show @test(abs(obj-4.17264)<1e-4)

################################################################################
#### Clear data
################################################################################
theProb=nothing;
pitprops=nothing;
B=nothing;
################################################################################
######## Test #2: Wine data set.
###############################################################################
wine=[[ 0.659062    , 0.0856113  , 0.0471152    ,  -0.841093   ,  3.13988   , 0.146887   ,  0.192033   , -0.0157543   , 0.0635175   ,  1.02828    ,-0.0133134   , 0.0416978   ,   164.567		];
[   0.0856113  ,  1.24802   ,  0.050277    ,    1.07633   ,  -0.87078  , -0.234338  ,  -0.45863   ,   0.0407334  , -0.141147   ,   0.644838  , -0.143326   , -0.292447   ,    -67.5489   ];
[   0.0471152  ,  0.050277  ,  0.0752646   ,    0.406208  ,   1.12294  ,  0.0221456 ,   0.0315347 ,   0.00635847 ,  0.00151558 ,   0.164654  , -0.00468215 ,  0.000761836,     19.3197   ];
[  -0.841093   ,  1.07633   ,  0.406208    ,   11.1527    ,  -3.97476  , -0.671149  ,  -1.17208   ,   0.150422   , -0.377176   ,   0.145024  , -0.209118   , -0.656234   ,   -463.355    ];
[   3.13988    , -0.87078   ,  1.12294     ,   -3.97476   , 203.989    ,  1.91647   ,   2.79309   ,  -0.455563   ,  1.93283    ,   6.62052   ,  0.180851   ,  0.669308   ,   1769.16     ];
[   0.146887   , -0.234338  ,  0.0221456   ,   -0.671149  ,   1.91647  ,  0.39169   ,   0.54047   ,  -0.0350451  ,  0.219373   ,  -0.0799975 ,  0.0620389  ,  0.311021   ,     98.1711   ];
[   0.192033   , -0.45863   ,  0.0315347   ,   -1.17208   ,   2.79309  ,  0.54047   ,   0.997719  ,  -0.066867   ,  0.373148   ,  -0.399169  ,  0.124082   ,  0.558262   ,    155.447    ];
[  -0.0157543  ,  0.0407334 ,  0.00635847  ,    0.150422  ,  -0.455563 , -0.0350451 ,  -0.066867  ,   0.0154886  , -0.0260599  ,   0.0401205 , -0.00747118 , -0.0444692  ,    -12.2036   ];
[   0.0635175  , -0.141147  ,  0.00151558  ,   -0.377176  ,   1.93283  ,  0.219373  ,   0.373148  ,  -0.0260599  ,  0.327595   ,  -0.0335039 ,  0.0386646  ,  0.210933   ,     59.5543   ];
[   1.02828    ,  0.644838  ,  0.164654    ,    0.145024  ,   6.62052  , -0.0799975 ,  -0.399169  ,   0.0401205  , -0.0335039  ,   5.37445   , -0.276506   , -0.705813   ,    230.767    ];
[  -0.0133134  , -0.143326  , -0.00468215  ,   -0.209118  ,   0.180851 ,  0.0620389 ,   0.124082  ,  -0.00747118 ,  0.0386646  ,  -0.276506  ,  0.052245   ,  0.0917662  ,     17.0002   ];
[   0.0416978  , -0.292447  ,  0.000761836 ,   -0.656234  ,   0.669308 ,  0.311021  ,   0.558262  ,  -0.0444692  ,  0.210933   ,  -0.705813  ,  0.0917662  ,  0.504086   ,     69.9275   ];
[ 164.567      ,-67.5489    , 19.3197      , -463.355     ,1769.16     , 98.1711    , 155.447     , -12.2036     , 59.5543     , 230.767     , 17.0002     , 69.9275     ,  99166.7      ];]
wine=reshape(wine, (13,13))
B=sqrt(wine);
theProb=problem(B, wine)
println("Running wine data set k=5")
obj, xVal, timetoBound, timetoConverge, timeOut, explored, toPrint=branchAndBound(theProb, 5)

println("Testing wine data set k=5")
@show @test(abs(obj-99201.31)<1.0)

println("Running wine data set k=10")
obj, xVal, timetoBound, timetoConverge, timeOut, explored, toPrint=branchAndBound(theProb, 10)

println("Testing wine data set k=10")
@show @test(abs(obj-99201.78)<1.0)
################################################################################
#### Clear data
################################################################################
theProb=nothing;
wine=nothing;
B=nothing;
################################################################################
######## Test #3: Norm wine data set.
###############################################################################

normwine=[[1.0       ,  0.0943969 ,  0.211545   , -0.310235  ,  0.270798  ,  0.289101 ,   0.236815 , -0.155929 ,  0.136698    , 0.546364   ,-0.0717472 ,  0.0723432 ,   0.64372   ];
[0.0943969 ,  1.0       ,  0.164045   ,  0.2885    , -0.0545751 , -0.335167 ,  -0.411007 ,  0.292977 , -0.220746    , 0.248985   ,-0.561296  , -0.36871   ,  -0.192011  ];
[0.211545  ,  0.164045  ,  1.0        ,  0.443367  ,  0.286587  ,  0.12898  ,   0.115077 ,  0.18623  ,  0.00965194  , 0.258887   ,-0.0746669 ,  0.00391123,   0.223626  ];
[-0.310235 ,   0.2885   ,   0.443367  ,   1.0      ,  -0.0833331,  -0.321113,   -0.35137 ,   0.361922,  -0.197327   ,  0.018732  , -0.273955 ,  -0.276769 ,   -0.440597 ];
[0.270798  , -0.0545751 ,  0.286587   , -0.0833331 ,  1.0       ,  0.214401 ,   0.195784 , -0.256294 ,  0.236441    , 0.19995    , 0.0553982 ,  0.0660039 ,   0.393351  ];
[0.289101  , -0.335167  ,  0.12898    , -0.321113  ,  0.214401  ,  1.0      ,   0.864564 , -0.449935 ,  0.612413    ,-0.0551364  , 0.433681  ,  0.699949  ,   0.498115  ];
[0.236815  , -0.411007  ,  0.115077   , -0.35137   ,  0.195784  ,  0.864564 ,   1.0      , -0.5379   ,  0.652692    ,-0.172379   , 0.543479  ,  0.787194  ,   0.494193  ];
[-0.155929 ,   0.292977 ,   0.18623   ,   0.361922 ,  -0.256294 ,  -0.449935,   -0.5379  ,   1.0     ,  -0.365845   ,  0.139057  , -0.26264  ,  -0.50327  ,   -0.311385 ];
[0.136698  , -0.220746  ,  0.00965194 , -0.197327  ,  0.236441  ,  0.612413 ,   0.652692 , -0.365845 ,  1.0         ,-0.0252499  , 0.295544  ,  0.519067  ,   0.330417  ];
[0.546364  ,  0.248985  ,  0.258887   ,  0.018732  ,  0.19995   , -0.0551364,  -0.172379 ,  0.139057 , -0.0252499   , 1.0        ,-0.521813  , -0.428815  ,   0.3161    ];
[-0.0717472,  -0.561296 ,  -0.0746669 ,  -0.273955 ,   0.0553982,   0.433681,    0.543479,  -0.26264 ,   0.295544   , -0.521813  ,  1.0      ,   0.565468 ,    0.236183 ];
[0.0723432 , -0.36871   ,  0.00391123 , -0.276769  ,  0.0660039 ,  0.699949 ,   0.787194 , -0.50327  ,  0.519067    ,-0.428815   , 0.565468  ,  1.0       ,   0.312761  ];
[0.64372   , -0.192011  ,  0.223626   , -0.440597  ,  0.393351  ,  0.498115 ,   0.494193 , -0.311385 ,  0.330417    , 0.3161     , 0.236183  ,  0.312761  ,   1.0	];]
normwine=reshape(normwine, (13,13))
B=sqrt(normwine);
theProb=problem(B, normwine)
println("Running normwine data set k=5")
obj, xVal, timetoBound, timetoConverge, timeOut, explored, toPrint=branchAndBound(theProb, 5)

println("Testing wine data set k=5")
@show @test(abs(obj-3.43978)<1e-4)

println("Running normwine data set k=10")
obj, xVal, timetoBound, timetoConverge, timeOut, explored, toPrint=branchAndBound(theProb, 10)

println("Testing wine data set k=10")
@show @test(abs(obj-4.59429)<1e-4)
################################################################################
#### Clear data
################################################################################
theProb=nothing;
normwine=nothing;
B=nothing;

################################################################################
######## Test #4: miniBooNE data set.
###############################################################################
miniboone=load("data/miniBoone.jld",  "miniBooNE")
B=sqrt(miniboone)
theProb=problem(B, miniboone)
println("Running miniBooNE data set k=5")
obj, xVal, timetoBound, timetoConverge, timeOut, explored, toPrint=branchAndBound(theProb, 5)

println("Testing miniBooNE data set k=5")
@show @test(abs(obj-1.96827e9)<1e7) # source: Berk+B paper
                                                    # Large numbers so numerical precision won't be great, hence tolerance is loose.
println("Running miniBooNE data set k=10")
obj, xVal, timetoBound, timetoConverge, timeOut, explored, toPrint=branchAndBound(theProb, 10)

println("Testing miniBooNE data set k=10")
@show @test(abs(obj-1.96827e9)<1e7)
################################################################################
#### Clear data
################################################################################
spca_miniboone=nothing;
miniboone=nothing;
B=nothing;

################################################################################
######## Test #4: normMiniBooNE data set.
###############################################################################
normminiboone=load("data/miniBoone.jld",  "normMiniBooNE")
B=sqrt(normminiboone)
theProb=problem(B, normminiboone)
println("Running normMiniBooNE data set k=5")
obj, xVal, timetoBound, timetoConverge, timeOut, explored, toPrint=branchAndBound(theProb, 5)

println("Testing normMiniBooNE data set k=5")
@show @test(abs(obj-5.0000)<1e-2)

println("Running normMiniBooNE data set k=10")
obj, xVal, timetoBound, timetoConverge, timeOut, explored, toPrint=branchAndBound(theProb, 10)

println("Testing normMiniBooNE data set k=10")
@show @test(abs(obj-9.99999)<1e-2)
################################################################################
#### Clear data
################################################################################
theProb=nothing;
normminiboone=nothing;
B=nothing;
################################################################################
######## Test #6: Communities data set.
###############################################################################

communities=load("data/communities.jld",  "communities");
B=sqrt(communities);
theProb=problem(B, communities)

println("Running communities data set k=5");
obj, xVal, timetoBound, timetoConverge, timeOut, explored, toPrint=branchAndBound(theProb, 5)
@show explored, timetoConverge
println("Testing communities data set k=5")
@show @test(abs(obj-0.2771)<1e-2)

println("Running communities data set k=10")
obj, xVal, timetoBound, timetoConverge, timeOut, explored, toPrint=branchAndBound(theProb, 10)
@show explored, timetoConverge
println("Testing communities data set k=10")
@show @test(abs(obj-0.4452)<1e-2)

################################################################################
#### Clear data
################################################################################
theProb=nothing;
communities=nothing;
B=nothing;
################################################################################
######## Test #7: Norm Communities data set.
###############################################################################
normcommunities=load("data/communities.jld", "normCommunities")
B=sqrt(normcommunities);
theProb=problem(B, normcommunities)
println("Running norm communities data set k=5")
obj, xVal, timetoBound, timetoConverge, timeOut, explored, toPrint=branchAndBound(theProb, 5)
@show explored, timetoConverge
println("Testing norm communities data set k=5")

@show @test(abs(obj-4.86051)<1e-2)

println("Running norm communities data set k=10")
obj, xVal, timetoBound, timetoConverge, timeOut, explored, toPrint=branchAndBound(theProb, 10)
@show explored, timetoConverge
println("Testing norm communities data set k=10")
@show @test(abs(obj-8.8236)<1e-2)
################################################################################
#### Clear data
################################################################################
spca_normcommunities=nothing;
normcommunities=nothing;
B=nothing;
end;
################################################################################
##### End test
################################################################################

max_stars_repo_name: ryancorywright/Optimal-SPCA

max_issues_repo_name: ryancorywright/Optimal-SPCA

max_forks_repo_name: ryancorywright/Optimal-SPCA


text:
module cmap_helpers

using Pandas, PyCall
PyCmap = pyimport("pycmap")
#cmap = PyCmap.API(token="your-own-API-key")
cmap = PyCmap.API()

"""
    get(t::String,v::String)

Retrieve variable v from CMAP table t. Return along with
meta-data, position and time, in a Dict.
"""
function get(t::String,v::String)
    df=Pandas.DataFrame(cmap.get_dataset(t))
    me=Pandas.DataFrame(cmap.get_metadata(t,v))
    x=Dict("Variable" => v,"Unit" => deepcopy(values(me["Unit"])[1]),
    "Long_Name" => deepcopy(values(me["Long_Name"])[1]),
    "Data_Source" => deepcopy(values(me["Data_Source"])[1]),
    "lon" => deepcopy(values(df[:lon])),
    "lat" => deepcopy(values(df[:lat])),
    "time" => deepcopy(values(df[:time])),
    "val" => deepcopy(values(df[Symbol(v)])))
    return x
end

"""
    tables(ListName::String)

List of CMAP tables for subset `ListName`.
"""
    function tables(ListName::String)
        list0=[];
        if ListName=="G3"
            list0=["tblKM1906_Gradients3","tblKM1906_Gradients3_uway_optics","tblKM1906_Gradients3_uwayCTD"]
            #"tblKM1906_Gradients3_uw_tsg"
        elseif ListName=="G2"
            list0=["tblMGL1704_Gradients2_CTD","tblMGL1704_Gradients2_uway_optics"]
        elseif ListName=="G1"
            list0=["tblKOK1606_Gradients1_CTD","tblKOK1606_Gradients1_uway_optics"]
        elseif ListName=="MoreG1"
            list0=["tblKOK1606_Gradients1_Nutrients","tblKOK1606_Gradients1_Dissolved_Gasses",
            "tblKOK1606_Gradients1_TargetedMetabolites","tblKOK1606_Gradients1_Diazotroph"]
        elseif ListListNameId=="MoreG2"
            list0=["tblMGL1704_Gradients2_Nutrients","tblMGL1704_Gradients2_Diazotroph",
            "tblMGL1704_Gradients2_TargetedMetabolites","tblMGL1704_Gradients2_Trace_Metals"]
        end
        return list0
    end

end

module cbiomes_helpers

using MeshArrays, MITgcmTools
pth="../samples/gradients/"
!isdir("$pth") ? mkdir("$pth") : nothing

"""
    myinterp(pth::String,v::String,lon,lat)

Retrieve variable v from model output (`v*.".bin" binary file) and
interpolate using coefficients provided in `c*.".mat" file.
Return result as a Dict along with lon & lat.
"""
function myinterp(pth::String,v::String,lon,lat)

    gridpath="$pth"*"GRID_LLC90/"
    gitpath="https://github.com/gaelforget/GRID_LLC90"
    !isdir(gridpath) ? run(`git clone $gitpath $gridpath`) : nothing
    γ=GridSpec("LatLonCap",gridpath)
    Γ=GridLoad(γ)

    μ=γ.read(gridpath*"$v"*".bin",MeshArray(γ,Float32))
    msk=1.0 .+ 0.0 * mask(view(Γ["hFacC"],:,1),NaN,0.0)
    (f,i,j,w)=InterpolationFactors(Γ,vec(lon),vec(lat))
    μ=Interpolate(μ.*msk,f,i,j,w)

    return Dict("val" => μ, "lon" => lon, "lat" => lat, "Data_Source" => "ECCOv4r2 (Gael Forget)")
end

end

max_stars_repo_name: gaelforget/MarineEcosystemNotebooks

max_forks_repo_name: gaelforget/MarineEcosystemNotebooks


text:
# TODO: Function to sync optimiser state. Will become easier once Optimisers.jl becomes
#       the standard

"""
broadcast_parameters(model; root_rank::Integer = 0, blocking_communication::Bool = true)
broadcast_parameters(ps::Params; root_rank::Integer = 0, blocking_communication::Bool = true)

Sync the parameters of the model across all processes.
"""
broadcast_parameters(model; kwargs...) = broadcast_parameters(params(model); kwargs...)

function broadcast_parameters(ps::Params; root_rank::Integer=0, blocking_communication::Bool=true)
    @assert 0 <= root_rank <= total_workers() - 1 "Valid `root_rank` Range: [0, $(total_workers() - 1)]"
    if blocking_communication
        requests = Vector{Request}(undef, length(ps))
        for (i, p) in enumerate(ps)
            _, request = Ibcast!(p, root_rank, MPI.COMM_WORLD)
            requests[i] = request
        end
        Waitall!(requests)
    else
        Bcast!.(ps, root_rank, MPI.COMM_WORLD)
    end
    return
end


text:
# module tempmod # for redefining structure easily while developing
	mutable struct Tex
		eachline::Vector{<:AbstractString}
		istarget::Vector{Bool} # Default is all true. See outer constructor.
		targetrange::Vector{UnitRange{Int64}}

	end
	function Tex(eachline::Vector{<:AbstractString})
	targetrange = map(x -> 1:length(x), eachline);
	istarget =  .!isempty.(targetrange); # 1:0 is empty. Which means "" is not the target.

	Tex(eachline,istarget,targetrange);
	# require more complex function to also auto generate correct targetrange?
	end

# end

"""
`refresh!(TEX::Tex)` refresh `TEX.istarget` according to `TEX.targetrange`.
This could be useless.
"""
function refresh!(TEX0::Tex)
	TEX0.istarget .= .!isempty.(TEX0.targetrange); # 1:0 is empty. Which means "" is not the target.
end

"""
```
texpath1 = joinpath(pwd(),"main1.tex")
texpath0 = joinpath(pwd(),"main0.tex")
pushback(texpath0, texpath1);
```
Push the contents of `tex1` back into the original `tex0`, replacing the original contents while preserving comments and begin-ends sections in `tex0`.
These contents are those not belong to
- section (e.g. `\\section{}`, `\\figure{}`, etc.)
- comments (e.g. `% hello I'm comment`)
"""
function pushback(texpath0, texpath1)
	# should be equivalent to `tex0 = open(x->readlines(x), texpath0)`;
	tex0 = open(texpath0) do file
		readlines(file);
	end
	tex1 = open(texpath1) do file
		readlines(file);
	end

	TEX0 = Tex(tex0);
	TEX1 = Tex(tex1);

	# criteria:
	# before \begin{document}, all lines should be neglected
	# all sections (e.g. \begin{xxx}...\end{xxx}) should be ignored except
	# all lines begin with `\` despite `\cite`

	# Expression for the criteria that defines our target sentences.
	expr1 = r"^([A-Za-z\[\]\(\)\$\"`]|[1-9]|,|;|\\cite|\s+\w)[^%]*(?=%?)";
	# We target lines who start with:
	# 	[A-Za-z\[\]\(\)\$] denotes
	# 		- normal characters (A-Z and a-z)
	# 		- parenthesis "(", ")", "[", "]"
	# 		- variable or inline equation, e.g. "$F = ma$ the Newton's 2nd law blablabla"
	# and
	# 		- inline citation may be just at the first of a line ("\cite{blablabla2021}")
	# 		- anything starts with one or more white space and words (\s+\w)
	# and anything that is not a comment "[^%]*"
	# that may ended with one or no comment "(?=%?)".

	# https://regex101.com/r/x0HEer/1
	pushback_rm!(TEX0, expr1);
	pushback_rm!(TEX1, expr1);

	if sum(TEX0.istarget) != sum(TEX1.istarget)
		outputpath0_a = pushback_genfile1(texpath0, TEX0);
		outputpath0_b = pushback_genfile1(texpath1, TEX1);
		openit(dirname(outputpath0_b));
		fname1_b = basename(outputpath0_b);
		fname1_a = basename(outputpath0_a);
		fname0_a = basename(texpath0);
		fname0_b = basename(texpath1);

		error("""Inconsistent in the number of lines. Please edit on `$fname1_b` to make the lines consistent with `$fname1_a` (e.g. use vscode compare), where adding dummy new lines might be necessary (e.g. valid: "(CHECKPOINT: nothing here)"; invalid: "\\n"). Then, run the same function but with modified later file (e.g. `pushback("$(fname0_a)","$(fname1_b)")`). After that, it is suggested to compare $fname1_b with the original $fname0_b""");
	end

	newstr = TEX1.eachline[TEX1.istarget];
	for (i, TF) in enumerate(TEX0.istarget)
		if TF
			str0v = split(TEX0.eachline[i],""); # split into vector of chars
			rnga = TEX0.targetrange[i]; # range of non-comment
			# str0a = join(str0v[rnga]); # old content without comment
			str0v[rnga] .= ""; # remove old content but reserve comment.
			str0c = join(str0v); # comment only
			str1a = popfirst!(newstr); # new content
			TEX0.eachline[i] = str1a*str0c; # finish replacing old content by new one despite old comment.
		end
	end
	outputpath0_b = pushback_genfile2(texpath0, TEX0);
	fnm0 = basename(texpath0);
	fnm1 = basename(texpath1);
	fnm2 = basename(outputpath0_b);
	println("Content successfully combined. Since a few specific types of contents won't be successfully pushed back (e.g. changes in title and subtitle), please check the results by: ");
	println("1. Comaring $fnm2 with $fnm1 (old) to check if contents were correctly inserted. (There shouldn't be any differences except comments and begin-end sections)")
	println("2. Comparing $fnm2 with $fnm0 for english editing viewing.");
	openit(dirname(outputpath0_b));
end

function pushback_genfile1(texpath0,TEX0::Tex)
		fname0 = basename(texpath0);
		strvec0 = getext(fname0);
		strvec0[1] = strvec0[1]*"_contentonly";
		fname0_a = join(strvec0)
		dir_0 = dirname(texpath0);
		outputpath0_b = joinpath(dir_0,fname0_a);
		open(outputpath0_b,"w") do io
			for (i, TF) in enumerate(TEX0.istarget)
				if TF
					range_i = TEX0.targetrange[i];
					println(io, TEX0.eachline[i][range_i]);
				end
			end
		end
		return outputpath0_b
end

function pushback_genfile2(texpath0,TEX0::Tex)
		fname0 = basename(texpath0);
		strvec0 = getext(fname0);
		strvec0[1] = strvec0[1]*"_combined";
		fname0_a = join(strvec0)
		dir_0 = dirname(texpath0);
		outputpath0_b = joinpath(dir_0,fname0_a);
		open(outputpath0_b,"w") do io
			for i in eachindex(TEX0.eachline)
				println(io, TEX0.eachline[i]);
			end
		end
		return outputpath0_b
end

function pushback_rm!(TEX0::Tex,expr1::Regex)
	for (i,mt) in enumerate(match.(expr1, TEX0.eachline))
		if !isnothing(mt)
			TEX0.targetrange[i] = mt.offset:length(mt.match);
			TEX0.istarget[i] = true; # this could be unnecessary
		else # nothing matched
			TEX0.targetrange[i] = 1:0;
			TEX0.istarget[i] = false;
		end
	end
	removesection!(TEX0, "figure");
	removesection!(TEX0, "table");
	removesection!(TEX0, "equation");
end

"""
removecomment!(tex)
- Remove all comments (text after the `%`) in the `.tex` file.
- `tex` should be an 1d array of strings, with each element the texts of a line.
- Please check again the result (e.g. check with git compare).

# Example:
```
tex = open(inputfpath_tex) do file
	readlines(file);
end
removecomment!(tex)
```
"""
function removecomment!(tex::Vector{String})
	ind2rm = Int64[];
	for i = 1:length(tex)
		onelncmt = match( r"^%",tex[i]);# oneline comment
		if isnothing(onelncmt)
			m = match(r"%", tex[i])
			if ~isnothing(m)
				newstr = tex[i][1:m.offset-1]
				tex[i] = newstr;
			end
		else # it begins with "%" therefore the whole line should be removed
			push!(ind2rm, i);
		end
	end
	deleteat!(tex, ind2rm);
	return nothing
end

"""
removesection!(tex::Vector{<:AbstractString}, sectionname)
- Remove section(s) in the `.tex` file.
- In which, every section should begin with the first line being `\\begin{sectionname[*]}` and the last line `\\end{sectionname[*]}` (where the asterisk denotes the width of the section (e.g. table, figure) is two column width).

For example, `removesection!(tex, "figure")` remove both

```
\\begin{figure}
<Figure here>
\\end{figure}
```

and

```
\\begin{figure*}
<Figure here>
\\end{figure*}
```
# WARNING:
The `\\begin...` or `\\end...` must be at the first of a line. For example, match missed when ` \\begin...`.
"""
function removesection!(tex::Vector{<:AbstractString}, sectionname)
	id2rm = removesection_secrange(tex,sectionname);
    deleteat!(tex, id2rm);
    return nothing
end

"""
removesection!(TEX::Tex, sectionname)
Modify the `Tex` object where the field `istarget` is set to false for the corresponding section specified by `sectionname`. This function only change the indicators, not any content will be removed.
"""
function removesection!(TEX::Tex, sectionname)
	tex = TEX.eachline;
	id2rm = removesection_secrange(tex,sectionname);
	for i in id2rm
		TEX.targetrange[i] = 1:0; # don't use broadcast
		TEX.istarget[i] = false;
	end
    return nothing
end

function removesection_secrange(tex,sectionname)
    begins = findall(.~isnothing.(match.(Regex("^(\\\\begin{$(sectionname)[\\*]*})"), tex)));
    ends = findall(.~isnothing.(match.(Regex("^(\\\\end{$(sectionname)[\\*]*})"), tex)));
    id2rm = Int64[];
    for (id0, id1) in (begins .=> ends)
        id2rm = vcat(id2rm, collect(id0:id1));
    end
	return id2rm
end


text:
export HomogeneousModel

"""
    HomogeneousModel{I,O} <: FixedTimeModel{I,O}

A dynamic model with a fixed time step and same transition model at every time point..
The constructor is called with the initial sfunc and transition sfuncs.
The constructor is also called with an optional dt (defaults to 1).
"""
struct HomogeneousModel{I,J,O} <: TimelessFixedTimeModel{I,J,O}
    initial :: SFunc{I,O}
    transition :: SFunc{J,O}
    dt :: Number
end
function HomogeneousModel(init::SFunc{I,O}, trans::SFunc{J,O}, dt = 1) where {I,J,O}
    HomogeneousModel{I,J,O}(init, trans, dt)
end

get_initial(m::HomogeneousModel) = m.initial

get_transition(m::HomogeneousModel) = m.transition

get_dt(m::HomogeneousModel) = m.dt


text:
module TestExpressionsIndexColumnNames

import Test: @testset, @test
import Volcanito: index_column_names, ColumnName

test_cases = (
    (
        input = (ColumnName(:a, false), ColumnName(:b, false)),
        output = Dict(ColumnName(:a, false) => 1, ColumnName(:b, false) => 2),
    ),
    (
        input = (
            ColumnName(:a, false),
            ColumnName(:b, false),
            ColumnName(:c, false),
        ),
        output = Dict(
            ColumnName(:a, false) => 1,
            ColumnName(:b, false) => 2,
            ColumnName(:c, false) => 3,
        ),
    ),
)

@testset "index_column_names"  begin
    for case in test_cases
        @test index_column_names(case.input) == case.output
    end
end

end


text:
### Singletons ###
abstract type ASetConstant <: AConstantSymbol end

# abstract type AΩ <: ASetConstant end
# Ω = Term(Symbol("Ω"), AΩ)

# ∅ = Term(Symbol("∅"), ASetConstant)
@new_constant ASetConstant AΩ Ω Sets
@new_constant ASetConstant A∅ ∅ Sets
@new_constant ASetConstant AStatements Statements Sets
@new_constant ASetConstant ASets Sets Sets

### Functions ###
abstract type ASetFunction <: AFunctionSymbol end

@new_function(ASetFunction, ASetUnion, SetUnion, (Sets, Sets), Sets, 2, ∪, false)
@new_function(ASetFunction, ASetIntersection, SetIntersection, (Sets, Sets), Sets, 2, ∩, false)
@new_function(ASetFunction, ASetDifference, SetDifference, (Sets, Sets), Sets, 2, -, false)

@new_function(ASetFunction, AConditionSet, ConditionSet, (Variable, Statements), Sets, 2)
ConditionSet(x::Term, condition::Term) = Term(:ConditionSet, ConditionSet, (x, condition))  # {x | condition}
Base.show(io::IO, self::Term{<:AConditionSet}) = print(io, "$(Meta.parse("$(data(self))($(args(self)[1]), $(args(self)[2]))"))")


text:
__precompile__()
"""
Reachability analysis of high-dimensional affine systems with nondeterministic inputs.

## Introduction

This module provides an interface to set-based reachability algorithms for affine
systems with nondeterministic inputs. It implements several reachability
algorithms using decomposition into low-dimensional subsystems. The code is optimized
for high-dimensional linear systems.

ReachSets.jl relies in the following modules:

|Module name|Features|
|-----------|--------|
|`LazySets.jl`|Set-based representations and operations in the lazy framework. Support vector and support function calculus.|
|`MathematicalSystems.jl`|Representation and transformation of affine systems with non-deterministic inputs.|

The following dependencies are optional and provide additional functionality:

|Optional module name | Features|
|------|------|
|`Expokit.jl` | Fast computation of matrix exponential of very large and very sparse matrices.|
|`Polyhedra.jl`, `CDDLib.jl` | Julia Polyhedra library and the CDD backend.|

The visualization routines are deferred to the corresponding `LazySets.jl`, such as
`plot_polygon`, implemented in `Polygon.jl`.

This module is accompanied with a set of unit tests, available in the `test` folder.

## Available algorithms

The dictionary ``available_algorithms`` specifies the registered algorithms. Each
algorithm is registered with its name, and a dictionary value which any number
of entries. Most common entries are:

|Key|Description|
|-----|----|
|``'func'``|the corresponding function's handle|
|``'is_explicit'`` |if true, the output contains explicitly computed polygons; otherwise, the lazy-set framework is returned|

The algorithms implemented in `ReachSets/` are summarized in the following
table.

|Function name|Description|
|-------------|------|
|``lazy``     | Lazy computation (all blocks)|
|``explicit`` | Explicit computation|

For additional information on other algorithms see the documentation in their
corresponding source files.

[1] Decomposing Reachability Computations for Affine Systems.
Sergiy Bogomolov, M. F., Goran Frehse, Andreas Podelski, Christian Schilling,
Frédéric Viry. 2017, to be submitted.

## Projections

The final goal of any reachability computation is a sequence of points (vertices)
describing the evolution of one or more variables in time, as a function of
other variables, or a linear combination of them. This information is required
for plotting the output in two or three dimensions.

This module provides the function ``project_reach`` for the purpose of projecting
a high-dimensional explicit or lazy set computation into a low number of variables
efficiently. Additionally, the time variable can be passed (use `0`), and an arbitrary
linear combination of state variables can be handled as well.

## Helper functions and macros

This module provides commonly used

|Macro or function name| Description|
|-----|-----|
|``@filename_to_png``|Convert the current script name's suffix to ".png"|
|``@block_id v``| Return the block number associated to a given variable.|

AUTHORS:

- Marcelo Forets
- Christian Schilling
- Frederic Viry
"""
module ReachSets

using LazySets, MathematicalSystems, Expokit, ..Utils, ProgressMeter

import LazySets.Approximations:symmetric_interval_hull,
                               symmetric_interval_hull_parallel,
                               decompose,
                               decompose_parallel,
                               overapproximate,
                               box_approximation
import Reachability.tocc

# ========================================
# Discretize and compute bloating factors
# ========================================
include("discretize.jl")

export discretize

# =======================================
# Mapping from input to output variables
# =======================================
include("inout_map_reach.jl")

export inout_map_reach

# ========================
# Reachability Algorithms
# ========================

# dictionary of registered algorithms
available_algorithms = Dict{String, Dict{String, Any}}()

# "explicit" backends
include("reach_blocks.jl")
push!(available_algorithms, "explicit_blocks"=>Dict("func"=>reach_blocks!,
                                                    "is_explicit"=>true))

# "explicit" backends
include("reach_blocks_parallel.jl")
push!(available_algorithms, "explicit_blocks_parallel"=>Dict("func"=>reach_blocks_parallel!,
                                                        "is_explicit"=>true))

include("reach_blocks_wrapping_effect.jl")
push!(available_algorithms,
      "wrap"=>Dict("func"=>reach_blocks_wrapping_effect!,
                   "is_explicit"=>true))

export available_algorithms

# =========================
# External reach interface
# =========================
include("reach.jl")

export reach

# ==============================================
# Projection of the reach set in two dimensions
# ==============================================
include("project_reach.jl")

export project_reach

end # module ReachSets


text:
module ScikitLearn_

#> for all Supervised models:
import MLJBase
import MLJBase: @mlj_model, metadata_model,
    _process_model_def, _model_constructor, _model_cleaner
using ScientificTypes
using Tables

#> for all classifiers:
using CategoricalArrays
using Tables

# NOTE: legacy code for SVM models does not use the @sk_model macro.

#> import package:
import ..ScikitLearn: @sk_import
import ..ScikitLearn
@sk_import svm: SVC
@sk_import svm: NuSVC
@sk_import svm: LinearSVC
@sk_import svm: SVR
@sk_import svm: NuSVR
@sk_import svm: LinearSVR

include("svm.jl")

const Option{T} = Union{Nothing, T}

const SKLM = ((ScikitLearn.Skcore).pyimport("sklearn.linear_model"))
const SKGP = ((ScikitLearn.Skcore).pyimport("sklearn.gaussian_process"))
const SKEN = ((ScikitLearn.Skcore).pyimport("sklearn.ensemble"))
const SKDU = ((ScikitLearn.Skcore).pyimport("sklearn.dummy"))
const SKNB = ((ScikitLearn.Skcore).pyimport("sklearn.naive_bayes"))
const SKNE = ((ScikitLearn.Skcore).pyimport("sklearn.neighbors"))
const SKDA = ((ScikitLearn.Skcore).pyimport("sklearn.discriminant_analysis"))
# const SKNN = ((ScikitLearn.Skcore).pyimport("sklearn.neural_network"))

"""
_skmodel_fit_reg

Called as part of [`@sk_reg`](@ref), returns the expression corresponing to the `fit` method
for the ScikitLearn regression model.
"""
function _skmodel_fit_reg(modelname, params)
        quote
                function MLJBase.fit(model::$modelname, verbosity::Int, X, y)
                        # set X and y into a format that can be processed by sklearn
                        Xmatrix   = MLJBase.matrix(X)
                        yplain    = y
                        targnames = nothing
                        # check if it's a multi-target regression case, in that case keep
                        # track of the names of the target columns so that the prediction
                        # can be named accordingly
                        if Tables.istable(y)
                           yplain    = MLJBase.matrix(y)
                           targnames = Tables.schema(y).names
                        end
                        # Call the parent constructor from Sklearn.jl named $modelname_
                        skmodel = $(Symbol(modelname, "_"))($((Expr(:kw, p, :(model.$p)) for p in params)...))
                        fitres  = ScikitLearn.fit!(skmodel, Xmatrix, yplain)
                        # TODO: we may want to use the report later on
                        report  = NamedTuple()
                        # the first nothing is so that we can use the same predict for reg and clf
                        return ((fitres, nothing, targnames), nothing, report)
                end
        end
end

"""
_skmodel_fit_clf

Called as part of [`@sk_clf`](@ref), returns the expression corresponing to the `fit` method
for the ScikitLearn classifier model.
"""
function _skmodel_fit_clf(modelname, params)
        quote
                function MLJBase.fit(model::$modelname, verbosity::Int, X, y)
                        Xmatrix = MLJBase.matrix(X)
                        yplain  = MLJBase.int(y)
                        skmodel = $(Symbol(modelname, "_"))($((Expr(:kw, p, :(model.$p)) for p in params)...))
                        fitres  = ScikitLearn.fit!(skmodel, Xmatrix, yplain)
                        # TODO: we may want to use the report later on
                        report  = NamedTuple()
                        # pass y[1] for decoding in predict method, first nothing is targnames
                        return ((fitres, y[1], nothing), nothing, report)
                end
        end
end

"""
_skmodel_predict

Called as part of [`@sk_model`](@ref), returns the expression corresponing to the `predict` method
for the ScikitLearn model (for a deterministic model)
"""
function _skmodel_predict(modelname)
        quote
                function MLJBase.predict(model::$modelname, (fitresult, y1, targnames), Xnew)
                        Xmatrix = MLJBase.matrix(Xnew)
                        preds   = ScikitLearn.predict(fitresult, Xmatrix)
                        if isa(preds, Matrix)
                                # only regressors are possibly multitarget;
                                # build a table with the appropriate column names
                                return MLJBase.table(preds, names=targnames)
                        end
                        if y1 !== nothing
                                # if it's a classifier)
                                return preds |> MLJBase.decoder(y1)
                        end
                        return preds
                end
        end
end

"""
_skmodel_predict_prob

Same as `_skmodel_predict` but with probabilities. Note that only classifiers are probabilistic
in sklearn so that we always decode.
"""
function _skmodel_predict_prob(modelname)
        quote
                # there are no multi-task classifiers in sklearn
                function MLJBase.predict(model::$modelname, (fitresult, y1, _), Xnew)
                        Xmatrix = MLJBase.matrix(Xnew)
                        # this is an array of size n x c with rows that sum to 1
                        preds   = ScikitLearn.predict_proba(fitresult, Xmatrix)
                        classes = MLJBase.classes(y1)
                        return [MLJBase.UnivariateFinite(classes, preds[i, :]) for i in 1:size(Xmatrix,1)]
                end
        end
end

# --------------------------------------------------------
# functions to help have a different macro for regressors
# and classifiers, these functions are just there to avoid
# duplication of code
# --------------------------------------------------------
function _sk_constructor(ex)
        # similar to @mlj_model
    ex, modelname, params, defaults, constraints = _process_model_def(ex)
        # keyword constructor
    const_ex = _model_constructor(modelname, params, defaults)
        # associate the constructor with the definition of the struct
    push!(ex.args[3].args, const_ex)
        # cleaner
    clean_ex = _model_cleaner(modelname, defaults, constraints)
        # return
        return modelname, params, clean_ex, ex
end

function _sk_finalize(m, clean_ex, fit_ex, ex)
        # call a different predict based on whether probabilistic or deteterministic
        if ex.args[2].args[2] == :(MLJBase.Probabilistic)
                predict_ex = _skmodel_predict_prob(m)
        else
                predict_ex = _skmodel_predict(m)
        end
    esc(
                quote
                        # Base.@__doc__ $ex
                export $m
                $ex
                $fit_ex
                $clean_ex
                $predict_ex
                MLJBase.load_path(::Type{<:$m})           = "MLJModels.ScikitLearn_.$(MLJBase.name($m))"
                MLJBase.package_name(::Type{<:$m})    = "ScikitLearn"
                MLJBase.package_uuid(::Type{<:$m})    = "3646fa90-6ef7-5e7e-9f22-8aca16db6324"
                MLJBase.is_pure_julia(::Type{<:$m})   = false
                MLJBase.package_url(::Type{<:$m})     = "https://github.com/cstjean/ScikitLearn.jl"
                MLJBase.package_license(::Type{<:$m}) = "BSD"
            end
        )
end

"""
macro sk_reg(ex)

Helper macro for defining interfaces of ScikitLearn regression models. Struct fields require a type
annotation and a default value as in the example below. Constraints for parameters (fields) are
introduced as for field3 below. The constraint must refer to the parameter as `arg`. If the used
parameter does not meet a constraint the default value is used.

@sk_reg mutable struct SomeRegression <: MLJBase.Deterministic
    field1::Int = 1
    field2::Any = nothing
    field3::Float64 = 0.5::(0 < arg < 0.8)
end

MLJBase.fit and MLJBase.predict methods are also produced. See also [`@sk_clf`](@ref)
"""
macro sk_reg(ex)
        modelname, params, clean_ex, ex = _sk_constructor(ex)
        fit_ex = _skmodel_fit_reg(modelname, params)
        _sk_finalize(modelname, clean_ex, fit_ex, ex)
end

"""
macro sk_clf(ex)

Same as [`@sk_reg`](@ref) but for classifiers.
"""
macro sk_clf(ex)
        modelname, params, clean_ex, ex = _sk_constructor(ex)
        fit_ex = _skmodel_fit_clf(modelname, params)
        _sk_finalize(modelname, clean_ex, fit_ex, ex)
end

include("linear-regressors.jl")
include("linear-classifiers.jl")
include("gaussian-process.jl")
include("ensemble.jl")
include("discriminant-analysis.jl")
include("misc.jl")

end # module