text:
include("precond.jl")
include("stepsize.jl")
include("stan.jl")

######################
### Mutable states ###
######################

mutable struct ThreePhaseState{T<:Integer}
    n           :: T
    window_size :: T
    next_window :: T
end

################
### Adapterers ###
################

abstract type CompositeAdapter <: AbstractAdapter end

struct NaiveCompAdapter <: CompositeAdapter
    pc  :: PreConditioner
    ssa :: StepSizeAdapter
end

function getss(tp::CompositeAdapter)
    return getss(tp.ssa)
end

# Acknowledgement: this adaption settings is mimicing Stan's 3-phase adaptation.
struct ThreePhaseAdapter{T<:Integer} <: CompositeAdapter
    n_adapts    :: T
    init_buffer :: T
    term_buffer :: T
    pc          :: PreConditioner
    ssa         :: StepSizeAdapter
    state       :: ThreePhaseState{T}
end

function get_threephase_params(::Nothing)
    init_buffer, term_buffer, window_size = 75, 50, 25
    next_window = init_buffer + window_size - 1
    return init_buffer, term_buffer, window_size, next_window
end

function ThreePhaseAdapter(spl::Sampler{<:AdaptiveHamiltonian}, ϵ::Real, dim::Integer)
    # Diagonal pre-conditioner
    # pc = UnitPreConditioner()
    pc = DiagPreConditioner(dim)
    # pc = DensePreConditioner(dim)
    # Dual averaging for step size
    ssa = DualAveraging(spl, spl.info[:adapt_conf], ϵ)
    # Window parameters
    init_buffer, term_buffer, window_size, next_window = get_threephase_params(spl.info[:adapt_conf])
    threephasestate = ThreePhaseState(0, window_size, next_window)
    return ThreePhaseAdapter(spl.alg.n_adapts, init_buffer, term_buffer, pc, ssa, threephasestate)
end

# Ref: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/windowed_adaptation.hpp
function in_adaptation(tp::ThreePhaseAdapter)
    return (tp.state.n >= tp.init_buffer) &&
           (tp.state.n < tp.n_adapts - tp.term_buffer) &&
           (tp.state.n != tp.n_adapts)
end

function is_windowend(tp::ThreePhaseAdapter)
    return (tp.state.n == tp.state.next_window) &&
           (tp.state.n != tp.n_adapts)
end

function compute_next_window!(tp::ThreePhaseAdapter)
    if ~(tp.state.next_window == tp.n_adapts - tp.term_buffer - 1)
        tp.state.window_size *= 2
        tp.state.next_window = tp.state.n + tp.state.window_size
        if ~(tp.state.next_window == tp.n_adapts - tp.term_buffer - 1)
            next_window_boundary = tp.state.next_window + 2 * tp.state.window_size
            if (next_window_boundary >= tp.n_adapts - tp.term_buffer)
                tp.state.next_window = tp.n_adapts - tp.term_buffer - 1
            end
        end
    end
end

function adapt!(tp::ThreePhaseAdapter, stats::Real, θ; adapt_ϵ=false, adapt_M=false)
    if tp.state.n < tp.n_adapts
        tp.state.n += 1
        if tp.state.n == tp.n_adapts
            if adapt_ϵ
                tp.ssa.state.ϵ = exp(tp.ssa.state.x_bar)
            end
            @info " Adapted ϵ = $(getss(tp)), std = $(string(tp.pc)); $(tp.state.n) iterations is used for adaption."
        else
            if adapt_ϵ
                is_updateμ = is_windowend(tp)# || tp.state.n == tp.n_adapts
                adapt!(tp.ssa, stats, is_updateμ)
            end

            # Ref: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/hmc/nuts/adapt_diag_e_nuts.hpp
            if adapt_M
                is_addsample, is_updatestd = in_adaptation(tp), is_windowend(tp)
                adapt!(tp.pc, θ, is_addsample, is_updatestd)
            end

            # If window ends, compute next window
            is_windowend(tp) && compute_next_window!(tp)
        end
    end
end


text:
#: Utils

export conjugate, proxy_objfun, proj_objfun
export no_constraints, indicator

# Scaled version of proximable/projectionable functions

struct ScaledProximableFun{T,N}<:ProximableFunction{T,N}
    c::Real
    g::ProximableFunction{T,N}
end

proxy!(y::AbstractArray{CT,N}, λ::T, g::ScaledProximableFun{CT,N}, x::AbstractArray{CT,N}) where {T<:Real,N,CT<:RealOrComplex{T}} = proxy!(y, λ*g.c, g.g, x)
project!(y::AbstractArray{CT,N}, ε::T, g::ScaledProximableFun{CT,N}, x::AbstractArray{CT,N}) where {T<:Real,N,CT<:RealOrComplex{T}} = project!(y, g.c/ε, g.g, x)

# LinearAlgebra

Base.:*(c::T, g::ProximableFunction{CT,N}) where {T<:Real,N,CT<:RealOrComplex{T}} = ScaledProximableFun{CT,N}(c, g)
Base.:*(c::T, g::ScaledProximableFun{CT,N}) where {T<:Real,N,CT<:RealOrComplex{T}} = ScaledProximableFun{CT,N}(c*g.c, g.g)
Base.:/(g::ProximableFunction{CT,N}, c::T) where {T<:Real,N,CT<:RealOrComplex{T}} = ScaledProximableFun{CT,N}(T(1)/c, g)
Base.:/(g::ScaledProximableFun{CT,N}, c::T) where {T<:Real,N,CT<:RealOrComplex{T}} = ScaledProximableFun{CT,N}(g.c/c, g.g)

# Conjugation of proximable functions

struct ConjugateProxFun{T,N}<:ProximableFunction{T,N}
    g::ProximableFunction{T,N}
end

function proxy!(y::AbstractArray{CT,N}, λ::T, g::ConjugateProxFun{CT,N}, x::AbstractArray{CT,N}) where {T<:Real,N,CT<:RealOrComplex{T}}
    proxy!(y/λ, T(1)/λ, g.g, x)
    x .= y-λ*x
    return x
end

conjugate(g::ProximableFunction{T,N}) where {T,N} = ConjugateProxFun{T,N}(g)
conjugate(g::ConjugateProxFun) = g.g

# Constraint sets

## No constraints

struct NoConstraints{T,N}<:ProjectionableSet{T,N} end

no_constraints(T::DataType, N::Int64) = NoConstraints{T,N}()

project!(x::AbstractArray{T,N}, ::NoConstraints{T,N}, y::Array{T,N}) where {T,N} = (y .= x)

## Sub-level sets with proximable functions

"""
Constraint set C = {x:g(x)<=ε}
"""
struct SubLevelSet{T,N}<:ProjectionableSet{T,N}
    g::ProximableFunction{T,N}
    ε::Real
end

Base.:≤(g::ProximableFunction{CT,N}, ε::T) where {T<:Real,N,CT<:RealOrComplex{T}} = SubLevelSet{CT,N}(g, ε)

project!(x::AbstractArray{T,N}, C::SubLevelSet{T,N}, y::AbstractArray{T,N}) where {T,N} = project!(x, C.ε, C.g, y)
project!(x::AbstractArray{T,N}, C::SubLevelSet{T,N}, y::AbstractArray{T,N}, opt::OptimOptions) where {T,N} = project!(x, C.ε, C.g, y, opt)

## Indicator function

"""
Indicator function δ_C(x) = {0, if x ∈ C; ∞, otherwise} for convex sets C
"""
struct IndicatorFunction{T,N}<:ProximableFunction{T,N}
    C::ProjectionableSet{T,N}
end

indicator(C::ProjectionableSet{T,N}) where {T,N} = IndicatorFunction{T,N}(C)

proxy!(y::AbstractArray{CT,N}, ::T, δ::IndicatorFunction{CT,N}, x::AbstractArray{CT,N}) where {T<:Real,N,CT<:RealOrComplex{T}} = project!(y, δ.C, x)
proxy!(y::AbstractArray{CT,N}, ::T, δ::IndicatorFunction{CT,N}, x::AbstractArray{CT,N}, opt::OptimOptions) where {T<:Real,N,CT<:RealOrComplex{T}} = project!(y, δ.C, x, opt)

# Proximable + linear operator

struct WeightedProximableFun{T,N1,N2}<:ProximableFunction{T,N1}
    g::ProximableFunction{T,N2}
    A::AbstractLinearOperator{T,N1,N2}
end

Base.:∘(g::ProximableFunction{T,N2}, A::AbstractLinearOperator{T,N1,N2}) where {T,N1,N2} = WeightedProximableFun{T,N1,N2}(g, A)

function proxy!(y::AbstractArray{CT,N1}, λ::T, g::WeightedProximableFun{CT,N1,N2}, x::AbstractArray{CT,N1}, opt::OptimOptions) where {T<:Real,N1,N2,CT<:RealOrComplex{T}}

    # Objective function (dual problem)
    f = leastsquares_misfit(adjoint(g.A), y/λ)+conjugate(g.g)/λ

    # Minimization (dual variable)
    p0 = similar(y, range_size(g.A)); p0 .= T(0)
    p = minimize(f, p0, opt)

    # Dual to primal solution
    return x .= y-λ*(adjoint(g.A)*p)

end

function project!(y::AbstractArray{CT,N1}, ε::T, g::WeightedProximableFun{CT,N1,N2}, x::AbstractArray{CT,N1}, opt::OptimOptions) where {T<:Real,N1,N2,CT<:RealOrComplex{T}}

    # Objective function (dual problem)
    f = leastsquares_misfit(adjoint(g.A), y)+conjugate(indicator(g.g ≤ ε))

    # Minimization (dual variable)
    p0 = similar(y, range_size(g.A)); p0 .= T(0)
    p = minimize(f, p0, opt)

    # Dual to primal solution
    return x .= y-adjoint(g.A)*p

end

(g::WeightedProximableFun{T,N1,N2})(x::AbstractArray{T,N1}) where {T,N1,N2} = g.g(g.A*x)

Flux.gpu(g::WeightedProximableFun{T,N1,N2}) where {T,N1,N2} = WeightedProximableFun{T,N1,N2}(g.g, gpu(g.A))
Flux.cpu(g::WeightedProximableFun{T,N1,N2}) where {T,N1,N2} = WeightedProximableFun{T,N1,N2}(g.g, cpu(g.A))

# Proximable function evaluation

struct ProxyObjFun{T,N} <: DifferentiableFunction{T,N}
    λ::Real
    g::ProximableFunction{T,N}
    opt::Union{Nothing,OptimOptions}
end

proxy_objfun(λ::T, g::ProximableFunction{CT,N}) where {T<:Real,N,CT<:RealOrComplex{T}} = ProxyObjFun{CT,N}(λ, g, nothing)
proxy_objfun(λ::T, g::ProximableFunction{CT,N}, opt::OptimOptions) where {T<:Real,N,CT<:RealOrComplex{T}} = ProxyObjFun{CT,N}(λ, g, opt)

function grad!(f::ProxyObjFun{T,N}, y::AbstractArray{T,N}, g::Union{Nothing,AbstractArray{T,N}}) where {T,N}
    f.opt === nothing ? (x = proxy(y, f.λ, f.g)) : (x = proxy(y, f.λ, f.g, f.opt))
    g !== nothing && (g .= y-x)
    return T(0.5)*norm(x-y)^2+f.λ*f.g(x)
end

struct ProjObjFun{T,N} <: DifferentiableFunction{T,N}
    ε::Real
    g::ProximableFunction{T,N}
    opt::Union{Nothing,OptimOptions}
end

proj_objfun(ε::T, g::ProximableFunction{CT,N}) where {T<:Real,N,CT<:RealOrComplex{T}} = ProjObjFun{CT,N}(ε, g, nothing)
proj_objfun(ε::T, g::ProximableFunction{CT,N}, opt::OptimOptions) where {T<:Real,N,CT<:RealOrComplex{T}} = ProjObjFun{CT,N}(ε, g, opt)

function grad!(f::ProjObjFun{T,N}, y::AbstractArray{T,N}, g::Union{Nothing,AbstractArray{T,N}}) where {T,N}
    f.opt === nothing ? (x = project(y, f.ε, f.g)) : (x = project(y, f.ε, f.g, f.opt))
    g !== nothing && (g .= y-x)
    return T(0.5)*norm(x-y)^2
end

# Minimizable type utils

struct DiffPlusProxFun{T,N}<:MinimizableFunction{T,N}
    f::DifferentiableFunction{T,N}
    g::ProximableFunction{T,N}
end

Base.:+(f::DifferentiableFunction{T,N}, g::ProximableFunction{T,N}) where {T,N} = DiffPlusProxFun{T,N}(f, g)

ext: jl

max_stars_repo_name: grizzuti/FastSolversForWeightedTV

max_issues_repo_name: grizzuti/FastSolversForWeightedTV

max_forks_repo_name: grizzuti/FastSolversForWeightedTV


text:
##
# # A simple two-dimensional optimization problem

# ## Loading ClimateMARGO.jl
using ClimateMARGO # Julia implementation of the MARGO model
using PyPlot # A basic plotting package

# Loading submodules for convenience
using ClimateMARGO.Models
using ClimateMARGO.Utils
using ClimateMARGO.Diagnostics
using ClimateMARGO.Optimization

# ## Loading the default MARGO configuration
params = deepcopy(ClimateMARGO.IO.included_configurations["default"])

# Slightly increasing the discount rate to 3% to be more comparable with other models
params.economics.ρ = 0.03

# ## Reducing the default problem's dimensionality from ``4N`` to ``2``.

# The default optimization problem consists of optimizing the values of each of the 4 controls for N timesteps, giving a total problem size of:
4*length(t(params))

# ### Modifying `ClimateModelParameters`
# Thanks to MARGO's flexibility, we can get rid of two of the controls and set the other two to have constant values, effectively reducing the problem size to 2.

# First, we have to remove initial conditions on the control variables, which would conflict with the constant-value constraint
params.economics.mitigate_init = nothing
params.economics.remove_init = nothing
params.economics.geoeng_init = nothing
params.economics.adapt_init = nothing

