text:
runmodel("2DEG")

#ForcingMult = 1.39; fertilization = 0.59;
#ForcingMult = 1.390131560400832; fertilization = 0.5877118470187332;

#annualEmissions = getscenario("2DEG");
#p = ClimateParams(0.01, 0.8, 0.0, 0.0, 1.0, 2.0, 278.0);
#calibrateforcing!(p);
#calibratefertilization!(annualEmissions, p);
annualEmissions = getscenario("RCP45");
results, p, rcp = solveclimate(annualEmissions);
printresults(2010:10:2100, results, p, annualEmissions, rcp);

@code_warntype ClimateCalculator.historicdata()
@code_warntype getscenario("2DEG")
# next line (and similar) requires change of function declaration to "global function"
#@code_warntype ClimateCalculator_refactor.forcingerror(ForcingMult, lambda)
@code_warntype calibrateforcing!(p, "RCP45")
#@code_warntype ClimateCalculator_refactor.concentrationerror(fertilization, lambda, ForcingMult, annualEmissions)
@code_warntype calibratefertilization!(annualEmissions, p, "RCP45")
@code_warntype ClimateCalculator.initclimate(p, 2010, false)
@code_warntype solveclimate(annualEmissions, p)
@code_warntype solveclimate(annualEmissions, p, 2010, 2100, true, "RCP45")
@code_warntype solveclimate(annualEmissions, p, 2010, 2100, false, "RCP45")
@code_warntype ClimateCalculator_refactor.printresults(lambda, ForcingMult, fertilization, annualEmissions, results, false)
@code_warntype runmodel("2DEG",0.8; timestep=0.1, firstyear=2010, lastyear=2100, usecache=true, output=true);

using BenchmarkTools
annualEmissions = getscenario("RCP45");
results, p, rcp = solveclimate(annualEmissions);
p.timestep = 0.1
@benchmark solveclimate(annualEmissions, p)

using BenchmarkTools
annualEmissions = getscenario("RCP45");
results, p, rcp = solveclimate(annualEmissions, usecache=false);
p.timestep = 0.1
@benchmark solveclimate(annualEmissions, p, 2010, 2100, false)

@benchmark solveclimate(annualEmissions, timestep=0.1)
@benchmark solveclimate(annualEmissions, p, firstyear=2010, lastyear=2100)
@benchmark runmodel("2DEG", output=false)

lambda, annualEmissions = getscenario("2DEG");
p = ClimateParams(0.1, lambda, 0.0, 0.0, 1.0, 2.0, 278.0);
calibrateforcing!(p);
calibratefertilization!(annualEmissions, p);
s = ClimateState();
ClimateCalculator.init_temperatures!(s);
ClimateCalculator.init_oceancarbon!(s, p);
ClimateCalculator.init_biocarbon!(s, p);
ClimateCalculator.init_concentration_and_forcing!(s, p);
@code_warntype ClimateCalculator.init_temperatures!(s)
@code_warntype ClimateCalculator.oceancarbon_init!(s, p)
@code_warntype ClimateCalculator.biocarbon_init!(s, p)
@code_warntype ClimateCalculator.concentration_RF_init!(s, p)
@code_warntype ClimateCalculator.backgroundCO2_iteration(0.0, equilibriumCO2)
@code_warntype ClimateCalculator.backgroundCO2_solve(equilibriumCO2)
@code_warntype ClimateCalculator.calc_netfluxOcean(CO2concentration, 0.0, Temp[1])

const GAS3 = [:CO2, :CH4, :N2O]
t = 2017.55
const y0 = 1765
iyear(t,f) = f(Int,t) - y0 + 1
iyear(t) = iyear(t,round)
interpolate(t, annualvector) =  (1-(t-floor(t))) * annualvector[iyear(t,floor)] +
                                    (t-floor(t)) * annualvector[iyear(t,ceil)]
interpolate(t, annualdict, gas) = (1-(t-floor(t))) * annualdict[iyear(t,floor)][gas] +
                                      (t-floor(t)) * annualdict[iyear(t,ceil)][gas]
emissions = Dict(   :CO2 => interpolate(t, annualEmissions[:CO2]),
                    :CH4 => interpolate(t, annualEmissions[:CH4]),
                    :N2O => interpolate(t, annualEmissions[:N2O]))
ClimateCalculator.carbonbalance!(emissions, s, p)
ClimateCalculator.oceancarbon!(s, p)
ClimateCalculator.biocarbon!(s, p)
ClimateCalculator.concentrations!(emissions, s, p)
ClimateCalculator.radiativeforcing!(t, s, p)
ClimateCalculator.temperatures!(s, p)
@code_warntype interpolate(t, annualEmissions[:CO2])
@code_warntype interpolate(t, [Dict(:CO2 => 23.0)], :CO2)
@code_warntype ClimateCalculator.carbonbalance!(emissions, s, p)
@code_warntype ClimateCalculator.oceancarbon!(s, p)
@code_warntype ClimateCalculator.biocarbon!(s, p)
@code_warntype ClimateCalculator.concentrations!(emissions, s, p)
@code_warntype ClimateCalculator.radiativeforcing!(t, s, p)
@code_warntype ClimateCalculator.temperatures!(s, p)

@benchmark ClimateCalculator.carbonbalance!(emissions, s, p)
@benchmark ClimateCalculator.oceancarbon!(s, p)
@benchmark ClimateCalculator.biocarbon!(s, p)
@benchmark ClimateCalculator.concentrations!(emissions, s, p)
@benchmark ClimateCalculator.radiativeforcing!(t, s, p)
@benchmark ClimateCalculator.temperatures!(s, p)

results = [solveclimate(annualEmissions, p, timestep=[0.25,0.1,0.01,0.001][i]) for i=1:4];
pp=plot();for i=1:4 pp=plot!(1765:2500, get.(getfield.(results[i], :Concentration),:CO2,0.0)  );end;display(pp)
pp=plot();for i=1:4 pp=plot!(1765:2500, getfield.(results[i], :Temp_global)  );end;display(pp)

max_stars_repo_name: niclasmattsson/ClimateCalculator

max_issues_repo_name: niclasmattsson/ClimateCalculator

max_forks_repo_name: niclasmattsson/ClimateCalculator


text:
using CircuitConstants
using Documenter

DocMeta.setdocmeta!(CircuitConstants, :DocTestSetup, :(using CircuitConstants); recursive=true)

makedocs(;
    modules=[CircuitConstants],
    authors="Glen Hertz <glen.hertz@gmail.com> and contributors",
    repo="https://github.com/GlenHertz/CircuitConstants.jl/blob/{commit}{path}#{line}",
    sitename="CircuitConstants.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://glenhertz.github.io/CircuitConstants.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/GlenHertz/CircuitConstants.jl",
    devbranch="main",
)

ext: jl

max_stars_repo_name: GlenHertz/CircuitConstants

max_issues_repo_name: GlenHertz/CircuitConstants

max_forks_repo_name: GlenHertz/CircuitConstants


text:
module App

import HTTP, Sockets

include("Config.jl")
include("State.jl")
include("Router.jl")

function julia_main()::Cint
    try
        real_main()
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

function real_main()
    config = Config.getconfig()
    state = State.SharedState(config)
    router = Router.HttpRouter(state)

    host = "0.0.0.0"
    port = 8080
    server = Sockets.listen(Sockets.InetAddr(parse(Sockets.IPAddr, host), port))
    p = @async HTTP.serve(router, host, port; server = server)
    @info "Service is running..."

    # close server which will stop HTTP.serve
    Base.exit_on_sigint(false)
    try
        wait(p)
    catch e
        if isa(e, InterruptException)
            @warn "Caught SIGINT"
            close(server)
            wait(p)
        else
            @warn "Caught SIGTERM"
            kill(p) # SIGTERM
            rethrow()
        end
    end
    return
end

end # module

max_stars_repo_name: DevOpsDemoTF/DevOpsDemo-template-Julia

max_issues_repo_name: DevOpsDemoTF/DevOpsDemo-template-Julia

max_forks_repo_name: DevOpsDemoTF/DevOpsDemo-template-Julia


text:
# Add PATH (adapt!)
push!(LOAD_PATH, "/Volumes/cfranken/code/gitHub/LSM-SPAM/src/Leaf/");
push!(LOAD_PATH, "/Volumes/cfranken/code/gitHub/LSM-SPAM/src/Utils/");

using Leaf
using MathTools
using LeafPhotosynthesis
using BenchmarkTools
using Plots

pyplot()

# Create a leaf structure
l = leaf_params{Float32}();
l2 = leaf_params{Float32}();
# Create a Flux structure
f = LeafPhotosynthesis.fluxes{Float32}();

# initialize some reasonable values
f.je = 100;
f.gbc = 100;
f.gbv = 100;
f.ceair=1500;
f.eair = 1500;
f.APAR = 500;

l.Kn = 2.44
LeafPhotosynthesis.LeafPhotosynthesis(f,l,298.0,f.APAR);

# Test speed
@benchmark LeafPhotosynthesis.LeafPhotosynthesis(f,l,rand(Float32,1)[1]+300, 1000)

apar = Float32[]
gpp  = Float32[]
apar_ = Float32[]
gpp_  = Float32[]
f.cair = 400
l.vcmax25 = 90
l.jmax25 = 90*1.9
l2.vcmax25 = 50
l2.jmax25 = 50*1.9
for f.APAR in 0:10:1000
    LeafPhotosynthesis.LeafPhotosynthesis(f,l,298.0,f.APAR);
    push!(apar, f.APAR)
    push!(gpp, f.an)
    LeafPhotosynthesis.LeafPhotosynthesis(f,l2,298.0,f.APAR);
    push!(apar_, f.APAR)
    push!(gpp_, f.an)
    #println(f.APAR, " ", f.ag, " ", f.an)
end

plot(apar, gpp, lw=3, label="Vcmax25=90", ylabel="An", xlabel="PAR")
plot!(apar, gpp_, lw=3, label="Vcmax25=50")

# Test T dependence
T_ = Float32[]
ag_T_  = Float32[]
an_T_  = Float32[]
vcmax_T_  = Float32[]
electron_co2_T_  = Float32[]
gamma_T_  = Float32[]
f.APAR = 500
for T in 260:1:320
    LeafPhotosynthesis.LeafPhotosynthesis(f,l,T,f.APAR);
    push!(T_, T)
    push!(ag_T_, f.ag)
    push!(an_T_, f.an)
    push!(vcmax_T_, l.vcmax)
    push!(electron_co2_T_, l.CO2_per_electron)
    push!(gamma_T_ , l.Γstar)
end

plot(T_.-273.15, ag_T_, lw=3, label="Ag; Vcmax25=90, PAR=500", xlabel="T (degC)")
plot!(T_.-273.15, an_T_, lw=3, label="An; Vcmax25=90, PAR=500")
plot!(T_.-273.15, vcmax_T_/10, lw=3, label="Vcmax/10; Vcmax25=90, PAR=500")
plot!(T_.-273.15, electron_co2_T_*100, lw=3, label="CO2 per electron*100; Vcmax25=90, PAR=500")
plot!(T_.-273.15, gamma_T_/10, lw=3, label="CO2 Compensation Point/10; Vcmax25=90, PAR=500")

tmax = 5000 # Time in seconds here
l = leaf_params{Float32}()
l.gstyp=0;
f = LeafPhotosynthesis.fluxes{Float32}()
l.dynamic_state = true
#l.gstyp = 1
f.cair = 400
l.vcmax25 = 90
l.jmax25 = 90*1.9
uu = zeros(tmax,10)
uu[1,2] = 0.07 # Start gs
uu[1,1] = 0.5  # Start Kn
f.APAR = 500
l.rdleaf = 0.0

for c = 2:1:tmax
    l.gs = uu[c-1,2]
    l.Kn = uu[c-1,1]

    uu[c-1,3]=l.ϕs
    uu[c-1,4]=l.Ci
    uu[c-1,5]=f.φ
    uu[c-1,6]=l.Kp
    uu[c-1,8]=f.an
    uu[c-1,9]=l.Fm
    uu[c-1,7]=f.APAR
    #l.Kp = 4
    LeafPhotosynthesis.LeafPhotosynthesis(f,l,298.0,f.APAR)
    uu[c,1] = uu[c-1,1]+(l.Kn_ss-l.Kn)/0.2*1/60
    uu[c,2] = uu[c-1,2]+(l.gs_ss-l.gs)/15*1/60
    #println(l.Kn_ss, " ",  l.Ci, " ", f.φ)
end
t = range(1,tmax,step=1)

@show f.φ
@show l.ϕs
@show l.CO2_per_electron
@show l.Ci
@show l.Γstar

# Compute steady state solution
l2.vcmax25 = 90
l2.jmax25 = 90*1.9
l2.gstyp=0
l2.dynamic_state=false
LeafPhotosynthesis.LeafPhotosynthesis(f,l2,298.0, f.APAR)

println(f.φ)
println(l2.ϕs)
println(l2.CO2_per_electron)
println(l2.Ci)
println(l2.Γstar)

#p1 = plot(t[2:end-1]/60,uu[2:end-1,1],ylabel="Kn",lw=3, label="Kn") # Make a line plot
#p2 = plot(t[2:end-1]/60,uu[2:end-1,6],ylabel="Kp",lw=3, label="Kp") # Make a scatter plot
p3 = plot(t[2:end-1]/60,uu[2:end-1,3],ylabel="Fs Yield",lw=3, label="Fs Yield")
plot!([t[2],t[end-1]]/60,[l2.ϕs,l2.ϕs ],label="Fs Yield Steady State" )
p4 = plot(t[2:end-1]/60,uu[2:end-1,4],ylabel="Ci",lw=3, label="Ci")
plot!([t[2],t[end-1]]/60,[l2.Ci,l2.Ci ],label="Ci Steady State" )
#p5 = plot(t[2:end-1]/60,uu[2:end-1,7],ylabel="PAR",lw=3, label="PAR")# Four histograms each with 10 points? Why not!
p6 = plot(t[2:end-1]/60,uu[2:end-1,8],ylabel="An",lw=3, label="An", xlabel="Time (min)")# Four histograms each with 10 points? Why not!
plot!([t[2],t[end-1]]/60,[f.an,f.an ],label="An Steady State" )
p7 = plot(t[2:end-1]/60,uu[2:end-1,2],ylabel="gs",lw=3, label="gs", xlabel="Time (min)")
plot!([t[2],t[end-1]]/60,[l.gs,l.gs ],label="gs Steady State" )
#p8 = plot(t[2:end-1]/60,uu[2:end-1,5],ylabel="PSII yield",lw=3, label="PSII yield")

plot(p3,p4,p6,p7,layout=(2,2))

using DelimitedFiles
PAR  = readdlm("/Users/cfranken/PAR.dat");
length(PAR)

tmax = length(PAR) # Time in seconds here
t = range(1,tmax,step=1)
l = leaf_params{Float32}()
f = LeafPhotosynthesis.fluxes{Float32}()
l.dynamic_state = true
l.gstyp = 0
l.g1=4
f.cair = 400
l.vcmax25 = 90
l.jmax25 = 90*1.9
uu = zeros(tmax,10)
uu[1,2] = 0.07
uu[1,1] = 0.5
f.APAR = 50
l.rdleaf = 0.0

#apar[1:500]=0
#apar[501:100]=0
for c = 2:1:tmax
    #println(apar[c])
    f.APAR = PAR[c]*0.5
    f.APAR = max(10,f.APAR)

       l.gs = uu[c-1,2]
       l.Kn = uu[c-1,1]

       uu[c-1,3]=l.ϕs
       uu[c-1,4]=l.Ci
       uu[c-1,5]=f.φ
       uu[c-1,6]=l.Kp
       uu[c-1,8]=f.an
       uu[c-1,9]=l.Fm′
       uu[c-1,7]=f.APAR
       #l.Kp = 4
       LeafPhotosynthesis.LeafPhotosynthesis(f,l,298.0,f.APAR)
    if (l.Kn_ss-l.Kn) > 0
        tau_k = 1
    else
        tau_k = 1
    end
       uu[c,1] = uu[c-1,1]+(l.Kn_ss-l.Kn)/tau_k*1/60
       uu[c,2] = uu[c-1,2]+(l.gs_ss-l.gs)/20*1/60
    #println(l.Kn_ss, " ",  l.Ci, " ", f.φ, " ", l.gs_ss, " ", f.APAR)
