text:
# Examining results from the search for the eighth planet
# without the Earth-like core-mass-fraction constraint, with
# random initial conditions.

using JLD2; using PyPlot; using DelimitedFiles

elements_7planet = readdlm("../../../data/elements_noprior_students.txt",',')

fnames = readdlm("../../../data/v13/file_list.txt");

size(fnames)

# Minimum chi-square with 7 planets:
chi7 = 1396.205; chimin = 1396.205
chi_all = Float64[]; p8_all = Float64[]; m8_all = Float64[]; t08_all = Float64[]; ecc8_all = Float64[]
m8_nonzero = Float64[]; chi_nonzero = Float64[]
mb_all = Float64[]; mc_all = Float64[]; md_all = Float64[]; me_all = Float64[]; mf_all = Float64[]; mg_all = Float64[]; mh_all = Float64[]
count = 0
elements_all = zeros(9,7,1651)
iall = 1
for fname in fnames
    @load string("../../../data/v13/",fname) chi_grid elements_grid
    for j=1:100
        if chi_grid[j] != 0.0 && chi_grid[j] < chi7 && elements_grid[9,1,j] > 0
            push!(m8_nonzero,elements_grid[9,1,j])
            push!(chi_nonzero,chi_grid[j])
            #println(fname," ",j," ",elements_grid[2:9,1:5,j]," ",chi_grid[j])
        end
        if chi_grid[j] != 0.0
            global count +=1
        end
        ecc8 = sqrt(elements_grid[9,4,j]^2+elements_grid[9,5,j]^2)
        if chi_grid[j] != 0.0 && chi_grid[j] < chi7 && elements_grid[9,1,j] > 0.0 # && ecc8 < 0.01
          #println(fname," ",j," ",chi_grid[j]," ",chi7-chi_grid[j]," ",elements_grid[2:9,1:5,j])
          if chi_grid[j] < chimin
                global chimin = chi_grid[j]
                println(fname," j: ",j," chi_min: ",chimin," ",elements_grid[2:9,1:5,j])
          end
          push!(chi_all,chi_grid[j])
          push!(p8_all,elements_grid[9,2,j])
          push!(t08_all,elements_grid[9,3,j])
          #push!(ecc8_all,sqrt(elements_grid[9,4,j]^2+elements_grid[9,5,j]^2))
          push!(ecc8_all,ecc8)
          push!(m8_all,elements_grid[9,1,j])
          push!(mb_all,elements_grid[2,1,j]);  push!(mc_all,elements_grid[3,1,j]);
          push!(md_all,elements_grid[4,1,j]);  push!(me_all,elements_grid[5,1,j])
          push!(mf_all,elements_grid[6,1,j]);  push!(mg_all,elements_grid[7,1,j])
          push!(mh_all,elements_grid[8,1,j])
          elements_all[:,:,iall] .= elements_grid[:,:,j]
          global iall += 1
        end
    end
end
nall = iall - 1

# Okay, so the results show that the minimum chi-square
# is 1362.7, which is about ~33 less than a fit with
# seven planets with CMF-earth constraint.  So, not a
# huge improvement, and this is close to a 2:1 resonance,
# but not near a Laplace resonance.
# which breaks the pattern seen with the inner 7.  Also, the
# eccentricity is pretty large (~0.08), and the mass is fairly
# small (1.24x10^{-5} ~ 0.37 M_earth.  BIC expected is ~5 log(447)
# ~ 30.5, while AIC is ~10, so either way, not a great improvement.
# It's interesting that the best-fit values are similar to the
# CMF-earth-like constraint case.  But, it would be more
# convincing if the delta-chi-square were significant.  Either
# way we can use this for a photometric search for transits,
# as well as an estimate of the systematic errors.

pygui(false)

if ~@isdefined(CGS)
    include("../../../src/CGS.jl")
    using Main.CGS
end

fac_mass = 0.09*MSUN/MEARTH

fig,axes = subplots(1,3,sharey="row",constrained_layout="True",figsize=(6,4))
ax = axes[1]
ax.plot(p8_all,chi7 .-chi_all,".",alpha=0.5)
ax.plot(p8_all[ecc8_all .< 0.01],chi7 .-chi_all[ecc8_all .< 0.01],"o")

ax.set_xlabel(L"$P_i$ [days]")
ax.plot([20.0,45.0],[30.5,30.5])
ax.axis([20,45,0,40])
ax.text([30],[29],L"$\Delta$ BIC = 0")
ax.set_ylabel(L"$\Delta 2 \log\mathcal{L}$")

ax = axes[2]
ax.plot(ecc8_all,chi7 .-chi_all,".",alpha=0.5)
ax.plot(ecc8_all[ecc8_all .< 0.01],chi7 .-chi_all[ecc8_all .< 0.01],"o")

ax.set_xlabel(L"$e_i$")
ax.plot([0.0,0.35],[30.5,30.5])
ax.axis([0,0.35,0,40])
ax.text([0.1],[29],L"$\Delta$ BIC = 0")
#ax.set_ylabel(L"$\Delta 2\log\mathcal{L}$")

ax = axes[3]
ax.plot(m8_all*fac_mass,chi7 .-chi_all,".",alpha=0.5)
ax.plot(m8_all[ecc8_all .< 0.01]*fac_mass,chi7 .-chi_all[ecc8_all .< 0.01],"o")

ax.plot([0.0,1.0],[30.5,30.5])
ax.text([0.5],[29],L"$\Delta$ BIC = 0")
ax.axis([0.0,1.0,0,40])
ax.set_xlabel(L"$M_i/M_\oplus$")
#ax.set_ylabel(L"$\Delta 2\log\mathcal{L}$")
subplots_adjust(wspace=0)
#tight_layout()
savefig("../Planet_i_properties.pdf",bbox_inches="tight")

#savefig("Planet_i_mass.pdf",bbox_inches="tight")

using Statistics
(median(chi_all)-minimum(chi_all))/std(chi_all)

using Statistics; using Printf
m7 =[1.3774,1.3109,0.3885, 0.6933, 1.0410, 1.3237, 0.3264]
m7_opt = elements_7planet[2:8,1]*fac_mass
s7 =[0.0593,0.0450, 0.0074, 0.0128, 0.0155, 0.0171, 0.0186]

println("b ",@sprintf("%6.3f",median(mb_all)*fac_mass-m7[1]),"+-",@sprintf("%5.3f",s7[1]),"+-",@sprintf("%5.3f",std(mb_all)*fac_mass))
println("c ",@sprintf("%6.3f",median(mc_all)*fac_mass-m7[2]),"+-",@sprintf("%5.3f",s7[2]),"+-",@sprintf("%5.3f",std(mc_all)*fac_mass))
println("d ",@sprintf("%6.3f",median(md_all)*fac_mass-m7[3]),"+-",@sprintf("%5.3f",s7[3]),"+-",@sprintf("%5.3f",std(md_all)*fac_mass))
println("e ",@sprintf("%6.3f",median(me_all)*fac_mass-m7[4]),"+-",@sprintf("%5.3f",s7[4]),"+-",@sprintf("%5.3f",std(me_all)*fac_mass))
println("f ",@sprintf("%6.3f",median(mf_all)*fac_mass-m7[5]),"+-",@sprintf("%5.3f",s7[5]),"+-",@sprintf("%5.3f",std(mf_all)*fac_mass))
println("g ",@sprintf("%6.3f",median(mg_all)*fac_mass-m7[6]),"+-",@sprintf("%5.3f",s7[6]),"+-",@sprintf("%5.3f",std(mg_all)*fac_mass))
println("h ",@sprintf("%6.3f",median(mh_all)*fac_mass-m7[7]),"+-",@sprintf("%5.3f",s7[7]),"+-",@sprintf("%5.3f",std(mh_all)*fac_mass))

fig,axes = subplots(2,4,sharey="row",figsize=(10,10))
ax = axes[1]
ax.plot(mb_all.*fac_mass,chi7 .-chi_all,".",label="b")
ax.plot(mb_all[ecc8_all .< 0.01].*fac_mass,chi7 .-chi_all[ecc8_all .< 0.01],".",label="low ecc")
ax.fill_between([m7[1]-s7[1],m7[1]+s7[1]],[0,0],[40,40],alpha=0.2,color="C1")
#ax.fill_between([m7_opt[1]-s7[1],m7_opt[1]+s7[1]],[0,0],[40,40],alpha=0.2,color="C2")
ax.set_xlabel(L"$M_b/M_\oplus$")
ax.set_ylabel(L"$\Delta 2\log \mathcal{L}$")
ax.legend()

ax = axes[2]
ax.plot(mc_all.*fac_mass,chi7 .-chi_all,".",label="c")
ax.plot(mc_all[ecc8_all .< 0.01].*fac_mass,chi7 .-chi_all[ecc8_all .< 0.01],".",label="low ecc")
ax.fill_between([m7[2]-s7[2],m7[2]+s7[2]],[0,0],[40,40],alpha=0.2,color="C1")
#ax.fill_between([m7_opt[2]-s7[2],m7_opt[2]+s7[2]],[0,0],[40,40],alpha=0.2,color="C2")

ax.legend()
ax.set_xlabel(L"$M_c/M_\oplus$")
ax.set_ylabel(L"$\Delta 2\log \mathcal{L}$")

ax = axes[3]
ax.plot(md_all.*fac_mass,chi7 .-chi_all,".",label="d")
ax.plot(md_all[ecc8_all .< 0.01].*fac_mass,chi7 .-chi_all[ecc8_all .< 0.01],".",label="low ecc")
ax.fill_between([m7[3]-s7[3],m7[3]+s7[3]],[0,0],[40,40],alpha=0.2,color="C1")
#ax.fill_between([m7_opt[3]-s7[3],m7_opt[3]+s7[3]],[0,0],[40,40],alpha=0.2,color="C2")

ax.legend()
ax.set_xlabel(L"$M_d/M_\oplus$")
#ax.set_ylabel(L"$\Delta 2\log \mathcal{L}$")

ax = axes[4]
ax.plot(me_all.*fac_mass,chi7 .-chi_all,".",label="e")
ax.plot(me_all[ecc8_all .< 0.01].*fac_mass,chi7 .-chi_all[ecc8_all .< 0.01],".",label="low ecc")
ax.fill_between([m7[4]-s7[4],m7[4]+s7[4]],[0,0],[40,40],alpha=0.2,color="C1")
#ax.fill_between([m7_opt[4]-s7[4],m7_opt[4]+s7[4]],[0,0],[40,40],alpha=0.2,color="C2")

ax.legend()
ax.set_xlabel(L"$M_e/M_\oplus$")
#ax.set_ylabel(L"$\Delta 2\log \mathcal{L}$")

ax = axes[5]
ax.plot(mf_all.*fac_mass,chi7 .-chi_all,".",label="f")
ax.plot(mf_all[ecc8_all .< 0.01].*fac_mass,chi7 .-chi_all[ecc8_all .< 0.01],".",label="low ecc")
ax.fill_between([m7[5]-s7[5],m7[5]+s7[5]],[0,0],[40,40],alpha=0.2,color="C1")
#ax.fill_between([m7_opt[5]-s7[5],m7_opt[5]+s7[5]],[0,0],[40,40],alpha=0.2,color="C2")

ax.legend()
ax.set_xlabel(L"$M_f/M_\oplus$")
#ax.set_ylabel(L"$\Delta 2\log \mathcal{L}$")

ax = axes[6]
ax.plot(mg_all.*fac_mass,chi7 .-chi_all,".",label="g")
ax.plot(mg_all[ecc8_all .< 0.01].*fac_mass,chi7 .-chi_all[ecc8_all .< 0.01],".",label="low ecc")
ax.fill_between([m7[6]-s7[6],m7[6]+s7[6]],[0,0],[40,40],alpha=0.2,color="C1")
#ax.fill_between([m7_opt[6]-s7[6],m7_opt[6]+s7[6]],[0,0],[40,40],alpha=0.2,color="C2")

ax.legend()
ax.set_xlabel(L"$M_g/M_\oplus$")
#ax.set_ylabel(L"$\Delta 2\log \mathcal{L}$")

ax = axes[7]
ax.plot(mh_all.*fac_mass,chi7 .-chi_all,".",label="h")
ax.plot(mh_all[ecc8_all .< 0.01].*fac_mass,chi7 .-chi_all[ecc8_all .< 0.01],".",label="low ecc")
ax.fill_between([m7[7]-s7[7],m7[7]+s7[7]],[0,0],[40,40],alpha=0.2,color="C1")
#ax.fill_between([m7_opt[7]-s7[7],m7_opt[7]+s7[7]],[0,0],[40,40],alpha=0.2,color="C2")

ax.legend()
ax.set_xlabel(L"$M_h/M_\oplus$")
#ax.set_ylabel(L"$\Delta 2\log \mathcal{L}$")
ax = axes[8]; ax.axis("off"); tight_layout()
subplots_adjust(wspace=0)
savefig("T1_8planet_mass_comarison.png",bbox_inches="tight")

m7 .- m7_opt

plot(mc_all.*fac_mass,chi_all,"."); elements_7planet[3,1]*fac_mass

plot(mg_all.*fac_mass,chi_all,"."); elements_7planet[7,1]*fac_mass

