text:
# Autogenerated wrapper script for normaliz_jll for i686-linux-gnu-cxx03
export libnormaliz, normaliz

using GMP_jll
using MPFR_jll
using FLINT_jll
using nauty_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("normaliz")
JLLWrappers.@declare_library_product(libnormaliz, "libnormaliz.so.3")
JLLWrappers.@declare_executable_product(normaliz)
function __init__()
    JLLWrappers.@generate_init_header(GMP_jll, MPFR_jll, FLINT_jll, nauty_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libnormaliz,
        "lib/libnormaliz.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        normaliz,
        "bin/normaliz",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#! format: off
struct RangeReserve <: AbstractReservesFormulation end
struct StepwiseCostReserve <: AbstractReservesFormulation end
struct RampReserve <: AbstractReservesFormulation end

############################### Reserve Variables #########################################

get_variable_sign(_, ::Type{<:PSY.Reserve}, ::AbstractReservesFormulation) = NaN
############################### ActiveServiceVariable, Reserve #########################################

get_variable_binary(::ActiveServiceVariable, ::Type{<:PSY.Reserve}, ::AbstractReservesFormulation) = false
get_variable_upper_bound(::ActiveServiceVariable, ::PSY.Reserve, d::PSY.Component, _) = PSY.get_max_active_power(d)
get_variable_upper_bound(::ActiveServiceVariable, ::PSY.Reserve, d::PSY.Storage, _) =  PSY.get_output_active_power_limits(d).max
get_variable_lower_bound(::ActiveServiceVariable, ::PSY.Reserve, ::PSY.Component, _) = 0.0

############################### ServiceRequirementVariable, ReserveDemandCurve ################################

get_variable_binary(::ServiceRequirementVariable, ::Type{<:PSY.ReserveDemandCurve}, ::AbstractReservesFormulation) = false
get_variable_upper_bound(::ServiceRequirementVariable, ::PSY.ReserveDemandCurve, d::PSY.Component, ::AbstractReservesFormulation) = PSY.get_max_active_power(d)
get_variable_lower_bound(::ServiceRequirementVariable, ::PSY.ReserveDemandCurve, ::PSY.Component, ::AbstractReservesFormulation) = 0.0

#! format: on
################################## Reserve Requirement Constraint ##########################
function service_requirement_constraint!(
    optimization_container::OptimizationContainer,
    service::SR,
    ::ServiceModel{SR, T},
) where {SR <: PSY.Reserve, T <: AbstractReservesFormulation}
    parameters = model_has_parameters(optimization_container)
    initial_time = model_initial_time(optimization_container)
    @debug initial_time
    time_steps = model_time_steps(optimization_container)
    name = PSY.get_name(service)
    constraint =
        get_constraint(optimization_container, make_constraint_name(REQUIREMENT, SR))
    reserve_variable = get_variable(optimization_container, name, SR)
    use_slacks = get_services_slack_variables(optimization_container.settings)

    ts_vector = get_time_series(optimization_container, service, "requirement")

    use_slacks && (slack_vars = reserve_slacks(optimization_container, name))

    requirement = PSY.get_requirement(service)
    if parameters
        container = get_parameter_container(
            optimization_container,
            UpdateRef{SR}(SERVICE_REQUIREMENT, "requirement"),
        )
        param = get_parameter_array(container)
        multiplier = get_multiplier_array(container)
        for t in time_steps
            param[name, t] = add_parameter(optimization_container.JuMPmodel, ts_vector[t])
            multiplier[name, t] = requirement
            if use_slacks
                resource_expression = sum(reserve_variable[:, t]) + slack_vars[t]
            else
                resource_expression = sum(reserve_variable[:, t])
            end
            constraint[name, t] = JuMP.@constraint(
                optimization_container.JuMPmodel,
                resource_expression >= param[name, t] * requirement
            )
        end
    else
        for t in time_steps
            constraint[name, t] = JuMP.@constraint(
                optimization_container.JuMPmodel,
                sum(reserve_variable[:, t]) >= ts_vector[t] * requirement
            )
        end
    end
    return
end

function service_requirement_constraint!(
    optimization_container::OptimizationContainer,
    service::SR,
    ::ServiceModel{SR, T},
) where {SR <: PSY.StaticReserve, T <: AbstractReservesFormulation}
    initial_time = model_initial_time(optimization_container)
    @debug initial_time
    time_steps = model_time_steps(optimization_container)
    name = PSY.get_name(service)
    constraint =
        get_constraint(optimization_container, make_constraint_name(REQUIREMENT, SR))
    reserve_variable = get_variable(optimization_container, name, SR)
    use_slacks = get_services_slack_variables(optimization_container.settings)

    use_slacks && (slack_vars = reserve_slacks(optimization_container, name))

    requirement = PSY.get_requirement(service)
    for t in time_steps
        resource_expression = JuMP.GenericAffExpr{Float64, JuMP.VariableRef}()
        JuMP.add_to_expression!(resource_expression, sum(reserve_variable[:, t]))
        if use_slacks
            resource_expression += slack_vars[t]
        end
        constraint[name, t] = JuMP.@constraint(
            optimization_container.JuMPmodel,
            resource_expression >= requirement
        )
    end

    return
end

function cost_function!(
    optimization_container::OptimizationContainer,
    service::SR,
    ::ServiceModel{SR, T},
) where {SR <: PSY.Reserve, T <: AbstractReservesFormulation}
    reserve = get_variable(optimization_container, PSY.get_name(service), SR)
    for r in reserve
        JuMP.add_to_expression!(
            optimization_container.cost_function,
            r,
            DEFAULT_RESERVE_COST,
        )
    end
    return
end

function service_requirement_constraint!(
    optimization_container::OptimizationContainer,
    service::SR,
    ::ServiceModel{SR, StepwiseCostReserve},
) where {SR <: PSY.ReserveDemandCurve}
    initial_time = model_initial_time(optimization_container)
    @debug initial_time
    time_steps = model_time_steps(optimization_container)
    name = PSY.get_name(service)
    constraint =
        get_constraint(optimization_container, make_constraint_name(REQUIREMENT, SR))
    reserve_variable = get_variable(optimization_container, name, SR)
    requirement_variable = get_variable(optimization_container, SERVICE_REQUIREMENT, SR)

    for t in time_steps
        constraint[name, t] = JuMP.@constraint(
            optimization_container.JuMPmodel,
            sum(reserve_variable[:, t]) >= requirement_variable[name, t]
        )
    end

    return
end

_get_ramp_limits(::PSY.Component) = nothing
_get_ramp_limits(d::PSY.ThermalGen) = PSY.get_ramp_limits(d)
_get_ramp_limits(d::PSY.HydroGen) = PSY.get_ramp_limits(d)

function _get_data_for_ramp_limit(
    optimization_container::OptimizationContainer,
    service::SR,
    contributing_devices::U,
) where {
    SR <: PSY.Reserve,
    U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
} where {D <: PSY.Component}
    time_frame = PSY.get_time_frame(service)
    resolution = model_resolution(optimization_container)
    if resolution > Dates.Minute(1)
        minutes_per_period = Dates.value(Dates.Minute(resolution))
    else
        @warn("Not all formulations support under 1-minute resolutions. Exercise caution.")
        minutes_per_period = Dates.value(Dates.Second(resolution)) / 60
    end
    lenght_contributing_devices = length(contributing_devices)
    idx = 0
    data = Vector{ServiceRampConstraintInfo}(undef, lenght_contributing_devices)

    for d in contributing_devices
        name = PSY.get_name(d)
        non_binding_up = false
        non_binding_down = false
        ramp_limits = _get_ramp_limits(d)
        if !(ramp_limits === nothing)
            p_lims = PSY.get_active_power_limits(d)
            max_rate = abs(p_lims.min - p_lims.max) / time_frame
            if (ramp_limits.up >= max_rate) & (ramp_limits.down >= max_rate)
                @debug "Generator $(name) has a nonbinding ramp limits. Constraints Skipped"
                continue
            else
                idx += 1
            end
            ramp = (up = ramp_limits.up * time_frame, down = ramp_limits.down * time_frame)
            data[idx] = ServiceRampConstraintInfo(name, ramp)
        end
    end
    if idx < lenght_contributing_devices
        deleteat!(data, (idx + 1):lenght_contributing_devices)
    end
    return data
end

function ramp_constraints!(
    optimization_container::OptimizationContainer,
    service::SR,
    contributing_devices::U,
    ::ServiceModel{SR, T},
) where {
    SR <: PSY.Reserve{PSY.ReserveUp},
    T <: AbstractReservesFormulation,
    U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
} where {D <: PSY.Component}
    initial_time = model_initial_time(optimization_container)
    data = _get_data_for_ramp_limit(optimization_container, service, contributing_devices)
    service_name = PSY.get_name(service)
    if !isempty(data)
        service_upward_rateofchange!(
            optimization_container,
            data,
            make_constraint_name(RAMP, SR),
            make_variable_name(service_name, SR),
            service_name,
        )
    else
        @warn "Data doesn't contain contributing devices with ramp limits for service $service_name, consider adjusting your formulation"
    end
    return
end

function ramp_constraints!(
    optimization_container::OptimizationContainer,
    service::SR,
    contributing_devices::U,
    ::ServiceModel{SR, T},
) where {
    SR <: PSY.Reserve{PSY.ReserveDown},
    T <: AbstractReservesFormulation,
    U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
} where {D <: PSY.Component}
    initial_time = model_initial_time(optimization_container)
    data = _get_data_for_ramp_limit(optimization_container, service, contributing_devices)
    service_name = PSY.get_name(service)
    if !isempty(data)
        # TODO: `make_constraint_name` to build unique constraint name using  service_name
        service_downward_rateofchange!(
            optimization_container,
            data,
            make_constraint_name(RAMP, SR),
            make_variable_name(service_name, SR),
            service_name,
        )
    else
        @warn "Data doesn't contain contributing devices with ramp limits for service $service_name, consider adjusting your formulation"
    end
    return
end

function AddCostSpec(
    ::Type{T},
    ::Type{StepwiseCostReserve},
    optimization_container::OptimizationContainer,
) where {T <: PSY.Reserve}
    return AddCostSpec(;
        variable_type = ServiceRequirementVariable,
        component_type = T,
        has_status_variable = false,
        has_status_parameter = false,
        variable_cost = PSY.get_variable,
        start_up_cost = nothing,
        shut_down_cost = nothing,
        fixed_cost = nothing,
        sos_status = SOSStatusVariable.NO_VARIABLE,
    )
end

function add_to_cost!(
    optimization_container::OptimizationContainer,
    spec::AddCostSpec,
    service::SR,
    component_name::String,
) where {SR <: PSY.Reserve}
    time_steps = model_time_steps(optimization_container)
    use_forecast_data = model_uses_forecasts(optimization_container)
    if !use_forecast_data
        error("StepwiseCostReserve is only supported with forecast")
    end
    variable_cost_forecast =
        get_time_series(optimization_container, service, "variable_cost")
    variable_cost_forecast = map(PSY.VariableCost, variable_cost_forecast)
    for t in time_steps
        variable_cost!(
            optimization_container,
            spec,
            component_name,
            variable_cost_forecast[t],
            t,
        )
    end
    return
end

function cost_function!(
    optimization_container::OptimizationContainer,
    service::SR,
    model::ServiceModel{SR, StepwiseCostReserve},
) where {SR <: PSY.ReserveDemandCurve}
    spec = AddCostSpec(SR, model.formulation, optimization_container)
    @debug SR, spec
    add_to_cost!(optimization_container, spec, service, PSY.get_name(service))
    return
end

function modify_device_model!(
    devices_template::Dict{Symbol, DeviceModel},
    service_model::ServiceModel{<:PSY.Reserve, <:AbstractReservesFormulation},
    contributing_devices::Vector{<:PSY.Device},
)
    device_types = unique(typeof.(contributing_devices))
    for dt in device_types
        for (device_model_name, device_model) in devices_template
            # add message here when it exists
            get_component_type(device_model) != dt && continue
            service_model in device_model.services && continue
            push!(device_model.services, service_model)
        end
    end

    return
end

function include_service!(
    constraint_info::T,
    services,
    ::ServiceModel{SR, <:AbstractReservesFormulation},
) where {
    T <: Union{AbstractRangeConstraintInfo, AbstractRampConstraintInfo},
    SR <: PSY.Reserve{PSY.ReserveUp},
}
    for (ix, service) in enumerate(services)
        push!(
            constraint_info.additional_terms_ub,
            make_constraint_name(PSY.get_name(service), SR),
        )
    end
    return
end

function include_service!(
    constraint_info::T,
    services,
    ::ServiceModel{SR, <:AbstractReservesFormulation},
) where {
    T <: Union{AbstractRangeConstraintInfo, AbstractRampConstraintInfo},
    SR <: PSY.Reserve{PSY.ReserveDown},
}
    for (ix, service) in enumerate(services)
        push!(
            constraint_info.additional_terms_lb,
            make_constraint_name(PSY.get_name(service), SR),
        )
    end
    return
end

function include_service!(
    constraint_info::T,
    services,
    ::ServiceModel{SR, RampReserve},
) where {T <: AbstractRampConstraintInfo, SR <: PSY.Reserve{PSY.ReserveDown}}
    return
end

function include_service!(
    constraint_info::ReserveRangeConstraintInfo,
    services,
    ::ServiceModel{SR, RampReserve},
) where {SR <: PSY.Reserve{PSY.ReserveUp}}
    for (ix, service) in enumerate(services)
        # Should this be make_variable_name ?
        name = make_constraint_name(PSY.get_name(service), SR)
        push!(constraint_info.additional_terms_up, name)
        set_time_frame!(constraint_info, (name => PSY.get_time_frame(service)))
    end
    return
end

function include_service!(
    constraint_info::ReserveRangeConstraintInfo,
    services,
    ::ServiceModel{SR, RampReserve},
) where {SR <: PSY.Reserve{PSY.ReserveDown}}
    for (ix, service) in enumerate(services)
        # Should this be make_variable_name ?
        name = make_constraint_name(PSY.get_name(service), SR)
        push!(constraint_info.additional_terms_dn, name)
        set_time_frame!(constraint_info, (name => PSY.get_time_frame(service)))
    end
    return
end

function add_device_services!(
    constraint_info::T,
    device::D,
    model::DeviceModel,
) where {
    T <: Union{AbstractRangeConstraintInfo, AbstractRampConstraintInfo},
    D <: PSY.Device,
}
    for service_model in get_services(model)
        if PSY.has_service(device, service_model.component_type)
            services = (
                s for s in PSY.get_services(device) if isa(s, service_model.component_type)
            )
            @assert !isempty(services)
            include_service!(constraint_info, services, service_model)
        end
    end
    return
end

function add_device_services!(
    constraint_info::T,
    device::D,
    model::DeviceModel{D, BatteryAncillaryServices},
) where {
    T <: Union{AbstractRangeConstraintInfo, AbstractRampConstraintInfo},
    D <: PSY.Storage,
}
    return
end

function add_device_services!(
    constraint_info::ReserveRangeConstraintInfo,
    device::D,
    model::DeviceModel{D, BatteryAncillaryServices},
) where {D <: PSY.Storage}
    for service_model in get_services(model)
        if PSY.has_service(device, service_model.component_type)
            services = (
                s for s in PSY.get_services(device) if isa(s, service_model.component_type)
            )
            @assert !isempty(services)
            include_service!(constraint_info, services, service_model)
        end
    end
    return
end

function add_device_services!(
    constraint_data_in::AbstractRangeConstraintInfo,
    constraint_data_out::AbstractRangeConstraintInfo,
    device::D,
    model::DeviceModel{D, <:AbstractStorageFormulation},
) where {D <: PSY.Storage}
    for service_model in get_services(model)
        if PSY.has_service(device, service_model.component_type)
            services = (
                s for s in PSY.get_services(device) if isa(s, service_model.component_type)
            )
            @assert !isempty(services)
            if service_model.component_type <: PSY.Reserve{PSY.ReserveDown}
                for service in services
                    push!(
                        constraint_data_in.additional_terms_ub,
                        make_constraint_name(
                            PSY.get_name(service),
                            service_model.component_type,
                        ),
                    )
                end
            elseif service_model.component_type <: PSY.Reserve{PSY.ReserveUp}
                for service in services
                    push!(
                        constraint_data_out.additional_terms_ub,
                        make_constraint_name(
                            PSY.get_name(service),
                            service_model.component_type,
                        ),
                    )
                end
            end
        end
    end
    return
end

ext: jl


text:
using Test

using CSV
using DataFrames

using FlightMechanicsSimulator

@testset "aero.jl" begin

    α_test = LinRange(-10, 45, 20)
    β_test = LinRange(-30, 30, 20)

    de_test = LinRange(-25, 25, 20)
    da_test = LinRange(-21.5, 21.5, 20)
    dr_test = LinRange(-30, 30, 20)

    df = DataFrame(CSV.File("data/damp.csv"))
    for case in eachrow(df)
        rv1 = F16Stevens.damp(case.alpha)
        @test isapprox(rv1, Array(case[2:end]), atol = 1.0e-14)
    end

    df = DataFrame(CSV.File("data/cx.csv"))
    for case in eachrow(df)
        rv1 = F16Stevens.CX(case.alpha, case.de)
        @test isapprox(rv1, case.cx, atol = 1.0e-15)
    end

    df = DataFrame(CSV.File("data/cy.csv"))
    for case in eachrow(df)
        rv1 = F16Stevens.CY(case.beta, case.da, case.dr)
        @test isapprox(rv1, case.cy, atol = 1.0e-15)
    end

    df = DataFrame(CSV.File("data/cz.csv"))
    for case in eachrow(df)
        rv1 = F16Stevens.CZ(case.alpha, case.beta, case.de)
        @test isapprox(rv1, case.cz, atol = 1.0e-15)
    end

    df = DataFrame(CSV.File("data/cm.csv"))
    for case in eachrow(df)
        rv1 = F16Stevens.CM(case.alpha, case.de)
        @test isapprox(rv1, case.cm, atol = 1.0e-15)
    end

    df = DataFrame(CSV.File("data/aero_coeffs.csv"))
    for case in eachrow(df)
        rv1 = F16Stevens.CL(case.alpha, case.beta)
        @test isapprox(rv1, case.cl, atol = 1.0e-15)

        rv1 = F16Stevens.CN(case.alpha, case.beta)
        @test isapprox(rv1, case.cn, atol = 1.0e-15)

        rv1 = F16Stevens.DLDA(case.alpha, case.beta)
        @test isapprox(rv1, case.dlda, atol = 1.0e-15)

        rv1 = F16Stevens.DLDR(case.alpha, case.beta)
        @test isapprox(rv1, case.dldr, atol = 1.0e-15)

        rv1 = F16Stevens.DNDA(case.alpha, case.beta)
        @test isapprox(rv1, case.dnda, atol = 1.0e-15)

        rv1 = F16Stevens.DNDR(case.alpha, case.beta)
        @test isapprox(rv1, case.dndr, atol = 1.0e-15)
    end
end


text:
const p1 = Proj{1}()
const p2 = Proj{2}()
const p3 = Proj{3}()
const p4 = Proj{4}()

# add = PrimRec(p1, Comb(succ, p1))
const add = PrimRec(p1, succ(p1))

# mult = PrimRec(zro, Comb(add, p1, p3))
const mult = PrimRec(zro, add(p1, p3))

const pred = PrimRec(zro, p2)

# cosubrev = PrimRec(p1, Comb(pred, p1))
const cosubrev = PrimRec(p1, pred(p1))

# cosub = Comb(cosubrev, p2, p1)
const cosub = cosubrev(p2, p1)

# dff = Comb(add, cosub, cosubrev)
const dff = add(cosub, cosubrev)

# const1 = Comb(succ, zro)
const const1 = succ(zro)

# sgn = PrimRec(zro, Comb(const1, p1))
const sgn = PrimRec(zro, const1(p1))

# nsgn = PrimRec(const1, Comb(zro, p1))
const nsgn = PrimRec(const1, zro(p1))

# eq = Comb(nsgn, dff)
const eq = nsgn(dff)

# neq = Comb(sgn, dff)
const neq = sgn(dff)

# remainder = PrimRec(
#     zro,
#     Comb(
#         mult,
#         Comb(succ, p1),
#         Comb(neq, p1, Comb(pred, p3)),
#     ),
# )
const remainder = PrimRec(zro, mult(succ(p1), neq(p1, pred(p3))))

# ge = Comb(nsgn, cosubrev)
const ge = nsgn(cosubrev)

ifel = PrimRec(p2, p3)

# square = Comb(mult, p1, p1)
const square = mult(p1, p1)

# floor_sqrt = PrimRec(
#     zro,
#     Comb(ifel, Comb(eq, Comb(square, Comb(succ, p1)), Comb(succ, p2)), Comb(succ, p1), p1),
# )
const floor_sqrt = PrimRec(zro, ifel(eq(square(succ(p1)), succ(p2)), succ(p1), p1))

# traingle = PrimRec(zro, Comb(add, p1, p2))
const traingle = PrimRec(zro, add(p1, p2))

# pair = Comb(add, p1, Comb(traingle, Comb(add, p1, Comb(succ, p2))))
const pair = add(p1, traingle(add(p1, succ(p2))))

# fst = PrimRec(zro, Comb(mult, Comb(neq, Comb(pair, p1, zro), p2), Comb(succ, p1)))
const fst = PrimRec(zro, mult(neq(pair(p1, zro), p2), succ(p1)))

# scd = PrimRec(
#     zro,
#     Comb(
#         ifel,
#         Comb(eq, p2, Comb(pair, Comb(fst, p2), zro)),
#         Comb(succ, Comb(fst, p2)),
#         Comb(pred, p1),
#     ),
# )
const scd = PrimRec(zro, ifel(eq(p2, pair(fst(p2), zro)), succ(fst(p2)), pred(p1)))

# fib_pair = PrimRec(
#     Comb(pair, zro, const1),
#     Comb(pair, Comb(scd, p1), Comb(add, Comb(fst, p1), Comb(scd, p1))),
# )
const fib_pair = PrimRec(pair(zro, const1), pair(scd(p1), add(fst(p1), scd(p1))))

# fib = Comb(fst, fib_pair)
const fib = fst(fib_pair)


text:
function test_energy(mesh, sbp, eqn, opts)
  energy_final = calcNorm(eqn, eqn.q_vec)
  q_initial = zeros(eqn.q_vec)
  ICfunc = AdvectionEquationMod.ICDict[opts["IC_name"]]
  ICfunc(mesh, sbp, eqn, opts, q_initial)
  energy_initial = calcNorm(eqn, q_initial)

  @test  abs(energy_initial - energy_final)  < 1e-12
end

"""
  Test energy stability in 2 and 3 dimensions.
"""
function test_energy_serial()
  @testset "----- Testing Energy Stability -----" begin

    start_dir = pwd()
    cd(dirname(@__FILE__))

    fname = "input_vals_periodic.jl"
    cd("./2dp1")
    mesh, sbp, eqn, opts = solvePDE(fname)
    test_energy(mesh, sbp, eqn, opts)

    cd("../3dp1")
    mesh, sbp, eqn, opts = solvePDE(fname)
    test_energy(mesh, sbp, eqn, opts)

    cd(start_dir)

  end

  return nothing
end

#test_energy_serial()
add_func1!(AdvectionTests, test_energy_serial, [TAG_SHORTTEST])


text:
using BinaryBuilder

cuda_version = v"10.1.243"   # NOTE: could be less specific

sources_linux_x64 = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.0.3.33/10.1_20200825/cudnn-10.1-linux-x64-v8.0.3.33.tgz",
                  "4752ac6aea4e4d2226061610d6843da6338ef75a93518aa9ce50d0f58df5fb07")
]
sources_linux_ppc64le = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.0.3.33/10.1_20200825/cudnn-10.1-linux-ppc64le-v8.0.3.33.tgz",
                  "c546175f6ec86a11ee8fb9ab5526fa8d854322545769a87d35b1a505992f89c3")
]
sources_windows = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.0.3.33/10.1_20200825/cudnn-10.1-windows10-x64-v8.0.3.33.zip",
                  "9749e8a8f1dc0cfbb8a70d91825d0c5864426b44c3371ada1e077d8bd86b129f")
]