# ### Instanciating the `ClimateModel`
# Now that we've finished changing parameter values, we can create our MARGO model instance:
m = ClimateModel(params)

# ### Modifying keyword arguments for `optimize_controls!`
# We can make the controls constant in time by asserting a maximum deployment rate of zero
max_slope = Dict("mitigate"=>0., "remove"=>0., "geoeng"=>0., "adapt"=>0.);

# Since we make the controls constant in time, we have to let them turn on immediately by removing any deployment delays
delay_deployment = Dict("mitigate"=>0., "remove"=>0., "geoeng"=>0., "adapt"=>0.);

# Now we get rid of geoengineering and adaptation options by setting their maximum deployment fraction to zero
max_deployment = Dict("mitigate"=>1.0, "remove"=>1.0, "geoeng"=>0., "adapt"=>0.);

# ### Run the optimization problem once to test
@time optimize_controls!(m, obj_option = "net_benefit", max_slope=max_slope, max_deployment=Dict("mitigate"=>1.0, "remove"=>0., "geoeng"=>0., "adapt"=>0.), delay_deployment=delay_deployment);

# ### Visualizing the results
fig, axes = ClimateMARGO.Plotting.plot_state(m);
axes[3].legend(loc="upper left")
axes[4].legend(loc="lower left")
axes[5].legend(loc="lower right")
axes[6].legend(loc="upper right")
axes[6].set_ylim(0,2.5)
savefig("figures/1D_optimization_state.png", dpi=150, bbox_inches="tight")
gcf()

##
# ## Comparing the two-dimensional optimization with the brute-force parameter sweep method

# ### Parameter sweep

# In the brute-force approach, we sweep through all possible values of **M**itigation and **R**emoval to map out the objective functions and constraints and visually identify the "optimal" solution in this 2-D space.
Ms = 0.:0.005:1.0;
Rs = 0.:0.005:1.0;

# We will also consider four different temperature thresholds and visualize these constraints in the 2-D space
temp_goals = [1.5, 2.0, 3., 4.0]

control_cost = zeros(length(Rs), length(Ms)) .+ 0. #
net_benefit = zeros(length(Rs), length(Ms)) .+ 0. #
max_temp = zeros(length(Rs), length(Ms)) .+ 0. # stores the maximum temperature acheived for each combination
min_temp = zeros(length(Rs), length(Ms)) .+ 0. # stores the minimum temperature acheived for each combination
optimal_controls = zeros(2, length(temp_goals), 2) # to hold optimal values computed using JuMP

for (o, option) = enumerate(["temp", "net_benefit"])
    for (i, M) = enumerate(Ms)
        for (j, R) = enumerate(Rs)
            m.controls.mitigate = zeros(size(t(m))) .+ M
            m.controls.remove = zeros(size(t(m))) .+ R
            if minimum(c(m, M=true, R=true)) <= 100.
                continue
            end
            control_cost[j, i] = net_present_cost(m, discounting=true, M=true, R=true)
            net_benefit[j, i] = net_present_benefit(m, discounting=true, M=true, R=true)
            max_temp[j, i] = maximum(T(m, M=true, R=true))
            min_temp[j, i] = minimum(T(m, M=true, R=true))
        end
    end
    for (g, temp_goal) = enumerate(temp_goals)
        optimize_controls!(m, obj_option = option, temp_goal = temp_goal, max_slope=max_slope, max_deployment=max_deployment, delay_deployment=delay_deployment);
        optimal_controls[:, g, o] = [m.controls.mitigate[1], m.controls.remove[1]]
    end
end

##
# ### Visualizing the one-dimensional mitigation optimization problem
# In the limit of zero-carbon dioxide removal, we can recover the 1D mitigation optimization problem from the 2D one.
col = ((1., 0.8, 0.), (0.8, 0.5, 0.), (0.7, 0.2, 0.), (0.6, 0., 0.),)

fig = figure(figsize=(14,5))
ax = subplot(1,2,1)
ind1 = argmin(abs.(max_temp[1,:] .- 1.5))
ind2 = argmin(abs.(max_temp[1,:] .- 2.))
ind3 = argmin(abs.(max_temp[1,:] .- 3.))
ind4 = argmin(abs.(max_temp[1,:] .- 4.))
plot(Ms, control_cost[1,:], "k-", lw=2)
plot(Ms[ind1], control_cost[1,ind1], "o", color=col[1], markersize=10)
plot(Ms[ind2], control_cost[1,ind2], "o", color=col[2], markersize=10)
plot(Ms[ind3], control_cost[1,ind3], "o", color=col[3], markersize=10)
plot(Ms[ind4], control_cost[1,ind4], "o", color=col[4], markersize=10)
minM1 = Ms[ind1]
minM2 = Ms[ind2]
minM3 = Ms[ind3]
minM4 = Ms[ind4]
ylims = ax.get_ylim()
fill_between([minM2,minM1], [ylims[1], ylims[1]], [ylims[2],ylims[2]], color=col[1], alpha=0.2)
fill_between([minM3,minM2], [ylims[1], ylims[1]], [ylims[2],ylims[2]], color=col[2], alpha=0.2)
fill_between([minM4,minM3], [ylims[1], ylims[1]], [ylims[2],ylims[2]], color=col[3], alpha=0.2)
fill_between([0,minM4], [ylims[1], ylims[1]], [ylims[2],ylims[2]], color=col[4], alpha=0.2)
axvline(minM1, color=col[1])
axvline(minM2, color=col[2])
axvline(minM3, color=col[3])
axvline(minM4, color=col[4])
ylim(ylims)
xlim(0,1)
xlabel("Emissions mitigation level [% reduction]")
xticks(0.:0.2:1.0, ["0%", "20%", "40%", "60%", "80%", "100%"])
ylabel("Net present cost of controls [trillion USD]")
plot([], [], "o", color="grey", label="lowest cost", markersize=10.)
legend(loc="upper left")

ax = subplot(1,2,2)
plot(Ms, net_benefit[1,:], "k-", lw=2)
ind = argmax(net_benefit[1,:])
plot(Ms[ind], net_benefit[1,ind], "ko", markersize=10, label="most benefits")
ylims = ax.get_ylim()
axvline(minM1, color=col[1])
axvline(minM2, color=col[2])
axvline(minM3, color=col[3])
axvline(minM4, color=col[4])
ylim(ylims)
xlim(0,1)
xlabel("Emissions mitigation level [% reduction]")
xticks(0.:0.2:1.0, ["0%", "20%", "40%", "60%", "80%", "100%"])
ylabel("Net present benefits, relative to baseline [trillion USD]")
legend(loc="upper left")
savefig("figures/1D_optimization.png", dpi=150, bbox_inches="tight")
gcf()

##
# ### Visualizing the two-dimensional optimization problem

figure(figsize=(14, 5))

o = 1
subplot(1,2,o)
pcolor(Ms, Rs, control_cost, cmap="Greys", vmin=0., vmax=150.)
cbar = colorbar(label="Net present cost of controls [trillion USD]")
contour(Ms, Rs, control_cost, levels=[25, 50, 75], colors="k", linewidths=0.85, alpha=0.4)

grid(true, color="k", alpha=0.25)
temp_mask = ones(size(max_temp))
temp_mask[max_temp .<= 4.0] .= NaN
contourf(Ms, Rs, temp_mask, cmap="Reds", alpha=1.0, vmin=0., vmax=2.)

temp_mask = ones(size(min_temp))
temp_mask[min_temp .> 0.] .= NaN
contourf(Ms, Rs, temp_mask, cmap="Blues", alpha=1.0, vmin=0., vmax=2.5)
contour(Ms, Rs, min_temp, colors="darkblue", levels=[0.], linewidths=2.5)
plot([], [], "o", color="grey", label="lowest cost", markersize=10.)
plot([], [], "-", color="darkblue", label=latexstring("\$\\min(T)=\$", 0), lw=2.5)

for (g, temp_goal) = enumerate(temp_goals)
    contour(Ms, Rs, max_temp, colors=[col[g]], levels=[temp_goal], linewidths=2.5)
    plot(optimal_controls[1,g,o], optimal_controls[2,g,o], "o", color=col[g], markersize=10.)
    plot([], [], "-", color=col[g], label=latexstring("\$\\max(T)=\$", temp_goal), lw=2.5)
end
legend(loc="upper left")

xlabel("Emissions mitigation level [% reduction]")
xticks(0.:0.2:1.0, ["0%", "20%", "40%", "60%", "80%", "100%"])
ylabel(L"CO$_{2e}$ removal rate [% of present-day emissions]")
yticks(0.:0.2:1.0, ["0%", "20%", "40%", "60%", "80%", "100%"])
annotate(L"$T > 4\degree$C", (0.02, 0.04), xycoords="axes fraction", color="darkred", fontsize=13)
annotate(L"$T < 0\degree$C", (0.74, 0.66), xycoords="axes fraction", color="darkblue", fontsize=13)
title("Cost-effectiveness analysis")

o = 2
subplot(1,2,o)
q = pcolor(Ms, Rs, net_benefit, cmap="Greys_r", vmin=-75, vmax=75)
cbar = colorbar(label="Net present benefits, relative to baseline [trillion USD]", extend="both")
contour(Ms, Rs, net_benefit, levels=[0, 25, 50], colors="k", linewidths=0.85, alpha=0.4)

grid(true, color="k", alpha=0.25)
temp_mask = ones(size(min_temp))
temp_mask[(min_temp .> 0.)] .= NaN
contourf(Ms, Rs, temp_mask, cmap="Blues", alpha=1.0, vmin=0., vmax=2.5)
contour(Ms, Rs, min_temp, colors="darkblue", levels=[0.], linewidths=2.5)
plot([], [], "-", color="darkblue")

plot(optimal_controls[1,1,2], optimal_controls[2,1,2], "o", color="k", markersize=10., label="most benefits")
for (g, temp_goal) = enumerate(temp_goals)
    contour(Ms, Rs, max_temp, colors=[col[g]], levels=[temp_goal], linewidths=2.5)
end
legend()

xlabel("Emissions mitigation level [% reduction]")
xticks(0.:0.2:1.0, ["0%", "20%", "40%", "60%", "80%", "100%"])
ylabel(L"CO$_{2e}$ removal rate [% of present-day emissions]")
yticks(0.:0.2:1.0, ["0%", "20%", "40%", "60%", "80%", "100%"])
annotate(L"$T < 0\degree$C", (0.74, 0.66), xycoords="axes fraction", color="darkblue", fontsize=13)
title("Cost-benefit analysis")
savefig("figures/2D_optimization.png", dpi=150, bbox_inches="tight")
gcf()


text:
# This file is a part of BitArraynge.jl, licensed under the MIT License (MIT).
@assert true


text:
# Our developed libraries
using CombinatorialSpaces
using CombinatorialSpaces.ExteriorCalculus

using Catlab
using Catlab.Present
using Catlab.Graphics
using Catlab.CategoricalAlgebra
using Catlab.Programs
using LinearAlgebra

# Julia community libraries
using MeshIO
using CairoMakie
using DifferentialEquations
using Logging: global_logger
using TerminalLoggers: TerminalLogger
global_logger(TerminalLogger())

using Revise
using Decapodes.Simulations
using Decapodes.Examples
using Decapodes.Diagrams
using Decapodes.Schedules
using Decapodes.Debug
using Decapodes.OpenDiagrams

using CombinatorialSpaces: volume
using SparseArrays

draw_dwd(dwd) = begin
  new_res = deepcopy(dwd)
  new_res.diagram[:outer_out_port_type] .= nothing
  new_res.diagram[:outer_in_port_type] .= nothing
  new_res.diagram[:out_port_type] .= nothing
  new_res.diagram[:in_port_type] .= nothing
  to_graphviz(new_res)
end
draw_equations(eq) = begin
  to_graphviz(eq,
  node_labels=true, prog="neato",
  node_attrs=Dict(:shape=>"oval"),
  graph_attrs=Dict(:nodesep=>"4.0"))
end

######################
# Functions for Curl #
######################
include("special_ops.jl")
include("deca_schema.jl")
include("physics.jl")
include("complex_operators.jl")
include("constants.jl")

###############
# Import Mesh #
###############

s1 = EmbeddedDeltaSet2D("../su2_mesh_square_small_51.stl")
s = EmbeddedDeltaSet2D{Bool, Point{3, Float64}}()
copy_parts!(s, s1)
sd = dual(s);
if ⋆(Val{1}, sd)[1,1] < 0.0
  orient_component!(s, 1, false)
  sd = dual(s);
end

include("boundaries.jl")
include("functions.jl")

c_objs = fill(288.15, nv(s))
c_objs[∂ₒ₀] .= 350.0
#c_objs[∂ₗ₀] .= t2e * 461.04
#c_objs[∂ᵣ₀] .= t2e * 115.26
velocity(p) = [3.402, 0.0, 0.0]
gravity(p) = [0.0,0.0,0.0]
v = ♭(sd, DualVectorField(velocity.(sd[triangle_center(sd),:dual_point]))).data;
g = ♭(sd, DualVectorField(gravity.(sd[triangle_center(sd),:dual_point]))).data;
p = [density for p in s[:point]] * (288.15 * R₀)

@show ⋆(Val{1},sd)[1,1]

exp_dwd = Examples.expand_dwd(new_dwd, rules)
new_dwd = deepcopy(exp_dwd)
Examples.zip_dwd!(new_dwd)
Examples.contract_matrices!(new_dwd, funcs)

func, _ = gen_sim(new_dwd, funcs, sd; autodiff=false, params=[:G]);

tstep = 2.0

