text:
# Mutation operators
# ==================

# Isotropic mutation operator y' := y + σ(N_1(0, 1), ..., N_N(0, 1))
function isotropic(recombinant::T, s::S) where {T <: Vector, S <: Strategy}
    vals = randn(length(recombinant)) * s[:σ]
    recombinant += vals
    return recombinant
end

# Anisotropic mutation operator y' := y + σ(N_1(0, 1), ..., N_N(0, 1))
function anisotropic(recombinant::T, s::S) where {T <: Vector, S <: Strategy}
    @assert length(s[:σ]) == length(recombinant) "Sigma parameters must be defined for every dimension of objective parameter"
    vals = randn(length(recombinant)) .* s[:σ]
    recombinant += vals
    return recombinant
end

# Strategy mutation operators
# ===========================

# Isotropic strategy mutation σ' := σ exp(τ N(0, 1))
function isotropicSigma(s::S) where {S <: Strategy}
    @assert :σ ∈ keys(s) && :τ ∈ keys(s) "Strategy must have parameters: σ, τ"
    return strategy(σ = s[:σ] * exp(s[:τ]*randn()), τ = s[:τ])
end

# Anisotropic strategy mutation σ' := exp(τ0 N_0(0, 1))(σ_1 exp(τ N_1(0, 1)), ..., σ_N exp(τ N_N(0, 1)))
function anisotropicSigma(s::S) where {S <: Strategy}
    @assert :σ ∈ keys(s) && :τ ∈ keys(s) && :τ0 ∈ keys(s) "Strategy must have parameters: σ, τ0, τ"
    @assert isa(s[:σ], Vector) "Sigma must be a vector of parameters"
    #σ = exp(s[:τ0]*randn())*exp(s[:τ]*randn(length(s[:σ])))
    σ = exp.(s[:τ]*randn(length(s[:σ])))
    return strategy(σ = σ, τ = s[:τ], τ0 = s[:τ0])
end

# Genetic mutations
# =================

# Binary mutations
# ----------------
function flip(recombinant::Vector{Bool})
    s = length(recombinant)
    pos = rand(1:s)
    recombinant[pos] = !recombinant[pos]
    return recombinant
end

# Real valued mutation
# Mühlenbein, H. and Schlierkamp-Voosen, D.: Predictive Models for the Breeder Genetic Algorithm: I. Continuous Parameter Optimization. Evolutionary Computation, 1 (1), pp. 25-49, 1993.
# --------------------
function domainrange(valrange::Vector, m::Int = 20)
    prob = 1.0 / m
    function mutation(recombinant::T) where {T <: Vector}
        d = length(recombinant)
        @assert length(valrange) == d "Range matrix must have $(d) columns"
        δ = zeros(m)
        for i in 1:length(recombinant)
            for j in 1:m
                δ[j] = (rand() < prob) ? δ[j] = 2.0^(-j) : 0.0
            end
            if rand() > 0.5
                recombinant[i] += sum(δ)*valrange[i]
            else
                recombinant[i] -= sum(δ)*valrange[i]
            end
        end
        return recombinant
    end
    return mutation
end

# Combinatorial mutations (applicable to binary vectors)
# ------------------------------------------------------
function inversion(recombinant::T) where {T <: Vector}
    l = length(recombinant)
    from, to = rand(1:l, 2)
    from, to = from > to ? (to, from)  : (from, to)
    l = round(Int,(to - from)/2)
    for i in 0:(l-1)
        swap!(recombinant, from+i, to-i)
    end
    return recombinant
end

function insertion(recombinant::T) where {T <: Vector}
    l = length(recombinant)
    from, to = rand(1:l, 2)
    val = recombinant[from]
    deleteat!(recombinant, from)
    return insert!(recombinant, to, val)
end

function swap2(recombinant::T) where {T <: Vector}
    l = length(recombinant)
    from, to = rand(1:l, 2)
    swap!(recombinant, from, to)
    return recombinant
end

function scramble(recombinant::T) where {T <: Vector}
    l = length(recombinant)
    from, to = rand(1:l, 2)
    from, to = from > to ? (to, from)  : (from, to)
    diff = to - from + 1
    if diff > 1
        patch = recombinant[from:to]
        idx = randperm(diff)
        # println("$(from)-$(to), P: $(patch), I: $(idx)")
        for i in 1:(diff)
            recombinant[from+i-1] = patch[idx[i]]
        end
    end
    return recombinant
end

function shifting(recombinant::T) where {T <: Vector}
    l = length(recombinant)
    from, to, where = sort(rand(1:l, 3))
    patch = recombinant[from:to]
    diff = where - to
    if diff > 0
        # move values after tail of patch to the patch head position
        #println([from, to, where, diff])
        for i in 1:diff
            recombinant[from+i-1] = recombinant[to+i]
        end
        # place patch values in order
        start = from + diff
        for i in 1:length(patch)
            recombinant[start+i-1] = patch[i]
        end
    end
    return recombinant
end

# Utils
# =====
function swap!(v::T, from::Int, to::Int) where {T <: Vector}
    val = v[from]
    v[from] = v[to]
    v[to] = val
end

function mutationwrapper(gamutation::Function)
    wrapper(recombinant::T, s::S) where {T <: Vector, S <: Strategy} =  gamutation(recombinant)
    return wrapper
end


text:
# Test that the land model still runs, even with the lowest/simplest
# version of soil (prescribed heat and prescribed water - no state
# variables)
using MPI
using OrderedCollections
using StaticArrays
using Test

using CLIMAParameters
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

using ClimateMachine
using ClimateMachine.Land
using ClimateMachine.Land.SoilWaterParameterizations
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.DGMethods: BalanceLaw, LocalGeometry
using ClimateMachine.MPIStateArrays
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.VariableTemplates
using ClimateMachine.SingleStackUtils
using ClimateMachine.BalanceLaws:
    BalanceLaw, Prognostic, Auxiliary, Gradient, GradientFlux, vars_state

@testset "Prescribed Models" begin
    ClimateMachine.init()
    FT = Float64

    function init_soil_water!(land, state, aux, localgeo, time) end

    soil_water_model = PrescribedWaterModel()
    soil_heat_model = PrescribedTemperatureModel()
    soil_param_functions = nothing
    m_soil = SoilModel(soil_param_functions, soil_water_model, soil_heat_model)
    sources = ()
    m = LandModel(
        param_set,
        m_soil;
        source = sources,
        init_state_prognostic = init_soil_water!,
    )

    N_poly = 5
    nelem_vert = 10

    # Specify the domain boundaries
    zmax = FT(0)
    zmin = FT(-1)

    driver_config = ClimateMachine.SingleStackConfiguration(
        "LandModel",
        N_poly,
        nelem_vert,
        zmax,
        param_set,
        m;
        zmin = zmin,
        numerical_flux_first_order = CentralNumericalFluxFirstOrder(),
    )

    t0 = FT(0)
    timeend = FT(60)
    dt = FT(1)

    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        ode_dt = dt,
    )
    mygrid = solver_config.dg.grid
    Q = solver_config.Q
    aux = solver_config.dg.state_auxiliary

    ClimateMachine.invoke!(solver_config;)

    t = ODESolvers.gettime(solver_config.solver)
    state_vars = SingleStackUtils.get_vars_from_nodal_stack(
        mygrid,
        Q,
        vars_state(m, Prognostic(), FT),
    )
    aux_vars = SingleStackUtils.get_vars_from_nodal_stack(
        mygrid,
        aux,
        vars_state(m, Auxiliary(), FT),
    )
    #Make sure it runs, and that there are no state variables, and only "x,y,z" as aux.
    @test t == timeend
    @test size(Base.collect(keys(aux_vars)))[1] == 3
    @test size(Base.collect(keys(state_vars)))[1] == 0
end

ext: jl


text:
# Autogenerated wrapper script for JonkerVolgenant_jll for powerpc64le-linux-gnu
export bipartite_assignement

JLLWrappers.@generate_wrapper_header("JonkerVolgenant")
JLLWrappers.@declare_library_product(bipartite_assignement, "bipartite_assignement.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        bipartite_assignement,
        "lib/bipartite_assignement.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function approximation_term(
    ::MinimumBudget{SpanningTreeInstance{T, Maximise}, U},
    ::LagrangianAlgorithm,
) where {T, U}
    return NaN
end
function approximation_ratio(
    ::MinimumBudget{SpanningTreeInstance{T, Maximise}, U},
    ::LagrangianAlgorithm,
) where {T, U}
    return NaN
end

function _st_prim_budgeted_lagrangian(
    i::MinimumBudget{SpanningTreeInstance{T, Maximise}, U},
    λ::Float64,
) where {T, U}
    # Solve the subproblem for one value of the dual multiplier λ:
    #     l(λ) = \max_{x spanning tree} (rewards + λ weights) x - λ budget.
    sti_rewards = Dict{Edge{T}, Float64}(
        e => i.instance.rewards[e] + λ * i.weights[e] for
        e in keys(i.instance.rewards)
    )
    sti = SpanningTreeInstance(i.instance.graph, sti_rewards)
    sti_sol = solve(sti, PrimAlgorithm())
    return _budgeted_spanning_tree_compute_value(sti, sti_sol.variables) -
           λ * i.min_budget,
    sti_sol.variables
end

function solve(
    i::MinimumBudget{SpanningTreeInstance{T, Maximise}, U},
    ::LagrangianAlgorithm;
    ε::Float64,
) where {T, U}
    # Approximately solve the problem \min_{l ≥ 0} l(λ), where
    #     l(λ) = \max_{x spanning tree} (rewards + λ weights) x - λ budget.
    # This problem is the Lagrangian dual of the budgeted maximum spanning-tree problem:
    #     \max_{x spanning tree} rewards x  s.t.  weights x >= budget.
    # This algorithm provides no guarantee on the optimality of the solution.

    # Initial set of values for λ. The optimum is guaranteed to be contained in this interval.
    weights_norm_inf = maximum(values(i.weights)) # Maximum weight.
    m = nv(i.instance.graph) - 1 # Maximum number of items in a solution. Easy to compute for a spanning tree!
    λmax = Float64(weights_norm_inf * m + 1)

    λlow = 0.0
    λhigh = λmax

    # Perform a golden-ratio search.
    while (λhigh - λlow) > ε
        λmidlow = λhigh - (λhigh - λlow) / MathConstants.φ
        λmidhigh = λlow + (λhigh - λlow) / MathConstants.φ
        vmidlow, _ = _st_prim_budgeted_lagrangian(i, λmidlow)
        vmidhigh, _ = _st_prim_budgeted_lagrangian(i, λmidhigh)

        if vmidlow < vmidhigh
            λhigh = λmidhigh
            vhigh = vmidhigh
        else
            λlow = λmidlow
            vlow = vmidlow
        end
    end

    vlow, stlow = _st_prim_budgeted_lagrangian(i, λlow)
    vhigh, sthigh = _st_prim_budgeted_lagrangian(i, λhigh)
    if vlow < vhigh
        return BudgetedSpanningTreeLagrangianSolution(
            i,
            stlow,
            λlow,
            vlow,
            λmax,
        )
    else
        return BudgetedSpanningTreeLagrangianSolution(
            i,
            sthigh,
            λhigh,
            vhigh,
            λmax,
        )
    end
end


text:
using JudiLing
using CSV
using Test
using DataFrames