end

step = 10
xl = (60*750,60*1400)
ranger = xl[1]:step:xl[2]

p1 = plot(t[ranger]/60,uu[ranger,1],ylabel="Kn",lw=3, title="Kn", legend=false) # Make a line plot
p2 = plot(t[ranger]/60,uu[ranger,6],ylabel="Kp",lw=3, title="Kp", legend=false) # Make a scatter plot
p3 = plot(t[ranger]/60,uu[ranger,3],ylabel="Fs Yield",lw=3, title="Fs Yield", legend=false)
p4 = plot(t[ranger]/60,uu[ranger,4],ylabel="Ci",lw=3, title="Ci", legend=false,  ylim=(220,330))
p5 = plot(t[ranger]/60,uu[ranger,7],ylabel="PAR",lw=3, title="PAR", legend=false)# Four histograms each with 10 points? Why not!
p6 = plot(t[ranger]/60,uu[ranger,8],ylabel="An",lw=3, title="An", legend=false)# Four histograms each with 10 points? Why not!
p7 = plot(t[ranger]/60,uu[ranger,2],ylabel="gs'",lw=3, title="gs", legend=false)
p8 = plot(t[ranger]/60,uu[ranger,5],ylabel="PSII yield",lw=3, title="PSII yield", legend=false)

plot(p1,p2,p3,p4,p5,p6,p7,p8,layout=(4,2))

tmax = length(PAR) # Time in seconds here
t = range(1,tmax,step=1)
l = leaf_params{Float32}()
f = LeafPhotosynthesis.fluxes{Float32}()
l.dynamic_state = false
l.g1=4
l.gstyp = 0
f.cair = 400
l.vcmax25 = 90
l.jmax25 = 90*1.9
uus = zeros(tmax,10)
uus[1,2] = 0.07
uus[1,1] = 0.5
f.APAR = 50
l.rdleaf = 0.0

#apar[1:500]=0
#apar[501:100]=0
for c = 2:1:tmax
    #println(apar[c])
    f.APAR = PAR[c]*0.5
    f.APAR = max(10,f.APAR)
    LeafPhotosynthesis.LeafPhotosynthesis(f,l,298.0,f.APAR)
    uus[c,2]=l.gs
    uus[c,1]=l.Kn_ss

    uus[c,3]=l.ϕs
    uus[c,4]=l.Ci
    uus[c,5]=f.φ
    uus[c,6]=l.Kp
    uus[c,8]=f.an
    uus[c,9]=l.Fm′
    uus[c,7]=f.APAR

end

step = 30
xl = (60*750,60*1400)
ranger = xl[1]:step:xl[2]

p1 = plot(t[ranger]/60,uus[ranger,1],ylabel="Kn",lw=3, title="Kn", legend=false) # Make a line plot
p2 = plot(t[ranger]/60,uus[ranger,6],ylabel="Kp",lw=3, title="Kp", legend=false) # Make a scatter plot
p3 = plot(t[ranger]/60,uus[ranger,3],ylabel="Fs Yield",lw=3, title="Fs Yield", legend=false)
p4 = plot(t[ranger]/60,uus[ranger,4],ylabel="Ci",lw=3, title="Ci", legend=false,  ylim=(220,330))
p5 = plot(t[ranger]/60,uus[ranger,7],ylabel="PAR",lw=3, title="PAR", legend=false)# Four histograms each with 10 points? Why not!
p6 = plot(t[ranger]/60,uus[ranger,8],ylabel="An",lw=3, title="An", legend=false)# Four histograms each with 10 points? Why not!
p7 = plot(t[ranger]/60,uus[ranger,2],ylabel="gs'",lw=3, title="gs", legend=false)
p8 = plot(t[ranger]/60,uus[ranger,5],ylabel="PSII yield",lw=3, title="PSII yield", legend=false)

plot(p1,p2,p3,p4,p5,p6,p7,p8,layout=(4,2))

step = 30
xl = (60*750,60*1400)
ranger = xl[1]:step:xl[2]

p1 = plot(t[ranger]/60,uus[ranger,1],ylabel="Kn",lw=3, title="Kn", legend=false) # Make a line plot
p1 = plot!(t[ranger]/60,uu[ranger,1],ylabel="Kn",lw=3, title="Kn", legend=false) # Make a line plot

p2 = plot(t[ranger]/60,uus[ranger,6],ylabel="Kp",lw=3, title="Kp", legend=false) # Make a scatter plot
p2 = plot!(t[ranger]/60,uu[ranger,6],ylabel="Kp",lw=3, title="Kp", legend=false) # Make a scatter plot

p3 = plot(t[ranger]/60,uus[ranger,3],ylabel="Fs Yield",lw=3, title="Fs Yield", legend=false)
p3 = plot!(t[ranger]/60,uu[ranger,3],ylabel="Fs Yield",lw=3, title="Fs Yield", legend=false)

p4 = plot(t[ranger]/60,uus[ranger,4],ylabel="Ci",lw=3, title="Ci", legend=false,  ylim=(220,330))
p4 = plot!(t[ranger]/60,uu[ranger,4],ylabel="Ci",lw=3, title="Ci", legend=false,  ylim=(220,330))
#
p5 = plot(t[ranger]/60,uus[ranger,7],ylabel="PAR",lw=3, title="PAR", legend=false)# Four histograms each with 10 points? Why not!
p5 = plot!(t[ranger]/60,uu[ranger,7],ylabel="PAR",lw=3, title="PAR", legend=false)# Four histograms each with 10 points? Why not!
#
p6 = plot(t[ranger]/60,uus[ranger,8],ylabel="An",lw=3, title="An", legend=false)# Four histograms each with 10 points? Why not!
p6 = plot!(t[ranger]/60,uu[ranger,8],ylabel="An",lw=3, title="An", legend=false)# Four histograms each with 10 points? Why not!
#
p7 = plot(t[ranger]/60,uus[ranger,2],ylabel="gs'",lw=3, title="gs", legend=false)
p7 = plot!(t[ranger]/60,uu[ranger,2],ylabel="gs'",lw=3, title="gs", legend=false)
#
p8 = plot(t[ranger]/60,uus[ranger,5],ylabel="PSII yield",lw=3, title="PSII yield", legend=false)
p8 = plot!(t[ranger]/60,uu[ranger,5],ylabel="PSII yield",lw=3, title="PSII yield", legend=false)

#p2 = plot(t[ranger]/60,uu[ranger,6],ylabel="Kp",lw=3, title="Kp", legend=false) # Make a scatter plot
#p3 = plot(t[ranger]/60,uu[ranger,3],ylabel="Fs Yield",lw=3, title="Fs Yield", legend=false)
#p4 = plot(t[ranger]/60,uu[ranger,4],ylabel="Ci",lw=3, title="Ci", legend=false,  ylim=(220,330))
#p5 = plot(t[ranger]/60,uu[ranger,7],ylabel="PAR",lw=3, title="PAR", legend=false)# Four histograms each with 10 points? Why not!
#p6 = plot(t[ranger]/60,uu[ranger,8],ylabel="An",lw=3, title="An", legend=false)# Four histograms each with 10 points? Why not!
#p7 = plot(t[ranger]/60,uu[ranger,2],ylabel="gs'",lw=3, title="gs", legend=false)
#p8 = plot(t[ranger]/60,uu[ranger,5],ylabel="PSII yield",lw=3, title="PSII yield", legend=false)

plot(p1,p2,p3,p4,p5,p6,p7,p8,layout=(4,2))

tmax = 6100 # Time in seconds here
t = range(1,tmax,step=1)
l = leaf_params{Float32}()
f = LeafPhotosynthesis.fluxes{Float32}()
l.dynamic_state = true
#l.gstyp = 1
f.cair = 400
l.vcmax25 = 90
l.jmax25 = 90*1.9
uu = zeros(tmax,10)
uu[1,2] = 0.07
uu[1,1] = 0.5
f.APAR = 50
l.rdleaf = 0.0
apar = zeros(tmax)
apar[:].=100
for i in 1:8
    #println(i," ", (i-1)*500+1, " ", (i)*500, " ", i*100 )
    apar[(i-1)*500+1:(i)*500].=i*150
end
for i in 9:12
    println(i," ", (i-1)*500+1, " ", (i)*500, " ", 1500-(i-9)*500 )
    apar[(i-1)*500+1:(i)*500].=1500-(i-9)*700
end
#apar[1:500]=0
#apar[501:100]=0
for c = 2:1:tmax
    #println(apar[c])
    f.APAR = apar[c]

       l.gs = uu[c-1,2]
       l.Kn = uu[c-1,1]

       uu[c-1,3]=l.ϕs
       uu[c-1,4]=l.Ci
       uu[c-1,5]=f.φ
       uu[c-1,6]=l.Kp
       uu[c-1,8]=f.an
       uu[c-1,9]=l.Fm′
       uu[c-1,7]=f.APAR
       #l.Kp = 4
       LeafPhotosynthesis.LeafPhotosynthesis(f,l,298.0,f.APAR)
    if (l.Kn_ss-l.Kn) > 0
        tau_k = 0.1
    else
        tau_k = 1
    end
       uu[c,1] = uu[c-1,1]+(l.Kn_ss-l.Kn)/tau_k*1/60
       uu[c,2] = uu[c-1,2]+(l.gs_ss-l.gs)/20*1/60
       #println(l.Kn_ss, " ",  l.Ci, " ", f.φ)
end

p1 = plot(t[2:end-1]/60,uu[2:end-1,1],ylabel="Kn",lw=3, title="Kn", legend=false) # Make a line plot
p2 = plot(t[2:end-1]/60,uu[2:end-1,6],ylabel="Kp",lw=3, title="Kp", legend=false) # Make a scatter plot
p3 = plot(t[2:end-1]/60,uu[2:end-1,3],ylabel="Fs Yield",lw=3, title="Fs Yield", legend=false)
p4 = plot(t[2:end-1]/60,uu[2:end-1,4],ylabel="Ci",lw=3, title="Ci", legend=false)
p5 = plot(t[2:end-1]/60,uu[2:end-1,7],ylabel="PAR",lw=3, title="PAR", legend=false)# Four histograms each with 10 points? Why not!
p6 = plot(t[2:end-1]/60,uu[2:end-1,8],ylabel="An",lw=3, title="An", legend=false)# Four histograms each with 10 points? Why not!
p7 = plot(t[2:end-1]/60,uu[2:end-1,9],ylabel="Fm'",lw=3, title="Fm'", legend=false)
p8 = plot(t[2:end-1]/60,uu[2:end-1,5],ylabel="PSII yield",lw=3, title="PSII yield", legend=false)

plot(p1,p2,p3,p4,p5,p6,p7,p8,layout=(4,2))

using Statistics
PAR2 = PAR;
# 10 minutes
wind = 60*10

for i in 1:length(PAR2)

    mini = max(1,i-wind)
    maxi = min(length(PAR2),i+wind)
    PAR2[i] = mean(PAR[mini:maxi])
end

tmax = length(PAR) # Time in seconds here
t = range(1,tmax,step=1)
l = leaf_params{Float32}()
f = LeafPhotosynthesis.fluxes{Float32}()
l.dynamic_state = true
l.gstyp = 0
l.g1=3.5
f.cair = 400
l.vcmax25 = 90
l.jmax25 = 90*1.9
uu = zeros(tmax,10)
uu[1,2] = 0.07
uu[1,1] = 0.5
f.APAR = 50
l.rdleaf = 0.0

#apar[1:500]=0
#apar[501:100]=0
for c = 2:1:tmax
    #println(apar[c])
    f.APAR = PAR2[c]*0.5
    f.APAR = max(10,f.APAR)

       l.gs = uu[c-1,2]
       l.Kn = uu[c-1,1]

       uu[c-1,3]=l.ϕs
       uu[c-1,4]=l.Ci
       uu[c-1,5]=f.φ
       uu[c-1,6]=l.Kp
       uu[c-1,8]=f.an
       uu[c-1,9]=l.Fm′
       uu[c-1,7]=f.APAR
       #l.Kp = 4
       LeafPhotosynthesis.LeafPhotosynthesis(f,l,298.0,f.APAR)
    if (l.Kn_ss-l.Kn) > 0
        tau_k = 1
    else
        tau_k = 1
    end
       uu[c,1] = uu[c-1,1]+(l.Kn_ss-l.Kn)/tau_k*1/60
       uu[c,2] = uu[c-1,2]+(l.gs_ss-l.gs)/20*1/60
    #println(l.Kn_ss, " ",  l.Ci, " ", f.φ, " ", l.gs_ss, " ", f.APAR)
end

tmax = length(PAR) # Time in seconds here
t = range(1,tmax,step=1)
l = leaf_params{Float32}()
f = LeafPhotosynthesis.fluxes{Float32}()
l.dynamic_state = false
l.g1=3.5
l.gstyp = 0
f.cair = 400
l.vcmax25 = 90
l.jmax25 = 90*1.9
uus = zeros(tmax,10)
uus[1,2] = 0.07
uus[1,1] = 0.5
f.APAR = 50
l.rdleaf = 0.0

#apar[1:500]=0
#apar[501:100]=0
for c = 2:1:tmax
    #println(apar[c])
    f.APAR = PAR2[c]*0.5
    f.APAR = max(10,f.APAR)
    LeafPhotosynthesis.LeafPhotosynthesis(f,l,298.0,f.APAR )
    uus[c,2]=l.gs
    uus[c,1]=l.Kn_ss

    uus[c,3]=l.ϕs
    uus[c,4]=l.Ci
    uus[c,5]=f.φ
    uus[c,6]=l.Kp
    uus[c,8]=f.an
    uus[c,9]=l.Fm′
    uus[c,7]=f.APAR

end

step = 30
xl = (60*750,60*1400)
ranger = xl[1]:step:xl[2]

p1 = plot(t[ranger]/60,uus[ranger,1],ylabel="Kn",lw=3, title="Kn", legend=false) # Make a line plot
p1 = plot!(t[ranger]/60,uu[ranger,1],ylabel="Kn",lw=3, title="Kn", legend=false) # Make a line plot

p2 = plot(t[ranger]/60,uus[ranger,6],ylabel="Kp",lw=3, title="Kp", legend=false) # Make a scatter plot
p2 = plot!(t[ranger]/60,uu[ranger,6],ylabel="Kp",lw=3, title="Kp", legend=false) # Make a scatter plot

p3 = plot(t[ranger]/60,uus[ranger,3],ylabel="Fs Yield",lw=3, title="Fs Yield", legend=false)
p3 = plot!(t[ranger]/60,uu[ranger,3],ylabel="Fs Yield",lw=3, title="Fs Yield", legend=false)