# Wow, this case has a *very* small mass for planet i,
# yet improves on the chi-square from our 7-planet fit.
# That's not good...  but it's a very high eccentricity.
@load "../../../data/T1_CMFearth_run094.jld2" chi_grid elements_grid

# Note to future self:  this cell below (between the multi-line commens)
# took about 5 hours to run.
# The results are saved in T1_8planet_noprior_efree.jld2, but I should
# have saved specific variables. To get back the interesting quantities
# in the future without having to rerun these simulations (which took about 5 hours),
# just carry out:

@load "../../../data/T1_8planet_noprior_efree.jld2" ecosfree esinfree

# The rest of the necessary quantities were computed in the cells above,
# which are quick to run.

# Now, let's compute the free eccentricities of these planets
# and see if we can find a relation to the initial eccentricities:

#=
include("/Users/ericagol/Observing/Spitzer/DDT2019/Campaign04/v09/compute_elements.jl")
include("/Users/ericagol/Software/NbodyGradient/src/ttv.jl")

using Statistics

function normalize(phi)
while phi < 0
  phi += 2pi
end
while phi > 2pi
  phi -= 2pi
end
return phi
end

function offset(phi,nt)
for i=1:nt-1
  while phi[i+1] > (phi[i]+3.0)
    phi[i+1] -= pi
  end
  while phi[i+1] < (phi[i]-3.0)
    phi[i+1] += pi
  end
end
return phi
end

@load "../../../data/T1_CMFearth_run001.jld2"
const NDIM  = 3
const third = 1.0/3.0
const alpha0 = 0.0

ecosfree = zeros(n-1,nall)
esinfree = zeros(n-1,nall)
tt1 = zeros(n,1160)
rstar = 1e12
tmax = 1000.0
i1 = 1; i2 = 16668
elements  = zeros(n,7)

# I'm going to compute the free eccentricities in the
# same manner as I did for the v09/Hyak/stability_10k/ case
# /Users/ericagol/Observing/Spitzer/DDT2019/Campaign04/v09/Hyak/stability_10k/compute_free_eccentricity_10k.jl
# so that I can see if I can use the same transformation in
# the eight planet case as I did in the seven planet case:

for iall=1:nall
  # Insert into the orbital element array:
  elements .= elements_all[:,:,iall]
  fileout = string("T1_xv_8planet_",lpad(iall,5,"0"),".txt")
  rm(fileout,force=true)
  @time dq0 = ttv_elements!(n,t0,h,tmax,elements,tt1,count1,0.0,0,0,rstar;fout=fileout,iout=1)
  # Now, compute the time-dependent elements:
  t,eloft = compute_elements(fileout,elements[:,1],n)
  # Now, measure free eccentricities of the planets:
  for ip=1:n-1
    ecosfree[ip,iall] = mean(eloft[i1:i2,3,ip].*cos.(eloft[i1:i2,4,ip]))
    esinfree[ip,iall] = mean(eloft[i1:i2,3,ip].*sin.(eloft[i1:i2,4,ip]))
  end
  # Need to remove the file since these are large:
  rm(fileout,force=true)
end

# Now save the results:

#@save "T1_8planet_noprior_efree.jld2"

=#

ecosfree[:,1]

esinfree[:,1]

#@load "/Users/ericagol/Observing/Spitzer/DDT2019/Campaign04/v09/Hyak/stability_10k/efree_osc_coeff.jld2"
@load "../../../data/efree_osc_coeff.jld2"

# Convert free eccentricities to elements, and compare:
ecos_osc = zeros(7,nall); esin_osc = zeros(7,nall)
for iall = 1:nall, ip=1:7
  ecos_osc[ip,iall] = c_ecos_f2o[ip,1] + c_ecos_f2o[ip,2]*ecosfree[ip,iall] + c_ecos_f2o[ip,3]*esinfree[ip,iall]
  esin_osc[ip,iall] = c_esin_f2o[ip,1] + c_esin_f2o[ip,2]*ecosfree[ip,iall] + c_esin_f2o[ip,3]*esinfree[ip,iall]
  #println(iall," ",ip," ",elements_all[ip+1,4:5,iall]," ",ecos_osc[ip,iall]," ",esin_osc[ip,iall]," ",
  #      elements_all[ip+1,4,iall]-ecos_osc[ip,iall]," ",elements_all[ip+1,5,iall]-esin_osc[ip,iall])
end

using PyPlot

# Make plots of the residuals from the prediction for the osculating elements from
# the free elements (based on the fit to seven planets from
# v09/Hyak/stability_10k/plot_eccentricity_correlations.jl
# efree_osc_coeff.jld2 in that same directory.  Looks pretty good
# - predicted values agree to ~10^{-4} rms.  Worst agreement is for h,
# not surprisingly.

pname = ["b","c","d","e","f","g","h"]
fig,axes = subplots(1,2)
ax = axes[1]
for ip=1:7
  ax.plot(elements_all[ip+1,4,:],ecos_osc[ip,:] .-elements_all[ip+1,4,:],".",label=pname[ip])
  println(ip," ecos: ",std(ecos_osc[ip,:] .-elements_all[ip+1,4,:]))
end
ax.legend()
ax.axis([-0.015,0.015,-0.015,0.015])
ax = axes[2]
for ip=1:7
  ax.plot(elements_all[ip+1,5,:],esin_osc[ip,:] .-elements_all[ip+1,5,:],".",label=pname[ip])
  println(ip," esin: ",std(esin_osc[ip,:] .-elements_all[ip+1,5,:]))
end
ax.legend()
ax.axis([-0.015,0.015,-0.015,0.015])

# Plot chi-square of 8-planet fit vs. residuals of forecast
# osculating eccentricity vectors (based on free eccentricity
# vector).  Some of the lower values of chi-square have the
# significant outliers.

for ip=7:7
  plot(ecos_osc[ip,:] .-elements_all[ip+1,4,:],chi_all,".")
  plot(esin_osc[ip,:] .-elements_all[ip+1,5,:],chi_all,".",label=pname[ip])
end
legend()

# Plot residuals of the prediction vs. orbital period of outer planet.
# Looks most significant near ~39 days (~2:1 period ratio, or p=q=1):
for ip=7:7
  plot(ecos_osc[ip,:] .-elements_all[ip+1,4,:],elements_all[9,2,:]./elements_all[8,2,:],".")
  plot(esin_osc[ip,:] .-elements_all[ip+1,5,:],elements_all[9,2,:]./elements_all[8,2,:],".",label=pname[ip])
end
legend()

# Plot residuals of the prediction vs. eccentricity of outer planet.
# Looks most significant near small eccentricity:
for ip=7:7
  plot(ecos_osc[ip,:] .-elements_all[ip+1,4,:],elements_all[9,4,:],".")
  plot(esin_osc[ip,:] .-elements_all[ip+1,5,:],elements_all[9,5,:],".",label=pname[ip])
end
legend()


text:
# A higher level stepper, defined as a wrapper around another stepper.

#TODO: how about having an DenseStepper <: AbstractWrapper <: AbstractStepper?
immutable DenseStepper <: AbstractStepper
    solver::Solver
end

solve(ode::ExplicitODE,
      stepper::DenseStepper,
      options::Options) = Solver(ode,stepper,options)

dense(sol::Solver) = solve(sol.ode, DenseStepper(sol), sol.options)

"""

The state of the dense stepper

- s0, s1: Previous steps, used to produce interpolated output
- solver_state: The state of the associated solver
- ytmp: work array

"""
type DenseState{T,S} <: AbstractState
    s0::Step{T,S}
    s1::Step{T,S}
    last_tout::T
    first_step
    solver_state::AbstractState
    # used for storing the interpolation result
    ytmp::S
    solver_done
end

function start(s::Solver{DenseStepper})
    # extract the real solver
    solver = s.stepper.solver
    t0  = solver.ode.t0
    y0  = solver.ode.y0
    dy0 = deepcopy(y0)
    solver.ode.F!(t0,y0,dy0)
    step0 = Step(t0,deepcopy(y0),deepcopy(dy0))
    step1 = Step(t0,deepcopy(y0),deepcopy(dy0))
    solver_state = start(solver)
    ytmp = deepcopy(y0)
    return DenseState(step0, step1, t0-1, true, solver_state, ytmp, false)
end

# m3: I think it would be nice to factor out the dense-output and
# root-finding into its own function.  That way it could be used also
# independently of the dense-output iterator.  Also, it would make
# this next function more compact.

# pwl: I agree, but then the problem is that once you decouple them
# you would lose the opprotunity to detect the roots with each step.
function next(s::Solver{DenseStepper}, state::DenseState)

    # m3: I'm not 100% sure what happens here.  I would implement it like so:
    # Initialize in `start`: calculate next t1, y1 and also hold onto IC in t0,y0,
    # set state.last_t=1
    #
    # in next have this loop
    # for t in tspan[state.last_t:end]
    #     if t>t1
    #        make new t1, y1, move old t1, y1 into t0, y0
    #     end
    #     make dense output at t
    #     find events
    #     state.last_t += 1
    #     return ((t, y), state)
    # end

    # pwl: @m3 this is basically what happens here:-), although I'm
    # not using the index of tspan anywhere explicitly.

    solver = s.stepper.solver

    # these guys store the intermediate steps we make
    s0, s1 = state.s0, state.s1
    t0, t1 = s0.t, s1.t

    # assuming the last output was done at state.last_tout set the
    # t_goal to the next larger time from tspan.  Strong inequality
    # below is crucial, otherwise we would be selecting the same step
    # every time.
    tspan  = s.options.tspan
    t_goal = tspan[findfirst(t->(t>state.last_tout), tspan)]

    # Keep computing new steps (i.e. new pairs (t0,t1)) until we reach
    # t0 < t_goal <= t1, then we use interpolation to get the value at
    # t_goal.  Unless points==:all, then we break the while loop after
    # making the first step.
    while t_goal > t1

        # s1 stores the last succesfull step, the new step is stored
        # in s0

        if done(solver, state.solver_state)
            warn("The iterator was exhausted before the dense output completed.")
            # prevents calling done(..) twice
            state.solver_done = true
            # TODO: deepcopy?
            # Return whatever we got as the last step
            return ((s0.t,s0.y[:]),state)
        else
            # at this point s0 is updated with the new step, "s2" if you will
            ((s0.t,s0.y[:]), state.solver_state) = next(solver, state.solver_state)
        end

        # swap s0 and s1
        s0, s1 = s1, s0
        # and times
        t0, t1 = s0.t, s1.t

        # update the state accordingly
        state.s0, state.s1 = s0, s1

        # we haven't reached t_goal yet (t1<t_goal) but the option
        # points==:all calls for an output at every successful
        # step.
        if s.options.points == :all && t1 < t_goal-eps(t1)
            state.last_tout = t1
            return ((t1,s1.y),state)
        end

        if s.options.stopevent(t1,s1.y)
            break
        end

    end

    # at this point we have t0 < t_goal < t1 so we can apply the
    # interpolation to get a value of the solution at t_goal

    # TODO: is this necessary?  The solver should store the value of dy.
    solver.ode.F!(t0,s0.y,s0.dy)
    solver.ode.F!(t1,s1.y,s1.dy)

    if s.options.stopevent(t1,s1.y)
        function stopfun(t)
            hermite_interp!(state.ytmp,t,s0,s1)
            res = Int(s.options.stopevent(t,state.ytmp))
            return 2*res-1      # -1 if false, +1 if true
        end
        t_goal = findroot(stopfun, [s0.t,s1.t], s.options.roottol)
        # state.ytmp is already overwriten to the correct result as a
        # side-effect of calling stopfun
    else
        hermite_interp!(state.ytmp,t_goal,s0,s1)
    end

    # update the last output time
    state.last_tout = t_goal

    return ((t_goal,state.ytmp),state)

end

function done(s::Solver{DenseStepper}, state::DenseState)

    return (
            state.solver_done ||
            state.last_tout >= s.options.tspan[end] ||
            s.options.stopevent(state.s1.t,state.s1.y)
            )
end

function hermite_interp!(y,t,step0::Step,step1::Step)
    # For dense output see Hairer & Wanner p.190 using Hermite
    # interpolation. Updates y in-place.
    #
    # f_0 = f(x_0 , y_0) , f_1 = f(x_0 + h, y_1 )
    # this is O(3). TODO for higher order.

    y0,  y1  = step0.y, step1.y
    dy0, dy1 = step0.dy, step1.dy

    if t == step0.t
        copy!(y,y0)
    elseif t == step1.t
        copy!(y,y1)
    else
        dt       = step1.t-step0.t
        theta    = (t-step0.t)/dt
        for i=1:length(y0)
            y[i] = ((1-theta)*y0[i] + theta*y1[i] + theta*(theta-1) *
                    ((1-2*theta)*(y1[i]-y0[i]) + (theta-1)*dt*dy0[i] + theta*dt*dy1[i]) )
        end
    end
end

ext: jl


text:
using Revise

using Tiger
import Tiger.build_base_model

using JuMP, Gurobi
using SparseArrays

import Base.convert