@testset "path finding for latin" begin
    latin_train =
        DataFrame(CSV.File(joinpath("data", "latin_mini.csv")))
    cue_obj_train = JudiLing.make_cue_matrix(
        latin_train,
        grams = 3,
        target_col = :Word,
        tokenized = false,
        keep_sep = false,
    )

    latin_val = latin_train[101:150, :]
    cue_obj_val = JudiLing.make_cue_matrix(
        latin_val,
        cue_obj_train,
        grams = 3,
        target_col = :Word,
        tokenized = false,
        keep_sep = false,
    )

    n_features = size(cue_obj_train.C, 2)

    S_train, S_val = JudiLing.make_S_matrix(
        latin_train,
        latin_val,
        ["Lexeme"],
        ["Person", "Number", "Tense", "Voice", "Mood"],
        ncol = n_features,
    )

    G_train = JudiLing.make_transform_matrix(S_train, cue_obj_train.C)

    Chat_train = S_train * G_train
    Chat_val = S_val * G_train

    F_train = JudiLing.make_transform_matrix(cue_obj_train.C, S_train)

    Shat_train = cue_obj_train.C * F_train
    Shat_val = cue_obj_val.C * F_train

    A = cue_obj_train.A

    max_t = JudiLing.cal_max_timestep(latin_train, latin_val, :Word)

    res_train, gpi_train = JudiLing.learn_paths(
        latin_train,
        latin_train,
        cue_obj_train.C,
        S_train,
        F_train,
        Chat_train,
        A,
        cue_obj_train.i2f,
        cue_obj_train.f2i,
        gold_ind = cue_obj_train.gold_ind,
        Shat_val = Shat_train,
        check_gold_path = true,
        max_t = max_t,
        max_can = 10,
        grams = 3,
        threshold = 0.1,
        tokenized = false,
        sep_token = "_",
        keep_sep = false,
        target_col = :Word,
        issparse = :dense,
        verbose = false,
    )

    res_val, gpi_val = JudiLing.learn_paths(
        latin_train,
        latin_val,
        cue_obj_train.C,
        S_val,
        F_train,
        Chat_val,
        A,
        cue_obj_train.i2f,
        cue_obj_train.f2i,
        gold_ind = cue_obj_val.gold_ind,
        Shat_val = Shat_val,
        check_gold_path = true,
        max_t = max_t,
        max_can = 10,
        grams = 3,
        threshold = 0.1,
        is_tolerant = true,
        tolerance = 0.1,
        max_tolerance = 0,
        tokenized = false,
        sep_token = "-",
        keep_sep = false,
        target_col = :Word,
        issparse = :dense,
        verbose = false,
    )

    res_train = JudiLing.learn_paths(
        latin_train,
        latin_train,
        cue_obj_train.C,
        S_train,
        F_train,
        Chat_train,
        A,
        cue_obj_train.i2f,
        cue_obj_train.f2i,
        max_t = max_t,
        max_can = 10,
        grams = 3,
        threshold = 0.1,
        tokenized = false,
        sep_token = "_",
        keep_sep = false,
        target_col = :Word,
        issparse = :dense,
        verbose = false,
    )

    res_val = JudiLing.learn_paths(
        latin_train,
        latin_val,
        cue_obj_train.C,
        S_val,
        F_train,
        Chat_val,
        A,
        cue_obj_train.i2f,
        cue_obj_train.f2i,
        max_t = max_t,
        max_can = 10,
        grams = 3,
        threshold = 0.1,
        is_tolerant = true,
        tolerance = 0.1,
        max_tolerance = 0,
        tokenized = false,
        sep_token = "-",
        keep_sep = false,
        target_col = :Word,
        issparse = :dense,
        verbose = false,
    )

    res_train = JudiLing.learn_paths(
        latin_train,
        cue_obj_train,
        S_train,
        F_train,
        Chat_train,
        verbose = false,
    )

    res_train, ts_train = JudiLing.learn_paths(
        latin_train,
        latin_train,
        cue_obj_train.C,
        S_train,
        F_train,
        Chat_train,
        A,
        cue_obj_train.i2f,
        cue_obj_train.f2i,
        max_t = max_t,
        max_can = 10,
        grams = 3,
        threshold = 0.1,
        tokenized = false,
        sep_token = "_",
        keep_sep = false,
        target_col = :Word,
        issparse = :dense,
        check_threshold_stat = true,
        verbose = false,
    )

    res_val, ts_val = JudiLing.learn_paths(
        latin_train,
        latin_val,
        cue_obj_train.C,
        S_val,
        F_train,
        Chat_val,
        A,
        cue_obj_train.i2f,
        cue_obj_train.f2i,
        max_t = max_t,
        max_can = 10,
        grams = 3,
        threshold = 0.1,
        is_tolerant = true,
        tolerance = 0.1,
        max_tolerance = 0,
        tokenized = false,
        sep_token = "-",
        keep_sep = false,
        target_col = :Word,
        issparse = :dense,
        check_threshold_stat = true,
        verbose = false,
    )

    res_train, gpi_train, ts_train = JudiLing.learn_paths(
        latin_train,
        latin_train,
        cue_obj_train.C,
        S_train,
        F_train,
        Chat_train,
        A,
        cue_obj_train.i2f,
        cue_obj_train.f2i,
        gold_ind = cue_obj_train.gold_ind,
        Shat_val = Shat_train,
        check_gold_path = true,
        max_t = max_t,
        max_can = 10,
        grams = 3,
        threshold = 0.1,
        tokenized = false,
        sep_token = "_",
        keep_sep = false,
        target_col = :Word,
        issparse = :dense,
        check_threshold_stat = true,
        verbose = false,
    )

    res_val, gpi_val, ts_val = JudiLing.learn_paths(
        latin_train,
        latin_val,
        cue_obj_train.C,
        S_val,
        F_train,
        Chat_val,
        A,
        cue_obj_train.i2f,
        cue_obj_train.f2i,
        gold_ind = cue_obj_val.gold_ind,
        Shat_val = Shat_val,
        check_gold_path = true,
        max_t = max_t,
        max_can = 10,
        grams = 3,
        threshold = 0.1,
        is_tolerant = true,
        tolerance = 0.1,
        max_tolerance = 0,
        tokenized = false,
        sep_token = "-",
        keep_sep = false,
        target_col = :Word,
        issparse = :dense,
        check_threshold_stat = true,
        verbose = false,
    )
end


text:
## C3V Symmetry
# Symmetric stretch
"""
Susceptibility for tensor elements xxz and yyz.
"""
function susceptibility_c3v_ss_xxz(θ, R; β_ccc=1.0, N=1.0)
    1/2 * N * β_ccc * ( (1+R) * cos(θ) - (1-R) * cos(θ)^3 )
end

"""
Susceptibility for tensor elements xzx, zxx, yzy and zyy.
"""
function susceptibility_c3v_ss_xzx(θ, R; β_ccc=1.0, N=1.0)
    1/2 * N * β_ccc * (1-R) * ( cos(θ) - cos(θ)^3 )
end

"""
Susceptibility for tensor element zzz
"""
function susceptibility_c3v_ss_zzz(θ, R; β_ccc=1.0, N=1.0)
    N * β_ccc * ( R * cos(θ) + (1-R) * cos(θ)^3 )
end

# Asymmetric stretch
"""
Susceptibility for tensor elements xxz and yyz.
"""
function susceptibility_c3v_as_xxz(θ; β_aca=1.0, N=1.0)
    -N * β_aca * (cos(θ) - cos(θ)^3)
end

"""
Susceptibility for tensor elements xzx, zxx, yzy and zyy.
"""
function susceptibility_c3v_as_xzx(θ; β_aca=1.0, N=1.0)
    N * β_aca * cos(θ)^3
end

"""
Susceptibility for tensor element zzz
"""
function susceptibility_c3v_as_zzz(θ; β_aca=1.0, N=1.0)
    2N * β_aca * (cos(θ) - cos(θ)^3)
end

"""
    susceptibility(θ::Real, t::Symbol; pointgroup, mode, R, beta, N; <kwargs>)

Compute the susceptibility for angle `θ` and tensor element `t`.
Tensor elements important for specific polarization combinations:
* ssp: `:yyz`
* sps: `:yzy`
* pss: `:zyy`
* ppp: `:xxz`, `:xzx`, `:zxx` and `:zzz`

# Arguments
* `θ`: title angle of the moiety
* `t`: tensor element (see above)

# Keyword Arguments
* `pointgroup::Symbol`: point group of the moiety (possible values: `:c3v`)
* `mode::Symbol`: vibrational mode (possible values: `:ss` (symmetric stretch), `:as` (asymmetric stretch))
* `R::Real`: the hyperpolarizability ratio \$R = \\beta_{aac} / \\beta_{ccc} = \\beta_{bbc} / \\beta_{ccc}\$. Has only to be provided for point group `:c3v` and mode `:ss`.
* `beta::Real`: hyperpolarizability
    * ``\\beta_{ccc}`` for `:c3v` and `:ss`
    * ``\\beta_{aca} = \\beta_{bcb}`` for `:c3v` and `:as`
* `N`::Real: number of oscillators

| Molecule                  	| Group  	| R   	|
|:--------------------------	|:-------	|:----	|
| methanol                  	| methyl 	| 1.7 	|
| ethanol and longer X-1-ol 	| methyl 	| 3.4 	|
| acetone                   	| methyl 	| 1.9 	|
| DMSO                      	| methyl 	| 2.3 	|
"""
function susceptibility(θ::Real, t::Symbol;
        pointgroup::Union{Nothing,Symbol}=nothing, mode::Union{Nothing,Symbol}=nothing,
        R::Union{Nothing,Real}=nothing, beta::Real=1.0, N::Real=1.0)
    # Run checks
    isnothing(pointgroup) && error("no pointgroup provided")
    isnothing(mode) && error("no mode provided")
    pointgroup in [:c3v] || error("unknown point group $pointgroup")
    mode in [:ss, :as] || error("unknown mode $mode")
    t in [:xxz, :yyz, :xzx, :zxx, :yzy, :zyy, :zzz] || error("unknown tensor element $t")
    # For the point group of c3v we have to provide an R value
    pointgroup in [:c3v] && mode == :ss && isnothing(R) && error("You have to provide the keyword argument R.")

    if pointgroup == :c3v
        if mode == :ss
            if t in [:xxz, :yyz]
                χ = susceptibility_c3v_ss_xxz(θ, R; β_ccc=beta, N=N)
            elseif t in [:xzx, :zxx, :yzy, :zyy]
                χ = susceptibility_c3v_ss_xzx(θ, R; β_ccc=beta, N=N)
            elseif t in [:zzz]
                χ = susceptibility_c3v_ss_zzz(θ, R; β_ccc=beta, N=N)
            end
        elseif mode == :as
            if t in [:xxz, :yyz]
                χ = susceptibility_c3v_as_xxz(θ; β_aca=beta, N=N)
            elseif t in [:xzx, :zxx, :yzy, :zyy]
                χ = susceptibility_c3v_as_xzx(θ; β_aca=beta, N=N)
            elseif t in [:zzz]
                χ = susceptibility_c3v_as_zzz(θ; β_aca=beta, N=N)
            end
        end
    end
    return χ
end

"""
    effective_susceptibility_ssp(n1, n2, n′, β, χ_yyz)

Compute the effective susceptibility in ssp polarization.

# Arguments
* `n1`: refractive indices of the medium in which the waves are propagating towards the interface
* `n2`: refractive indices of the other medium
* `n′`: refractive index of the interface for \$\\omega_2\$
* `β`:  incident or reflected angles
* `χ_yyz`: susceptibility

The arguments `n1`, `n2`, and `β` have to be arrays where the first element
corresponds to the property for the sum frequency, the second to \$\\omega_1\$
and the third to \$\\omega_2\$.

https://doi.org/10.1080/01442350500225894
"""
function effective_susceptibility_ssp(
        n1::Array{<:Number},
        n2::Array{<:Number},
        n′::Array{<:Number},
        β::Array{<:Number},
        χ_yyz::Number
    )
    χ2_eff = fresnel_y(n1[1], n2[1], β[1]) *
             fresnel_y(n1[2], n2[2], β[2]) *
             fresnel_z(n1[3], n2[3], n′[3], β[3]) *
             sin(β[3]) * χ_yyz
end

"""
    effective_susceptibility_ssp(s::SFGAnalysis.Setup, χ_yyz)

Compute the effective susceptibility in ssp polarization.

# Arguments
* `s`: setup object
* `χ_yyz`: susceptibility
"""
function effective_susceptibility_ssp(s::Setup, χ_yyz::Number)
    n1 = [s.Ω.n1, s.Ω1.n1, s.Ω2.n1]
    n2 = [s.Ω.n2, s.Ω1.n2, s.Ω2.n2]
    n′ = [s.Ω.n′, s.Ω1.n′, s.Ω2.n′]
    β = [s.Ω.β, s.Ω1.β, s.Ω2.β]
    effective_susceptibility_ssp(n1, n2, n′, β, χ_yyz)
end

"""
    effective_susceptibility_sps(n1, n2, n′, β, χ_yyz)

Compute the effective susceptibility in sps polarization.

# Arguments
* `n1`: refractive indices of the medium in which the waves are propagating towards the interface
* `n2`: refractive indices of the other medium
* `n′`: refractive index of the interface for \$\\omega_2\$
* `β`:  incident or reflected angles
* `χ_yyz`: susceptibility

The arguments `n1`, `n2`, and `β` have to be arrays where the first element
corresponds to the property for the sum frequency, the second to \$\\omega_1\$
and the third to \$\\omega_2\$.

https://doi.org/10.1080/01442350500225894
"""
function effective_susceptibility_sps(
        n1::Array{<:Number},
        n2::Array{<:Number},
        n′::Array{<:Number},
        β::Array{<:Number},
        χ_yzy::Number
    )
    χ2_eff = fresnel_y(n1[1], n2[1], β[1]) *
             fresnel_z(n1[2], n2[2], n′[2], β[2]) *
             fresnel_y(n1[3], n2[3], β[3]) *
             sin(β[2]) * χ_yzy
end

"""
    effective_susceptibility_sps(s::SFGAnalysis.Setup, χ_yzy)

Compute the effective susceptibility in sps polarization.

# Arguments
* `s`: setup object
* `χ_yzy`: susceptibility
"""
function effective_susceptibility_sps(s::Setup, χ_yyz::Number)
    n1 = [s.Ω.n1, s.Ω1.n1, s.Ω2.n1]
    n2 = [s.Ω.n2, s.Ω1.n2, s.Ω2.n2]
    n′ = [s.Ω.n′, s.Ω1.n′, s.Ω2.n′]
    β = [s.Ω.β, s.Ω1.β, s.Ω2.β]
    effective_susceptibility_sps(n1, n2, n′, β, χ_yyz)
end

# pss
"""
    effective_susceptibility_pss(n1, n2, n′, β, χ_yyz)

Compute the effective susceptibility in pss polarization.

# Arguments
* `n1`: refractive indices of the medium in which the waves are propagating towards the interface
* `n2`: refractive indices of the other medium
* `n′`: refractive index of the interface for \$\\omega_2\$
* `β`:  incident or reflected angles
* `χ_yyz`: susceptibility

The arguments `n1`, `n2`, and `β` have to be arrays where the first element
corresponds to the property for the sum frequency, the second to \$\\omega_1\$
and the third to \$\\omega_2\$.

https://doi.org/10.1080/01442350500225894
"""
function effective_susceptibility_pss(
        n1::Array{<:Number},
        n2::Array{<:Number},
        n′::Array{<:Number},
        β::Array{<:Number},
        χ_zyy::Number
    )
    χ2_eff = fresnel_z(n1[1], n2[1], n′[1], β[1]) *
             fresnel_y(n1[2], n2[2], β[2]) *
             fresnel_y(n1[3], n2[3], β[3]) *
             sin(β[1]) * χ_zyy
end