include("../common.jl")

max_stars_repo_name: mortenpi/Yggdrasil

max_issues_repo_name: mortenpi/Yggdrasil

max_forks_repo_name: mortenpi/Yggdrasil


text:
using GNSSTools
using PyPlot
pygui(true)

t_length = 1.  # second
f_s = 5e6  # Hz
v_0 = L1_freq  # Hz
h₋₂, h₋₁, h₀, h₁, h₂ = h_parms_tcxo[4]

noise_h₋₂ = generate_phase_noise(t_length, f_s, v_0, [h₋₂, 0,  0,  0,  0])
noise_h₋₁ = generate_phase_noise(t_length, f_s, v_0, [0, h₋₁,  0,  0,  0])
noise_h₀  = generate_phase_noise(t_length, f_s, v_0, [0,   0, h₀,  0,  0])
noise_h₁  = generate_phase_noise(t_length, f_s, v_0, [0,   0,  0, h₁,  0])
noise_h₂  = generate_phase_noise(t_length, f_s, v_0, [0,   0,  0,  0, h₂])

t = range(0, t_length, length=length(noise_h₋₂))
title_font_size = 10

fig = figure(figsize=(6.5,1.2))
ax1 = fig.add_subplot(1,5,1)
ax1.plot(t, real.(noise_h₋₂), "k")
title("h₋₂", fontsize=title_font_size)
tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
ax2 = fig.add_subplot(1,5,2)
ax2.plot(t, real.(noise_h₋₁), "k")
title("h₋₁", fontsize=title_font_size)
tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
ax3 = fig.add_subplot(1,5,3)
ax3.plot(t, real.(noise_h₀), "k")
title("h₀", fontsize=title_font_size)
tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
ax4 = fig.add_subplot(1,5,4)
ax4.plot(t, real.(noise_h₁), "k")
title("h₁", fontsize=title_font_size)
tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
ax5 = fig.add_subplot(1,5,5)
ax5.plot(t, real.(noise_h₂), "k")
title("h₂", fontsize=title_font_size)
tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
subplots_adjust(hspace=0.3, wspace=0.4, top=0.81, left=0.08, right=0.93, bottom=0.06)
savefig("figures/ch2_oscillator_noise.pdf", dpi=300)


text:
module Func

function clip_grads(grads,max_norm)
    total_norm = 0

    for i in 1:length(grads)
        total_norm += sum(grads[i].^2)
    end
    total_norm = sqrt(total_norm)

    rate = max_norm / (total_norm + 1e-6)
    if rate < 1
        for i in 1:length(grads)
            grads[i] .*= rate
        end
    end

end

end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# https://discourse.julialang.org/t/help-to-get-my-slow-julia-code-to-run-as-fast-as-rust-java-lisp/65741

# %%
# download dictionary
if !isfile("dictionary.txt")
    dictionary_url = "https://raw.githubusercontent.com/renatoathaydes/prechelt-phone-number-encoding/julia/dictionary.txt"
    Downloads.download(dictionary_url, "dictionary.txt")
end

# %%
# https://gist.github.com/jonathanBieler/de37a190d590297fc6b5d0ffee3c18dc

#=
# Port of Peter Norvig's Common Lisp program from http://norvig.com/java-lisp.html.
#
# - Julia version: 1.6.2
# - Author: Renato Athaydes
# - Date: 2021-07-24
=#
const emptyStrings = String[]

function printTranslations(num, digits, start=1, words=String[])
    if start > length(digits)
       return println(num, ": ", join(words, " "))
    end
    foundWord = false
    n = BigInt(1)
    for i in start:length(digits)
        n = n * 10 + nthDigit(digits, i)
        for word in get(dict, n, emptyStrings)
            foundWord = true
            printTranslations(num, digits, i + 1, [words; word])
        end
    end
    if !foundWord &&
        !(!isempty(words) && length(words[end]) == 1 && isdigit(words[end][begin]))
        printTranslations(num, digits, start + 1, [words; string(nthDigit(digits, start))])
    end
end

