text:
@defcomp AdaptationCosts begin
    region = Index()

    y_year_0 = Parameter(unit="year")
    y_year = Parameter(index=[time], unit="year")
    gdp = Parameter(index=[time, region], unit="\$M")
    cf_costregional = Parameter(index=[region], unit="none") # first value should be 1.

    automult_autonomouschange = Parameter(unit="none")
    impmax_maximumadaptivecapacity = Parameter(index=[region], unit="driver")
    #tolerability parameters
    plateau_increaseintolerableplateaufromadaptation = Parameter(index=[region], unit="driver")
    pstart_startdateofadaptpolicy = Parameter(index=[region], unit="year")
    pyears_yearstilfulleffect = Parameter(index=[region], unit="year")
    impred_eventualpercentreduction = Parameter(index=[region], unit= "%")
    istart_startdate = Parameter(index=[region], unit = "year")
    iyears_yearstilfulleffect = Parameter(index=[region], unit= "year")

    cp_costplateau_eu = Parameter(unit="%GDP/driver")
    ci_costimpact_eu = Parameter(unit="%GDP/%driver")

    atl_adjustedtolerablelevel = Variable(index=[time, region]) # Unit depends on instance (degreeC or m)
    imp_adaptedimpacts = Variable(index=[time, region], unit="%")

    # Mostly for debugging
    autofac_autonomouschangefraction = Variable(index=[time], unit="none")
    acp_adaptivecostplateau = Variable(index=[time, region], unit="\$million")
    aci_adaptivecostimpact = Variable(index=[time, region], unit="\$million")

    ac_adaptivecosts = Variable(index=[time, region], unit="\$million")

    function run_timestep(p, v, d, tt)

        # Hope (2009), p. 21, equation -5
        auto_autonomouschangepercent = (1 - p.automult_autonomouschange^(1/(p.y_year[end] - p.y_year_0)))*100 # % per year
        v.autofac_autonomouschangefraction[tt] = (1 - auto_autonomouschangepercent/100)^(p.y_year[tt] - p.y_year_0) # Varies by year

        for rr in d.region
            #calculate adjusted tolerable level and max impact based on adaptation policy
            if (p.y_year[tt] - p.pstart_startdateofadaptpolicy[rr]) < 0
                v.atl_adjustedtolerablelevel[tt,rr]= 0
            elseif ((p.y_year[tt]-p.pstart_startdateofadaptpolicy[rr])/p.pyears_yearstilfulleffect[rr])<1.
                v.atl_adjustedtolerablelevel[tt,rr]=
                    ((p.y_year[tt]-p.pstart_startdateofadaptpolicy[rr])/p.pyears_yearstilfulleffect[rr]) *
                    p.plateau_increaseintolerableplateaufromadaptation[rr]
            else
                v.atl_adjustedtolerablelevel[tt,rr] = p.plateau_increaseintolerableplateaufromadaptation[rr]
            end

            if (p.y_year[tt]- p.istart_startdate[rr]) < 0
                v.imp_adaptedimpacts[tt,rr] = 0
            elseif ((p.y_year[tt]-p.istart_startdate[rr])/p.iyears_yearstilfulleffect[rr]) < 1
                v.imp_adaptedimpacts[tt,rr] =
                    (p.y_year[tt]-p.istart_startdate[rr])/p.iyears_yearstilfulleffect[rr]*
                    p.impred_eventualpercentreduction[rr]
            else
                v.imp_adaptedimpacts[tt,rr] = p.impred_eventualpercentreduction[rr]
            end

            # Hope (2009), p. 25, equations 1-2
            cp_costplateau_regional = p.cp_costplateau_eu * p.cf_costregional[rr]
            ci_costimpact_regional = p.ci_costimpact_eu * p.cf_costregional[rr]

            # Hope (2009), p. 25, equations 3-4
            v.acp_adaptivecostplateau[tt, rr] = v.atl_adjustedtolerablelevel[tt, rr] * cp_costplateau_regional * p.gdp[tt, rr] * v.autofac_autonomouschangefraction[tt] / 100
            v.aci_adaptivecostimpact[tt, rr] = v.imp_adaptedimpacts[tt, rr] * ci_costimpact_regional * p.gdp[tt, rr] * p.impmax_maximumadaptivecapacity[rr] * v.autofac_autonomouschangefraction[tt] / 100

            # Hope (2009), p. 25, equation 5
            v.ac_adaptivecosts[tt, rr] = v.acp_adaptivecostplateau[tt, rr] + v.aci_adaptivecostimpact[tt, rr]
        end
    end
end

function addadaptationcosts_sealevel(model::Model)
    adaptationcosts = add_comp!(model, AdaptationCosts, :AdaptiveCostsSeaLevel)

    # Sea Level-specific parameters
    setdistinctparameter(model, :AdaptiveCostsSeaLevel, :impmax_maximumadaptivecapacity, readpagedata(model, "data/impmax_sealevel.csv"))
    setdistinctparameter(model, :AdaptiveCostsSeaLevel, :plateau_increaseintolerableplateaufromadaptation, readpagedata(model, "data/sealevel_plateau.csv"))
    setdistinctparameter(model, :AdaptiveCostsSeaLevel, :pstart_startdateofadaptpolicy, readpagedata(model, "data/sealeveladaptstart.csv"))
    setdistinctparameter(model, :AdaptiveCostsSeaLevel, :pyears_yearstilfulleffect, readpagedata(model, "data/sealeveladapttimetoeffect.csv"))
    setdistinctparameter(model, :AdaptiveCostsSeaLevel, :impred_eventualpercentreduction, readpagedata(model, "data/sealevelimpactreduction.csv"))
    setdistinctparameter(model, :AdaptiveCostsSeaLevel, :istart_startdate, readpagedata(model, "data/sealevelimpactstart.csv"))
    setdistinctparameter(model, :AdaptiveCostsSeaLevel, :iyears_yearstilfulleffect, readpagedata(model, "data/sealevelimpactyearstoeffect.csv"))
    setdistinctparameter(model, :AdaptiveCostsSeaLevel, :cp_costplateau_eu, 0.0233333333333333)
    setdistinctparameter(model, :AdaptiveCostsSeaLevel, :ci_costimpact_eu, 0.00116666666666667)

    return adaptationcosts
end

function addadaptationcosts_economic(model::Model)
    adaptationcosts = add_comp!(model, AdaptationCosts, :AdaptiveCostsEconomic)

    # Economic-specific parameters
    setdistinctparameter(model, :AdaptiveCostsEconomic, :impmax_maximumadaptivecapacity, readpagedata(model, "data/impmax_economic.csv"))
    setdistinctparameter(model, :AdaptiveCostsEconomic, :plateau_increaseintolerableplateaufromadaptation, readpagedata(model, "data/plateau_increaseintolerableplateaufromadaptationM.csv"))
    setdistinctparameter(model, :AdaptiveCostsEconomic, :pstart_startdateofadaptpolicy, readpagedata(model, "data/pstart_startdateofadaptpolicyM.csv"))
    setdistinctparameter(model, :AdaptiveCostsEconomic, :pyears_yearstilfulleffect, readpagedata(model, "data/pyears_yearstilfulleffectM.csv"))
    setdistinctparameter(model, :AdaptiveCostsEconomic, :impred_eventualpercentreduction, readpagedata(model, "data/impred_eventualpercentreductionM.csv"))
    setdistinctparameter(model, :AdaptiveCostsEconomic, :istart_startdate, readpagedata(model, "data/istart_startdateM.csv"))
    setdistinctparameter(model, :AdaptiveCostsEconomic, :iyears_yearstilfulleffect, readpagedata(model, "data/iyears_yearstilfulleffectM.csv"))
    setdistinctparameter(model, :AdaptiveCostsEconomic, :cp_costplateau_eu, 0.0116666666666667)
    setdistinctparameter(model, :AdaptiveCostsEconomic, :ci_costimpact_eu, 0.0040000000)

    return adaptationcosts
end

function addadaptationcosts_noneconomic(model::Model)
    adaptationcosts = add_comp!(model, AdaptationCosts, :AdaptiveCostsNonEconomic)

    # Non-economic-specific parameters
    setdistinctparameter(model, :AdaptiveCostsNonEconomic, :impmax_maximumadaptivecapacity, readpagedata(model, "data/impmax_noneconomic.csv"))
    setdistinctparameter(model, :AdaptiveCostsNonEconomic, :plateau_increaseintolerableplateaufromadaptation, readpagedata(model, "data/plateau_increaseintolerableplateaufromadaptationNM.csv"))
    setdistinctparameter(model, :AdaptiveCostsNonEconomic, :pstart_startdateofadaptpolicy, readpagedata(model, "data/pstart_startdateofadaptpolicyNM.csv"))
    setdistinctparameter(model, :AdaptiveCostsNonEconomic, :pyears_yearstilfulleffect, readpagedata(model, "data/pyears_yearstilfulleffectNM.csv"))
    setdistinctparameter(model, :AdaptiveCostsNonEconomic, :impred_eventualpercentreduction, readpagedata(model, "data/impred_eventualpercentreductionNM.csv"))
    setdistinctparameter(model, :AdaptiveCostsNonEconomic, :istart_startdate, readpagedata(model, "data/istart_startdateNM.csv"))
    setdistinctparameter(model, :AdaptiveCostsNonEconomic, :iyears_yearstilfulleffect, readpagedata(model, "data/iyears_yearstilfulleffectNM.csv"))
    setdistinctparameter(model, :AdaptiveCostsNonEconomic, :cp_costplateau_eu, 0.0233333333333333)
    setdistinctparameter(model, :AdaptiveCostsNonEconomic, :ci_costimpact_eu, 0.00566666666666667)

    return adaptationcosts
end