"""
    effective_susceptibility_pss(s::SFGAnalysis.Setup, χ_zyy)

Compute the effective susceptibility in pss polarization.

# Arguments
* `s`: setup object
* `χ_zyy`: susceptibility
"""
function effective_susceptibility_pss(s::Setup, χ_yyz::Number)
    n1 = [s.Ω.n1, s.Ω1.n1, s.Ω2.n1]
    n2 = [s.Ω.n2, s.Ω1.n2, s.Ω2.n2]
    n′ = [s.Ω.n′, s.Ω1.n′, s.Ω2.n′]
    β = [s.Ω.β, s.Ω1.β, s.Ω2.β]
    effective_susceptibility_pss(n1, n2, n′, β, χ_yyz)
end

"""
    effective_susceptibility_ppp(n1, n2, n′, β, χ)

Compute the effective susceptibility in ppp polarization.

# Arguments
* `n1`: refractive indices of the medium in which the waves are propagating towards the interface
* `n2`: refractive indices of the other medium
* `n′`: refractive indices of the interface
* `β`:  incident or reflected angles
* `χ`: susceptibilities as an array `[χ_xxz, χ_xzx, χ_zxx, χ_zzz]`

The arguments `n1`, `n2`, `n′` and `β` have to be arrays where the first element
corresponds to the property for the sum frequency, the second to \$\\omega_1\$
and the third to \$\\omega_2\$.

https://doi.org/10.1080/01442350500225894
"""
function effective_susceptibility_ppp(n1::Array{<:Number},
                n2::Array{<:Number}, n′::Array{<:Number},
                β::Array{<:Number}, χ::Array{<:Number})
    χ_xxz, χ_xzx, χ_zxx, χ_zzz = χ
    χ2_eff = -fresnel_x(n1[1], n2[1], β[1]) *
             fresnel_x(n1[2], n2[2], β[2]) *
             fresnel_z(n1[3], n2[3], n′[3], β[3]) *
             cos(β[1]) * cos(β[2]) * sin(β[3]) * χ_xxz -
             fresnel_x(n1[1], n2[1], β[1]) *
             fresnel_z(n1[2], n2[2], n′[2], β[2]) *
             fresnel_x(n1[3], n2[3], β[3]) *
             cos(β[1]) * sin(β[2]) * cos(β[3]) * χ_xzx +
             fresnel_z(n1[1], n2[1], n′[1], β[1]) *
             fresnel_x(n1[2], n2[2], β[2]) *
             fresnel_x(n1[3], n2[3], β[3]) *
             sin(β[1]) * cos(β[2]) * cos(β[3]) * χ_zxx +
             fresnel_z(n1[1], n2[1], n′[1], β[1]) *
             fresnel_z(n1[2], n2[2], n′[2], β[2]) *
             fresnel_z(n1[3], n2[3], n′[3], β[3]) *
             sin(β[1]) * sin(β[2]) * sin(β[3]) * χ_zzz
end

"""
    effective_susceptibility_ppp(s::SFGAnalysis.Setup, χ)

# Arguments
* `s`: setup object
* `χ`: susceptibilities as an array `[χ_xxz, χ_xzx, χ_zxx, χ_zzz]`
"""
function effective_susceptibility_ppp(s::Setup, χ::Array{<:Number})
    n1 = [s.Ω.n1, s.Ω1.n1, s.Ω2.n1]
    n2 = [s.Ω.n2, s.Ω1.n2, s.Ω2.n2]
    n′ = [s.Ω.n′, s.Ω1.n′, s.Ω2.n′]
    β = [s.Ω.β, s.Ω1.β, s.Ω2.β]
    effective_susceptibility_ppp(n1, n2, n′, β, χ)
end

"""
    effective_susceptibility(θ::Real, s::Setup, p::Symbol; <kwargs>)

Compute the effective susceptibility for a moiety with tilt angle `θ` for a `Setup` `s`
in polarization `p`.

# Arguments
* `θ`: tilt angle of the moiety
* `s`: `Setup` object
* `p`: polarization combination (`:ppp` or `:ssp`)

# Keyword Arguments
* `pointgroup::Symbol`: point group of the moiety (possible values: `:c3v`)
* `mode::Symbol`: vibrational mode (possible values: `:ss` (symmetric stretch), `:as` (asymmetric stretch))
* `R::Real`: the hyperpolarizability ratio \$R = \\beta_{aac} / \\beta_{ccc} = \\beta_{bbc} / \\beta_{ccc}\$. Has only to be provided for point group `:c3v` and mode `:ss`.
* `beta::Real`: hyperpolarizability
    * ``\\beta_{ccc}`` for `:c3v` and `:ss`
    * ``\\beta_{aca} = \\beta_{bcb}`` for `:c3v` and `:as`
* `N::Real`: number of oscillators

"""
function effective_susceptibility(θ::Real, s::Setup, p::Symbol;
    pointgroup::Union{Nothing,Symbol}=nothing,
    mode::Union{Nothing,Symbol}=nothing,
    R::Union{Nothing,Real}=nothing,
    beta::Real=1.0, N::Real=1.0)

    p in [:ssp, :ppp, :sps, :pss] || error("unknown polarization combination $p")

    if p == :ppp
        χ = susceptibility.(θ, [:xxz, :xzx, :zxx, :zzz];
            pointgroup=pointgroup,
            mode=mode, R=R,
            beta=beta, N=N
        )
        χ_eff = effective_susceptibility_ppp(s, χ)
    elseif p == :ssp
        χ = susceptibility(θ, :yyz;
            pointgroup=pointgroup,
            mode=mode, R=R,
            beta=beta, N=N
        )
        χ_eff = effective_susceptibility_ssp(s, χ)
    elseif p == :sps
        χ = susceptibility(θ, :yzy;
            pointgroup=pointgroup,
            mode=mode, R=R,
            beta=beta, N=N
        )
        χ_eff = effective_susceptibility_sps(s, χ)
    elseif p == :pss
        χ = susceptibility(θ, :zyy;
            pointgroup=pointgroup,
            mode=mode, R=R,
            beta=beta, N=N
        )
        χ_eff = effective_susceptibility_pss(s, χ)
    end
    χ_eff
end

"""
    effective_susceptibility(d::Distribution, s::Setup, p::Symbol; <kwargs>)

Same as above. In this case the tilt angle `θ` is expressed by a `Distribution` `d`.
"""
function effective_susceptibility(d::Distribution, s::Setup, p::Symbol;
    pointgroup::Union{Nothing,Symbol}=nothing,
    mode::Union{Nothing,Symbol}=nothing,
    R::Union{Nothing,Real}=nothing,
    beta::Real=1.0, N::Real=1.0,
    samples::Int=360)

    sus = θ -> effective_susceptibility(θ, s, p; pointgroup=pointgroup, mode=mode, R=R, beta=beta, N=N)
    integral, error = quadgk(θ -> sus(θ) * pdf(d, θ), -π, π, rtol=1e-8)
    integral
end


text:
using GapFindFill, Test, Gurobi
include("loadGAMS.jl")

println("======Tesing Gap Finding===")
println("======loading data===")
#*************SETS********************************
reactionList = read_GAMS_format_file("gapfind/reactions.txt")
@assert(length(reactionList) == 2383)
# foreach(println, reactionList)
revRxnList = read_GAMS_format_file("gapfind/reversible_reactions.txt")
@assert(length(revRxnList) == 852)
# foreach(println, revRxnList)
compoundList = read_GAMS_format_file("gapfind/compounds.txt")
@assert(length(compoundList) == 1668)
# foreach(println, compoundList)
cytCompoundList = read_GAMS_format_file("gapfind/cytosolic_compounds.txt")
@assert(length(cytCompoundList) == 951)
# foreach(println, cytCompoundList)
extCompoundList = read_GAMS_format_file("gapfind/extracellular_compounds.txt")
@assert(length(extCompoundList) == 299)
# foreach(println, extCompoundList)

#*****************PARAMETERS*****************************
fluxUBList = read_GAMS_format_file("gapfind/upperbound_on_fluxes.txt")
@assert(length(fluxUBList) == 2383)
# foreach(println, fluxUBList)
fluxLBList = read_GAMS_format_file("gapfind/lowerbound_on_fluxes.txt")
@assert(length(fluxLBList) == 2383)
# foreach(println, fluxLBList)
stoiList = read_GAMS_format_file("gapfind/S_matrix.txt")
@assert(length(stoiList) == 9329)
# foreach(println, stoiList)

#****************PARAMETERS PREPROCESSING***************************
# map rxn_name to index
rxnDict = Dict{String, Integer}()
for (id, name) in enumerate(reactionList)
  rxnDict[name] = id
end
# map compound_name to index
compoundDict = Dict{String, Integer}()
for (id, name) in enumerate(compoundList)
  compoundDict[name] = id
end
noRxn = length(reactionList)
noCompound = length(compoundList)
# isRev -- isRev[i] is true if rxn in rxnDict(rxn, i) is reversible
isRev = zeros(Bool, noRxn)
for el in revRxnList
  isRev[rxnDict[el]] = true
end
@assert(length(revRxnList) == sum(isRev))
# isCyt -- isCyt[i] is true if cp in compoundDict(cp, i) is in cytosol
isCyt = zeros(Bool, noCompound)
for el in cytCompoundList
  isCyt[compoundDict[el]] = true
end
@assert(length(cytCompoundList) == sum(isCyt))
# isExt -- isExt[i] is true if cp in compoundDict(cp, i) is in extracellular
isExt = zeros(Bool, noCompound)
for el in extCompoundList
  isExt[compoundDict[el]] = true
end
@assert(length(extCompoundList) == sum(isExt))
# get S matrix from stoiList
stoiMatrix = zeros(Float64, noCompound, noRxn)
for elm in stoiList
  stoiMatrix[compoundDict[elm[1]],rxnDict[elm[2]]] = elm[3]
end
println("size of S matrix: ", size(stoiMatrix))
# get upperbound on fluxes
fluxUB = zeros(Float64, noRxn)
for (name, val) in fluxUBList
  fluxUB[rxnDict[name]] = val
end
println("big-M in upperbound: ", maximum(fluxUB))
# get lowerbound on fluxes
fluxLB = zeros(Float64, noRxn)
for (name, val) in fluxLBList
  fluxLB[rxnDict[name]] = val
end
println("big-M in lowerbound: ", minimum(fluxLB))

println("======calling find_gaps===")
#*********Calling function*********************
model, objVal, status, solveTime, xValue = find_gaps(isRev, isCyt, isExt,
  stoiMatrix, fluxLB, fluxUB, solver=Gurobi)
println("======results report===")
println("Best objective: $(objVal)")
println("Solve Time by \"time()\": $(solveTime)")
println("Solve Status: $(status)")

#***************VERIFICATION*********************
println("======VERIFICATION===")
stdResultList = read_GAMS_format_file("gapfind/sample_results.txt")
stdResult = Set(stdResultList)
@assert(length(stdResult) == (noCompound - objVal))
println(sum(xValue))
@assert(objVal == sum(xValue))
try
  for (id, val) in enumerate(xValue)
    if val == 0
      @assert(in(compoundList[id], stdResult))
    end
  end
  println("Consistent results")
  println("======Gap Finding Test Passed===")
catch e
  println("Inconsistent results")
end


text:
module JuliaUtils

using LinearAlgebra
using ForwardDiff

export linearls, predict
export newton

include("LinearLS.jl")
include("newton.jl")

end


text:
#For running on a server, need to flush out the print calls
function message(x...)
    println(x...)
    flush(stdout)
end

####################################################################
# Data Structures
####################################################################

#Collection of variables to pass along to different functions in the FGES algorithm.
struct ParseData{M, F<:AbstractFloat}
    data::M #orginal data
    normAugData::M #standardized by the mean and std of each column, appended with ones column at end
    scatterMat::M # (covariance matrix not scaled by the number of observations)
    numFeatures::Int #number of columns
    numObservations::Int #number of rows
    penalty::F

    function ParseData(data::M, scatterMat, penalty) where M

        #Determine the data type of the data inputted
        baseType = eltype(data)

        #Get the dimensions of the input data
        numObservations, numFeatures = size(data)

        #Copy the data and standardize each column
        normAugData = copy(data)
        normAugData .-= mean(normAugData, dims=1) #subtract mean
        normAugData ./= std(normAugData, dims=1) #divide by standard deviation

        #Augment a column of ones on the end for linear regression
        normAugData = [normAugData ones(baseType, numObservations)]

        #If no scatter matrix is provided, calculate in construction
        if isnothing(scatterMat)
            scatterMat = normAugData'normAugData
        end

        #Ensure that the penalty is the same type as the data (e.g. Float32)
        penalty = baseType(penalty)

        return new{M, baseType}(data, normAugData, scatterMat, numFeatures, numObservations, penalty)
    end
end

#Simple structure to hold the current edge, a subset of neighbors, and a score change
Base.@kwdef mutable struct Step{A,B}
    edge::Edge = Edge(1,2,true)
    subset::Vector{A} = A[]
    Δscore::B = zero(B)
end

####################################################################
# Base Function Overloads
####################################################################

#Print method to display the parsed data
function show(io::IO, d::ParseData{T}) where T
    print(io, "$(d.numObservations)×$(d.numFeatures) ParseData{$(T)}")
end

#Print method to display the Step
function show(io::IO, newStep::Step{A,B}) where {A,B}
    print(io, "Edge: $(newStep.edge), Subset: $(newStep.subset), Δscore: $(newStep.Δscore)")
end

#This could be dangerous. The @memoize macro has to check if ParseData is the same argument. We only ever define one immutable ParseData, so it will never change. Any equality check should be true.
==(a::T, b::T) where T <: ParseData = true

####################################################################
# Main Entry point for the Algorithm
####################################################################