indicator_name(varname) = "I__" * varname
function add_indicator(model::Model, var::VariableRef)
    ind = @variable(model, base_name=indicator_name(name(var)), binary=true)
    @constraint(model, ind => {var == 0.0})
    return ind
end

function add_indicator(model::Model, var::String)
    add_indicator(model, variable_by_name(model, var))
end

function create_primal_dual(cobra::Tiger.CobraModel; bound_vars::Vector{String}=[])
    c = cobra.c
    lb = cobra.lb
    ub = cobra.ub
    n = size(cobra.S,2)

    ieq = cobra.csense .== '='
    Aeq = cobra.S[ieq,:]
    beq = cobra.b[ieq]

    leq = cobra.csense .== '<'
    Ale = cobra.S[leq,:]
    ble = cobra.b[leq]

    geq = cobra.csense .== '>'
    Age = cobra.S[geq,:]
    bge = cobra.b[geq]

    A = [Ale; -Age]
    b = [ble; -bge]

    model = build_base_model(cobra)
    x = variable_by_name.(model, cobra.vars)

    meq = size(Aeq,1)
    mle = size(A,1)
    @variable(model, lambda[1:meq])
    @variable(model, 0 <= mu[1:mle])
    @variable(model, 0 <= muU[1:n])
    @variable(model, 0 <= muL[1:n])

    # z = muU .* ub
    @variable(model, 0 <= z[1:n])

    @constraint(model, Aeq'*lambda + A'*mu - muL + muU .== c)

    # strict duality
    @constraint(model, c[1:n]'*x[1:n] == lambda'*beq + mu'*b - muL'*lb + z'*ub)

    # bounding constraints
    nbound = length(bound_vars)
    @variable(model, Inds[1:nbound], binary=true)
    set_name.(Inds, "I__" .* bound_vars)

    for i = 1:n
        if name(x[i]) in bound_vars
            Ind = variable_by_name(model, "I__" * name(x[i]))
            @constraint(model, z[i] == muU[i])
            @constraint(model, !Ind => {muU[i] == 0.0})
            @constraint(model, !Ind => {z[i] == 0.0})
            @constraint(model, !Ind => {x[i] <= 0.0})
        else
            @constraint(model, z[i] == muU[i])
        end
    end

    @objective(model, Max, 0)

    return model, Inds
end

function find_conditional_media(cobra::Tiger.CobraModel, ko; min_wt_growth=1.0, max_ko_growth=0.0, exchanges=nothing, max_kos=0, optimize=true)
    lp = StandardLP(cobra)
    n = length(lp.vars)

    model = Model(Gurobi.Optimizer)

    # add WT condition
    v_wt = @variable(model, [1:n])
    set_name.(v_wt, lp.vars .* "_WT")
    set_lower_bound.(v_wt, lp.lb)
    set_upper_bound.(v_wt, lp.ub)
    @constraint(model, lp.A * v_wt .<= lp.b)
    @constraint(model, lp.Aeq * v_wt .== lp.beq)

    # find the WT growth
    @objective(model, Max, cobra.c'*v_wt)
    optimize!(model)
    wt_growth = objective_value(model)

    # add KO condition with primal-dual expansion
    ko_idx = findfirst(x->x==ko, lp.vars)
    ko_lb = copy(lp.lb)
    ko_ub = copy(lp.ub)
    ko_lb[ko_idx] = 0.0
    ko_ub[ko_idx] = 0.0
    v_ko = @variable(model, [1:n])
    set_name.(v_ko, lp.vars .* "_KO")
    set_lower_bound.(v_ko, ko_lb)
    set_upper_bound.(v_ko, ko_ub)
    @constraint(model, lp.A * v_ko .<= lp.b)
    @constraint(model, lp.Aeq * v_ko .== lp.beq)

    @variable(model, lambda[1:size(lp.Aeq,1)])
    @variable(model, mu[1:size(lp.A,1)] >= 0)
    @variable(model, muU[1:n] >= 0)
    @variable(model, muL[1:n] >= 0)

    @variable(model, I[1:n], binary=true)

    if isnothing(exchanges)
        ex_names = get_exchange_rxns(cobra, objective=false)
    else
        ex_names = exchanges
    end
    Ie = [x in ex_names for x in lp.vars]

    gene_names = cobra.genes
    Ig = [x in gene_names for x in lp.vars]

    fix.(I[.!(Ie .| Ig)], 1)

    @constraint(model, lp.Aeq'*lambda + lp.A'*mu - muL + muU .== -lp.sense*lp.c)

    @constraint(model, v_wt[Ie] .<= lp.ub[Ie].*I[Ie])
    @constraint(model, v_ko[Ie] .<= ko_ub[Ie].*I[Ie])
    if max_kos > 0
        @constraint(model, v_wt[Ig] .<= lp.ub[Ig].*I[Ig])
        @constraint(model, v_ko[Ig] .<= ko_ub[Ig].*I[Ig])
        @constraint(model, sum(Ig) - sum(I[Ig]) <= max_kos)
    end
    #@constraint(model, muI[Ie] .<= M*(1 .- I[Ie]))

    # strong duality
    @constraint(model, lambda'*lp.beq + mu'*lp.b - muL'*ko_lb + muU'*(ko_ub.*I) == -lp.sense*lp.c'*v_ko)

    # constraint WT to grow
    @constraint(model, lp.c' * v_wt >= min_wt_growth * wt_growth)

    # constraint KO to not grow
    @constraint(model, lp.c' * v_ko <= max_ko_growth * wt_growth)

    if max_kos > 0
        @objective(model, Max, sum(I[Ig]) - sum(I[Ie]))
    else
        @objective(model, Min, sum(I[Ie]))
    end

    if !optimize
        media = nothing
        kos = nothing
    else
        set_optimizer_attribute(model, "IntegralityFocus", 1)
        optimize!(model)
        media = exchanges[isapprox.(value.(I[Ie]), 1.0, atol=1e-4)]
        kos = gene_names[isapprox.(value.(I[Ig]), 0.0, atol=1e-4)]
    end

    return model, media, kos
end

# =================== TESTING ===================

# cobra = read_cobra("/Users/jensen/Dropbox/research/tiger.jl/Tiger/test/cobra_or.mat", "cobra")
# #cobra = read_cobra("/Users/jensen/Dropbox/research/tiger.jl/Tiger/test/cobra_small.mat", "cobra")

# cobra = extend_cobra_cnf(cobra, ub=1.0)
# model = find_conditional_media(cobra, "g6a", min_wt_growth=0.1)
# optimize!(model)

cobra = read_cobra("/Users/jensen/Dropbox/repos/COBRA_models/iSMU.mat", "iSMU")
cobra.lb[cobra.lb .> 0.0] .= 0.0  # remove NGAM
cobra = extend_cobra_cnf(cobra, ub=1000)

# we only want the exchanges with positive fluxes; these are
# the media exchanges in iSMU
ex_names = get_exchange_rxns(cobra, objective=false)
is_ex = [x in ex_names for x in cobra.vars]
ingredients = ex_names[nonzeros(cobra.S[:,is_ex]) .> 0]
is_ingredient = [x in ingredients for x in cobra.vars]

cobra.ub[is_ingredient] .= 10
cobra.ub[cobra.ub .>  1000] .=  1000
cobra.lb[cobra.lb .< -1000] .= -1000

optimize!(build_base_model(StandardLP(cobra)))

model, media, kos = find_conditional_media(cobra, "SMU_308", min_wt_growth=0.1, exchanges=ingredients)
#optimize!(model)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "wxWidgets"
version = v"3.1.5"

version_mm = "$(version.major).$(version.minor)"
version_no_sep = "$(version.major)$(version.minor)$(version.patch)"

gen_libnames(lib) = ["libwx_gtk3u_$(lib)-$(version_mm)",
                    "wxmsw$(version_no_sep)u_$(lib)_gcc_custom",
                    "libwx_osx_cocoau_$(lib)-$(version_mm)-x86_64-apple-darwin14",
                    "libwx_osx_cocoau_$(lib)-$(version_mm)-aarch64-apple-darwin20"]

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/wxWidgets/wxWidgets/archive/refs/tags/v$version.tar.gz", "e8fd5f9fbff864562aa4d9c094f898c97f5e1274c90f25beb0bfd5cb61319dea"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""

cd $WORKSPACE/srcdir/wxWidgets-*

if [[ "${target}" == *-linux-musl* ]]; then

    #derived from https://github.com/wxWidgets/wxWidgets/commit/f4eae4df2a256d0c1bab1657438c2de449cdcd67?branch=f4eae4df2a256d0c1bab1657438c2de449cdcd67&diff=split
    #this was merged on master, so probably can get rid of next release
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/musl-add-locale-functions.patch

    #help find zlib for some reason
    export CPPFLAGS="-I${includedir}"

    # Delete libexpat to prevent it from being picked up by mistake
    rm /usr/lib/libexpat.so*

elif [[ "${target}" == *-freebsd* ]]; then

    #help find libpng for some reason
    export CPPFLAGS="-I${includedir}"

elif [[ "${target}" == *-apple* ]]; then

    #see https://trac.wxwidgets.org/ticket/19159 for issue tracker and https://github.com/wxWidgets/wxWidgets/pull/2354 for source of patch
    #this has since been merged into master, unsure about future status?
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/apple-add-bridging-commit.patch

    #fix missing symbols error - Undefined symbols for architecture x86_64: "___isPlatformVersionAtLeast". I think something to do with cross-compiling and we don't have access to XCode libraries??
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/osx-disable-builtin-platform-check.patch

fi

CONFIGURE_FLAGS=(--prefix=${prefix}
                --build=${MACHTYPE}
                --host=${target}
                --includedir=${includedir}
                --libdir=${libdir}
                --with-libiconv=ON
                --with-libcurl=ON
                --with-libpng=sys
                --with-libjpeg=sys
                --with-libtiff=sys
                --with-zlib=sys
                --with-expat=sys
                --disable-tests
                )

#wxWidgets has a couple of variants available (see "Supported Platforms" section https://www.wxwidgets.org/about/).
#On Unix platforms, we'll build the wxGTK variant, on MinGW platforms we'll build the wxMSW variant, and we'll build the OSX/OSCocoa variant on mac platforms.

if [[ "${target}" == *-apple* ]]; then
    CONFIGURE_FLAGS+=(--with-osx_cocoa)
    CONFIGURE_FLAGS+=(--with-macosx-version-min=10.12)
elif [[ "${target}" == *-mingw* ]]; then
    CONFIGURE_FLAGS+=(--with-msw)
else
    CONFIGURE_FLAGS+=(--with-gtk=3)
fi

./configure ${CONFIGURE_FLAGS[@]}

# Override insane script which tries to code sign **all** libraries in
# ${libdir}, seriously??
echo > change-install-names

make -j${nproc}
make install

install_license docs/preamble.txt docs/licence.txt docs/licendoc.txt docs/gpl.txt docs/lgpl.txt docs/xserver.txt
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms(; experimental = true))
filter!(p -> arch(p) != "armv6l", platforms)

# The products that we will ensure are always built
#[linux_naming_scheme, mingw_naming_scheme, macos_naming_scheme]
#see https://github.com/wxWidgets/wxWidgets/blob/master/docs/contributing/about-platform-toolkit-and-library-names.md for more info on naming
products = [
    ExecutableProduct("wxrc-$(version_mm)", :wxrc),
    LibraryProduct(["libwx_baseu-$(version_mm)","wxbase$(version_no_sep)u_gcc_custom", "libwx_baseu-$(version_mm)-x86_64-apple-darwin14", "libwx_baseu-$(version_mm)-aarch64-apple-darwin20"], :baseu),
    LibraryProduct(["libwx_baseu_net-$(version_mm)","wxbase$(version_no_sep)u_net_gcc_custom", "libwx_baseu_net-$(version_mm)-x86_64-apple-darwin14", "libwx_baseu_net-$(version_mm)-aarch64-apple-darwin20"], :baseu_net),
    LibraryProduct(["libwx_baseu_xml-$(version_mm)","wxbase$(version_no_sep)u_xml_gcc_custom", "libwx_baseu_xml-$(version_mm)-x86_64-apple-darwin14", "libwx_baseu_xml-$(version_mm)-aarch64-apple-darwin20"], :baseu_xml),
    LibraryProduct(gen_libnames("aui"), :aui),
    #LibraryProduct(gen_libnames("adv"), :adv),
    LibraryProduct(gen_libnames("core"), :core),
    LibraryProduct(gen_libnames("html"), :html),
    LibraryProduct(gen_libnames("propgrid"), :propgrid),
    LibraryProduct(gen_libnames("qa"), :qa),
    LibraryProduct(gen_libnames("ribbon"), :ribbon),
    LibraryProduct(gen_libnames("richtext"), :richtext),
    LibraryProduct(gen_libnames("stc"), :stc),
    LibraryProduct(gen_libnames("xrc"), :xrc)

]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="libpng_jll", uuid="b53b4c65-9356-5827-b1ea-8c7a1a84506f"))
    Dependency(PackageSpec(name="JpegTurbo_jll", uuid="aacddb02-875f-59d6-b918-886e6ef4fbf8"))
    Dependency(PackageSpec(name="Libtiff_jll", uuid="89763e89-9b03-5906-acba-b20f662cd828"))
    Dependency(PackageSpec(name="LibCURL_jll", uuid="deac9b47-8bc7-5906-a0fe-35ac56dc84c0"))
    Dependency(PackageSpec(name="GTK3_jll", uuid="77ec8976-b24b-556a-a1bf-49a033a670a6"))
    Dependency(PackageSpec(name="Zlib_jll", uuid = "83775a58-1f1d-513f-b197-d71354ab007a"))
    BuildDependency(PackageSpec(name="Xorg_xorgproto_jll", uuid="c4d99508-4286-5418-9131-c86396af500b"))
    Dependency("Expat_jll"; compat="2.2.10")
]

# Build the tarballs, and possibly a `build.jl` as well.
#wxMSW fails on gcc4 and 5, wxGTK works on everything
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6",preferred_gcc_version=v"6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
@testset "Result" begin
    fm = _build_dmip_formulation()
    let state = _CurrentState()
        result = @inferred Cerberus.Result(state, CONFIG)
        @test result.primal_bound == Inf
        @test result.dual_bound == -Inf
        @test result.termination_status == Cerberus.EARLY_TERMINATION
        @test result.total_node_count == 0
        @test result.total_simplex_iters == 0
        @test result.total_model_builds == 0
        @test result.total_warm_starts == 0
    end

    let state = _CurrentState(primal_bound = 12.4)
        result = @inferred Cerberus.Result(state, CONFIG)
        @test result.primal_bound == 12.4
        @test result.dual_bound == -Inf
        @test result.termination_status == Cerberus.EARLY_TERMINATION
        @test result.total_node_count == 0
        @test result.total_simplex_iters == 0
        @test result.total_model_builds == 0
        @test result.total_warm_starts == 0
    end
end

@testset "End-to-end" begin
    for (
        warm_start_strategy,
        model_reuse_strategy,
        model_builds,
        warm_starts,
    ) in [
        (Cerberus.NO_WARM_STARTS, Cerberus.NO_MODEL_REUSE, 3, 0),
        (Cerberus.NO_WARM_STARTS, Cerberus.REUSE_MODEL_ON_DIVES, 2, 0),
        (Cerberus.NO_WARM_STARTS, Cerberus.USE_SINGLE_MODEL, 1, 0),
        (Cerberus.WARM_START_WHEN_BACKTRACKING, Cerberus.NO_MODEL_REUSE, 3, 1),
        (
            Cerberus.WARM_START_WHEN_BACKTRACKING,
            Cerberus.REUSE_MODEL_ON_DIVES,
            2,
            1,
        ),
        (
            Cerberus.WARM_START_WHEN_BACKTRACKING,
            Cerberus.USE_SINGLE_MODEL,
            1,
            1,
        ),
        (Cerberus.WARM_START_WHENEVER_POSSIBLE, Cerberus.NO_MODEL_REUSE, 3, 2),
        (
            Cerberus.WARM_START_WHENEVER_POSSIBLE,
            Cerberus.REUSE_MODEL_ON_DIVES,
            2,
            2,
        ),
        (
            Cerberus.WARM_START_WHENEVER_POSSIBLE,
            Cerberus.USE_SINGLE_MODEL,
            1,
            2,
        ),
    ]
        config = Cerberus.AlgorithmConfig(
            warm_start_strategy = warm_start_strategy,
            model_reuse_strategy = model_reuse_strategy,
            silent = true,
        )
        fm = _build_dmip_formulation()
        result = Cerberus.optimize!(fm, config)
        @test result.primal_bound ≈ 1 / 21
        @test result.dual_bound ≈ 1 / 21
        @test result.best_solution ≈ [1.0, 20 / 21, 0.0]
        @test result.termination_status == Cerberus.OPTIMAL
        @test result.total_node_count == 3
        @test result.total_simplex_iters == 0
        @test result.total_elapsed_time_sec > 0
        @test result.total_model_builds == model_builds
        @test result.total_warm_starts == warm_starts
    end
end


text:
"
It is possible to have matrix product (operators / states) that are also periodic in the vertical direction
For examples, as fix points of statmech problems
These should be represented as respectively MultiLine{<:InfiniteMPO} / Multiline{<:InfiniteMPS}
"
struct Multiline{T}
    data::PeriodicArray{T,1}
end

Base.length(t::Multiline) = prod(size(t));
Base.size(t::Multiline) = (length(t.data),length(t.data[1]));
Base.size(t::Multiline,i) = size(t)[i];
Base.getindex(t::Multiline,i) = t.data[i];
Base.copy(t::Multiline) = Multiline(map(copy,t.data));
Multiline(t::AbstractArray) = Multiline(PeriodicArray(t));

Base.convert(::Vector,t::Multiline) = t.data.data;
Base.convert(::PeriodicArray,t::Multiline) = t.data;

Base.convert(::Multiline,v::AbstractArray) = Multiline(v);

#--- implementation of MPSMultiline
const MPSMultiline = Multiline{<:InfiniteMPS}

function MPSMultiline(pspaces::AbstractArray{S,2},Dspaces::AbstractArray{S,2};kwargs...) where S
    MPSMultiline(map(zip(circshift(Dspaces,(0,-1)),pspaces,Dspaces)) do (D1,p,D2)
        TensorMap(rand,Defaults.eltype,D1*p,D2)
    end; kwargs...)
end

function MPSMultiline(data::AbstractArray{T,2};kwargs...) where T<:GenericMPSTensor
    Multiline(PeriodicArray(map(1:size(data,1)) do row
        InfiniteMPS(data[row,:];kwargs...)
    end))
end
function MPSMultiline(data::AbstractArray{T,2},Cinit::AbstractArray{O,1};kwargs...) where {T<:GenericMPSTensor,O<:MPSBondTensor}
    Multiline(PeriodicArray(map(1:size(data,1)) do row
        InfiniteMPS(data[row,:],Cinit[row];kwargs...);
    end))
end
function Base.getproperty(psi::MPSMultiline,prop::Symbol)
    if prop == :AL
        return ALView(psi)
    elseif prop == :AR
        return ARView(psi)
    elseif prop == :AC
        return ACView(psi)
    elseif prop == :CR
        return CRView(psi)
    else
        return getfield(psi,prop)
    end
end

for f in (:l_RR, :l_RL, :l_LL, :l_LR)
    @eval $f(t::MPSMultiline,i,j = 1) = $f(t[i],j)
end

for f in (:r_RR, :r_RL, :r_LR,:r_LL)
    @eval $f(t::MPSMultiline,i,j = size(t,2)) = $f(t[i],j)
end

TensorKit.dot(a::MPSMultiline,b::MPSMultiline;kwargs...) = sum(dot.(a.data,b.data;kwargs...))

Base.convert(::Type{MPSMultiline},st::InfiniteMPS) = Multiline([st]);
Base.convert(::Type{InfiniteMPS},st::MPSMultiline) = st[1];
Base.eltype(t::MPSMultiline) = eltype(t[1]);
virtualspace(t::MPSMultiline,i::Int,j::Int) = virtualspace(t[i],j);


text:
using NetCDF

function transp_potential!(M, transp_file, time)
    tt = Float64.(ncread(transp_file,"TIME"))
    ind = argmin(abs.(tt .- 1e-3*time))
    psi = Float64.(ncread(transp_file,"PLFLX")[:,ind])
    epot = Float64.(ncread(transp_file,"EPOTNC")[:,ind])
    phi_itp = LinearInterpolation(psi,epot,extrapolation_bc=Flat())
    phi = phi_itp.(M.psi)
    phi_itp = CubicSplineInterpolation(M.psi,phi,extrapolation_bc=Flat())
    M.phi = phi_itp
    nothing
end


text:
s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => true)
@testset "test ac polar pf" begin
    # @testset "3-bus case" begin
    #     result = run_acdcpf("../test/data/case3.m", ACPPowerModel, ipopt_solver, setting = Dict("output" => Dict("branch_flows" => true)))
    #
    #     @test result["termination_status"] == LOCALLY_SOLVED
    #     @test isapprox(result["objective"], 0; atol = 1e-2)
    #
    #     @test isapprox(result["solution"]["gen"]["2"]["pg"], 1.600063; atol = 1e-3)
    #     @test isapprox(result["solution"]["gen"]["3"]["pg"], 0.0; atol = 1e-3)
    #
    #     @test isapprox(result["solution"]["bus"]["1"]["vm"], 1.10000; atol = 1e-3)
    #     @test isapprox(result["solution"]["bus"]["1"]["va"], 0.00000; atol = 1e-3)
    #     @test isapprox(result["solution"]["bus"]["2"]["vm"], 0.92617; atol = 1e-3)
    #     @test isapprox(result["solution"]["bus"]["3"]["vm"], 0.90000; atol = 1e-3)
    #
    #     @test isapprox(result["solution"]["dcline"]["1"]["pf"],  0.10; atol = 1e-5)
    #     @test isapprox(result["solution"]["dcline"]["1"]["pt"], -0.10; atol = 1e-5)
    #     @test isapprox(result["solution"]["dcline"]["1"]["qf"], -0.403045; atol = 1e-5)
    #     @test isapprox(result["solution"]["dcline"]["1"]["qt"],  0.0647562; atol = 1e-5)
    # end
    @testset "5-bus ac dc case" begin
        result = run_acdcpf("../test/data/case5_acdc.m", ACPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 0; atol = 1e-2)

        @test isapprox(result["solution"]["gen"]["1"]["pg"], 1.3494; atol = 1e-3)
        @test isapprox(result["solution"]["gen"]["2"]["pg"], 0.40; atol = 1e-3)

        @test isapprox(result["solution"]["bus"]["1"]["vm"], 1.06; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["1"]["va"], 0.00000; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["2"]["vm"], 1.00; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["3"]["vm"], 0.995; atol = 1e-3)

        @test isapprox(result["solution"]["convdc"]["2"]["pgrid"], -0.1954; atol = 1e-3)
        @test isapprox(result["solution"]["convdc"]["3"]["pdc"], 0.36421; atol = 1e-3)
        @test isapprox(result["solution"]["busdc"]["1"]["vm"], 1.008; atol = 1e-3)

    end
    @testset "5-bus ac dc case with 2 seperate ac grids" begin
        result = run_acdcpf("../test/data/case5_2grids.m", ACPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 0; atol = 1e-2)

        @test isapprox(result["solution"]["gen"]["1"]["pg"], 1.9326; atol = 1e-3)
        @test isapprox(result["solution"]["gen"]["4"]["pg"], 0.40; atol = 1e-3)

        @test isapprox(result["solution"]["bus"]["1"]["vm"], 1.06; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["1"]["va"], 0.00000; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["3"]["vm"], 0.987; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["7"]["va"], -0.0065; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["10"]["vm"], 0.972; atol = 1e-3)

        @test isapprox(result["solution"]["convdc"]["1"]["pgrid"], 0.6; atol = 1e-3)
        @test isapprox(result["solution"]["convdc"]["2"]["pdc"], 0.56872; atol = 1e-3)
        @test isapprox(result["solution"]["busdc"]["1"]["vm"], 1.015; atol = 1e-3)
    end
    # REMOVED for TRAVIS, otherwise case ok
    # @testset "24-bus rts ac dc case with three zones" begin
    #     result = run_acdcpf("../test/data/case24_3zones_acdc.m", ACPPowerModel, ipopt_solver; setting = s)
    #
    #     @test result["termination_status"] == LOCALLY_SOLVED
    #     @test isapprox(result["objective"], 0; atol = 1e-2)
    #
    #     @test isapprox(result["solution"]["gen"]["65"]["pg"], 1.419; atol = 1e-3)
    #     @test isapprox(result["solution"]["gen"]["65"]["qg"], -1.29648; atol = 1e-3)
    #
    #
    #     @test isapprox(result["solution"]["bus"]["101"]["vm"], 1.035; atol = 1e-3)
    #     @test isapprox(result["solution"]["bus"]["101"]["va"], -0.1389; atol = 1e-3)
    #     @test isapprox(result["solution"]["bus"]["205"]["vm"], 1.0318; atol = 1e-3)
    #     @test isapprox(result["solution"]["bus"]["301"]["vm"], 1.026266; atol = 1e-3)
    #
    #     @test isapprox(result["solution"]["convdc"]["2"]["pgrid"], -0.753; atol = 1e-3)
    #     @test isapprox(result["solution"]["convdc"]["3"]["pdc"], -1.37301; atol = 1e-3)
    #     @test isapprox(result["solution"]["busdc"]["5"]["vm"], 1.01731; atol = 1e-3)
    # end