text:
"""
First update is phi = 0, so BFGS for hessian
"""
function first_update(X, p, q)
    X + p*p'/dot(p, q) - X*q*q'*X/(q'*X*q)
end
function second_update(X, p, q)
    dqp = dot(q, p)
    (I - p*q'/dqp)*X*(I - q*p'/dqp) + p*p'/dqp
end
function first_update!(A, p, q)
    A .+= p*p'/dot(p, q) - A*(q*q')*A/(q'*A*q)
end
function second_update!(A, p, q)
    dqp = dot(q, p)
    # so right now, we just skip the update if dqp is zero
    # but we might do something else here
    if dqp > 0.0
        A .= (I - p*q'/dqp)*A*(I - q*p'/dqp) + p*p'/dqp
    end
    A
end

ext: jl


text:
#!/usr/local/bin/julia
#
const PKGS = ["Distributions","PyPlot","Plots","FFTW","DSP","OrdinaryDiffEq","Sundials","Roots",
              "Sundials","DifferentialEquations","Calculus","ForwardDiff","SymPy","QuadGK",
              "HCubature","JuMP","Clp","GLPK","Optim","SimJulia","ResumableFunctions"]

println("Installing required packages ...")

using Pkg
ipks = Pkg.installed()
for p in PKGS
    println("> $p")
    !haskey(ipks,p) && Pkg.add(p)
end
println("Done!")


text:
struct Identity{M,N,S} <: TaskMap{M,N,S} end
struct IdentityT{M,N,S} <: TaskMapT{M,N,S}
    base_map::Identity{M,N,S}
end
IdentityT{M,N,S}() where {M,N,S} = IdentityT{M,N,S}(Identity{M,N,S}())

task_map_emb(::EmbRep, ::EmbRep, xm, task_map::Identity{ℝ{m},ℝ{m}}) where m = xm

task_map_emb(::EmbRep, ::EmbRep, xme, task_map::Identity{𝕊{m},𝕊{m}}) where m = xme
domain_coord_rep(::Identity{𝕊{n},𝕊{n}}) where n = EmbRep()
codomain_coord_rep(::Identity{𝕊{n},𝕊{n}}) where n = EmbRep()

task_map_chart(::ChartRep, ::ChartRep, xm, task_map::Identity{SO3,SO3}, CM, CN) = xm
domain_coord_rep(::Identity{SO3,SO3}) = ChartRep()
codomain_coord_rep(::Identity{SO3,SO3}) = ChartRep()


text:
module TestPutGet

using StaticStorages
using StaticStorages.Implementations: DEFAULT_BUCKET_KEY, LOCAL_BUCKETS_NAME
using Test

const SPECIAL_BUCKET_KEY = StaticStorages.BucketKey()

const KEY1 = StaticStorages.put!(@__MODULE__, gensym(:value))
const KEY2 = StaticStorages.put!(@__MODULE__, [])
const KEY3 = StaticStorages.put!(@__MODULE__, SPECIAL_BUCKET_KEY, [])

const VALUE1 = StaticStorages.get(KEY1)
const VALUE2 = StaticStorages.get(KEY2)
const VALUE3 = StaticStorages.get(SPECIAL_BUCKET_KEY, KEY3)

function test_get_hit()
    @test StaticStorages.get(KEY1) === VALUE1
    @test StaticStorages.get(KEY2) === VALUE2
    @test StaticStorages.get(SPECIAL_BUCKET_KEY, KEY3) === VALUE3
end

function test_get_miss()
    @test StaticStorages.get(SPECIAL_BUCKET_KEY, KEY1) === nothing
    @test StaticStorages.get(SPECIAL_BUCKET_KEY, KEY2) === nothing
    @test StaticStorages.get(KEY3) === nothing
end

macro local_buckets()
    esc(LOCAL_BUCKETS_NAME)
end

local_buckets() = @local_buckets
local_storages() = get(local_buckets(), DEFAULT_BUCKET_KEY, nothing)

function test_local_buckets()
    @test local_storages()[KEY1] === VALUE1
    @test local_storages()[KEY2] === VALUE2
    @test local_buckets()[SPECIAL_BUCKET_KEY][KEY3] === VALUE3
end

function test_getbucket()
    @test StaticStorages.getbucket(SPECIAL_BUCKET_KEY) == Dict(KEY3 => VALUE3)
end

end  # module


text:
module Plug
using LightXML, Printf
using DifferentialEquations, Sundials
include("Constants.jl")
include("Utils.jl")
include("SurfaceReactions.jl")

"""
Plug flow
"""
function plug(input_file::AbstractString, lib_dir::AbstractString, sens= false)

    #locate the lib directory

    thermo_file = lib_dir*"therm.dat"

    local gasphase = Array{String,1}

    xmldoc = parse_file(input_file)
    xmlroot = root(xmldoc)

    #Get gasphase from the xml input
    gasphase = get_collection_from_xml(xmlroot,"gasphase")
    thermo_obj = SurfaceReactions.IdealGas.create_thermo(gasphase,thermo_file)

    #Get the molefractions
    mole_fracs = get_molefraction_from_xml(xmlroot, thermo_obj.molwt, gasphase)

    #Convert mole fractions to mass fractions
    mass_fracs = zeros( length(gasphase))
    molefrac_to_massfrac!(mass_fracs,mole_fracs,thermo_obj.molwt)

    #Read the inlet temperature
    local T = get_value_from_xml(xmlroot,"T")

    #Read the velocity
    local u = get_value_from_xml(xmlroot,"u")

    #Read the pressure
    local p = get_value_from_xml(xmlroot,"p")

    #Read the length of the reactor
    local l = get_value_from_xml(xmlroot,"length")

    #Read the diameter of the reactor
    local dia = get_value_from_xml(xmlroot,"dia")
    local Ac = π*dia^2/4.0
    local As_per_unit_length = π*dia

    #Read the wall temperature
    local Tw = get_value_from_xml(xmlroot,"Tw")

    #Read the surface area per unit length
    local cat_geom = get_value_from_xml(xmlroot,"cat-geom-factor")

    #Read the temperature condition
    isothermal = lowercase(get_text_from_xml(xmlroot,"isothermal")) == "true" ? true : false

    #Get the mechanism file from xml
    local mech_file = get_text_from_xml(xmlroot,"surface_mech")
    mech_file = lib_dir*"/"*mech_file

    #Create the mechanism definition
    md = SurfaceReactions.compile_mech(mech_file,thermo_obj,gasphase)

    #create the soln vector
    soln = copy(mass_fracs) #mass fractions
    push!(soln,density(mole_fracs,thermo_obj.molwt,T,p)*u) #mass flux ρu
    if !isothermal
        push!(soln,T)
    end

    #file output streams for saving the data
    g_stream = open("gas_profile.dat","w")
    s_stream = open("surf_profile.dat","w")
    os_streams = (g_stream,s_stream)
    geometry = (Ac,As_per_unit_length,cat_geom)

    create_header(g_stream,["z","T","p","u","rho"],gasphase)
    create_header(s_stream,"z","T",md.sm.species)
    #Set up the problem
    t_span = (0,l)
    covg = md.sm.si.ini_covg
    n_species = length(gasphase)+length(md.sm.species)
    rate = zeros(n_species)
    #Create the state object
    state = SurfaceReactions.State(mole_fracs,covg,T,p,rate)
    #this step is to get the steady state coverage before starting the plug flow integration
    t, state = SurfaceReactions.calculate_ss_molar_production_rates!(state,thermo_obj,md,1.0)
    params = (state,thermo_obj,md, geometry, os_streams)

    prob = ODEProblem(residual!,soln,t_span,params)
    if sens == true
        return (params, prob, t_span)
    end
    cb = FunctionCallingCallback(write_out_put)
    sol = solve(prob, CVODE_BDF(linear_solver=:GMRES), reltol=1e-6, abstol=1e-10, save_everystep=false,callback=cb);

    close(s_stream)
    close(g_stream)

    return sol.retcode
end

"""
residual function that defined the governing equations
"""
function residual!(du,u,p,t)
    #unpack the parameters
    state,thermo_obj,md,geom,os_stream = p
    Ac,Aspul,cat_geom = geom

    #Convert massfractions to molefractions
    ng = length(state.mole_frac)
    massfracs = u[1:ng]

    massfrac_to_molefrac!(state.mole_frac,massfracs,thermo_obj.molwt)

    #Calculate the molar production rates
    SurfaceReactions.calculate_ss_molar_production_rates!(state,thermo_obj,md,1.0)

    #species residual
    du[1:ng] = (state.s_rate[1:ng] .* thermo_obj.molwt) * cat_geom*(Aspul/Ac)/u[ng+1]

    #mass flux
    du[ng+1] = sum(state.s_rate[1:ng] .* thermo_obj.molwt)*cat_geom*Aspul/Ac
end

function write_out_put(u,t,integrator)
    state = integrator.p[1]
    thermo_obj = integrator.p[2]
    g_stream, s_stream = integrator.p[5]

    d = density(state.mole_frac,thermo_obj.molwt,state.T,state.p)
    vel = u[length(state.mole_frac)+1]/d
    write_to_file(g_stream,t,state.T,state.p,vel,d,state.mole_frac)
    write_to_file(s_stream,t,state.T,state.covg)
    @printf("%.4e\n", t)
end

end


text:
abstract type AbstractSlackBridge{T, VF, ZS, F, S} <: AbstractBridge end

function MOIB.added_constrained_variable_types(::Type{<:AbstractSlackBridge{T, VF, ZS, F, S}}) where {T, VF, ZS, F, S}
    return [(S,)]
end
function MOIB.added_constraint_types(::Type{<:AbstractSlackBridge{T, VF, ZS, F}}) where {T, VF, ZS, F}
    return [(F, ZS)]
end

function MOI.get(::AbstractSlackBridge{T, VF, ZS, F},
                 ::MOI.NumberOfConstraints{F, ZS}) where {T, VF, ZS, F}
    return 1
end
function MOI.get(::AbstractSlackBridge{T, VF, ZS, F, S},
                 ::MOI.NumberOfConstraints{VF, S}) where {T, VF, ZS, F, S}
    return 1
end
function MOI.get(bridge::AbstractSlackBridge{T, VF, ZS, F},
                 ::MOI.ListOfConstraintIndices{F, ZS}) where {T, VF, ZS, F}
    return [bridge.equality]
end
function MOI.get(bridge::AbstractSlackBridge{T, VF, ZS, F, S},
                 ::MOI.ListOfConstraintIndices{VF, S}) where {T, VF, ZS, F, S}
    return [bridge.slack_in_set]
end

# Indices
function MOI.delete(model::MOI.ModelLike, bridge::AbstractSlackBridge)
    MOI.delete(model, bridge.equality)
    MOI.delete(model, bridge.slack)
    return
end

# Attributes, Bridge acting as a constraint
function MOI.supports(
    ::MOI.ModelLike,
    ::Union{MOI.ConstraintPrimalStart, MOI.ConstraintDualStart},
    ::Type{<:AbstractSlackBridge})

    return true
end
function MOI.get(model::MOI.ModelLike,
                 attr::Union{MOI.ConstraintPrimal, MOI.ConstraintPrimalStart},
                 bridge::AbstractSlackBridge)
    # due to equality, slack should have the same value as original affine function
    return MOI.get(model, attr, bridge.slack_in_set)
end
function MOI.set(model::MOI.ModelLike,
                 attr::MOI.ConstraintPrimalStart,
                 bridge::AbstractSlackBridge,
                 value)
    if bridge isa ScalarSlackBridge
        MOI.set(model, MOI.VariablePrimalStart(), bridge.slack, value)
    else
        MOI.set.(model, MOI.VariablePrimalStart(), bridge.slack, value)
    end
    MOI.set(model, attr, bridge.slack_in_set, value)
    MOI.set(model, attr, bridge.equality, zero(value))
end
function MOI.get(model::MOI.ModelLike,
                 a::Union{MOI.ConstraintDual, MOI.ConstraintDualStart},
                 bridge::AbstractSlackBridge)
    # The dual constraint on slack (since it is free) is
    # -dual_slack_in_set + dual_equality = 0 so the two duals are
    # equal and we can return either one of them.
    return MOI.get(model, a, bridge.slack_in_set)
end
function MOI.set(model::MOI.ModelLike,
                 attr::MOI.ConstraintDualStart,
                 bridge::AbstractSlackBridge,
                 value)
    # As the slack appears `+slack` in `slack_in_set` and `-slack` in equality,
    # giving `value` to both will cancel it out in the Lagrangian.
    # Giving `value` to `bridge.equality` will put the function in the
    # lagrangian as expected.
    MOI.set(model, attr, bridge.slack_in_set, value)
    MOI.set(model, attr, bridge.equality, value)
end

function MOI.modify(model::MOI.ModelLike, bridge::AbstractSlackBridge,
                    change::MOI.AbstractFunctionModification)
    MOI.modify(model, bridge.equality, change)
end

function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintSet,
                 bridge::AbstractSlackBridge{T, VF, ZS, F, S}, change::S) where {T, VF, ZS, F, S}
    MOI.set(model, MOI.ConstraintSet(), bridge.slack_in_set, change)
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintFunction,
                 bridge::AbstractSlackBridge)
    return MOIU.remove_variable(MOI.get(model, attr, bridge.equality), bridge.slack)
end
function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintSet,
                 bridge::AbstractSlackBridge)
    return MOI.get(model, attr, bridge.slack_in_set)
end

# scalar version

"""
    ScalarSlackBridge{T, F, S}

The `ScalarSlackBridge` converts a constraint `G`-in-`S` where `G` is a function different
from `SingleVariable` into the constraints `F`-in-`EqualTo{T}` and `SingleVariable`-in-`S`.
`F` is the result of subtracting a `SingleVariable` from `G`.
Typically `G` is the same as `F`, but that is not mandatory.
"""
struct ScalarSlackBridge{T, F, S} <: AbstractSlackBridge{T, MOI.SingleVariable, MOI.EqualTo{T}, F, S}
    slack::MOI.VariableIndex
    slack_in_set::CI{MOI.SingleVariable, S}
    equality::CI{F, MOI.EqualTo{T}}
end
function bridge_constraint(::Type{ScalarSlackBridge{T, F, S}}, model,
                           f::MOI.AbstractScalarFunction, s::S) where {T, F, S}
    slack, slack_in_set = MOI.add_constrained_variable(model, s)
    new_f = MOIU.operate(-, T, f, MOI.SingleVariable(slack))
    equality = MOI.add_constraint(model, new_f, MOI.EqualTo(zero(T)))
    return ScalarSlackBridge{T, F, S}(slack, slack_in_set, equality)
end

# start allowing everything (scalar)
function MOI.supports_constraint(
    ::Type{ScalarSlackBridge{T}},
    F::Type{<:MOI.AbstractScalarFunction},
    ::Type{<:MOI.AbstractScalarSet}) where {T}
    return MOIU.is_coefficient_type(F, T)
end
# then restrict (careful with ambiguity)
MOI.supports_constraint(::Type{ScalarSlackBridge{T}},
                        ::Type{<:MOI.SingleVariable},
                        ::Type{<:MOI.EqualTo}) where {T} = false
MOI.supports_constraint(::Type{ScalarSlackBridge{T}},
                        ::Type{<:MOI.SingleVariable},
                        ::Type{<:MOI.AbstractScalarSet}) where {T} = false
MOI.supports_constraint(::Type{ScalarSlackBridge{T}},
                        ::Type{<:MOI.AbstractScalarFunction},
                        ::Type{<:MOI.EqualTo}) where {T} = false
function concrete_bridge_type(::Type{<:ScalarSlackBridge{T}},
                              F::Type{<:MOI.AbstractScalarFunction},
                              S::Type{<:MOI.AbstractScalarSet}) where T
    F2 = MOIU.promote_operation(-, T, F, MOI.SingleVariable)
    return ScalarSlackBridge{T, F2, S}
end

# Attributes, Bridge acting as a model
MOI.get(b::ScalarSlackBridge, ::MOI.NumberOfVariables) = 1
MOI.get(b::ScalarSlackBridge, ::MOI.ListOfVariableIndices) = [b.slack]

# Attributes, Bridge acting as a constraint
function MOI.get(model::MOI.ModelLike, ::MOI.ConstraintBasisStatus,  bridge::ScalarSlackBridge)
    MOI.get(model, MOI.ConstraintBasisStatus(), bridge.slack_in_set)
end

function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintFunction,
                 bridge::ScalarSlackBridge{T, F, S}, func::F) where {T, F, S}
    new_func = MOIU.operate(-, T, func, MOI.SingleVariable(bridge.slack))
    MOI.set(model, MOI.ConstraintFunction(), bridge.equality, new_func)
end

# vector version

"""
    VectorSlackBridge{T, F, S}

The `VectorSlackBridge` converts a constraint `G`-in-`S` where `G` is a function different
from `VectorOfVariables` into the constraints `F`in-`Zeros` and `VectorOfVariables`-in-`S`.
`F` is the result of subtracting a `VectorOfVariables` from `G`.
Tipically `G` is the same as `F`, but that is not mandatory.
"""
struct VectorSlackBridge{T, F, S} <: AbstractSlackBridge{T, MOI.VectorOfVariables, MOI.Zeros, F, S}
    slack::Vector{MOI.VariableIndex}
    slack_in_set::CI{MOI.VectorOfVariables, S}
    equality::CI{F, MOI.Zeros}
end
function bridge_constraint(::Type{VectorSlackBridge{T, F, S}}, model,
                           f::MOI.AbstractVectorFunction, s::S) where {T, F, S}
    d = MOI.dimension(s)
    slack, slack_in_set = MOI.add_constrained_variables(model, s)
    new_f = MOIU.operate(-, T, f, MOI.VectorOfVariables(slack))
    equality = MOI.add_constraint(model, new_f, MOI.Zeros(d))
    return VectorSlackBridge{T, F, S}(slack, slack_in_set, equality)
end

function MOI.supports_constraint(
    ::Type{VectorSlackBridge{T}},
    F::Type{<:MOI.AbstractVectorFunction},
    ::Type{<:MOI.AbstractVectorSet}) where {T}
    return MOIU.is_coefficient_type(F, T)
end
MOI.supports_constraint(::Type{VectorSlackBridge{T}},
                        ::Type{<:MOI.VectorOfVariables},
                        ::Type{<:MOI.Zeros}) where {T} = false
MOI.supports_constraint(::Type{VectorSlackBridge{T}},
                        ::Type{<:MOI.AbstractVectorFunction},
                        ::Type{<:MOI.Zeros}) where {T} = false
MOI.supports_constraint(::Type{VectorSlackBridge{T}},
                        ::Type{<:MOI.VectorOfVariables},
                        ::Type{<:MOI.AbstractVectorSet}) where {T} = false