"""
    fges(data; debug=false)
Compute a causal graph for the given observed data.
"""
function fges(data; penalty = 1.0, scatterMat=nothing, debug=false)

    #Parse the inputted data
    dataParsed = ParseData(data, scatterMat, penalty)

    #Create an empty graph with one node for each feature
    g = PDAG(dataParsed.numFeatures)

    debug && message("Start forward search")
    #Perform the forward search
    Search!(g, dataParsed, Insert!, debug)

    debug && message("Start backward search")
    #Perform the backward search
    Search!(g, dataParsed, Delete!, debug)

    #Return the graph
    return g
end

####################################################################
# Insert and Delete Operators
####################################################################

function Insert!(g::PDAG, newStep::Step)
    edge = newStep.edge
    T = newStep.subset

    #Add a directed edge x→y
    addedge!(g,edge)

    #Orient all edges incident into child node
    y = edge.child
    for t ∈ T
        orientedge!(g,t,y) #t→y
    end

    return nothing
end

function Delete!(g::PDAG, newStep::Step)
    edge = newStep.edge
    H = newStep.subset

    #remove directed and unidrected edges (x→y and x-y)
    remedge!(g, edge)

    #Orient all vertices in H toward x and y
    x,y = edge.parent, edge.child
    for h ∈ H
        orientedge!(g,x,h) #x→h
        orientedge!(g,y,h) #y→h
    end

    return nothing
end

####################################################################
# Forward and Backward search
####################################################################

function Search!(g, dataParsed::ParseData{Matrix{A}}, operator, debug) where A

    #Create a container to hold information about the next step
    newStep = Step{Int,A}()

    #Continually add/remove edges to the graph until the score stops increasing
    while true

        #Get the new best step (depends on if we're inserting or deleting)
        if operator == Insert!
            findNextEquivClass!(newStep, dataParsed, g, findBestInsert, debug)
        else
            findNextEquivClass!(newStep, dataParsed, g, findBestDelete, debug)
        end

        #If the score did not improve...
        if newStep.Δscore ≤ zero(A)
            #...stop searching
            break
        end

        #Use the insert or delete operator update the graph
        operator(g, newStep)

        #Covert the PDAG to a complete PDAG
        #undirect all edges unless they participate in a v-structure
        graphVStructure!(g)
        #Apply the 4 Meek rules to orient some edges in the graph
        meekRules!(g)

        #Reset the score
        newStep.Δscore = zero(A)

    end

    return nothing
end

function findNextEquivClass!(newStep, dataParsed, g, findBestOperation::Function, debug)

    #Parallelization with synced threads (to avoid race conditions)
    @sync begin
        #Loop through all possible node combinations
        for (x,y) in allpairs(vertices(g))
            #Spawn a new thread
            Threads.@spawn begin
                #Check if there is an edge between two vertices
                hasEdge = isadjacent(g,x,y)

                #Skip over adjacent edges if we're trying to insert and vice versa
                if findBestOperation == findBestInsert ? !hasEdge : hasEdge
                    #For this pair of nodes, the following function will:
                        #(1) check if a valid operator exists
                        #(2) score all valid operators and return the one with the highest score
                        #findBestOperation is either "findBestInsert" or "findBestDelete"
                    (bestSubset, bestScore) = findBestOperation(dataParsed, g, x, y, debug)

                    #If the best valid operator was better than any previously found...
                    if bestScore > newStep.Δscore
                        #...update newStep
                        newStep.edge = Edge(x,y,true)
                        newStep.subset = bestSubset
                        newStep.Δscore = bestScore

                        debug && message(newStep)
                    end
                end

            end
        end
    end
end

function findBestInsert(dataParsed::ParseData{Matrix{A}}, g, x, y, debug) where A
    #Calculate two (possibly empty) sets of nodes
    # NAxy: any nodes that are undirected neighbors of y and connected to x by any edge
    # Txy: any subset of the undirected neighbors of y not connected to x
    Tyx = calcT(g, Edge(x,y,true))
    NAyx = calcNAyx(g, Edge(x,y,true))

    #Ceate two containers to hold the best found score and best subset of Tyx
    bestScore = zero(A)
    bestT = Vector{Int}()

    #Keep a list of invalid sets
    invalid = Vector{Vector{Int}}()

    #Loop through all possible subsets of Tyx
    for T in powerset(Tyx)
        if checkSupersets(T,invalid)
            NAyxT = NAyx ∪ T
            if isclique(g, NAyxT) && isblocked(g, y, x, NAyxT)

                #Score the valid Insert
                PAy = parents(g,y)
                PAy⁺ = PAy ∪ x
                newScore = score(dataParsed, NAyxT ∪ PAy⁺, y, debug) - score(dataParsed, NAyxT ∪ PAy, y, debug)

                #Save the new score if it was better than any previous
                if newScore > bestScore
                    bestT = T
                    bestScore = newScore
                end
            end
        else
            #Record that the subset T is invalid
            push!(invalid,T)
        end
    end

    return (bestT, bestScore)
end

#Check if the set T is contained in any invalid set
function checkSupersets(T,invalid)
    for i ∈ invalid
        if i ⊆ T
            return false
        end
    end
    return true
end

function findBestDelete(dataParsed::ParseData{Matrix{A}}, g, x, y, debug) where A
    #Calculate two (possibly empty) sets of nodes
    # NAxy: any nodes that are undirected neighbors of y and connected to x by any edge
    # Hyx: any subset of the undirected neighbors of y that are connected to x
    NAyx = calcNAyx(g, Edge(x,y,true))
    Hyx = NAyx

    #Ceate two containers to hold the best found score and best subset of Tyx
    bestScore = zero(A)
    bestH = Vector{Int}()

    #Loop through all possible subsets of Hyx
    for H in powerset(Hyx)
        #Calculate NAyx \ {H}
        NAyx_H = setdiff(NAyx,H)

        #Check if the operator is valid
        if isclique(g, NAyx_H)

            #Score the valid operator
            PAy = parents(g,y)
            PAy⁻ = setdiff(PAy, x)
            newScore = score(dataParsed, NAyx_H ∪ PAy⁻, y, debug) - score(dataParsed, NAyx_H ∪ PAy, y, debug)

            if newScore > bestScore
                bestH = H
                bestScore = newScore
            end
        end
    end

    return (bestH, bestScore)
end

####################################################################
# Scoring function
####################################################################

@memoize LRU(maxsize=10_000_000) function score(dataParsed::ParseData{Matrix{A}}, nodeParents, node, debug) where A
    #debug && message("Parents: $(nodeParents), Child: $(node)")
    #Unpack some variables from the dataParsed structure
    n = A(dataParsed.numObservations) #convert datatype
    scatterMat = dataParsed.scatterMat
    data = dataParsed.normAugData
    p = dataParsed.penalty

    #The last column of dataParsed.normAugData is all ones which is required for a linear regression with an intercept. If there are no node parents, model the child node with just the intercept, else use the parents and the intercept
    if isempty(nodeParents)
        parentsAndIncept = [dataParsed.numFeatures+1]
    else
        parentsAndIncept = [nodeParents; dataParsed.numFeatures+1]
    end

    #To calculate the score we need a mean-squared error which we can get by regessing the the child node onto the parents
    #Create variables for a linear regression y = X*b

    #Use views to avoid creating copies of the data in memory
        # X is the design matrix, augmented with a column of ones at the end
        # X is also been standardized so mean(columns)=0 and std(column)=1
        # y is data from the child node being tested
        # We're using the "normal equation" approach to calculate the linear regression
        # Instead of X\y we use (XᵀX)\(Xᵀy) which is slightly less accurate but in general twice as fast
        # XᵀX is precomputed in the scatter matrix so we need only to take slices it
    @views begin
        Xᵀy = scatterMat[parentsAndIncept,node]

        XᵀX = scatterMat[parentsAndIncept,parentsAndIncept]

        y = data[:,node]
        X = data[:,parentsAndIncept]
    end

    #Perform a linear regression
    b = XᵀX \ Xᵀy

    #I think the above is still using pivoted QR factorization, but to truly use cholesky decomposition I would have to use cholesky(Hermitian(XᵀX)). However, cholesky(⋅) doesn't have a method for SubArrays and not using a view is more costly.
    #Anyways, maybe this could be solved with LinearSolver.jl

    #Get the estimation
    ŷ = X*b

    #Next we want to calculate the log likelihood that these are the parents of our node
    # score = log(P(data|Model)) ≈ -BIC/2
    # because we're only comparing log likelihoods we'll ignore the 1/2 factor
    # when P(⋅) is Guassian, log(P(data|Model)) takes the form:
    # -k⋅log(n) - n⋅log(mse)
    # k is the number of free parameters and mse is mean squared error
    k = length(parentsAndIncept) #includes the intercept
    mse = sum(x->x^2, y-ŷ) / n

    #return the final score we want to maximize (which is technically -2BIC)
    return -p*k*log(n) - n*log(mse)
end

function statusUpdate(g::PDAG, phase)

    numNodes = g.nv
    maxEdges = numNodes*(numNodes - 1) ÷ 2

    aveDegree = round(degreeAverage(g), sigdigits=3)

    currentTime = Dates.format(now(), dateformat"m/d/yyyy II:MM p")

    message("------------------------")
    message("Updated at $(currentTime)")
    message("Phase: $(phase)")

    percentEdges = round(ne(g)/maxEdges, sigdigits=3)
    message("Edges $(ne(g)) / Total $(maxEdges) ($(percentEdges)%)")
    message("Average Degree: $(aveDegree)")
end


text:
MOI.canget(instance::ECOSOptimizer, ::MOI.TerminationStatus) = true
function MOI.get(instance::ECOSOptimizer, ::MOI.TerminationStatus)
    flag = instance.sol.ret_val
    if flag == ECOS.ECOS_OPTIMAL
        MOI.Success
    elseif flag == ECOS.ECOS_PINF
        MOI.Success
    elseif flag == ECOS.ECOS_DINF  # Dual infeasible = primal unbounded, probably
        MOI.Success
    elseif flag == ECOS.ECOS_MAXIT
        MOI.IterationLimit
    elseif flag == ECOS.ECOS_OPTIMAL + ECOS.ECOS_INACC_OFFSET
        m.solve_stat = MOI.AlmostSuccess
    else
        m.solve_stat = MOI.OtherError
    end
end

MOI.canget(instance::ECOSOptimizer, ::MOI.ObjectiveValue) = true
MOI.get(instance::ECOSOptimizer, ::MOI.ObjectiveValue) = instance.sol.objval

function MOI.canget(instance::ECOSOptimizer, ::MOI.PrimalStatus)
    instance.sol.ret_val != ECOS.ECOS_PINF
end
function MOI.get(instance::ECOSOptimizer, ::MOI.PrimalStatus)
    flag = instance.sol.ret_val
    if flag == ECOS.ECOS_OPTIMAL
        MOI.FeasiblePoint
    elseif flag == ECOS.ECOS_PINF
        MOI.InfeasiblePoint
    elseif flag == ECOS.ECOS_DINF  # Dual infeasible = primal unbounded, probably
        MOI.InfeasibilityCertificate
    elseif flag == ECOS.ECOS_MAXIT
        MOI.UnknownResultStatus
    elseif flag == ECOS.ECOS_OPTIMAL + ECOS.ECOS_INACC_OFFSET
        m.solve_stat = MOI.NearlyFeasiblePoint
    else
        m.solve_stat = MOI.OtherResultStatus
    end
end
# Swapping indices 2 <-> 3 is an involution (it is its own inverse)
const reorderval = orderval
function MOI.canget(instance::ECOSOptimizer, ::Union{MOI.VariablePrimal, MOI.ConstraintPrimal}, ::Type{<:MOI.Index})
    instance.sol.ret_val != ECOS.ECOS_PINF
end
function MOI.get(instance::ECOSOptimizer, ::MOI.VariablePrimal, vi::VI)
    instance.sol.primal[vi.value]
end
MOI.get(instance::ECOSOptimizer, a::MOI.VariablePrimal, vi::Vector{VI}) = MOI.get.(instance, a, vi)
setconstant(instance::ECOSOptimizer, offset, ::CI{<:MOI.AbstractFunction, <:MOI.EqualTo}) = instance.cone.eqsetconstant[offset]
setconstant(instance::ECOSOptimizer, offset, ::CI) = instance.cone.ineqsetconstant[offset]
_unshift(instance::ECOSOptimizer, offset, value, ::CI) = value
_unshift(instance::ECOSOptimizer, offset, value, ci::CI{<:MOI.AbstractScalarFunction, <:MOI.AbstractScalarSet}) = value + setconstant(instance, offset, ci)
function MOI.get(instance::ECOSOptimizer, ::MOI.ConstraintPrimal, ci::CI{<:MOI.AbstractFunction, MOI.Zeros})
    rows = constrrows(instance, ci)
    zeros(length(rows))
end
function MOI.get(instance::ECOSOptimizer, ::MOI.ConstraintPrimal, ci::CI{<:MOI.AbstractFunction, <:MOI.EqualTo})
    offset = constroffset(instance, ci)
    setconstant(instance, offset, ci)
end
function MOI.get(instance::ECOSOptimizer, ::MOI.ConstraintPrimal, ci::CI{<:MOI.AbstractFunction, S}) where S <: MOI.AbstractSet
    offset = constroffset(instance, ci)
    rows = constrrows(instance, ci)
    _unshift(instance, offset, scalecoef(rows, reorderval(instance.sol.slack[offset + rows], S), false, S), ci)
end

function MOI.canget(instance::ECOSOptimizer, ::MOI.DualStatus)
    instance.sol.ret_val != ECOS.ECOS_DINF
end
function MOI.get(instance::ECOSOptimizer, ::MOI.DualStatus)
    flag = instance.sol.ret_val
    if flag == ECOS.ECOS_OPTIMAL
        MOI.FeasiblePoint
    elseif flag == ECOS.ECOS_PINF
        MOI.InfeasibilityCertificate
    elseif flag == ECOS.ECOS_DINF  # Dual infeasible = primal unbounded, probably
        MOI.InfeasiblePoint
    elseif flag == ECOS.ECOS_MAXIT
        MOI.UnknownResultStatus
    elseif flag == ECOS.ECOS_OPTIMAL + ECOS.ECOS_INACC_OFFSET
        m.solve_stat = MOI.NearlyFeasiblePoint
    else
        m.solve_stat = MOI.OtherResultStatus
    end