end

@testset "test dc pf" begin
    # @testset "3-bus case" begin
    #     result = run_acdcpf("../test/data/case3.m", DCPPowerModel, ipopt_solver)
    #
    #     @test result["termination_status"] == LOCALLY_SOLVED
    #     @test isapprox(result["objective"], 0; atol = 1e-2)
    #
    #     @test isapprox(result["solution"]["gen"]["1"]["pg"], 1.54994; atol = 1e-3)
    #
    #     @test isapprox(result["solution"]["bus"]["1"]["va"],  0.00000; atol = 1e-5)
    #     @test isapprox(result["solution"]["bus"]["2"]["va"],  0.09147654582; atol = 1e-5)
    #     @test isapprox(result["solution"]["bus"]["3"]["va"], -0.28291891895; atol = 1e-5)
    # end
    @testset "5-bus ac dc case" begin
        result = run_acdcpf("../test/data/case5_acdc.m", DCPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 0; atol = 1e-2)

        @test isapprox(result["solution"]["gen"]["1"]["pg"], 1.2169; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["5"]["va"], -0.09289; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["3"]["va"], -0.0826; atol = 1e-3)
        @test isapprox(result["solution"]["convdc"]["2"]["pgrid"], -0.2831; atol = 1e-3)
        @test isapprox(result["solution"]["convdc"]["3"]["pdc"], -0.3385; atol = 1e-3)
    end
    @testset "24-bus rts ac dc case with three zones" begin
        result = run_acdcpf("../test/data/case24_3zones_acdc.m", DCPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 0; atol = 1e-2)
        @test isapprox(result["solution"]["gen"]["65"]["pg"], 1.419; atol = 1e-3)

        @test isapprox(result["solution"]["bus"]["119"]["va"], 0.17208; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["224"]["va"], 0.07803; atol = 1e-3)

        @test isapprox(result["solution"]["convdc"]["2"]["pgrid"], -0.753; atol = 1e-3)
        @test isapprox(result["solution"]["convdc"]["3"]["pdc"], 1.43579; atol = 1e-3)
    end