function concrete_bridge_type(::Type{<:VectorSlackBridge{T}},
                              F::Type{<:MOI.AbstractVectorFunction},
                              S::Type{<:MOI.AbstractVectorSet}) where T
    F2 = MOIU.promote_operation(-, T, F, MOI.VectorOfVariables)
    return VectorSlackBridge{T, F2, S}
end

# Attributes, Bridge acting as a model
MOI.get(b::VectorSlackBridge, ::MOI.NumberOfVariables) = length(b.slack)
MOI.get(b::VectorSlackBridge, ::MOI.ListOfVariableIndices) = b.slack

# Attributes, Bridge acting as a constraint
function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintFunction,
                 bridge::VectorSlackBridge{T, F, S}, func::F) where {T, F, S}
    new_func = MOIU.operate(-, T, func, MOI.VectorAffineFunction{T}(MOI.VectorOfVariables(bridge.slack)))
    MOI.set(model, MOI.ConstraintFunction(), bridge.equality, new_func)
end


text:
using Test

using Example2018

## Polynomials form a ring
# 1. Addition is associative
# 2. Additive identity
# 3. Additive inverse
# 4. Additive commutes
## 5. Associative
# 5. Multiplication is ...

@test 2==2
@test 1==1

z = Polynomial([0])
p = Polynomial([1,2,3])
## Need to create function to compare polynomials, like substract them and find out it is a zero polynomial
@test p + z == p
@test z + p == p
@test p + (-p) == z

### 1,2,3 polynomial is not strong, now, let's introduce a test set

@testset "Polynomials form a ring" begin
z = Polynomial([0])

    for n in 1:100
        ##p = Polynomial([1,2,3])
        deg  = rand(0:100)
        p = Polynomial(rand(-1000000:100000, deg+1))
        @test p + z == p
        @test z + p == p
        @test p + (-p) == z
    end
end


text:
#!/usr/bin/env julia6
using ReadGlobal

@inline function energydistr_scalar(t11::T,t22::T,t33::T,t12::T,t13::T,t23::T,
    d11::T,d22::T,d33::T,d12::T,d13::T,d23::T,
    w12::T,w13::T,w23::T) where {T<:Number}

    #ɛ = t11 * d11 + t22 * d22 + t33 * d33 + 2.0*(t12 * d12 + t13*d13 + t23*d23)
    ɛ = muladd(t11, d11 ,muladd(t22, d22, muladd(t33, d33, 2.0*muladd(t12, d12, muladd(t13, d13, t23*d23)))))

    #dir11 = (d11*t11+d12*t12+d13*t13-t12*w12-t13*w13) - ɛ/3.
    dir11 = muladd(d11, t11, muladd(d12, t12, muladd(d13, t13, -muladd(t12, w12, muladd(t13, w13, ɛ/3)))))

    #dir22 = (d12*t12+d22*t22+d23*t23+t12*w12-t23*w23) - ɛ/3.
    dir22 = muladd(d12, t12, muladd(d22, t22, muladd(d23, t23, muladd(t12, w12, -muladd(t23, w23, ɛ/3)))))

    #dir33 = (d13*t13+d23*t23+d33*t33+t13*w13+t23*w23) - ɛ/3.
    dir33 = muladd(d13, t13, muladd(d23, t23, muladd(d33, t33, muladd(t13, w13, muladd(t23, w23, - ɛ/3)))))

    #dir12 = 0.5*(d11*t12+d12*t11+d12*t22+d13*t23+d22*t12+d23*t13+t11*w12-t13*w23-t22*w12-t23*w13)
    dir12 = 0.5*muladd(d11, t12, muladd(d12, t11, muladd(d12, t22, muladd(d13, t23,
               muladd(d22, t12, muladd(d23, t13, muladd(t11, w12, -muladd(t13, w23, muladd(t22, w12, t23*w13)))))))))

    #dir13 = 0.5*(d11*t13+d12*t23+d13*t11+d13*t33+d23*t12+d33*t13+t11*w13+t12*w23-t23*w12-t33*w13)
    dir13 = 0.5*muladd(d11, t13, muladd(d12, t23, muladd(d13, t11, muladd(d13, t33,
               muladd(d23, t12, muladd(d33, t13, muladd(t11, w13, muladd(t12, w23, -muladd(t23, w12, t33*w13)))))))))

    #dir23 = 0.5*(d12*t13+d13*t12+d22*t23+d23*t22+d23*t33+d33*t23+t12*w13+t13*w12+t22*w23-t33*w23)
    dir23 = 0.5*muladd(d12, t13, muladd(d13, t12, muladd(d22, t23, muladd(d23, t22,
               muladd(d23, t33, muladd(d33, t23, muladd(t12, w13, muladd(t13, w12, muladd(t22, w23, -t33*w23)))))))))
    return ε, dir11, dir22, dir33, dir12, dir13, dir23
end

function energydistr_part!(t11::T,t22::T,t33::T,t12::T,t13::T,t23::T,
                        d11::T,d22::T,d33::T,d12::T,d13::T,d23::T,
                        w12::T,w13::T,w23::T,
                        dir11::T,dir22::T,dir33::T,dir12::T,dir13::T,dir23::T,
                        ɛ::T) where {T<:Vector{<:Number}}

    Threads.@threads for i in 1:length(d11)

      @inbounds ɛ[i], dir11[i], dir22[i], dir33[i], dir12[i], dir13[i], dir23[i] = energydistr_scalar(t11[i],t22[i],t33[i],t12[i],t13[i],t23[i],d11[i],d22[i],d33[i],d12[i],d13[i],d23[i],w12[i],w13[i],w23[i])

    end
    return nothing
end

function energydistr()
  nx, ny, nz, xs, ys, zs = getdimsize()

  dim = (nx,ny,nz)

  place = split(pwd(),"/")
  if ismatch(r"^Model_",place[end])
    N = place[end-1][2:end]
    Fil = place[end-2][1]

    ft11 = "Tm11"
    ft22 = "Tm22"
    ft33 = "Tm33"
    ft12 = "Tm12"
    ft13 = "Tm13"
    ft23 = "Tm23"

    fd11 = "../$(Fil)D11_N$N"
    fd22 = "../$(Fil)D22_N$N"
    fd33 = "../$(Fil)D33_N$N"
    fd12 = "../$(Fil)D12_N$N"
    fd13 = "../$(Fil)D13_N$N"
    fd23 = "../$(Fil)D23_N$N"

    fw12 = "../$(Fil)W12_N$N"
    fw13 = "../$(Fil)W13_N$N"
    fw23 = "../$(Fil)W23_N$N"

  else

    N = place[end][2:end]
    Fil = place[end-1][1]

    ft11 = "$(Fil)T11_N$N"
    ft22 = "$(Fil)T22_N$N"
    ft33 = "$(Fil)T33_N$N"
    ft12 = "$(Fil)T12_N$N"
    ft13 = "$(Fil)T13_N$N"
    ft23 = "$(Fil)T23_N$N"

    fd11 = "$(Fil)D11_N$N"
    fd22 = "$(Fil)D22_N$N"
    fd33 = "$(Fil)D33_N$N"
    fd12 = "$(Fil)D12_N$N"
    fd13 = "$(Fil)D13_N$N"
    fd23 = "$(Fil)D23_N$N"

    fw12 = "$(Fil)W12_N$N"
    fw13 = "$(Fil)W13_N$N"
    fw23 = "$(Fil)W23_N$N"

  end

  fɛ = "Pr"

  fdir11 = "dir11"
  fdir22 = "dir22"
  fdir33 = "dir33"
  fdir12 = "dir12"
  fdir13 = "dir13"
  fdir23 = "dir23"

  doinchunks(energydistr_part!,
             input=(ft11,ft22,ft33,ft12,ft13,ft23,fd11,fd22,fd33,fd12,fd13,fd23,fw12,fw13,fw23),
             output=(fdir11,fdir22,fdir33,fdir12,fdir13,fdir23,fɛ))

end

energydistr()


text:
module BioMetaheuristics

import RandomNumbers.Xorshifts
using Printf: @printf
using StaticArrays
using Random
import Base: isless

# Benchmark functions module
export TestFunctions
include(joinpath("TestFunctions", "TestFunctions.jl"))

# Types
export Particle, Population, OptimizationResults, PopulationBase, Metaheuristic,
TrajectoryBase
include(joinpath("types", "population.jl"))
include(joinpath("types", "solvers.jl"))

# Some utilities
include("utils.jl")

# Solvers
export PSO, SimulatedAnnealing, GeneralSimulatedAnnealing
include(joinpath("solvers", "pso.jl"))
include(joinpath("solvers", "simulated_annealing.jl"))

# Common API
export optimize
include("optimize.jl")

end


text:
using Tokens
using Test

include("testpacked31.jl")
include("testliterals.jl")
include("testflytoken.jl")
include("testtoken.jl")
include("testioshared.jl")
include("testvector.jl")
include("testtree.jl")
include("testlexer.jl")

@testset "Tokens.jl" begin
    # Write your own tests here.
    testpacked31()
    testliterals()
    testflytoken()
    testtoken()
    testshared()
    testvector()
    testtree()
    testlexer()
end


text:
using Test
using MPI

if get(ENV,"JULIA_MPI_TEST_ARRAYTYPE","") == "CuArray"
    import CUDA
    ArrayType = CUDA.CuArray
else
    ArrayType = Array
end

MPI.Init()

comm = MPI.COMM_WORLD
size = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)

dst = mod(rank+1, size)
src = mod(rank-1, size)

N = 32

send_mesg = ArrayType{Float64}(undef,N)
recv_mesg = ArrayType{Float64}(undef,N)
recv_mesg_expected = ArrayType{Float64}(undef,N)

fill!(send_mesg, Float64(rank))
fill!(recv_mesg_expected, Float64(src))

rreq = MPI.Irecv!(recv_mesg, src,  src+32, comm)
sreq = MPI.Isend(send_mesg, dst, rank+32, comm)

stats = MPI.Waitall!([sreq, rreq])
@test rreq isa MPI.Request
@test sreq isa MPI.Request
@test MPI.Get_source(stats[2]) == src
@test MPI.Get_tag(stats[2]) == src+32
@test recv_mesg == recv_mesg_expected

(done, stats) = MPI.Testall!([sreq, rreq])
@test done

if rank == 0
    MPI.send(send_mesg, dst, rank+32, comm)
    recv_mesg = recv_mesg_expected
elseif rank == size-1
    (recv_mesg, _) = MPI.recv(src, src+32, comm)
else
    (recv_mesg, _) = MPI.recv(src, src+32, comm)
    MPI.send(send_mesg, dst, rank+32, comm)
end
@test recv_mesg == recv_mesg_expected

if rank == 0
    MPI.Send(Float64(rank), dst, rank+32, comm)
    recv_val = Float64(src)
elseif rank == size-1
    (recv_val, _) = MPI.Recv(Float64, src, src+32, comm)
else
    (recv_val, _) = MPI.Recv(Float64, src, src+32, comm)
    MPI.Send(Float64(rank), dst, rank+32, comm)
end
@test recv_val == Float64(src)

rreq = MPI.Irecv!(recv_mesg, src,  src+32, comm)
sreq = MPI.Isend(send_mesg, dst, rank+32, comm)

req_arr = [sreq,rreq]
inds = MPI.Waitsome(req_arr)
for i in inds
    @test MPI.Test(req_arr[i])
end

rreq = MPI.Irecv!(recv_mesg, src,  src+32, comm)
MPI.Cancel!(rreq)
MPI.Wait(rreq)
@test rreq.buffer == nothing

GC.gc()

# ---------------------
# MPI_Sendrecv function
# ---------------------
#
# send datatype
# ---------------------
# We test this function by executing a left shift of the leftmost element in a 1D
# cartesian topology with periodic boundary conditions.
#
# Assume we have two processors, the data will look like this
# proc 0 | proc 1
#  0 0 0 |  1 1 1
#
# After the shift the data will contain
# proc 0 | proc 1
#  0 0 1 |  1 1 0
#
# init data
comm_rank = MPI.Comm_rank(comm)
comm_size = MPI.Comm_size(comm)
a = Float64[comm_rank, comm_rank, comm_rank]

# construct cartesian communicator with 1D topology
comm_cart = MPI.Cart_create(comm, 1, Cint[comm_size], Cint[1], false)

# get source and dest ranks using Cart_shift
src_rank, dest_rank = MPI.Cart_shift(comm_cart, 0, -1)

# execute left shift using subarrays
MPI.Sendrecv!(@view(a[1]), dest_rank, 0,
              @view(a[3]), src_rank,  0, comm_cart)

@test a == [comm_rank, comm_rank, (comm_rank+1) % comm_size]

# send elements from a buffer
# ---------------------------
a = Float64[comm_rank, comm_rank, comm_rank]
b = Float64[       -1,        -1,        -1]
MPI.Sendrecv!(@view(a[1:2]), dest_rank, 1,
              @view(b[1:2]), src_rank,  1, comm_cart)

@test b == [(comm_rank+1) % comm_size, (comm_rank+1) % comm_size, -1]

# send entire buffer
# ---------------------------
a = Float64[comm_rank, comm_rank, comm_rank]
b = Float64[       -1,        -1,        -1]
MPI.Sendrecv!(a, dest_rank, 2,
              b,  src_rank, 2, comm_cart)

@test b == [(comm_rank+1) % comm_size, (comm_rank+1) % comm_size, (comm_rank+1) % comm_size]

@test MPI.Waitall!(MPI.Request[]) == MPI.Status[]

MPI.Finalize()
# @test MPI.Finalized()


text:
module DiffEqBenchmarks

using Weave, Pkg, IJulia, InteractiveUtils, Markdown

repo_directory = joinpath(@__DIR__,"..")