function loadDictionary(file)::Dict{BigInt, Vector{String}}
    local dict = Dict{BigInt, Vector{String}}()
    for word in eachline(file)
        push!(get!(dict, wordToNumber(word)) do; String[] end, word)
    end
    dict
end

function nthDigit(digits::String, i::Int64)::UInt
    UInt(digits[i]) - UInt('0')
end

function charToDigit(ch::Char)::UInt
    ch = lowercase(ch)
    ch == 'e' && return 0
    ch in ['j', 'n', 'q'] && return 1
    ch in ['r', 'w', 'x'] && return 2
    ch in ['d', 's', 'y'] && return 3
    ch in ['f', 't'] && return 4
    ch in ['a', 'm'] && return 5
    ch in ['c', 'i', 'v'] && return 6
    ch in ['b', 'k', 'u'] && return 7
    ch in ['l', 'o', 'p'] && return 8
    ch in ['g', 'h', 'z'] && return 9
    throw(DomainError(ch, "Not a letter"))
end

function wordToNumber(word::String)::BigInt
    n = BigInt(1)
    for ch in word
        if isletter(ch) && isascii(ch)
            n = n * 10 + charToDigit(ch)
        end
    end
    n
end

# dict = open(isempty(ARGS) ? "tests/words.txt" : ARGS[begin]) do file
#     loadDictionary(file)
# end

# open(length(ARGS) < 2 ? "tests/numbers.txt" : ARGS[begin+1]) do file
#     for num in eachline(file)
#         printTranslations(num, filter(isdigit, num))
#     end
# end

# %%
# original printTranslations added the arguments io and dict
function printTranslations_original(io, dict, num, digits, start=1, words=String[])
    if start > length(digits)
       return println(io, num, ": ", join(words, " "))
    end
    foundWord = false
    n = BigInt(1)
    for i in start:length(digits)
        n = n * 10 + nthDigit(digits, i)
        for word in get(dict, n, emptyStrings)
            foundWord = true
            printTranslations_original(io, dict, num, digits, i + 1, [words; word])
        end
    end
    if !foundWord &&
        !(!isempty(words) && length(words[end]) == 1 && isdigit(words[end][begin]))
        printTranslations_original(io, dict, num, digits, start + 1, [words; string(nthDigit(digits, start))])
    end
end

# translate numbers by printTranslations with dict and print the result to io
function translate(io::IO, printTranslations, dict, numbers)
    (num -> printTranslations(io, dict, num, filter(isdigit, num))).(numbers)
end

# %%
# generate the test data `numbers` of length n = 10^6

using Random

const _chars = ('0', '1', '2', '3', '4', '5', '6', '7', '8', '9', '/', '-')
const _maxlen = 50
const _allowempty = false

function randphone(; chars=_chars, maxlen=_maxlen, allowempty=_allowempty)
    while true
        phone = randstring(chars, rand(1:maxlen))
        (allowempty || occursin(r"[0-9]", phone)) && return phone
    end
end

function randphone(n; chars=_chars, maxlen=_maxlen, allowempty=_allowempty)
    [randphone(; chars, maxlen, allowempty) for _ in 1:n]
end

n = 10^6
println("numbers = randphone($n):")
@time numbers = randphone(n)
println()

# %%
# load dictionary
println("open(loadDictionary, \"dictionary.txt\"):")
@time dict = open(loadDictionary, "dictionary.txt")
println()

# %%
# translate numbers with dict and print the result to io
io = IOBuffer()
println("translate by printTranslations_original:")
@time translate(io, printTranslations_original, dict, numbers)
println()
result_original = String(take!(io))
close(io)

# %%
# https://gist.github.com/jonathanBieler/de37a190d590297fc6b5d0ffee3c18dc

#=
# Port of Peter Norvig's Common Lisp program from http://norvig.com/java-lisp.html.
#
# - Julia version: 1.6.2
# - Author: Renato Athaydes
# - Date: 2021-07-24
=#
# const emptyStrings = String[]

function printTranslations(io, dict, num, digits, start=1, words=String[])
    if start > length(digits)
       return println(io, num, ": ", join(words, " "))
    end
    foundWord = false
    n = BigInt(1)
    for i in start:length(digits)
        Base.GMP.MPZ.mul_si!(n, 10)
        Base.GMP.MPZ.add_ui!(n, nthDigit(digits, i))
        for word in get(dict, n, emptyStrings)
            foundWord = true
            printTranslations(io, dict, num, digits, i + 1, [words; word])
        end
    end
    if !foundWord &&
        !(!isempty(words) && length(words[end]) == 1 && isdigit(words[end][begin]))
        printTranslations(io, dict, num, digits, start + 1, [words; string(nthDigit(digits, start))])
    end
end

function loadDictionary(file)::Dict{BigInt, Vector{String}}
    local dict = Dict{BigInt, Vector{String}}()
    for word in eachline(file)
        push!(get!(dict, wordToNumber(word)) do; String[] end, word)
    end
    dict
end

function nthDigit(digits::String, i::Int64)
    UInt(digits[i]) - UInt('0')
end

function charToDigit(ch::Char)
    ch = lowercase(ch)
    ch == 'e' && return 0
    ch in ('j', 'n', 'q') && return 1
    ch in ('r', 'w', 'x') && return 2
    ch in ('d', 's', 'y') && return 3
    ch in ('f', 't') && return 4
    ch in ('a', 'm') && return 5
    ch in ('c', 'i', 'v') && return 6
    ch in ('b', 'k', 'u') && return 7
    ch in ('l', 'o', 'p') && return 8
    ch in ('g', 'h', 'z') && return 9
    throw(DomainError(ch, "Not a letter"))
end

function wordToNumber(word::String)
    n = BigInt(1)
    for ch in word
        if isletter(ch) && isascii(ch)
            Base.GMP.MPZ.mul_si!(n, 10)
            Base.GMP.MPZ.add_ui!(n, charToDigit(ch))
        end
    end
    n
end

# patch in method to add integer to BigInt in-place
@eval Base.GMP.MPZ begin
    add_ui!(x::BigInt, a::BigInt, b) = (ccall((:__gmpz_add_ui, :libgmp), Cvoid, (mpz_t, mpz_t, Clong), x, a, b); x)
    add_ui!(x::BigInt, b) = add_ui!(x, x, b)
end

function main()

    path = "/Users/jbieler/Downloads/tmp/prechelt-phone-number-encoding/"

    dict = open("$(path)dictionary.txt") do file
        loadDictionary(file)
    end
    io = IOBuffer()
    open("$(path)input.txt") do file
        for num in eachline(file)
            printTranslations(io, dict, num, filter(isdigit, num))
        end
    end

end

# main() #make sure it's compiled before timing it
# @time main()

# %%
io = IOBuffer()
println("translate by printTranslations of jonathanBieler:")
@time translate(io, printTranslations, dict, numbers)
println()
result_jonathanBieler = String(take!(io))
close(io)

# %%
@show result_original == result_jonathanBieler;

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
module BFParser

import REPL
using REPL.TerminalMenus

struct ParseError <: Exception
    msg::AbstractString
    stacktrace::Base.StackTraces.StackTrace
end
struct ErrorWrapper <: Exception
    err::Exception
    backtraces::Vector
end

abstract type AbstractBF end
tosource(t::AbstractBF) = sprint(tosource, t)

tosource(terms::Vector{AbstractBF}) = sprint(tosource, terms)
function tosource(io::IO, terms::Vector{AbstractBF})
    for line in terms
        tosource(io, line)
    end
end

# BF Script
struct BF <: AbstractBF
    lines::Vector{AbstractBF}
end
tosource(io::IO, t::BF) = tosource(io, t.lines)

# BF Loop
struct BFLE <: AbstractBF
    lines::Vector{AbstractBF}
    lineno::Int
end
_lineno(t::BFLE) = t.lineno
function tosource(io::IO, t::BFLE)
    print(io, '[')
    tosource(io, t.lines)
    print(io, ']')
end

# BF Increment Pointer
struct BFI{N} <: AbstractBF
    lineno::Int
end
_lineno(t::BFI) = t.lineno
function tosource(io::IO, ::BFI{N}) where N
    for _=1:N
        print(io, '>')
    end
end

# BF Decrement Pointer
struct BFD{N} <: AbstractBF
    lineno::Int
end
_lineno(t::BFD) = t.lineno
function tosource(io::IO, ::BFD{N}) where N
    for _=1:N
        print(io, '<')
    end
end

# BF Increment Value
struct BFP{N} <: AbstractBF
    lineno::Int
end
_lineno(t::BFP) = t.lineno
function tosource(io::IO, ::BFP{N}) where N
    for _=1:N
        print(io, '+')
    end
end

# BF Decrement Value
struct BFM{N} <: AbstractBF
    lineno::Int
end
_lineno(t::BFM) = t.lineno
function tosource(io::IO, ::BFM{N}) where N
    for _=1:N
        print(io, '-')
    end
end

# BF Input
struct BFR <: AbstractBF
    lineno::Int
end
_lineno(t::BFR) = t.lineno
tosource(io::IO, ::BFR) = print(io, ',')

# BF Output
struct BFW <: AbstractBF
    lineno::Int
end
_lineno(t::BFW) = t.lineno
tosource(io::IO, ::BFW) = print(io, '.')

function parseBF(::Type{BF}, src::String)
    _lines = AbstractBF[]
    _stack = Tuple{Int, Vector{AbstractBF}}[]
    re = r"\[|\]|<+|>+|\++|-+|,|\."
    for m in eachmatch(re, src)
        if m.match == "["
            push!(_stack, (m.offset, _lines))
            _lines = AbstractBF[]
        elseif m.match == "]"
            isempty(_stack) && throw(ParseError("syntax: ']' mismatch", [Base.StackTraces.StackFrame("]", src, m.offset)]))
            _lnum, __lines = pop!(_stack)
            bf = BFLE(_lines, _lnum)
            _lines = __lines
            push!(_lines, bf)
        elseif startswith(m.match, '>')
            push!(_lines, BFI{length(m.match)}(m.offset))
        elseif startswith(m.match, '<')
            push!(_lines, BFD{length(m.match)}(m.offset))
        elseif startswith(m.match, '+')
            push!(_lines, BFP{length(m.match)}(m.offset))
        elseif startswith(m.match, '-')
            push!(_lines, BFM{length(m.match)}(m.offset))
        elseif m.match == ","
            push!(_lines, BFR(m.offset))
        elseif m.match == "."
            push!(_lines, BFW(m.offset))
        end
    end
    if !isempty(_stack)
        _src = "[" * tosource(_lines)
        stacktrace = [Base.StackTraces.StackFrame("[", _src, 1)]
        _fn = _src
        _lnum, _lines = pop!(_stack)
        while !isempty(_stack)
            _src = "[" * tosource(_lines) * _fn
            __lnum, __lines = pop!(_stack)
            push!(stacktrace, Base.StackTraces.StackFrame(_fn, _src, _lnum - __lnum + 1))
            _fn = _src
            _lnum = __lnum
            _lines = __lines
        end
        _src = tosource(_lines) * _fn
        push!(stacktrace, Base.StackTraces.StackFrame(_fn, _src, _lnum))
        throw(ParseError("syntax: '[' mismatch", stacktrace))
    end
    BF(_lines)
end

# BF terms which has length
const BFN = Union{BFI, BFD, BFP, BFM}

# BF terms which has no length
const BFS = Union{BFLE, BFR, BFW}

# BF terms
const BFT = Union{BFN, BFS}

offset(t1::BFT, t2::BFT) = _lineno(t1) - _lineno(t2)

# BF Interpreter
mutable struct BFInterpreter{I <: IO, O <: IO}
    tape::Vector{UInt8}
    ptr::Int
    inio::I
    outio::O
    (::Type{BFInterpreter})(inio::I, outio::O, n::Int=30000) where {I <: IO, O <: IO} =
        new{I, O}(zeros(UInt8, n), 1, inio, outio)
end
(::Type{BFInterpreter})(input::AbstractString, n::Int=30000) =
    BFInterpreter(IOBuffer(input), IOBuffer(), n)
(::Type{BFInterpreter})(input::AbstractVector{UInt8}, n::Int=30000) =
    BFInterpreter(IOBuffer(input), IOBuffer(), n)
(::Type{BFInterpreter})(io::IO, n::Int=30000) =
    BFInterpreter(stdin, io, n)
(::Type{BFInterpreter})(n::Int=30000) =
    BFInterpreter(stdin, stdout, n)

# BF Execute Increment Pointer
function Base.run(bfi::BFInterpreter, ::BFI{N}) where {N}
    bfi.ptr += N
end

# BF Execute Decrement Pointer
function Base.run(bfi::BFInterpreter, ::BFD{N}) where {N}
    bfi.ptr -= N
end

# BF Execute Increment Value
function Base.run(bfi::BFInterpreter, ::BFP{N}) where {N}
    try
        bfi.tape[bfi.ptr] += N % UInt8
    catch ex
        isa(ex, BoundsError) && throw(ErrorWrapper(ex, Base.StackTraces.backtrace()))
        rethrow()
    end
end

# BF Execute Decrement Value
function Base.run(bfi::BFInterpreter, ::BFM{N}) where {N}
    try
        bfi.tape[bfi.ptr] -= N % UInt8
    catch ex
        isa(ex, BoundsError) && throw(ErrorWrapper(ex, Base.StackTraces.backtrace()))
        rethrow()
    end
end

# BF Execute Input
function Base.run(bfi::BFInterpreter, ::BFR)
    try
        bfi.tape[bfi.ptr] = eof(bfi.inio) ? 0xff : read(bfi.inio, UInt8)
    catch ex
        isa(ex, EOFError) && return 0xff
        isa(ex, BoundsError) && throw(ErrorWrapper(ex, Base.StackTraces.backtrace()))
        rethrow()
    end
end

# BF Execute Output
function Base.run(bfi::BFInterpreter, ::BFW)
    try
        write(bfi.outio, bfi.tape[bfi.ptr])
    catch ex
        isa(ex, BoundsError) && throw(ErrorWrapper(ex, Base.StackTraces.backtrace()))
        rethrow()
    end
end