using JSON
for i in 1:10
  GC.gc()
  prev_loc = "res_$(i-1)"
  loc = "res_$i"
  mkpath(loc)
  u0 = vcat(c_objs, p, v)
  if i == 1
    pre_cond = open("../cht_sims/final_cond_v2.json", "r") do f
        JSON.parse(f)
    end
    u0[(1:ne(s)) .+ 2*nv(s)] .= pre_cond[(nv(s)+1):(nv(s)+ne(s))]
  else
    pre_cond = open("$prev_loc/sim_res.json", "r") do f
        JSON.parse(f)
    end
    u0 .= pre_cond["$tstep"]
    u0[∂ₑ₀₊] .= c_objs[∂ₑ₀₊]
    u0[∂ₑ₀₊ .+ nv(s)] .= p[∂ₑ₀₊]
    u0[∂ₑ₁₊ .+ (2*nv(s))] .= v[∂ₑ₁₊]
  end
  ##############
  # Simulation #
  ##############
  dt = 0.01
  prob = ODEProblem(func, u0, (0.0, tstep))
  sol1 = solve(prob, Tsit5(), progress=true, progress_steps=1, dtmax=#=4e-5=#8e-4, saveat=dt, p=g)

  #=@show extrema(sol1.t[2:(end-1)] .- sol1.t[1:(end-2)])
  fig, ax, ob = plot(sol1.t[2:(end-1)] .- sol1.t[1:(end-2)])
  save("$(loc)/timesteps.png", fig)=#

  #################
  # Write Results #
  #################
  t_range = 1:nv(s)
  ρ_range = (1:nv(s)) .+ nv(s)
  v_range = (1:ne(s)) .+ (2 * nv(s))

  using JSON
  res = Dict{Float64, Vector}()
  for i in 0.0:dt:sol1.t[end]
      res[i] = sol1(i)
  end
  open("$loc/sim_res.json", "w") do f
      JSON.print(f, res)
  end

  using WriteVTK
  for i in 1:length(dt:dt:sol1.t[end])
    tris = triangle_vertices(s)
    cells = [MeshCell(VTKCellTypes.VTK_TRIANGLE, [tris[1][t],tris[2][t],tris[3][t] ]) for t in 1:ntriangles(s)];
    x =  [p[1] for p in s[:point]]
    y =  [p[2] for p in s[:point]]
    z =  [p[3] for p in s[:point]]
    vtkfile = vtk_grid("$(loc)/su2_coupled_v1_$(lpad(i, 4, "0"))", x, y, z, cells)
    vel_form = ♯(sd, sol1(dt*i)[v_range], CombinatorialSpaces.DiscreteExteriorCalculus.PPSharp())

    vtkfile["temperature", VTKPointData()] = e2t * sol1(dt*i)[t_range]
    vtkfile["pressure", VTKPointData()] = sol1(dt*i)[ρ_range]
    vtkfile["vel", VTKPointData()] = vel_form
    vtkfile["v_mag", VTKPointData()] = sqrt.(abs.(1e-7 .+ inv_hodge_star(Val{0}, sd)*pd_wedge(Val{(1,1)}, sd, sol1(dt*i)[v_range], ⋆(Val{1}, sd) * sol1(dt*i)[v_range])))
    vtk_save(vtkfile)
  end

  fig = Figure()
  ax, ob = mesh(fig[1,1], s, color = e2t * sol1[end][t_range], colormap=:seismic)#, colorrange=(115 + 100,462 - 100))
  xl, xr = (0.0, 11.0)
  yl, yr = (-1, 1)
  asp = 11 / 2.0
  xlims!(ax, (xl,xr))
  ylims!(ax, (yl,yr))
  ax.aspect = AxisAspect(asp)
  Colorbar(fig[1,2], ob)
  fig
  save("$(loc)/heat_prof.png", fig)

  times = range(0.01, sol1.t[end]/2, length=150)
  colors = [sol1(t)[t_range] for t in times]
  fig = Figure()
  axis, ob = mesh(fig[1,1], s, color=colors[1], colormap=:seismic)
  axis.aspect = AxisAspect(asp)
  Colorbar(fig[1,2], ob)
  framerate = 30
  fig

  xlims!(axis, (xl, xr))
  ylims!(axis, (yl, yr))
  record(fig, "$(loc)/conc_flow.gif", collect(1:length(collect(times))); framerate = framerate) do i
    ob.color = colors[i]
  end

  times = range(0.01, sol1.t[end], length=150)
  colors = [sol1(t)[ρ_range] for t in times]

  fig = Figure()
  axis, ob = mesh(fig[1,1], s, color=colors[1],
                                     colorrange=(minimum(vcat(colors...)),maximum(vcat(colors...))))

  xlims!(axis, (xl, xr))
  ylims!(axis, (yl, yr))
  axis.aspect = AxisAspect(asp)
  Colorbar(fig[1,2], ob)
  framerate = 30
  fig

  record(fig, "$(loc)/press_flow.gif", collect(1:length(collect(times))); framerate = framerate) do i
    ob.color = colors[i]
  end

  times = range(0, sol1.t[end], length=150)
  colors = [ sol1(t)[ρ_range] ./ (R₀ * sol1(t)[t_range]) for t in times]

  figure, axis, ob = mesh(s, color=colors[1],
                                     colorrange=(minimum(vcat(colors...)),maximum(vcat(colors...))))

  xlims!(axis, (xl, xr))
  ylims!(axis, (yl, yr))
  axis.aspect = AxisAspect(asp)
  framerate = 30

  record(figure, "$(loc)/density_flow.gif", collect(1:length(collect(times))); framerate = framerate) do i
    ob.color = colors[i]
  end

  magnitudes(t) = sqrt.(abs.(1e-4 .+ inv_hodge_star(Val{0}, sd)*pd_wedge(Val{(1,1)}, sd, sol1(t)[v_range], ⋆(Val{1}, sd) * sol1(t)[v_range])))
  times = range(0.01, sol1.t[end], length=150)
  colors = [magnitudes(t) for t in times]
  fig = Figure()
  axis, ob = mesh(fig[1,1], s, color=colors[3], colorrange=(minimum(vcat(colors...)),maximum(vcat(colors...))))

  xlims!(axis, (xl, xr))
  ylims!(axis, (yl, yr))
  axis.aspect = AxisAspect(asp)
  Colorbar(fig[1,2], ob)
  framerate = 30
  fig
  record(fig, "$(loc)/vel_flow.gif", collect(1:length(collect(times))); framerate = framerate) do i
    ob.color = colors[i]
  end
end

max_stars_repo_name: AlgebraicJulia/DECAPODES-Benchmarks

max_issues_repo_name: AlgebraicJulia/DECAPODES-Benchmarks

max_forks_repo_name: AlgebraicJulia/DECAPODES-Benchmarks


text:
using Random
using POMDPSimulators

@testset "trace_generation.jl" begin
    pp = simple_pp()
    mg = PMDPs.PMDPg(pp)
    h = PMDPs.simulate_trace(mg, MersenneTwister(12))
    @test isa(h, SimHistory)
end


text:
function nearest_neighbor(grd::Grid, Ns)

    _, proj_j, _ = findnz(grd.Proj)
    dK = proj_j[end] - proj_j[1]

    is = Int64[]
    js = Int64[]
    dats = Int64[]

    for j=1:grd._Ny
        for i=1:grd._Nx

            r = grd.r0 + grd.dr*[i-1, j-1]

            k_cart = i + (j-1)*grd._Nx
            if k_cart > proj_j[1]
                k_approx = Int64(ceil((k_cart - proj_j[1]) * grd.Nk / dK))
            else
                k_approx = 1
            end

            kmax = minimum([k_approx+Ns, grd.Nk])
            kmin = maximum([k_approx-Ns, 1])
            dists = [norm(r - grd.points[:,i]) for i=kmin:kmax]

            try
                mindist, _knearest = findmin(dists)
                k_nearest = [kmin:kmax...][_knearest]
                is = [is; k_cart]
                js = [js; k_nearest]
                dats = [dats; 1]
            catch
                println(k_approx)
                println(kmin)
                println(kmax)
            end
        end
    end

    return sparse(is, js, dats, grd._Nx*grd._Ny, grd.Nk)
end


text:
include("test_Dog.jl")
include("test_Dogs.jl")
include("test_Sheeps.jl")
include("test_DenseSheeps.jl")
include("test_Simulation.jl")


text:
using GraphSignals
using CUDA
using FillArrays
using GraphLaplacians
using LightGraphs
using SimpleWeightedGraphs
using Test

tests = [
    "featuredgraph",
    "graph",
    "linalg",
    "simplegraph",
    "weightedgraph",
]

if CUDA.functional()
    push!(tests, "cuda")
end

@testset "GraphSignals.jl" begin
    for t in tests
        include("$(t).jl")
    end
end


text:
###############################################################
# MultiValue Type
###############################################################

"""
Type representing a multi-dimensional value
"""
abstract type MultiValue{S,T,N,L} <: Number end

@inline Base.Tuple(arg::MultiValue) = arg.data

# Custom type printing

function show(io::IO,v::MultiValue)
  print(io,v.data)
end

function show(io::IO,::MIME"text/plain",v:: MultiValue)
  print(io,typeof(v))
  print(io,v.data)
end

###############################################################
# Other constructors and conversions implemented for more generic types
###############################################################

change_eltype(::Type{<:Number},::Type{T}) where {T} = T
change_eltype(::Number,::Type{T2}) where {T2} = change_eltype(Number,T2)

num_components(::Type{<:Number}) = 1
num_components(::Number) = num_components(Number)

function n_components(a)
  msg = "Function n_components has been removed, use num_components instead"
  error(msg)
end

function data_index(::Type{<:MultiValue},i...)
  @abstractmethod
end


text:
abstract type AbstractWorkingData{NG,T} end

struct WorkingData{NG,T} <: AbstractWorkingData{NG,T}
    inverse_wisharts::Vector{InverseWishart{T}}
    individual_probs::ResizableMatrix{T,NG}
    groups::Groups{NG}
end
struct WorkingDataUnifCache{NG,T} <: AbstractWorkingData{NG,T}
    inverse_wisharts::Vector{InverseWishart{T}}
    individual_probs::ResizableMatrix{T,NG}
    uniform_probs::Vector{T}
    groups::Groups{NG}
end
function WorkingData(N, ::Type{T}, ::Val{NG}) where {T,NG}
    WorkingData(
        Vector{InverseWishart{T}}(undef, NG), # inverse_wisharts
        ResizableMatrix{T,NG}(undef, N), # individual_probs
        Groups{NG}(undef, N) # groups
    )
end
function WorkingDataUnifCache(N, ::Type{T}, ::Val{NG}) where {T,NG}
    WorkingData(
        Vector{InverseWishart{T}}(undef, NG), # inverse_wisharts
        ResizableMatrix{T,NG}(undef, N), # individual_probs
        Vector{T}(undef, N), #uniform_probs
        Groups{NG}(undef, N) # groups
    )
end
function Base.resize!(wd::WorkingData{NG,T}, N) where {NG,T}
    W = LoopVectorization.pick_vector_width(T)
    # r = ((N - 1) & (W - 1)) + 1
    # N_cap = N + W - r
    # sizehint!(wd.individual_probs, N + W)
    # sizehint!(wd.groups, N + W)
    resize!(wd.individual_probs, N)
    resize!(wd.groups, N)
    wd
end

function Base.resize!(wd::WorkingDataUnifCache, N)
    resize!(wd.individual_probs, N)
    resize!(wd.uniform_probs, N)
    resize!(wd.groups, N)
    wd
end

function WorkDataChains(N, ::Type{T}, ::Val{NG}, chains) where {T,NG}
    workdatachains = Vector{WorkingData{NG,T}}(undef, chains)
    @threads for i ∈ 1:chains
        workdatachains[i] = WorkingData(N, T, Val(NG))
    end
    workdatachains
end


text:
#!/usr/bin/env julia
# -*- coding: UTF-8 -*-
# __julia-version__ = 1.7.2
# __author__        = "Max Kannenberg"
# __copyright__     = "2022"
# __license__       = "ISC"

#########################
## literature the driving dynamics equations are based on:
##
## @incollection{Bruenger:2014, % Chapter 4
##   author    = {Brünger, Olaf and Dahlhaus, Elias},
##   year      = {2014},
##   title     = {Running Time Estimation},
##   pages     = {65--90},
##   booktitle = {Railway Timetabling \& Operations.},
##   editora   = {Hansen, Ingo A.},
##   editorb   = {Pachl, Jörn},
##   isbn      = {978-3-777-10462-1},
##   publisher = {Eurailpress DVV Media Group},
## }
## @Book{Wende:2003,
##   author    = {Wende, Dietrich},
##   date      = {2003},
##   title     = {Fahrdynamik des Schienenverkehrs},
##   isbn      = {978-3-322-82961-0},
##   publisher = {Springer-Verlag},
## }
#########################

approxLevel = 6
v00 = 100/3.6     # velocity factor (in m/s)

## calculate forces

#TODO: replace the ? ? ?
"""
    calcTractionUnitResistance(v, train)

Calculate the vehicle resistance for the traction unit of the `train` dependend on the velocity `v`.

...
# Arguments
- `v::AbstractFloat`: the current velocity in m/s.
- `train::Train`: ? ? ?
...

# Examples
```julia-repl
julia> calcTractionUnitResistance(30.0, ? ? ?)
? ? ?
```
"""
function calcTractionUnitResistance(v::AbstractFloat, train::Train)
    # equation is based on [Wende:2003, page 151]
    f_Rtd0 = train.f_Rtd0 # coefficient for basic resistance due to the traction units driving axles (in ‰)
    f_Rtc0 = train.f_Rtc0 # coefficient for basic resistance due to the traction units carring axles (in ‰)
    F_Rt2  = train.F_Rt2  # coefficient for air resistance of the traction units (in N)
    m_td   = train.m_td   # mass on the traction unit's driving axles (in kg)
    m_tc   = train.m_tc   # mass on the traction unit's carrying axles (in kg)

    F_R_tractionUnit = f_Rtd0/1000 * m_td * g + f_Rtc0/1000 * m_tc * g + F_Rt2 * ((v + Δv_air) /v00)^2   # vehicle resistance of the traction unit (in N)   # /1000 because of the unit ‰
    # TODO: use calcForceFromCoefficient? F_R_tractionUnit = calcForceFromCoefficient(f_Rtd0, m_td) + calcForceFromCoefficient(f_Rtc0, m_tc) + F_Rt2 * ((v + Δv_air) /v00)^2       # vehicle resistance of the traction unit (in N)
    return F_R_tractionUnit
    #TODO: same variable name like in the rest of the tool? return R_traction
    #TODO: just one line? return train.f_Rtd0/1000*train.m_td*g+train.f_Rtc0/1000*train.m_tc*g+train.F_Rt2*((v+train.Δv_air)/v00)^2    # /1000 because of the unit ‰
end #function calcTractionUnitResistance

"""
TODO
calculate and return the wagons vehicle resistance dependend on the velocity
"""
function calcWagonsResistance(v::AbstractFloat, train::Train)
    # equation is based on a combination of the equations of Strahl and Sauthoff [Wende:2003, page 153] with more detailled factors (Lehmann, page 135)
    f_Rw0  = train.f_Rw0  # coefficient for basic resistance of the set of wagons (consist)  (in ‰)
    f_Rw1  = train.f_Rw1  # coefficient for the consists resistance to rolling (in ‰)
    f_Rw2  = train.f_Rw2  # coefficient fo the consistsr air resistance (in ‰)
    m_w    = train.m_w    # mass of the set of wagons (consist)  (in kg)

    F_R_wagons = m_w *g *(f_Rw0/1000 + f_Rw1/1000 *v /v00 + f_Rw2/1000 * ((v + Δv_air) /v00)^2)     # vehicle resistance of the wagons (in N)      # /1000 because of the unit ‰
# TODO: use calcForceFromCoefficient?    F_R_wagons = calcForceFromCoefficient(f_Rw0, m_w) + calcForceFromCoefficient(f_Rw1, m_w) *v /v00 + calcForceFromCoefficient(f_Rw2, m_w) * ((v + Δv_air) /v00)^2     # vehicle resistance of the wagons (in N)
    return F_R_wagons
end #function calcWagonsResistance

function calcForceFromCoefficient(f_R::Real, m::Real)
    # equation is based on [Wende:2003, page 8]

    # f_R: specific resistance (in ‰)
    # m: vehicle's mass (in kg)

    F_R = f_R /1000 *m *g     # Resisting Force (in N)  # /1000 because of the unit ‰
    return F_R
end #function calcForceFromCoefficient

function calcAcceleration(F_T::Real, F_R::Real, m_train::Real, ξ_train::Real)
    # equation is based on [Bruenger:2014, page 72] with a=dv/dt

    # F_T: tractive effort (in N)
    # F_R: resisting forces (in N)
    # m_train: train's mass (in kg)
    # ξ_train: train's rotation mass factor (without unit)

    a = (F_T - F_R) /m_train /ξ_train      # acceleration (in m/s)
    return a
end #function calcAcceleration

function calc_Δs_with_Δt(Δt::Real, a_prev::Real, v_prev::Real)
    # equation is based on [Wende:2003, page 37]

    # Δt: time step (in s)
    # a_prev: acceleration from previous data point
    # v_prev: velocitiy from previous data point
    Δs = Δt * (2*v_prev + Δt*a_prev) /2        # step size (in m)
    return Δs
end #function calc_Δs_with_Δt

function calc_Δs_with_Δv(Δv::Real, a_prev::Real, v_prev::Real)
    # equation is based on [Wende:2003, page 37]

    # Δv: velocity step (in m/s)
    # a_prev: acceleration from previous data point
    # v_prev: velocitiy from previous data point
    Δs = ((v_prev + Δv)^2 - v_prev^2)/2/a_prev      # step size (in m)
    return Δs
end #function calc_Δs_with_Δv

function calc_Δt_with_Δs(Δs::Real, a_prev::Real, v_prev::Real)
    # equation is based on [Wende:2003, page 37]

    # Δs: distance step (in m)
    # a_prev: acceleration from previous data point
    # v_prev: velocitiy from previous data point

    Δt = sign(a_prev) *sqrt((v_prev /a_prev)^2 + 2 *Δs /a_prev) - v_prev /a_prev       # step size (in m/s)
    return Δt
end #function calc_Δt_with_Δs

function calc_Δt_with_Δv(Δv::Real, a_prev::Real)
    # equation is based on [Wende:2003, page 37]

    # Δv: velocity step (in m/s)
    # a_prev: acceleration from previous data point
    Δt = Δv /a_prev        # step size (in s)
    return Δt
end #function calc_Δt_with_Δv

function calc_Δt_with_constant_v(Δs::Real, v::Real)
    # equation is based on [Wende:2003, page 37]

    # Δs: distance step (in m)
    # v: constant velocity (in m/s)
    Δt = Δs /v        # step size (in s)
    return Δt
end #function calc_Δt_with_constant_v

function calc_Δv_with_Δs(Δs::Real, a_prev::Real, v_prev::Real)
    # equation is based on [Wende:2003, page 37]

    # Δs: distance step (in m)
    # a_prev: acceleration from previous data point
    # v_prev: velocitiy from previous data point
    Δv = sqrt(v_prev^2 + 2*Δs*a_prev) - v_prev      # step size (in m/s)
    return Δv
end #function calc_Δv_with_Δs

function calc_Δv_with_Δt(Δt::Real, a_prev::Real)
    # equation is based on [Wende:2003, page 37]

    # Δt: time step (in s)
    # a_prev: acceleration from previous data point
    Δv = Δt * a_prev        # step size (in m/s)
    return Δv
end #function calc_Δv_with_Δt

function calc_ΔW(F_T_prev::Real, Δs::Real)
    # equation is based on [Wende:2003, page 17]

    # F_T_prev: tractive force from previous data point
    # Δs: distance step
    ΔW = F_T_prev * Δs      # mechanical work in this step (in Ws)
    return ΔW
end #function calc_ΔW

function calc_ΔE(ΔW::Real)
    # simplified equation
    # TODO!
    # ΔW: mechanical work in this step (in Ws)
    ΔE = ΔW                 # energy consumption in this step (in Ws)
    return ΔE
end #function calc_ΔW

function calcBrakingDistance(v_start::Real, v_end::Real, a_braking::Real)
    # equation is based on [Wende:2003, page 37]

    # v_start: velocity at the start of braking (in m/s)
    # v_end: target velocity at the end of braking (in m/s)
    # a_braking: constant braking acceleration (in m/s^2)
    s_braking = (v_end^2 - v_start^2) /2 /a_braking             # braking distance (in m)
    # TODO: also possible: calc_Δs_with_Δv(v_end-v_start, a_braking, v_start)
#    return max(0.0, ceil(s_braking, digits=approxLevel))         # ceil is used to be sure that the train stops at s_exit in spite of rounding errors
    return max(0.0, ceil(s_braking, digits=approxLevel +1))         # ceil is used to be sure that the train stops at s_exit in spite of rounding errors
end #function calcBrakingDistance

function calcBrakingStartVelocity(v_end::Real, a_braking::Real, s_braking::Real)
    # equation is based on [Wende:2003, page 37]

    # v_end: target velocity at the end of braking (in m/s)
    # a_braking: constant braking acceleration (in m/s^2)
    # s_braking: braking distance (in Ws)
    v_start = sqrt(v_end^2 - 2*a_braking *s_braking)          # braking start velocity (in m/s)
#    return floor(v_start, digits=approxLevel)
    return floor(v_start, digits=approxLevel +1)
end #function calcBrakingStartVelocity

function calcBrakingAcceleration(v_start::Real, v_end::Real, s_braking::Real)
    # equation is based on [Wende:2003, page 37]

    # v_start: braking start velocity (in m/s)
    # v_end: target velocity at the end of braking (in m/s)
    # s_braking: braking distance (in Ws)
    a_braking = (v_end^2 - v_start^2) /2 /s_braking       # constant braking acceleration (in m/s^2)
    return a_braking
end #function calcBrakingAcceleration


text:
module GmshDiscreteModelsTests

using GridapGmsh
using Gridap.Geometry

mshfile = joinpath(@__DIR__,"t1.msh")
model = GmshDiscreteModel(mshfile)
test_discrete_model(model)

mshfile = joinpath(@__DIR__,"twoTetraeder.msh")
model = GmshDiscreteModel(mshfile; renumber=true)
test_discrete_model(model)

mshfile = joinpath(@__DIR__,"..","demo","demo.msh")
model = GmshDiscreteModel(mshfile)
test_discrete_model(model)

mshfile = joinpath(@__DIR__,"square.msh")
model = GmshDiscreteModel(mshfile)
test_discrete_model(model)

mshfile = joinpath(@__DIR__,"cube.msh")
model = GmshDiscreteModel(mshfile)
test_discrete_model(model)

end # module


text:
"""
Run two coarse rising thermal bubble simulations and make sure that when
restarting from a checkpoint, the restarted simulation matches the non-restarted
simulation numerically.
"""
function run_thermal_bubble_checkpointer_tests()
    Nx, Ny, Nz = 16, 16, 16
    Lx, Ly, Lz = 100, 100, 100
    Δt = 6

    true_model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), ν=4e-2, κ=4e-2)

    # Add a cube-shaped warm temperature anomaly that takes up the middle 50%
    # of the domain volume.
    i1, i2 = round(Int, Nx/4), round(Int, 3Nx/4)
    j1, j2 = round(Int, Ny/4), round(Int, 3Ny/4)
    k1, k2 = round(Int, Nz/4), round(Int, 3Nz/4)
    true_model.tracers.T.data[i1:i2, j1:j2, k1:k2] .+= 0.01

    checkpointed_model = deepcopy(true_model)

    time_step!(true_model, 10, Δt)

    checkpointer = Checkpointer(dir=".", prefix="test_", frequency=5, padding=1)
    push!(checkpointed_model.output_writers, checkpointer)

    # Checkpoint should be saved as "test_model_checkpoint_5.jld" after the
    # 5th iteration.
    time_step!(checkpointed_model, 5, Δt)

    # Remove all knowledge of the checkpointed model.
    checkpointed_model = nothing

    restored_model = restore_from_checkpoint("test_model_checkpoint_5.jld")

    time_step!(restored_model, 5, Δt)

    # Now the true_model and restored_model should be identical.
    @test all(restored_model.velocities.u.data .≈ true_model.velocities.u.data)
    @test all(restored_model.velocities.v.data .≈ true_model.velocities.v.data)
    @test all(restored_model.velocities.w.data .≈ true_model.velocities.w.data)
    @test all(restored_model.tracers.T.data    .≈ true_model.tracers.T.data)
    @test all(restored_model.tracers.S.data    .≈ true_model.tracers.S.data)
    @test all(restored_model.G.Gu.data         .≈ true_model.G.Gu.data)
    @test all(restored_model.G.Gv.data         .≈ true_model.G.Gv.data)
    @test all(restored_model.G.Gw.data         .≈ true_model.G.Gw.data)
    @test all(restored_model.G.GT.data         .≈ true_model.G.GT.data)
    @test all(restored_model.G.GS.data         .≈ true_model.G.GS.data)