function weave_file(folder,file,build_list=(:script,:html,:pdf,:github,:notebook))
  println("File: $file")
  tmp = joinpath(repo_directory,"benchmarks",folder,file)
  args = Dict{Symbol,String}(:folder=>folder,:file=>file)
  if :script ∈ build_list
    println("Building Script")
    dir = joinpath(repo_directory,"script",folder)
    isdir(dir) || mkdir(dir)
    tangle(tmp;out_path=dir)
  end
  if :html ∈ build_list
    println("Building HTML")
    dir = joinpath(repo_directory,"html",folder)
    isdir(dir) || mkdir(dir)
    weave(tmp,doctype = "md2html",out_path=dir,args=args,fig_ext=".svg")
  end
  if :pdf ∈ build_list
    println("Building PDF")
    dir = joinpath(repo_directory,"pdf",folder)
    isdir(dir) || mkdir(dir)
    weave(tmp,doctype="md2pdf",out_path=dir,args=args)
  end
  if :github ∈ build_list
    println("Building Github Markdown")
    dir = joinpath(repo_directory,"markdown",folder)
    isdir(dir) || mkdir(dir)
    weave(tmp,doctype = "github",out_path=dir,args=args)
  end
  if :notebook ∈ build_list
    println("Building Notebook")
    dir = joinpath(repo_directory,"notebook",folder)
    isdir(dir) || mkdir(dir)
    Weave.convert_doc(tmp,joinpath(dir,file[1:end-4]*".ipynb"))
  end
end

function weave_all()
  for folder in readdir(joinpath(repo_directory,"benchmarks"))
    folder == "test.jmd" && continue
    weave_folder(folder)
  end
end

function weave_folder(folder)
  for file in readdir(joinpath(repo_directory,"benchmarks",folder))
    println("Building $(joinpath(folder,file)))")
    try
      weave_file(folder,file)
    catch
    end
  end
end

function bench_footer(folder=nothing, file=nothing)
    display(md"""
    ## Appendix

    These benchmarks are a part of the DiffEqBenchmarks.jl repository, found at: <https://github.com/JuliaDiffEq/DiffEqBenchmarks.jl>
    """)
    if folder !== nothing && file !== nothing
        display(Markdown.parse("""
        To locally run this tutorial, do the following commands:
        ```
        using DiffEqBenchmarks
        DiffEqBenchmarks.weave_file("$folder","$file")
        ```
        """))
    end
    display(md"Computer Information:")
    vinfo = sprint(InteractiveUtils.versioninfo)
    display(Markdown.parse("""
    ```
    $(vinfo)
    ```
    """))

    ctx = Pkg.API.Context()
    pkgs = Pkg.Display.status(Pkg.API.Context(), use_as_api=true);

    display(md"""
    Package Information:
    """)

    io = IOBuffer()
    for pkg in pkgs
        ver = pkg.new.ver === nothing ? "" : pkg.new.ver
        println(io, "[$(pkg.uuid)] $(pkg.name) $ver")
    end
    proj = String(take!(io))
    md = """
    ```
    Status: `$(ctx.env.project_file)`
    $(chomp(proj))
    ```
    """
    display(Markdown.parse(md))
end

function open_notebooks()
  Base.eval(Main, Meta.parse("import IJulia"))
  path = joinpath(repo_directory,"notebook")
  IJulia.notebook(;dir=path)
end

end


text:
import ClimaCore.Geometry, LinearAlgebra
import ClimaCore:
    Fields,
    Domains,
    Topologies,
    Meshes,
    DataLayouts,
    Operators,
    Geometry,
    Spaces

using OrdinaryDiffEq: OrdinaryDiffEq, ODEProblem, solve, SSPRK33

import Logging
import TerminalLoggers
Logging.global_logger(TerminalLoggers.TerminalLogger())
const FT = Float64

domain = Domains.IntervalDomain(
    Geometry.ZPoint{FT}(0.0),
    Geometry.ZPoint{FT}(4pi),
    boundary_names = (:left, :right),
)
mesh = Meshes.IntervalMesh(domain; nelems = 30)

cspace = Spaces.CenterFiniteDifferenceSpace(mesh)
fspace = Spaces.FaceFiniteDifferenceSpace(cspace)

zc = Fields.coordinate_field(cspace)
u = sin.(zc.z)
p = Geometry.WVector.(zeros(Float64, fspace))

Y = Fields.FieldVector(u = u, p = p)

function tendency!(dY, Y, _, t)
    u = Y.u
    p = Y.p

    du = dY.u
    dp = dY.p

    ∂f = Operators.GradientC2F(
        left = Operators.SetValue(0.0),
        right = Operators.SetValue(0.0),
    )
    ∂c = Operators.DivergenceF2C()

    @. dp = -Geometry.WVector(∂f(u))
    @. du = -∂c(p)

    return dY
end

@show tendency!(similar(Y), Y, nothing, 0.0)

# Solve the ODE operator
Δt = 0.01
prob = ODEProblem(tendency!, Y, (0.0, 4 * pi))
sol = solve(
    prob,
    SSPRK33(),
    dt = Δt,
    saveat = 10 * Δt,
    progress = true,
    progress_message = (dt, u, p, t) -> t,
);

ENV["GKSwstype"] = "nul"
using ClimaCorePlots, Plots
Plots.GRBackend()

dirname = "wave"
path = joinpath(@__DIR__, "output", dirname)
mkpath(path)

anim = Plots.@animate for u in sol.u
    Plots.plot(u.u, xlim = (-1, 1), label = "u", legend = true)
end
Plots.mp4(anim, joinpath(path, "wave.mp4"), fps = 10)

Plots.png(
    Plots.plot(sol.u[end].u, label = "u", legend = true),
    joinpath(path, "wave_end.png"),
)

function linkfig(figpath, alt = "")
    # buildkite-agent upload figpath
    # link figure in logs if we are running on CI
    if get(ENV, "BUILDKITE", "") == "true"
        artifact_url = "artifact://$figpath"
        print("\033]1338;url='$(artifact_url)';alt='$(alt)'\a\n")
    end
end

linkfig(
    relpath(joinpath(path, "wave_end.png"), joinpath(@__DIR__, "../..")),
    "Wave End",
)


text:
# for testing code that is not exported or included in Laplacians.jl

#============================================================#
# Test Lex by running code from its notebook
#=
fn = "$(@__DIR__)/../src/lex.jl"
println(fn)

include(fn)

setLexDebugFlag(true)

n = 10
Pn = pathGraph(n)

isTerm = zeros(Bool, n)
isTerm[1] = true
isTerm[n] = true

initVal = zeros(n)
initVal[n] = 1.0

# inf-minimizer
infMinVolt = CompInfMin(Pn, isTerm, initVal)
println(infMinVolt)
println(MaxEdgeGrad(Pn, infMinVolt))

# lex-minimizer
lexMinVolt = CompLexMin(Pn, isTerm, initVal)

println(lexMinVolt)

n = 100
G = chimera(n,1,ver=RandomV06.V06)
isTerm = zeros(Bool, n)
# arbitrary terminal values
isTerm[1] = true
isTerm[5] = true
isTerm[11] = true
isTerm[18] = true

initVal = zeros(Float64, n)
initVal[1] = 0.0
initVal[5] = 13
initVal[11] = 7
initVal[18] = 11

infMinVolt = CompInfMin(G, isTerm, initVal)
println(infMinVolt)
println(MaxEdgeGrad(G, infMinVolt))

lexMinVolt = simIterLex(500, G, isTerm, initVal)
println(lexMinVolt)
println(MaxEdgeGrad(G, lexMinVolt))
println(checkLex(G, isTerm, initVal, lexMinVolt))

simIterLexTest()
=#
#============================================================#
# test isotonicIPM

include("$(@__DIR__)/../src/isotonicIPM.jl")

n = 100;
A = grid2(10,10)
A = triu(A);
v = randn(n)+collect(1:n)/2;
isoTestValues = collect(1:n);
newOrder = randperm(n)
permMat = (sparse(I,n,n))[1:n,newOrder];
A = permMat*A*permMat';
v = permMat*v;
@time (x,acc,itercount) = isotonicIPM(A,v);


text:
#Variables

#Entero
entero = 999
println(entero)
println(typeof(entero))

#Flotante
flotante = 23.45
println(flotante)
println(typeof(flotante))

#String
cadenaCaracteres = "Dibit Learning"
println(cadenaCaracteres)
println(typeof(cadenaCaracteres))

#Caracter
caracter = 'K'
println(caracter)
println(typeof(caracter))

#Bool
booleano = true
println(booleano)
println(typeof(booleano))

#Caracteres especiales UTF-8
letrasRaras = "人人生而自由，在尊严和权利上一律平等。"
println(letrasRaras)
println(typeof(letrasRaras))

#Unicode names are allowed (UTF-8)
生而自由 = "Aprender es bonito"
println(生而自由)
println(typeof(生而自由))

#Constantes
const anioNacimiento = "2020"
println(anioNacimiento)
println(typeof(anioNacimiento))

#const es para el tipo de dato
anioNacimiento = "2021"
println(anioNacimiento)
println(typeof(anioNacimiento))

#¿Se puede redefinir el tipo de dato?
entero = "Ahora es string"
println(entero)
println(typeof(entero))

#anioNacimiento = 2021 No se puede dar
#const evita que se cambie el tipo de dato de una variable

#¿32-bit o 64-bit
println(Sys.WORD_SIZE)

#nombre::String = "Agreda"

#Fuera del alcance de una función
#variable::[tipo de dato] = [dato]
#Optimización
function prueba()
    nombre::String = "Agreda"
    println(nombre)
    println(typeof(nombre))
end

prueba()


text:
# I want to try a pure MC simulation that doesn't rely on a neural network and
# see how things
using Statistics: mean, quantile
using StatsBase: countmap
using Game2048
using Game2048: bitboard_to_array, DIRS
using Serialization: serialize, deserialize
using Base.Threads
using BSON
using BenchmarkTools
using Flux.Data: DataLoader
using Flux: throttle
using Flux: softmax, onecold
# make a nn to study these games
using Flux, CUDA
using IterTools
CUDA.allowscalar(false)

### load all results are analysis
full_results = mapreduce(deserialize, vcat, readdir("c:/data/game2048/raw/", join=true));

# countmap(map(x->maximum(x.final_state), full_results))
#   5  => 48
#   6  => 223
#   7  => 1120
#   11 => 21269
#   10 => 106756
#   9  => 41368
#   12 => 1
#   8  => 7115

# how many are bigger than 11
good_results = filter(x->maximum(x.final_state) >= 11, full_results);
length(good_results)
mapreduce(x->x.counter, +, good_results)

bad_results = filter(x->maximum(x.final_state) <=9, full_results);
mapreduce(x->x.counter, +, bad_results)

results_for_modelling = vcat(good_results, bad_results);

function make_x_for_value_modelling(results_for_modelling)
    # preallocate the memory
    res = Array{Int8}(undef, 4, 4, mapreduce(x->x.counter, +, results_for_modelling))
    i = 1
    for result in results_for_modelling
        for state in result.states
            res[:, :, i] .= bitboard_to_array(state)
            i += 1
        end
    end
    res
end

function make_y_for_value_modelling(results_for_modelling)
    # preallocate the memory
    res = Array{Int32}(undef, mapreduce(x->x.counter, +, results_for_modelling))
    i = 1
    for result in results_for_modelling
        val = sum(2 .<< bitboard_to_array(result.final_state))
        for _ in 1:result.counter
            res[i] = val
            i += 1
        end
    end
    res
end

@time x = reshape(make_x_for_value_modelling(results_for_modelling), 4, 4, 1, :);
@time y = make_y_for_value_modelling(results_for_modelling);

value_model = Chain(
    Conv((2, 2), 1=>8, relu),
    Conv((2, 2), 8=>8, relu),
    flatten,
    Dense(32, 1)
) #|> gpu

ly = length(y)
using StatsBase
idx1 = sort(sample(1:ly, round(Int, ly * 0.50), replace=false))
# idx2 = sort(sample(setdiff(1:ly, idx1), round(Int, ly * 0.1), replace=false))

xtrain = (x[:, :, :, idx1] ./ Float32(11)) #|> cu;
ytrain = Float32.(reshape(log.(y[idx1]), 1, :)) #|> cu;

using Statistics
ytrain = (ytrain .- mean(ytrain))/std(ytrain)

testidx = sample(setdiff(1:ly, idx1), round(Int, ly * 0.10), replace=false)
xtest = x[:, :, :, testidx] ./ Float32(11) #|> cu;
ytest = Float32.(reshape(log.(y[testidx]), 1, :)) #|> cu;
ytest = (ytest .- mean(ytest))/std(ytest)

loss(x, y) = Flux.mse(value_model(x), y; agg=sum)

loss(xtest, ytest)

using StatsPlots

using DataFrames
px = value_model(xtest)
df = DataFrame(px = reshape(px, :), ytest =reshape(ytest, :))

sort!(df, :ytest)

cb = Flux.throttle(()->println(loss(xtest, ytest)), 8)

opt= ADAM()

dl = DataLoader((xtrain, ytrain); batchsize=32);

@time Flux.train!(loss, Flux.params(value_model), dl, opt; cb)


text:
########################################################################

"""
$(TYPEDSIGNATURES)

Total magnetic field (2D) for a set of polygonal bodies defined by their corners. Takes into account both induced and remnant magnetization.
 Based on Talwani & Heitzler (1964), the default algorithm in Mag2Dpoly package.
"""
function tmagpolybodies2D(xzobs::Array{<:Real,2},Jinds::Vector{MagnetizVector},Jrems::Vector{MagnetizVector},
                    northxax::Real,bodies::MagPolygBodies2D)

    nbody = length(bodies.bo)
    tmag = zeros(eltype(Jinds[1].mod),size(xzobs,1))
    for i=1:nbody
        tmag .+= tmagpoly2D(xzobs,Jinds[i],Jrems[i],northxax,bodies.bo[i])
    end
    return tmag
end
###################################################################################

"""
$(TYPEDSIGNATURES)

Total magnetic field (2D) for a set of polygonal bodies defined by their corners. Takes into account both induced and remnant magnetization.
Generic version containing four different algorithm formulations `forwardtype`, passed as a string:
  - "talwani"      --> Talwani & Heitzler (1964)
  - "talwani_red"  --> Talwani & Heitzler (1964) rederived from Kravchinsky et al. 2019
  - "krav"         --> Kravchinsky et al. (2019) rectified by Ghirotto et al. (2020)
  - "wonbev"       --> Won & Bevis (1987)
"""
function tmagpolybodies2Dgen(xzobs::Array{<:Real,2},Jinds::Vector{MagnetizVector},Jrems::Vector{MagnetizVector},
                       northxax::Real,bodies::MagPolygBodies2D,forwardtype::String)

    nbody = length(bodies.bo)
    tmag = zeros(eltype(Jinds[1].mod),size(xzobs,1))
    for i=1:nbody
        tmag .+= tmagpoly2Dgen(xzobs,Jinds[i],Jrems[i],northxax,bodies.bo[i],forwardtype)
    end
    return tmag