# BF Execute Loop
function Base.run(bfi::BFInterpreter, bfle::BFLE)
    try
        while bfi.tape[bfi.ptr] != 0
            for bf in bfle.lines
                run(bfi, bf)
            end
        end
    catch ex
        isa(ex, BoundsError) && throw(ErrorWrapper(ex, Base.StackTraces.backtrace()))
        rethrow()
    end
end

# BF Execute Script
function Base.run(bfi::BFInterpreter, bfs::BF)
    for bf in bfs.lines
        run(bfi, bf)
    end
end

end # module


text:
# The code to generate Fig 5
# The relevant statistics reported in the paper are calculated at the end of the
# file
using PyPlot
using SurNoR_2020
using Statistics
using MAT
using HypothesisTests

PyPlot.svg(true)
rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["svg.fonttype"] = "none"
rcParams["pdf.fonttype"] = 42

# ------------------------------------------------------------------------------
# Training Algorithms
# ------------------------------------------------------------------------------
for Part_ind = 1:4  # For 4 different repetitions of the optimization
    # Main alternative algorithms for Fig 5
    Agents_MBQ0S    = SurNoR_2020.Func_MBQ0S_CV_train(Part_ind=Part_ind)        # MB + S + OI
    Agents_MBSExp   = SurNoR_2020.Func_MBSExp_CV_train(Part_ind=Part_ind)       # MB + S + U
    Agents_MB_Leaky = SurNoR_2020.Func_MB_Leaky_CV_train(Part_ind=Part_ind)     # MB + N
    Agents_MBS      = SurNoR_2020.Func_MBS_CV_train(Part_ind=Part_ind)          # MB + S + N
    Agents_MFQ0     = SurNoR_2020.Func_MFQ0_CV_train(Part_ind=Part_ind)         # MF + OI
    Agents_MFSExp   = SurNoR_2020.Func_MFSExp_CV_train(Part_ind=Part_ind)       # MF + S + U
    Agents_MFN      = SurNoR_2020.Func_MFN_CV_train(Part_ind=Part_ind)          # MF + N
    Agents_MFNS     = SurNoR_2020.Func_MFNS_CV_train(Part_ind=Part_ind)         # MF + S + N
    Agents_HybQ0S   = SurNoR_2020.Func_HybQ0S_CV_train(Part_ind=Part_ind)       # Hyb + S + OI
    Agents_HybSExp  = SurNoR_2020.Func_HybSExp_CV_train(Part_ind=Part_ind)      # Hyb + S + U
    Agents_HybN     = SurNoR_2020.Func_HybN_CV_train(Part_ind=Part_ind)         # Hyb + N
    Agents_SurNoR   = SurNoR_2020.Func_SurNoR_CV_train(Part_ind=Part_ind)       # SurNoR (Hyb + S + N)

    # Control algorithms
    Agents_HybSTrapDet           = SurNoR_2020.Func_HybSTrapDet_CV_train(Part_ind=Part_ind)         # Binary Novelty 1
    Agents_HybSTrapDet2          = SurNoR_2020.Func_HybSTrapDet2_CV_train(Part_ind=Part_ind)        # Binary Novelty 2
    Agents_SurNoR_NovInExploit   = SurNoR_2020.Func_SurNoR_NovInExploit_CV_train(Part_ind=Part_ind) # SurNoR with β_N ≧ 0 also in Epi > 1

    if Part_ind == 1
        global Agents_Tot_Alg = cat(Agents_MBQ0S, Agents_MBSExp, Agents_MB_Leaky, Agents_MBS,
                                    Agents_MFQ0, Agents_MFSExp, Agents_MFN, Agents_MFNS,
                                    Agents_HybQ0S, Agents_HybSExp, Agents_HybN, Agents_SurNoR,
                                    Agents_HybSTrapDet, Agents_HybSTrapDet2, Agents_SurNoR_NovInExploit,  dims=4)

    else
        Agents_Tot_Alg_temp = cat(Agents_MBQ0S, Agents_MBSExp, Agents_MB_Leaky, Agents_MBS,
                                  Agents_MFQ0, Agents_MFSExp, Agents_MFN, Agents_MFNS,
                                  Agents_HybQ0S, Agents_HybSExp, Agents_HybN, Agents_SurNoR,
                                  Agents_HybSTrapDet, Agents_HybSTrapDet2, Agents_SurNoR_NovInExploit, dims=4)

        Agents_Tot_Alg = cat(Agents_Tot_Alg, Agents_Tot_Alg_temp, dims=5)
    end
end
SurNoR_ind = 12   # index corresponding to SurNoR

# Extracting Log-likelihoods
Sub_Num = 12
Alg_Num = size(Agents_Tot_Alg)[4]
Part_Num = size(Agents_Tot_Alg)[5]
LogL_Matrix = zeros((Sub_Num,5,2,Alg_Num,Part_Num))
LogL_Unif_Matrix = zeros((Sub_Num,5,2,Alg_Num,Part_Num))

for Sub = 1:Sub_Num
    for Epi = 1:5
        for Block = 1:2
            for Alg = 1:Alg_Num
                for Part_ind = 1:Part_Num
                    LogL_Matrix[Sub,Epi,Block,Alg,Part_ind] =
                            Agents_Tot_Alg[Sub,Block,Epi,Alg,Part_ind].LogL
                    LogL_Unif_Matrix[Sub,Epi,Block,Alg,Part_ind] =
                            Agents_Tot_Alg[Sub,Block,Epi,Alg,Part_ind].LogL_Unif
                end
            end
        end
    end
end

LogL_Matrix_Epi1 = LogL_Matrix[:,1,1,:,:]
LogL_Unif_Matrix_Epi1 = LogL_Unif_Matrix[:,1,1,:,:]

LogL_Matrix_Epi1B2 = LogL_Matrix[:,1,2,:,:]
LogL_Unif_Matrix_Epi1B2 = LogL_Unif_Matrix[:,1,2,:,:]

LogL_Sum = sum(LogL_Matrix,dims=1)
LogL_Unif_Sum = sum(LogL_Unif_Matrix,dims=1)

LogL_Sum_Epi1 = sum(LogL_Matrix_Epi1,dims=1)
LogL_Unif_Sum_Epi1 = sum(LogL_Unif_Matrix_Epi1,dims=1)

LogL_Sum_Epi1B2 = sum(LogL_Matrix_Epi1B2,dims=1)
LogL_Unif_Sum_Epi1B2 = sum(LogL_Unif_Matrix_Epi1B2,dims=1)

dLogL_Sum_Block = sum(LogL_Sum - LogL_Unif_Sum,dims=2)
dLogL_Sum_Block_Epi1 = LogL_Sum_Epi1 - LogL_Unif_Sum_Epi1
dLogL_Sum_Block_Epi1B2 = LogL_Sum_Epi1B2 - LogL_Unif_Sum_Epi1B2

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Fig 5A: Log Evidences
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
Alg_Num_plot = 12
x = 1:Alg_Num_plot
Colors = ["#69C6E5","#45A2C1","#167CB7","#0F5680",
          "#E56969","#C14545","#B71616","#6E0D0D",
          "#D7A8F0","#C785EA","#B159E0","#8723BE"]

Legend = ["MB+S+OI", "MB+S+U", "Leak+N", "MB+S+N",
          "MF+OI" , "MF+S+U", "MF+N", "MF+S+N",
          "Hyb+S+OI", "Hyb+S+U", "Hyb+N", "SurNoR"]

# All data
y = dLogL_Sum_Block[1,1,1,:,:] .+ dLogL_Sum_Block[1,1,2,:,:]
dy = std(y,dims=2)/sqrt(Part_Num)
y = mean(y,dims=2)
fig = figure(figsize=(6,6.5)); ax = gca()
for i = 1:Alg_Num_plot
   ax.bar(x[i],y[i],color = Colors[i])
   ax.errorbar(x[i],y[i],yerr=dy[i],color="k",linewidth=1,drawstyle="steps",capsize=3)
end
title("Delta Log Evidence - w.r.t. RC - all epi. for both blocks")
ax.set_xticks(1:Alg_Num_plot)
ax.set_ylim([450,1300])
ax.set_xticklabels(Legend,rotation=90)
ax.set_ylabel("delta Log Evidence");
fig.subplots_adjust(bottom = 0.2, left = 0.15)

# Episode 1 of block 1
y = dLogL_Sum_Block_Epi1
dy = std(y,dims=3)/sqrt(Part_Num)
y = mean(y,dims=3)
fig = figure(figsize=(6,6.5)); ax = gca()
for i = 1:Alg_Num_plot
   ax.bar(x[i],y[i],color = Colors[i])
   ax.errorbar(x[i],y[i],yerr=dy[i],color="k",linewidth=1,drawstyle="steps",capsize=3)
end
title("Delta Log Evidence - w.r.t. RC - B1E1")
ax.set_xticks(1:Alg_Num_plot)
ax.set_ylim([-20,400])
ax.set_xticklabels(Legend,rotation=90)
title("Delta Log Evidence - w.r.t. RC - B1E1")
fig.subplots_adjust(bottom = 0.2, left = 0.15)

# Episode 1 of block 2
y = dLogL_Sum_Block_Epi1B2
dy = std(y,dims=3)/sqrt(Part_Num)
y = mean(y,dims=3)
fig = figure(figsize=(6,6.5)); ax = gca()
for i = 1:Alg_Num_plot
   ax.bar(x[i],y[i],color = Colors[i])
   ax.errorbar(x[i],y[i],yerr=dy[i],color="k",linewidth=1,drawstyle="steps",capsize=5)
end
title("Delta Log Evidence - w.r.t. RC - B2E1")
ax.set_xticks(1:Alg_Num_plot)
ax.set_ylim([-10,200])
ax.set_xticklabels(Legend,rotation=90)
ax.set_ylabel("delta Log Evidence");
fig.subplots_adjust(bottom = 0.2, left = 0.15)

# ------------------------------------------------------------------------------
# Save LogL as .MAT for Bayesian Model Selection (KE Stephan, et al. 2009)
# ------------------------------------------------------------------------------
LogL_MAT = mean(LogL_Matrix,dims=5)
LogL_MAT = cat(LogL_MAT[:,:,:,1:Alg_Num_plot,1],
               LogL_Unif_Matrix[:,:,:,1,1], dims=4)

file = matopen("./data/Log_evidences/LogL_Matrix.mat", "w")
write(file, "LogL_Matrix", LogL_MAT)
close(file)

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Fig 5B: Expected posterior
# WARNING: Before running this part, in principle, one needs to run the MATLAB
#          code "./data/Log_evidences/LogL_BMS.m" to use package SPM12 and
#          produce the results and save the file "./data/Log_evidences/BMS_Data.mat"
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
BMS = matread("./data/Log_evidences/BMS_Data.mat")

Colors = ["#69C6E5","#45A2C1","#167CB7","#0F5680",
          "#E56969","#C14545","#B71616","#6E0D0D",
          "#D7A8F0","#C785EA","#B159E0","#8723BE", "#474848"]

Legend = ["MB+S+OI", "MB+S+U", "Leak+N", "MB+S+N",
          "MF+OI" , "MF+S+U", "MF+N", "MF+S+N",
          "Hyb+S+OI", "Hyb+S+U", "Hyb+N", "SurNoR", "RC"]

fig = figure(figsize=(7,6)); ax = gca()
for i = 1:13
   ax.bar(i,BMS["BMS"]["exp_r"][i], color=Colors[i])
   if (i==SurNoR_ind)
        ax.text(i,BMS["BMS"]["exp_r"][i]+0.01,
                string("ϕ = ",string(round(BMS["BMS"]["pxp"][i],digits=2))),
                fontsize=12, horizontalalignment="center", rotation=0)
   end
end
ax.plot([0,14],[1,1] ./ 13, linestyle="dashed",linewidth=1,color="k")
title("Posterior Probabilities for Different Models")
ax.set_xticks(1:13)
ax.set_xticklabels(Legend,rotation=90)
ax.set_yticks(0:0.1:1)
ax.set_yticklabels(0:0.1:1)
ax.set_ylabel("Model Posterior Probability");
ax.set_ylim([0,0.8])
ax.set_xlim([0.25,13.75])
fig.subplots_adjust(bottom = 0.2, left = 0.15)

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Fig 5C: Accuracy rates
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Reading the accuracy rates
Sub_Num = 12
Acc_Matrix = zeros((Sub_Num,5,2,Alg_Num,Part_Num))      # To save accuracy rates
NCor_Matrix = zeros((Sub_Num,5,2,Alg_Num,Part_Num))     # To save number of correctly predicted actions
NTot_Matrix = zeros((Sub_Num,5,2,Alg_Num,Part_Num))     # To save total number of actions

for Sub = 1:Sub_Num
    for Epi = 1:5
        for Block = 1:2
            for Alg = 1:Alg_Num
                for Part_ind = 1:Part_Num
                    Acc_Matrix[Sub,Epi,Block,Alg,Part_ind] =
                            Agents_Tot_Alg[Sub,Block,Epi,Alg,Part_ind].Acc_Rate
                    NCor_Matrix[Sub,Epi,Block,Alg,Part_ind] =
                            sum(Agents_Tot_Alg[Sub,Block,Epi,Alg,Part_ind].Sort_Dist_Data[4,:])
                    NTot_Matrix[Sub,Epi,Block,Alg,Part_ind] =
                            length(Agents_Tot_Alg[Sub,Block,Epi,Alg,Part_ind].Est_Act)
                end
            end
        end
    end
end

Acc_Mean_4Parts = mean(Acc_Matrix,dims=1)
Acc_Std_4Parts = std(Acc_Matrix,dims=1)/sqrt(Sub_Num)

Acc_Mean = mean(Acc_Mean_4Parts,dims=5)
Acc_Std = sqrt.(mean(Acc_Std_4Parts,dims=5).^2 .+ var(Acc_Mean_4Parts,dims=5)./4)

# Computing entropy
Average_Entropy_Matrix = zeros((Sub_Num,5,2,Alg_Num,Part_Num))
for Sub = 1:Sub_Num
    for Epi = 1:5
        for Block = 1:2
            for Alg = 1:Alg_Num
                for Part_ind = 1:Part_Num
                    temp = Agents_Tot_Alg[Sub,Block,Epi,Alg,Part_ind].P_Action
                    temp = - temp .* log.(temp)
                    Average_Entropy_Matrix[Sub,Epi,Block,Alg,Part_ind] = mean(sum(temp,dims=1))
                end
            end
        end
    end