end

"""
Run a coarse thermal bubble simulation and save the output to NetCDF at the
10th time step. Then read back the output and test that it matches the model's
state.
"""
function run_thermal_bubble_netcdf_tests()
    Nx, Ny, Nz = 16, 16, 16
    Lx, Ly, Lz = 100, 100, 100
    Δt = 6

    model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), ν=4e-2, κ=4e-2)

    # Add a cube-shaped warm temperature anomaly that takes up the middle 50%
    # of the domain volume.
    i1, i2 = round(Int, Nx/4), round(Int, 3Nx/4)
    j1, j2 = round(Int, Ny/4), round(Int, 3Ny/4)
    k1, k2 = round(Int, Nz/4), round(Int, 3Nz/4)
    model.tracers.T.data[i1:i2, j1:j2, k1:k2] .+= 0.01

    nc_writer = NetCDFOutputWriter(dir=".", prefix="test_", frequency=10, padding=1)
    push!(model.output_writers, nc_writer)

    time_step!(model, 10, Δt)

    u = read_output(nc_writer, "u", 10)
    v = read_output(nc_writer, "v", 10)
    w = read_output(nc_writer, "w", 10)
    T = read_output(nc_writer, "T", 10)
    S = read_output(nc_writer, "S", 10)

    @test all(u .≈ model.velocities.u.data)
    @test all(v .≈ model.velocities.v.data)
    @test all(w .≈ model.velocities.w.data)
    @test all(T .≈ model.tracers.T.data)
    @test all(S .≈ model.tracers.S.data)
end


text:
using BranchAndBound
using JuMP
using Ipopt

const BB = BranchAndBound

# continuous relaxation of MINLP
m = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))

# The following will be supported in the later version of JuMP.
# @variable(m, x >= 0, Int)
# @variable(m, y >= 0, Int)
# undo_relax = JuMP.relax_integrality(m)

@variable(m, x >= 0) # assume integrality
@variable(m, y >= 0) # assume integrality

@objective(m, Min, -x - y)
@constraint(m, -2*x*y + 2*y >= 1)
@constraint(m, -8*x + 10*y <= 13)

# Initialize BNB tree
tree = BB.initialize_tree(m)

# User-defined branch function
function BB.branch!(tree::BB.AbstractTree, node::BB.AbstractNode)
    @info " Node id $(node.id), status $(node.solution_status), bound $(node.bound)"
    if node.bound >= tree.best_incumbent
        @info " Fathomed by bound"
    elseif node.solution_status in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
        if !isinteger(node.solution[x])
            @info " Branch at x = $(node.solution[x]), x >= $(ceil(node.solution[x])), x <= $(floor(node.solution[x]))"
            child1 = BB.create_child_node_with_lb(node, x, ceil(node.solution[x]))
            child2 = BB.create_child_node_with_ub(node, x, floor(node.solution[x]))
            BB.push!(tree, child1)
            BB.push!(tree, child2)
        elseif !isinteger(node.solution[y])
            @info " Branch at y = $(node.solution[y]), y >= $(ceil(node.solution[y])), y <= $(floor(node.solution[y]))"
            child1 = BB.create_child_node_with_lb(node, y, ceil(node.solution[y]))
            child2 = BB.create_child_node_with_ub(node, y, floor(node.solution[y]))
            BB.push!(tree, child1)
            BB.push!(tree, child2)
        elseif node.bound < tree.best_incumbent
            @info " New incumbent bound: $(node.bound)"
            tree.best_incumbent = node.bound
        end
    else
        @info " Fathomed by solution status: $(node.solution_status)"
    end