end

#TODO check for minimum loss solutions
# @testset "test soc pf" begin
#     @testset "5-bus ac dc case" begin
#         result = run_acdcpf("../test/data/case5_acdc.m", SOCWRPowerModel, ipopt_solver)
#
#         @test result["termination_status"] == LOCALLY_SOLVED
#         @test isapprox(result["objective"], 0; atol = 1e-2)
#
#         @test isapprox(result["solution"]["gen"]["1"]["pg"], 12.258; atol = 1e-3)
#         @test isapprox(result["solution"]["gen"]["1"]["qg"], 13.264; atol = 1e-3)
#
#         @test isapprox(result["solution"]["bus"]["1"]["vm"], 1.06; atol = 1e-3)
#         @test isapprox(result["solution"]["bus"]["1"]["va"], 0.00000; atol = 1e-3)
#         @test isapprox(result["solution"]["bus"]["3"]["vm"], 0.987; atol = 1e-3)
#         @test isapprox(result["solution"]["bus"]["7"]["va"], -0.0065; atol = 1e-3)
#         @test isapprox(result["solution"]["bus"]["10"]["vm"], 0.972; atol = 1e-3)
#
#         @test isapprox(result["solution"]["convdc"]["1"]["pgrid"], 0.6; atol = 1e-3)
#         @test isapprox(result["solution"]["convdc"]["2"]["pdc"], 0.56872; atol = 1e-3)
#         @test isapprox(result["solution"]["busdc"]["1"]["vm"], 1.015; atol = 1e-3)
#     end
#     @testset "24-bus rts ac dc case with three zones" begin
#         result = run_acdcpf("../test/data/case24_3zones_acdc.m", SOCWRPowerModel, ipopt_solver)
#
#         @test result["termination_status"] == LOCALLY_SOLVED
#         @test isapprox(result["objective"], 0; atol = 1e-2)
#
#         @test isapprox(result["solution"]["gen"]["65"]["pg"], 1.419; atol = 1e-3)
#         @test isapprox(result["solution"]["gen"]["65"]["qg"], 11.1026; atol = 1e-3)
#
#
#         @test isapprox(result["solution"]["bus"]["101"]["vm"], 1.035; atol = 1e-3)
#         @test isapprox(result["solution"]["bus"]["205"]["vm"], 0.88996; atol = 1e-3)
#         @test isapprox(result["solution"]["bus"]["301"]["vm"], 1.0202; atol = 1e-3)
#
#         @test isapprox(result["solution"]["convdc"]["2"]["pgrid"], -0.753; atol = 1e-3)
#         @test isapprox(result["solution"]["convdc"]["3"]["pdc"], -1.3671; atol = 1e-3)
#         @test isapprox(result["solution"]["busdc"]["5"]["vm"], 1.0004; atol = 1e-3)
#     end
# end
#
#
#
# @testset "test soc distflow pf_bf" begin
#     @testset "5-bus ac dc case" begin
#         result = run_acdcpf("../test/data/case5_acdc.m", SOCDFPowerModel, ipopt_solver)
#
#         @test result["termination_status"] == LOCALLY_SOLVED
#         @test isapprox(result["objective"], 0; atol = 1e-2)
#
#         @test result["solution"]["gen"]["1"]["pg"] >= 1.480
#
#         @test isapprox(result["solution"]["gen"]["2"]["pg"], 1.600063; atol = 1e-3)
#         @test isapprox(result["solution"]["gen"]["3"]["pg"], 0.0; atol = 1e-3)
#
#         @test isapprox(result["solution"]["bus"]["1"]["vm"], 1.09999; atol = 1e-3)
#         @test isapprox(result["solution"]["bus"]["2"]["vm"], 0.92616; atol = 1e-3)
#         @test isapprox(result["solution"]["bus"]["3"]["vm"], 0.89999; atol = 1e-3)
#
#         @test isapprox(result["solution"]["dcline"]["1"]["pf"],  0.10; atol = 1e-4)
#         @test isapprox(result["solution"]["dcline"]["1"]["pt"], -0.10; atol = 1e-4)
#     end
#     @testset "24-bus rts ac dc case with three zones" begin
#         result = run_acdcpf("../test/data/case24_3zones_acdc.m", SOCDFPowerModel, ipopt_solver)
#
#         @test result["termination_status"] == LOCALLY_SOLVED
#         @test isapprox(result["objective"], 0; atol = 1e-2)
#     end
# end
#
#
# @testset "test sdp pf" begin
#     #=
#     #seems to be having an issue on linux (04/02/18)
#     @testset "5-bus ac dc case" begin
#         result = run_pf("../test/data/case5_acdc.m", SDPWRMPowerModel, scs_solver)
#
#         @test result["termination_status"] == :Optimal
#         @test isapprox(result["objective"], 0; atol = 1e-2)
#
#         @test result["solution"]["gen"]["1"]["pg"] >= 1.480
#
#         @test isapprox(result["solution"]["gen"]["2"]["pg"], 1.600063; atol = 1e-3)
#         @test isapprox(result["solution"]["gen"]["3"]["pg"], 0.0; atol = 1e-3)
#
#         @test isapprox(result["solution"]["bus"]["1"]["vm"], 1.09999; atol = 1e-3)
#         @test isapprox(result["solution"]["bus"]["2"]["vm"], 0.92616; atol = 1e-3)
#         @test isapprox(result["solution"]["bus"]["3"]["vm"], 0.89999; atol = 1e-3)
#
#         @test isapprox(result["solution"]["dcline"]["1"]["pf"],  0.10; atol = 1e-4)
#         @test isapprox(result["solution"]["dcline"]["1"]["pt"], -0.10; atol = 1e-4)
#     end
#     =#
#     @testset "24-bus rts ac dc case with three zones" begin
#         result = run_acdcpf("../test/data/case24_3zones_acdc.m", SDPWRMPowerModel, scs_solver)
#
#         @test result["termination_status"] == :Optimal
#         @test isapprox(result["objective"], 0; atol = 1e-2)
#     end
# end


text:
using BinDeps

export download_datasets

function download_datasets(startfile=1)
    datafolder = getconfig("datafolder")
    if !isfile(datafolder)
        mkpath(datafolder)
    end

    WDPA_filename = "WDPA_WDOECM_search_feadb1b4f30799a6dc3ad0b16116d3530ec4a477898f4e10e097e2030e167128_shp"

    # tuples of (dataset_name, filename, url)
    datasets = [
        # https://globalwindatlas.info/api/gis/global/wind-speed/10
        # https://globalwindatlas.info/api/gis/global/wind-speed/50
        # https://globalwindatlas.info/api/gis/global/wind-speed/100
        # https://globalwindatlas.info/api/gis/global/wind-speed/150
        # https://globalwindatlas.info/api/gis/global/wind-speed/200
        ("Global Wind Atlas", "Global Wind Atlas v3 - 100m wind speed.tif",
            "https://chalmersuniversity.box.com/shared/static/wfr6dm9bcmj0mcqtdn0uimhg0otd4ht1.tif"),
        ("WDPA (protected areas):", "WDPA.zip",
            "https://d1gam3xoknrgr2.cloudfront.net/current/WDPA_WDOECM_search_feadb1b4f30799a6dc3ad0b16116d3530ec4a477898f4e10e097e2030e167128_shp.zip"),
            # "https://chalmersuniversity.box.com/shared/static/wn1kznvy7qh1issqcxdlsq64kgtkaayi.zip"),
        ("GADM (global administrative areas)", "gadm36.zip",
            "https://biogeo.ucdavis.edu/data/gadm3.6/gadm36_shp.zip"),
        ("NUTS (administrative areas in Europe)", "nuts-2016-01m.shp.zip",
            "https://ec.europa.eu/eurostat/cache/GISCO/distribution/v2/nuts/download/ref-nuts-2016-01m.shp.zip"),
        ("Land Cover", "Landcover - USGS MODIS.tif",
            "https://chalmersuniversity.box.com/shared/static/zm8zdkv1mz0wns0u77afkna95wbo7ggo.tif"),
        ("ETOPO1 Topography", "ETOPO1_Ice_c_geotiff.zip",
            "https://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/ice_surface/cell_registered/georeferenced_tiff/ETOPO1_Ice_c_geotiff.zip"),
        ("1-km downscaled spatial population scenarios (Gao et al), file 1/3", "temp_ssp1.tar",
            "http://www.cgd.ucar.edu/iam/modeling/data/SSP1_1km_netcdf.tar"),
        ("1-km downscaled spatial population scenarios (Gao et al), file 2/3", "temp_ssp2.tar",
            "http://www.cgd.ucar.edu/iam/modeling/data/SSP2_1km_netcdf.tar"),
        ("1-km downscaled spatial population scenarios (Gao et al), file 3/3", "temp_ssp3.tar",
            "http://www.cgd.ucar.edu/iam/modeling/data/SSP3_1km_netcdf.tar"),
        ("Global population & GDP", "temp_popgdp.zip",
            "https://osf.io/7jv3n/download"),
        ("WRI Global Power Plant Database", "WRI - Global Power Plant Database v1.10.zip",
            "https://chalmersuniversity.box.com/shared/static/ss6gycw7hf10e1fiicbxl5rgk08q5xr9.zip"),
            # switch to the official v1.2 link later (some plant cleanup is hardcoded for v1.1)
            # "http://datasets.wri.org/dataset/540dcf46-f287-47ac-985d-269b04bea4c6/resource/c240ed2e-1190-4d7e-b1da-c66b72e08858/download/globalpowerplantdatabasev120"),
        ("Time zone shape file", "timezones-with-oceans.shapefile.zip",
            "https://github.com/evansiroky/timezone-boundary-builder/releases/download/2019b/timezones-with-oceans.shapefile.zip"),
        ("Average monthly wind speeds 1979-2019", "era5monthlywind.h5",
            "https://chalmersuniversity.box.com/shared/static/otvb5nq0lz5ntqx0e65kocos2afo7gas.h5"),
        ("Average monthly solar insolation 1979-2019", "era5monthlysolar.h5",
            "https://chalmersuniversity.box.com/shared/static/jlpspmp9ou96hk7xno46rf79wg3d5hqc.h5"),
        ("Various smaller datasets", "Various_smaller_datasets.zip",
            "https://chalmersuniversity.box.com/shared/static/w3pmx4xhgorgd6jejv23gn4ycsnza8s6.zip")
    ]

    for (i, datasetinfo) in enumerate(datasets)
        i < startfile && continue
        name, filename, url = datasetinfo

        println("\nDownloading dataset $i: $name")
        download_progressbar(url, joinpath(datafolder, filename))
    end

    println("\nDownloads complete.")

    for datasetinfo in datasets
        name, filename, url = datasetinfo
        foldername, extension = splitext(filename)
        fullpath = joinpath(datafolder, filename)

        if extension in [".zip", ".tar"] && isfile(fullpath)
            println("\nUnpacking archive: $filename")
            unpack(fullpath, joinpath(datafolder, foldername), extension)
        end
    end

    function renameWDPAfiles(WDPAfolder)
        for filename in readdir(WDPAfolder)
            newname = replace(filename, WDPA_filename => "WDPA-shapefile")
            if newname != filename
                mv(joinpath(WDPAfolder, filename), joinpath(WDPAfolder, newname))
            end
        end
    end

    if startfile <= 2
        renameWDPAfiles(joinpath(datafolder, "WDPA"))
        for i = 0:2
            foldername = "WDPA-shapefile$i"
            println("\nUnpacking archive: $foldername.zip")
            unpack(joinpath(datafolder, "WDPA", "$foldername.zip"),
                    joinpath(datafolder, "WDPA", foldername), ".zip")
            renameWDPAfiles(joinpath(datafolder, "WDPA", foldername))
        end
    end

    if startfile <= 4
        filename = "NUTS_RG_01M_2016_4326_LEVL_3.shp.zip"
        println("\nUnpacking archive: $filename")
        unpack(joinpath(datafolder, "nuts-2016-01m.shp", filename), joinpath(datafolder, "nuts2016-level3"), ".zip")
        rm(joinpath(datafolder, "nuts-2016-01m.shp"), force=true, recursive=true)
    end

    println("\n\n\nCleaning up...")

    if startfile <= 6
        mv(joinpath(datafolder, "ETOPO1_Ice_c_geotiff", "ETOPO1_Ice_c_geotiff.tif"), joinpath(datafolder, "ETOPO1_Ice_c_geotiff.tif"))
        rm(joinpath(datafolder, "ETOPO1_Ice_c_geotiff"))
    end
    if startfile <= 9
        sspfolder = joinpath(datafolder, "SSP_1km")
        !isdir(sspfolder) && mkdir(sspfolder)
        for ssp = 1:3
            for y = 2010:10:2100
                sspfile = joinpath(datafolder, "temp_ssp$ssp", "SSP$(ssp)_1km", "ssp$(ssp)_total_$y.nc4")
                isfile(sspfile) && mv(sspfile, joinpath(sspfolder, "ssp$(ssp)_total_$y.nc4"))
            end
            rm(joinpath(datafolder, "temp_ssp$ssp"), force=true, recursive=true)
        end
    end
    if startfile <= 10
        mv(joinpath(datafolder, "temp_popgdp", "data"), joinpath(datafolder, "global_population_and_gdp"))
        rm(joinpath(datafolder, "temp_popgdp"))
    end
    if startfile <= 11
        mv(joinpath(datafolder, "WRI - Global Power Plant Database v1.10"), joinpath(datafolder, "tempWRI"))
        mv(joinpath(datafolder, "tempWRI", "WRI - Global Power Plant Database v1.10"),
            joinpath(datafolder, "WRI - Global Power Plant Database v1.10"))
        rm(joinpath(datafolder, "tempWRI"))
    end
    if startfile <= 15
        mixeddir = joinpath(datafolder, "Various_smaller_datasets")
        for file in readdir(mixeddir)
            mv(joinpath(mixeddir, file), joinpath(datafolder, file))
        end
        rm(mixeddir)
    end

    for datasetinfo in datasets
        name, filename, url = datasetinfo
        foldername, extension = splitext(filename)

        if extension in [".zip", ".tar"]
            fullpath = joinpath(datafolder, filename)
            isfile(fullpath) && rm(fullpath)
        end
    end
    println("\n\nDone.")