end

Average_Entropy_Mean_4Parts = mean(Average_Entropy_Matrix,dims=1)
Average_Entropy_Std_4Parts = std(Average_Entropy_Matrix,dims=1)/sqrt(Sub_Num)

Average_Entropy_Mean = mean(Average_Entropy_Mean_4Parts,dims=5)
Average_Entropy_Std = sqrt.(mean(Average_Entropy_Std_4Parts,dims=5).^2 .+ var(Average_Entropy_Mean_4Parts,dims=5)./4)

Alg_plot = SurNoR_ind  # Choosing SurNoR

fig, (ax1, ax2) = subplots(1, 2, figsize=(10,6.5))

x = 1:0.5:3

y = 100 * Acc_Mean[1,:,1,Alg_plot,1]
dy = 100 * Acc_Std[1,:,1,Alg_plot,1]
ax1.bar(x, y, width=0.3, align="center", color = "#8723BE")
ax1.errorbar(x,y[:],yerr=dy[:],color="k",linewidth=1,capsize=5,linestyle ="")

y = 100 * Acc_Mean[1,:,2,Alg_plot,1]
dy = 100 * Acc_Std[1,:,2,Alg_plot,1]
ax2.bar(x, y, width=0.3, align="center", color = "#8723BE")
ax2.errorbar(x,y[:],yerr=dy[:],color="k",linewidth=1,capsize=5,linestyle ="")

ax1.set_xticks(x)
ax1.set_xticklabels(["1", "2", "3", "4", "5"])
ax2.set_xticks(x)
ax2.set_xticklabels(["1", "2", "3", "4", "5"])
ax1.set_yticks([0,25,50,75,100])
ax1.set_yticklabels(["0%","25%", "50%", "75%", "100%"])
ax2.set_yticks([])
ax2.set_yticklabels([])

ax1.set_ylim(0,100); ax2.set_ylim(0,100);
ax1.set_xlim(0.5,3.5); ax2.set_xlim(0.5,3.5);
ax1.set_xlabel("Episode"); ax2.set_xlabel("Episode")
ax1.set_ylabel("Accuracy Rate of SurNoR in Predicting Actions");
ax1.set_title("1st Block"); ax2.set_title("2nd Block")

ax12 = ax1.twinx()
y = Average_Entropy_Mean[1,:,1,Alg_plot,1]
dy = Average_Entropy_Std[1,:,1,Alg_plot,1]
ax12.errorbar(x, y[:], yerr=dy[:],color="#2D2D2D", linestyle="dashed",capsize=5)

ax22 = ax2.twinx()
y = Average_Entropy_Mean[1,:,2,Alg_plot,1]
dy = Average_Entropy_Std[1,:,2,Alg_plot,1]
ax22.errorbar(x, y[:], yerr=dy[:],color="#2D2D2D", linestyle="dashed",capsize=5)

ax12.set_ylim([0,1.3])
ax22.set_ylim([0,1.3])
ax12.set_yticks([])
ax12.set_yticklabels([])

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Reporting the accuracy rates of different phases of SurNoR
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Reporting the accuracy after the first time finding the goal
NTot_AG = sum(sum(NTot_Matrix[:,2:end,:,SurNoR_ind:SurNoR_ind,:],dims=2),dims=3) .+
          NTot_Matrix[:,1:1,2:2,SurNoR_ind:SurNoR_ind,:]
NCor_AG = sum(sum(NCor_Matrix[:,2:end,:,SurNoR_ind:SurNoR_ind,:],dims=2),dims=3) .+
          NCor_Matrix[:,1:1,2:2,SurNoR_ind:SurNoR_ind,:]
Acc_AG = NCor_AG ./ NTot_AG

Acc_AG_Mean_4Parts = mean(Acc_AG,dims=1)
Acc_AG_Std_4Parts = std(Acc_AG,dims=1)/sqrt(Sub_Num)

Acc_AG_Mean = mean(Acc_AG_Mean_4Parts,dims=5)
Acc_AG_Std = sqrt.(mean(Acc_AG_Std_4Parts,dims=5).^2 .+ var(Acc_AG_Mean_4Parts,dims=5)./4)

println("Accuracy rate after the 1st time finiding the goal:")
@show Acc_AG_Mean[1];
println("SEM of the accuracy rate after the 1st time finiding the goal:")
@show Acc_AG_Std[1];
println("Total number of actions after the 1st time finiding the goal:")
@show sum(NTot_AG)/Part_Num

# Reporting the accuracy before the first time finding the goal
NTot_BG = sum(sum(NTot_Matrix[:,1:1,1:1,SurNoR_ind:SurNoR_ind,:],dims=2),dims=3)
NCor_BG = sum(sum(NCor_Matrix[:,1:1,1:1,SurNoR_ind:SurNoR_ind,:],dims=2),dims=3)
Acc_BG = NCor_BG ./ NTot_BG

Acc_BG_Mean_4Parts = mean(Acc_BG,dims=1)
Acc_BG_Std_4Parts = std(Acc_BG,dims=1)/sqrt(Sub_Num)

Acc_BG_Mean = mean(Acc_BG_Mean_4Parts,dims=5)
Acc_BG_Std = sqrt.(mean(Acc_BG_Std_4Parts,dims=5).^2 .+ var(Acc_BG_Mean_4Parts,dims=5)./4)

println("Accuracy rate before the 1st time finiding the goal (B1E1):")
@show Acc_BG_Mean[1];
println("SEM of the accuracy rate before the 1st time finiding the goal (B1E1):")
@show Acc_BG_Std[1];
println("Total number of actions before the 1st time finiding the goal (B1E1):")
@show sum(NTot_BG)/Part_Num

# Reporting the accuracy of all actions
NTot_All = sum(sum(NTot_Matrix[:,:,:,SurNoR_ind:SurNoR_ind,:],dims=2),dims=3)
NCor_All = sum(sum(NCor_Matrix[:,:,:,SurNoR_ind:SurNoR_ind,:],dims=2),dims=3)
Acc_All = NCor_All ./ NTot_All

Acc_All_Mean_4Parts = mean(Acc_All,dims=1)
Acc_All_Std_4Parts = std(Acc_All,dims=1)/sqrt(Sub_Num)

Acc_All_Mean = mean(Acc_All_Mean_4Parts,dims=5)
Acc_All_Std = sqrt.(mean(Acc_All_Std_4Parts,dims=5).^2 .+ var(Acc_All_Mean_4Parts,dims=5)./4)

println("Total accuracy rate:")
@show Acc_All_Mean[1];
println("SEM of the total accuracy rate:")
@show Acc_All_Std[1];
println("Total number of actions:")
@show sum(NTot_All)/Part_Num

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Reporting the accuracy rates of other exploration strategies
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
OI_Exp_index = [1, 5, 9]
U_Exp_index = [2, 6, 10]

Acc_Matrix_Mean4Parts = mean(Acc_Matrix,dims=5)

println("********************************************")
println("OI:")
println("********************************************")
println("SurNoR vs MB+S+OI:")
@show OneSampleTTest(Acc_Matrix_Mean4Parts[:,1,1,SurNoR_ind,1],Acc_Matrix_Mean4Parts[:,1,1,OI_Exp_index[1],1])
println("--------------------------------------------")
println("SurNoR vs MF+OI:")
@show OneSampleTTest(Acc_Matrix_Mean4Parts[:,1,1,SurNoR_ind,1],Acc_Matrix_Mean4Parts[:,1,1,OI_Exp_index[2],1])
println("--------------------------------------------")
println("SurNoR vs Hyb+S+OI:")
@show OneSampleTTest(Acc_Matrix_Mean4Parts[:,1,1,SurNoR_ind,1],Acc_Matrix_Mean4Parts[:,1,1,OI_Exp_index[3],1])
println("********************************************")
println("OI:")
println("********************************************")
println("SurNoR vs MB+S+U:")
@show OneSampleTTest(Acc_Matrix_Mean4Parts[:,1,1,SurNoR_ind,1],Acc_Matrix_Mean4Parts[:,1,1,U_Exp_index[1],1])
println("--------------------------------------------")
println("SurNoR vs MF+S+U:")
@show OneSampleTTest(Acc_Matrix_Mean4Parts[:,1,1,SurNoR_ind,1],Acc_Matrix_Mean4Parts[:,1,1,U_Exp_index[2],1])
println("--------------------------------------------")
println("SurNoR vs Hyb+S+U:")
@show OneSampleTTest(Acc_Matrix_Mean4Parts[:,1,1,SurNoR_ind,1],Acc_Matrix_Mean4Parts[:,1,1,U_Exp_index[3],1])
println("--------------------------------------------")

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Reporting the difference in log-evidence between SurNoR and control algorithms
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
y = dLogL_Sum_Block[1,1,1,:,:] .+ dLogL_Sum_Block[1,1,2,:,:]
dy = std(y,dims=2)/sqrt(Part_Num)
y = mean(y,dims=2)

println("--------------------------------------------")
println("SurNoR vs. Binary Novelty 1:")
println("Mean:")
@show (y[SurNoR_ind] - y[SurNoR_ind+1])
println("SEM:")
@show sqrt(dy[SurNoR_ind]^2 + dy[SurNoR_ind+1]^2)

println("--------------------------------------------")
println("SurNoR vs. Binary Novelty 2:")
println("Mean:")
@show (y[SurNoR_ind] - y[SurNoR_ind+2])
println("SEM:")
@show sqrt(dy[SurNoR_ind]^2 + dy[SurNoR_ind+2]^2)

println("--------------------------------------------")
println("SurNoR vs. SurNoR with β_N ≧ 0 also in Epi > 1:")
println("Mean:")
@show (y[SurNoR_ind] - y[SurNoR_ind+3])
println("SEM:")
@show sqrt(dy[SurNoR_ind]^2 + dy[SurNoR_ind+3]^2)

# Only 1st episode
y = dLogL_Sum_Block_Epi1
dy = std(y,dims=3)/sqrt(Part_Num)
y = mean(y,dims=3)

println("--------------------------------------------")
println("--------------------------------------------")
println("SurNoR vs. Binary Novelty 1 in B1E1:")
println("Mean:")
@show (y[SurNoR_ind] - y[SurNoR_ind+1])
println("SEM:")
@show sqrt(dy[SurNoR_ind]^2 + dy[SurNoR_ind+1]^2)

println("--------------------------------------------")
println("SurNoR vs. Binary Novelty 2 in B1E1:")
println("Mean:")
@show (y[SurNoR_ind] - y[SurNoR_ind+2])
println("SEM:")
@show sqrt(dy[SurNoR_ind]^2 + dy[SurNoR_ind+2]^2)


text:
using Rankings, Test

@testset "Литвак (1982), p.78" begin
@test  ⦸(ranking"a₂,a₄,a₁,a₃",
         ranking"a₁,a₃∼a₄,a₂",
         ranking"a₂∼a₃,a₄,a₁",
         ranking"a₃,a₂,a₁∼a₄") ==
[
  0   2  0   1;
 -2   0  1  -2;
  0  -1  0  -1;
 -1   2  1   0
]
end

@testset "Bury, Wagner (2008), p.154" begin
@test  ranking"O₂,O₃,O₄,O₁,O₅" ⦸

       ranking"(O₂,O₃),O₄,(O₁,O₅)" ⦸

       ranking"O₂,(O₁,O₃,O₄),O₅" ⦸

       ranking"O₂,(O₁,O₃,O₄,O₅)" ==

[

  0  4   2   2  -2;

 -4  0  -3  -4  -4;

 -2  3   0  -2  -3;

 -2  4   2   0  -3;

  2  4   3   3   0

]
end

@testset "Bredereck (2009), p.21" begin
@test  ranking"a≻b≻c≻d" ⦸

       ranking"a≻b≻d≻c" ⦸

       ranking"a≻c≻b≻d" ⦸

       ranking"b≻a≻c≻d" ⦸

       ranking"b≻a≻d≻c" ==

[

 0  -1  -5  -5;

 1   0  -3  -5;

 5   3   0  -1;

 5   5   1   0

]
end

@testset "Muravyov, Tarakanov (2011)" begin
@test  ranking"A K D N J H E∼L P I B∼O F R M T Q S G C" ⦸

       ranking"T I C A Q F H N M G B O D R P L E K J S" ⦸

       ranking"R N T P∼I C G∼D K O H J Q S∼M B E A F L" ⦸

       ranking"H B M A Q N D O T L I G R S E∼C J F K P" ==

[

 0   0   0  -2  -2  -4  -2  0   0  -2  -2  -4   0  -2  -2  -2  -2  -2  -2   0;

 0   0   0   0  -2  -2   0  4   2   0   0  -2   0   2  -1   0   0  -2  -2   0;

 0   0   0   0  -1  -2   0  0   4  -2  -2   0   0   2   0   0   0   2   0   4;

 2   0   0   0  -4  -2  -1  0   0  -4  -2  -4   0   2  -2  -2   0  -2  -4   0;

 2   2   1   4   0  -2   2  4   2   0   0   1   2   4   2   0   2   2   0   2;

 4   2   2   2   2   0   0  2   4   2   0   0   0   2   2   0   2   0   0   2;

 2   0   0   1  -2   0   0  2   4  -2  -2   0   2   4   0   0   2   0  -2   4;

 0  -4   0   0  -4  -2  -2  0   0  -2   0  -4  -4   0  -2  -2  -2  -2  -4   0;

 0  -2  -4   0  -2  -4  -4  0   0  -2  -2   0  -2   2  -2  -1  -2  -2  -4   2;

 2   0   2   4   0  -2   2  2   2   0   2   0   0   4   2   0   0   2  -2   2;

 2   0   2   2   0   0   2  0   2  -2   0   0   0   2   0   0   0   2  -2   2;

 4   2   0   4  -1   0   0  4   0   0   0   0   2   4   2   0   2   0  -2   2;

 0   0   0   0  -2   0  -2  4   2   0   0  -2   0   2   0   0   0   0  -3   0;

 2  -2  -2  -2  -4  -2  -4  0  -2  -4  -2  -4  -2   0  -4  -4   0  -2  -4  -2;

 2   1   0   2  -2  -2   0  2   2  -2   0  -2   0   4   0   0   0  -2  -4   0;

 2   0   0   2   0   0   0  2   1   0   0   0   0   4   0   0   0   2  -2   2;

 2   0   0   0  -2  -2  -2  2   2   0   0  -2   0   0   0   0   0   0  -4   2;

 2   2  -2   2  -2   0   0  2   2  -2  -2   0   0   2   2  -2   0   0  -4   0;

 2   2   0   4   0   0   2  4   4   2   2   2   3   4   4   2   4   4   0   4;

 0   0  -4   0  -2  -2  -4  0  -2  -2  -2  -2   0   2   0  -2  -2   0  -4   0

]
@test  ranking"AIFDHCBEJG" ⦸

       ranking"BHDAFCIEGJ" ⦸

       ranking"DIEBCJAHGF" ⦸

       ranking"CJHADIFGBE" ⦸

       ranking"ICAJHDGBFE" ==