end

function BB.termination(tree::BB.AbstractTree)
    @info "Tree nodes: processed $(length(tree.processed)), left $(length(tree.nodes)), total $(tree.node_counter), best bound $(tree.best_bound), best incumbent $(tree.best_incumbent)"
    if BB.isempty(tree)
        @info "Completed the tree search"
        return true
    end
    if length(tree.processed) >= 20
        @info "Reached node limit"
        return true
    end
    return false
end

BB.run(tree)


text:
include("google53_gates.jl")

export Lattice53, rand_google53

struct Lattice53
    labels::Matrix{Int}
end

function Lattice53(;nbits::Int=53)
    config = ones(Bool, 5, 12)
    config[end,2:2:end] .= false
    config[1, 7] = false
    labels = zeros(Int, 5, 12)
    k = 0
    for (i,c) in enumerate(config)
        if c
            k += 1
            labels[i] = k
            k>=nbits && break
        end
    end
    return Lattice53(labels)
end

nbits(lattice::Lattice53) = maximum(lattice.labels)

function Base.getindex(lattice::Lattice53, i, j)
    1<=i<=size(lattice.labels, 1) && 1<=j<=size(lattice.labels, 2) ? lattice.labels[i,j] : 0
end
upperleft(lattice::Lattice53,i,j) = lattice[i-j%2,j-1]
lowerleft(lattice::Lattice53,i,j) = lattice[i+(j-1)%2,j-1]
upperright(lattice::Lattice53,i,j) = lattice[i-j%2,j+1]
lowerright(lattice::Lattice53,i,j) = lattice[i+(j-1)%2,j+1]

function pattern53(lattice::Lattice53, chr::Char)
    res = Tuple{Int,Int}[]
    # i0, di, j0, dj and direction
    di = 1 + (chr>'D')
    dj = 2 - (chr>'D')
    j0 = 1 + min(dj-1, mod(chr-'A',2))
    direction = 'C'<=chr<='F' ? lowerright : upperright
    for j=j0:dj:12
        i0 = chr>'D' ? mod((chr-'D') + (j-(chr>='G'))÷2, 2) : 1
        for i = i0:di:5
            src = lattice[i, j]
            dest = direction(lattice, i, j)
            src!=0 && dest !=0 && push!(res, (src, dest))
        end
    end
    return res
end

function print_lattice53(lattice, pattern)
    for i_=1:10
        i = (i_+1)÷2
        for j=1:12
            if i_%2 == j%2 && lattice[i,j]!=0
                print(" ∘  ")
            else
                print("    ")
            end
        end
        println()
        for j=1:12
            if i_%2 == j%2 && lattice[i,j]!=0
                hasll = (lowerleft(lattice, i, j), lattice[i,j]) in pattern
                haslr = (lattice[i,j], lowerright(lattice, i, j)) in pattern
                print(hasll ? "/ " : "  ")
                print(haslr ? " \\" : "  ")
            else
                print("    ")
            end
        end
        println()
    end
end

"""
    rand_google53(depth::Int; seed=1, nbits=53) -> AbstactBlock

random google supremacy circuit with 53 qubits.
"""
function rand_google53(depth::Int; seed=1, nbits::Int=53)
    Random.seed!(seed)
    c = chain(nbits)
    lattice = Lattice53(nbits=nbits)
    k = 0
    for pattern in Iterators.cycle(['A', 'B', 'C', 'D', 'C', 'D', 'A', 'B'])
        push!(c, rand_google53_layer(lattice, pattern))
        k += 1
        k>=depth && break
    end
    return c
end

function rand_google53_layer(lattice, pattern)
    nbit = nbits(lattice)
    chain(chain([put(nbit, i=>rand([SqrtW, SqrtX, SqrtY])) for i=1:nbit]),
        chain([entangler_google53(nbit,i,j) for (i,j) in pattern53(lattice, pattern)])
        )
end


text:
using CSV
using DataFrames
using Missings
using Query

filename = Pkg.dir("Multilane", "data/baseline_curve_Thursday_22_Mar_18_41.csv")
data = CSV.read(filename, nullable=true)

data[:time_to_lane] = data[:steps_to_lane].*data[:dt]
data[:brakes_per_step] = data[:nb_brakes]./data[:n_steps]

combined = by(data, [:solver, :lambda]) do df
    n = size(df, 1)
    return DataFrame(reached_lane=sum(df[:terminal].=="lane")/n,
                     had_brakes=sum(df[:nb_brakes].>=1)/n,
                    )

    # return DataFrame(ttl=mean(df[:time_to_lane]),
    #                  ttl_sem=std(df[:time_to_lane])./sqrt(size(df,1)),
    #                  bps=mean(df[:brakes_per_step]),
    #                  bps_sem=std(df[:brakes_per_step])./sqrt(size(df,1)),
    #                  nb_missing=sum(df[:steps_to_lane].==1000)
    #                 )
end

@show combined

# selected = @from r in data begin
#     @where r.lambda == 2
#     @select r.time_to_lane
#     @collect
# end


text:
# Eric S. Tellez <eric.tellez@infotec.mx>

#Pkg.add("GZip")
#Pkg.add("Glob")
#Pkg.add("JSON")

using ArgParse
import JSON
import Glob

include("textmodel.jl")
include("io.jl")

function delete(fun, input, keys)
    itertweets(input) do tweet
        for key in keys
            delete!(tweet, key)
        end

        fun(tweet)
    end
end

function select(fun, input, keys)
    itertweets(input) do tweet
        a = Dict()
        for key in keys
            a[key] = getkeypath(tweet, key)
        end

        fun(a)
    end
end

function main()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--delete"
        help = "deletes all pairs with the keys listed here"
        arg_type = String
        default = ""
        "--select"
        help = "selects all pairs with the keys listed here"
        arg_type = String
        default = ""
    end

    args = parse_args(ARGS, s)

    if length(args["delete"]) > 0
        delete(STDIN, split(args["delete"], ',')) do tweet
            write(STDOUT, JSON.json(tweet), '\n')
        end
    elseif length(args["select"]) > 0
        select(STDIN, split(args["select"], ',')) do tweet
            write(STDOUT, JSON.json(tweet), '\n')
        end
    end
end

if length(ARGS) > 0
    main()
end

max_stars_repo_name: INGEOTEC/PatternAnalysisTools

max_issues_repo_name: INGEOTEC/PatternAnalysisTools

max_forks_repo_name: INGEOTEC/PatternAnalysisTools


text:
const td = mktempdir()
flush_td() = (isdir(td) && rm(td; recursive=true); mkdir(td))
F.FOLDER_PATH[] = td

fd2html_td(e)  = fd2html(e; dir=td)
fd2html_tdv(e) = F.fd2html_v(e; dir=td)

F.def_GLOBAL_VARS!()
F.def_GLOBAL_LXDEFS!()

mkdir(F.PATHS[:libs])
# copying _libs/katex in the F.PATHS[:libs] so that it can be used in testing
# the js_prerender_math
cp(joinpath(dirname(dirname(pathof(Franklin))), "test", "_libs", "katex"), joinpath(F.PATHS[:libs], "katex"))

@testset "Set vars" begin
    d = F.PageVars(
    	"a" => 0.5 => (Real,),
    	"b" => "hello" => (String, Nothing))
    F.set_vars!(d, ["a"=>"5", "b"=>"nothing"])

    @test d["a"].first == 5
    @test d["b"].first === nothing

    @test_throws F.PageVariableError F.set_vars!(d, ["a"=> "sqrt(-1)"])

    # assigning new variables

    F.set_vars!(d, ["blah"=>"1"])
    @test d["blah"].first == 1
end

@testset "Def+coms" begin # see #78
    st = raw"""
        @def title = "blah" <!-- comment -->
        @def hasmath = false
        etc
        """
    m = F.convert_md(st)
    @test F.locvar("title") == "blah"
    @test F.locvar("hasmath") == false
end


text:
function ufpair(rt, data :: Data)

	@unpack sig6, eps4, eps24 = data

	r6 = (rt[3])^3
	r4 = (rt[3])^2

	u = eps4*(sig6/r6 - 1.0)/r6

	f = eps24 * (r6 - 2*sig6) / (r6*r4*r4)

	fx = f * rt[1]
	fy = f * rt[2]

	return u, (-fx, -fy)
end


text:
@testset "invert" begin
    # arrays of arrays
    a = [[1,2,3],[4,5,6]]
    @test @inferred(invert(a))::Vector{Vector{Int}} == [[1,4],[2,5],[3,6]]

    b = [[0,0],[0,0],[0,0]]
    @test invert!(b, a) == [[1,4],[2,5],[3,6]]

    # Issue #14
    a3 = [[1 2; 3 4],[5 6; 7 8]]
    result3 = Matrix{Vector{Int}}(undef,2,2) # literal syntax is `hvcat`... :(
    result3[1, 1] = [1, 5]
    result3[1, 2] = [2, 6]
    result3[2, 1] = [3, 7]
    result3[2, 2] = [4, 8]
    @test @inferred(invert(a3))::Matrix{Vector{Int}} == result3

    # tuples of tuples
    tt = ((true, 2.0, 3f0), ('d', 0x05, 6))
    @test @inferred(invert(tt)) === ((true, 'd'), (2.0, 0x05), (3f0, 6))

    # tuples / arrays
    ta = ([1, 2], [3.1, 4.1])
    at = [(1, 3.1), (2, 4.1)]
    @test @inferred(invert(ta))::Vector{Tuple{Int, Float64}} == at
    @test @inferred(invert(at))::Tuple{Vector{Int}, Vector{Float64}} == ta

    # named tuples
    nn = (a = (x = true, y = 2), b = (x = 3.0, y = 4.0f0))
    @test @inferred(invert(nn)) === (x = (a = true, b = 3.0), y = (a = 2, b = 4.0f0))

    # tuples / named tuples
    tn = ((a = true, b = 2), (a = 3.0, b = 4.0f0))
    nt = (a = (true, 3.0), b = (2, 4.0f0))
    @test @inferred(invert(nt)) === tn
    @test @inferred(invert(tn)) === nt

    # arrays / named tuples
    na = (a = [1,2,3], b = [2.0, 4.0, 6.0])
    an = [(a = 1, b = 2.0), (a = 2, b = 4.0), (a = 3, b = 6.0)]
    @test @inferred(invert(na))::typeof(an) == an
    @test @inferred(invert(an))::typeof(na) == na
end


text:
# H = Σⱼ XⱼXⱼ₊₁ + YⱼYⱼ₊₁
function ITensors.MPO(::Model"xx", s)
  N = length(s)
  os = OpSum()
  for n in 1:(N - 1)
    os .+= 1, "X", n, "X", n + 1
    os .+= 1, "Y", n, "Y", n + 1
  end
  return splitblocks(linkinds, MPO(os, s))
end

# H = X₁X₂ + Y₁Y₂
function ITensors.OpSum(::Model"xx", n1, n2)
  opsum = OpSum()
  for j in 1:(N - 1)
    opsum .+= 0.5, "S+", j, "S-", j + 1
    opsum .+= 0.5, "S-", j, "S+", j + 1
  end
  return opsum
end

# H = X₁X₂ + Y₁Y₂
# XXX: use `op` instead of `ITensor`
function ITensors.ITensor(::Model"xx", s1::Index, s2::Index)
  # op("X", s1) * op("X", s2) + op("Y", s1) * op("Y", s2)
  opsum = OpSum()
  n = 1
  opsum += "X", n, "X", n + 1
  opsum += "Y", n, "Y", n + 1
  return prod(MPO(opsum, [s1, s2]))
end

function reference(::Model"xx", ::Observable"energy"; N=∞)
  isinf(N) && return -4 / π
  # Exact eigenvalues of uniform symmetric tridiagonal matrix
  λ(k) = cos(k * π / (N + 1))
  return 4 * sum(k -> min(λ(k), 0.0), 1:N) / N
end


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

name = "libsigcpp"
version = v"3.0.7"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://download.gnome.org/sources/libsigc++/3.0/libsigc%2B%2B-$(version).tar.xz",
                  "bfbe91c0d094ea6bbc6cbd3909b7d98c6561eea8b6d9c0c25add906a6e83d733")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libsigc++*/
mkdir meson
cd meson
meson --cross-file=${MESON_TARGET_TOOLCHAIN}
ninja -j${nproc}
ninja install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms(; experimental=true))

# The products that we will ensure are always built
products = [
    LibraryProduct(["libsigc-3.0", "libsigc-3"], :libsigc)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version = v"7.1.0")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
@enum HtlParserState STATE_DATA STATE_TAG_OPEN STATE_END_TAG_OPEN STATE_TAG_NAME STATE_BEFORE_ATTRIBUTE_NAME STATE_AFTER_ATTRIBUTE_NAME STATE_ATTRIBUTE_NAME STATE_BEFORE_ATTRIBUTE_VALUE STATE_ATTRIBUTE_VALUE_DOUBLE_QUOTED STATE_ATTRIBUTE_VALUE_SINGLE_QUOTED STATE_ATTRIBUTE_VALUE_UNQUOTED STATE_AFTER_ATTRIBUTE_VALUE_QUOTED STATE_SELF_CLOSING_START_TAG STATE_COMMENT_START STATE_COMMENT_START_DASH STATE_COMMENT STATE_COMMENT_LESS_THAN_SIGN STATE_COMMENT_LESS_THAN_SIGN_BANG STATE_COMMENT_LESS_THAN_SIGN_BANG_DASH STATE_COMMENT_LESS_THAN_SIGN_BANG_DASH_DASH STATE_COMMENT_END_DASH STATE_COMMENT_END STATE_COMMENT_END_BANG STATE_MARKUP_DECLARATION_OPEN STATE_RAWTEXT STATE_RAWTEXT_LESS_THAN_SIGN STATE_RAWTEXT_END_TAG_OPEN STATE_RAWTEXT_END_TAG_NAME