p4 = plot(t[ranger]/60,uus[ranger,4],ylabel="Ci",lw=3, title="Ci", legend=false,  ylim=(220,330))
p4 = plot!(t[ranger]/60,uu[ranger,4],ylabel="Ci",lw=3, title="Ci", legend=false,  ylim=(220,330))
#
p5 = plot(t[ranger]/60,uus[ranger,7],ylabel="PAR",lw=3, title="PAR", legend=false)# Four histograms each with 10 points? Why not!
p5 = plot!(t[ranger]/60,uu[ranger,7],ylabel="PAR",lw=3, title="PAR", legend=false)# Four histograms each with 10 points? Why not!
#
p6 = plot(t[ranger]/60,uus[ranger,8],ylabel="An",lw=3, title="An", legend=false)# Four histograms each with 10 points? Why not!
p6 = plot!(t[ranger]/60,uu[ranger,8],ylabel="An",lw=3, title="An", legend=false)# Four histograms each with 10 points? Why not!
#
p7 = plot(t[ranger]/60,uus[ranger,2],ylabel="gs'",lw=3, title="gs", legend=false)
p7 = plot!(t[ranger]/60,uu[ranger,2],ylabel="gs'",lw=3, title="gs", legend=false)
#
p8 = plot(t[ranger]/60,uus[ranger,5],ylabel="PSII yield",lw=3, title="PSII yield", legend=false)
p8 = plot!(t[ranger]/60,uu[ranger,5],ylabel="PSII yield",lw=3, title="PSII yield", legend=false)

#p2 = plot(t[ranger]/60,uu[ranger,6],ylabel="Kp",lw=3, title="Kp", legend=false) # Make a scatter plot
#p3 = plot(t[ranger]/60,uu[ranger,3],ylabel="Fs Yield",lw=3, title="Fs Yield", legend=false)
#p4 = plot(t[ranger]/60,uu[ranger,4],ylabel="Ci",lw=3, title="Ci", legend=false,  ylim=(220,330))
#p5 = plot(t[ranger]/60,uu[ranger,7],ylabel="PAR",lw=3, title="PAR", legend=false)# Four histograms each with 10 points? Why not!
#p6 = plot(t[ranger]/60,uu[ranger,8],ylabel="An",lw=3, title="An", legend=false)# Four histograms each with 10 points? Why not!
#p7 = plot(t[ranger]/60,uu[ranger,2],ylabel="gs'",lw=3, title="gs", legend=false)
#p8 = plot(t[ranger]/60,uu[ranger,5],ylabel="PSII yield",lw=3, title="PSII yield", legend=false)

plot(p1,p2,p3,p4,p5,p6,p7,p8,layout=(4,2))

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
# Test data are generated using abelfunctions and Sage 8.6

ϵ = 1e-3

@testset "Genus 1" begin
    z = [0.85746943+0.30772689im];
    τ = [0.63443932+0.49919024im];
    Y = imag(τ);
    x = real(z);
    y = convert(Array{Float64}, imag(z));
    Yinv = inv.(Y);
    y0 = Yinv.*y;
    R = RiemannMatrix(τ, siegel=false);
    @test Theta.oscillatory_part(R, x, y0, [0]) ≈ -0.257+0.197im atol=ϵ
    @test theta(z, R) ≈ -0.466+0.358im atol=ϵ
    @test theta(z, R, derivs=[[1]]) ≈ 1.724+9.922im atol=ϵ
    @test theta(z, R, derivs=[[1],[1]]) ≈ 67.973-10.642im atol=ϵ
    @test theta(z, R, char=[[1],[0]]) ≈ -1.770-1.225im atol=ϵ
    @test theta(z, R, char=[[1],[1]]) ≈ -1.669+0.298im atol=ϵ
    @test theta(z, R, char=[[1],[0]], derivs=[[1]]) ≈ -2.368+6.931im atol=ϵ
    @test theta(z, R, char=[[1],[1]], derivs=[[1]]) ≈ -0.488+6.871im atol=ϵ
    @test symplectic_transform(siegel_transform(τ)[1], τ) ≈ siegel_transform(τ)[2]
end

@testset "Genus 2" begin
    z = [0.81149300+0.27027128im; 0.77132834+0.26619567im];
    τ = [0.71106237+1.20021283im 0.57281731+0.89762698im; 0.57281731+0.89762698im 0.22079146+0.68617488im];
    Y = imag(τ);
    x = real(z);
    y = convert(Array{Float64}, imag(z));
    Yinv = inv(Y);
    y0 = Yinv*y
    R = RiemannMatrix(τ, siegel=false);
    @test Theta.oscillatory_part(R, x, y0, [0,0]) ≈ 0.5894+0.3593im atol=ϵ
    @test theta(z, R) ≈ 1.700489+1.03657im atol=ϵ
    @test theta(z, R, derivs=[[1,0]]) ≈ -13.73322+34.5327im atol=ϵ
    @test theta(z, R, derivs=[[0,1]]) ≈ 14.71768-51.3656im atol=ϵ
    @test theta(z, R, derivs=[[1,0], [1,0]]) ≈ -1859.654+898.206im atol=ϵ
    @test theta(z, R, derivs=[[1,0], [0,1]]) ≈ 2533.2922-1397.9940im atol=ϵ
    @test theta(z, R, char=[[1,0],[1,1]]) ≈ -5.91176+4.90798im atol=ϵ
    @test theta(z, R, char=[[1,0],[0,1]]) ≈ -0.68197+0.13723im atol=ϵ
    @test theta(z, R, char=[[1,0],[1,1]], derivs=[[1,0]]) ≈ -33.55188-104.14870im atol=ϵ
    @test theta(z, R, char=[[1,0],[0,1]], derivs=[[0,1]]) ≈ 20.42986+95.749276im atol=ϵ
    @test symplectic_transform(siegel_transform(τ)[1], τ) ≈ siegel_transform(τ)[2]
end

@testset "Genus 3" begin
    z = [0.76394595+0.5283701im; 0.03744967+0.0898654im; 0.65786408+0.77005086im];
    τ = [0.85625765+1.13705995im 0.14147041+0.43763364im 0.64991508+1.07829521im; 0.14147041+0.43763364im 0.28765738+0.60567227im 0.27742587+0.38069182im; 0.64991508+1.07829521im 0.27742587+0.38069182im 0.23093880+1.81289325im];
    Y = imag(τ);
    x = real(z);
    y = convert(Array{Float64}, imag(z));
    Yinv = inv(Y);
    y0 = Yinv*y
    R = RiemannMatrix(τ, siegel=false);
    @test Theta.oscillatory_part(R, x, y0, [0,0,0]) ≈ 0.11045-0.28985im atol=ϵ
    @test theta(z, R) ≈ 0.34318-0.900595im atol=ϵ
    @test theta(z, R, derivs=[[1,0,0]]) ≈ -11.647091+6.51591im atol=ϵ
    @test theta(z, R, derivs=[[0,1,0]]) ≈ 5.57557-4.17216im atol=ϵ
    @test theta(z, R, derivs=[[0,0,1]]) ≈ 3.68145+1.71192im atol=ϵ
    @test theta(z, R, derivs=[[1,0,0], [1,0,0]]) ≈ 7.8837+62.5345im atol=ϵ
    @test theta(z, R, derivs=[[1,0,0], [0,1,0]]) ≈ -9.2622-32.7841im atol=ϵ
    @test theta(z, R, derivs=[[0,0,1], [0,0,1]]) ≈ 9.8031-19.6301im atol=ϵ
    @test theta(z, R, char=[[0,0,1],[0,1,1]]) ≈ -1.345015+0.31343im atol=ϵ
    @test theta(z, R, char=[[1,1,0],[0,1,0]]) ≈ 1.039441-0.95682im atol=ϵ
    @test theta(z, R, char=[[0,0,1],[0,1,1]], derivs=[[0,1,0]]) ≈ 4.045043+3.521485im atol=ϵ
    @test symplectic_transform(siegel_transform(τ)[1], τ) ≈ siegel_transform(τ)[2]
end

@testset "Genus 4" begin
    z = [0.04134584+0.40910551im; 0.20972589+0.90269823im; 0.39996195+0.42432923im; 0.73063375+0.49945621im];
    τ = [0.95870734+0.73587725im 0.22092477+0.76863646im 0.53877459+0.87577267im 0.68177023+0.867436im; 0.22092477+0.76863646im 0.98812562+1.79674905im 0.54859032+1.10626215im 0.63310305+1.30158981im; 0.53877459+0.87577267im 0.54859032+1.10626215im 0.50173043+1.27729044im 0.49163557+1.33147334im; 0.68177023+0.867436im 0.63310305+1.30158981im 0.49163557+1.33147334im 0.35312207+1.60745975im];
    Y = imag(τ);
    x = real(z);
    y = convert(Array{Float64}, imag(z));
    Yinv = inv(Y);
    y0 = Yinv*y
    R = RiemannMatrix(τ, siegel=false);
    @test Theta.oscillatory_part(R, x, y0, [0,0,0,0]) ≈ -1.35525+1.05162im atol=ϵ
    @test theta(z, R) ≈ -8.24515+6.39791im atol=ϵ
    @test theta(z, R, derivs=[[1,0,0,0]]) ≈ 62.08697+26.90192im atol=ϵ
    @test theta(z, R, derivs=[[0,1,0,0]]) ≈ 21.74519+37.036989im atol=ϵ
    @test theta(z, R, derivs=[[0,0,1,0]]) ≈ -75.17050+14.28415im atol=ϵ
    @test theta(z, R, derivs=[[0,0,0,1]]) ≈ 32.02015-36.54835im atol=ϵ
    @test theta(z, R, derivs=[[1,0,0,0], [1,0,0,0]]) ≈ -133.4336-555.3991im atol=ϵ
    @test theta(z, R, derivs=[[1,0,0,0], [0,1,0,0]]) ≈ 4.08109-225.3261im atol=ϵ
    @test theta(z, R, derivs=[[0,1,0,0], [0,0,1,0]]) ≈ 132.01459+365.6160im atol=ϵ
    @test theta(z, R, derivs=[[0,0,0,1], [0,0,0,1]]) ≈ 89.8953-145.9657im atol=ϵ
    @test theta(z, R, char=[[0,1,0,1],[0,1,0,0]]) ≈ 15.51204+1.32669im atol=ϵ
    @test theta(z, R, char=[[0,1,0,0],[0,0,0,1]]) ≈ 1.852505-1.75216im atol=ϵ
    @test theta(z, R, char=[[0,1,0,1],[0,1,0,0]], derivs=[[0,0,1,0]]) ≈ 28.8916+90.992105im atol=ϵ
    @test symplectic_transform(siegel_transform(τ)[1], τ) ≈ siegel_transform(τ)[2]
end

@testset "Genus 5" begin
    z = [0.76931119+0.61774429im; 0.74693159+0.23398854im; 0.68457901+0.49648231im; 0.59332210+0.25122584im; 0.02494457+0.27172437im];
    τ = [0.17012929+3.36021399im 0.42793573+1.65263414im 0.55503983+2.1393021im 0.49182605+2.17345553im 0.93865305+1.81009273im; 0.42793573+1.65263414im 0.87402978+1.33203047im 0.53489381+0.68238397im 0.57601957+1.20509707im 0.50028888+1.11605622im; 0.55503983+2.1393021im 0.53489381+0.68238397im 0.21882206+2.28319153im 0.74421992+1.51352273im 0.41855281+0.87400937im; 0.49182605+2.17345553im 0.57601957+1.20509707im 0.74421992+1.51352273im 0.51070570+1.62305859im 0.52380319+1.15073893im; 0.93865305+1.81009273im 0.50028888+1.11605622im 0.41855281+0.87400937im 0.52380319+1.15073893im 0.79948032+1.38925759im];
    Y = imag(τ);
    x = real(z);
    y = convert(Array{Float64}, imag(z));
    Yinv = inv(Y);
    y0 = Yinv*y
    R = RiemannMatrix(τ, siegel=false);
    @test Theta.oscillatory_part(R, x, y0, [0,0,0,0,0]) ≈ -0.13999-0.20302im atol=ϵ
    @test theta(z, R) ≈ -0.43584-0.632066im atol=ϵ
    @test theta(z, R, derivs=[[1,0,0,0,0]]) ≈ -16.240575+7.856723im atol=ϵ
    @test theta(z, R, derivs=[[0,1,0,0,0]]) ≈ 14.447895+1.37822im atol=ϵ
    @test theta(z, R, derivs=[[0,0,1,0,0]]) ≈ 14.232547+0.23889im atol=ϵ
    @test theta(z, R, derivs=[[0,0,0,1,0]]) ≈ -5.509064-15.47362im atol=ϵ
    @test theta(z, R, derivs=[[0,0,0,0,1]]) ≈ 0.796003+6.15536im atol=ϵ
    @test theta(z, R, derivs=[[1,0,0,0,0], [1,0,0,0,0]]) ≈ 150.0687+120.7466im atol=ϵ
    @test theta(z, R, derivs=[[0,0,1,0,0], [0,0,0,1,0]]) ≈ -33.6935+132.9402im atol=ϵ
    @test theta(z, R, derivs=[[0,0,1,0,0], [0,1,0,0,0]]) ≈ -15.0943-138.2480im atol=ϵ
    @test theta(z, R, derivs=[[0,0,0,1,0], [0,0,0,1,0]]) ≈ 330.7400-145.8428im atol=ϵ
    @test theta(z, R, derivs=[[0,0,0,0,1], [0,1,0,0,0]]) ≈ 127.4712+47.0716im atol=ϵ
    @test theta(z, R, char=[[0,1,0,1,1],[0,1,1,0,0]]) ≈ 1.03805-0.46603im atol=ϵ
    @test theta(z, R, char=[[1,1,1,0,1],[0,1,1,1,0]]) ≈ -1.41336+4.24756im atol=ϵ
    @test theta(z, R, char=[[1,1,1,0,1],[0,1,1,1,0]], derivs=[[1,0,0,0,0]]) ≈ -1.49522 + 5.56114im atol=ϵ
    @test symplectic_transform(siegel_transform(τ)[1], τ) ≈ siegel_transform(τ)[2]
end

@testset "Azygetic" begin
    @test check_azygetic([[[1,0,1,0], [1,0,1,0]], [[0,0,0,1], [1,0,0,0]], [[0,0,1,1], [1,0,1,1]]]) == true
end


text:
#!/usr/bin/env julia

using Test
using QuadOsc
using QuadGK
using SpecialFunctions

@testset "SeriesAccelerations Cohen-Villegas-Zagier" begin
	# Check first few coefficients as given in Cohen et al 2000
	ak = [1]
	A = QuadOsc.accel_cohen_villegas_zagier(ak)
	@test A == 2/3

	ak = [1, 0]
	A = QuadOsc.accel_cohen_villegas_zagier(ak)
	@test A ≈ 16/17
	ak = [0, -1]
	A = QuadOsc.accel_cohen_villegas_zagier(ak)
	@test A ≈ -8/17

	ak = [1, 0, 0]
	A = QuadOsc.accel_cohen_villegas_zagier(ak)
	@test A ≈ 98/99
	ak = [0, -1, 0]
	A = QuadOsc.accel_cohen_villegas_zagier(ak)
	@test A ≈ -80/99
	ak = [0, 0, 1]
	A = QuadOsc.accel_cohen_villegas_zagier(ak)
	@test A ≈ 32/99

	ak = [1, 0, 0, 0]
	A = QuadOsc.accel_cohen_villegas_zagier(ak)
	@test A ≈ 576/577
	ak = [0, -1, 0, 0]
	A = QuadOsc.accel_cohen_villegas_zagier(ak)
	@test A ≈ -544/577
	ak = [0, 0, 1, 0]
	A = QuadOsc.accel_cohen_villegas_zagier(ak)
	@test A ≈ 384/577
	ak = [0, 0, 0, -1]
	A = QuadOsc.accel_cohen_villegas_zagier(ak)
	@test A ≈ -128/577

	# Ex 1 in Cohen et al 2000
	# Note: we actually disagree with the answer given in the paper. What's
	# going on? Since all other tests pass, we probably solving a different
	# problem than them.
	fn(n) = (-1)^(n+1) * loggamma(1 + 1 / (2*n - 1))
	println("n an\tsum(ak)\tA")
	for n=1:20
		ak = @. fn(1:n)
		lnA = QuadOsc.accel_cohen_villegas_zagier(ak)
		A = exp(lnA)
		Adirect = exp(sum(ak))
		println("$n $(ak[end])\t$Adirect\t$A")
	end

	# Ex 3 in Cohen et al 2000
	riemann_zeta(s; n=20) = begin
		k = 0:n-1
		ak = @. (-1)^k / (k + 1)^s / (1 - 2^(1 - s))
		return QuadOsc.accel_cohen_villegas_zagier(ak)
	end
	@test riemann_zeta(1/2) ≈ -1.4603545088 rtol=0 atol=1e-10
	@test riemann_zeta(-1+im) ≈ 0.0168761517 - 0.1141564804*im rtol=0 atol=1e-10