end

function download_progressbar(url::AbstractString, filename::AbstractString)
    println("Downloading to $filename...")
    curl = Base.find_curl()
    if curl == nothing
        # if curl is not on the system, fall back to a download without progress bar
        download(url, filename)
    else
        download_curl(curl, url, filename)
    end
end

function download_curl(curl_exe::AbstractString, url::AbstractString, filename::AbstractString)
    run(`$curl_exe --progress-bar -g -L -f --retry 5 -o $filename $url`, wait=true)
    return filename
end

function unpack(inputfilename, outputpath, extension)
    mkdir(outputpath)
    run(BinDeps.unpack_cmd(inputfilename, outputpath, extension, ""))
end

max_stars_repo_name: niklasjakobsson/GlobalEnergyGIS

max_issues_repo_name: niklasjakobsson/GlobalEnergyGIS

max_forks_repo_name: niklasjakobsson/GlobalEnergyGIS


text:
struct SimpleGame
    γ  # discount factor
    ℐ  # agents
    𝒜  # joint action space
    R  # joint reward function
end


text:
function pseudocharge(gridpos, Ls, atoms, YukawaK, epsil0)
# Create the pseudo-charge for Coulomb interaction
# $\rho_a = \sum_{N_{\text{atoms}}} - Z_j/\sqrt{2 \pi \sigma_i}$

Nsglb = length(gridpos);
Natoms = atoms.Natoms;

Z     = atoms.Z;
sigma = atoms.sigma;
R     = atoms.R;

# Calculate the total pseudo-charge
# NOTE: The pseudo-charge should not extend over twice the domain size!
rhoa = zeros(Nsglb);

    for j=1:Natoms
        d = R[j] .- gridpos;
        d = d - round.(Integer, d/Ls)*Ls;
        # Note: Z has minus sign in front of it!
        rhoa = rhoa - Z[j]./sqrt(2*pi*sigma[j]^2) .* (exp.(-0.5*(d./sigma[j]).^2 ));
    end

drhoa = zeros(Nsglb, Natoms);

    for j=1:Natoms
        d = R[j] .- gridpos;
        d = d - round.(Integer, d/Ls)*Ls;
        drhoa[:,j] = drhoa[:,j] - Z[j] ./sqrt(2*pi*sigma[j]^2)./sigma[j]^2 .* (exp.(-0.5*(d./sigma[j]).^2)).*d;
    end
    return rhoa, drhoa
end

max_stars_repo_name: Forgotten/modelscf

max_issues_repo_name: Forgotten/modelscf

max_forks_repo_name: Forgotten/modelscf


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

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

function __solve(prob::OptimizationProblem, opt::Evolutionary.AbstractOptimizer, data = DEFAULT_DATA;
                 cb = (args...) -> (false), maxiters = nothing,
                 progress = false, kwargs...)
    local x, cur, state

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

    cur, state = iterate(data)

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

    if !(isnothing(maxiters)) && maxiters <= 0.0
        error("The number of maxiters has to be a non-negative and non-zero number.")
    elseif !(isnothing(maxiters))
        maxiters = convert(Int, maxiters)
    end

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

    t0 = time()

    result = Evolutionary.optimize(_loss, prob.u0, opt, !isnothing(maxiters) ? Evolutionary.Options(;iterations = maxiters, callback = _cb, kwargs...)
                                                                        : Evolutionary.Options(;callback = _cb, kwargs...))
    t1 = time()

    SciMLBase.build_solution(prob, opt, Evolutionary.minimizer(result), Evolutionary.minimum(result); original=result)
end


text:
### A Pluto.jl notebook ###
# v0.12.18

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ e0e61a2a-5b42-11eb-3e21-a5af0c3ba8fb
begin
	import Pkg
	Pkg.add("InteractiveUtils")
 using Interact
using Images
	Pkg.add("WebIO")
	using WebIO
end

# ╔═╡ 53e94c18-5684-11eb-21e9-01510c841f58
md"""
This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
"""

# ╔═╡ fc21549e-5683-11eb-07bf-51233c3f5e9b

md"""
## My Notebook $~~~~~~~$ Jan 13 20201

For no brake spaces use tildes for

spaces between two dollar signs.

Leave a blank line for a line feed.

Math symbols are rendered between two $ signs.

[LaTex Help Html Link](https://kapeli.com/cheat_sheets/LaTeX_Math_Symbols.docset/Contents/Resources/Documents/index)

$\alpha$

$ \alpha = \beta $

$\mathscr David~ Marc~ Tan~ Creti$

$\mathfrak David~\mathfrak Marc~\mathfrak Tan~\mathfrak Creti$

$\mathcal David Marc Tan Creti$

$$\mathbb Greetings. $$

$$This~display~has~spaces~made~with~tildes!$$

---
*help*
__me2__
_me_

---

Hello from between horizontal rulers.

---

``\LaTeX\ \ \ \ A = \pi r^2   ``

``   3^6 pi  e^i\pi   ``
``                    ``
```math
f(a) = \frac{1}{2\pi}\int_{0}^{2\pi} (\alpha+R\cos(\theta))d\theta
```

```
  Here is a
   code block
    beween triple
     back      ticks.

```

[^1] Note one

[^Barringer] Reference Article

# Level One
## Level Two
### Level Three
#### Level Four
##### Level Five
###### Level Six

Here's a quote:
> Julia is a high-level, high-performance dynamic programming language for
> technical computing, with syntax that is familiar to users of other
> technical computing environments.

  + item one   two spaces before
    + two
      + three
        + four
          + five
    +
  + item two
  + item three

  1)
  2)
  3)
---
!!! note "Note"
    This is the content of the note.

!!! warning "Beware!"
    And this is another one.
!!! hint "MyHint"
    This is the content of the hint.

!!! help "Help!"
    How may we help you?

!!! red "Advice"
    A word to the wise is sufficient!
!!! whatever "Additional Help"
    Try again

    and again!

!!! divider ""

```javascript
var s = "JavaScript syntax highlighting";
    alert(s);
```

|----Column One---|----Column Two----|----Column Three----|
|:----------:      | :----------      |:------------:      |
| *Row* 2         | **Row** 2        |``3B``Column        |

"""

# ╔═╡ 58484c68-5bb9-11eb-2cea-6df30770eaaa
md"""
!["A pretty tiger"](https://upload.wikimedia.org/wikipedia/commons/5/56/Tiger.50.jpg)
"""

# ╔═╡ a634c51e-5bb9-11eb-185b-7777c0b5a778
md"""
!["A pretty tiger"](https://thumbs2.imgbox.com/c6/6f/T13ZjubW_t.jpg)
"""

# ╔═╡ 86e17588-5bbf-11eb-01c5-f19a6581d3c9
pwd()*"/Tiger.50.jpg"

# ╔═╡ 090ececa-5bfc-11eb-2c32-dbac0bd527b0
md"""
!["A pretty tiger"](https://drive.google.com/file/d/1osFq5AIooG8YaYRkKnuhg_qvEa4Arima/view?usp=sharing)
"""

# ╔═╡ 46304b82-5c0e-11eb-0444-97b86948cdfd
md"""
[Tiger](https://www.dropbox.com/s/anoj6agrz3fnebr/CornBowel1.jpg?dl=0)
"""

# ╔═╡ c402c28e-5c1c-11eb-1f89-e5dd0316bf75
dog = load("./Tiger.50.jpg")

# ╔═╡ e63d1158-5c1f-11eb-2cf7-0f14c4309fb5
philip_file = download("https://drive.google.com/file/d/1gMBxAVh9NzNsQztfQqD2aQgN5Mg8baUM/view?usp=download x.jpg")

# ╔═╡ f243c06e-5c1f-11eb-15e4-3506e26e1d9a
load(philip_file)

# ╔═╡ 3315fed2-59f4-11eb-065c-5f4f4aa3fd85
Text(sprint(versioninfo))

# ╔═╡ b61675e2-5a4b-11eb-1983-b96c71b4e0d1
myVariable = 65

# ╔═╡ 2ad843c6-5b00-11eb-3092-51f414ed4efe
myString = "My variable is"

# ╔═╡ 9cb89e92-59f6-11eb-14b4-0dc23f5c914d
s2 = "$myString: $myVariable"

# ╔═╡ 02b028e4-5b02-11eb-0fcf-9d8d07ff2817
s2

# ╔═╡ 6c34f1e4-5b04-11eb-21a6-7d61a928dba5
md"""
$myString
-->
$myVariable

"""

# ╔═╡ d69de236-5b41-11eb-1dbe-df913335e3f0
#textarea(hint="Enter name"; rows=3)

# ╔═╡ 5cd77430-5b9d-11eb-0661-eb3cb1b47151
Pkg.installed()

# ╔═╡ 7cc8cea0-5ba8-11eb-183a-2139b05c710b

# ╔═╡ Cell order:
# ╠═53e94c18-5684-11eb-21e9-01510c841f58
# ╠═e0e61a2a-5b42-11eb-3e21-a5af0c3ba8fb
# ╠═fc21549e-5683-11eb-07bf-51233c3f5e9b
# ╠═58484c68-5bb9-11eb-2cea-6df30770eaaa
# ╠═a634c51e-5bb9-11eb-185b-7777c0b5a778
# ╠═86e17588-5bbf-11eb-01c5-f19a6581d3c9
# ╠═090ececa-5bfc-11eb-2c32-dbac0bd527b0
# ╠═46304b82-5c0e-11eb-0444-97b86948cdfd
# ╠═c402c28e-5c1c-11eb-1f89-e5dd0316bf75
# ╠═e63d1158-5c1f-11eb-2cf7-0f14c4309fb5
# ╠═f243c06e-5c1f-11eb-15e4-3506e26e1d9a
# ╠═3315fed2-59f4-11eb-065c-5f4f4aa3fd85
# ╠═b61675e2-5a4b-11eb-1983-b96c71b4e0d1
# ╠═2ad843c6-5b00-11eb-3092-51f414ed4efe
# ╠═9cb89e92-59f6-11eb-14b4-0dc23f5c914d
# ╠═02b028e4-5b02-11eb-0fcf-9d8d07ff2817
# ╠═6c34f1e4-5b04-11eb-21a6-7d61a928dba5
# ╠═d69de236-5b41-11eb-1dbe-df913335e3f0
# ╠═5cd77430-5b9d-11eb-0661-eb3cb1b47151
# ╠═7cc8cea0-5ba8-11eb-183a-2139b05c710b

max_stars_repo_name: paradocs/hello-world

max_issues_repo_name: paradocs/hello-world

max_forks_repo_name: paradocs/hello-world


text:
include("euler/euler.jl")
using .BigIntegers: BigInteger
using .Numbers: digits_sum
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