is_alpha(ch) = 'A' <= ch <= 'Z' || 'a' <= ch <= 'z'
is_space(ch) = ch in ('\t', '\n', '\f', ' ')
normalize(s) = replace(replace(s, "\r\n" => "\n"), "\r" => "\n")
nearby(x,i) = i+10>length(x) ? x[i:end] : x[i:i+8] * "…"

"""
    interpolate(args)::Expr

Take an interweaved set of Julia expressions and strings, tokenize the
strings according to the HTML specification [1], wrapping the
expressions with wrappers based upon the escaping context, and returning
an expression that combines the result with an `Result` wrapper.

For these purposes, a `Symbol` is treated as an expression to be
resolved; while a `String` is treated as a literal string that won't be
escaped. Critically, interpolated strings to be escaped are represented
as an `Expr` with `head` of `:string`.

There are tags, "script" and "style" which are rawtext, in these cases
there is no escaping, and instead raise an exception if the appropriate
ending tag is in substituted content.

[1] https://html.spec.whatwg.org/multipage/parsing.html#tokenization
"""
function interpolate(args)
    state = STATE_DATA
    parts = Union{String,Expr}[]
    attribute_start = attribute_end = 0
    element_start = element_end = 0
    buffer_start = buffer_end = 0
    attribute_tag = nothing
    element_tag = nothing
    state_tag_is_open = false

    function choose_tokenizer()
        if state_tag_is_open
            if element_tag in (:style, :xmp, :iframe, :noembed,
                               :noframes, :noscript, :script)
                return STATE_RAWTEXT
            end
        end
        return STATE_DATA
    end

    args = [a for a in args if a != ""]

    for j in 1:length(args)
        input = args[j]
        if !isa(input, String)
            if state == STATE_DATA || state == STATE_COMMENT
                push!(parts, :(content($(esc(input)))))
            elseif state == STATE_RAWTEXT
                if :script === element_tag
                    push!(parts, :(ScriptTag($(esc(input)))))
                elseif :style === element_tag
                    push!(parts, :(StyleTag($(esc(input)))))
                else
                    throw(DomainError(element_tag,
                      "Only script and style rawtext tags are supported."))
                end
            elseif state == STATE_BEFORE_ATTRIBUTE_VALUE
                state = STATE_ATTRIBUTE_VALUE_UNQUOTED
                # rewrite previous string to remove ` attname=`
                @assert parts[end] isa String
                name = parts[end][attribute_start:attribute_end]
                parts[end] = parts[end][1:(attribute_start-2)]
                attribute = normalize_attribute_name(name)
                push!(parts, :(attribute_pair($attribute, $(esc(input)))))
                # peek ahead to ensure we have a delimiter
                if j < length(args)
                    next = args[j+1]
                    if next isa String && !occursin(r"^[\s+\/>]", next)
                        msg = "$(name)=$(nearby(next,1))"
                        throw(DomainError(msg, "Unquoted attribute " *
                          "interpolation is limited to a single component"))
                    end
                end
            elseif state == STATE_ATTRIBUTE_VALUE_UNQUOTED
                throw(DomainError(input, "Unquoted attribute " *
                  "interpolation is limited to a single component"))
            elseif state == STATE_ATTRIBUTE_VALUE_SINGLE_QUOTED ||
                   state == STATE_ATTRIBUTE_VALUE_DOUBLE_QUOTED
                @assert parts[end] isa String
                name = parts[end][attribute_start:attribute_end]
                push!(parts, :(attribute_value($(esc(input)))))
            elseif state == STATE_BEFORE_ATTRIBUTE_NAME ||
                   state == STATE_AFTER_ATTRIBUTE_NAME
                # strip space before interpolated element pairs
                if parts[end] isa String
                    if parts[end][end] == ' '
                       parts[end] = parts[end][1:length(parts[end])-1]
                    end
                end
                # ensure a space between this and next attribute
                if j < length(args)
                    next = args[j+1]
                    if next isa String && !occursin(r"^[\s+\/>]", next)
                        args[j+1] = " " * next
                    end
                end
                append!(parts, rewrite_inside_tag(input))
            else
                throw("unexpected binding $(state)")
            end
        else
            input = normalize(input)
            inputlength = lastindex(input)
            i = 1
            while i <= inputlength
                ch = input[i]

                if state == STATE_DATA
                    if ch === '<'
                        state = STATE_TAG_OPEN
                    end

                elseif state == STATE_RAWTEXT
                    if ch === '<'
                        state = STATE_RAWTEXT_LESS_THAN_SIGN
                    end

                elseif state == STATE_TAG_OPEN
                    if ch === '!'
                        state = STATE_MARKUP_DECLARATION_OPEN
                    elseif ch === '/'
                        state = STATE_END_TAG_OPEN
                    elseif is_alpha(ch)
                        state = STATE_TAG_NAME
                        state_tag_is_open = true
                        element_start = i
                        i = prevind(input, i)
                    elseif ch === '?'
                        # this is an XML processing instruction, with
                        # recovery production called "bogus comment"
                        throw(DomainError(nearby(input, i-1),
                          "unexpected question mark instead of tag name"))
                    else
                        throw(DomainError(nearby(input, i-1),
                          "invalid first character of tag name"))
                    end

                elseif state == STATE_END_TAG_OPEN
                    @assert !state_tag_is_open
                    if is_alpha(ch)
                        state = STATE_TAG_NAME
                        i = prevind(input, i)
                    elseif ch === '>'
                        state = STATE_DATA
                    else
                        throw(DomainError(nearby(input, i-1),
                          "invalid first character of tag name"))
                    end

                elseif state == STATE_TAG_NAME
                    if isspace(ch) || ch === '/' || ch === '>'
                        if state_tag_is_open
                            element_tag = Symbol(lowercase(
                                            input[element_start:element_end]))
                            element_start = element_end = 0
                        end
                        if isspace(ch)
                            state = STATE_BEFORE_ATTRIBUTE_NAME
                            # subordinate states use state_tag_is_open flag
                        elseif ch === '/'
                            state = STATE_SELF_CLOSING_START_TAG
                            state_tag_is_open = false
                        elseif ch === '>'
                            state = choose_tokenizer()
                            state_tag_is_open = false
                        end
                    else
                        if state_tag_is_open
                            element_end = i
                        end
                    end

                elseif state == STATE_BEFORE_ATTRIBUTE_NAME
                    if is_space(ch)
                        nothing
                    elseif ch === '/' || ch === '>'
                        state = STATE_AFTER_ATTRIBUTE_NAME
                        i = prevind(input, i)
                    elseif ch in  '='
                        throw(DomainError(nearby(input, i-1),
                          "unexpected equals sign before attribute name"))
                    else
                        state = STATE_ATTRIBUTE_NAME
                        attribute_start = i
                        attribute_end = nothing
                        i = prevind(input, i)
                    end

                elseif state == STATE_ATTRIBUTE_NAME
                    if is_space(ch) || ch === '/' || ch === '>'
                        state = STATE_AFTER_ATTRIBUTE_NAME
                        i = prevind(input, i)
                    elseif ch === '='
                        state = STATE_BEFORE_ATTRIBUTE_VALUE
                    elseif ch in ('"', '\"', '<')
                        throw(DomainError(nearby(input, i-1),
                          "unexpected character in attribute name"))
                    else
                        attribute_end = i
                    end

                elseif state == STATE_AFTER_ATTRIBUTE_NAME
                    if is_space(ch)
                        nothing
                    elseif ch === '/'
                        state = STATE_SELF_CLOSING_START_TAG
                    elseif ch === '='
                        state = STATE_BEFORE_ATTRIBUTE_VALUE
                    elseif ch === '>'
                        state = choose_tokenizer()
                        state_tag_is_open = false
                    else
                        state = STATE_ATTRIBUTE_NAME
                        attribute_start = i
                        attribute_end = nothing
                        i = prevind(input, i)
                    end

                elseif state == STATE_BEFORE_ATTRIBUTE_VALUE
                    if is_space(ch)
                        nothing
                    elseif ch === '"'
                        attribute_tag = input[attribute_start:attribute_end]
                        state = STATE_ATTRIBUTE_VALUE_DOUBLE_QUOTED
                    elseif ch === '\''
                        attribute_tag = input[attribute_start:attribute_end]
                        state = STATE_ATTRIBUTE_VALUE_SINGLE_QUOTED
                    elseif ch === '>'
                        throw(DomainError(nearby(input, i-1),
                          "missing attribute value"))
                    else
                        state = STATE_ATTRIBUTE_VALUE_UNQUOTED
                        i = prevind(input, i)
                    end

                elseif state == STATE_ATTRIBUTE_VALUE_DOUBLE_QUOTED
                    if ch === '"'
                        state = STATE_AFTER_ATTRIBUTE_VALUE_QUOTED
                        attribute_tag = nothing
                    end

                elseif state == STATE_ATTRIBUTE_VALUE_SINGLE_QUOTED
                    if ch === '\''
                        state = STATE_AFTER_ATTRIBUTE_VALUE_QUOTED
                        attribute_tag = nothing
                    end

                elseif state == STATE_ATTRIBUTE_VALUE_UNQUOTED
                    if is_space(ch)
                        state = STATE_BEFORE_ATTRIBUTE_NAME
                    elseif ch === '>'
                        state = choose_tokenizer()
                        state_tag_is_open = false
                    elseif ch in ('"', '\'', "<", "=", '`')
                        throw(DomainError(nearby(input, i-1),
                          "unexpected character in unquoted attribute value"))
                    end

                elseif state == STATE_AFTER_ATTRIBUTE_VALUE_QUOTED
                    if is_space(ch)
                        state = STATE_BEFORE_ATTRIBUTE_NAME
                    elseif ch === '/'
                        state = STATE_SELF_CLOSING_START_TAG
                    elseif ch === '>'
                        state = choose_tokenizer()
                        state_tag_is_open = false
                    else
                        throw(DomainError(nearby(input, i-1),
                          "missing whitespace between attributes"))
                    end

                elseif state == STATE_SELF_CLOSING_START_TAG
                    if ch === '>'
                        state = STATE_DATA
                    else
                        throw(DomainError(nearby(input, i-1),
                          "unexpected solidus in tag"))
                    end

                elseif state == STATE_MARKUP_DECLARATION_OPEN
                    if ch === '-' && input[i + 1] == '-'
                        state = STATE_COMMENT_START
                        i = nextind(input, i)
                    elseif startswith(input[i:end], "DOCTYPE")
                        throw("DOCTYPE not supported")
                    elseif startswith(input[i:end], "[CDATA[")
                        throw("CDATA not supported")
                    else
                        throw(DomainError(nearby(input, i-1),
                          "incorrectly opened comment"))
                    end

                elseif state == STATE_COMMENT_START
                    if ch === '-'
                        state = STATE_COMMENT_START_DASH
                    elseif ch === '>'
                        throw(DomainError(nearby(input, i-1),
                          "abrupt closing of empty comment"))
                    else
                        state = STATE_COMMENT
                        i = prevind(input, i)
                    end

                elseif state == STATE_COMMENT_START_DASH
                    if ch === '-'
                        state = STATE_COMMENT_END
                    elseif ch === '>'
                        throw(DomainError(nearby(input, i-1),
                          "abrupt closing of empty comment"))
                    else
                        state = STATE_COMMENT
                        i = prevind(input, i)
                    end

                elseif state == STATE_COMMENT
                    if ch === '<'
                        state = STATE_COMMENT_LESS_THAN_SIGN
                    elseif ch === '-'
                        state = STATE_COMMENT_END_DASH
                    end

                elseif state == STATE_COMMENT_LESS_THAN_SIGN
                    if ch === '!'
                        state = STATE_COMMENT_LESS_THAN_SIGN_BANG
                    elseif ch === '<'
                        nothing
                    else
                        state = STATE_COMMENT
                        i = prevind(input, i)
                    end

                elseif state == STATE_COMMENT_LESS_THAN_SIGN_BANG
                    if ch == '-'
                        state = STATE_COMMENT_LESS_THAN_SIGN_BANG_DASH
                    else
                        state = STATE_COMMENT
                        i = prevind(input, i)
                    end

                elseif state == STATE_COMMENT_LESS_THAN_SIGN_BANG_DASH
                    if ch == '-'
                        state = STATE_COMMENT_LESS_THAN_SIGN_BANG_DASH_DASH
                    else
                        state = STATE_COMMENT_END
                        i = prevind(input, i)
                    end

                elseif state == STATE_COMMENT_LESS_THAN_SIGN_BANG_DASH_DASH
                    if ch == '>'
                        state = STATE_COMMENT_END
                        i = prevind(input, i)
                    else
                        throw(DomainError(nearby(input, i-1),
                          "nested comment"))
                    end

                elseif state == STATE_COMMENT_END_DASH
                    if ch === '-'
                        state = STATE_COMMENT_END
                    else
                        state = STATE_COMMENT
                        i = prevind(input, i)
                    end

                elseif state == STATE_COMMENT_END
                    if ch === '>'
                        state = STATE_DATA
                    elseif ch === '!'
                        state = STATE_COMMENT_END_BANG
                    elseif ch === '-'
                        nothing
                    else
                        state = STATE_COMMENT
                        i = prevind(input, i)
                    end

                elseif state == STATE_COMMENT_END_BANG
                    if ch === '-'
                        state = STATE_COMMENT_END_DASH
                    elseif ch === '>'
                        throw(DomainError(nearby(input, i-1),
                          "nested comment"))
                    else
                        state = STATE_COMMENT
                        i = prevind(input, i)
                    end

                elseif state == STATE_RAWTEXT_LESS_THAN_SIGN
                    if ch === '/'
                        state = STATE_RAWTEXT_END_TAG_OPEN
                    elseif ch === '!' && element_tag == :script
                        # RAWTEXT differs from SCRIPT here
                        throw("script escape or comment is not implemented")
                    else
                        state = STATE_RAWTEXT
                        # do not "reconsume", even though spec says so
                    end

                elseif state == STATE_RAWTEXT_END_TAG_OPEN
                    if is_alpha(ch)
                        state = STATE_RAWTEXT_END_TAG_NAME
                        buffer_start = i
                        i = prevind(input, i)
                    else
                        state = STATE_RAWTEXT
                        i = prevind(input, i)
                    end

                elseif state == STATE_RAWTEXT_END_TAG_NAME
                    if is_alpha(ch)
                        buffer_end = i
                    elseif ch in ('/', '>') || is_space(ch)
                        # test for "appropriate end tag token"
                        current = input[buffer_start:buffer_end]
                        if Symbol(lowercase(current)) == element_tag
                            if ch === '/'
                                state = STATE_SELF_CLOSING_START_TAG
                            elseif ch === '>'
                                state = STATE_DATA
                            else
                                state = STATE_BEFORE_ATTRIBUTE_NAME
                            end
                            continue
                        else
                            state = STATE_RAWTEXT
                        end
                    else
                        state = STATE_RAWTEXT
                    end

                else
                    throw("unhandled state transition")
                end

                i = nextind(input, i)
            end
            push!(parts, input)
        end
    end
    parts = Expr[(x isa String ? :(Bypass($x)) : x) for x in parts]
    return Expr(:call, :Result, parts...)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 7c80e210-3d63-11eb-18de-ef8be75be533