end

@testset "SeriesAccelerations Wynn's epsilon" begin
	# Ex 1 in Cohen et al 2000
	# Note: we actually disagree with the answer given in the paper. What's
	# going on? Since all other tests pass, we probably solving a different
	# problem than them.
	fn(n) = (-1)^(n+1) * loggamma(1 + 1 / (2*n - 1))
	println("n an\tsum(ak)\tA")
	for n=1:20
		ak = @. fn(1:n)
		lnA = QuadOsc.accel_wynn_eps(ak)
		A = exp(lnA)
		Adirect = exp(sum(ak))
		println("$n $(ak[end])\t$Adirect\t$A")
	end

	# Ex 3 in Cohen et al 2000
	riemann_zeta(s; n=20) = begin
		k = 0:n-1
		ak = @. (-1)^k / (k + 1)^s / (1 - 2^(1 - s))
		return QuadOsc.accel_wynn_eps(ak)
	end
	@test riemann_zeta(1/2) ≈ -1.4603545088 rtol=0 atol=1e-10
	@test riemann_zeta(-1+im) ≈ 0.0168761517 - 0.1141564804*im rtol=0 atol=1e-10
end

function test_quadosc()
	integrand(x) = x == 0.0 ? 0.0 : x * besselj0(x) / (1 + x^2)
	nterm = 15
	a = 0.0
	b = 0.0
	sum = 0.0
	ak = Float64[]
	ans = 0.0
	correct = 0.4210244382
	println("N\tSum(direct)\tSum(Cohen+)\tSum(Wynn)")
	for n in 0:nterm
		b += pi
		s, E = quadgk(integrand, a, b)
		push!(ak, s)
		a = b
		sum += s
		ans = QuadOsc.accel_cohen_villegas_zagier(ak)
		ans2 = QuadOsc.accel_wynn_eps(ak)
		println("$n\t$sum\t$ans\t$ans2")
	end
	directans, E = quadosc(integrand, 0, Inf, n -> n * pi)
	println("ans          = $ans")
	println("directans    = $directans +- $E")
	println("correct + dI = $(correct) + $(directans - correct)")
	@test abs(ans - correct) < 1e-10
	@test abs(ans - directans) < 1e-10
	@test 0 <= E <= 1e-10
end

function test2()
	thisans, E = quadosc(x -> sphericalbesselj(0, x), 0, Inf, n -> n * pi)
	println("int sphbes0: $thisans +- $E =? $(pi/2)")
	println("I          : $(pi/2) + $(thisans - pi/2)")
	@test abs(thisans - pi/2) < E
	@test 0 <= E <= 1e-10
end

function test_int_Jν()
	for ℓ=0:100
		Jl_zeros(n,ℓ) = (n + ℓ/2) * π  # just an estimate
		I, E = quadosc(t -> besselj(ℓ+0.5, t), 0, Inf, n -> Jl_zeros(n,ℓ))
		@show ℓ,I,E
		@test I ≈ 1
		@test abs(I - 1) <= E
		@test 0 <= E <= 2*sqrt(eps(I))
	end
end

function typetest()
	integrand(x::Float64) = (trunc(Int, x) % 2 == 1) ? -1/x : 1/x
	@time quadosc(integrand, 1, Inf, n -> n)
	@time quadosc(integrand, 1, Inf, n -> n)
end

@testset "QuadOsc" begin
	test_quadosc()
	test2()
	test_int_Jν()
	typetest()
end


text:
# TODO: limit types for arguments of f
struct Predicate
  f::Function
end

struct Mapping
  f::Function
end

(p::Predicate)(x,y) = p.f(x,y)
(m::Mapping)(x) = m.f(x)

# ====================  Auxiliary functions ====================

"""
    check_for_equivalence_relation(f::Predicate, dom)

Given a Predicate `f` over a domain `dom`, return `true` if `f`
is a equivalence relation, `false` otherwise.
One can create a symmetric equivalency relation with [`make_symmetric(f)`](@ref)
TODO: check transitivity.

# Examples
```
julia> check_for_equivalence_relation(Predicate((x,y) -> x == -y), -5:10)
true
julia> check_for_equivalence_relation(Predicate((x,y) -> x == y+2), -1:10)
false
```
"""
function check_for_equivalence_relation(f::Predicate, dom)
  for (i,el1) in enumerate(dom)
    for el2 in dom
      !(f(el1,el2) == f(el2,el1)) && return false
    end
  end
  #TODO: check transitivity
  return true
end


text:
struct InterpolatingPolynomial{N,B,T}
    coeffs::Matrix{T}
    basis::B
    function InterpolatingPolynomial(
        coeffs::Matrix{T},
        basis::B,
    ) where {B<:AbstractBasis{dim}} where {dim,T}

        N, NF = size(coeffs)
        @assert number_of_basis_functions(basis) == NF
        new{N,B,T}(coeffs, basis)

    end
end

function InterpolatingPolynomial(
    T::Type{<:Number},
    N::Z,
    basis::B,
) where {Z<:Integer} where {B<:AbstractBasis{dim,NF}} where {dim,NF}

    coeffs = zeros(T, N, NF)
    return InterpolatingPolynomial(coeffs, basis)

end

function InterpolatingPolynomial(
    N::Z,
    basis::B,
) where {Z<:Integer} where {B<:AbstractBasis}
    T = typeof(0.0)
    return InterpolatingPolynomial(T, N, basis)
end

function InterpolatingPolynomial(
    N::Z,
    dim::Z,
    order::Z;
    start = -1.0,
    stop = 1.0,
) where {Z<:Integer}
    T = typeof(0.0)
    basis = LagrangeTensorProductBasis(dim, order, start = start, stop = stop)
    return InterpolatingPolynomial(T, N, basis)
end

function Base.show(io::IO, poly::InterpolatingPolynomial{N,B,T}) where {N,B,T}
    p = order(poly)
    dim = dimension(poly)
    NF = number_of_basis_functions(poly.basis)
    basistype = typeof(poly.basis)
    str =
        "InterpolatingPolynomial\n\tDimension  : $dim\n\t" *
        "Basis Type : $basistype\n\t" *
        "Order      : $p\n\t" *
        "Num. Funcs.: $NF\n\t" *
        "DOFs/Func. : $N"
    print(io, str)
end

function basis(poly::InterpolatingPolynomial)
    return poly.basis
end

function dimension(poly::InterpolatingPolynomial)
    return dimension(poly.basis)
end

function order(poly::InterpolatingPolynomial)
    return order(poly.basis)
end

function update!(P::InterpolatingPolynomial, coeffs)
    P.coeffs[:] = coeffs[:]
end

################################################################################
function (P::InterpolatingPolynomial{1,B,T})(x) where {B<:AbstractBasis{1},T}
    return ((P.coeffs)*(P.basis(x)))[1]
end

function (P::InterpolatingPolynomial{1,B,T})(
    x::V,
) where {B<:AbstractBasis{1},T} where {V<:AbstractVector}
    return P(x[1])
end
################################################################################

################################################################################
function (P::InterpolatingPolynomial{1,B,T})(x, y) where {B<:AbstractBasis{2},T}
    return ((P.coeffs)*(P.basis(x, y)))[1]
end

function (P::InterpolatingPolynomial{1,B,T})(
    x::V,
) where {B<:AbstractBasis{2},T} where {V<:AbstractVector}
    @assert length(x) == 2
    return ((P.coeffs)*(P.basis(x[1], x[2])))[1]
end

function (P::InterpolatingPolynomial{N,B,T})(
    x,
    y,
) where {B<:AbstractBasis{2},N,T}
    return ((P.coeffs) * (P.basis(x, y)))
end

function (P::InterpolatingPolynomial{N,B,T})(
    x::V,
) where {B<:AbstractBasis{2},N,T} where {V<:AbstractVector}
    @assert length(x) == 2
    return ((P.coeffs) * (P.basis(x[1], x[2])))
end
################################################################################

function gradient(
    P::InterpolatingPolynomial{1,B,T},
    x,
) where {B<:AbstractBasis{1},T}
    return ((P.coeffs) * (gradient(P.basis, x)))
end

function gradient(
    P::InterpolatingPolynomial{1,B,T},
    x,
    y,
) where {B<:AbstractBasis{2},T}
    return ((P.coeffs) * (gradient(P.basis, x, y)))
end

function gradient(
    P::InterpolatingPolynomial{1,B,T},
    x::V,
) where {B<:AbstractBasis{2},T} where {V<:AbstractVector}
    return ((P.coeffs) * (gradient(P.basis, x)))
end

function hessian(
    P::InterpolatingPolynomial{1,B,T},
    x,
) where {B<:AbstractBasis{1},T}
    return ((P.coeffs) * (hessian(P.basis, x)))
end

function hessian(
    P::InterpolatingPolynomial{1,B,T},
    x,
    y,
) where {B<:AbstractBasis{2},T}

    return P.coeffs * hessian(P.basis, x, y)
end

function hessian(
    P::InterpolatingPolynomial{1,B,T},
    x::V,
) where {B<:AbstractBasis{2},T} where {V<:AbstractVector}
    return ((P.coeffs) * (hessian(P.basis, x)))
end


text:
function eosshow(io::IO, ::MIME"text/plain", model::EoSModel)
    print(io, typeof(model))
    model.lengthcomponents == 1 && println(io, " with 1 component:")
    model.lengthcomponents > 1 && println(io, " with ", model.lengthcomponents, " components:")
    for i in model.icomponents
        print(io, " \"", model.components[i], "\"")
        println(io)
    end
    print(io,"Contains parameters: ")
    firstloop = true
    for fieldname in fieldnames(typeof(model.params))
        firstloop == false && print(io, ", ")
        print(io, fieldname)
        firstloop = false
    end
end
function eosshow(io::IO, model::EoSModel)
    print(io, typeof(model))
    firstloop = true
    print(io, "(")
    for i in model.icomponents
        firstloop == false && print(io, ", ")
        print(io, "\"", model.components[i], "\"")
        firstloop = false
    end
    print(io, ")")
end

function Base.show(io::IO, ::MIME"text/plain", model::GCSAFTModel)
    print(io, typeof(model))
    model.lengthcomponents == 1 && println(io, " with 1 component:")
    model.lengthcomponents > 1 && println(io, " with ", model.lengthcomponents, " components:")
    for i in model.icomponents
        print(io, " \"", model.components[i], "\": ")
        firstloop = true
        for k in 1:length(model.allcomponentgroups[i])
            firstloop == false && print(io, ", ")
            print(io, "\"", model.allcomponentgroups[i][k], "\" => ", model.allcomponentngroups[i][k])
            firstloop = false
        end
        println(io)
    end
    print(io, "Contains parameters: ")
    firstloop = true
    for fieldname in fieldnames(typeof(model.params))
        firstloop == false && print(io, ", ")
        print(io, fieldname)
        firstloop = false
    end
end

function Base.show(io::IO, model::GCSAFTModel)
    print(io, typeof(model))
    firstloop = true
    for i in model.icomponents
        firstloop == false && print(io, ", ")
        print(io, "\"", model.components[i], "\"")
        firstloop = false
    end
    print(io, ")")
end


text:
mutable struct ROCK2ConstantCache{T,T2,zType} <: OrdinaryDiffEqConstantCache
  ms::SVector{46, Int}
  fp1::SVector{46, T}
  fp2::SVector{46, T}
  recf::Vector{T2}
  zprev::zType
  mdeg::Int
  deg_index::Int
  start::Int
  min_stage::Int
  max_stage::Int