compute(n::Int, p::Int)::Int = digits_sum(BigInteger(n) ^ p)

compute(2, 15)

compute(2, 1_000)

@benchmark compute(2, 1_000)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
# Copyright (c) 2020: Matthew Wilhelm & Matthew Stuber.
# This work is licensed under the Creative Commons Attribution-NonCommercial-
# ShareAlike 4.0 International License. To view a copy of this license, visit
# http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative
# Commons, PO Box 1866, Mountain View, CA 94042, USA.
#############################################################################
# Dynamic Bounds Base
# The underlying abstraction layer for DynamicBounds.
# See https://github.com/PSORLab/DynamicBoundsBase.jl
#############################################################################
# src/problem_types/ODERelaxProb.jl
# Defines structures to hold parametric ODEs and access functions.
#############################################################################

"""
AbstractODERelaxProblem

Abstract type for problems used to construct relaxations of parametric ODEs.
"""
abstract type AbstractODERelaxProblem <: AbstractDERelaxProblem end
const AODERP = AbstractODERelaxProblem

"""
$(TYPEDEF)

A structure used to hold a parametric ODEs problem.

$(TYPEDFIELDS)
"""
mutable struct ODERelaxProb{F,JX,JP,xType,K} <: AODERP
    "Right-hand side function."
    f::F
    "Jacobian of rhs w.r.t x."
    Jx!::JX
    "Jacobian of rhs w.r.t p."
    Jp!::JP
    "Initial condition function."
    x0::xType
    "Lower constant state bound storage."
    xL::Vector{Float64}
    "Upper constant state bound storage."
    xU::Vector{Float64}
    "Time span of integration."
    tspan::Tuple{Float64, Float64}
    "Support points: Points 'x' values will be queried."
    tsupports::Vector{Float64}
    "Point to relax/bound in over the decision space."
    p::Vector{Float64}
    "Lower bounds of the decision space."
    pL::Vector{Float64}
    "Upper bounds of the decision space."
    pU::Vector{Float64}
    "Problem has a user-defined Jacobian of the rhs w.r.t x"
    user_Jx::Bool
    "Problem has a user-defined Jacobian of the rhs w.r.t p"
    user_Jp::Bool
    "Constant user state bounds have been set"
    user_state_bnd::Bool
    "Variable user state bounds have been set"
    variable_state_bnd::Bool
    "State space dimension"
    nx::Int
    "Decision space dimension"
    np::Int
    "Optional keywork argument: Indicates the rhs function and the initial
    condition function take an additional argument `param`. For example,
    it is now `f!(dx, x, p, param)`."
    params::Vector{Float64}
    "Storage for polyhedral constraint, if any"
    polyhedral_constraint::Union{PolyhedralConstraint, Nothing}
    "Storage for the constant state bounds, if any"
    constant_state_bounds::Union{ConstantStateBounds, Nothing}
    "The support set if used"
    support_set::SupportSet{Float64}
    "Additional keyword arguments"
    kwargs::K
end

"""
ODERelaxProb(f, tspan, x0, pL, pU; kwargs...)

Constructor for basic parametric ODE problem. The `f` is rhs function. The integration time span is `tspan`. The initial
condition (which may be a function of p) is `x0`. The lower decision variable bounds are given by `pL` and the upper
decision variable bounds are given by `pU`. Other inputs are set via keyword arguments.
"""
function ODERelaxProb(f::F, tspan::Tuple{Float64, Float64}, x0::xType,
                      pL::Vector{Float64}, pU::Vector{Float64}; params::Vector{Float64} = Float64[],
                      xL::Vector{Float64} = Float64[], xU::Vector{Float64} = Float64[],
                      Jx! = nothing, Jp! = nothing, kwargs...) where {F,xType,tType,pBnd}

    if haskey(kwargs,:p)
        p::Vector{Float64} = kwargs[:p]
    else
        p = 0.5*(pL + pU)
    end
    support_set = haskey(kwargs, :support_set) ? kwargs[:support_set] : SupportSet()

    np::Int = length(p)

    @assert((xL !== nothing && xU !== nothing) || (xL === nothing && xU === nothing),
            "Both state bounds must be set or neither must be set")
    @assert((pL !== nothing && pU !== nothing), "Box constraints on parameters must be defined.")
    @assert(length(pL) === np === length(pU), "The parameter vector and bounds must have the same length.")

    nx::Int64 = haskey(kwargs, :nx) ? kwargs[:nx] : length(x0(pL))
    tsupports::Vector{Float64} =  haskey(kwargs, :tsupports) ? kwargs[:tsupports] : Float64[]

    polyhedral = haskey(kwargs, :polyhedral_constraint) ?  kwargs[:polyhedral_constraint] : nothing
    constantstate = haskey(kwargs, :constant_state_bounds) ? kwargs[:constant_state_bounds] : nothing

    user_state_bnd::Bool = false
    variable_state_bnd = !isempty(xL) && !isempty(xU)
    user_state_bnd |= variable_state_bnd

    if !isempty(xL) && !isempty(xU)
        user_state_bnd |= true
        @assert length(xL) == length(xU)
    end

    # compute Jacobians
    x0s = zeros(Real, nx)
    Jx0 = zeros(Real, nx, nx)
    Jp0 = zeros(Real, nx, np)
    user_Jx!::Bool = Jx! !== nothing
    user_Jp!::Bool = Jp! !== nothing
    if user_Jx! && user_Jp!
        JxWrap = wrapfun_iip(Jx!, (Jx0, x0s, p, tspan[1],))
        JpWrap = wrapfun_iip(Jx!, (Jp0, x0s, p, tspan[1],))
        fwrap = ODEFunction(f; jac = JxWrap, paramjac = JpWrap)
    elseif user_Jx!
        JxWrap = wrapfun_iip(Jx!, (Jx0, x0s, p, tspan[1],))
        fwrap = ODEFunction(f; jac = JxWrap)
    elseif user_Jp!
        JpWrap = wrapfun_iip(Jp!, (Jp0, x0s, p, tspan[1],))
        fwrap = ODEFunction(f; paramjac = JpWrap)
    else
        fwrap = ODEFunction(f)
    end

    return ODERelaxProb(fwrap, Jx!, Jp!, x0, xL, xU,
                        tspan, tsupports, p, pL, pU, user_Jx!, user_Jp!,
                        user_state_bnd, variable_state_bnd, nx, np, params,
                        polyhedral, constantstate, support_set, kwargs)
end

supports(::ODERelaxProb, ::HasStateBounds) = true
supports(::ODERelaxProb, ::HasConstantStateBounds) = true
supports(::ODERelaxProb, ::HasVariableStateBounds) = true
supports(::ODERelaxProb, ::HasUserJacobian) = true
supports(::ODERelaxProb, ::ConstantStateBounds) = true
supports(::ODERelaxProb, ::PolyhedralConstraint) = true
supports(::ODERelaxProb, ::SupportSet) = true

get(x::ODERelaxProb, t::HasStateBounds)::Bool = x.user_state_bnd
get(x::ODERelaxProb, t::HasConstantStateBounds)::Bool = ~(isempty(x.xL) && isempty(x.xU))
function get(x::ODERelaxProb, t::HasVariableStateBounds)::Bool
    return x.variable_state_bnd
end
function get(x::ODERelaxProb, t::HasUserJacobian)::Bool
    return x.user_Jx || x.user_Jp
end

function get(x::ODERelaxProb, t::ConstantStateBounds)
    return x.constant_state_bounds
end
function set!(x::ODERelaxProb, bnds::ConstantStateBounds)
    x.constant_state_bounds = bnds
    return
end

function get(x::ODERelaxProb, invariant::PolyhedralConstraint)
    return x.polyhedral_constraint
end
function set!(x::ODERelaxProb, invariant::PolyhedralConstraint)
    x.polyhedral_constraint = invariant
    return
end

function get(x::ODERelaxProb, support_set::SupportSet{Float64})
    return x.support_set
end
function set!(x::ODERelaxProb, support_set::SupportSet{Float64})
    x.support_set = support_set
    return
end

function get(x::ODERelaxProb, t::ConstantParameterValue)
    return x.params[t.i]
end
function set!(x::ODERelaxProb, t::ConstantParameterValue, v)
    x.params[t.i] = v
    return
end
function getall(x::ODERelaxProb, t::ConstantParameterValue)
    return x.params
end
function getall!(out, x::ODERelaxProb, t::ConstantParameterValue)
    out .= x.params
    return
end
function setall!(x::ODERelaxProb, t::ConstantParameterValue, v)
    x.params .= v
    return
end

mutable struct ODELocalIntegrator{N}
    problem
    ode_problem
    sensitivity_problem
    integrator
    p::Vector{Float64}
    pduals::Vector{Dual{Nothing,Float64,N}}
    x0::Vector{Float64}
    x::Matrix{Float64}
    dxdp::Vector{Matrix{Float64}}
    x0duals::Vector{Dual{Nothing,Float64,N}}
    user_t::Vector{Float64}
    integrator_t::Vector{Float64}
    local_t_dict_flt::Dict{Float64,Int64}
    local_t_dict_indx::Dict{Int64,Int64}
    abs_tol::Float64
    rel_tol::Float64
    function ODELocalIntegrator{N}(prob::ODERelaxProb, integrator; kwargs...) where N
        d = new()
        d.abs_tol = 1E-9
        d.rel_tol = 1E-8
        d.problem = prob
        d.integrator = integrator
        d.ode_problem = ODEProblem(prob.f, zeros(Float64, prob.nx),
                                   prob.tspan, prob.p)
        d.sensitivity_problem = ODEForwardSensitivityProblem(prob.f,
                                                             zeros(Float64, prob.nx),
                                                             prob.tspan,
                                                             prob.p)
        d.x0 = zeros(Float64, prob.nx)
        dxdp = Matrix{Float64}[]
        for i = 1:prob.np
            push!(dxdp, zeros(Float64, prob.nx, length(prob.support_set.s)))
        end
        d.x = zeros(prob.nx, length(prob.support_set.s))
        d.dxdp = dxdp
        d.p = copy(prob.p)
        d.pduals = seed_duals(Val(N),prob.p)
        d.x0duals = fill(Dual{Nothing}(0.0,
                                      single_seed(Partials{N, Float64}, Val(1))),
                                      (prob.nx,))
        support_set = get(prob, SupportSet())
        d.user_t = copy(support_set.s)
        d.integrator_t = copy(support_set.s)
        d.local_t_dict_flt = Dict{Float64,Int64}()
        d.local_t_dict_indx = Dict{Int64,Int64}()
        for (i,s) in enumerate(d.user_t)
            d.local_t_dict_flt[s] = i
            d.local_t_dict_indx[i] = i
        end
        return d
    end
end
function ODELocalIntegrator(prob::ODERelaxProb, integrator)
    ODELocalIntegrator{prob.np}(prob, integrator)
end

function integrate!(::Val{true}, d::AbstractODERelaxIntegrator, p::ODERelaxProb)

    local_prob_storage = get(d, LocalIntegrator())
    np = get(d, ParameterNumber())
    nx = get(d, StateNumber())
    if size(local_prob_storage.x0, 1) != nx + nx*np
        local_prob_storage.x0 = zeros(nx + nx*np)
    end
    for i = 1:nx
        local_prob_storage.x0[i] = local_prob_storage.x0duals[i].value
        for j = 1:np
            local_prob_storage.x0[(nx + j + (i-1)*np)] = local_prob_storage.x0duals[i].partials[j]
        end
    end

    local_prob_storage.sensitivity_problem = remake(local_prob_storage.sensitivity_problem,
                                            u0 = local_prob_storage.x0,
                                            p = local_prob_storage.p)

    solution = solve(local_prob_storage.sensitivity_problem,
                     local_prob_storage.integrator,
                     saveat = local_prob_storage.user_t,
                     abstol = local_prob_storage.abs_tol,
                     tstops = local_prob_storage.user_t,
                     reltol = local_prob_storage.rel_tol)

    x, dxdp = extract_local_sensitivities(solution)

    new_length = size(x, 2)
    prior_length = length(local_prob_storage.integrator_t)
    if new_length == prior_length
        local_prob_storage.integrator_t .= solution.t
    else
        local_prob_storage.x = zeros(nx, new_length)
        for i = 1:np
            local_prob_storage.dxdp[i] = zeros(nx, new_length)
        end
        local_prob_storage.integrator_t = solution.t
    end
    local_prob_storage.x .= x
    for i = 1:np
        local_prob_storage.dxdp[i] .= dxdp[i]
    end

    return solution.t
end