# ╔═╡ ce3d3820-3d62-11eb-1a66-55f477b0a1e8
import Pkg; Pkg.add("Plots")

# ╔═╡ 51cd5580-3d63-11eb-2c0a-b55067ce10f4
1 + 1

# ╔═╡ Cell order:
# ╠═7c80e210-3d63-11eb-18de-ef8be75be533
# ╟─ce3d3820-3d62-11eb-1a66-55f477b0a1e8
# ╠═51cd5580-3d63-11eb-2c0a-b55067ce10f4

max_stars_repo_name: groundf/course-julia

max_issues_repo_name: groundf/course-julia

max_forks_repo_name: groundf/course-julia


text:
using PyPlot

pospart(x) = ( x > 0 ? x : zero(x) )

function K2(x, y, a, b)
    pt = [a, (a+b)/2, b]
    ℓ2 = (x-pt[1])*(x-pt[3])/((pt[2]-pt[1])*(pt[2]-pt[3]))
    ℓ3 = (x-pt[1])*(x-pt[2])/((pt[3]-pt[1])*(pt[3]-pt[2]))
    Qπ = pospart(pt[2]-y)^2 * ℓ2 / 2 + (b-y)^2 * ℓ3 / 2
    return Qπ - pospart(x-y)^2/2
end

a = -1.0
b = 1.0
N = 200
x = LinRange(a, b, N+1)
y = LinRange(a, b, N+1)
z = Float64[ K2(x[k], y[j], a, b) for j=1:N+1, k=1:N+1 ]
fig = figure(1)
surf(x, y, z, cstride=4, rstride=4, cmap="jet", linewidth=0.25)
xlabel(L"x")
ylabel(L"y")
ax = gca()
ax[:view_init](elev=40,azim=-100)
savefig("Quadratic_PeanoK_3d.pdf")

figure(2)
contour(x, y, z, 12)
colorbar()
xlabel(L"x")
ylabel(L"y")
grid(true)
savefig("Quadratic_PeanoK_contour.pdf")

max_stars_repo_name: billmclean/ComputationalMathsNotes

max_issues_repo_name: billmclean/ComputationalMathsNotes

max_forks_repo_name: billmclean/ComputationalMathsNotes


text:
struct Edge{T} <: AbstractLinePrimitive{T}
    a::CartesianPoint{T}
    b::CartesianPoint{T}
end

direction(e::Edge) = e.b - e.a

Line(e::Edge{T}) where {T} = Line{T}(e.a, direction(e))

function distance(pt::CartesianPoint{T}, e::Edge{T}) where {T}
    return if (pt - e.b) ⋅ direction(e) >= 0
        norm(pt - e.b)
    elseif (pt - e.a) ⋅ direction(e) <= 0
        norm(pt - e.a)
    else
        distance(pt, Line(e))
    end
end

function sample(e::Edge{T}; n = 2)::Vector{CartesianPoint{T}} where {T}
    xs = range(zero(T), stop = one(T), length = n)
    pts = Vector{CartesianPoint{T}}(undef, n)
    dir = direction(e)
    for i in eachindex(pts)
        pts[i] = e.a + xs[i] .* dir
    end
    pts
end
function sample(es::Vector{Edge{T}}; n = 2)::Vector{CartesianPoint{T}} where {T}
    pts = Vector{CartesianPoint{T}}(undef, length(es) * n)
    for i in eachindex(es)
        pts[(i-1)*n+1:i*n] = sample(es[i], n = n)
    end
    pts
end


text:
# define macro for statically AD-ing through activation function
#macro frule()

# Define activation functions
@inline ℒ2(x::T) where T = T(2)*cosh(x)
@inline ∂logℒ2(x) = tanh(x)
@inline logℒ2(x::T) where T<:Real =
    (abs(x)<= T(12.0) ? log(cosh(x)) : abs(x) - log(2one(x)))
@inline logℒ2(x::Complex{T}) where T<:Real =
    logℒ2(real(x)) + log(cos(imag(x)) + im * tanh(real(x)) * sin(imag(x)))

@inline ℒ(x) = one(x) + exp(x)
@inline ∂logℒ(x) = one(x)/(one(x)+exp(-x))
#const ∂logℒ = NNlib.σ
#NNlib.σ(x::Complex) = one(x)/(one(x)+exp(-x))
@inline logℒ(x) = log1p(exp(x))#log(one(x) + exp(x))
#const logℒ = NNlib.softplus
#NNlib.softplus(x::Complex) = log1p(exp(x))#log(one(x) + exp(x))

@inline fwd_der(f::typeof(logℒ)) = ∂logℒ
@inline fwd_der(f::typeof(logℒ), x) = ∂logℒ(x)

@inline fwd_der(f::typeof(logℒ2)) = ∂logℒ
@inline fwd_der(f::typeof(logℒ2), x) = ∂logℒ2(x)

@inline fwd_der(f::typeof(identity)) = identity
@inline fwd_der(f::typeof(identity), x::T) where T = one(T)

export af_softplus, af_logcosh
"""
    af_softplus(x) = log( 1 + exp(x) )

Classic softplus (https://en.wikipedia.org/wiki/Sigmoid_function) activation function.
"""
const af_softplus = logℒ

# TODO use NNlib implementation
"""
  logcosh(x)

  Return log(cosh(x)) which is computed in a (somewhat) numerically stable way.
"""
const af_logcosh = logℒ2


text:
abstract Event

immutable StreamStartEvent <: Event
    start_mark::Mark
    end_mark::Mark
    encoding::AbstractString
end

immutable StreamEndEvent <: Event
    start_mark::Mark
    end_mark::Mark
end

immutable DocumentStartEvent <: Event
    start_mark::Mark
    end_mark::Mark
    explicit::Bool
    version::(@compat Union{AbstractString, Void})
    tags::(@compat Union{Dict{AbstractString, AbstractString}, Void})

    function DocumentStartEvent(start_mark::Mark,end_mark::Mark,
                                explicit::Bool, version=nothing,
                                tags=nothing)
        new(start_mark, end_mark, explicit, version, tags)
    end
end

immutable DocumentEndEvent <: Event
    start_mark::Mark
    end_mark::Mark
    explicit::Bool
end

immutable AliasEvent <: Event
    start_mark::Mark
    end_mark::Mark
    anchor::(@compat Union{AbstractString, Void})
end

immutable ScalarEvent <: Event
    start_mark::Mark
    end_mark::Mark
    anchor::(@compat Union{AbstractString, Void})
    tag::(@compat Union{AbstractString, Void})
    implicit::Tuple
    value::AbstractString
    style::(@compat Union{Char, Void})
end

immutable SequenceStartEvent <: Event
    start_mark::Mark
    end_mark::Mark
    anchor::(@compat Union{AbstractString, Void})
    tag::(@compat Union{AbstractString, Void})
    implicit::Bool
    flow_style::Bool
end

immutable SequenceEndEvent <: Event
    start_mark::Mark
    end_mark::Mark
end

immutable MappingStartEvent <: Event
    start_mark::Mark
    end_mark::Mark
    anchor::(@compat Union{AbstractString, Void})
    tag::(@compat Union{AbstractString, Void})
    implicit::Bool
    flow_style::Bool
end

immutable MappingEndEvent <: Event
    start_mark::Mark
    end_mark::Mark
end


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations
gamma = 5/3
equations = CompressibleEulerEquations3D(gamma)

"""
    initial_condition_blob(x, t, equations::CompressibleEulerEquations3D)

The blob test case taken from
- Agertz et al. (2006)
  Fundamental differences between SPH and grid methods
  [arXiv: astro-ph/0610051](https://arxiv.org/abs/astro-ph/0610051)
"""
function initial_condition_blob(x, t, equations::CompressibleEulerEquations3D)
  # blob test case, see Agertz et al. https://arxiv.org/pdf/astro-ph/0610051.pdf
  # other reference: https://arxiv.org/pdf/astro-ph/0610051.pdf
  # change discontinuity to tanh
  # typical domain is rectangular, we change it to a square
  # resolution 128^3, 256^3
  # domain size is [-20.0,20.0]^3
  # gamma = 5/3 for this test case
  R = 1.0 # radius of the blob
  # background density
  rho = 1.0
  Chi = 10.0 # density contrast
  # reference time of characteristic growth of KH instability equal to 1.0
  tau_kh = 1.0
  tau_cr = tau_kh / 1.6 # crushing time
  # determine background velocity
  v1 = 2 * R * sqrt(Chi) / tau_cr
  v2 = 0.0
  v3 = 0.0
  Ma0 = 2.7 # background flow Mach number Ma=v/c
  c = v1 / Ma0 # sound speed
  # use perfect gas assumption to compute background pressure via the sound speed c^2 = gamma * pressure/density
  p = c * c * rho / equations.gamma
  # initial center of the blob
  inicenter = [-15, 0, 0]
  x_rel = x - inicenter
  r = sqrt(x_rel[1]^2 + x_rel[2]^2 + x_rel[3]^2)
  # steepness of the tanh transition zone
  slope = 2
  # density blob
  rho = rho + (Chi - 1) * 0.5 * (1 + (tanh(slope * (r + R)) - (tanh(slope *(r - R)) + 1)))
  # velocity blob is zero
  v1 = v1 - v1 * 0.5 * (1 + (tanh(slope *(r + R)) - (tanh(slope *(r - R)) + 1)))
  return prim2cons(SVector(rho, v1, v2, v3, p), equations)
end
initial_condition = initial_condition_blob

volume_flux = flux_ranocha
solver = DGSEM(polydeg=3, surface_flux=flux_hllc,
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

coordinates_min = (-20.0, -20.0, -20.0)
coordinates_max = ( 20.0,  20.0,  20.0)

refinement_patches = (
  (type="box", coordinates_min=(-20.0, -10.0, -10.0), coordinates_max=(-10.0, 10.0, 10.0)),
  (type="box", coordinates_min=(-20.0,  -5.0,  -5.0), coordinates_max=(-10.0,  5.0,  5.0)),
  (type="box", coordinates_min=(-17.0,  -2.0,  -2.0), coordinates_max=(-13.0,  2.0,  2.0)),
  (type="box", coordinates_min=(-17.0,  -2.0,  -2.0), coordinates_max=(-13.0,  2.0,  2.0)),
)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=2,
                refinement_patches=refinement_patches,
                n_cells_max=100_000,)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 2.5)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 200
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=200,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

amr_indicator = IndicatorLöhner(semi,
                                variable=Trixi.density)
amr_controller = ControllerThreeLevel(semi, amr_indicator,
                                      base_level=1,
                                      med_level =0, med_threshold=0.1, # med_level = current level
                                      max_level =6, max_threshold=0.3)
amr_callback = AMRCallback(semi, amr_controller,
                           interval=3,
                           adapt_initial_condition=false,
                           adapt_initial_condition_only_refine=true)

stepsize_callback = StepsizeCallback(cfl=1.7)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution,
                        amr_callback, stepsize_callback)

stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=(1.0e-4, 1.0e-4),
                                                     variables=(Trixi.density, pressure))

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(stage_limiter!, williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
"""
    metnum.jl

Modul ini berisi metode untuk menyelesaikan

1. Akar pertaksamaan tak-linear: `bisection`, `regulaFalsi`, `fixpoint`, `newtonRaphson`, `secant`
2. SPL langsung: `backsub`, `elimGaussNonPivoting`, `elimGaussWithPivoting`, `LUtanpaP`, `LUdenganP`
3. SPL iteratif: `jacobi`, `gaussSeidel`, `rekons` ,`conGrad`
4. Interpolas: `lagrange`, `newton`
5. Regresi: `reglin`, `regpower`, `regpoly`
6. Turunan: `bedaPusat`, `richardson`
7. Integral: `komptrap`, `kompsimp`, `rekursif`, `romberg`, `adaptif`
8. PDB dengan MNA: `euler`, `heun`, `taylor`, `rungekutta`, `rungekuttasistem`, `rkf45`
9. PDB dengan MNB: `linearshooting`, `findiff`
10. PDP: `gelombang`, `panas`, `dirichlet`
"""
module metnum

using LinearAlgebra
using Statistics
using Dates
using IJulia

export regpoly, regpower, regulaFalsi, rekons, rekursif, richardson, rkf45
export romberg, rungekutta, rungekuttasistem, secant
export taylor, adaptif, backsub, bedaPusat, bisection
export conGrad, dirichlet, elimGaussNonPivoting, elimGaussWithPivoting
export euler, findiff, fixpoint, gaussSeidel, gelombang, heun, jacobi
export kompsimp, komptrap, lagrange, linearshooting, LUdenganP, LUtanpaP
export newton, newtonRaphson, panas, reglin, mulaipraktikum

include("regpoly.jl")
include("regpower.jl")
include("regulaFalsi.jl")
include("rekons.jl")
include("rekursif.jl")
include("richardson.jl")
include("rkf45.jl")
include("romberg.jl")
include("rungekutta.jl")
include("rungekuttasistem.jl")
include("secant.jl")
include("taylor.jl")
include("adaptif.jl")
include("backsub.jl")
include("bedaPusat.jl")
include("bisection.jl")
include("conGrad.jl")
include("dirichlet.jl")
include("elimGaussNonPivoting.jl")
include("elimGaussWithPivoting.jl")
include("euler.jl")
include("findiff.jl")
include("fixpoint.jl")
include("gaussSeidel.jl")
include("gelombang.jl")
include("heun.jl")
include("jacobi.jl")
include("kompsimp.jl")
include("komptrap.jl")
include("lagrange.jl")
include("linearshooting.jl")
include("LUdenganP.jl")
include("LUtanpaP.jl")
include("newton.jl")
include("newtonRaphson.jl")
include("panas.jl")
include("reglin.jl")
include("mulaipraktikum.jl")

end


text:
include("utils.jl")
include("../Applications/SP.jl")

n = 20
k = 2
Γ = 3
inst_id = 1
heuristicmode = "DUAL"

instance = "$(inst_id)_$n.txt"
data = read_data(instance)
K = 1:k
N = 1:data.n

# The following excutions avoid warming-up issues
const STARTTIME = time()
@info "Executing functions to warm up ..."
oldstd = stdout
redirect_stdout(open("/dev/null", "w"))
scenario_generation(false,false)
redirect_stdout(oldstd) # recover original stdout
@info "warm-up finished!"

# Now comes the instance that is tested, passed by the arguments
n = parse(Int, ARGS[1])
k = parse(Int, ARGS[2])
inst_id = parse(Int, ARGS[3])
Γ = parse(Int, ARGS[4])

instance = "$(inst_id)_$n.txt"
@info "Solving instance $inst_id with k=$k and Γ=$Γ and ϵ=$ϵ ..."
data = read_data(instance)
K = 1:k
N = 1:data.n

outputfilename = "/../results/SP_scenario_generation_$(n)_$(k)_$(Γ).csv"
reset_timer!(to::TimerOutput)
STARTTIME = time()
runtime = @elapsed lb, ub, solution, gap, NUMBER_GENSC, NUMITER, NUMBER_GENSOL = scenario_generation(false,true)
outputfile = open(string(dirname(@__FILE__),outputfilename), "a")
output = "$n $inst_id $k $Γ $NUMBER_GENSC $NUMBER_GENSOL $ϵ $runtime $NUMITER $lb $ub $gap $(TimerOutputs.time(to["heuristic"])/10^9) $(TimerOutputs.time(to["gen_sol"])/10^9)"
println(output)
println()
println(outputfile, output)
close(outputfile)
show(to)

# outputfilename = "/../results/SP_rc_generation_$(n)_$(k)_$(Γ).csv"
# reset_timer!(to::TimerOutput)
# STARTTIME = time()
# runtime = @elapsed lb, ub, incumbent, gap, it = rcg()
# outputfile = open(string(dirname(@__FILE__),outputfilename), "a")
# output = "$n $inst_id $k $Γ $runtime $it $lb $ub $gap"
# println(output)
# println()
# println(outputfile, output)
# close(outputfile)
# show(to)

max_stars_repo_name: mjposs/min-max-min

max_issues_repo_name: mjposs/min-max-min

max_forks_repo_name: mjposs/min-max-min


text:
__precompile__(false)

module MagneticLocSuchowiak

using Plots
using DataFrames
using CSV
using MLJ
using Distances
using Statistics

@load XGBoostRegressor
@load NeuralNetworkRegressor
@load KNNRegressor
@load NuSVR
@load RandomForestRegressor pkg=DecisionTree
#@load LGBMRegressor

include("graphs.jl")
include("datareader.jl")
include("models.jl")

export createfullgraph, createcomparisongraphs, readdata, fitknn, testmachines, fitforest, fitnusvr, fitcustommodel

end

max_stars_repo_name: SebastianSuchowiak/MagneticLocSuchowiak

max_issues_repo_name: SebastianSuchowiak/MagneticLocSuchowiak

max_forks_repo_name: SebastianSuchowiak/MagneticLocSuchowiak


text:
"""
Do not distribute the work at all.
Perform all computations one after another.
"""
struct Serial <: DistributionStrategy end

function step!(new, old, ::Serial)
    m, n = size(old)
    I = (2:m-1,2:n-1)
    prepare!(old)
    life_rule!(new, old, I, I)
end


text:
"""
    PedestrianLatLonAccel
Pedestrian walking action. Acceleration in the Frenet frame, along with
desired lane after crossing the street.
"""
struct PedestrianLatLonAccel
    a_lat::Float64
    a_lon::Float64
    lane_des::Lane
end

function AutomotiveDrivingModels.propagate(ped::Vehicle,
                                           action::PedestrianLatLonAccel,
                                           roadway::Roadway,
                                           ΔT::Float64=0.1
                                           )
    state = propagate(ped, LatLonAccel(action.a_lat, action.a_lon), roadway, ΔT)
    roadproj = proj(state.posG, roadway[action.lane_des.tag], roadway, move_along_curves=false)
    retval = VehicleState(Frenet(roadproj, roadway), roadway, state.v)
    return retval
end

max_stars_repo_name: intelligent-control-lab/Autoenv

max_issues_repo_name: intelligent-control-lab/Autoenv

max_forks_repo_name: intelligent-control-lab/Autoenv


text:
using Revise, Test, ForwardDiff, Parameters, Setfield, Plots, LinearAlgebra
using BifurcationKit
const BK = BifurcationKit

norminf(x) = norm(x, Inf)
####################################################################################################
function TMvf!(dz, z, p, t)
	@unpack J, α, E0, τ, τD, τF, U0 = p
	E, x, u = z
	SS0 = J * u * x * E + E0
	SS1 = α * log(1 + exp(SS0 / α))
	[
		(-E + SS1) / τ
		(1.0 - x) / τD - u * x * E
		(U0 - u) / τF +  U0 * (1.0 - u) * E
	]
end

TMvf(z, p) = TMvf!(similar(z), z, p, 0)

# we group the differentials together
jet  = BK.getJet(TMvf, matrixfree=false)

par_tm = (α = 1.5, τ = 0.013, J = 3.07, E0 = -2.0, τD = 0.200, U0 = 0.3, τF = 1.5, τS = 0.007) #2.87
z0 = [0.238616, 0.982747, 0.367876 ]

opts_br = ContinuationPar(pMin = -10.0, pMax = -0.9, ds = 0.04, dsmax = 0.125, nInversion = 8, detectBifurcation = 3, maxBisectionSteps = 25, nev = 3)
	opts_br = @set opts_br.newtonOptions.verbose = false
	br, = continuation(jet[1], jet[2], z0, par_tm, (@lens _.E0), opts_br;
	printSolution = (x, p) -> (E = x[1], x = x[2], u = x[3]),
	tangentAlgo = BorderedPred(),
	plot = true, normC = norminf)

plot(br, plotfold=false, markersize=4, legend=:topleft)
####################################################################################################
hopfpt = computeNormalForm(jet..., br, 4)

# newton parameters
optn_po = NewtonPar(verbose = true, tol = 1e-8,  maxIter = 8)

# continuation parameters
opts_po_cont = ContinuationPar(dsmax = 0.1, ds= 0.0001, dsmin = 1e-4, pMax = 0., pMin=-5., maxSteps = 110, newtonOptions = (@set optn_po.tol = 1e-7), nev = 3, precisionStability = 1e-8, detectBifurcation = 0, plotEveryStep = 20, saveSolEveryStep=1)

# arguments for periodic orbits
args_po = (	recordFromSolution = (x, p) -> begin
		xtt = BK.getPeriodicOrbit(p.prob, x, @set par_tm.E0 = p.p)
		return (max = maximum(xtt[1,:]),
				min = minimum(xtt[1,:]),
				period = getPeriod(p.prob, x, @set par_tm.E0 = p.p))
	end,
	plotSolution = (x, p; k...) -> begin
		xtt = BK.getPeriodicOrbit(p.prob, x, @set par_tm.E0 = p.p)
		@show size(xtt[:,:]) maximum(xtt[1,:])
		plot!(xtt.t, xtt[1,:]; label = "E", k...)
		plot!(xtt.t, xtt[2,:]; label = "x", k...)
		plot!(xtt.t, xtt[3,:]; label = "u", k...)
		plot!(br; subplot = 1, putspecialptlegend = false)
		end,
	normC = norminf)

Mt = 200 # number of sections
	br_potrap, utrap = continuation(
	jet..., br, 4, opts_po_cont,
	PeriodicOrbitTrapProblem(M = Mt);
	jacobianPO = :Dense,
	verbosity = 2,	plot = true,
	args_po...,
	)

plot(br, br_potrap, markersize = 3)
	plot!(br_potrap.param, br_potrap.min, label = "")
####################################################################################################
# based on collocation
hopfpt = computeNormalForm(jet..., br, 4)

# newton parameters
optn_po = NewtonPar(verbose = true, tol = 1e-8,  maxIter = 25)

# continuation parameters
opts_po_cont = ContinuationPar(dsmax = 0.1, ds= -0.001, dsmin = 1e-4, pMax = 0., pMin=-5., maxSteps = 150, newtonOptions = (@set optn_po.tol = 1e-7), nev = 3, precisionStability = 1e-5, detectBifurcation = 0, plotEveryStep = 40, saveSolEveryStep=1)

br_pocoll, ucoll, = @time continuation(
	jet..., br, 4, opts_po_cont,
	PeriodicOrbitOCollProblem(20, 5);
	tangentAlgo = BorderedPred(),
	updateSectionEveryStep = 1,
	verbosity = 2,	plot = false,
	args_po...,
	plotSolution = (x, p; k...) -> begin
		xtt = BK.getPeriodicOrbit(p.prob, x, p.p)
		plot!(xtt.t, xtt[1,:]; label = "", marker =:d, markersize = 1.5, k...)
		plot!(br; subplot = 1, putspecialptlegend = false)

	end,
	callbackN = BK.cbMaxNorm(1000.),
	finaliseSolution = (z, tau, step, contResult; k...) -> begin
		prob = k[:prob]
		newt, err = BK.computeError(prob, z.u; verbosity = 3, par = BK.setParam(contResult, z.p))#, K = 100)
		return true
	end,

	)

plot(br, br_pocoll, markersize = 3)
	# plot!(br_pocoll.param, br_pocoll.min, label = "")
	# plot!(br, br_potrap, markersize = 3)
	# plot!(br_potrap.param, br_potrap.min, label = "", marker = :d)

####################################################################################################
# idem with Standard shooting
using DifferentialEquations#, TaylorIntegration

# this is the ODEProblem used with `DiffEqBase.solve`
probsh = ODEProblem(TMvf!, copy(z0), (0., 1000.), par_tm; abstol = 1e-10, reltol = 1e-9)

opts_po_cont = ContinuationPar(dsmax = 0.1, ds= -0.0001, dsmin = 1e-4, pMax = 0., pMin=-5., maxSteps = 120, newtonOptions = (@set optn_po.tol = 1e-6), nev = 3, precisionStability = 1e-8, detectBifurcation = 0, plotEveryStep = 10, saveSolEveryStep=0)

br_posh, = @time continuation(jet...,
	br, 4, opts_po_cont,
	# this is where we tell that we want Standard Shooting
	# with 15 time sections
	ShootingProblem(15, probsh, Rodas4(), parallel = true);
	# this to help branching
	δp = 0.0005,
	# deflation helps not converging to an equilibrium instead of a PO
	usedeflation = true,
	# this linear solver is specific to ODEs
	# it is computed using AD of the flow and
	# updated inplace
	jacobianPO = :autodiffDense,
	# we update the section along the branches
	updateSectionEveryStep = 2,
	# regular continuation parameters
	verbosity = 2,	plot = true,
	args_po...)

plot(br_posh, br, markersize=3)
	# plot(br, br_potrap, br_posh, markersize=3)
####################################################################################################
# idem with Poincaré shooting
@assert 1==0 "There is a PB with DE!!! I am getting NaN :("

function TMvfExtended!(dz, z, p, t)
	# we write the first part
	TMvf!(dz, z, p, t)
	dz[4:end] .= @views ForwardDiff.derivative(t -> TMvf(z[1:3] .+ t .* z[4:end], p), 0)
end

probmono = ODEProblem(TMvfExtended!, vcat(z0, z0), (0., 1000.), par_tm; abstol = 1e-10, reltol = 1e-9)

opts_po_cont = ContinuationPar(dsmax = 0.02, ds= 0.001, dsmin = 1e-6, pMax = 0., pMin=-5., maxSteps = 200, newtonOptions = NewtonPar(optn_po;tol = 1e-6, maxIter=15), nev = 3, precisionStability = 1e-8, detectBifurcation = 0, plotEveryStep = 1, saveSolEveryStep = 1)

br_popsh, = @time continuation(
	jet..., br, 4,
	# arguments for continuation
	opts_po_cont,
	# this is where we tell that we want Poincaré Shooting
	# PoincareShootingProblem(Mt, probsh, Rodas5(), probmono, Rodas4P(), parallel = false);
	PoincareShootingProblem(10, probsh, Rodas4P2(), parallel = false);
	ampfactor = 1.0, δp = 0.005,
	# usedeflation = true,
	jacobianPO = :autodiffDenseAnalytical,
	# jacobianPO = :autodiffDense,
	linearAlgo = MatrixBLS(),
	updateSectionEveryStep = 2,
	verbosity = 2,	plot = true,
	args_po...,
	callbackN = BK.cbMaxNorm(1e2),
	recordFromSolution = (x, p) -> (return (max = getMaximum(p.prob, x, @set par_tm.E0 = p.p), period = getPeriod(p.prob, x, @set par_tm.E0 = p.p))),
	normC = norminf)

plot(br, br_popsh, markersize=3)