end

###################################################################################

"""
$(TYPEDSIGNATURES)

Total magnetic field (2D) for a polygon defined by its corners. Takes into account both induced and remnant magnetization.
Generic version containing four different algorithm formulations `forwardtype`, passed as a string:
  - "talwani"      --> Talwani & Heitzler (1964)
  - "talwani_red"  --> Talwani & Heitzler (1964) rederived from Kravchinsky et al. 2019
  - "krav"         --> Kravchinsky et al. (2019) rectified by Ghirotto et al. (2020)
  - "wonbev"       --> Won & Bevis (1987)
"""
function tmagpoly2Dgen(xzobs::Array{<:Real,2},Jind::MagnetizVector,Jrem::MagnetizVector,
                       northxax::Real,body::BodySegments2D,forwardtype::String)

    ## LOOPING on segments MUST be in ANTI-CLOCKWISE order
    aclockw = checkanticlockwiseorder(body)

    if aclockw==false
        error("tmagforwardpoly2D(): vertices *not* ordered anticlockwise. Aborting.")
    end

    ##---------------------
    ## Get the angles

    ## check modules
    @assert Jind.mod >= 0.0
    @assert Jrem.mod >= 0.0

    ## `northxax` is the angle between geographic north and the positive x axis
    @assert 0.0 <= northxax <= 360.0
    Cnorth = deg2rad(northxax)

    ## check angles
    @assert -90.0 <= Jind.Ideg <= 90.0
    @assert -90.0 <= Jrem.Ideg <= 90.0

    @assert -180.0 <= Jind.Ddeg <= 180.0
    @assert -180.0 <= Jrem.Ddeg <= 180.0

    # check right forwardtype
    if forwardtype != "talwani" && forwardtype != "talwani_red" && forwardtype != "krav" && forwardtype != "wonbev"
        error("tmagforwardpoly2D(): [forwardtype] must be 'talwani' or 'talwani_red' or 'krav' or 'wonbev'")
    end

    # deg to rad
    Iind = deg2rad(Jind.Ideg)
    Dind = deg2rad(Jind.Ddeg)
    Irem = deg2rad(Jrem.Ideg)
    Drem = deg2rad(Jrem.Ddeg)

    # Calculation of Jx and Jz only for case != wonbev
    if forwardtype != "wonbev"
        Jtotx,_,Jtotz = magcomp(Jind.mod,Iind,Dind,Jrem.mod,Irem,Drem,Cnorth)
    end

    ##-------------------------------------
    ## Loop on observation points and segments
    nobs = size(xzobs,1)
    totfield = Vector{typeof(Jind.mod)}(undef,nobs)

    ## loop on observation points
    for iob=1:nobs

        xo = xzobs[iob,1]
        zo = xzobs[iob,2]

        ## loop on segments
        tsum = 0.0
        for ise=1:body.nsegm

            x1 = body.ver1[ise,1]-xo
            z1 = body.ver1[ise,2]-zo
            x2 = body.ver2[ise,1]-xo
            z2 = body.ver2[ise,2]-zo

            if forwardtype == "talwani"
                tsum += tmagtalwani(x1,z1,x2,z2,Jtotx,Jtotz,Iind,Dind,Cnorth)

            elseif forwardtype == "talwani_red"
                tsum += tmagtalwanired(x1,z1,x2,z2,Jtotx,Jtotz,Iind,Dind,Cnorth)

            elseif forwardtype == "krav"
                tsum += tmagkrav(x1,z1,x2,z2,Jtotx,Jtotz,Iind,Dind,Cnorth)

            elseif forwardtype == "wonbev"
                tsum += tmagwonbev(x1,z1,x2,z2,Jind.mod,Jrem.mod,Iind,Dind,Irem,Drem,Cnorth)

            end
         end

        # convert field from A/m to nT
        totfield[iob] = convert_H_to_B_nT( tsum )
    end

    return totfield
end

########################################################################

"""
$(TYPEDSIGNATURES)

Total magnetic field (2D) for a polygon defined by its corners. Takes into account both induced and remnant magnetization.
 Based on Talwani & Heitzler (1964), the default algorithm in Mag2Dpoly package.
"""
function tmagpoly2D(xzobs::Array{<:Real,2},Jind::MagnetizVector,Jrem::MagnetizVector,
                    northxax::Real,body::BodySegments2D)

    ## LOOPING on segments MUST be in ANTI-CLOCKWISE order
    aclockw = checkanticlockwiseorder(body)

    if aclockw==false
        error("tmagforwardpoly2D(): vertices *not* ordered anticlockwise. Aborting.")
    end

    ##---------------------
    ## Get the angles

    ## check modules
    @assert Jind.mod >= 0.0
    @assert Jrem.mod >= 0.0

    ## `northxax` is the angle between geographic north and the positive x axis
    @assert 0.0 <= northxax <= 360.0
    Cnorth = deg2rad(northxax)

    ## check angles
    @assert -90.0 <= Jind.Ideg <= 90.0
    @assert -90.0 <= Jrem.Ideg <= 90.0

    @assert -180.0 <= Jind.Ddeg <= 180.0
    @assert -180.0 <= Jrem.Ddeg <= 180.0

    # deg to rad
    Iind = deg2rad(Jind.Ideg)
    Dind = deg2rad(Jind.Ddeg)
    Irem = deg2rad(Jrem.Ideg)
    Drem = deg2rad(Jrem.Ddeg)

    # Calculation of Jx and Jz only for case != wonbev
    Jtotx,_,Jtotz = magcomp(Jind.mod,Iind,Dind,Jrem.mod,Irem,Drem,Cnorth)

    ##-------------------------------------
    ## Loop on observation points and segments
    nobs = size(xzobs,1)
    totfield = Vector{typeof(Jind.mod)}(undef,nobs)

    ## loop on observation points
    for iob=1:nobs

        xo = xzobs[iob,1]
        zo = xzobs[iob,2]

        ## loop on segments
        tsum = 0.0
        for ise=1:body.nsegm

            x1 = body.ver1[ise,1]-xo
            z1 = body.ver1[ise,2]-zo
            x2 = body.ver2[ise,1]-xo
            z2 = body.ver2[ise,2]-zo

            tsum += tmagtalwani(x1,z1,x2,z2,Jtotx,Jtotz,Iind,Dind,Cnorth)
         end

        # convert field from A/m to nT
        totfield[iob] = convert_H_to_B_nT( tsum )
    end

    return totfield
end

#######################################################################################
"""
$(TYPEDSIGNATURES)

Check whether the polygonal body has segments ordered anticlockwise.
"""
function checkanticlockwiseorder(body::BodySegments2D)
    ## https://stackoverflow.com/questions/1165647/how-to-determine-if-a-list-of-polygon-points-are-in-clockwise-order
    ## https://www.element84.com/blog/determining-the-winding-of-a-polygon-given-as-a-set-of-ordered-points
    #
    # Check direction (anti)clockwise for a reference
    #   system like the following:
    #
    #   z
    #  /\
    #  |      2
    #  |   1     3
    #  |      4
    #  |
    #  -------------> x
    #
    encarea2=0.0
    for ise=1:body.nsegm
        x1 = body.ver1[ise,1]
        z1 = body.ver1[ise,2]
        x2 = body.ver2[ise,1]
        z2 = body.ver2[ise,2]
        encarea2 += (x2-x1)*(z2+z1)
    end
    #
    # anticlockwise -> encarea2 < 0.0
    # clockwise -> encarea2 > 0.0
    if encarea2<0.0
        anticlockw=true
    else
        anticlockw=false
    end
    #
    # The reference system for the magnetic anomaly functions
    #   is reversed in z:
    #
    #  -------------> x
    #  |
    #  |      4
    #  |   1     3
    #  |      2
    #  \/
    #  z
    #
    # so, consequently, we flip the direction of
    # clockwise/anticlockwise:   !(anticlockw)
    return !anticlockw
end

#####################################################################################

"""
$(TYPEDSIGNATURES)

 Total magnetic field (2D) for a line segment. Formulas from Talwani & Heitzler (1964).
"""
function tmagtalwani(x1::Real,z1::Real,x2::Real,z2::Real,
                     Jx::Real,Jz::Real,Iind::Real,Dind::Real,C::Real)

    # Quantities for error definitions
    small = 1e4*eps(typeof(Jx))
    anglelim = 0.995*π

    #--------------
    x21 = x2-x1
    z21 = z2-z1
    s = sqrt(x21^2+z21^2)

    # Return 0 if two corners are too close
    if s < small
        return 0.0
    end

    #-----------------------

    # Get the angles
    θ1 = atan(z1,x1)
    θ2 = atan(z2,x2)

    # If z21 = 0.0 no contribution
    if z21 != 0.0
        g = -x21/z21
    else
        return 0.0
    end

    ϕ = acot(g)

    θdiff = θ2-θ1
    # In the case polygon sides cross the x axis
    if θdiff < -π
        θdiff = θdiff + 2.0*π
    elseif θdiff > π
        θdiff = θdiff - 2.0*π
    end

    #--------------------------------------

    # Error if a corner is too close to the observation point (calculation continues)
    # and the corner are slightly moved away
    if x1 < small && z1 < small
        x1 = small
        z1 = small
        @warn "A corner is too close to an observation point (calculation continues)"
    end

    if x2 < small && z2 < small
        x2 = small
        z2 = small
        @warn "A corner is too close to an observation point (calculation continues)"
    end

    r1 = sqrt(x1^2+z1^2)
    r2 = sqrt(x2^2+z2^2)

    flog = log(r2)-log(r1)

    # Error if the side is too close to the observation point (calculation continues)
    if abs(θdiff) > anglelim
        @warn "A polygon side is too close to an observation point (calculation continues)"
    end

    #-------------------------------------

    # vertical component
    V = 2.0*sin(ϕ) * (Jx*( (θdiff)*cos(ϕ) + sin(ϕ)*flog) -
                      Jz*( (θdiff)*sin(ϕ) - cos(ϕ)*flog) )

    # horizonatal component
    H = 2.0*sin(ϕ) * (Jx*( (θdiff)*sin(ϕ) - cos(ϕ)*flog)+
                      Jz*( (θdiff)*cos(ϕ) + sin(ϕ)*flog) )

    # Divided by 4π to take into account algorithm formulation in emu units
    totfield = (1.0/(4.0*π)) * (H*cos(Iind)*cos(C-Dind) + V*sin(Iind))

    return totfield
end

###########################################################################

"""
$(TYPEDSIGNATURES)

 Total magnetic field (2D) for a line segment. Formulas from Kravchinsky et al (2019) rectified by Ghirotto et al. (2021).
"""
function tmagkrav(x1::Real,z1::Real,x2::Real,z2::Real,
                  Jtotx::Real,Jtotz::Real,Iind::Real,Dind::Real,Cnorth::Real)

    # Quantities for errors definitions
    small = 1e4*eps(typeof(Jtotx))
    anglelim = 0.995*π

    #--------------
    x21 = x2-x1
    z21 = z2-z1
    tmpγ = sqrt(x21^2+z21^2)

    # Return 0 if two corners are too close
    if tmpγ < small
        return 0.0
    end

    #--------------------

    # check if den ≠ 0.0
    if tmpγ!=0.0
        γx = x21 / tmpγ
        γz = z21 / tmpγ
    else
        return 0.0
    end

    # if the segment is horizontal it provides no contribution!
    if z21==0.0
        return 0.0
    end

    #------------
    g = x21/z21

    if x1 >= g*z1
        δ = 1.0
    elseif x1 < g*z1
        δ = -1.0
    end

    #--------------------
    # Get the angles
    α1 = atan(δ*(z1+g*x1),(x1-g*z1))
    α2 = atan(δ*(z2+g*x2),(x2-g*z2))

    #In the case polygon sides cross the x axis
    αdiff = α2 - α1

    if αdiff < -π
        αdiff = αdiff + 2.0*π
    elseif αdiff > π
        αdiff = αdiff - 2.0*π
    end

    #----------------------------------

    # Error if a corner is too close to the observation point (calculation continues)
    # and the corner are slightly moved away
    if x1 < small && z1 < small
        x1 = small
        z1 = small
        @warn  "A corner is too close to an observation point (calculation continues)"
    end

    if x2 < small && z2 < small
        x2 = small
        z2 = small
        @warn "A corner is too close to an observation point (calculation continues)"
    end

    r1 = sqrt(x1^2+z1^2)
    r2 = sqrt(x2^2+z2^2)

    lor21 = log(r2)-log(r1)

    # Error if the side is too close to the observation point (calculation continues)
    if abs(αdiff) > anglelim
        @warn "A polygon side is too close to an observation point (calculation continues)"
    end

    #------------------------------------

    P = γz*γx*lor21 + δ*(γz^2)*(αdiff)
    Q = (γz^2)*lor21 - δ*γx*γz*(αdiff)

    ## horizonatl and vertical field components
    H = 1.0/(2.0*pi) * (Jtotz*Q + Jtotx*P)
    V = 1.0/(2.0*pi) * (Jtotx*Q - Jtotz*P)

    ## total field anomaly
    totfield = V*sin(Iind)+H*cos(Iind)*cos(Cnorth-Dind)

    return totfield
end

########################################################################################