[

  0  -1   1  -1  -3  -5  -5  -1  -1  -1;

  1   0   1   3  -3  -1  -1   1   3  -1;

 -1  -1   0   1  -3  -1  -5  -1   1  -5;

  1  -3  -1   0  -5  -3  -5   1  -1  -1;

  3   3   3   5   0   3  -1   3   5  -1;

  5   1   1   3  -3   0  -1   3   3   1;

  5   1   5   5   1   1   0   5   5   3;

  1  -1   1  -1  -3  -3  -5   0   1   1;

  1  -3  -1   1  -5  -3  -5  -1   0  -3;

  1   1   5   1   1  -1  -3  -1   3   0

]
@test  ranking"J∼K A∼B∼C∼D∼E∼F∼G∼H∼I∼L∼M" ⦸

       ranking"J A∼B∼C∼D∼E∼F∼G∼H∼I∼K∼L∼M" ⦸

       ranking"J∼K A∼B∼C∼D∼E∼F∼G∼H∼I∼L∼M" ⦸

       ranking"I∼J∼K A∼B∼C∼D∼E∼F∼G∼H∼L∼M" ⦸

       ranking"H∼I∼J∼K L∼A∼B∼C∼D∼E∼F∼G∼M" ⦸

       ranking"H∼I∼J A∼B∼C∼D∼E∼F∼G∼K∼L∼M" ⦸

       ranking"K∼L A∼B∼C∼D∼E∼F∼G∼H∼I∼J∼M" ⦸

       ranking"I∼J∼K A∼B∼C∼D∼E∼F∼G∼H∼L∼M" ⦸

       ranking"B∼C A∼D∼E∼F∼G∼H∼I∼J∼K∼L∼M" ⦸

       ranking"H∼I A∼B∼C∼D∼E∼F∼G∼J∼K∼L∼M" ==

[

  0   1   1   0   0   0   0   3   5  7   6   1   0;

 -1   0   0  -1  -1  -1  -1   2   4  6   5   0  -1;

 -1   0   0  -1  -1  -1  -1   2   4  6   5   0  -1;

  0   1   1   0   0   0   0   3   5  7   6   1   0;

  0   1   1   0   0   0   0   3   5  7   6   1   0;

  0   1   1   0   0   0   0   3   5  7   6   1   0;

  0   1   1   0   0   0   0   3   5  7   6   1   0;

 -3  -2  -2  -3  -3  -3  -3   0   2  4   3  -2  -3;

 -5  -4  -4  -5  -5  -5  -5  -2   0  2   1  -4  -5;

 -7  -6  -6  -7  -7  -7  -7  -4  -2  0  -1  -6  -7;

 -6  -5  -5  -6  -6  -6  -6  -3  -1  1   0  -5  -6;

 -1   0   0  -1  -1  -1  -1   2   4  6   5   0  -1;

  0   1   1   0   0   0   0   3   5  7   6   1   0

]
end

@testset "Farnoud, Milenkovic (2014), p.6427" begin
@test  ranking"1423" ⦸

       ranking"1432" ⦸

       ranking"2314" ⦸

       ranking"4231" ⦸

       ranking"3241" ==

[

  0   1   1  -1;

 -1   0  -1   1;

 -1   1   0   1;

  1  -1  -1   0

]
@test  ranking"54132" ⦸

       ranking"15423" ⦸

       ranking"43512" ⦸

       ranking"13452" ⦸

       ranking"42531" ⦸

       ranking"12534" ⦸

       ranking"24351" ==

[

  0  -3  -1  1   1;

  3   0  -1  3   1;

  1   1   0  3   1;

 -1  -3  -3  0  -1;

 -1  -1  -1  1   0

]

@test  ranking"123" ⦸

       ranking"123" ⦸

       ranking"321" ⦸

       ranking"213" ==

[

 0  0  -2;

 0  0  -2;

 2  2   0

]
end

@testset "Бухарин, Марышев, Дивуева (2014), p.55" begin
@test  ranking"a₃,a₅,a₂,a₄,a₁" ⦸

       ranking"a₃,a₅,a₄,a₂,a₁" ⦸

       ranking"a₃,a₅,a₂∼a₄,a₁" ==

[

  0   3  3   3   3;

 -3   0  3   0   3;

 -3  -3  0  -3  -3;

 -3   0  3   0   3;

 -3  -3  3  -3   0

]

end

@testset "Болтенков, Куваева, Червоненко (2018), p.55" begin
@test  ranking"c≻e≻d≻b≻a" ⦸

       ranking"c≻e≻d≻b≻a" ⦸

       ranking"c≻e≻b≻d≻a" ⦸

       ranking"c≻e≻d≻b≻a" ⦸

       ranking"e≻c≻d≻a≻b" ⦸

       ranking"e≻d≻c≻a≻b" ⦸

       ranking"e≻c≻d≻a≻b" ⦸

       ranking"b≻e≻d≻c≻a" ⦸

       ranking"d≻b≻e≻c≻a" ==

[

  0   3   9   9  9;

 -3   0   5   5  5;

 -9  -5   0  -3  1;

 -9  -5   3   0  7;

 -9  -5  -1  -7  0

]

end

@testset "Muravyov, Baranov, Emelyanova (2019)" begin
@test  ranking"A∼C,B" ⦸

       ranking"C,A∼B" ⦸

       ranking"B∼C,A" ==

[

  0   0  2;

  0   0  2;

 -2  -2  0

]
@test  ranking"A B F D (C E)" ⦸

       ranking"D E A B C F" ⦸

       ranking"B (E A)(C D) F" ⦸

       ranking"(F C) D B (A E)" ⦸

       ranking"C D (B F) E A" ==

[

  0   1  -1  1   1  -1;

 -1   0  -1  1  -3  -2;

  1   1   0  0   0  -2;

 -1  -1   0  0  -3  -1;

 -1   3   0  3   0   1;

  1   2   2  1  -1   0

]

end


text:
line = readlines(a)
line = line[21:end]
nodes = Array{Float64}[]
#read nodes
println("reading nodes...\n")
global i=1
while true
    global i
    l = line[i]
    l = split(l)
    arr = [parse(Float64,k) for k in l]
    if length(arr)==1 & (arr[1]== -1.)
        break
    end
    i = i+2
    push!(nodes,arr)
end
line = line[i+2:end]
#read elements
println("reading elements....\n")
elem = Array{Int64}[]
edge = Array{Int64}[]

i = 1
while true

    global i
    l = line[i]
    l = split(l)
    arr = [parse(Float64,k) for k in l]
    if length(arr)==1 & (arr[1]==-1)
        break
    elseif arr[2]==22
        i = i+3
        l = line[i-1]
        l = split(l)
        arr = [parse(Float64,k) for k in l]
        push!(edge,arr)
    elseif arr[2]==42
        l = line[i+1]
        l = split(l)
        arr = [parse(Float64,k) for k in l]
        push!(elem,arr)
        i = i+2
    end

end
#read boundaries
println("reading boundaries...\n")

function splitToInt(l)
    l = split(l)
    l = [parse(Int64,k) for k in l]
    return l
end
i = i+3 #record 1 of 2467 format
n = splitToInt(line[i])[end] #shows how many elems are to read
line = line[i+2:end]
bc = Array{Int64}[]
global n = div(n,2)+mod(n,2)

while true
global i = 1
global line,n

println("splitting boundaries \n")
while true
    global i,line,n
    if (i>n)
        break
    end
    if tryparse(Float64, split(line[i])[1]) == nothing
        break
    end
    arr = splitToInt(line[i])
    push!(bc,[arr[2]])
    if length(arr)>5
        push!(bc,[arr[6]])
    end
    i = i+1
end
    push!(bc,[-1])
    arr = splitToInt(line[i])
    if arr[1]==-1
        break
    end
    n = arr[end]
    n = div(n,2)+mod(n,2)
    line = line[i+2:end]
end

println("making the arrays \n")
using Flux:batch
nodes = Transpose(batch(nodes))
edge = Transpose(batch(edge))
elem = Transpose(batch(elem))
elem = elem[:,[1 3 5 2 4 6]][:,1,:]
bc = Transpose(batch(bc))

max_stars_repo_name: billyon/first

max_issues_repo_name: billyon/first

max_forks_repo_name: billyon/first


text:
"""
Basic initialization
"""
function initialize(user_props;grid_dims=(3,3),seed=0)
    # pulses = []
    # for i in 1:args[:q]
    #     push!(pulses,[])
    # end

    # Global model props
    default_props = Dict(
        :ticks => 0,# # time unit
        :pkt_id => 0,
        :amsg_id =>0,
        :ofmsg_id=>0,
        :capacity_factor=>1.2, #default capacity factor of pkts processed per tick (.2 extra is to have always room for management msgs/pkts.)
        :mapping_ctl_ntw => Dict{Int64,Int64}(), # mapping between (Ctl) Agent and SimNE
        :mapping_ntw_sne => Dict{Int64,Int64}(), #mapping btwn the underlying network and the corresponding simNE agent
        :ntw_links_msgs=>Dict{Tuple{Int,Int},Vector{Vector{OFMessage}}}(),
        :ntw_links_delays =>Dict{Tuple{Int,Int},Int}(),
        :ntw_links_capacity=>Dict{Tuple{Int,Int},Int}(),
        :state_trj => Vector{ModelState}(),
        :base_ntw_graph => user_props[:ntw_graph],
        :dropped_nodes => Vector{Tuple{Int,Int}}()
    )
    #For G6:
    #prob_eq_queries_cycle: 0.2
    #clear_cache_graph_freq: 30 or 15
    #query_cycle: 20 or 10
    # requires redundancy of queries + longer time graph cache (> query cycle)
    #For G7/G4:
    #prob_eq_queries_cycle: 0.82
    #clear_cache_graph_freq: 10
    #query_cycle: 10

    Random.seed!(seed)
    props = merge(default_props,user_props)
    #space = GridSpace(grid_dims, moore=true)
    space = GraphSpace(props[:ntw_graph])
    agent_types = Union{SimNE,Agent}
    model = ABM(agent_types, space;
    scheduler = Agents.Schedulers.randomly #random_activation
    , properties = props)
    init_model!(model)
    #create
    create_agents!(model)
    model
end

"""
    Create simulated asset agents
"""
function create_sim_asset_agents!(model)
    a_params = Dict()

    # create SimNE
    for i in 1:nv(model.properties[:ntw_graph])
        #next_fire = rand(0:0.2:model.:Τ)
        id = nextid(model)

        deterioration = id in model.init_sne_params.ids ? model.init_sne_params.deterioration[first(indexin(id,model.init_sne_params.ids))] : model.deterioration

        mnt = @match model.mnt_policy begin
            1 => MaintenanceInfoPreventive(deterioration,model)
            2 => MaintenanceInfoPredictive(deterioration,model)
            _ => MaintenanceInfoCorrective(deterioration,model)
        end

        capacity_factor =  id in model.init_sne_params.ids ? model.init_sne_params.capacity_factor[first(indexin(id,model.init_sne_params.ids))] : model.capacity_factor
        max_q_capacity =  capacity_factor * model.pkt_per_tick
        a_params[:pkt_per_tick] = max_q_capacity
        a = add_agent_pos!(
                SimNE(id,i,a_params,max_q_capacity,mnt),model
            )

        set_prop!(model.ntw_graph, id, :eid, id )
        # create initial mapping btwn network and SimNE
        model.mapping_ntw_sne[i] = i
    end

    set_indexing_prop!(model.ntw_graph, :eid)
end

"""
    Create control agents
"""
function create_control_agents!(model::ABM)

    a_params::Dict{Symbol,Any} = Dict(:pkt_per_tick=>model.pkt_per_tick)
    mnt::MaintenanceInfo = @match model.mnt_policy begin
                1 => MaintenanceInfoPreventive(model)
                2 => MaintenanceInfoPredictive(model)
                _ => MaintenanceInfoCorrective(model)
            end

    if model.ctrl_model == CENTRALISED
        id = nextid(model)
        a = add_agent_pos!(
                Agent(id,1,mnt,Array{Float64}(undef,0,nv(model.properties[:ntw_graph])),a_params),model
            )
        #set_props!(model.ctl_graph, 1, Dict(:eid => 1, :aid => id) )
         ##assign controller to SimNE
        for j in 1:nv(model.properties[:ntw_graph])
            set_control_agent!(j,id,model)
        end
    else
        for i in 1:nv(model.properties[:ctl_graph])
            id = nextid(model)
            a = add_agent_pos!(
                    Agent(id,i,mnt,Array{Int64}(undef,0,1),a_params),model
                )
            set_props!(model.ctl_graph, i, Dict(:eid => i, :aid => id) )
            ##assign controller to SimNE
            #one-to-one mapping
            set_control_agent!(i,id,model)
        end
        set_indexing_prop!(model.ctl_graph, :aid)
    end

end

"""
Simplest create agents
"""
function create_agents!(model)
    create_sim_asset_agents!(model)

    create_control_agents!(model)

    log_info("Total agents created: $(length(allagents(model))), seed: $(model.seed)")

    init_agents!(model)
end