function integrate!(::Val{false}, d::AbstractODERelaxIntegrator, p::ODERelaxProb)

    local_prob_storage = get(d, LocalIntegrator())
    np = get(d, ParameterNumber())
    nx = get(d, StateNumber())

    if size(local_prob_storage.x0, 1) != nx
        local_prob_storage.x0 = zeros(nx)
    end
    for i = 1:nx
        local_prob_storage.x0[i] = local_prob_storage.x0duals[i].value
    end

    local_prob_storage.ode_problem = remake(local_prob_storage.ode_problem,
                                            u0 = local_prob_storage.x0,
                                            p = local_prob_storage.p)

    solution = solve(local_prob_storage.ode_problem,
                     local_prob_storage.integrator,
                     saveat = local_prob_storage.user_t,
                     abstol = local_prob_storage.abs_tol,
                     tstops = local_prob_storage.user_t,
                     reltol = local_prob_storage.rel_tol)

    x = solution.u

    new_length = length(x)
    prior_length = length(local_prob_storage.integrator_t)
    if new_length == prior_length
        local_prob_storage.integrator_t .= solution.t
    else
        local_prob_storage.x = zeros(nx, new_length)
        local_prob_storage.integrator_t = solution.t
    end
    for i = 1:new_length
        local_prob_storage.x[:,i] .= x[i]
    end

    return solution.t
end

function integrate!(d::AbstractODERelaxIntegrator, p::ODERelaxProb)

    local_prob_storage = get(d, LocalIntegrator())::ODELocalIntegrator

    getall!(local_prob_storage.p, d, ParameterValue())
    local_prob_storage.pduals .= seed_duals(Val(length(local_prob_storage.p)), local_prob_storage.p)
    local_prob_storage.x0duals = p.x0(d.local_problem_storage.pduals)

    solution_t = integrate!(Val(get(d, LocalSensitivityOn())), d, p)

    empty!(local_prob_storage.local_t_dict_flt)
    empty!(local_prob_storage.local_t_dict_indx)

    for (tindx, t) in enumerate(solution_t)
        local_prob_storage.local_t_dict_flt[t] = tindx
    end

    if !isempty(local_prob_storage.user_t)
        next_support_time = local_prob_storage.user_t[1]
        supports_left = length(local_prob_storage.user_t)
        loc_count = 1
        for (tindx, t) in enumerate(solution_t)
            if t == next_support_time
                local_prob_storage.local_t_dict_indx[loc_count] = tindx
                loc_count += 1
                supports_left -= 1
                if supports_left > 0
                    next_support_time = local_prob_storage.user_t[loc_count]
                end
            end
        end
    end

    return nothing
end

function get_val_loc_local(t::AbstractODERelaxIntegrator, index::Int64, time::Float64)

    local_prob_storage = get(t, LocalIntegrator())::ODELocalIntegrator

    (index <= 0 && time == -Inf) && error("Must set either index or time.")
    if index > 0
        return local_prob_storage.local_t_dict_indx[index]
    end
    local_prob_storage.local_t_dict_flt[time]
end

function get(out::Vector{Float64}, t::AbstractODERelaxIntegrator, v::Value)
    local_prob_storage = get(t, LocalIntegrator())::ODELocalIntegrator
    val_loc = get_val_loc_local(t, v.index, v.time)
    out .= local_prob_storage.x[:, val_loc]
    return
end

function getall!(out::Array{Float64,2}, t::AbstractODERelaxIntegrator, v::Value)
    local_prob_storage = get(t, LocalIntegrator())::ODELocalIntegrator
    copyto!(out, local_prob_storage.x)
    return
end

function getall!(out::Vector{Array{Float64,2}}, t::AbstractODERelaxIntegrator, g::Gradient{Nominal})
    local_prob_storage = get(t, LocalIntegrator())::ODELocalIntegrator
    for i = 1:get(t, ParameterNumber())
        copyto!(out[i], local_prob_storage.dxdp[i])
    end
    return
end


text:
# ---
# title: 1224. Maximum Equal Frequency
# id: problem1224
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Hash Table
# link: <https://leetcode.com/problems/maximum-equal-frequency/description/>
# hidden: true
# ---
#
# Given an array `nums` of positive integers, return the longest possible length
# of an array prefix of `nums`, such that it is possible to remove **exactly
# one** element from this prefix so that every number that has appeared in it
# will have the same number of occurrences.
#
# If after removing one element there are no remaining elements, it's still
# considered that every appeared number has the same number of ocurrences (0).
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [2,2,1,1,5,3,3,5]
#     Output: 7
#     Explanation: For the subarray [2,2,1,1,5,3,3] of length 7, if we remove nums[4]=5, we will get [2,2,1,1,3,3], so that each number will appear exactly twice.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1,1,1,2,2,2,3,3,3,4,4,4,5]
#     Output: 13
#
#
# **Example 3:**
#
#
#
#     Input: nums = [1,1,1,2,2,2]
#     Output: 5
#
#
# **Example 4:**
#
#
#
#     Input: nums = [10,2,8,9,3,8,1,5,2,3,7,6]
#     Output: 8
#
#
#
#
# **Constraints:**
#
#   * `2 <= nums.length <= 10^5`
#   * `1 <= nums[i] <= 10^5`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
mutable struct LQSolver{M <: AbstractLinearQuadraticModel, S <: AbstractMathProgSolver}
    lqmodel::M
    optimsolver::S

    function (::Type{LQSolver})(model::JuMP.Model,optimsolver::AbstractMathProgSolver)
        lqmodel = LinearQuadraticModel(optimsolver)
        solver = new{typeof(lqmodel),typeof(optimsolver)}(lqmodel,optimsolver)
        loadproblem!(solver.lqmodel,loadLP(model)...)
        setvartype!(lqmodel, model.colCat)
        return solver
    end
end

function (solver::LQSolver)(x₀::AbstractVector)
    if applicable(setwarmstart!,solver.lqmodel,x₀)
        try
            setwarmstart!(solver.lqmodel,x₀)
        catch
            warn("Could not set warm start for some reason")
        end
    end
    optimize!(solver.lqmodel)

    optimstatus = status(solver)
    !(status(solver) ∈ [:Optimal,:Infeasible,:Unbounded]) && error("LP could not be solved, returned status: ", optimstatus)

    return nothing
end

function getsolution(solver::LQSolver)
    optimstatus = status(solver)
    if optimstatus == :Optimal
        return getsolution(solver.lqmodel)
    else
        error("LP was not solved to optimality, returned status: ", optimstatus)
    end
end

function getobjval(solver::LQSolver)
    optimstatus = status(solver)
    if optimstatus == :Optimal
        return getobjval(solver.lqmodel)
    elseif optimstatus ==:Infeasible
        return Inf
    elseif status == :Unbounded
        solver.obj = -Inf
    else
        error("LP could not be solved, returned status: ", optimstatus)
    end
end

function getredcosts(solver::LQSolver)
    cols = numvar(solver.lqmodel)
    optimstatus = status(solver)
    if optimstatus == :Optimal
        return try
            getreducedcosts(solver.lqmodel)[1:cols]
        catch
            fill(NaN, cols)
        end
    else
        error("LP was not solved to optimality. Return status: ", optimstatus)
    end
end

function getduals(solver::LQSolver)
    rows = numconstr(solver.lqmodel)
    optimstatus = status(solver)
    if optimstatus == :Optimal
        return try
            getconstrduals(solver.lqmodel)[1:rows]
        catch
            fill(NaN, rows)
        end
    elseif optimstatus ==:Infeasible
        return try
            getinfeasibilityray(solver.lqmodel)[1:rows]
        catch
            fill(NaN, rows)
        end
    else
        error("LP was not solved to optimality, and the model was not infeasible. Returned status: ", optimstatus)
    end
end

status(solver::LQSolver) = MathProgBase.SolverInterface.status(solver.lqmodel)

function loadLP(m::JuMP.Model)
    l = m.colLower
    u = m.colUpper

    # Build objective
    # ==============================
    c = JuMP.prepAffObjective(m)
    c *= m.objSense == :Min ? 1 : -1

    # Build constraints
    # ==============================
    # Non-zero row indices
    I = Vector{Int}()
    # Non-zero column indices
    J = Vector{Int}()
    # Non-zero values
    V = Vector{Float64}()
    # Lower constraint bound
    lb = zeros(Float64,length(m.linconstr))
    # Upper constraint bound
    ub = zeros(Float64,length(m.linconstr))

    for (i,constr) in enumerate(m.linconstr)
        coeffs = constr.terms.coeffs
        vars = constr.terms.vars
        @inbounds for (j,var) = enumerate(vars)
            if var.m == m
                push!(I,i)
                push!(J,var.col)
                push!(V,coeffs[j])
            end
        end
        lb[i] = constr.lb
        ub[i] = constr.ub
    end
    A = sparse(I,J,V,length(m.linconstr),m.numCols)

    return A,l,u,c,lb,ub,:Min
end

function getVariableDuals(solver::LQSolver,d::AbstractVector)
    l = getvarLB(solver.lqmodel)
    v = zero(l)
    u = getvarUB(solver.lqmodel)
    w = zero(u)
    for i = 1:length(l)
        if l[i] > -Inf
            if u[i] == Inf
                v[i] = d[i]
            else
                if solver.x[i] == l[i]
                    v[i] = d[i]
                elseif solver.x[i] == u[i]
                    w[i] = d[i]
                else
                    continue
                end
            end
        else
            if u[i] < Inf
                w[i] = d[i]
            end
        end
    end
    return v,w
end

function getConstrDuals(solver::LQSolver)
    λ = getconstrduals(solver.lqmodel)
    Ax = getconstrsolution(solver.lqmodel)

    lb = getconstrLB(solver.lqmodel)
    λl = zero(lb)
    ub = getconstrUB(solver.lqmodel)
    λu = zero(ub)
    for i = 1:length(lb)
        if lb[i] > -Inf
            if ub[i] == Inf
                λl[i] = λ[i]
            else
                if Ax[i] == lb[i]
                    λl[i] = λ[i]
                elseif Ax[i] == ub[i]
                    λu = λ[i]
                else
                    continue
                end
            end
        else
            if ub[i] < Inf
                λu[i] = λ[i]
            end
        end
    end
    return λl,λu
end


text:
# include("general.jl")
# include("reducers.jl")

function response_vanilla(mu::Array{<:Real}, omega::Array{<:Real}, lambda::Array{<:Real}, i::Int64, j::Int64, epsilon::Float64)::Array{Float64}
    lambda_response = copy(lambda)
    lambda_response[i] = (mu[i] * omega[i] + mu[j] * omega[j] - epsilon * omega[j]) / (omega[i] + omega[j])
    lambda_response[j] = lambda_response[i] + epsilon
    return copy(lambda_response)
end

function response_spectral(mu::Array{<:Real}, omega::Array{<:Real}, lambda::Array{<:Real}, i::Int64, j::Int64, epsilon::Float64, laplacian::Array{<:Real}, R::Real)::Array{Float64}
    lambda_vanilla = response_vanilla(mu, omega, lambda, i, j, epsilon)
end

function response_oracle()
    # call either reponse_vanilla or response_spectral
end

max_stars_repo_name: tomaskocak/BAI

max_issues_repo_name: tomaskocak/BAI

max_forks_repo_name: tomaskocak/BAI


text:
# load code
include("SDE_su2.jl")
include("BSE_su2_s.jl")
include("BSE_su2_t.jl")
include("BSE_su2_u.jl")


text:
# Copyright 2022 RelationalAI, Inc.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

import ConfParser
using Base: Filesystem

mutable struct Config
    region::Union{String,Nothing}
    scheme::Union{String,Nothing}
    host::Union{String,Nothing}
    port::Union{String,Nothing}
    credentials::Union{Credentials,Nothing}
end

function _load_stanza(fname::AbstractString, profile::AbstractString)
    fname = Filesystem.expanduser(fname)
    conf = ConfParser.ConfParse(fname)
    ConfParser.parse_conf!(conf)
    return ConfParser.retrieve(conf, profile)
end

# Note, ConfParser returns Dict{String => Vector{T}}
_get_value(d::Dict, k::String) = first(get(d, k, (nothing,)))

function load_config(; fname = nothing, profile = nothing)::Config
    isnothing(fname) && (fname = "~/.rai/config")
    isnothing(profile) && (profile = "default")
    stanza = _load_stanza(fname, profile)
    region = _get_value(stanza, "region")
    scheme = _get_value(stanza, "scheme")
    host = _get_value(stanza, "host")
    port = _get_value(stanza, "port")
    client_id = _get_value(stanza, "client_id")
    client_secret = _get_value(stanza, "client_secret")
    credentials = nothing
    if !isnothing(client_id) && !isnothing(client_secret)
        client_credentials_url = _get_value(stanza, "client_credentials_url")
        credentials = ClientCredentials(client_id, client_secret, client_credentials_url)
    end
    return Config(region, scheme, host, port, credentials)
end

max_stars_repo_name: RelationalAI/rai-sdk-julia

max_issues_repo_name: RelationalAI/rai-sdk-julia

max_forks_repo_name: RelationalAI/rai-sdk-julia


text:
# Assign ranks to values but keep the rank the same if the values are within
# tolerance of each other.
function assign_ranks_within_tolerance(values; by = (x) -> x, tolerance = 1e-5, rev = false)

  perm = sortperm(values, by = by, rev = rev)
  ranked = Any[]
  rank = 1
  prev = by(values[perm[1]])
  num_of_this_rank = 0

  for i in 1:length(perm)
    r = values[perm[i]]
    v = by(r)
    if abs(prev - v) > tolerance
      rank += num_of_this_rank
      num_of_this_rank = 1
    else
      num_of_this_rank += 1
    end
    push!(ranked, (rank, r, v))
    prev = v
  end

  ranked

end