end
function MOI.canget(instance::ECOSOptimizer, ::MOI.ConstraintDual, ::Type{<:CI})
    instance.sol.ret_val != ECOS.ECOS_DINF
end
_dual(instance, ci::CI{<:MOI.AbstractFunction, <:ZeroCones}) = instance.sol.dual_eq
_dual(instance, ci::CI) = instance.sol.dual_ineq
function MOI.get(instance::ECOSOptimizer, ::MOI.ConstraintDual, ci::CI{<:MOI.AbstractFunction, S}) where S <: MOI.AbstractSet
    offset = constroffset(instance, ci)
    rows = constrrows(instance, ci)
    scalecoef(rows, reorderval(_dual(instance, ci)[offset + rows], S), false, S)
end

MOI.canget(instance::ECOSOptimizer, ::MOI.ResultCount) = true
MOI.get(instance::ECOSOptimizer, ::MOI.ResultCount) = 1


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

using ValueShapes
using Test

using Random
using ElasticArrays
import TypedTables

@testset "abstract_value_shape" begin
    @testset "default_datatype" begin
        @test @inferred(ValueShapes.default_datatype(Integer)) == Int
        @test @inferred(ValueShapes.default_datatype(AbstractFloat)) == Float64
        @test @inferred(ValueShapes.default_datatype(Real)) == Float64
        @test @inferred(ValueShapes.default_datatype(Complex)) == Complex{Float64}
    end
end


text:
"""
    EncodedTableRow{M, N} <: Block

Block for processed rows having a tuple of M categorical and
N continuous value collections.
"""
struct EncodedTableRow{M, N, T} <: Block
    catcols::NTuple{M}
    contcols::NTuple{N}
    categorydict::T
end

function EncodedTableRow(catcols, contcols, categorydict)
    EncodedTableRow{length(catcols), length(contcols)}(catcols, contcols, categorydict)
end

function checkblock(::EncodedTableRow{M, N}, x::Tuple{Vector, Vector}) where {M, N}
    length(x[1]) == M && length(x[2]) == N
end

function showblock!(io, ::ShowText, block::EncodedTableRow, obs)
    print(io, "EncodedTableRow(...)")
end

# ## Encoding

"""
    TabularPreprocessing <: Encoding

Encodes a `TableRow` by applying the following preprocessing steps:
- [`DataAugmentation.NormalizeRow`](#) (for normalizing a row of data for continuous columns)
- [`DataAugmentation.FillMissing`](#) (for filling missing values)
- [`DataAugmentation.Categorify`](#) (for label encoding categorical columns,
    which can be later used for indexing into embedding matrices)
or a sequence of these transformations.
"""
struct TabularPreprocessing{T} <: Encoding
	tfms::T
end

TabularPreprocessing(td::Datasets.TableDataset) = TabularPreprocessing(gettransforms(td))

function encodedblock(::TabularPreprocessing, block::TableRow)
    EncodedTableRow(block.catcols, block.contcols, block.categorydict)
end

function encode(tt::TabularPreprocessing, _, block::TableRow, row)
    columns = Tables.columnnames(row)
    usedrow = NamedTuple(filter(
            x -> x[1] ∈ block.catcols || x[1] ∈ block.contcols,
            collect(zip(columns, row))
        ))
    tfmrow = DataAugmentation.apply(
        tt.tfms,
        DataAugmentation.TabularItem(usedrow, keys(usedrow))
    ).data
    catvals = collect(map(col -> tfmrow[col], block.catcols))
    contvals = collect(map(col -> tfmrow[col], block.contcols))
    (catvals, contvals)
end

function setup(::Type{TabularPreprocessing}, block::TableRow, data::TableDataset)
    return TabularPreprocessing(gettransforms(data, block.catcols, block.contcols))
end

# ## `blockmodel`

"""
    blockmodel(inblock::TableRow{M, N}, outblock::Union{Continuous, OneHotTensor{0}}, backbone=nothing) where {M, N}

Contruct a model for tabular classification or regression. `backbone` should be a
NamedTuple of categorical, continuous, and a finalclassifier layer, with
the first two taking in batches of corresponding row value matrices.
"""

"""
    blockmodel(::EncodedTableRow, ::OneHotTensor[, backbone])

Create a model for tabular classification. `backbone` should be named tuple
`(categorical = ..., continuous = ...)`. See [`TabularModel`](#) for more info.
"""
function blockmodel(inblock::EncodedTableRow, outblock::OneHotTensor{0}, backbone)
    TabularModel(
        backbone.categorical,
        backbone.continuous,
        Dense(100, length(outblock.classes))
    )
end

"""
    blockmodel(::EncodedTableRow, ::Continuous[, backbone])

Create a model for tabular regression. `backbone` should be named tuple
`(categorical = ..., continuous = ...)`. See [`TabularModel`](#) for more info.
"""
function blockmodel(inblock::EncodedTableRow, outblock::Continuous, backbone)
    TabularModel(
        backbone.categorical,
        backbone.continuous,
        Dense(100, outblock.size)
    )
end

function blockbackbone(inblock::EncodedTableRow{M, N}) where {M, N}
    embedszs = _get_emb_sz(collect(map(col->length(inblock.categorydict[col]), inblock.catcols)))
    catback = tabular_embedding_backbone(embedszs)
    contback = tabular_continuous_backbone(N)
    return (categorical = catback, continuous = contback)
end

# ## Utilities

"""
The helper functions defined below can be used for quickly constructing a dictionary,
which will be required for creating various tabular transformations available in DataAugmentation.jl.

These functions assume that the table in the TableDataset object td has Tables.jl columnaccess interface defined.
"""
function gettransformdict(td, ::Type{DataAugmentation.NormalizeRow}, cols)
    dict = Dict()
    map(cols) do col
        vals = skipmissing(Tables.getcolumn(td.table, col))
        dict[col] = (Statistics.mean(vals), Statistics.std(vals))
    end
    dict
end

function gettransformdict(td, ::Type{DataAugmentation.FillMissing}, cols)
    dict = Dict()
    map(cols) do col
        vals = skipmissing(Tables.getcolumn(td.table, col))
        dict[col] = Statistics.median(vals)
    end
    dict
end

function gettransformdict(td, ::Type{DataAugmentation.Categorify}, cols)
    dict = Dict()
    map(cols) do col
        vals = Tables.getcolumn(td.table, col)
        dict[col] = unique(vals)
    end
    dict
end

"""
    getcoltypes(td::Datasets.TableDataset)

Returns the categorical and continuous columns present in a `TableDataset`.
"""
function getcoltypes(td::Datasets.TableDataset)
    schema = Tables.schema(td.table)

    contcols = Tuple(name for (name, T) in zip(schema.names, schema.types)
        if T <: Union{<:Number, <:Union{Missing, <:Number}})

    catcols = Tuple(name for name in schema.names if !(name in contcols))
    catcols, contcols
end

"""
    gettransforms(td::Datasets.TableDataset)

Returns a composition of basic tabular transformations constructed
for the given TableDataset.
"""
function gettransforms(td::TableDataset, catcols, contcols)
    normstats = gettransformdict(td, DataAugmentation.NormalizeRow, contcols)
    fmvals = gettransformdict(td, DataAugmentation.FillMissing, contcols)
    catdict = gettransformdict(td, DataAugmentation.Categorify, catcols)
    normalize = DataAugmentation.NormalizeRow(normstats, contcols)
    categorify = DataAugmentation.Categorify(catdict, catcols)
    fm = DataAugmentation.FillMissing(fmvals, contcols)

    return fm |> normalize |> categorify
end

gettransforms(td::TableDataset) = gettransforms(td, getcoltypes(td)...)

# ## Tests

@testset "TabularPreprocessing [encoding]" begin
    cols = [:col1, :col2, :col3, :col4, :col5]
    vals = [1, 2, 3, "a", "x"]
    row = NamedTuple(zip(cols, vals))

    catcols = (:col4, :col5)
    contcols = (:col1, :col2, :col3)

    col1_mean, col1_std = 10, 100
    col2_mean, col2_std = 100, 10
    col3_mean, col3_std = 15, 1

    normdict = Dict(
        :col1 => (col1_mean, col1_std),
        :col2 => (col2_mean, col2_std),
        :col3 => (col3_mean, col3_std)
    )

    tfm = TabularPreprocessing(
        DataAugmentation.NormalizeRow(normdict, contcols)
    )

    block = TableRow(
        catcols,
        contcols,
        Dict(:col4=>["a", "b"], :col5=>["x", "y", "z"])
    )

    testencoding(tfm, block, row)
    testencoding(setup(TabularPreprocessing, block, TableDataset(DataFrame([row, row]))), block, row)
end

@testset "blockbackbone" begin
    @test_nowarn FastAI.blockbackbone(EncodedTableRow((:x,), (:y,), Dict(:x => [1, 2])))
end


text:
# Command-Line 2048 Example
module CLI2048

using Termbox
using Compat
if VERSION >= v"0.7-"
    using Random
    using Printf
end

const Grid = Matrix

function Base.show(io::IO, g::Grid)
    n = size(g, 1)
    for row in 1:n
        for col in 1:(n-1)
            print(io, g[row, col], '\t')
        end
        print(io, g[row, n], '\n')
    end
end

newGrid(n::Int = 4, T = Int) = zeros(T, n, n)

function randomGrid(n::Int = 4, numNZ::Integer = 2, NonzeroElements::Vector{T} = [2,4]) where {T}
    randomGrid(newGrid(n, T), numNZ, NonzeroElements)
end
function randomGrid(grid::Grid{T}, numNZ::Integer = 2, NonzeroElements::Vector{T} = [2,4]) where {T}
    nGrid = length(grid)
    for ind in randperm(nGrid)[1:numNZ]
        grid[ind] = rand(NonzeroElements)
    end
    grid
end
zeroElements(g::Grid) = findall([gv == 0 for gv in g])
function randomGrid!(grid::Grid{T}, numNZ::Integer = 2, NonzeroElements::Vector{T} = [2,4]) where {T}
    inds = shuffle(zeroElements(grid))[1:numNZ]
    for ind in inds
        grid[ind] = rand(NonzeroElements)
    end
end

#TODO: Move asserts into some kind of test function
v1 = newGrid(4, Int)
@assert length(v1) == 4*4
@assert sum(v1) == 0
v1Rand = randomGrid(v1)
@assert size(v1Rand) == size(v1)
@assert sum(v1Rand) >= 2*2
v2 = randomGrid()
@assert sum(v2) >= 2*2
v3 = randomGrid(4, 3, [2,])
@assert sum(v3) == 3*2

# Functions to complete a single move in one direction
# NOTE: Output from below can be generated ahead of time and cached.
function compactLeft!(v::AbstractVector{T}, n) where {T}
    row = zeros(T, n)
    p = 1
    for (i, vi) in enumerate(v)
        if vi != 0
            row[p] = vi
            p += 1
        end
    end
    for i in eachindex(v)
        v[i] = row[i]
    end
end
function mergeLeft!(v::AbstractVector{T}, diffRef) where {T}
    n = length(v)
    compactLeft!(v, n)
    for i in 1:(n-1)
        vi = v[i]
        vip = v[i + 1]
        if vip == vi
            diffRef[] += vi
            v[i] *= 2
            v[i+1] = 0
        end
    end
    compactLeft!(v, n)
    v
end

function move!(g::Grid, dir::Val{:left}, scoreDiffRef)
    n = size(g, 1)
    for row in 1:n
        mergeLeft!(view(g, row, :), scoreDiffRef)
    end
end
function move!(g::Grid, dir::Val{:right}, scoreDiffRef)
    n = size(g, 1)
    for row in 1:n
        reverse!(mergeLeft!(reverse!(view(g, row, :)), scoreDiffRef))
    end
end
function move!(g::Grid, dir::Val{:up}, scoreDiffRef)
    n = size(g, 1)
    for col in 1:n
        mergeLeft!(view(g, :, col), scoreDiffRef)
    end
end
function move!(g::Grid, dir::Val{:down}, scoreDiffRef)
    n = size(g, 1)
    for col in 1:n
        reverse!(mergeLeft!(reverse!(view(g, :, col)), scoreDiffRef))
    end
end

# Termbox GUI
function putchar(k, i, c)
    tb_change_cell(k, i, c, 0x0001, 0x0000)
end

#TODO: Refactor interface below into boxes
function tb_print(g::Grid, numWidth::Int = 4)
    tb_clear()
    n = size(g, 1)
    rowCounter = 1
    for row in 1:n
        for i in 1:(n*(numWidth+4)+1)
            putchar(i, rowCounter, '#')
        end
        rowCounter += 1
        colCounter = 1
        for col in 1:n
            putchar(colCounter, rowCounter, '#')
            colCounter += 2 + numWidth + 2
            putchar(colCounter, rowCounter, '#')
        end
        rowCounter += 1
        colCounter = 1
        for col in 1:n
            colCounter += 2
            for (i, c) in enumerate(rpad(sprint(show, g[row, col]), numWidth))
                putchar(colCounter, rowCounter, c)
                colCounter += 1
            end
            colCounter += 2
            putchar(colCounter, rowCounter, '#')
        end
        rowCounter += 1
        colCounter = 1
        for col in 1:n
            putchar(colCounter, rowCounter, '#')
            colCounter += 2 + numWidth + 2
            putchar(colCounter, rowCounter, '#')
        end
        rowCounter += 1
    end
    for i in 1:(n*(numWidth+4)+1)
        putchar(i, rowCounter, '#')
    end