"""
$(TYPEDSIGNATURES)

 Total magnetic field (2D) for a ribbon. Talwani & Heitzler (1964) modified by Kravchinsky et al. (2019).
"""
function tmagtalwanired(x1::Real,z1::Real,x2::Real,z2::Real,
                        Jx::Real,Jz::Real,Iind::Real,Dind::Real,C::Real)

    # Quantities for errors definitions
    small = 1e4*eps(typeof(Jx))
    anglelim = 0.995*π

    #--------------
    x21 = x2-x1
    z21 = z2-z1
    s = sqrt(x21^2+z21^2)

    # Return 0 if two corners are too close
    if s < small
        return 0.0
    end

    #-----------------------

    # if the segment is horizontal it provides no contribution!
    if z21 != 0.0
        g = -x21/z21
    else
        return 0.0
    end

    ϕ = acot(g)

    den1 = x1+z1*cot(ϕ)
    den2 = x2+z2*cot(ϕ)
    num1 = z1-x1*cot(ϕ)
    num2 = z2-x2*cot(ϕ)

    # Controls on signs of atan argument (abs in den1 and den2)
    #-----------------------
    if den1 < 0.0
        den1 = -den1
        δ = -1.0
        θ1 = atan(num1,den1)
    else
        δ = 1.0
        θ1 = atan(num1,den1)
    end

    if den2 < 0.0
        den2 = -den2
        θ2 = atan(num2,den2)
    else
        θ2 = atan(num2,den2)
    end
    #-----------------------

    # In the case polygon sides cross the x axis
    θdiff = θ2-θ1

    if θdiff < -π
        θdiff = θdiff + 2.0*π
    elseif θdiff > π
        θdiff = θdiff - 2.0*π
    end

    #-------------------------------------------------

    # Error if a corner is too close to the observation point (calculation continues)
    # and the corner are slightly moved away
    if x1 < small && z1 < small
        x1 = small
        z1 = small
        @warn "A corner is too close to an observation point (calculation continues)\n"
    end

    if x2 < small && z2 < small
        x2 = small
        z2 = small
        @warn "A corner is too close to an observation point (calculation continues)"
    end

    r1 = sqrt(x1^2+z1^2)
    r2 = sqrt(x2^2+z2^2)

    flog = log(r2)-log(r1)

    # Error if the side is too close to the observation point (calculation continues)
    if abs(θdiff) > anglelim
        @warn "A polygon side is too close to an observation point (calculation continues)"
    end

    #----------------------------------------------------

    # vertical component
    V = 2.0*sin(ϕ) * (Jx * (δ*(θdiff)*cos(ϕ) + sin(ϕ)*flog)-
                    Jz * (δ*(θdiff)*sin(ϕ) - cos(ϕ)*flog) )

    # horizontal component
    H = 2.0*sin(ϕ) * (Jx * (δ*(θdiff)*sin(ϕ) - cos(ϕ)*flog)+
                    Jz * (δ*(θdiff)*cos(ϕ) + sin(ϕ)*flog) )

    ## total field anomaly divided by 4π to take into account algorithm formulation in emu units
    totfield = (1.0/(4.0*π)) * (H*cos(Iind)*cos(C-Dind) + V*sin(Iind))

    return totfield
end

#######################################################################################################################

"""
$(TYPEDSIGNATURES)

 Total magnetic field (2D) for a line segment. Formulas from Won & Bevis (1987).
"""
function tmagwonbev(x1::Real,z1::Real,x2::Real,z2::Real,
                    modJind::Real,modJrem::Real,Iind::Real,Dind::Real,
                    Irem::Real,Drem::Real,C::Real)

    # β is angle among North and profle direction
    βi = Dind - C + π/2
    βr = Drem - C + π/2

    #-------------------

    # Quantities for errors definitions
    small = 1e4*eps(typeof(modJrem))
    anglelim = 0.995*π

    #----------------------
    x21 = x2-x1
    z21 = z2-z1
    R = sqrt(x21^2+z21^2)

    # Return 0 if two corners are too close
    if R < small
        return 0.0
    end

    #------------------------

    # Error if a corner is too close to the observation point (calculation continues)
    # and the corner are slightly moved away
    if x1 < small && z1 < small
        x1 = small
        z1 = small
        @warn "A corner is too close to an observation point (calculation continues)"
    end

    if x2 < small && z2 < small
        x2 = small
        z2 = small
        @warn "A corner is too close to an observation point (calculation continues)"
    end

    r1 = sqrt(x1^2+z1^2)
    r2 = sqrt(x2^2+z2^2)

    lor21 = log(r2) - log(r1)

    θ1 = atan(z1,x1)
    θ2 = atan(z2,x2)

    # In the case polygon sides cross the x axis
    if sign(z1) != sign(z2)
        test = x1*z2 - x2*z1
        if test > 0.0
            if z1 >= 0.0
                θ2 = θ2 + 2π
            end
        elseif test < 0.0
            if z2 >= 0.0
                θ1 = θ1 + 2π
            end
        else
            return 0.0
        end
    end

    # Error if the side is too close to the observation point (calculation continues)
    θdiff = θ1-θ2
    if abs(θdiff) > anglelim
        @warn "A polygon side is too close to an observation point (calculation continues)"
    end

    #-----------------------

    P = (1/R^2)*(x1*z2 - x2*z1)*(((x1*x21 - z1*z21)/(r1^2))-
                                 ((x2*x21 - z2*z21)/(r2^2)))

    Q = (1/R^2)*(x1*z2 - x2*z1)*(((x1*z21 + z1*x21)/(r1^2))-
                                 ((x2*z21 + z2*x21)/(r2^2)))

    if x21 != 0.0

        g = z21/x21
        derZz = ((x21^2)/(R^2))*((θdiff) + g*lor21) - P
        derZx = -((x21*z21)/(R^2))*((θdiff) + g*lor21) + Q
        derXz = -((x21^2)/(R^2))*(g*(θdiff) - lor21) + Q
        derXx = ((x21*z21)/(R^2))*(g*(θdiff) - lor21) + P

    else

        derZz = -P
        derZx = -((z21^2)/(R^2))*lor21 + Q
        derXz = Q
        derXx = ((z21^2)/(R^2))*(θdiff) + P

    end

    # Magnetic strenght components due to induced magnetization
    ΔHzind = 2.0*modJind*(sin(Iind)*derZz + sin(βi)*cos(Iind)*derZx)
    ΔHxind = 2.0*modJind*(sin(Iind)*derXz + sin(βi)*cos(Iind)*derXx)

    # Magnetic strenght components due to remnant magnetization
    ΔHzrem = 2.0*modJrem*(sin(Irem)*derZz + sin(βr)*cos(Irem)*derZx)
    ΔHxrem = 2.0*modJrem*(sin(Irem)*derXz + sin(βr)*cos(Irem)*derXx)

    ΔHztot = ΔHzind + ΔHzrem
    ΔHxtot = ΔHxind + ΔHxrem

    ## total field anomaly divided by 4π to take into account algorithm formulation in emu units
    ΔHtot = -(1.0/(4.0*π))*(ΔHztot*sin(Iind) + ΔHxtot*sin(βi)*cos(Iind))

    return ΔHtot
end

###################################################################################


text:
# Tests for Cost Russell DEA Models
@testset "CostRussellDEAModel" begin

    # Test using Book data
    X = [2 2; 1 4; 4 1; 4 3; 5 5; 6 1; 2 5; 1.6 8];
    Y = [1; 1; 1; 1; 1; 1; 1; 1];
    W = [1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1];

    # Cost Russell
    deacostrus = deacostrussell(X, Y, W)

    @test typeof(deacostrus) == CostRussellDEAModel

    @test nobs(deacostrus) == 8
    @test ninputs(deacostrus) == 2
    @test noutputs(deacostrus) == 1
    @test ismonetary(deacostrus) == false

    @test efficiency(deacostrus) == efficiency(deacostrus, :Economic)
    @test efficiency(deacostrus, :Economic)   ≈ [0.0; 0.5; 0.5; 0.5; 0.6; 1.5; 0.75; 1.75]
    @test efficiency(deacostrus, :Technical)  ≈ [0.0; 0.0; 0.0; 5/12; 0.6; 1/6; 0.35; 7/16]
    @test efficiency(deacostrus, :Allocative) ≈ [0.0; 0.5; 0.5; 1/12; 0.0; 4/3; 0.4; 21/16]
    @test normfactor(deacostrus) ≈ [4.0; 2.0; 2.0; 6.0; 10.0; 2.0; 4.0; 3.2]

    # With monetary = true
    deacostrusmon = deacostrussell(X, Y, W, monetary = true)

    @test ismonetary(deacostrusmon) == true

    @test efficiency(deacostrusmon) ≈ efficiency(deacostrus) .* normfactor(deacostrus)
    @test efficiency(deacostrusmon, :Technical) ≈ efficiency(deacostrus, :Technical) .* normfactor(deacostrus)
    @test efficiency(deacostrusmon, :Allocative) ≈ efficiency(deacostrus, :Allocative) .* normfactor(deacostrus)

    # With CRS (Same result in this example)
    deacostruscrs = deacostrussell(X, Y, W, rts = :CRS)

    @test efficiency(deacostruscrs, :Economic)   ≈ [0.0; 0.5; 0.5; 0.5; 0.6; 1.5; 0.75; 1.75]
    @test efficiency(deacostruscrs, :Technical)  ≈ [0.0; 0.0; 0.0; 5/12; 0.6; 1/6; 0.35; 7/16]
    @test efficiency(deacostruscrs, :Allocative) ≈ [0.0; 0.5; 0.5; 1/12; 0.0; 4/3; 0.4; 21/16]

    # Print
    show(IOBuffer(), deacostrus)

    # Test errors
    @test_throws DimensionMismatch deacostrussell([1; 2 ; 3], [4 ; 5], [1; 1; 1]) #  Different number of observations
    @test_throws DimensionMismatch deacostrussell([1; 2; 3], [4; 5; 6], [1; 2; 3; 4]) # Different number of observation in prices
    @test_throws DimensionMismatch deacostrussell([1 1; 2 2; 3 3 ], [4; 5; 6], [1 1 1; 2 2 2; 3 3 3]) # Different number of input prices and inputs
    @test_throws ArgumentError deacostrussell([1; 2; 3], [4; 5; 6], [1; 2; 3], rts = :Error) # Invalid returns to scale

end


text:
@testset "Creating an LP" begin
    model = GLPK.Optimizer()
    x̂ = generate_lp(model,10,5)

    MOI.set(model, MOI.Silent(), true)
    MOI.optimize!(model)

    @test MOI.get(model, MOI.TerminationStatus()) in [MOI.LOCALLY_SOLVED, MOI.OPTIMAL]
end

@testset "Creating a convex QP" begin
    model = MOI.instantiate(OSQP.Optimizer, with_bridge_type=Float64)
    x̂ = generate_qp(model,10,5,5)

    MOI.set(model, MOI.Silent(), true)
    MOI.optimize!(model)

    @test MOI.get(model, MOI.TerminationStatus()) in [MOI.LOCALLY_SOLVED, MOI.OPTIMAL]
end


text:
using Test
using Scruff
using Scruff.Utils
using Scruff.RTUtils
using Scruff.Models
using Scruff.SFuncs
using Scruff.Operators
import Scruff.Algorithms: VE, ve, infer, probability, greedy_order, unconnected_neighbors, cost, copy_graph, eliminate