end
@cache struct ROCK2Cache{uType,rateType,uNoUnitsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  uᵢ₋₁::uType
  uᵢ₋₂::uType
  tmp::uType
  atmp::uNoUnitsType
  fsalfirst::rateType
  k::rateType
  constantcache::ROCK2ConstantCache
end

function alg_cache(alg::ROCK2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  constantcache = ROCK2ConstantCache(uEltypeNoUnits, uEltypeNoUnits, u)
  uᵢ₋₁ = similar(u)
  uᵢ₋₂ = similar(u)
  tmp = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  fsalfirst = zero(rate_prototype)
  k = zero(rate_prototype)
  ROCK2Cache(u, uprev, uᵢ₋₁, uᵢ₋₂, tmp, atmp, fsalfirst, k, constantcache)
end

function alg_cache(alg::ROCK2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  ROCK2ConstantCache(uEltypeNoUnits, uEltypeNoUnits, u)
end

mutable struct ROCK4ConstantCache{T,T2,T3,T4,zType} <: OrdinaryDiffEqConstantCache
  ms::SVector{50, Int}
  fpa::Vector{T}
  fpb::Vector{T2}
  fpbe::Vector{T3}
  recf::Vector{T4}
  zprev::zType
  mdeg::Int
  deg_index::Int
  start::Int
  min_stage::Int
  max_stage::Int
end

@cache struct ROCK4Cache{uType,rateType,uNoUnitsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  uᵢ₋₁::uType
  uᵢ₋₂::uType
  uᵢ₋₃::uType
  tmp::uType
  atmp::uNoUnitsType
  fsalfirst::rateType
  k::rateType
  constantcache::ROCK4ConstantCache
end

function alg_cache(alg::ROCK4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  constantcache = ROCK4ConstantCache(real(uEltypeNoUnits), real(uEltypeNoUnits), u)
  uᵢ₋₁ = similar(u)
  uᵢ₋₂ = similar(u)
  uᵢ₋₃ = similar(u)
  tmp = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  fsalfirst = zero(rate_prototype)
  k = zero(rate_prototype)
  ROCK4Cache(u, uprev, uᵢ₋₁, uᵢ₋₂, uᵢ₋₃, tmp, atmp, fsalfirst, k, constantcache)
end

function alg_cache(alg::ROCK4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  ROCK4ConstantCache(real(uEltypeNoUnits), real(uEltypeNoUnits),u)
end

mutable struct RKCConstantCache{zType} <: OrdinaryDiffEqConstantCache
  #to match the types to call maxeig!
  zprev::zType
end
@cache struct RKCCache{uType,rateType,uNoUnitsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  gprev::uType
  gprev2::uType
  tmp::uType
  atmp::uNoUnitsType
  fsalfirst::rateType
  k::rateType
  constantcache::RKCConstantCache
end

function alg_cache(alg::RKC,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  constantcache = RKCConstantCache(u)
  gprev = similar(u)
  gprev2 = similar(u)
  tmp = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  fsalfirst = zero(rate_prototype)
  k = zero(rate_prototype)
  RKCCache(u, uprev, gprev, gprev2, tmp, atmp, fsalfirst, k, constantcache)
end

function alg_cache(alg::RKC,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  RKCConstantCache(u)
end

@cache mutable struct IRKCConstantCache{uType,rateType,N} <: OrdinaryDiffEqConstantCache
  minm::Int64
  zprev::uType
  nlsolver::N
  du₁::rateType
  du₂::rateType
end

@cache mutable struct IRKCCache{uType,rateType,uNoUnitsType,N} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  gprev::uType
  gprev2::uType
  fsalfirst::rateType
  f1ⱼ₋₁::rateType
  f1ⱼ₋₂::rateType
  f2ⱼ₋₁::rateType
  atmp::uNoUnitsType
  nlsolver::N
  du₁::rateType
  du₂::rateType
  constantcache::IRKCConstantCache
end

function alg_cache(alg::IRKC,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  γ, c = 1.0, 1.0
  J, W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = oopnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  zprev = u
  du₁ = rate_prototype; du₂ = rate_prototype
  IRKCConstantCache(50,zprev,nlsolver,du₁,du₂)
end

function alg_cache(alg::IRKC,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  γ, c = 1.0, 1.0
  J, W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = iipnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)

  gprev = similar(u)
  gprev2 = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  fsalfirst = zero(rate_prototype)
  zprev = similar(u)
  f1ⱼ₋₁ = zero(rate_prototype)
  f1ⱼ₋₂ = zero(rate_prototype)
  f2ⱼ₋₁ = zero(rate_prototype)
  du₁ = zero(rate_prototype)
  du₂ = zero(rate_prototype)
  constantcache = IRKCConstantCache(50,zprev,nlsolver,du₁,du₂)
  IRKCCache(u,uprev,gprev,gprev2,fsalfirst,f1ⱼ₋₁,f1ⱼ₋₂,f2ⱼ₋₁,atmp,nlsolver,du₁,du₂,constantcache)
end

mutable struct ESERK4ConstantCache{T, zType} <: OrdinaryDiffEqConstantCache
  ms::SVector{46, Int}
  Cᵤ::SVector{4, Int}
  Cₑ::SVector{4, Int}
  zprev::zType
  Bᵢ::Vector{T}
  mdeg::Int
  start::Int
  internal_deg::Int
end

@cache struct ESERK4Cache{uType,rateType,uNoUnitsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  uᵢ::uType
  uᵢ₋₁::uType
  uᵢ₋₂::uType
  Sᵢ::uType
  tmp::uType
  atmp::uNoUnitsType
  fsalfirst::rateType
  k::rateType
  constantcache::ESERK4ConstantCache
end

function alg_cache(alg::ESERK4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  constantcache = ESERK4ConstantCache(u)
  uᵢ = similar(u)
  uᵢ₋₁ = similar(u)
  uᵢ₋₂ = similar(u)
  Sᵢ   = similar(u)
  tmp = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  fsalfirst = zero(rate_prototype)
  k = zero(rate_prototype)
  ESERK4Cache(u, uprev, uᵢ, uᵢ₋₁, uᵢ₋₂, Sᵢ, tmp, atmp, fsalfirst, k, constantcache)
end

function alg_cache(alg::ESERK4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  ESERK4ConstantCache(u)
end

mutable struct ESERK5ConstantCache{T, zType} <: OrdinaryDiffEqConstantCache
  ms::SVector{49, Int}
  Cᵤ::SVector{5, Int}
  Cₑ::SVector{5, Int}
  zprev::zType
  Bᵢ::Vector{T}
  mdeg::Int
  start::Int
  internal_deg::Int
end

@cache struct ESERK5Cache{uType,rateType,uNoUnitsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  uᵢ::uType
  uᵢ₋₁::uType
  uᵢ₋₂::uType
  Sᵢ::uType
  tmp::uType
  atmp::uNoUnitsType
  fsalfirst::rateType
  k::rateType
  constantcache::ESERK5ConstantCache
end

function alg_cache(alg::ESERK5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  constantcache = ESERK5ConstantCache(u)
  uᵢ = similar(u)
  uᵢ₋₁ = similar(u)
  uᵢ₋₂ = similar(u)
  Sᵢ   = similar(u)
  tmp = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  fsalfirst = zero(rate_prototype)
  k = zero(rate_prototype)
  ESERK5Cache(u, uprev, uᵢ, uᵢ₋₁, uᵢ₋₂, Sᵢ, tmp, atmp, fsalfirst, k, constantcache)
end

function alg_cache(alg::ESERK5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  ESERK5ConstantCache(u)
end

mutable struct SERK2ConstantCache{T, zType} <: OrdinaryDiffEqConstantCache
  ms::SVector{11, Int}
  zprev::zType
  Bᵢ::Vector{T}
  mdeg::Int
  start::Int
  internal_deg::Int
end

@cache struct SERK2Cache{uType,rateType,uNoUnitsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  uᵢ₋₁::uType
  uᵢ₋₂::uType
  Sᵢ::uType
  tmp::uType
  atmp::uNoUnitsType
  fsalfirst::rateType
  k::rateType
  constantcache::SERK2ConstantCache
end

function alg_cache(alg::SERK2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  constantcache = SERK2ConstantCache(u)
  uᵢ₋₁ = similar(u)
  uᵢ₋₂ = similar(u)
  Sᵢ   = similar(u)
  tmp = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  fsalfirst = zero(rate_prototype)
  k = zero(rate_prototype)
  SERK2Cache(u, uprev, uᵢ₋₁, uᵢ₋₂, Sᵢ, tmp, atmp, fsalfirst, k, constantcache)
end

function alg_cache(alg::SERK2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  SERK2ConstantCache(u)
end


text:
println("I'm excited to learn Julia!")

my_variable = 42

print(my_variable)
print(" ")
println(typeof(my_variable))

🐱 = "Emoji variable name"

println(🐱)

# Comments use the # symbol

sum = 3 + 8
diff = 10 - 3
product = 20 * 2
quotient = 100 / 10
power = 10 ^ 3
mod = 101 % 7

# Type conversion
days = 365
days_float = float(days)
# Alternatives
# days_float = convert(Float64, days)
# days_float = convert(Float32, days)

@assert days == 365
@assert days_float == 365.0

max_stars_repo_name: AustinTSchaffer/DailyProgrammer

max_issues_repo_name: AustinTSchaffer/DailyProgrammer

max_forks_repo_name: AustinTSchaffer/DailyProgrammer


text:
###############################################################################
#
# Symbols
#
###############################################################################
"""
    latex_symbol(
                mid::String;
                sub::String = "",
                sup::String = "",
                presub::String = "",
                presup::String = "",
                option::String = "mathrm"
    )

Return the latex symbol string, given
- `mid` Center symbol, italic only when length>1
- `sub` Optional: subscript after the `mid`
- `sup` Optional: supscript after the `mid`
- `presub` Optional: subscript before the `mid`
- `presup` Optional: supscript before the `mid`
- `option` Optional: choose from `text` and `mathrm` (default)
"""
function latex_symbol(
            mid::String;
            sub::String = "",
            sup::String = "",
            presub::String = "",
            presup::String = "",
            option::String = "mathrm"
)
    # start with "\$ "
    symbol_string = "\$ ";

    # add the pre sub and sup
    symbol_string *= subscript(presub; option=option, standalone=false);
    symbol_string *= superscript(presup; option=option, standalone=false);

    # add the mid
    if length(mid) == 1
        symbol_string *= mid * " ";
    elseif length(mid) > 1
        symbol_string *= "\\$(option){" * mid * "} "
    end

    # add the sub and sup
    symbol_string *= subscript(sub; option=option, standalone=false);
    symbol_string *= superscript(sup; option=option, standalone=false);

    # end with "\$"
    symbol_string *= "\$";

    return symbol_string
end


text:
using Bukdu

struct PlugController <: ApplicationController
    conn::Conn
end

function index(::PlugController)
    "ok"
end

if PROGRAM_FILE == basename(@__FILE__)

pipeline(:web) do conn::Conn
    plug(Plug.CSRF.Protection, conn)
end

routes(:web) do
    get("/", PlugController, index)
end

Bukdu.start(8080)

Router.call(get, "/") #
# CLI.routes()

Base.JLOptions().isinteractive==0 && wait()

# Bukdu.stop()

end # if

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
using Test
using ThermoState
using ThermoState.Types
using ThermoState.QuickStates
using ThermoState: value
using Unitful

@testset "spec" begin
    #value
    @test value(spec(t=42))==42

    #specification
    @test specification(spec(spec"t",1)) == spec"t"

    #@spec_str keyword
    st = spec(t=1)
    st2 = spec(spec"t",1)
    @test st==st2

    #variable spec
    st2 = spec(spec"t",VariableSpec())
    @test st2(1) == st

    #units
    st2 = spec(spec"t",1u"K",false)
    @test st != st2

    st2 = spec(spec"t",1u"K")
    @test st == st2

end

@testset "mass consistency" begin
    #==
    test all possible combinations of:
    -mass, mass fractions,
    -moles, mol fractions,
    -mol numbers, mass numbers
    ==#

    n0 = [1.1,2.2,3.3] #mol
    m0 = [0.011,0.044,0.099] #mass
    mw0 = [10.0,20.0,30.0]u"g/mol" #mw, #testing for correct unit comversion
    xn0 = n0 ./ sum(n0)
    xm0 = m0 ./ sum(m0)
    mass0 = sum(m0)
    moles0 = sum(n0)
    _t = spec(t=1)
    _p = spec(p=2)
    _n = spec(n=n0)
    _m = spec(m=m0)
    _xm = spec(xm=xm0)
    _xn = spec(xn=xn0)
    _mass = spec(mass=mass0)
    _moles = spec(moles=moles0)

    a1 = state(mass=mass0,xm = xm0,t=1,p=2)
    @test mol_fraction(FromState(),a1,nothing,mw0)≈ xn0
    @test mass_fraction(FromState(),a1,nothing,mw0)≈ xm0
    @test mass_number(FromState(),a1,u"kg",mw0)≈ m0
    @test mol_number(FromState(),a1,u"mol",mw0)≈ n0
    @test mass(FromState(),a1,u"kg",mw0)≈ mass0
    @test moles(FromState(),a1,u"mol",mw0)≈ moles0

    a1 = state(_mass,_xm,_t,_p)
    @test mol_fraction(FromState(),a1,nothing,mw0)≈ xn0
    @test mass_fraction(FromState(),a1,nothing,mw0)≈ xm0
    @test mass_number(FromState(),a1,u"kg",mw0)≈ m0
    @test mol_number(FromState(),a1,u"mol",mw0)≈ n0
    @test mass(FromState(),a1,u"kg",mw0)≈ mass0
    @test moles(FromState(),a1,u"mol",mw0)≈ moles0

    a1 = state(moles=moles0,xm = xm0,t=1,p=2)
    @test mol_fraction(FromState(),a1,nothing,mw0)≈ xn0
    @test mass_fraction(FromState(),a1,nothing,mw0)≈ xm0
    @test mass_number(FromState(),a1,u"kg",mw0)≈ m0 #fail
    @test mol_number(FromState(),a1,u"mol",mw0)≈ n0 #fail
    @test mass(FromState(),a1,u"kg",mw0)≈ mass0 #fail
    @test moles(FromState(),a1,u"mol",mw0)≈ moles0

    a1 = state(_moles,_xm,_t,_p)
    @test mol_fraction(FromState(),a1,nothing,mw0)≈ xn0
    @test mass_fraction(FromState(),a1,nothing,mw0)≈ xm0
    @test mass_number(FromState(),a1,u"kg",mw0)≈ m0 #fail
    @test mol_number(FromState(),a1,u"mol",mw0)≈ n0 #fail
    @test mass(FromState(),a1,u"kg",mw0)≈ mass0 #fail
    @test moles(FromState(),a1,u"mol",mw0)≈ moles0

    a1 = state(moles=moles0,xn = xn0,t=1,p=2)
    @test mol_fraction(FromState(),a1,nothing,mw0)≈ xn0
    @test mass_fraction(FromState(),a1,nothing,mw0)≈ xm0
    @test mass_number(FromState(),a1,u"kg",mw0)≈ m0
    @test mol_number(FromState(),a1,u"mol",mw0)≈ n0
    @test mass(FromState(),a1,u"kg",mw0)≈ mass0
    @test moles(FromState(),a1,u"mol",mw0)≈ moles0

    a1 = state(_moles,_xn,_t,_p)
    @test mol_fraction(FromState(),a1,nothing,mw0)≈ xn0
    @test mass_fraction(FromState(),a1,nothing,mw0)≈ xm0
    @test mass_number(FromState(),a1,u"kg",mw0)≈ m0
    @test mol_number(FromState(),a1,u"mol",mw0)≈ n0
    @test mass(FromState(),a1,u"kg",mw0)≈ mass0
    @test moles(FromState(),a1,u"mol",mw0)≈ moles0

    a1 = state(mass=mass0,xm = xm0,t=1,p=2)
    @test mol_fraction(FromState(),a1,nothing,mw0)≈ xn0
    @test mass_fraction(FromState(),a1,nothing,mw0)≈ xm0
    @test mass_number(FromState(),a1,u"kg",mw0)≈ m0
    @test mol_number(FromState(),a1,u"mol",mw0)≈ n0
    @test mass(FromState(),a1,u"kg",mw0)≈ mass0
    @test moles(FromState(),a1,u"mol",mw0)≈ moles0

    a1 = state(_mass,_xm,_t,_p)
    @test mol_fraction(FromState(),a1,nothing,mw0)≈ xn0
    @test mass_fraction(FromState(),a1,nothing,mw0)≈ xm0
    @test mass_number(FromState(),a1,u"kg",mw0)≈ m0
    @test mol_number(FromState(),a1,u"mol",mw0)≈ n0
    @test mass(FromState(),a1,u"kg",mw0)≈ mass0
    @test moles(FromState(),a1,u"mol",mw0)≈ moles0

    a1 = state(m = m0,t=1,p=2)
    @test mol_fraction(FromState(),a1,nothing,mw0)≈ xn0
    @test mass_fraction(FromState(),a1,nothing,mw0)≈ xm0
    @test mass_number(FromState(),a1,u"kg",mw0)≈ m0
    @test mol_number(FromState(),a1,u"mol",mw0)≈ n0
    @test mass(FromState(),a1,u"kg",mw0)≈ mass0
    @test moles(FromState(),a1,u"mol",mw0)≈ moles0

    a1 = state(_m,_t,_p)
    @test mol_fraction(FromState(),a1,nothing,mw0)≈ xn0
    @test mass_fraction(FromState(),a1,nothing,mw0)≈ xm0
    @test mass_number(FromState(),a1,u"kg",mw0)≈ m0
    @test mol_number(FromState(),a1,u"mol",mw0)≈ n0
    @test mass(FromState(),a1,u"kg",mw0)≈ mass0
    @test moles(FromState(),a1,u"mol",mw0)≈ moles0

    a1 = state(n = n0,t=1,p=2)
    @test mol_fraction(FromState(),a1,nothing,mw0)≈ xn0
    @test mass_fraction(FromState(),a1,nothing,mw0)≈ xm0
    @test mass_number(FromState(),a1,u"kg",mw0)≈ m0
    @test mol_number(FromState(),a1,u"mol",mw0)≈ n0
    @test mass(FromState(),a1,u"kg",mw0)≈ mass0
    @test moles(FromState(),a1,u"mol",mw0)≈ moles0

    a1 = state(_n,_t,_p)
    @test mol_fraction(FromState(),a1,nothing,mw0)≈ xn0
    @test mass_fraction(FromState(),a1,nothing,mw0)≈ xm0
    @test mass_number(FromState(),a1,u"kg",mw0)≈ m0
    @test mol_number(FromState(),a1,u"mol",mw0)≈ n0
    @test mass(FromState(),a1,u"kg",mw0)≈ mass0
    @test moles(FromState(),a1,u"mol",mw0)≈ moles0

    mw1 = 34.3
    moles1 = 2.3
    mass1 = moles1*mw1*0.001
    _mass = spec(mass=mass1)
    _moles= spec(moles=moles1)
    a1 = state(t=1,p=2)
    @test mol_fraction(FromState(),a1,nothing,mw1)≈ [1.0]
    @test mass_fraction(FromState(),a1,nothing,mw1)≈ [1.0]
    @test mass_number(FromState(),a1,u"kg",mw1)≈ [0.001*mw1]
    @test mol_number(FromState(),a1,u"mol",mw1)≈ [1.0]
    @test mass(FromState(),a1,u"kg",mw1)≈ 0.001*mw1
    @test moles(FromState(),a1,u"mol",mw1)≈   1.0

    a1 = state(_t,_p)
    @test mol_fraction(FromState(),a1,nothing,mw1)≈ [1.0]
    @test mass_fraction(FromState(),a1,nothing,mw1)≈ [1.0]
    @test mass_number(FromState(),a1,u"kg",mw1)≈ [0.001*mw1]
    @test mol_number(FromState(),a1,u"mol",mw1)≈ [1.0]
    @test mass(FromState(),a1,u"kg",mw1)≈ 0.001*mw1
    @test moles(FromState(),a1,u"mol",mw1)≈   1.0

    a1 = state(t=1,p=2,mass=mass1)
    @test mol_fraction(FromState(),a1,nothing,mw1)≈ [1.0]
    @test mass_fraction(FromState(),a1,nothing,mw1)≈ [1.0]
    @test mass_number(FromState(),a1,u"kg",mw1)≈ [mass1]
    @test mol_number(FromState(),a1,u"mol",mw1)≈ [moles1]
    @test mass(FromState(),a1,u"kg",mw1)≈ mass1
    @test moles(FromState(),a1,u"mol",mw1)≈  moles1

    a1 = state(_t,_p,_mass)
    @test mol_fraction(FromState(),a1,nothing,mw1)≈ [1.0]
    @test mass_fraction(FromState(),a1,nothing,mw1)≈ [1.0]
    @test mass_number(FromState(),a1,u"kg",mw1)≈ [mass1]
    @test mol_number(FromState(),a1,u"mol",mw1)≈ [moles1]
    @test mass(FromState(),a1,u"kg",mw1)≈ mass1
    @test moles(FromState(),a1,u"mol",mw1)≈  moles1

    a1 = state(t=1,p=2,moles=moles1)
    @test mol_fraction(FromState(),a1,nothing,mw1)≈ [1.0]
    @test mass_fraction(FromState(),a1,nothing,mw1)≈ [1.0]
    @test mass_number(FromState(),a1,u"kg",mw1)≈ [mass1]
    @test mol_number(FromState(),a1,u"mol",mw1)≈ [moles1]
    @test mass(FromState(),a1,u"kg",mw1)≈ mass1
    @test moles(FromState(),a1,u"mol",mw1)≈  moles1

    a1 = state(_t,_p,_moles)
    @test mol_fraction(FromState(),a1,nothing,mw1)≈ [1.0]
    @test mass_fraction(FromState(),a1,nothing,mw1)≈ [1.0]
    @test mass_number(FromState(),a1,u"kg",mw1)≈ [mass1]
    @test mol_number(FromState(),a1,u"mol",mw1)≈ [moles1]
    @test mass(FromState(),a1,u"kg",mw1)≈ mass1
    @test moles(FromState(),a1,u"mol",mw1)≈  moles1

end

@testset "energy consistency" begin
    #==
    test all possible combinations an energy spec with::
    -mass, mass fractions,
    -moles, mol fractions,
    -mol numbers, mass numbers
    ==#

    n0 = [1.1,2.2,3.3] #mol
    m0 = [0.011,0.044,0.099] #mass
    mw0 = [10.0,20.0,30.0] #mw, #testing for correct unit comversion
    xn0 = n0 ./ sum(n0)
    xm0 = m0 ./ sum(m0)
    mass0 = sum(m0)
    moles0 = sum(n0)
    mol_h0 = 101.433
    total_h0 = mol_h0*moles0
    mass_h0 = total_h0/mass0
    a1 = state(n = n0,t=1,mol_h = mol_h0)
    @test mol_enthalpy(FromState(),a1,u"J/mol",mw0) ≈ mol_h0
    @test mass_enthalpy(FromState(),a1,u"J/kg",mw0) ≈ mass_h0
    @test total_enthalpy(FromState(),a1,u"J",mw0) ≈ total_h0

    a1 = state(n = n0,t=1,mass_h = mass_h0)
    @test mol_enthalpy(FromState(),a1,u"J/mol",mw0) ≈ mol_h0
    @test mass_enthalpy(FromState(),a1,u"J/kg",mw0) ≈ mass_h0
    @test total_enthalpy(FromState(),a1,u"J",mw0) ≈ total_h0

    a1 = state(n = n0,t=1,total_h = total_h0)
    @test mol_enthalpy(FromState(),a1,u"J/mol",mw0) ≈ mol_h0
    @test mass_enthalpy(FromState(),a1,u"J/kg",mw0) ≈ mass_h0
    @test total_enthalpy(FromState(),a1,u"J",mw0) ≈ total_h0

   #test entropy
   a1 = state(n = n0,t=1,mol_s = mol_h0)
   @test mol_entropy(FromState(),a1,u"J/(mol*K)",mw0) ≈ mol_h0
   @test mass_entropy(FromState(),a1,u"J/(kg*K)",mw0) ≈ mass_h0
   @test total_entropy(FromState(),a1,u"J/K",mw0) ≈ total_h0

   a1 = state(n = n0,t=1,mass_s = mass_h0)
   @test mol_entropy(FromState(),a1,u"J/(mol*K)",mw0) ≈ mol_h0
   @test mass_entropy(FromState(),a1,u"J/(kg*K)",mw0) ≈ mass_h0
   @test total_entropy(FromState(),a1,u"J/K",mw0) ≈ total_h0

   a1 = state(n = n0,t=1,total_s = total_h0)
   @test mol_entropy(FromState(),a1,u"J/(mol*K)",mw0) ≈ mol_h0
   @test mass_entropy(FromState(),a1,u"J/(kg*K)",mw0) ≈ mass_h0
   @test total_entropy(FromState(),a1,u"J/K",mw0) ≈ total_h0
end

@testset "volume consistency" begin

n0 = [1.1,2.2,3.3] #mol
m0 = [0.011,0.044,0.099] #mass
mw0 = [10.0,20.0,30.0]#mw, #testing for correct unit comversion
xn0 = n0 ./ sum(n0)
xm0 = m0 ./ sum(m0)
mass0 = sum(m0)
moles0 = sum(n0)
V0 = 400.0 #m3
mol_v0 = V0/moles0
mol_rho0 = 1/mol_v0
mass_v0 = V0/mass0
mass_rho0 = 1/mass_v0

a1 = state(n = n0,t=1,total_v = V0)
@test total_volume(FromState(),a1,u"m^3",mw0) ≈ V0
@test mol_volume(FromState(),a1,u"m^3/mol",mw0) ≈ mol_v0
@test mass_volume(FromState(),a1,u"m^3/kg",mw0) ≈ mass_v0
@test mol_density(FromState(),a1,u"mol/m^3",mw0) ≈ mol_rho0
@test mass_density(FromState(),a1,u"kg/m^3",mw0) ≈ mass_rho0

a1 = state(n = n0,t=1,mass_v = mass_v0)
@test total_volume(FromState(),a1,u"m^3",mw0) ≈ V0
@test mol_volume(FromState(),a1,u"m^3/mol",mw0) ≈ mol_v0
@test mass_volume(FromState(),a1,u"m^3/kg",mw0) ≈ mass_v0
@test mol_density(FromState(),a1,u"mol/m^3",mw0) ≈ mol_rho0
@test mass_density(FromState(),a1,u"kg/m^3",mw0) ≈ mass_rho0

a1 = state(n = n0,t=1,mol_v = mol_v0)
@test total_volume(FromState(),a1,u"m^3",mw0) ≈ V0
@test mol_volume(FromState(),a1,u"m^3/mol",mw0) ≈ mol_v0
@test mass_volume(FromState(),a1,u"m^3/kg",mw0) ≈ mass_v0
@test mol_density(FromState(),a1,u"mol/m^3",mw0) ≈ mol_rho0
@test mass_density(FromState(),a1,u"kg/m^3",mw0) ≈ mass_rho0

a1 = state(n = n0,t=1,mol_rho = mol_rho0)
@test total_volume(FromState(),a1,u"m^3",mw0) ≈ V0
@test mol_volume(FromState(),a1,u"m^3/mol",mw0) ≈ mol_v0
@test mass_volume(FromState(),a1,u"m^3/kg",mw0) ≈ mass_v0
@test mol_density(FromState(),a1,u"mol/m^3",mw0) ≈ mol_rho0
@test mass_density(FromState(),a1,u"kg/m^3",mw0) ≈ mass_rho0

a1 = state(n = n0,t=1,mass_rho = mass_rho0)
@test total_volume(FromState(),a1,u"m^3",mw0) ≈ V0
@test mol_volume(FromState(),a1,u"m^3/mol",mw0) ≈ mol_v0
@test mass_volume(FromState(),a1,u"m^3/kg",mw0) ≈ mass_v0
@test mol_density(FromState(),a1,u"mol/m^3",mw0) ≈ mol_rho0
@test mass_density(FromState(),a1,u"kg/m^3",mw0) ≈ mass_rho0

end

@testset "state points" begin
using ThermoState.StatePoints

@test pressure(FromState(),StandardConditions()) == 100000.0
@test pressure(FromState(),NormalConditions()) == 101325.0
@test pressure(FromState(),NormalBoilingPoint()) == 101325.0
@test temperature(FromState(),StandardConditions()) == 273.15
@test temperature(FromState(),NormalConditions()) == 293.15
end

@testset "state errors" begin

    #@test_throws Exception state(t=1,T=1) #broken
    mass0 = spec(spec"mass",1)
    moles0 = spec(spec"moles",1)
    t0 = spec(spec"t",1)
    p0 = spec(spec"p",1)

    @test_throws Exception state(mass0,moles0,t0,p0)
end

@testset "variable spec" begin
    λ = VariableSpec()
    p0 = 1
    t0 = 1
    moles0 = 1
    st0 = state(t=t0,p=p0,moles=moles0)
    st1 = state(t=λ,p=p0,moles=moles0)
    @test st1(t0) == st0
    st1 = state(t=λ,p=λ,moles=moles0)
    @test st1(t0,p0) == st0
    st1 = state(t=λ,p=λ,moles=λ)
    @test st1(t0,p0,moles0) == st0

    sp0 = spec(t=t0)
    sp1 = spec(t=λ)
    @test sp1(t0) == sp0
end

@testset "state_type" begin
    st = state(t=1,p=2)
    @test state_type(st) isa SinglePT
    @test QuickStates.pt() isa SinglePT

    st = state(t=1,p=2,n=rand(5))
    @test state_type(st) isa MultiPT
    @test QuickStates.ptx() isa MultiPT
    @test QuickStates.ptn() isa MultiPT

    st = state(v=1,p=2)
    @test state_type(st) isa SinglePV
    @test QuickStates.pv() isa SinglePV

    st = state(v=1,p=2,n=rand(5))
    @test state_type(st) isa MultiPV
    @test QuickStates.pvx() isa MultiPV
    @test QuickStates.pvn() isa MultiPV

    st = state(t=2,v=1)
    @test state_type(st) isa SingleVT
    @test QuickStates.vt() isa SingleVT

    st = state(t=2,v=1,n=rand(5))
    @test state_type(st) isa MultiVT
    @test QuickStates.vtx() isa MultiVT
    @test QuickStates.vtn() isa MultiVT

    st = state(h=2,p=1)
    @test state_type(st) isa SinglePH
    @test QuickStates.ph() isa SinglePH

    st = state(h=2,p=1,n=rand(5))
    @test state_type(st) isa MultiPH
    @test QuickStates.phx() isa MultiPH
    @test QuickStates.phn() isa MultiPH

    st = state(sat=true,p=1)
    @test state_type(st) isa SingleSatP
    @test QuickStates.sat_p() isa SingleSatP

    st = state(sat=true,p=1,n=rand(5))
    @test state_type(st) isa MultiSatP
    @test QuickStates.sat_px() isa MultiSatP
    @test QuickStates.sat_pn() isa MultiSatP

    st = state(sat=true,t=1)
    @test state_type(st) isa SingleSatT
    @test QuickStates.sat_t() isa SingleSatT

    st = state(sat=true,t=1,n=rand(5))
    @test state_type(st) isa MultiSatT
    @test QuickStates.sat_tx() isa MultiSatT
    @test QuickStates.sat_tn() isa MultiSatT

end

@testset "@to_units" begin
    xn = rand(5)
    xn = xn ./ sum(xn)
    st = state(sat=true,t=1u"K",xn=xn)
    @test (@to_units temperature(FromState(),st)) == 1u"K"
    @test (@to_units temperature(FromState(),st,u"°C")) ==-272.15u"°C"
    @test (@to_units mol_fraction(FromState(),st)) == xn

end


text:
@cdef :PyBytes_FromStringAndSize PyPtr (Ptr{Cchar}, Py_ssize_t)
@cdef :PyBytes_AsStringAndSize Cint (PyPtr, Ptr{Ptr{Cchar}}, Ptr{Py_ssize_t})

const PyBytes_Type__ref = Ref(PyNULL)
PyBytes_Type() = pyglobal(PyBytes_Type__ref, :PyBytes_Type)

PyBytes_Check(o) = Py_TypeCheckFast(o, Py_TPFLAGS_BYTES_SUBCLASS)
PyBytes_CheckExact(o) = Py_TypeCheckExact(o, PyBytes_Type())

PyBytes_From(s::Union{Vector{Cuchar},Vector{Cchar},String,SubString{String}}) =
    PyBytes_FromStringAndSize(pointer(s), sizeof(s))

PyBytes_AsString(o) = begin
    ptr = Ref{Ptr{Cchar}}()
    len = Ref{Py_ssize_t}()
    err = PyBytes_AsStringAndSize(o, ptr, len)
    ism1(err) && return ""
    Base.unsafe_string(ptr[], len[])
end

PyBytes_AsVector(o, ::Type{T} = UInt8) where {T} = begin
    T in (Int8, UInt8) || throw(MethodError(PyBytes_AsVector, (o, T)))
    ptr = Ref{Ptr{Cchar}}()
    len = Ref{Py_ssize_t}()
    err = PyBytes_AsStringAndSize(o, ptr, len)
    ism1(err) && return T[]
    copy(Base.unsafe_wrap(Vector{T}, Ptr{T}(ptr[]), len[]))
end

PyBytes_TryConvertRule_vector(o, ::Type{Vector{X}}) where {X} = begin
    v = PyBytes_AsVector(o, X)
    isempty(v) && PyErr_IsSet() && return -1
    return putresult(v)
end

PyBytes_TryConvertRule_string(o, ::Type{String}) = begin
    v = PyBytes_AsString(o)
    isempty(v) && PyErr_IsSet() && return -1
    return putresult(v)
end


text:
module BlockSystems

using LinearAlgebra
using DocStringExtensions
using ModelingToolkit
using ModelingToolkit: ODESystem, Differential
using ModelingToolkit: get_iv, get_eqs, get_states
using ModelingToolkit: rename, getname, renamespace, namespace_equation, namespace_equations, value, vars
using ModelingToolkit: equation_dependencies, asgraph, variable_dependencies, eqeq_dependencies, varvar_dependencies
using ModelingToolkit.SymbolicUtils: Symbolic, operation, arguments, istree
using ModelingToolkit.Symbolics: tosymbol
using Graphs

export AbstractIOSystem, IOBlock, IOSystem, get_iv, equations

include("utils.jl")

"""
abstract supertype for [`IOBlock`](@ref) and [`IOSystem`](@ref).
"""
abstract type AbstractIOSystem end

function Base.getproperty(sys::AbstractIOSystem, name::Symbol)
    if name ∈ fieldnames(typeof(sys))
        return getfield(sys, name)
    end
    syms = vcat(getfield(sys, :inputs),
                getfield(sys, :iparams),
                getfield(sys, :istates),
                getfield(sys, :outputs),
                getfield(sys, :removed_states))

    i = findfirst(x -> getname(x) == name, syms)
    if i !== nothing
        return rename(syms[i], renamespace(getfield(sys, :name), name))
    end
    throw(error("Variable $name does not exist"))
end

function namespace_property(ios::AbstractIOSystem, prop::Symbol)
    [rename(sym, renamespace(ios.name, getname(sym)))
     for sym in getproperty(ios, prop)]
end
namespace_inputs(ios::AbstractIOSystem)  = namespace_property(ios, :inputs)
namespace_iparams(ios::AbstractIOSystem) = namespace_property(ios, :iparams)
namespace_istates(ios::AbstractIOSystem) = namespace_property(ios, :istates)
namespace_outputs(ios::AbstractIOSystem) = namespace_property(ios, :outputs)
namespace_rem_states(ios::AbstractIOSystem) = namespace_property(ios, :removed_states)

"""
$(TYPEDEF)

A basic IOSystem which consists of a single ODESystem.
$(TYPEDFIELDS)
"""
struct IOBlock <: AbstractIOSystem
    name::Symbol
    inputs::Vector{Symbolic}
    iparams::Vector{Symbolic}
    istates::Vector{Symbolic}
    outputs::Vector{Symbolic}
    system::ODESystem
    removed_states::Vector{Symbolic}
    removed_eqs::Vector{Equation}

    function IOBlock(name, inputs, iparams, istates, outputs, odes, rem_eqs; warn)
        @check name == getname(odes) "Name of inner ODESystem does not match name of IOBlock"
        @checkwarn warn Set(inputs) ⊆ Set(parameters(odes)) "Inputs should be parameters. You may ignore this warning if you want to speciy an input which is not used in the eqs."
        @checkwarn warn Set(outputs) ⊆ Set(states(odes)) "Outputs should be variables. You may ignore this waring if you want to specify an output which is not used in the eqs (completly implicit)."
        @checkwarn warn Set(iparams) ⊆ Set(parameters(odes)) "iparams should be parameters of the eqs system"
        @checkwarn warn Set(istates) ⊆ Set(states(odes)) "istates should be variables of the eqs system"
        @check Set(parameters(odes)) ⊆ Set(inputs ∪ iparams) "params(eqs) ⊆ inputs ∪ iparams"
        @check Set(states(odes)) ⊆ Set(outputs ∪ istates) "states(eqs) ⊆ outputs ∪ istates"

        # lhs of the equation should be either differential of state/output
        for eq in equations(odes)
            (type, lhs_var) = eq_type(eq)
            if type === :diffeq
                @check lhs_var ∈ Set(istates ∪ outputs) "$eq : Lhs variable of diff eqs should be istate or output."
            elseif type === :explicit_algebraic
                @check lhs_var ∈ Set(istates ∪ outputs) "$eq : Lhs variable of explicit algebraic eqs should be istate or output."
            elseif type === :implicit_algebraic
                @checkwarn warn isnothing(lhs_var) "$eq : Please provide constraints in form `0 ~ f(...)`."
            end
        end

        additional_vars = get_variables([eq.rhs for eq in rem_eqs])
        all_vars = vcat(inputs, outputs, istates, iparams, get_iv(odes))
        @check uniquenames(all_vars) "There seem to be a name clashes between inputs, iparams istates and outputs!"
        all_vars = Set(all_vars)
        @checkwarn warn additional_vars ⊆ all_vars "Removed eqs should not contain new variables. This may be a bug since this feature is not fully implemented."

        if isempty(rem_eqs)
            rem_states = Symbolic[]
        else
            rem_states = lhs_var.(rem_eqs)
            @checkwarn warn isempty(rem_states ∩ all_vars) "removed states should not appear in in/out/is/ip equations"
            rem_rhs_vars = union([get_variables(eq.rhs) for eq in rem_eqs]...)
            if !(rem_rhs_vars ⊆ all_vars)
                remaining = setdiff(rem_rhs_vars, all_vars)
                @warn "RHS variables of removed eqs are not a subset of in/out/is/ip. Maybe $remaining got totally removed during simplification? This will lead to errors if you wan't to build functions for removed states."
            end
        end

        # TODO: check IOBlock assumptions in inner constructor
        # - each state is represented by one lhs (static or diff) or implicit algebraic equation?
        # - lhs only first order or algebraic
        new(name, inputs, iparams, istates, outputs, odes, rem_states, rem_eqs)
    end
end

ModelingToolkit.get_iv(block::IOBlock) = get_iv(block.system)
ModelingToolkit.equations(block::IOBlock) = equations(block.system)

function namespace_rem_eqs(iob::IOBlock)
    eqs = iob.removed_eqs
    isempty(eqs) && return Equation[]

    map(eq->namespace_equation(eq, iob.system), eqs)
end

"""
$(SIGNATURES)

Construct a new IOBlock for the given arguments.

```@example
using BlockSystems, ModelingToolkit
@parameters t i(t)
@variables x(t) o(t)
D = Differential(t)

iob = IOBlock([D(x) ~ i, o ~ x], [i], [o], name=:iob)
```
"""
function IOBlock(eqs::Vector{<:Equation}, inputs, outputs; name = gensym(:IOBlock), iv = nothing, warn=true)
    IOBlock(name, eqs, inputs, outputs, Equation[]; iv, warn)
end

function IOBlock(name, eqs, inputs, outputs, rem_eqs; iv=nothing, warn=true)
    os = ODESystem(eqs, iv; name = name)

    inputs = value.(inputs) # gets the inputs as `tern` type
    outputs = value.(outputs) # gets the outputs as `tern` type
    istates = setdiff(get_states(os), outputs)
    iparams = setdiff(parameters(os), inputs)

    IOBlock(name, inputs, iparams, istates, outputs, os, rem_eqs; warn)
end

"""
$(SIGNATURES)

Construct a new IOBlock based on an existing. Deep-copy all fields and assigns new name.
"""
function IOBlock(iob::IOBlock; name=gensym(:IOBlock), warn=true)
    cp = deepcopy(iob)
    odes = ODESystem(get_eqs(cp.system), get_iv(cp), name=name)
    IOBlock(name, cp.inputs, cp.iparams, cp.istates, cp.outputs, odes, cp.removed_eqs; warn)
end

"""
$(TYPEDEF)

A composite `IOSystem` which consists of multiple [`AbstractIOSystem`](@ref) which are connected via
a vector of namespaced pairs (`subsys1.out => subsys2.in`).

An `IOSystem` contains maps how to promote the namespaced variables of the subsystem to the new scope
   subsys1₊x(t) => x(t)
   subsys1₊y(t) => subsys1₊y(t)
   subsys2₊y(t) => subsys2₊y(t)

$(TYPEDFIELDS)
"""
struct IOSystem <: AbstractIOSystem
    name::Symbol
    inputs::Vector{Symbolic}
    iparams::Vector{Symbolic}
    istates::Vector{Symbolic}
    outputs::Vector{Symbolic}
    removed_states::Vector{Symbolic}
    connections::Vector{Pair{Symbolic, Symbolic}}
    namespace_map::Dict{Symbolic, Symbolic}
    systems::Vector{AbstractIOSystem}
end
# TODO: check IOSystem assumptions in inner constructor?
# - check that every variable found in the equations is referenced by namespace map

ModelingToolkit.get_iv(sys::IOSystem) = get_iv(first(sys.systems))

"""
$(SIGNATURES)

Construct a new IOSystem from various subsystems.
Arguments:
 - `cons`:

    The connections in the form `sub1.output => sub2.input`. It is also possible to use simple
    algebraic equations such as `sub1.o1 + sub2.o2 => sub3.input`.

 - `io_systems`: Vector of subsystems
 - `namespace_map`: Provide collection of custom namespace promotions / renamings
   i.e. sub1.input => voltage`. Variables without entry in the map will be
   promoted automatically. Automatic promotion means that the sub-namespace is
   removed whenever it is possible without naming conflicts. The map may contain
   inputs, outputs, istates, iparams and removed_states. The rhs of the map can be provided
   as as Symbol: `sub1.input => :newname`.
 - `outputs`: Per default, all of the subsystem outputs will become system outputs. However,
   by providing a list of variables as outputs *only these* will become outputs of the
   new system. All other sub-outputs will become internal states of the connected system
   (and might be optimized away in `connect_system`).
 - `name`: namespace
 - `autopromote=true`: enable/disable automatic promotion of variable names to system namespace
"""
function IOSystem(cons,
                  io_systems::Vector{<:AbstractIOSystem};
                  namespace_map = nothing,
                  outputs = :all,
                  name = gensym(:IOSystem),
                  autopromote = true)
    namespaces = [sys.name for sys in io_systems]
    allunique(namespaces) || throw(ArgumentError("Namespace collision in subsystems!"))

    ivs = unique([get_iv(sys) for sys in io_systems])
    length(ivs) == 1 || throw(ArgumentError("Multiple independent variables!"))

    nspcd_inputs = vcat([namespace_inputs(sys) for sys in io_systems]...)
    nspcd_iparams = vcat([namespace_iparams(sys) for sys in io_systems]...)
    nspcd_istates = vcat([namespace_istates(sys) for sys in io_systems]...)
    nspcd_outputs = vcat([namespace_outputs(sys) for sys in io_systems]...)
    nspcd_rem_states = vcat([namespace_rem_states(sys) for sys in io_systems]...)

    # check validity of provided connections
    @check vars(first.(cons)) ⊆ Set(nspcd_outputs) "First argument in connection may only contain outputs of subsystem."
    @check Set(last.(cons)) ⊆ Set(nspcd_inputs) "Second argument in connection needs to be input of subsystem."
    @check allunique(last.(cons)) "Multiple connections to same input!"

    # reduce the inputs to the open inputs
    open_inputs = setdiff(nspcd_inputs, last.(cons))

    # check validity of provided namespace maps
    namespace_map = fix_map_types(namespace_map)
    allowed_lhs = Set(open_inputs ∪ nspcd_iparams ∪ nspcd_istates ∪ nspcd_outputs ∪ nspcd_rem_states)
    @check keys(namespace_map) ⊆ Set(allowed_lhs) "namespace_map keys contain unknown symbols"
    @check uniquenames(values(namespace_map)) "naming conflict in rhs of user provided namespace promotions"

    # if the user provided a list of outputs, all other outputs become istates
    if outputs != :all
        outputs::Vector{Union{Symbol, Symbolic}} = value.(outputs)
        # check if outputs ∈ nspcd_outputs or referenced in rhs of namespace_map
        for (i, o) in enumerate(outputs)
            if o ∉ Set(nspcd_outputs)
                if o isa Symbol
                    key = findfirst(v->isequal(getname(v), o), namespace_map)
                else
                    key = findfirst(v->isequal(v, o), namespace_map)
                end
                # check if key references a output (namespace_map contains all)
                key ∉ Set(nspcd_outputs) && throw(ArgumentError("output $o ∉ outputs ∪ outputs_promoted"))
                outputs[i] = key
            end
        end
        former_outputs = setdiff(nspcd_outputs, outputs)
        nspcd_istates  = vcat(nspcd_istates, former_outputs)
        nspcd_outputs  = outputs |> collect
    end

    # auto promotion of the rest
    all_symbols = vcat(open_inputs, nspcd_iparams, nspcd_istates, nspcd_outputs, nspcd_rem_states)
    left_symbols = setdiff(all_symbols, keys(namespace_map))

    auto_promotions = create_namespace_promotions(collect(left_symbols), values(namespace_map), autopromote)

    namespace_map = merge(namespace_map, auto_promotions)

    # assert that there are no namespace clashes. this should be allways true!
    @assert uniquenames(keys(namespace_map)) "ERROR: complete namespace map lhs not unique"
    @assert uniquenames(values(namespace_map)) "ERROR: complete namespace map lhs not unique"

    nmspc_promote = ModelingToolkit.substituter(namespace_map)
    inputs  = nmspc_promote.(open_inputs)
    iparams = nmspc_promote.(nspcd_iparams)
    istates = nmspc_promote.(nspcd_istates)
    outputs = nmspc_promote.(nspcd_outputs)
    rem_states = nmspc_promote.(nspcd_rem_states)

    IOSystem(name,
             inputs, iparams, istates, outputs, rem_states,
             cons,
             namespace_map,
             io_systems)
end

export BlockSpec, fulfills

"""
    struct BlockSpec
    BlockSpec(in::Vector, out::Vector; in_strict=true, out_strict=false)

Block specification, defines which inputs/outputs an [`AbstractIOSystem`](@ref)
should have. Contains two vectors of `Symbols`. Can be initialized with Vectors
of `Symbols`, `Num` or `<:Symbolic`.

If `strict=true` the in/outputs *must equal* the specification. If `strict=false`
the block *must contain* the in/outputs from the specification.

Object is functor: call `(::BlockSpec)(ios)` to check whether `ios` fulfills
specification. See also [`fulfills`](@ref).

```
iob = IOBlock(...)
spec = BlockSpec([:uᵣ, :uᵢ], [:iᵣ, :iᵢ])
fulfills(iob, spec)
spec(iob)
```
"""
struct BlockSpec
    inputs::Vector{Symbol}
    outputs::Vector{Symbol}
    in_strict::Bool
    out_strict::Bool
end

BlockSpec(in::Vector, out::Vector; in_strict=true, out_strict=false) = BlockSpec(in, out, in_strict, out_strict)
BlockSpec(in::Vector{Num}, out::Vector{Num}; kwargs...) = BlockSpec(value.(in), value.(out); kwargs...)
BlockSpec(in::Vector{<:Symbolic}, out::Vector{<:Symbolic}; kwargs...) = BlockSpec(getname.(in), getname.(out); kwargs...)

(bs::BlockSpec)(io) = fulfills(io, bs)

"""
    fulfills(io, bs::BlockSpec)::Bool

Check whether `io` fulfills the given [`BlockSpec`](@ref).
"""
function fulfills(io, bs::BlockSpec)
    if bs.in_strict
        Set(bs.inputs) == Set(getname.(io.inputs))
    else
        Set(bs.inputs) ⊆ Set(getname.(io.inputs))
    end && if bs.out_strict
        Set(bs.outputs) == Set(getname.(io.outputs))
    else
        Set(bs.outputs) ⊆ Set(getname.(io.outputs))
    end
end

"""
    fix_map_types(map)

Creates Dict from `map`. Changes `Num`-types to `Symbolic`-types in the
user provided maps. Map can be Dict or  Array of Pairs.
If the rhs is of type `Symbol` it will be converted to a `Symbolic` type.
"""
function fix_map_types(map)
    dict = Dict(map)
    newdict = Dict{Symbolic, Symbolic}()
    for k in keys(dict)
        newkey = value(k)
        if dict[k] isa Symbol
            newval = rename(newkey, dict[k])
        else
            newval = value(dict[k])
        end
        newdict[newkey] = newval
    end
    newdict
end
fix_map_types(::Nothing) = Dict()

"""
    create_namepsace_promotions(syms, forbidden, autopromote)

Takes two lists of Symbols, `syms` and `forbidden`. Returns a Dict of namespace
promotions for each symbol in `syms`. Avoids name collisions inside the list as well
as with the `forbidden` symbols.
"""
function create_namespace_promotions(syms, forbidden, autopromote)
    dict = Dict()
    if autopromote
        promoted = remove_namespace.(syms)
        for i in 1:length(syms)
            forbidden_names = getname.(forbidden) ∪ getname.(promoted[[1:i-1;i+1:end]])
            dict[syms[i]] = getname(promoted[i])∈forbidden_names ? syms[i] : promoted[i]
        end
    else
        for i in 1:length(syms)
            dict[syms[i]] = syms[i]
        end
    end
    return dict
end

include("transformations.jl")
include("function_generation.jl")
include("visualization.jl")

end


text:
export AnchorNode

import ..Nodes.draw

# This is an opionated node. It will automatically inject a filter node
# between the given parent and this anchor node.

mutable struct AnchorNode <: Ranger.AbstractNode
    base::Nodes.NodeData
    transform::Nodes.TransformProperties{Float64}

    # Collection of nodes.
    children::Array{Ranger.AbstractNode,1}

    mesh::Geometry.Mesh

    color::Rendering.Palette

    function AnchorNode(world::Ranger.World, name::String, parent::Ranger.AbstractNode)
        o = new()

        # Automatically insert a filter between the parent and this anchor node.
        # filter uses the default behaviour of the transform filter.
        filter = Filters.TransformFilter(world, "AnchorTransformFilter", parent)
        # Add the filter as a child of parent.
        push!(parent.children, filter)
        # Now make the anchor a child of the filter.
        push!(filter.children, o)

        o.base = Nodes.NodeData(name, filter, world)
        o.transform = Nodes.TransformProperties{Float64}()
        o.children = Array{Ranger.AbstractNode,1}[]
        o.mesh = Geometry.Mesh()
        o.color = Rendering.White()

        # horizontal
        Geometry.add_vertex!(o.mesh, -0.5, -0.5)
        Geometry.add_vertex!(o.mesh, 0.5, 0.5)

        # vertical
        Geometry.add_vertex!(o.mesh, -0.5, 0.5)
        Geometry.add_vertex!(o.mesh, 0.5, -0.5)

        Geometry.build!(o.mesh)

        o
    end
end

function Nodes.draw(node::AnchorNode, context::Rendering.RenderContext)
    if Nodes.is_dirty(node)
        Rendering.transform!(context, node.mesh)
        Nodes.set_dirty!(node, false)
    end

    Rendering.set_draw_color(context, node.color)

    v1 = node.mesh.bucket[1]
    v2 = node.mesh.bucket[2]
    v3 = node.mesh.bucket[3]
    v4 = node.mesh.bucket[4]
    Rendering.draw_line(context, v1.x, v1.y, v2.x, v2.y)

    Rendering.draw_line(context, v3.x, v3.y, v4.x, v4.y);
end

# --------------------------------------------------------
# Grouping
# --------------------------------------------------------
function Nodes.get_children(node::AnchorNode)
    node.children
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
using Documenter, SQLite

makedocs(
    modules = [SQLite],
    sitename = "SQLite.jl",
    pages = ["Home" => "index.md"]
)

deploydocs(
    repo = "github.com/JuliaDatabases/SQLite.jl.git",
    target = "build",
)


text:
# AUTO GENERATED FILE - DO NOT EDIT

export fileexplorer

"""
    fileexplorer(;kwargs...)

A FileExplorer component.

Keyword arguments:
- `id` (String; optional)
- `selectedPath` (String; optional)
- `value` (Array; optional)
"""
function fileexplorer(; kwargs...)
        available_props = Symbol[:id, :selectedPath, :value]
        wild_props = Symbol[]
        return Component("fileexplorer", "FileExplorer", "dash_cool_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: raulrozza/dash-cool-components

max_issues_repo_name: raulrozza/dash-cool-components

max_forks_repo_name: raulrozza/dash-cool-components


text:
int_rules_4_3_12 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.3.12*(d+e*x)^m*tan(a+b*x+c*x^2)^n =#
    @apply_utils Antiderivative(tan(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~(n'), ~x) => Unintegrable(tan(~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x)
    @apply_utils Antiderivative(cot(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~(n'), ~x) => Unintegrable(cot(~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x)
    @apply_utils Antiderivative((~d + ~(e') * ~x) * tan(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => (-(~e) * log(cos(~a + ~b * ~x + ~c * (~x) ^ 2))) / (2 * ~c) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && EqQ(2 * ~c * ~d - ~b * ~e, 0)
    @apply_utils Antiderivative((~d + ~(e') * ~x) * cot(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => (~e * log(sin(~a + ~b * ~x + ~c * (~x) ^ 2))) / (2 * ~c) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && EqQ(2 * ~c * ~d - ~b * ~e, 0)
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) * tan(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => (-(~e) * log(cos(~a + ~b * ~x + ~c * (~x) ^ 2))) / (2 * ~c) + ((2 * ~c * ~d - ~b * ~e) / (2 * ~c)) * Antiderivative(tan(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && NeQ(2 * ~c * ~d - ~b * ~e, 0)
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) * cot(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => (~e * log(sin(~a + ~b * ~x + ~c * (~x) ^ 2))) / (2 * ~c) + ((2 * ~c * ~d - ~b * ~e) / (2 * ~c)) * Antiderivative(cot(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && NeQ(2 * ~c * ~d - ~b * ~e, 0)
    #= Antiderivative((~x)^(~m)*tan((~a')+(~b')*(~x)+(~c')*(~x)^2),~x) :=  -x^(m-1)*log(cos(a+b*x+c*x^2))/(2*c) - b/(2*c)*Antiderivative(x^(m-1)*tan(a+b*x+c*x^2),x) + (m-1)/(2*c)*Antiderivative(x^(m-2)*log(cos(a+b*x+c*x^2)),x) <--  FreeQ([a,b,c],x) && GtQ(m,1) =#
    #= Antiderivative((~x)^(~m)*cot((~a')+(~b')*(~x)+(~c')*(~x)^2),~x) :=  x^(m-1)*log(sin(a+b*x+c*x^2))/(2*c) - b/(2*c)*Antiderivative(x^(m-1)*cot(a+b*x+c*x^2),x) - (m-1)/(2*c)*Antiderivative(x^(m-2)*log(sin(a+b*x+c*x^2)),x) <--  FreeQ([a,b,c],x) && GtQ(m,1)=#
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * tan(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~(n'), ~x) => Unintegrable((~d + ~e * ~x) ^ ~m * tan(~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x)
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * cot(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~(n'), ~x) => Unintegrable((~d + ~e * ~x) ^ ~m * cot(~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x)
end

max_stars_repo_path: src/rules/4 Trig functions/4.3 Tangent/4.3.12 (d+e x)^m tan(a+b x+c x^2)^n.jl


text:
"""
Package de fonction pour la recherche linéaire d'un MPCCActif

liste des fonctions :
Armijo(ma::ActifMPCCmod.MPCC_actif,xj::Any,d::Any,hg::Any,gradft::Any,stepmax::Float64,slope::Float64;
                verbose :: Bool=false, kwargs...)
ArmijoWolfe(ma::ActifMPCCmod.MPCC_actif,xj::Any,d::Any,hg::Any,gradft::Any,stepmax::Float64,
                slope::Float64;verbose :: Bool=false, kwargs...)
"""
module LineSearch

using ActifMPCCmod

"""
Armijo : Backtracking line search
         1D Minimization
"""
function Armijo(ma::ActifMPCCmod.MPCC_actif,xj::Any,d::Any,hg::Any,stepmax::Float64,slope::Float64,old_alpha::Float64;
                verbose :: Bool=false, kwargs...)

 good_grad=false
 nbW=0
 nbk=0

 step=min(stepmax,1.0)

 ht=ActifMPCCmod.obj(ma,xj+step*d)

 #critère d'Armijo : f(x+alpha*d)-f(x)<=tau_a*alpha*grad f^Td
 while nbk<ma.paramset.ite_max_armijo && ht-hg>ma.paramset.tau_armijo*step*slope
  step*=ma.paramset.armijo_update #step=step_0*(1/2)^m
  ht=ActifMPCCmod.obj(ma,xj+step*d)
  nbk+=1
 end
 step=step

 return step,good_grad,ht,nbk,nbW,gradft
end

"""
ArmijoWolfe : Backtracking line search + amélioration si le pas initial peut être augmenté
         1D Minimization
"""
function ArmijoWolfe(ma::ActifMPCCmod.MPCC_actif,xj::AbstractVector,d::AbstractVector,hg::Float64,stepmax::Float64,slope::Float64,old_alpha::Float64;
                verbose :: Bool=true, kwargs...)

 good_grad=false
 nbW=0
 nbk=0
 step=min(stepmax,1.0)

 #First try to increase t to satisfy loose Wolfe condition
 xjp=xj+step*d
 ht=ActifMPCCmod.obj(ma,xjp)

 gradft=ActifMPCCmod.grad(ma,xjp)
 slope_t=dot(d,gradft)

 while (slope_t<ma.paramset.tau_wolfe*slope) && (ht-hg<=ma.paramset.tau_armijo*step*slope) && (nbW<ma.paramset.ite_max_wolfe) && step<stepmax

  step=min(step*ma.paramset.wolfe_update,stepmax) #on ne peut pas dépasser le pas maximal
  xjp=xj+step*d
  ht=ActifMPCCmod.obj(ma,xjp)
  gradft=ActifMPCCmod.grad(ma,xjp)
  slope_t=dot(d,gradft)

  nbW+=1
 end

 #In this case the problem is very likely to be unbounded
 if nbW==ma.paramset.ite_max_wolfe
@show step stepmax
  return step,good_grad,ht,nbk,nbW,gradft
 end

 #critère d'Armijo : f(x+alpha*d)-f(x)<=tau_a*alpha*grad f^Td
 while nbk<ma.paramset.ite_max_armijo && ht-hg>ma.paramset.tau_armijo*step*slope
  step*=ma.paramset.armijo_update #step=step_0*(1/2)^m
  xjp=xj+step*d
  ht=ActifMPCCmod.obj(ma,xjp)

  nbk+=1
 end

 if nbk==0
  good_grad=true
 end

 return step,good_grad,ht,nbk,nbW,gradft
end

"""
ArmijoWolfe : 'Newarmijo_wolfe' de JPD
Problème avec Hager et Zhang numerical trick
+ Ca n'a pas passé mon test... même sans Hager et Zhang trick...
"""
function ArmijoWolfeHZ(ma::ActifMPCCmod.MPCC_actif,xj::Any,d::Any,hg::Any,stepmax::Float64,
                slope::Float64,old_alpha::Float64;verbose :: Bool=false, kwargs...)

 # Perform improved Armijo linesearch.
 nbk = 0
 nbW = 0
 step = min(stepmax,1.0)

 #First try to increase t to satisfy loose Wolfe condition
 ht=ActifMPCCmod.obj(ma,xj+step*d)
 gradft=ActifMPCCmod.grad(ma,xj+step*d)
 slope_t=dot(d,gradft)
 while (slope_t<ma.paramset.tau_wolfe*slope) && (ht-hg<=ma.paramset.tau_armijo*step*slope) && (nbW<ma.paramset.ite_max_armijo) && step<stepmax
  step=min(step*ma.paramset.wolfe_update,stepmax) #on ne peut pas dépasser le pas maximal
  ht=ActifMPCCmod.obj(ma,xj+step*d)
  slope_t=dot(d,ActifMPCCmod.grad(ma,xj+step*d))

  nbW+=1
 end

# Hager & Zhang numerical trick
 hgoal = hg+slope*step*ma.paramset.tau_armijo
 fact=-0.8
 prec=1e-10
 Armijo = (ht <= hgoal) || ((ht <= hg + prec*abs(hg)) && (slope_t <= fact * slope))
 good_grad=true
 #critère d'Armijo : f(x+alpha*d)-f(x)<=tau_a*alpha*grad f^Td
 while nbk<ma.paramset.ite_max_armijo && !Armijo
  step*=ma.paramset.armijo_update #step=step_0*(1/2)^m
  ht=ActifMPCCmod.obj(ma,xj+step*d)
  hgoal = hg+slope*step*ma.paramset.tau_armijo

  Armijo=false
  good_grad=false
  if ht <= hgoal
   Armijo=true
  elseif ht <= hg+prec*abs(hg)
   gradft=ActifMPCCmod.grad(ma,xj+step*d)
   slope_t = dot(d,gradft)
   good_grad=true
   if slope_t <= fact*slope
    Armijo = true
   end
  end

  nbk+=1
 end

 good_grad=false #à enlever ?
 return step,good_grad,ht,nbk,nbW,gradft
end

#end of module
end


text:
using Interpolations
using Makie
using GeometryBasics
import GeometryBasics: Ngon
using StaticArrays
using FractalTools

# Construct interpolation data
f(x, y) = x^2 + y^2
rect = Rect2D(BigFloat(-1.), BigFloat(-1.), BigFloat(2.), BigFloat(2.))
msh2 = GeometryBasics.mesh(Tesselation(rect, (10, 10)))
pnts3d = [Point(pnt[1], pnt[2], f(pnt[1], pnt[2])) for pnt in msh2.position]

# Construct interpolant
fitp = FractalTools.interpolate(pnts3d)
x = getindex.(pnts3d, 1) |> unique
y = getindex.(pnts3d, 2) |> unique
z = getindex.(pnts3d, 3) |> val -> reshape(val, 10, 10) |> collect
iitp = Interpolations.interpolate((x, y), z, Gridded(Linear()))


text:
function StructDefinition{false}(spec::SpecUnion)
    name = struct_name(spec, spec.is_returnedonly)
    p = Dict(
        :category => :struct,
        :decl => :($name <: $(spec.is_returnedonly ? :ReturnedOnly : :(VulkanStruct{false}))),
        :fields => [
            :($(spec.is_returnedonly ? :data : :vks)::$(spec.name)),
        ],
    )
    StructDefinition{false}(spec, p)
end

function StructDefinition{true}(def::StructDefinition{false,SpecUnion})
    spec = def.spec
    name = struct_name(spec, true)
    p = Dict(
        :category => :struct,
        :decl => :($name <: HighLevelStruct),
        :fields => [
            :(data::$(spec.name)),
        ],
    )
    StructDefinition{true}(spec, p)
end

function constructors(def::StructDefinition{HL,SpecUnion}) where {HL}
    spec = def.spec
    is_high_level = HL || spec.is_returnedonly
    name = struct_name(spec, is_high_level)
    supertypes = supertype_union.(spec.types, is_high_level)
    sig_types = map(zip(spec.types, supertypes)) do (type, supertype)
        if count(==(supertype), supertypes) > 1
            type
        else
            @match supertype begin
                :Unsigned || :Signed => count(in((:Unsigned,:Signed)), supertypes) == 1 ? :Integer : supertype
                _ => supertype
            end
        end
    end

    map(zip(spec.types, sig_types, spec.fields)) do (type, sig_type, field)
        var = wrap_identifier(field)
        call = type in spec_unions.name ? :($var.data) : var
        p = Dict(
            :category => :function,
            :name => name,
            :args => [:($var::$sig_type)],
            :body => :(
                $name($(spec.name)(($call)))
            ),
            :short => true,
        )
        Constructor(def, p)
    end
end

function Constructor(T::StructDefinition{false, SpecUnion}, x::StructDefinition{true, SpecUnion})
    p = Dict(
        :category => :function,
        :name => name(T),
        :args => [:(x::$(name(x)))],
        :body => :($(name(T))(getfield(x, :data))),
        :short => true,
    )
    Constructor(T, p)
end

function GetProperty(def::StructDefinition{true, SpecUnion})
    spec = def.spec
    itr = map(spec.fields) do field
        newfield = wrap_identifier(field)
        (:(sym === $(QuoteNode(newfield))), :(x.data.$field))
    end

    first, rest = Iterators.peel(itr)
    body = Expr(:block)
    ex = Expr(:if, first...)
    push!(body.args, ex)
    for (cond, branch) in rest
        _ex = Expr(:elseif, cond, branch)
        push!(ex.args, _ex)
        ex = _ex
    end
    push!(ex.args, :(getfield(x, sym)))

    p = Dict(
        :category => :function,
        :name => :(Base.getproperty),
        :args => [
            :(x::$(name(def))),
            :(sym::Symbol),
        ],
        :body => body,
        :short => false,
    )
    GetProperty(def, p)
end

function supertype_union(type, is_high_level)
    @match type begin
        :UInt8 || :UInt16 || :UInt32 || :UInt64 => :Unsigned
        :Int8 || :Int16 || :Int32 || :Int64 => :Signed
        :Float8 || :Float16 || :Float32 || :Float64 => :AbstractFloat
        _ => is_high_level ? hl_type(type) : idiomatic_julia_type(type)
    end
end