end

function handle_char!(c, g)
    if c == Int('r')
        fill!(g, 0)
        randomGrid!(g)
        tb_clear()
        return true
    end
    if c == Int('q')
        return false
    end
    return true
end
function handle_key!(k, g)
    scorediff = Ref(0)
    if k == TB_KEY_ARROW_LEFT
        move!(g, Val(:left), scorediff)
    elseif k == TB_KEY_ARROW_RIGHT
        move!(g, Val(:right), scorediff)
    elseif k == TB_KEY_ARROW_UP
        move!(g, Val(:up), scorediff)
    elseif k == TB_KEY_ARROW_DOWN
        move!(g, Val(:down), scorediff)
    end
    if isempty(zeroElements(g))
        return -1
    end
    if scorediff[] >= 0
        randomGrid!(g, 1)
    end
    return scorediff[]
end
function main_interactive_loop(boardSize::Int = 4)
    do_loop = true
    score = 0
    # g = newGrid(4)
    # g[2, 1] = 2
    # g[2, 2] = 2

    g = randomGrid(boardSize)

    tb_print(g)
    tb_present()
    while do_loop
        ev = tb_event()
        while tb_poll_event(ev) != 0
            if ev.etype == TB_EVENT_KEY
                if ev.key == TB_KEY_ESC
                    do_loop = false
                    break
                elseif ev.key == 0
                    if !(handle_char!(ev.ch, g))
                        do_loop = false
                        break
                    end
                else
                    scorediff = handle_key!(ev.key, g)
                    if scorediff < 0
                        do_loop = false
                        break
                    end
                    score += scorediff
                end
            end #if
            tb_print(g)
            tb_present()
            tb_peek_event(ev, 0)
        end
    end
    score
end

function main_interactive(boardSize::Int = 4)
    ret = tb_init()
    if ret != 0
        @printf stderr "tb_init() failed with error code %d\n" ret
        return 1
    end
    score = 0
    try
        tb_clear()
        tb_select_input_mode(TB_INPUT_ESC)
        tb_select_output_mode(TB_OUTPUT_NORMAL)

        score += main_interactive_loop(boardSize)
    finally
        tb_shutdown()
    end
    @show score
    return 0
end

function main_wargames_loop(sleepTime::Real = 0.05, boardSize::Int = 4, maxMoves::Integer = 1000)
    score = 0
    moves = 0

    g = randomGrid(boardSize)

    tb_print(g)
    tb_present()
    while (moves < maxMoves)
        scorediff = handle_key!(rand([TB_KEY_ARROW_LEFT, TB_KEY_ARROW_RIGHT, TB_KEY_ARROW_DOWN, TB_KEY_ARROW_UP]), g)
        if scorediff < 0
            break
        end
        score += scorediff
        tb_print(g)
        tb_present()
        sleep(sleepTime)
        moves += 1
    end
    score
end

function main_wargames(sleepTime::Real = 0.05, boardSize::Int = 4)
    ret = tb_init()
    if ret != 0
        @printf stderr "tb_init() failed with error code %d\n" ret
        return 1
    end
    score = 0
    try
        tb_clear()
        tb_select_input_mode(TB_INPUT_ESC)
        tb_select_output_mode(TB_OUTPUT_NORMAL)

        score += main_wargames_loop(sleepTime, boardSize)
    finally
        tb_shutdown()
    end
    @show score
    return 0
end
# v2 = newGrid(4)
# v2[2, 1] = 2
# v2[2, 2] = 2
# println(v2)
# println(move!(v2, Val(:left)))
# println(v2)

end # module

# To use this interactively, run
# julia -L 2048Example.jl -e 'CLI2048.main_interactive()'
#
# To see a game played randomly (and quickly), run
# julia -L 2048Example.jl -e 'CLI2048.main_wargames()'
# The latter command is suitable for iteration on the code itself.


text:
@testset "runtestWarningCountdown" begin
    println("Testing runtestWarningCountdown")
    include("../src/bfgaWarningCountdown.jl")
    using .bfgaWarningCountdown

    include("../src/Bf.jl")
    using .BfInterpreter

    @time model = test_serial()
    #@time model = test_parallel()
    prems = model.population[length(model.population)]
    println("Generation : $(model.params.currentGeneration)")
    println("& fitness : $(model.specific_fitness.fitness(prems, model.instructionsSet))")
    #bfcode :: Array{Float64,1}  = prems.dna
    println("code :  $(prems.program)") #$(join(genesToBfInstructions(bfcode), "") ) ")

    for elt in _trainingExamples
        output, _ = execute(prems.program, elt , model.instructionsSet)
        @show join(output, "")
        @test length(output) >= l_length * (target_length +1)
    end

    # do the test_parallel after ...
    #test_parallel()
end


text:
const GLOBAL_ATOL = Ref(1e-8)

########### macro tools #############
startwithdot(sym::Symbol) = string(sym)[1] == '.'
startwithdot(sym::Expr) = false
startwithdot(sym) = false
function removedot(f)
    string(f)[1] == '.' || error("$f is not a broadcasting.")
    Symbol(string(f)[2:end])
end

"""
    get_argname(ex)

Return the argument name of a function argument expression, e.g. `x::Float64 = 4` gives `x`.
"""
function get_argname(fname)
    @smatch fname begin
        ::Symbol => fname
        :($x::$t) => x
        :($x::$t=$y) => x
        :($x=$y) => x
        :($x...) => :($x...)
        :($x::$t...) => :($x...)
        Expr(:parameters, args...) => fname
        _ => error("can not get the function name of expression $fname.")
    end
end

"""
    match_function(ex)

Analyze a function expression, returns a tuple of `(macros, function name, arguments, type parameters (in where {...}), statements in the body)`
"""
function match_function(ex)
    @smatch ex begin
        :(function $(fname)($(args...)) $(body...) end) ||
        :($fname($(args...)) = $(body...)) => (nothing, fname, args, [], body)
        Expr(:function, :($fname($(args...)) where {$(ts...)}), xbody) => (nothing, fname, args, ts, xbody.args)
        Expr(:macrocall, mcname, line, fdef) => ([mcname, line], match_function(fdef)[2:end]...)
        _ => error("must input a function, got $ex")
    end
end

"""
    rmlines(ex::Expr)

Remove line number nodes for pretty printing.
"""
rmlines(ex::Expr) = begin
    hd = ex.head
    if hd == :macrocall
        Expr(:macrocall, ex.args[1], nothing, rmlines.(ex.args[3:end])...)
    else
        tl = Any[rmlines(ex) for ex in ex.args if !(ex isa LineNumberNode)]
        Expr(hd, tl...)
    end
end
rmlines(@nospecialize(a)) = a

########### ordered dict ###############
struct MyOrderedDict{TK,TV}
    keys::Vector{TK}
    vals::Vector{TV}
end

function MyOrderedDict{K,V}() where {K,V}
    MyOrderedDict(K[], V[])
end

function Base.setindex!(d::MyOrderedDict, val, key)
    ind = findfirst(x->x===key, d.keys)
    if ind isa Nothing
        push!(d.keys, key)
        push!(d.vals, val)
    else
        @inbounds d.vals[ind] = val
    end
    return d
end

function Base.getindex(d::MyOrderedDict, key)
    ind = findfirst(x->x===key, d.keys)
    if ind isa Nothing
        throw(KeyError(ind))
    else
        return d.vals[ind]
    end
end

function Base.delete!(d::MyOrderedDict, key)
    ind = findfirst(x->x==key, d.keys)
    if ind isa Nothing
        throw(KeyError(ind))
    else
        deleteat!(d.vals, ind)
        deleteat!(d.keys, ind)
    end
end

Base.length(d::MyOrderedDict) = length(d.keys)

function Base.pop!(d::MyOrderedDict)
    k = pop!(d.keys)
    v = pop!(d.vals)
    k, v
end

Base.isempty(d::MyOrderedDict) = length(d.keys) == 0

########### broadcasting ###############

"""
    unzipped_broadcast(f, args...)

unzipped broadcast for arrays and tuples, e.g. `SWAP.([1,2,3], [4,5,6])` will do inplace element-wise swap, and return `[4,5,6], [1,2,3]`.
"""
unzipped_broadcast(f) = error("must provide at least one argument in broadcasting!")
function unzipped_broadcast(f, arg::AbstractArray; kwargs...)
    arg .= f.(arg)
end
function unzipped_broadcast(f, arg::Tuple; kwargs...)
    f.(arg)
end
@generated function unzipped_broadcast(f, args::Vararg{<:AbstractArray,N}; kwargs...) where N
    argi = [:(args[$k][i]) for k=1:N]
    quote
        for i = 1:same_length(args)
            ($(argi...),) = f($(argi...); kwargs...)
        end
        return args
    end
end
@generated function unzipped_broadcast(f, args::Vararg{<:Tuple,N}; kwargs...) where N
    quote
        same_length(args)
        res = map(f, args...)
        ($([:($getindex.(res, $i)) for i=1:N]...),)
    end
end

function same_length(args)
    length(args) == 0 && error("can not broadcast over an empty set of arguments.")
    l = length(args[1])
    for j=2:length(args)
        @assert l == length(args[j]) "length of arguments should be the same `$(length(args[j])) != $l`"
    end
    return l
end


text:
# !INPUT/OUTPUT PARAMETERS:
#   z : atomic number (in,real)
# !DESCRIPTION:
#   Computes an approximate nuclear charge radius from the atomic number $Z$.
#   The nuclear mass number, $A$, is estimated using
#   $$ A=4.467\times 10^{-3}Z^2+2.163 Z-1.168, $$
#   [D. Andrae in {\it Relativistic Electronic Structure Theory - Fundamentals}
#   {\bf 11}, 203 (2002)], and the nuclear charge radius can be determined from
#   $$ r=\left(r_0+\frac{r_1}{A^{2/3}}+\frac{r_2}{A^{4/3}}\right)A^{1/3}, $$
#   where $r_0=0.9071$, $r_1=1.105$ and $r_2=-0.548$ [I. Angeli, {\it Atomic
#   Data and Nuclear Data Tables} {\bf 87}, 185 (2004)].

function radnucl( z::Float64 )

    # coefficients for computing mass number
    c2 = 4.467e-3
    c1 = 2.163
    c0 = -1.168

    # coefficients for computing charge radius (fm)
    r0 = 0.9071
    r1 = 1.105
    r2 = -0.548
    # Bohr radius in SI units (CODATA 2018)
    br_si = 0.529177210903e-10

    za = abs(z)
    a = 1.0
    # approximate nuclear mass number
    if za <= 1.0
        a = 1.0
    else
        a = c2*za^2 + c1*za + c0
    end
    # approximate nuclear charge radius
    a13 = a^(1.0/3.0)
    a23 = a13^2
    a43 = a13*a
    return (r0 + r1/a23 + r2/a43)*a13 * 1.0e-15/br_si

end


text:
using NumericalMethodsforEngineers

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  f1(x::Float64, y::Vector{Float64}) = [y[2], 2.0*y[1]-3.0*y[2]+3.0*x^2]
  steps = 4
  h = 0.05

  x = 0.0
  y = [1.0, 0.0]

  r = Array{Float64,2}[]
  push!(r, euler(f1, x, y, steps, h))
  push!(r, modified_euler(f1, x, y, steps, h))
  push!(r, mid_point_euler(f1, x, y, steps, h))
  push!(r, runga_kutta_4(f1, x, y, steps, h))

  @test r[1][steps+1,:][:] == [0.2, 1.0272471874999998, 0.3254353125]
  @test r[2][steps+1,:][:] == [0.2, 1.0338953695013426, 0.3110140747716065]
  @test r[3][steps+1,:][:] == [0.2, 1.033871576010437, 0.3107074949932251]
  @test r[4][steps+1,:][:] == [0.2, 1.033615605146455, 0.3112553545583527]

end


text:
using DynamicHMCModels, MCMCChains
using Flux

ProjDir = @__DIR__
cd(ProjDir)

include(joinpath(@__DIR__, "LBA_functions.jl"))

Base.@kwdef struct LBAModel{T}
   data::T
   N::Int
   Nc::Int
end

function make_transformation(model::LBAModel)
    as((v=as(Array,asℝ₊,Nc),A=asℝ₊,k=asℝ₊,tau=asℝ₊))
end

N = 20
data=simulateLBA(;Nd=N,v=[1.0,1.5,2.0],A=.8,k=.2,tau=.4)
Nc = 3

#dist = LBA(ν=[1.0,1.5,2.0],A=.8,k=.2,τ=.4)
#data = rand(dist,N)

model = LBAModel(; data=data, N=N, Nc=Nc)

function (model::LBAModel)(θ)
    @unpack data=model
    @unpack v,A,k,tau=θ
    d=LBA(ν=v,A=A,k=k,τ=tau)
    minRT = minimum(x->x[2],data)
    logpdf(d,data)+sum(logpdf.(TruncatedNormal(0,3,0,Inf),v)) +
    logpdf(TruncatedNormal(.8,.4,0,Inf),A)+logpdf(TruncatedNormal(.2,.3,0,Inf),k)+
    logpdf(TruncatedNormal(.4,.1,0,minRT),tau)
end

d = [(c,r) for (c,r) in zip(data.choice,data.rt)]

p = LBAModel(d,N,Nc)
p((v=fill(.5,Nc),A=.8,k=.2,tau=.4))

# Use Flux for the gradient.
P = TransformedLogDensity(make_transformation(p), p)
∇P = ADgradient(:Flux, P)

# Sample from the posterior.
results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 1000)
posterior = P.transformation.(results.chain)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule nPth_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("nPth")
JLLWrappers.@generate_main_file("nPth", UUID("58695248-dc95-54d8-8aad-a18bc9fe19aa"))
end  # module nPth_jll