@testset "VE" begin

    @testset "vegraph" begin
        g = Graph()
        add_node!(g, 1, 2)
        add_node!(g, 2, 3)
        add_node!(g, 3, 2)
        add_node!(g, 4, 2)
        add_node!(g, 5, 5)
        add_node!(g, 6, 1)
        add_undirected!(g, 1, 3)
        add_undirected!(g, 2, 3)
        add_undirected!(g, 3, 4)
        add_undirected!(g, 3, 5)
        add_undirected!(g, 1, 4)
        @testset "Construction" begin
            @test g.nodes == [1, 2, 3, 4, 5, 6]
            @test g.sizes == Dict(1 => 2, 2 => 3, 3 => 2, 4 => 2, 5 => 5, 6 => 1)
            @test g.edges == Dict(1 => [3, 4], 2 => [3], 3 => [1, 2, 4, 5],
                                    4 => [3, 1], 5 => [3], 6 => [])
        end

        @testset "Unconnected neighbors" begin
            @test unconnected_neighbors(g, 3) == [(1,2), (1,5), (2,4), (2,5), (4,5)]
            @test unconnected_neighbors(g, 6) == []
            @test cost(g, 3) == 5
            @test cost(g, 6) == 0
        end

        @testset "Elimination" begin
            h = copy_graph(g)
            eliminate(h, 3)
            @test h.nodes == [1, 2, 4, 5, 6]
            @test h.sizes == Dict(1 => 2, 2 => 3, 4 => 2, 5 => 5, 6 => 1)
            @test h.edges == Dict(1 => [4, 2, 5], 2 => [1, 4, 5],
                                    4 => [1, 2, 5], 5 => [1, 2, 4], 6 => [])
        end

        @testset "Greedy elimination order" begin
            @testset "With all variables eliminated" begin
                ord = greedy_order(g)
                @test Set(ord) == Set([1,2,3,4,5,6])
                @test length(ord) == 6
                # Only 3 has unconnected neighbors.
                # It cannot be eliminated before at least 3 of its neighbors
                # have been eliminated.
                inds = indexin([1,2,4,5,3], ord)
                count = 0
                if inds[1] < inds[4] count += 1 end
                if inds[2] < inds[4] count += 1 end
                if inds[3] < inds[4] count += 1 end
                if inds[4] < inds[4] count += 1 end
                @test count >= 3
            end

            @testset "With uneliminated variables" begin
                ord = greedy_order(g, [5, 6])
                @test Set(ord) == Set([1,2,3,4])
                @test length(ord) == 4
                # 3 must be eliminated last because all the others
                # are a disconnected neighbor from 5
                @test ord[4] == 3
            end
        end
    end

    @testset "range" begin
        dn1 = Cat([1,2], [0.1, 0.9])
        dv1 = dn1()(:dv1)
        dn2 = Cat([1,2,3], [0.2, 0.3, 0.5])
        dv2 = dn2()(:dv2)
        dn3 = DiscreteCPT([1,2],
            Dict((1,1) => [0.3, 0.7], (2,1) => [0.6, 0.4],
                 (1,2) => [0.4, 0.6], (2,2) => [0.7, 0.3],
                 (1,3) => [0.5, 0.5], (2,3) => [0.8, 0.2]))
        dv3 = dn3()(:dv3)
        cn1 = Normal(-0.1, 1.0)
        cv1 = cn1()(:cv1)
        cn2 = CLG(Dict((1,) => ((1.5,), -0.3, 0.5), (2,) => ((0.7,), 0.4, 0.5)))
        cv2 = cn2()(:cv2)

        network = InstantNetwork(Variable[dv1,dv2,dv3,cv1,cv2], VariableGraph(dv3=>[dv1,dv2],cv2=>[dv3,cv1]))
        runtime = Runtime(network)
        ensure_all!(runtime, 0)
        order = topsort(get_initial_graph(network))

        iv = Vector{Int}()
        fv = Vector{Float64}()
        dr1 = Operators.support(dn1, (), 10, iv)
        dr2 = Operators.support(dn2, (), 10, iv)
        dr3 = Operators.support(dn3, (dr1, dr2), 10, iv)
        cr1 = Operators.support(cn1, (), 10, fv)
        cr2 = Operators.support(cn2, (dr3, cr1), 10, fv)

        dx1 = Operators.support(dn1, (), 20, dr1)
        dx2 = Operators.support(dn2, (), 20, dr2)
        dx3 = Operators.support(dn3, (dx1, dx2), 20, dr3)
        cx1 = Operators.support(cn1, (), 20, cr1)
        cx2 = Operators.support(cn2, (dx3, cx1), 20, cr2)

        @testset "Setting initial ranges" begin
            set_ranges!(runtime, Dict{Symbol, Score}(), 10)
            @test get_range(runtime, dv1) == dr1
            @test get_range(runtime, dv2) == dr2
            @test get_range(runtime, dv3) == dr3
            @test get_range(runtime, cv1) == cr1
            @test get_range(runtime, cv2) == cr2
        end

        @testset "Setting expanded ranges" begin
            set_ranges!(runtime, Dict{Symbol, Score}(), 20)
            @test get_range(runtime, dv1) == dx1
            @test get_range(runtime, dv2) == dx2
            @test get_range(runtime, dv3) == dx3
            @test get_range(runtime, cv1) == cx1
            @test get_range(runtime, cv2) == cx2
        end

        @testset "Ranges from previous instance" begin
            ensure_all!(runtime, 2)
            @test get_range(runtime, dv1) == dx1
        end

    end

    @testset "ve" begin
        dn1 = Cat([1,2], [0.1, 0.9])
        i11 = indexin(1, dn1.range)[1]
        i12 = indexin(2, dn1.range)[1]
        dv1 = dn1()(:dv1)
        dn2 = Cat([1,2,3], [0.2, 0.3, 0.5])
        i21 = indexin(1, dn2.range)[1]
        i22 = indexin(2, dn2.range)[1]
        i23 = indexin(3, dn2.range)[1]
        dv2 = dn2()(:dv2)
        dn3 = DiscreteCPT([1,2], Dict((1,1) => [0.3, 0.7], (2,1) => [0.4, 0.6],
                                (1,2) => [0.5, 0.5], (2,2) => [0.6, 0.4],
                                (1,3) => [0.7, 0.3], (2,3) => [0.8, 0.2]))
        i31 = indexin(1, dn3.sfs[1,1].range)[1]
        i32 = indexin(2, dn3.sfs[1,1].range)[1]
        dv3 = dn3()(:dv3)
        cn1 = Normal(-0.1, 1.0)
        cv1 = cn1()(:cv1)
        cn2 = CLG(Dict((1,) => ((1.5,), -0.3, 0.5), (2,) => ((0.7,), 0.4, 0.5)))
        cv2 = cn2()(:cv2)
        net1 = InstantNetwork(Variable[dv1,dv2,dv3], VariableGraph(dv3=>[dv1,dv2]))
        ord1 = topsort(get_initial_graph(net1))

        @testset "A discrete network" begin

            @testset "With one query variable and bounds" begin
                runtime = Runtime(net1)
                ensure_all!(runtime, 0)
                set_ranges!(runtime, Dict{Symbol, Score}(), 10)
                ((l,u), ids) = ve(runtime, ord1, [dv3]; bounds = true)
                pa = 0.1 * 0.2 * 0.3 + 0.1 * 0.3 * 0.5 + 0.1 * 0.5 * 0.7 +
                        0.9 * 0.2 * 0.4 + 0.9 * 0.3 * 0.6 + 0.9 * 0.5 * 0.8
                pb = 1 - pa
                @test length(l.keys) == 1
                @test l.keys[1] == ids[dv3]
                @test l.dims == (2,)
                @test length(l.entries) == 2
                @test isapprox(l.entries[i31], pa, atol = 0.0000001)
                @test isapprox(l.entries[i32], pb, atol = 0.0000001)
                @test length(u.keys) == 1
                @test u.keys[1] == ids[dv3]
                @test u.dims == (2,)
                @test length(u.entries) == 2
                @test isapprox(u.entries[i31], pa, atol = 0.0000001)
                @test isapprox(u.entries[i32], pb, atol = 0.0000001)
            end

            @testset "With one query variable and no bounds" begin
                runtime = Runtime(net1)
                ensure_all!(runtime, 0)
                set_ranges!(runtime, Dict{Symbol, Score}(), 10)
                (l,ids) = ve(runtime, ord1, [dv3]; bounds = false)
                pa = 0.1 * 0.2 * 0.3 + 0.1 * 0.3 * 0.5 + 0.1 * 0.5 * 0.7 +
                        0.9 * 0.2 * 0.4 + 0.9 * 0.3 * 0.6 + 0.9 * 0.5 * 0.8
                pb = 1 - pa
                @test length(l.keys) == 1
                @test l.keys[1] == ids[dv3]
                @test l.dims == (2,)
                @test length(l.entries) == 2
                @test isapprox(l.entries[i31], pa, atol = 0.0000001)
                @test isapprox(l.entries[i32], pb, atol = 0.0000001)
            end

            @testset "With disconnected variable" begin
                x = Cat([1,2], [0.5, 0.5])()(:x)
                y = Cat([1,2], [0.2, 0.8])()(:y)
                net = InstantNetwork(Variable[x,y], VariableGraph())
                run = Runtime(net)
                ensure_all!(run)
                ord = topsort(get_initial_graph(net))
                set_ranges!(run, Dict{Symbol, Score}(), 2)
                (l,ids) = ve(run, ord, [x]; bounds = false)
                @test l.entries == [0.5, 0.5]
            end

            @testset "With two query variables" begin
                runtime = Runtime(net1)
                ensure_all!(runtime, 0)
                set_ranges!(runtime, Dict{Symbol, Score}(), 10)
                (l,ids) = ve(runtime, ord1, [dv3, dv1]; bounds = false)
                ppa = 0.1 * 0.2 * 0.3 + 0.1 * 0.3 * 0.5 + 0.1 * 0.5 * 0.7
                ppb = 0.1 * 0.2 * 0.7 + 0.1 * 0.3 * 0.5 + 0.1 * 0.5 * 0.3
                pqa = 0.9 * 0.2 * 0.4 + 0.9 * 0.3 * 0.6 + 0.9 * 0.5 * 0.8
                pqb = 0.9 * 0.2 * 0.6 + 0.9 * 0.3 * 0.4 + 0.9 * 0.5 * 0.2
                @test l.dims == (2,2)
                @test length(l.keys) == 2
                @test length(l.entries) == 4
                k1 = l.keys[1]
                k2 = l.keys[2]
                @test k1 == ids[dv1] && k2 == ids[dv3] || k1 == ids[dv3] && k2 == ids[dv1]
                if k1 == ids[dv1]
                    @test isapprox(l.entries[(i11-1)*2 + i31], ppa, atol = 0.000001)
                    @test isapprox(l.entries[(i11-1)*2 + i32], ppb, atol = 0.000001)
                    @test isapprox(l.entries[(i12-1)*2 + i31], pqa, atol = 0.000001)
                    @test isapprox(l.entries[(i12-1)*2 + i32], pqb, atol = 0.000001)
                else
                    @test isapprox(l.entries[(i31-1)*2 + i11], ppa, atol = 0.000001)
                    @test isapprox(l.entries[(i31-1)*2 + i12], pqa, atol = 0.000001)
                    @test isapprox(l.entries[(i32-1)*2 + i11], ppb, atol = 0.000001)
                    @test isapprox(l.entries[(i32-1)*2 + i12], pqb, atol = 0.000001)
                end

            end

            @testset "with hard evidence" begin
                runtime = Runtime(net1)
                ensure_all!(runtime)
                set_ranges!(runtime, Dict{Symbol, Score}(), 10)
                inst1 = current_instance(runtime, dv1)
                post_evidence!(runtime, inst1, HardScore(1))
                (l, ids) = ve(runtime, ord1, [dv3]; bounds = false)
                p31 = 0.1 * (0.2 * 0.3 + 0.3 * 0.5 + 0.5 * 0.7)
                p32 = 0.1 * (0.2 * 0.7 + 0.3 * 0.5 + 0.5 * 0.3)
                @test isapprox(l.entries[i31], p31, atol = 0.000001)
                @test isapprox(l.entries[i32], p32, atol = 0.000001)
            end

            @testset "with soft evidence" begin
                runtime = Runtime(net1)
                ensure_all!(runtime)
                set_ranges!(runtime, Dict{Symbol, Score}(), 10)
                inst1 = current_instance(runtime, dv1)
                post_evidence!(runtime, inst1, SoftScore([1,2], [3.0, 5.0]))
                (l, ids) = ve(runtime, ord1, [dv3]; bounds = false)
                p31 = 0.1 * 3.0 * (0.2 * 0.3 + 0.3 * 0.5 + 0.5 * 0.7) +
                        0.9 * 5.0 * (0.2 * 0.4 + 0.3 * 0.6 + 0.5 * 0.8)
                p32 = 0.1 * 3.0 * (0.2 * 0.7 + 0.3 * 0.5 + 0.5 * 0.3) +
                        0.9 * 5.0 * (0.2 * 0.6 + 0.3 * 0.4 + 0.5 * 0.2)
                @test isapprox(l.entries[i31], p31, atol = 0.000001)
                @test isapprox(l.entries[i32], p32, atol = 0.000001)
            end

            @testset "with separable models" begin
                sf1 = Cat([1,2], [0.1, 0.9])
                z1 = sf1()(:z1)
                sf2 = DiscreteCPT([1,2], Dict((1,) => [0.2, 0.8], (2,) => [0.3, 0.7]))
                z2 = sf2()(:z2)
                sf3 = DiscreteCPT([1,2], Dict((1,) => [0.4, 0.6], (2,) => [0.5, 0.5]))
                z3 = sf3()(:z3)
                cpt1 = Dict((1,) => [0.6, 0.4], (2,) => [0.7, 0.3])
                cpt2 = Dict((1,) => [0.8, 0.2], (2,) => [0.9, 0.1])
                cpts :: SepCPTs = [cpt1, cpt2]
                sf4 = Separable([1,2], [0.75, 0.25], cpts)
                z4 = sf4()(:z4)
                net = InstantNetwork(Variable[z1,z2,z3,z4], VariableGraph(z2=>[z1],z3=>[z1],z4=>[z2,z3]))
                ord = topsort(get_initial_graph(net))

                aceg = 0.1 * 0.2 * 0.4 * (0.75 * 0.6 + 0.25 * 0.8)
                aceh = 0.1 * 0.2 * 0.4 * (0.75 * 0.4 + 0.25 * 0.2)
                acfg = 0.1 * 0.2 * 0.6 * (0.75 * 0.6 + 0.25 * 0.9)
                acfh = 0.1 * 0.2 * 0.6 * (0.75 * 0.4 + 0.25 * 0.1)
                adeg = 0.1 * 0.8 * 0.4 * (0.75 * 0.7 + 0.25 * 0.8)
                adeh = 0.1 * 0.8 * 0.4 * (0.75 * 0.3 + 0.25 * 0.2)
                adfg = 0.1 * 0.8 * 0.6 * (0.75 * 0.7 + 0.25 * 0.9)
                adfh = 0.1 * 0.8 * 0.6 * (0.75 * 0.3 + 0.25 * 0.1)
                bceg = 0.9 * 0.3 * 0.5 * (0.75 * 0.6 + 0.25 * 0.8)
                bceh = 0.9 * 0.3 * 0.5 * (0.75 * 0.4 + 0.25 * 0.2)
                bcfg = 0.9 * 0.3 * 0.5 * (0.75 * 0.6 + 0.25 * 0.9)
                bcfh = 0.9 * 0.3 * 0.5 * (0.75 * 0.4 + 0.25 * 0.1)
                bdeg = 0.9 * 0.7 * 0.5 * (0.75 * 0.7 + 0.25 * 0.8)
                bdeh = 0.9 * 0.7 * 0.5 * (0.75 * 0.3 + 0.25 * 0.2)
                bdfg = 0.9 * 0.7 * 0.5 * (0.75 * 0.7 + 0.25 * 0.9)
                bdfh = 0.9 * 0.7 * 0.5 * (0.75 * 0.3 + 0.25 * 0.1)

                @testset "without evidence" begin
                    run = Runtime(net)
                    ensure_all!(run)
                    set_ranges!(run, Dict{Symbol, Score}(), 10)
                    (l, ids) = ve(run, ord, [z4]; bounds = false)
                    g = aceg + acfg + adeg + adfg + bceg + bcfg + bdeg + bdfg
                    h = aceh + acfh + adeh + adfh + bceh + bcfh + bdeh + bdfh
                    es = normalize(l.entries)
                    r4 = get_range(run, current_instance(run, z4))
                    i1 = indexin(1, r4)[1]
                    i2 = indexin(2, r4)[1]
                    @test isapprox(es[i1], g)
                    @test isapprox(es[i2], h)
                end

                @testset "with evidence" begin
                    run = Runtime(net)
                    ensure_all!(run)
                    set_ranges!(run, Dict{Symbol, Score}(), 10)
                    inst4 = current_instance(run, z4)
                    post_evidence!(run, inst4, HardScore(2))
                    (l, ids) = ve(run, ord, [z1]; bounds = false)
                    ah = aceh + acfh + adeh + adfh
                    bh = bceh + bcfh + bdeh + bdfh
                    es = normalize(l.entries)
                    r1 = get_range(run, current_instance(run, z1))
                    i1 = indexin(1, r1)[1]
                    i2 = indexin(2, r1)[1]
                    @test isapprox(es[i1], ah / (ah + bh))
                    @test isapprox(es[i2], bh / (ah + bh))
                end

            end

        end

        @testset "Using the VE instant algorithm" begin
            @testset "Basic" begin
                v1 = Cat([:a,:b], [0.1, 0.9])()(:v1)
                v2 = DiscreteCPT([1,2], Dict((:a,) => [0.2, 0.8], (:b,) => [0.3, 0.7]))()(:v2)
                net = InstantNetwork(Variable[v1,v2], VariableGraph(v2 => [v1]))
                runtime = Runtime(net)
                default_initializer(runtime)
                alg = VE([v1,v2])
                infer(alg, runtime)
                i1 = current_instance(runtime, v1)
                i2 = current_instance(runtime, v2)
                @test isapprox(probability(alg, runtime, i1, :a), 0.1)
                @test isapprox(probability(alg, runtime, i1, :b), 0.9)
                @test isapprox(probability(alg, runtime, i2, 1), 0.1 * 0.2 + 0.9 * 0.3)
                @test isapprox(probability(alg, runtime, i2, 2), 0.1 * 0.8 + 0.9 * 0.7)
            end

            @testset "With placeholder" begin
                p1 = Placeholder{Symbol}(:p1)
                v2 = DiscreteCPT([1,2], Dict((:a,) => [0.2, 0.8], (:b,) => [0.3, 0.7]))()(:v2)
                net = InstantNetwork(Variable[v2], VariableGraph(v2 => [p1]), Placeholder[p1])
                runtime = Runtime(net)
                default_initializer(runtime, 10, Dict(p1.name => Cat([:a,:b], [0.1, 0.9])))
                alg = VE([v2])
                infer(alg, runtime)
                i2 = current_instance(runtime, v2)
                @test isapprox(probability(alg, runtime, i2, 1), 0.1 * 0.2 + 0.9 * 0.3)
                @test isapprox(probability(alg, runtime, i2, 2), 0.1 * 0.8 + 0.9 * 0.7)
            end

            @testset "With evidence" begin
                v1 = Cat([:a,:b], [0.1, 0.9])()(:v1)
                v2 = DiscreteCPT([1,2], Dict((:a,) => [0.2, 0.8], (:b,) => [0.3, 0.7]))()(:v2)
                net = InstantNetwork(Variable[v1,v2], VariableGraph(v2 => [v1]))
                runtime = Runtime(net)
                default_initializer(runtime)
                alg = VE([v1,v2])
                infer(alg, runtime, Dict{Symbol, Score}(:v2 => HardScore(2)))
                i1 = current_instance(runtime, v1)
                i2 = current_instance(runtime, v2)
                p1 = 0.1 * 0.8
                p2 = 0.9 * 0.7
                z = p1 + p2
                @test isapprox(probability(alg, runtime, i1, :a), p1 / z)
                @test isapprox(probability(alg, runtime, i1, :b), p2 / z)
                @test isapprox(probability(alg, runtime, i2, 1), 0.0)
                @test isapprox(probability(alg, runtime, i2, 2), 1.0)
            end
        end

    end