"""
    It advance model one step. Agents are progressed here to ensure that one action is perform in
    all agents before the next action is performed in all agents, instead of having one action
    executed in some agents and then some start with second action.
"""
function model_step!(model)
    init_step_state!(model)
    # log_info(model.ticks,"==========")
    log_info(".")
    for a in allagents(model)
        init_state!(a)
    end
    generate_traffic!(model)
    for e in edges(model.ntw_graph)
        ntw_link_step!((e.src,e.dst),model)
    end
    ctl_links_step!(model)

    # Run controllers first
    for a in allagents(model)
        if typeof(a) == Agent
            # kn = [ a.ntw_graph[i,:eid] for i=1:nv(a.ntw_graph) ]
            do_agent_step!(a,model)
        end
    end
    #Then run SimNEs
    for a in allagents(model)
        if typeof(a) == SimNE
            do_agent_step!(a,model)
        end
    end
    for a in allagents(model)
        pending_pkt_handler(a,model)
        clear_cache!(a,model)
        calculate_metrics_step!(a,model)
    end
    trigger_random_node_drops!(model)
end

"""
    Progress control agents one step.
"""
function agent_step!(a::Agent,model)
    #placeholder
end

"""
    Progress SimNE one step.
"""
function agent_step!(a::SimNE,model)
    #placeholder
end

"""
Simple run model function
"""
function run_model(n,properties; agent_data, model_data)
    seed = properties[:seed]
    model = initialize(properties;seed)

    # agent_shape(a::Agent) = :square
    agent_size(a) = 7
    # agent_size(a::Agent) = 7

    plot_scheduler = model.scheduler#by_type((SimNE,Agent), false)#

    plotkwargs = (
            # ac = agent_color, am = agent_shape, as = agent_size,
            # scheduler = plot_scheduler, aspect_ratio = 1
            # #Graph space properties
            # , method = :circular, node_size = 0.2
            # ,
            size = (1000, 600),
            showaxis = false
    )

    #p = plotabm(model; plotkwargs...)
    df = init_agent_dataframe(model,agent_data)
    df_m = init_model_dataframe(model,model_data)

        if model.animation

        anim = @animate for i in 0:n
                p = plotabm_networks(model; plotkwargs...)
                title!(p, "step $(i)",position=(10,1))
                #annotate!((1,1,Plots.text("step $(i)", 11, :black, :center)))
                step!(model, agent_step!,model_step!)
                collect_agent_data!(df, model, agent_data, i)
                collect_model_data!(df_m, model, model_data, i)
            end
        plot_label = model.run_label * "_anim"
        check_create_dir!(model.plots_dir)
        gif(anim, model.plots_dir * plot_label * ".gif", fps = 5)
    else
        for i in 0:n
            step!(model, agent_step!,model_step!)
            collect_agent_data!(df, model, agent_data, i)
            collect_model_data!(df_m, model, model_data, i)
        end
    end
    df, df_m
end

function init_agents!(model)
    # for a in sort(allagents(model))
    #     init_agent(a,model)
    # end

    #force ordered start

    ids = [a.id for a in allagents(model)]

    for id in sort(ids)
        init_agent!(getindex(model,id),model)
    end
end

"""
Initialise control agents
"""
function init_agent!(a::Agent,model)

    if model.ctrl_model != GraphModel(1)
        #Calculate sub graphs and init msg channels among agents
        nodes = [get_controlled_assets(a.id,model)...]
        sub_g = get_subgraph(model.ntw_graph,nodes,:eid)

        nodes = [a.id]
        ctl_sub_g = get_subgraph(model.ctl_graph,nodes,:aid)
        a.ntw_graph = sub_g
        a.base_ntw_graph = sub_g
        a.ctl_graph = ctl_sub_g
        a.params[:delay_ctl_link] = 2 # 1: no delay
        #Init vector of msgs
        a.msgs_links = init_array_vectors(AGMessage,a.params[:delay_ctl_link],degree(a.ctl_graph,to_local_vertex(a.ctl_graph,a.id,:aid)))
    else
        a.ntw_graph = model.ntw_graph
        a.base_ntw_graph = model.ntw_graph
    end

    ##
    if a.maintenance.policy !=  CorrectiveM
        schedule_event!(a,CTL_Event(4),a.maintenance.predictive_freq,Array{Int64,1}())
        if a.maintenance.policy ==  PredictiveM
            if Symbol("py_integration") in keys(model.properties)
                #conversion to py
                ajm_py =  model.py_integration.np.matrix(adjacency_matrix(a.ntw_graph))
                model.py_integration.opt_init.optimisation_initialisation( ajm_py,
                model.traffic_dist_params
                #[1,0.05]
                , model.mnt_bc_cost, model.mnt_bc_duration, model.mnt_wc_cost, model.mnt_wc_duration)
            end
        end
    end
end

function label_paths(time::Int64,paths::Array{YenState{Float64,Int64},1})
    npaths = Dict()
    if !isempty(paths)
        for path in paths
            if !isempty(path.paths)
                #obtain only the first path
                npaths[(first(first(path.paths)),last(first(path.paths)))] = [(time,first(path.dists),first(path.paths))]
            end
        end
    end
    return npaths
end

function label_path(path)

    return (first(path),last(path),path)
end

"""
Initialise sne agents
"""
function init_agent!(sne::SimNE,model)
    #print("Initialisation of SimNE agent $(a.id)")
    nbs = all_neighbors(model.ntw_graph,get_address(sne.id,model.ntw_graph))

    push_ep_entry!(sne,(0,"h$(sne.id)")) # link to a host of the same id

    for i in 1:size(nbs,1)
        push_ep_entry!(sne,(i,"s$(nbs[i])"))
        #install flows to host of the same id
        # install_flow!(Flow(sne.id,MRule("$i","$(nbs[i])","$(sne.id)"),[0],OFS_Action(1)),sne,model)
    end

    init_condition!(sne,model)
    init_maintenance!(sne,model)
end

"""
    Next agent message's id
"""
function next_amid!(model)
    model.amsg_id += 1
    return model.amsg_id
end

"""
    Next Open Flow message id
"""
function next_ofmid!(model)
    model.ofmsg_id += 1
    return model.ofmsg_id
end

"""
    Initialise state of model for a given step
"""
function init_step_state!(m::ABM)
    m.ticks += 1
    new_state = ModelState(m.ticks)#deepcopy(get_state(m))
    #for e in edges(m.ntw_graph)
    #    new_state.links_load[(e.src,e.dst)] = 0
           # read from links according to delay
   # put msgs in destination nodes for those where delay is over

        #m.ntw_links[k]
    #end
    push!(m.state_trj,new_state)
end

function get_state(m::ABM)
    return last(m.state_trj)
end

function set_active_flows(m::ABM,af::Vector{Tuple{Int,Int,Flow_Type}})
    current = pop!(m.state_trj)
    current.active_flows = af
    push!(m.state_trj,current_logger)
end

function set_state!(m::ABM,new_state::ModelState)
    pop!(m.state_trj)
    push!(m.state_trj,new_state)
end

function init_model!(m::ABM)
    for e in edges(m.ntw_graph)
        #all delays equal initially
        m.ntw_links_delays[(e.src,e.dst)] = 1
        #set capacity for each link, checking if any init param was passed for the links
        if (e.src,e.dst) in m.init_link_params.ids
            idx = first([i for i=1:length(m.init_link_params.ids) if m.init_link_params.ids[i] == (e.src,e.dst) ])
            m.ntw_links_capacity[(e.src,e.dst)] = m.init_link_params.capacities[idx]
        else
            m.ntw_links_capacity[(e.src,e.dst)] = m.link_capacity
        end
    end
end

function init_link_msg!(l::Tuple{Int,Int},m::ABM)
    link_queue = Vector{Vector{OFMessage}}()
    for q=1:m.ntw_links_delays[l]
        push!(link_queue,Vector{OFMessage}())
    end
    m.ntw_links_msgs[l] = link_queue
end

"""
This function progresses one step each link of the sne network.
"""
function ntw_link_step!(l::Tuple{Int,Int},model)
    # Done for each registered link
    if haskey(model.ntw_links_msgs,l)

        msgs = model.ntw_links_msgs[l]
        to_deliver = first(msgs)
        in_pkt_count = 0

        if !isempty(to_deliver)
            if length(msgs) > 1
                msgs = msgs[begin+1:end]
                push!(msgs,Vector{OFMessage}())
            else
                msgs = [ Vector{OFMessage}() ]
            end
            model.ntw_links_msgs[l] = msgs

            for msg in to_deliver
                #Does it need to check address?, I don't think so
                dst = msg.dpid == l[1] ? getindex(model,l[2]) : getindex(model,l[1])

                if length(dst.queue.data) < dst.queue.sz_max - Int(round(dst.queue.sz_max * .2))

                     #TODO make % a parameter. This the gap for OF control msgs.

                    put!(dst.queue,msg)
                    in_pkt_count = get_state(dst).in_pkt + 1
                    set_in_pkt!(dst,in_pkt_count)
                else
                    #if destination queue is full the msg/pkt is dropped
                    drop_packet!(dst)
                end
            end
        end
    end

end

"""
    Initialize an array of dimensions d1 x d2 that contains vectors of type T
"""
function init_array_vectors(T,d1,d2)
    arr = d2 >= 1 ? Array{Vector{T}}(undef,d1,d2) : Array{Vector{T}}(undef,d1)
    [ arr[i,j] = Vector{T}() for i=1:size(arr,1) for j=1:size(arr,2)]
    return arr
end

function ctl_links_step!(model)
    if model.ctrl_model != GraphModel(1)
        ctl_ags = filter(a->typeof(a) == Agent,Set(allagents(model)))
        ctl_link_step!.(ctl_ags)
    end

end

function ctl_link_step!(a::Agent)
    #Merge msgs from all senders to be processed
    a.msgs_in = vcat(a.msgs_links[1,:]...)
    a.msgs_links[1,:] = init_array_vectors(AGMessage,size(a.msgs_links,2),1)
    a.msgs_links = circshift(a.msgs_links,1)
end

function get_state_trj(m::ABM)::Vector{ModelState}
    return m.state_trj
end

"""
Return total of control messages exchanged by agents
"""
function get_ag_msg(model)
    return cumsum(sum.(eachcol([ [ s.in_ag_msg for s in a.state_trj ] for a in allagents(model) if typeof(a) == Agent ]))...)
end

function is_active_flow(f::Tuple{Int,Int},model)
    v = !isempty(filter(af->(af[1],af[2]) == f || (af[2],af[1]) == f  ,get_state(model).active_flows))
    return v
end


text:
using EffectiveWaves, Test
using LinearAlgebra

@testset "Check material definitions" begin
    include("../src/materials.jl")
    @test true
end

@testset "Examples from: Reflection from multi-species.., Proc.R.Soc.(2018)" begin

    include("../docs/src/examples/concrete/concrete_species.jl")
    include("../docs/src/examples/concrete/concrete_species_volfrac.jl")
    # include("../docs/src/examples/concrete/concrete_species_large-freq.jl") # takes longer

    # Takes too long
    include("../docs/src/examples/emulsion/fluid_species.jl")
    include("../docs/src/examples/emulsion/fluid_species_volfrac.jl")
    # include("../docs/src/examples/emulsion/fluid_species_large-freq.jl") # takes longer

    @test true
end

include("specialfunctions.jl")

include("types_constructors.jl")

include("strong_low_freq_effective.jl")
include("high_frequency_effective.jl")

include("large_vol_low_freq_effective.jl")
# A heavy test..
include("weak_scatterers_effective.jl")

# test does not run on Julia version < 0.7 due to differences in Optim versions
include("path_mesh_wavenumbers.jl")

include("numerical_integration.jl")
include("integrated_reflection.jl")
include("average_integrand_kernel.jl")

include("match_wave.jl")
include("wiener-hopf-reflection.jl")
include("wiener-hopf-wave.jl")


text:
export expected,
       m_variance,
       apply!

"""
    expected(Op::Mpo{T}, psi::Mps{T}) where T<:Number

Compute expectation value < psi|Op|psi >.
"""
function expected(Op::Mpo{T}, psi::Mps{T}) where T<:Number
    L = ones(T, 1, 1, 1)
    for i=1:psi.L
        L = prop_right3(L, psi.M[i], Op.W[i], psi.M[i])
    end
    return L[1, 1, 1]
end

"""
    norm_of_apply(Op::Mpo{T}, psi::Mps{T}) where T<:Number

Compute norm of the operator acting on the state: || Op|psi > ||^2.
"""
function norm_of_apply(Op::Mpo{T}, psi::Mps{T}) where T<:Number
    L = ones(T, 1, 1, 1, 1)
    for i=1:psi.L
        hOp = permutedims(conj(Op.W[i]), (1, 3, 2, 4))
        L = prop_right4(L, psi.M[i], Op.W[i], hOp, psi.M[i])
    end
    return real(L[1, 1, 1, 1])
end

"""
    m_variance(Op::Mpo{T}, psi::Mps{T}) where T<:Number

Measure variance: <psi|Op^2|psi> - |<psi|Op|psi>|^2, return a real number.
"""
function m_variance(Op::Mpo{T}, psi::Mps{T}) where T<:Number
    return norm_of_apply(Op, psi) - abs2(expected(Op, psi))
end

"""
    apply!(Op::Mpo{T}, psi::Mps{T}) where T<:Number

Apply an Mpo to an Mps: Op|psi>. Keep new state unnormalized.
"""
function apply!(Op::Mpo{T}, psi::Mps{T}) where T<:Number
    psi.L == Op.L || throw("Mpo and Mps have different lengths.")
    psi.d == Op.d || throw("Mpo and Mps have different physical dimension.")
    for i=1:psi.L
        @tensor Mi[l1, l2, s2, r1, r2] := psi.M[i][l1, s1, r1]*Op.W[i][l2, s1, s2, r2]
        psi.M[i] = reshape(Mi, size(psi.M[i], 1)*size(Op.W[i], 1), psi.d,
                               size(psi.M[i], 3)*size(Op.W[i], 4))
    end
    return psi
end


text:
module ElectronDisplay

export electrondisplay

using Electron, Base64, Markdown

import IteratorInterfaceExtensions, TableTraits, TableShowUtils