text:
# https://adventofcode.com/2020/day/7

example1 = "
light red bags contain 1 bright white bag, 2 muted yellow bags.
dark orange bags contain 3 bright white bags, 4 muted yellow bags.
bright white bags contain 1 shiny gold bag.
muted yellow bags contain 2 shiny gold bags, 9 faded blue bags.
shiny gold bags contain 1 dark olive bag, 2 vibrant plum bags.
dark olive bags contain 3 faded blue bags, 4 dotted black bags.
vibrant plum bags contain 5 faded blue bags, 6 dotted black bags.
faded blue bags contain no other bags.
dotted black bags contain no other bags.
"

struct Bag
    kind::AbstractString
    quantity::Integer

    Bag(kind, quantity) = new(kind, quantity)

    function Bag(string::AbstractString)
        m = match(r"(\d*?) ?([a-z ]+?) (?:bag|bags)", string)
        if isnothing(m)
            error("invalid bag description")
        end
        quantity, kind = m.captures

        if kind == "no other"
            return nothing
        end
        return new(String(kind), quantity == "" ? 1 : parse(Int, quantity))
    end
end

function rule_parser(rule)
    bag, contains = split(rule, " contain ", limit=2)
    bag = Bag(bag)
    contains = [Bag(elem) for elem in split(contains, ',')]
    return bag, filter(x -> !isnothing(x), contains)
end

function rules_to_revdict(input)
    bags = Dict()

    for row in split(input, '\n', keepempty=false)
        bag, contains = rule_parser(row)

        for elem in contains
            if isnothing(elem)
                continue
            end

            if elem.kind in keys(bags)
                push!(bags[elem.kind], bag.kind)
            else
                bags[elem.kind] = [bag.kind]
            end
        end
    end

    return bags
end

@assert rules_to_revdict(example1)["shiny gold"] == ["bright white", "muted yellow"]

function search_for_possible_bag_colors(bags, target="shiny gold")
    possible_colors = Set()
    to_search = [target]

    while !isempty(to_search)
        found = []
        for elem in to_search
            try
                # we already have those on the list, ignore in case there are cycles
                for x in bags[elem]
                    push!(found, x)
                    push!(possible_colors, x)
                end
            catch KeyError
            end

        end
        to_search = found
    end
    return possible_colors
end

function part1(input)
    bags = rules_to_revdict(input)
    possible_colors = search_for_possible_bag_colors(bags)
    return length(possible_colors)
end

@assert part1(example1) == 4

function rules_to_dict(input)
    rules = Dict{String,Array{Bag}}()
    for row in split(input, '\n', keepempty=false)
        bag, contains = rule_parser(row)
        rules[bag.kind] = contains
    end
    return rules
end

@assert rules_to_dict(example1)["shiny gold"] == Bag[Bag("dark olive", 1), Bag("vibrant plum", 2)]

function count_bags(rules, bag)
    if isempty(rules[bag.kind])
        return bag.quantity
    else
        n = 0
        for elem in rules[bag.kind]
            n += count_bags(rules, elem)
        end
        return bag.quantity + bag.quantity * n
    end
end

@assert count_bags(rules_to_dict(example1), Bag("faded blue", 1)) == 1
@assert count_bags(rules_to_dict(example1), Bag("faded blue", 3)) == 3
@assert count_bags(rules_to_dict(example1), Bag("vibrant plum", 1)) == 12
@assert count_bags(rules_to_dict(example1), Bag("vibrant plum", 2)) == 24

function part2(input)
    rules = rules_to_dict(input)
    return count_bags(rules, Bag("shiny gold", 1)) - 1
end

example2 = "
shiny gold bags contain 2 dark red bags.
dark red bags contain 2 dark orange bags.
dark orange bags contain 2 dark yellow bags.
dark yellow bags contain 2 dark green bags.
dark green bags contain 2 dark blue bags.
dark blue bags contain 2 dark violet bags.
dark violet bags contain no other bags.
"

@assert part2(example2) == 126

test = read("data/day-07.txt", String)
println("Part 1: $(result1 = part1(test))")
println("Part 2: $(result2 = part2(test))")

@assert result1 == 222
@assert result2 == 13264

max_stars_repo_name: twolodzko/advent-of-code

max_issues_repo_name: twolodzko/advent-of-code

max_forks_repo_name: twolodzko/advent-of-code


text:
# Functions related to telinfo.yml files

using YAML: load_file

"""
    telinfo(yamlfile)

Loads YAML file `yamlfile` using `dicttype=Dict{Symbol,Amy}`.
"""
telinfo(fname) = load_file(fname, dicttype=Dict{Symbol,Any})

"""
    antpos(telinfo_yml::AbstracyString)
    antpos(telinfo::AbstractDict{Symbol, Any})
    antpos(antennas::AbstractVector{<:AbstractDict{Symbol,Amy}})

Return antenna positions from `antennas` Vector as a 3xN matrix (i.e. one
antenna per columm).  The antenna ordering is sorted by the `number` of each
antenna.

# Examples

```julia
ti = telinfo("telinfo.yml") # => Dict{Symbol, Any}

ap = antpos("telinfo.yml")  # => Array{Float64,2}
ap = antpos(ti)             # => Array{Float64,2}
ap = antpos(ti[:antennas])  # => Array{Float64,2}
```

!!! note
    The antenna numbers are used for sorting, but they are NOT antenna indexes.
"""
function antpos(antennas::AbstractVector{<:AbstractDict{Symbol,Any}})
    ants = sort(antennas, by=a->a[:number])
    reduce(hcat, (a[:position] for a in ants))
end,

function antpos(ti::AbstractDict{Symbol,Any})
    antpos(ti[:antennas])
end,

function antpos(tiyml::AbstractString)
    antpos(telinfo(tiyml))
end


text:
using StaticArrays
using Parameters
using LaserTypes
using LinearAlgebra
using Statistics
using ThreadsX
using Plots
using LaTeXStrings
using RecursiveArrayTools
using QuadGK
using Serialization

function folder_setup(path)
    path = replace(path,"/" => "")

    if split(pwd(),"/")[end] == path
        println("folder ok")
        return
    end

    if !isdir(path)
        mkdir(path)
        println("making folder")
    end

    cd(path)
    println("moved in $path")
end

cd(@__DIR__)
folder_setup("output5")

const c = 137.036 # speed of light
const ω = 0.057; const T₀ = 2π/ω # angular frequency, period
const k = ω/c # wavenumber
const λ = c*T₀ # wavelength
const w₀ = 5 * λ # is the waist radius
const a₀ = 2.; const AA = a₀*c # (qA₀)/(mc) = a₀
const w = 10.; const τ = w/ω; # temporal pulse decay time given in terms of number of oscilations and frequency
const n = 5 # number of periods to integrate before and after pulse collides with particles
const π₀ = a₀*c # in atomic units, a₀mc has units of linear momentum and sets the scale for linear momentum transfered to the particle, not that m = 1 for our particles
const maxR = 3.25w₀ # maximum radius for distributing praticles

@with_kw struct InterPars{qq,mm,typo,rewts}
    q::qq = -1
    m::mm = 1
    type::typo
    roots::rewts = Float64[]
end

################################################################################
#Right Circularly Polarized

# p = 0
CPLGTYPE00 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = 0, p = 0,  ξx = (1. + 0im)/√2, ξy = (1im)/√2, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))
CPLGTYPE0_1 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = -1, p = 0,  ξx = (1. + 0im)/√2, ξy = (1im)/√2, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))
CPLGTYPE01 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = 1, p = 0,  ξx = (1. + 0im)/√2, ξy = (1im)/√2, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))

# p = 1
CPLGTYPE10 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = 0, p = 1,  ξx = (1. + 0im)/√2, ξy = (1im)/√2, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))
CPLGTYPE1_1 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = -1, p = 1,  ξx = (1. + 0im)/√2, ξy = (1im)/√2, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))
CPLGTYPE11 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = 1, p = 1,  ξx = (1. + 0im)/√2, ξy = (1im)/√2, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))

# p = 2
CPLGTYPE2_1 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = -1, p = 2,  ξx = (1. + 0im)/√2, ξy = (1im)/√2, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))
CPLGTYPE21 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = 1, p = 2,  ξx = (1. + 0im)/√2, ξy = (1im)/√2, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))
CPLGTYPE20 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = 0, p = 2,  ξx = (1. + 0im)/√2, ξy = (1im)/√2, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))
CPLGTYPE2_2 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = -2, p = 2,  ξx = (1. + 0im)/√2, ξy = (1im)/√2, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))
CPLGTYPE22 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = 2, p = 2,  ξx = (1. + 0im)/√2, ξy = (1im)/√2, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))

################################################################################
#Left Circularly Polarized

# p = 0
MCPLGTYPE00 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = 0, p = 0,  ξx = (1. + 0im)/√2, ξy = (-1im)/√2, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))
MCPLGTYPE0_1 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = -1, p = 0,  ξx = (1. + 0im)/√2, ξy = (-1im)/√2, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))
MCPLGTYPE01 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = 1, p = 0,  ξx = (1. + 0im)/√2, ξy = (-1im)/√2, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))

# p = 1
MCPLGTYPE10 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = 0, p = 1,  ξx = (1. + 0im)/√2, ξy = (-1im)/√2, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))
MCPLGTYPE1_1 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = -1, p = 1,  ξx = (1. + 0im)/√2, ξy = (-1im)/√2, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))
MCPLGTYPE11 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = 1, p = 1,  ξx = (1. + 0im)/√2, ξy = (-1im)/√2, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))

################################################################################
#Linearly Polarized

# p = 0
LPLGTYPE00 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = 0, p = 0, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))
LPLGTYPE0_1 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = -1, p = 0, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))
LPLGTYPE01 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = 1, p = 0, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))

# p = 1
LPLGTYPE10 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = 0, p = 1, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))
LPLGTYPE1_1 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = -1, p = 1, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))
LPLGTYPE11 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = 1, p = 1, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))

# p = 2
LPLGTYPE20 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = 0, p = 2, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))
LPLGTYPE21 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = 1, p = 2, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))
LPLGTYPE22 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = 2, p = 2, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))
LPLGTYPE2_1 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = -1, p = 2, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))
LPLGTYPE2_2 = LaguerreGaussLaser(c = c, m_q = 1., q = -1., λ = λ, a₀ = a₀, ϕ₀ = -π/2,
    w₀ = w₀, m = -2, p = 2, profile = ConstantProfile()) # GaussProfile(c = c, τ = w/ω, z₀ = 0.))

field_param = (
    CPLGTYPE0_1 = InterPars(type = CPLGTYPE0_1, roots = [0., 1.954]),
    CPLGTYPE00 = InterPars(type = CPLGTYPE00, roots = [0., 1.518]),
    CPLGTYPE01 = InterPars(type = CPLGTYPE01, roots = [0., 1.954]),
    CPLGTYPE1_1 = InterPars(type = CPLGTYPE1_1, roots = [0., 1., 2.566]),
    CPLGTYPE10 = InterPars(type = CPLGTYPE10, roots = [ 0., 0.707, 2.323]),
    CPLGTYPE11 = InterPars(type = CPLGTYPE11, roots = [0., 1., 2.566]),
    CPLGTYPE2_2 = InterPars(type = CPLGTYPE2_2, roots = [0., 1.000, 1.731, 3.176]),
    CPLGTYPE2_1 = InterPars(type = CPLGTYPE2_1, roots = [0., 0.797, 1.538, 3.005]),
    CPLGTYPE20 = InterPars(type = CPLGTYPE20, roots = [0., 0.541, 1.306, 2.807]),
    CPLGTYPE21 = InterPars(type = CPLGTYPE21, roots = [0., 0.797, 1.538, 3.005]),
    CPLGTYPE22 = InterPars(type = CPLGTYPE22, roots = [0., 1.000, 1.731, 3.176]),
    MCPLGTYPE0_1 = InterPars(type = MCPLGTYPE0_1, roots = [0., 1.954]),
    MCPLGTYPE00 = InterPars(type = MCPLGTYPE00, roots = [0., 1.518]),
    MCPLGTYPE01 = InterPars(type = MCPLGTYPE01, roots = [0., 1.954]),
    MCPLGTYPE1_1 = InterPars(type = MCPLGTYPE1_1, roots = [0., 1., 2.566]),
    MCPLGTYPE10 = InterPars(type = MCPLGTYPE10, roots = [0., 0.707, 2.323]),
    MCPLGTYPE11 = InterPars(type = MCPLGTYPE11, roots = [0., 1., 2.566]),
    LPLGTYPE0_1 = InterPars(type = LPLGTYPE0_1, roots = [0., 1.954]),
    LPLGTYPE00 = InterPars(type = LPLGTYPE00, roots = [0., 1.518]),
    LPLGTYPE01 = InterPars(type = LPLGTYPE01, roots = [0., 1.954]),
    LPLGTYPE1_1 = InterPars(type = LPLGTYPE1_1, roots = [0., 1., 2.566]),
    LPLGTYPE10 = InterPars(type = LPLGTYPE10, roots = [0., 0.707, 2.323]),
    LPLGTYPE11 = InterPars(type = LPLGTYPE11, roots = [0., 1., 2.566]),
    LPLGTYPE2_2 = InterPars(type = LPLGTYPE2_2, roots = [0., 1.000, 1.731, 3.176]),
    LPLGTYPE2_1 = InterPars(type = LPLGTYPE2_1, roots = [0., 0.797, 1.538, 3.005]),
    LPLGTYPE20 = InterPars(type = LPLGTYPE20, roots = [0., 0.541, 1.306, 2.807]),
    LPLGTYPE21 = InterPars(type = LPLGTYPE21, roots = [0., 0.797, 1.538, 3.005]),
    LPLGTYPE22 = InterPars(type = LPLGTYPE22, roots = [0., 1.000, 1.731, 3.176])
)