end


text:
# This proposal samples from prior and updates exogenous variables
# Using inverse transform

function stdproposal(qω, x::Member{<:PrimDist}, ω)
  if x in keys(ω.data)
    nothing
  else
    x_ = rand(qω, x.class)
    @show x => x_
  end
end

stdproposal(qω, x, ω) = nothing

function stdproposal(qω, x::Member{<:Distribution}, ω)
  @show x
  @show ω.data[x]
  if x in keys(ω)
    @show "here"
    @assert false
  else
    @show "not here"
    @assert false
  end
end


text:
using MLStyle

module Linq
    map(arr, f) = Base.map(f, arr)

    filter(arr, f) = Base.filter(f, arr)

    collect(arr) = Base.collect(arr)

    collect(arr, ::Type{T}) where T = Base.collect(T, arr)

    flat_map(arr, f) = Base.vcat(Base.map(f, arr)...)

    skip(arr) = Base.view(arr, 2:Base.length(arr))

    skip(arr, n) = Base.view(arr, n:Base.length(arr))

    len(arr) = Base.length(arr)

    drop(arr, n) = Base.view(arr, 1:(Base.length(arr) - n))

    sum(arr, f) = Base.sum(Base.map(f, arr))

    group_by(arr, f) = begin
        result = OrderedDict()
        for elt in arr
            push!(get!(result, f(elt)) do
                    []
                  end,
                  elt)

        end
        result
    end

    group_by(arr) = begin
        result = OrderedDict()
        for elt in arr
            push!(get!(result, elt) do
                    []
                  end,
                  elt)

        end
        result
    end

    any(arr, f) = Base.any(f, arr)
    any(arr) = Base.any(arr)

    all(arr, f) = Base.all(f, arr)
    all(arr) = Base.all(arr)

    enum(arr) = enumerate(arr)

    foldl(arr, f) = Base.foldl(f, arr)
    foldl(arr, f, init) = Base.foldl(f, arr, init=init)
    foldr(arr, f) = Base.foldr(f, arr)
    foldr(arr, f, init) = Base.foldr(f, arr, init=init)
    sort(arr) = Base.sort(arr)
    sort(arr, f) = Base.sort(arr, by=f)

    export dispatch
    """
    You can implement your own extension methods with this.

    e.g.
    import Linq: dispatch

    function dispatch(Val{:expected}, subject::Vector{T}) where T
        view(subject, 1:2:length(subject))
    end

    @assert [1, 3] == @linq [1, 2, 3, 4].expected
    """
    function dispatch
    end

end

function linq(expr)
    @match expr begin
        :($subject.$method($(args...))) =>
            let subject = linq(subject)
                if isdefined(Linq, method)
                    let method = getfield(Linq, method)
                        :($method($subject, $(args...)))
                    end
                else
                    let method = Val(method)
                        :($(Linq.dispatch)($method, $subject, $(args...)))
                    end
                end
            end
        :($subject.$method) =>
            let subject = linq(subject)
                if isdefined(Linq, method)
                    let method = getfield(Linq, method)
                        :($method($subject))
                    end
                else
                    let method = Val(method)
                        :($(Linq.dispatch)($method, $subject))
                    end
                end
            end
        _ => expr
    end
end

macro linq(expr)
    esc(linq(expr))
end

import .Linq: dispatch

dispatch(::Val{:str}, arr :: Vector{T}) where T = string(arr)

dispatch(::Val{:println}, s :: String) where T = println(s)

@linq [1, 2, 3].map(x -> 2x).str.println

max_stars_repo_name: thautwarm/MLStyle-Playground

max_issues_repo_name: cherichy/MLStyle-Playground

max_forks_repo_name: cherichy/MLStyle-Playground


text:
# This code is part of QuantumCircuits.
#
# (C) Copyright Rafał Pracht 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

module Execute

using PyCall

using QuantumCircuits.QCircuits.Math
using QuantumCircuits.QCircuits.QBase
using QuantumCircuits.QCircuits.QBase: QuantumDevice, tomatrix
using QuantumCircuits.QCircuits.Circuit
using QuantumCircuits.QCircuits.Circuit: toQiskit
using QuantumCircuits.QCircuits.Qiskit
using QuantumCircuits.QCircuits.Qiskit: qiskit
using QuantumCircuits.QCircuits.Gates: ParamT
using QuantumCircuits.QCircuits.Registers

export execute, QuantumSimulator, QiskitQuantum, @ket_str, @bra_str,
       loss_expected_zero_state, qderivative, qexecute

"Quantum Simulator"
struct QuantumSimulator <: QuantumDevice end

"Qiskit Quantum Device"
struct QiskitQuantum <: QuantumDevice
    backend::PyObject
end
QiskitQuantum() = QiskitQuantum(qiskit.providers.aer.QasmSimulator())

######################################################

# Initial stata
const _0 = [1+0im; 0]
const _1 = [0im; 1]

const qubit_map = Dict('0' => _0, '1' => _1)

#const _00 = kron(_0, _0)

function str2state(s)
    @assert all([c in ('0', '1') for c in s]) "Only 0 and 1 are allowed."

    # Reverse the bit order
    s = reverse(s)
    state = qubit_map[s[1]]

    if length(s) == 1
        return state
    end

    for d in s[2:end]
        state = kron(qubit_map[d], state)
    end

    return state
end

"Macro proces ket binary vector to quantum state"
macro ket_str(s)
    return str2state(s)
end

"Macro proces ket binary vector to quantum state"
macro bra_str(s)
    return dagger(str2state(s))
end

"Method convert state to probability"
state2probability(s) = abs.(s) .^ 2

"Function execute the quantum circuit on simulator"
function execute(device::QuantumSimulator, qc::QCircuit, params=nothing; prev_matrix=nothing)
    # The initial, all zeros vector
    s = vcat([1], zeros(2^qc.qubits-1))

    m = tomatrix(qc, params)
    if prev_matrix != nothing
        s = prev_matrix * s
    end

    p = state2probability(m * s)
    return qc.measures_matrix * p
end

function extractProbability(counts::Dict, qubits::Integer, shots::Integer)
    p = zeros(2^qubits)

    for (k, n) in counts
        p[parse(Int, k; base=2) + 1] = n/shots
    end

    return p
end

function execute(device::QiskitQuantum, qc::QiskitCircuit; shots=8192, debug=true)
    job = qiskit.execute(qc.qc, device.backend, shots=shots)
    debug && println("Job ID: $(job.job_id()) run single circuit.")
    #debug && qiskit.tools.monitor.job_monitor(job)

    res = job.result().get_counts()

    return extractProbability(res, length(qc.cbits), shots)
end

function execute(device::QuantumSimulator, circuits::Vector{QCircuit})
    results = Vector{Float64}[]
    for qc in circuits
        push!(results, execute(device, qc))
    end

    return results
end

function execute(device::QiskitQuantum, circuits::Vector{QiskitCircuit}; shots=8192, debug=true)
    jobs = [qc.qc for qc in circuits]
    job = qiskit.execute(jobs, device.backend, shots=shots)
    debug && println("Job ID: $(job.job_id()) run $(length(jobs)) circuits.")
    debug && qiskit.tools.monitor.job_monitor(job)

    # Get results
    job_results = job.result().get_counts()

    results = Vector{Float64}[]
    for (res, qc) in zip(job_results, circuits)
        push!(results, extractProbability(res, length(qc.cbits), shots))
    end

    return results
end

######################################################

function generate_mesuere_circuit(circuit::QCircuit, i)
    qr = QuantumRegister(length(circuit.vqubits), "q")
    cr = ClassicalRegister(length(circuit.vcbits))
    qc = QCircuit(qr, cr)

    measureQubits = [getid(q) for (q, b) in circuit.measures]

    for (i, v) in enumerate(reverse(string(i, base=2, pad=length(circuit.vcbits))))
        if v == '1'
            qc.x(measureQubits[i])
        end
    end

    for (q, c) in circuit.measures
        qc.measure(qr[getid(q)], cr[getid(c)])
    end
    return qc
end

function generate_mesuere_circuits(circuit::QCircuit)
    circuits = QiskitCircuit[]
    for i in 0:(2^length(circuit.vcbits)-1)
        push!(circuits, toQiskit(generate_mesuere_circuit(circuit, i)))
    end

    return circuits
end

function correctMeasures(correctMatrix, state)
    nstate = correctMatrix * state
    nstate = min.(max.(nstate, 0.0), 1.0)

    return nstate / sum(nstate)
end

######################################################

"Loss method to check if final state is zero."
function loss_expected_zero_state(state)
    return -log(real(state[1])+1e-32) + sum([real(v)^2 for v in state[2:end]])
end

######################################################

"Set the parameters to the cicquit and convert to qiskit"
function setAndConvert(qc, params::Vector{T}) where T <: ParamT
    setparameters!(qc, params)
    return toQiskit(qc)
end

"Function return the shift"
function getShift(n, idx)
    shift = zeros(ParamT, n)
    shift[idx] = π/2

    return shift
end

"Caclulate the jacobian on device."
function qjacobian(backend, qc, params::Vector{T}, runN=1, corrMes=true) where T <: ParamT
    @assert (length(params) * runN * 2 + runN) < 300 "To mamny run."
    circuits = QiskitCircuit[]

    if corrMes
        mc = generate_mesuere_circuits(qc)
        mcN = length(mc)
        append!(circuits, mc)
    end

    # function call
    rqc = setAndConvert(qc, params)
    for i in 1:runN
        push!(circuits, rqc)
    end

    # shift rule for all params
    n = length(params)
    for i in 1:n
        shift = getShift(n, i)
        slqc = setAndConvert(qc, params + shift)
        srqc = setAndConvert(qc, params - shift)
        for i in 1:runN
            push!(circuits, slqc)
        end
        for i in 1:runN
            push!(circuits, srqc)
        end
    end

    results = execute(backend, circuits)

    # correct the measures
    if corrMes
        corrMes = hcat(results[1:mcN]...)
        correctMeas = inv(corrMes)
        results = [correctMeasures(correctMeas, r) for r in results[(mcN+1):end]]
        #results = results[(mcN+1):end]
    end

    f_results = results[1:runN]
    df_results_tmp = results[(runN+1):end]

    if runN > 1
        df_results = Vector{Float64}[]
        df = zeros(ParamT, length(df_results_tmp[1]))
        for i in 1:length(df_results_tmp)
            df = df + df_results_tmp[i]
            if i % runN == 0
                push!(df_results, df / runN)
                df = zeros(ParamT, length(df_results_tmp[1]))
            end
        end
        push!(df_results, df / runN)
    else
        df_results = df_results_tmp
    end

    # proces result
    mat = (df_results[1] - df_results[2]) / 2
    for i in 2:n
        idx = 2 * i
        m = (df_results[idx - 1] - df_results[idx]) / 2
        mat = hcat(mat, m)
    end

    return sum(f_results)/runN, mat
end

"Caclulate loss value and it derivateive"
function qderivative(backend, qc, loss, params::Vector{T}, runN=1) where T <: ParamT
    y, jak = qjacobian(backend, qc, params, runN)
    der = transpose(jak) * loss'(y)
    return loss(y), der
end

"Caclulate the the value."
function qexecute(backend, qc, loss, params::Vector{T}, runN=1, corrMes=true) where T <: ParamT

    # set the parameters
    setparameters!(qc, params)
    # convert to qiskit
    qqc = toQiskit(qc)

    if corrMes
        circuits = generate_mesuere_circuits(qc)
        mcN = length(circuits)

        # add the orginal one
        push!(circuits, qqc)

        results = execute(backend, circuits)

        # correct the measures
        corrMes = hcat(results[1:mcN]...)
        correctMeas = inv(corrMes)

        ϕ = correctMeasures(correctMeas, results[end])
    else
        ϕ = execute(backend, qqc)
    end

    return loss(ϕ)
end

end  # module Execute