Base.@kwdef mutable struct ElectronDisplayConfig
    showable = electron_showable
    single_window::Bool = false
    focus::Bool = true
    max_json_bytes::Int = 2^20
end

"""
    setconfig(config; kwargs...)

Update a copy of `config` based on `kwargs`.
"""
setconfig(
    config::ElectronDisplayConfig;
    showable = config.showable,
    single_window::Bool = config.single_window,
    focus::Bool = config.focus,
    max_json_bytes::Int = config.max_json_bytes,
) =
    ElectronDisplayConfig(
        showable = showable,
        single_window = single_window,
        focus = focus,
        max_json_bytes = max_json_bytes,
    )

struct ElectronDisplayType <: Base.AbstractDisplay
    config::ElectronDisplayConfig
end

ElectronDisplayType() = ElectronDisplayType(CONFIG)
newdisplay(; config...) = ElectronDisplayType(setconfig(CONFIG; config...))

electron_showable(m, x) =
    m ∉ ("application/vnd.dataresource+json", "text/html", "text/markdown") &&
    showable(m, x)

"""
    ElectronDisplay.CONFIG

Configuration for ElectronDisplay.

* `showable`: A callable with signature `showable(mime::String,
  x::Any) :: Bool`.  This determines if object `x` is displayed by
  `ElectronDisplay`.  Default is to return `false` if `mime` is
  `"text/html"` or `"text/markdown"` and otherwise fallbacks to
  `Base.showable(mime, x)`.

* `single_window::Bool = false`: If `true`, reuse existing window for
  displaying a new content.  If `false` (default), create a new window
  for each display.

* `focus::Bool = true`: Focus the Electron window on `display` if `true`
  (default).

* `max_json_bytes::Int = $(ElectronDisplayConfig().max_json_bytes)`:
  Maximum size in bytes for which JSON representation is used.  Otherwise,
  convert visualization locally in a binary form before sending it to the
  Electron display.  Currently only Vega and Vega-Lite support this.
"""
const CONFIG = ElectronDisplayConfig()

const _window = Ref{Window}()

function _getglobalwindow()
    if !(isdefined(_window, 1) && _window[].exists)
        _window[] = Electron.Window(
            URI("about:blank"),
            options=Dict("webPreferences" => Dict("webSecurity" => false)))
    end
    return _window[]
end

function displayhtml(d::ElectronDisplayType, payload; options::Dict=Dict{String,Any}())
    if d.config.single_window
        w = _getglobalwindow()
        load(w, payload)
        showfun = get(options, "show", d.config.focus) ? "show" : "showInactive"
        run(w.app, "BrowserWindow.fromId($(w.id)).$showfun()")
        return w
    else
        options = Dict{String,Any}(options)
        get!(options, "show", d.config.focus)
        return Electron.Window(payload; options=options)
    end
end

displayhtmlbody(d::ElectronDisplayType, payload) =
    displayhtml(d, string(
        """
        <!doctype html>

        <meta name="viewport" content="width=device-width, initial-scale=1">
        <link rel="stylesheet" href="file:///$(asset("github-markdown-css", "github-markdown.css"))">

        <link rel="stylesheet" href="file://$(asset("katex-0.11.1", "katex.min.css"))" integrity="sha384-zB1R0rpPzHqg7Kpt0Aljp8JPLqbXI3bhnPWROx27a9N0Ll6ZP/+DiW/UqRcLbRjq" crossorigin="anonymous">

        """,
        payload,
        """

        """,
    ))

function Base.display(d::ElectronDisplayType, ::MIME{Symbol("text/html")}, x)
    html_page = repr("text/html", x)
    if occursin(r"<html\b"i, html_page)
        # Detect if object `x` rendered itself as a "standalone" HTML page.
        # If so, display it as-is:
        displayhtml(d, html_page)
    else
        # Otherwise, i.e., if `x` only produced an HTML fragment, apply
        # our default CSS to it:
        displayhtmlbody(d, html_page)
    end
end

Base.displayable(d::ElectronDisplayType, ::MIME{Symbol("text/html")}) = true

Base.display(d::ElectronDisplayType, ::MIME{Symbol("text/markdown")}, x) =
    displayhtmlbody(d, repr("text/html", asmarkdown(x)))

asmarkdown(x::Markdown.MD) = x
asmarkdown(x) = Markdown.parse(repr("text/markdown", x))

Base.displayable(d::ElectronDisplayType, ::MIME{Symbol("text/markdown")}) = true

function Base.display(d::ElectronDisplayType, ::MIME{Symbol("image/png")}, x)
    img = stringmime(MIME("image/png"), x)

    payload = string("")

    displayhtml(d, payload)
end

Base.displayable(d::ElectronDisplayType, ::MIME{Symbol("image/png")}) = true

function Base.display(d::ElectronDisplayType, ::MIME{Symbol("image/svg+xml")}, x)
    payload = stringmime(MIME("image/svg+xml"), x)

    displayhtml(d, payload)
end

Base.displayable(d::ElectronDisplayType, ::MIME{Symbol("image/svg+xml")}) = true

asset(url...) = replace(normpath(joinpath(@__DIR__, "..", "assets", url...)), "\\" => "/")

function Base.display(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.plotly.v1+json")}, x)
    payload = stringmime(MIME("application/vnd.plotly.v1+json"), x)

    html_page = """

    """

    displayhtml(d, html_page, options=Dict("webPreferences" => Dict("webSecurity" => false)))
end

Base.displayable(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.plotly.v1+json")}) = true

function Base.display(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.dataresource+json")}, x)
    payload = stringmime(MIME("application/vnd.dataresource+json"), x)

    html_page = """

        <link rel="stylesheet" href="file://$(asset("ag-grid", "ag-grid.css"))">
        <link rel="stylesheet" href="file://$(asset("ag-grid", "ag-theme-balham.css"))">

    """

    displayhtml(d, html_page, options=Dict("webPreferences" => Dict("webSecurity" => false)))
end

Base.displayable(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.dataresource+json")}) = true

function Base.display(d::ElectronDisplayType, x)
    showable = d.config.showable
    if showable("application/vnd.vegalite.v3+json", x)
        display(d,MIME("application/vnd.vegalite.v3+json"), x)
    elseif showable("application/vnd.vegalite.v2+json", x)
        display(d,MIME("application/vnd.vegalite.v2+json"), x)
    elseif showable("application/vnd.vega.v5+json", x)
        display(d,MIME("application/vnd.vega.v5+json"), x)
    elseif showable("application/vnd.vega.v4+json", x)
        display(d,MIME("application/vnd.vega.v4+json"), x)
    elseif showable("application/vnd.vega.v3+json", x)
        display(d,MIME("application/vnd.vega.v3+json"), x)
    elseif showable("application/vnd.plotly.v1+json", x)
        display(d,MIME("application/vnd.plotly.v1+json"), x)
    elseif showable("application/vnd.dataresource+json", x)
        display(d, "application/vnd.dataresource+json", x)
    elseif showable("image/svg+xml", x)
        display(d,"image/svg+xml", x)
    elseif showable("image/png", x)
        display(d,"image/png", x)
    elseif showable("text/html", x)
        display(d, "text/html", x)
    elseif showable("text/markdown", x)
        display(d, "text/markdown", x)
    else
        throw(MethodError(Base.display,(d,x)))
    end
end

"""
    electrondisplay([mime,] x; config...)

Show `x` in Electron window.  Use MIME `mime` if specified.  The keyword
arguments can be used to override [`ElectronDisplay.CONFIG`](@ref) without
mutating it.

# Examples
```julia
electrondisplay(@doc reduce; single_window=true, focus=false)
```
"""
electrondisplay(mime, x; config...) = display(newdisplay(; config...), mime, x)

struct DataresourceTableTraitsWrapper{T}
    source::T
end

function Base.show(io::IO, ::MIME"application/vnd.dataresource+json", source::DataresourceTableTraitsWrapper)
    TableShowUtils.printdataresource(io, IteratorInterfaceExtensions.getiterator(source.source))
end

Base.showable(::MIME"application/vnd.dataresource+json", dt::DataresourceTableTraitsWrapper) = true

struct CachedDataResourceString
    content::String
end

Base.show(io::IO, ::MIME"application/vnd.dataresource+json", source::CachedDataResourceString) = print(io, source.content)

Base.showable(::MIME"application/vnd.dataresource+json", dt::CachedDataResourceString) = true

function electrondisplay(x; config...)
    d = newdisplay(; showable=showable, config...)
    if TableTraits.isiterabletable(x)!==false
        if showable("application/vnd.dataresource+json", x)
            display(d, x)
        elseif TableTraits.isiterabletable(x)===true
            display(d, DataresourceTableTraitsWrapper(x))
        else
            try
                buffer = IOBuffer()
                TableShowUtils.printdataresource(buffer, IteratorInterfaceExtensions.getiterator(x))

                buffer_asstring = CachedDataResourceString(String(take!(buffer)))
                display(d, buffer_asstring)
            catch err
                display(d, x)
            end
        end
    else
        display(d, x)
    end
end

include("vega.jl")

function __init__()
    Base.Multimedia.pushdisplay(ElectronDisplayType())
end

end # module


text:
## --- Read Heuret traces

    # using Images, ImageView
    using FileIO: load
    using Images, Plots, Compat
    using StatGeochem

    img = load("TrenchSedimentThickness.tif")

    sedcolor = [RGB{N0f8}(0.208,0.325,0.643)
        RGB{N0f8}(0.29,0.608,0.835)
        RGB{N0f8}(0.063,0.722,0.925)
        RGB{N0f8}(0.388,0.78,0.78)
        RGB{N0f8}(0.631,0.831,0.651)
        RGB{N0f8}(0.776,0.878,0.639)
        RGB{N0f8}(0.976,0.929,0.51)
        RGB{N0f8}(1.0,0.792,0.392)
        RGB{N0f8}(0.98,0.667,0.29)
        RGB{N0f8}(0.949,0.4,0.404)
        RGB{N0f8}(0.937,0.224,0.369)
    ]

    # Trench sediment thickness from Heuret et al, 2012
    sedthick = [ 0.2
                 0.4
                 0.5
                 0.7
                 0.9
                 1.5
                 2.5
                 3.5
                 4.5
                 6.5
                 7.5]

    latimg = repmat(collect(linspace(-90,90,1125)),1,2250)
    lonimg = repmat(collect(linspace(-180,180,2250))',1125,1)

    t = img .== sedcolor[1]
    latlist = latimg[t]
    lonlist = lonimg[t]
    sedthicklist = ones(sum(t))*sedthick[1]
    for i=2:length(sedcolor)
        t = img .== sedcolor[i]
        latlist = vcat(latlist, latimg[t])
        lonlist = vcat(lonlist, lonimg[t])
        sedthicklist = vcat(sedthicklist, ones(sum(t))*sedthick[i])
    end

## --- Produce histogram of sed thickness vs trench length

    using StatsBase

    # Resample, to account uncertainty, weighting to account for change in length with latitude
    # Assume uncertainty is of the order of half the thickness
    trenchrepresentative = bsresample(sedthicklist,sedthicklist/8,100000,cos.(latlist))

    # Collect histogram
    binedges = 0:8
    trenchfilldist = fit(Histogram, trenchrepresentative[:,1], binedges).weights

    # Plot results
    trenchfilldist *= 65000/length(trenchrepresentative) # Normalize global trench length of about 65000 km
    h = plot(cntr(binedges),trenchfilldist,seriestype=:bar,bar_width=1,label="", framestyle=:box)
    plot!(h,xlims=(0,8),ylims=(0,4E4),xlabel="Trench sediment thickness (km)",ylabel="Trench length (km)", markerstrokecolor=:auto)
    savefig(h,"TrenchFillingExtent.pdf")
    display(h)

## --- Read eHf data with lat-lon constraints

    data = readdlm("Bataille2017S1.csv",'|')
    data = elementify(data)

    lat = data["latitude"]
    lon = data["longitude"]
    sed = Array{Float64}(size(lat))
    sed_sigma = Array{Float64}(size(lat))
    distance = Array{Float64}(size(lat))
    for i=1:length(lat)
        dist = arcdistance(lat[i],lon[i],latlist,lonlist)
        k = argmin(dist)
        sed[i] = sedthicklist[k]
        distance[i] = dist[k]
        # Assume uncertainty is of the order of half the thickness
        sed_sigma[i] = sedthicklist[k]/2
    end

    data["eHf"] = eHf(data["176Hf_177Hf_normalized"], data["176Lu_177Hf"], data["Age_UPb"])

## --- Correlate each zircon

    # Only look at samples younger than 100 Ma and within five arc degress of an arc
    t = (data["Age_UPb"] .< 100) .& (distance .< 5) .& (data["eHf"] .> -25)

    # Clean up data
    LuHf = data["176Lu_177Hf"][t]
    LuHfSigma = data["176Lu_177Hf±"][t]
    LuHfSigma[isnan.(LuHfSigma)] = nanmean(LuHfSigma)
    HfHf = data["176Hf_177Hf_normalized"][t]
    HfHfSigma = data["176Hf_177Hf±"][t]
    HfHfSigma[isnan.(HfHfSigma)] = nanmean(HfHfSigma)
    AgeUPb = data["Age_UPb"][t]
    Age_UPb_Sigma = 0.05*data["Age_UPb"][t]

    # Resample and plot
    # (c,m,el,eu) = bin_bsr_means(sed[t],data["eHf"][t],0,6,6,sed_sigma[t]/2,1000)
    (c,m,el,eu) = bin_bsr_eHf(sed[t],HfHf, LuHf, AgeUPb, 0,6,6, sed_sigma[t],HfHfSigma,LuHfSigma,Age_UPb_Sigma,10000)
    h = plot(c, m, yerror=(el,eu), seriestype=:scatter, label="", framestyle=:box, color=:darkred,markerstrokecolor=:auto)
    plot!(h, xlims=(0,6), xlabel="Present-day trench sediment thickness (km)", ylabel="Zircon epsilon Hf")
    savefig(h,"SedimentSubductionEHfCorrelation.pdf")
    display(h)

## ---

max_stars_repo_name: brenhinkeller/GreatUnconformity

max_issues_repo_name: brenhinkeller/GreatUnconformity

max_forks_repo_name: brenhinkeller/GreatUnconformity