const τi = - n * τ
const τf =   n * τ
const time_samples = 2000
const Δτ = range(0, stop = 2*n*τ/T₀, length = time_samples+1)

const domainXY = range(-maxR, stop = maxR, length = 800)/w₀
const domainR  = range(0., stop = maxR, length = 800)/w₀
const time = range(- n * τ, stop = n * τ, length = 1200)/T₀

rectangle(w, h, x, y) = Shape(x .+ [0,w,w,0], y .+ [0,0,h,h])

function add_info!(plt, p, isscatter = false)
    #there is either a bug or confusing behaviour with ylims and aspect_ratio
    xl, yl = (isscatter ? ((-5.658368025329326, 5.667117711980884), (-3.662004592945377, 3.613120699860083)) : (xlims(plt), ylims(plt)))
    frame_width = xl[2] - xl[1]
    frame_height =yl[2] - yl[1]
    box_ratio = (w = 0.14 + 0.05 * abs(p.type.ξy*√2), h = 0.31) # hardcoded, I don't know how to get display length

    # lower left corner coordinates
    box = (x = (0.99 - box_ratio.w)*frame_width + xl[1],
           y = (0.99 - box_ratio.h)*frame_height + yl[1],
           w = box_ratio.w*frame_width,
           h = box_ratio.h*frame_height)

    # compute text positions
    text_x = box.x + box.w/2 # always centered on the box
    text_y = box.y + box.h # starting from the top, mind the minus
    δy = box.h/8

    # find the exact strings to print, related to polarization
    sgn = (imag(p.type.ξy) < 0 ? "-" : "")
    xi = (p.type.ξy == 0 ? "\$\\xi_y=0\$" : "\$\\xi_y= $sgn\\mathrm{i}/\\sqrt{2}\$")

    annotate!(plt, [(text_x, text_y - 1δy, Plots.text("\$a_0=$a₀\$",       :center, 10, "computer modern")),
                    (text_x, text_y - 3δy, Plots.text(xi,                  :center, 10, "computer modern")),
                    (text_x, text_y - 5δy, Plots.text("\$p=$(p.type.p)\$", :center, 10, "computer modern")),
                    (text_x, text_y - 7δy, Plots.text("\$m=$(p.type.m)\$", :center, 10, "computer modern"))])

    # add a rectangle without changing the frame size
    plot!(plt, rectangle(box.w, box.h, box.x, box.y);
                        xlims = xl, ylims = yl,
                        linewidth = 1.5,
                        fillcolor = :white, label = false)
end

function add_rect!(plt, X, Y, roots; kwargs...)
        h = ylims(plt)[2]-ylims(plt)[1]
        y₀ = ylims(plt)[1]
        underlay = plot([],[], label = false)
        for i in 1:(length(roots) - 1)
            a = roots[i]
            b = roots[i+1]
            plot!(underlay,rectangle(b-a, h, a, y₀),
            color_palette = :Set1_4, seriescolor = i,
            fillalpha = .4, linecolor = nothing, label = false)
        end
        plot!(underlay, X, Y; kwargs...)
end

function polarization(p)
    if imag(p.type.ξy)*√2 == 1.
        println("Right Circular Polarization, p = $(p.type.p), m = $(p.type.m)")
        return "right-handedCP"
    end
    if imag(p.type.ξy)*√2 == -1.
        println("Left Circular Polarization, p = $(p.type.p), m = $(p.type.m)")
        return "left-handedMCP"
    end
    if imag(p.type.ξy)*√2 == 0.
        println("Linear Polarization, p = $(p.type.p), m = $(p.type.m)")
        return "linearLP"
    end
end

for (sim_name, p) in pairs(field_param)
    folder_setup("$sim_name")
    folder_setup("fields")
    pol = polarization(p)

    # useful function definitions

    fieldE(x,y,z,t) = E([x,y,z],t,p.type)
    fieldB(x,y,z,t) = B([x,y,z],t,p.type)

    function wenergy(x,y,t)
        w = fieldE(x,y,0,t)⋅fieldE(x,y,0,t) + c*c * fieldB(x,y,0,t)⋅fieldB(x,y,0,t)
        if isnan(w) || !isfinite(w)
            println("Your field definition may have issues.")
            0
        else
            w
        end
    end

    wenergy₀(x,y) = wenergy(x,y,0)
    wenergy_period(x,y) = quadgk(t -> wenergy(x,y,t), 0, T₀, rtol = 1e-8)[1]

    wρ(ρ) = quadgk(ϕ -> wenergy₀(ρ*cos(ϕ), ρ*sin(ϕ))/(2*π), 0, 2*π, rtol = 1e-8)[1]

    Eρ(ρ,i) = quadgk(ϕ -> fieldE(ρ*cos(ϕ), ρ*sin(ϕ),0,0)[i]/(2*π), 0, 2*π, rtol = 1e-8)[1]

    xyEx(x,y) = fieldE(x,y,0,0)[1]
    xyEz(x,y) = fieldE(x,y,0,0)[3]

    jz₀(x,y,t = 0) = ([x,y,0.]×(fieldE(x,y,0,t)×fieldB(x,y,0,t)))[3]
    jzρ(ρ) = quadgk(ϕ -> jz₀(ρ*cos(ϕ), ρ*sin(ϕ))/(2*π), 0, 2*π, rtol = 1e-8, atol = 1.e-5)[1]

    # end of definitions

    # plots start here

    plot_data = [plot([], [], label = false) for i in 1:9]

    # 1. Heatmaps

    common_settings = (fontfamily = "computer modern", dpi = 360, label = false)
    plot_settings   = (color_palette = :Set1_4, seriescolor = 4, linewidth = 2.5)

    #  (x,y) -> (jz₀(w₀*x, w₀*y) + 0.00001)/(wenergy₀(w₀*x, w₀*y) + 0.00001)
    jz_period(x, y) = quadgk(t -> jz₀(x, y, t), 0, T₀, rtol = 1e-8, atol = 1e-8)[1]

    # plot_data[8] = heatmap(domainXY, domainXY,
    #             (x,y) -> 2ω*jz_period(w₀*x, w₀*y)/wenergy_period(w₀*x, w₀*y); #(any([(abs(hypot(x,y)-root) < 0.05) for root in p.roots]) ? (p.type.m + imag(p.type.ξy)*√2)/ω : jz₀(w₀*x, w₀*y, T₀/4)/wenergy(w₀*x, w₀*y, T₀/4));
    #             xlabel = L"x/w_0", ylabel = L"y/w_0",
    #             seriescolor = :lajolla,
    #             aspect_ratio = 1, common_settings...)
    # add_info!(plot_data[8], p, true)
    # savefig(plot_data[8], "j_w$(imag(p.type.ξy)*√2)$(p.type.p)$(p.type.m).png")

    jz_(x,y) = real((
                x*(
                    (E([x, y, 0], p.type)×conj(B([x, y, 0], p.type)))[2] +
                    (conj(E([x, y, 0], p.type))×B([x, y, 0], p.type))[2]
                    )
                    -
                y*(
                    (E([x, y, 0], p.type)×conj(B([x, y, 0], p.type)))[1] +
                    (conj(E([x, y, 0], p.type))×B([x, y, 0], p.type))[1]
                    ))
                    )/4

    w_(x,y) = real((E([x, y, 0], p.type)⋅conj(E([x, y, 0], p.type)) + c^2*conj(B([x, y, 0], p.type))⋅B([x, y, 0], p.type)))/2

    rjw(x,y) = (w_(x,y) < 1e-8 ? p.type.m + imag(p.type.ξy)*√2 : ω*jz_(x,y)/w_(x,y))

    plot_data[9] = heatmap(domainXY, domainXY,
                (x,y) -> rjw(w₀*x, w₀*y);
                xlabel = L"x/w_0", ylabel = L"y/w_0",
                seriescolor = :lajolla,
                aspect_ratio = 1, common_settings...)
    add_info!(plot_data[9], p, true)
    savefig(plot_data[9], "j_w_$(imag(p.type.ξy)*√2)$(p.type.p)$(p.type.m).png")

    plot_data[1] = heatmap(domainXY, domainXY, (x,y) -> wenergy₀(w₀*x, w₀*y);
                xlabel = L"x/w_0", ylabel = L"y/w_0",
                seriescolor = :linear_kryw_5_100_c67_n256,
                aspect_ratio = 1, common_settings...)
    add_info!(plot_data[1], p, true)
    savefig(plot_data[1], "xy00Energy$(imag(p.type.ξy)*√2)$(p.type.p)$(p.type.m).png")
    println("Energy density for time = 0, z = 0, in the xy-plane")

    plot_data[2] = heatmap(domainXY, domainXY, (x,y) -> xyEx(w₀*x, w₀*y);
                xlabel = L"x/w_0", ylabel = L"y/w_0",
                seriescolor = :vik,
                clims = x -> maximum(abs.(extrema(x))).*(-1,1),
                aspect_ratio = 1, common_settings...)
    add_info!(plot_data[2], p, true)
    savefig(plot_data[2], "xy00Ex$(imag(p.type.ξy)*√2)$(p.type.p)$(p.type.m).png")
    println("Electric field on the x-axis for time = 0, z = 0, in the xy-plane")

    plot_data[3] = heatmap(domainXY, domainXY, (x,y) -> xyEz(w₀*x, w₀*y);
                xlabel = L"x/w_0", ylabel = L"y/w_0",
                seriescolor = :vik,
                clims = x -> maximum(abs.(extrema(x))).*(-1,1),
                aspect_ratio = 1, common_settings...)
    add_info!(plot_data[3], p, true)
    savefig(plot_data[3], "xy00Ez$(imag(p.type.ξy)*√2)$(p.type.p)$(p.type.m).png")
    println("Electric field on the z-axis for time = 0, z = 0, in the xy-plane")

    plot_data[4] = heatmap(domainXY, domainXY, (x,y) -> jz₀(w₀*x, w₀*y);
                xlabel = L"x/w_0", ylabel = L"y/w_0",
                seriescolor = :vik,
                clims = x -> maximum(abs.(extrema(x))).*(-1,1),
                aspect_ratio = 1, common_settings...)
    add_info!(plot_data[4], p, true)
    savefig(plot_data[4], "xy00jz$(imag(p.type.ξy)*√2)$(p.type.p)$(p.type.m).png")
    println("Field angular momentum density on the z-axis for time = 0, z = 0, in the xy-plane")

    # 2. Azimuthally averaged radial plots

    jzρ_ = [jzρ(r*w₀) for r in domainR]
    plot_data[5] = plot(domainR, jzρ_, label = false)
    plot_data[5] = add_rect!(plot_data[5], domainR, jzρ_, p.roots;
                xlabel = L"\rho/w_0", ylabel = L"\overline{j_z^{xy}}(\rho, t = 0)\quad\mathrm{[a.u.]}",
                common_settings..., plot_settings...)
    add_info!(plot_data[5], p)
    savefig(plot_data[5], "rjz$(imag(p.type.ξy)*√2)$(p.type.p)$(p.type.m).png")
    println("Azimuthally averaged angular momentum on the z-axis for time = 0, z = 0, in the xy")

    wρ_ = [wρ(r*w₀) for r in domainR]
    plot_data[6] = plot(domainR, wρ_, label = false)
    plot_data[6] = add_rect!(plot_data[6], domainR, wρ_, p.roots;
        xlabel = L"\rho/w_0", ylabel = L"\mathrm{Intensity}(\rho)\quad\mathrm{[a.u.]}",
        common_settings..., plot_settings...)
    add_info!(plot_data[6], p)
    savefig(plot_data[6], "rw$(imag(p.type.ξy)*√2)$(p.type.p)$(p.type.m).png")
    println("Azimuthally averaged energy density for time = 0, z = 0, in the xy")

    # 3. Time axis plots.

    tE_ = [fieldE(w₀/√2/2, w₀/√2/2, 0, ti*T₀)[1] for ti in time]
    plot_data[7] = plot(time, tE_;
         xlabel = L"t/T_0", ylabel = L"E_x({w_0}/{\sqrt{8}}, {w_0}/{\sqrt{8}}, 0, t)\quad\mathrm{[a.u.]}",
         common_settings..., plot_settings...)
    add_info!(plot_data[7], p)
    savefig(plot_data[7],"tE$(imag(p.type.ξy)*√2)$(p.type.p)$(p.type.m).png")
    println("Electric field on the x-axis in time for x = w₀/√2/2, y = w₀/√2/2, z = 0")

    serialize("plot_data_$(pol)_$(p.type.p)_$(p.type.m)", plot_data)
    cd("..")
    cd("..")
end

using Telegram
using Telegram.API: sendAnimation, sendMessage, sendSticker
using ConfigEnv
dotenv("../.env")

sendAnimation(animation = "CgACAgQAAxkBAAMWYI1V67io3ILX4lFH-PzAZmoGfugAAi0CAAIvno1SOEmh8Ay-BCYfBA", caption = "Hei, tu!", disable_notification = true)
sendMessage(text = "Fișierul `$(basename(@__FILE__))`, de pe $(ENV["HOSTNAME"]) a fost executat cu succes\\! 💯⌛🧐", parse_mode = "MarkdownV2")
sendSticker(sticker ="CAACAgIAAxkBAAMXYI1V9sq0HUcFb7jdWxkWC5bboqoAAg4AA-nYEygTpj1DX_hIHx8E", caption = "Totul a mers ok!", disable_notification = true)

max_stars_repo_name: toma-vlad/electrodynamic-simulations

max_issues_repo_name: toma-vlad/electrodynamic-simulations

max_forks_repo_name: toma-vlad/electrodynamic-simulations