text:
module Expand

using Genie.Renderer.Html, Genie.Renderer, Stipple, StippleUI, StippleUI.API

@reactive mutable struct ExpansionModel <: ReactiveModel
    expanded::R{Bool} = true
end

function ui(model)
    layout(
        list(
            bordered = true,
            class = "rounded-borders",
            [
                expansionitem(
                    @bind(:expanded),
                    default_opened = true,
                    expandseparator = true,
                    icon = "perm_identity",
                    label = "Account settings",
                    caption = "John Doe",
                    [
                        card([
                            cardsection(
                                "Lorem ipsum dolor sit amet, consectetur adipisicing elit. Quidem, eius reprehenderit eos corrupti
                    commodi magni quaerat ex numquam, dolorum officiis modi facere maiores architecto suscipit iste
                    eveniet doloribus ullam aliquid.",
                            ),
                        ]),
                    ],
                ),
                expansionitem(
                    @bind(:expanded),
                    default_opened = true,
                    expandseparator = true,
                    icon = "signal_wifi_off",
                    label = "Wifi settings",
                    [
                        card([
                            cardsection(
                                "Lorem ipsum dolor sit amet, consectetur adipisicing elit. Quidem, eius reprehenderit eos corrupti
                    commodi magni quaerat ex numquam, dolorum officiis modi facere maiores architecto suscipit iste
                    eveniet doloribus ullam aliquid.",
                            ),
                        ]),
                    ],
                ),
                expansionitem(
                    @bind(:expanded),
                    default_opened = true,
                    expandseparator = true,
                    icon = "drafts",
                    label = "Drafts",
                    [
                        card([
                            cardsection(
                                "Lorem ipsum dolor sit amet, consectetur adipisicing elit. Quidem, eius reprehenderit eos corrupti
                    commodi magni quaerat ex numquam, dolorum officiis modi facere maiores architecto suscipit iste
                    eveniet doloribus ullam aliquid.",
                            ),
                        ]),
                    ],
                ),
            ],
        ),
    )
end

function handlers(model)
    on(model.isready) do _
        push!(model)
    end
    model
end

expansionmodel = handlers(Stipple.init(ExpansionModel; transport = Genie.WebThreads))

route("/") do
    html(ui(expansionmodel), context = @__MODULE__)
end

isrunning(:webserver) || up()
# Genie.config.log_requests = false
end #module

max_stars_repo_name: michaelfliegner/BitempWebApp

max_issues_repo_name: michaelfliegner/BitempWebApp


text:
"""
# Freeze.jl

A minimal package to expose an interface for reasoning about and controlling the mutability
of indexable containers, via:

 * `issettable(x)` returns a `Bool` indicating whether `setindex!` works on `x`.
 * `freeze(x)` returns a container with the values contained within `x` that cannot be mutated via `setindex!`.
 * `thaw(x)` is the reverse of `freeze`, returning a container with the values contained with `x` that can be mutated via `setindex!`.

In the general case, the above might require copying the contents of `x`, but in practice we
expect the Julia compiler to be good at reasoning about the lifetime of objects via escape
analysis.
"""
module Freeze

using SparseArrays

export issettable, freeze, thaw
export FrozenArray

"""
    issettable(container)

Return `true` if `setindex!` works on `container`, or `false` otherwise.
"""
function issettable end

"""
    freeze(container)

Freeze an indexable container to make it immutable. Disables the usage of `setindex!`.

See also `thaw`.
"""
function freeze end

"""
    thaw(container)

Thaw an indexable container to make it mutable. Enables the usage of `setindex!`.

See also `freeze`.
"""
function thaw end

include("array.jl")

end # module

ext: jl


text:
@testset "Emulation Model Build" begin
    template = get_thermal_dispatch_template_network()
    c_sys5 = PSB.build_system(
        PSITestSystems,
        "c_sys5_uc";
        add_single_time_series=true,
        force_build=true,
    )

    model = EmulationModel(template, c_sys5; optimizer=GLPK_optimizer)
    @test build!(model; executions=10, output_dir=mktempdir(cleanup=true)) ==
          BuildStatus.BUILT
    @test run!(model) == RunStatus.SUCCESSFUL

    template = get_thermal_standard_uc_template()
    c_sys5_uc_re = PSB.build_system(
        PSITestSystems,
        "c_sys5_uc_re";
        add_single_time_series=true,
        force_build=true,
    )
    set_device_model!(template, RenewableDispatch, RenewableFullDispatch)
    model = EmulationModel(template, c_sys5_uc_re; optimizer=GLPK_optimizer)

    @test build!(model; executions=10, output_dir=mktempdir(cleanup=true)) ==
          BuildStatus.BUILT
    @test run!(model) == RunStatus.SUCCESSFUL
    @test !isempty(collect(readdir(PSI.get_recorder_dir(model))))
end

@testset "Emulation Model initial_conditions test for ThermalGen" begin
    ######## Test with ThermalStandardUnitCommitment ########
    template = get_thermal_standard_uc_template()
    c_sys5_uc_re = PSB.build_system(
        PSITestSystems,
        "c_sys5_uc_re";
        add_single_time_series=true,
        force_build=true,
    )
    set_device_model!(template, RenewableDispatch, RenewableFullDispatch)
    model = EmulationModel(template, c_sys5_uc_re; optimizer=HiGHS_optimizer)
    @test build!(model; executions=10, output_dir=mktempdir(cleanup=true)) ==
          BuildStatus.BUILT
    check_duration_on_initial_conditions_values(model, ThermalStandard)
    check_duration_off_initial_conditions_values(model, ThermalStandard)
    @test run!(model) == RunStatus.SUCCESSFUL

    ######## Test with ThermalMultiStartUnitCommitment ########
    template = get_thermal_standard_uc_template()
    c_sys5_uc = PSB.build_system(
        PSITestSystems,
        "c_sys5_pglib";
        add_single_time_series=true,
        force_build=true,
    )
    set_device_model!(template, ThermalMultiStart, ThermalMultiStartUnitCommitment)
    model = EmulationModel(template, c_sys5_uc; optimizer=HiGHS_optimizer)
    @test build!(model; executions=1, output_dir=mktempdir(cleanup=true)) ==
          BuildStatus.BUILT

    check_duration_on_initial_conditions_values(model, ThermalStandard)
    check_duration_off_initial_conditions_values(model, ThermalStandard)
    check_duration_on_initial_conditions_values(model, ThermalMultiStart)
    check_duration_off_initial_conditions_values(model, ThermalMultiStart)
    @test run!(model) == RunStatus.SUCCESSFUL

    ######## Test with ThermalCompactUnitCommitment ########
    template = get_thermal_standard_uc_template()
    c_sys5_uc = PSB.build_system(
        PSITestSystems,
        "c_sys5_pglib";
        add_single_time_series=true,
        force_build=true,
    )
    set_device_model!(template, ThermalMultiStart, ThermalCompactUnitCommitment)
    model = EmulationModel(template, c_sys5_uc; optimizer=HiGHS_optimizer)
    @test build!(model; executions=1, output_dir=mktempdir(cleanup=true)) ==
          BuildStatus.BUILT
    check_duration_on_initial_conditions_values(model, ThermalStandard)
    check_duration_off_initial_conditions_values(model, ThermalStandard)
    check_duration_on_initial_conditions_values(model, ThermalMultiStart)
    check_duration_off_initial_conditions_values(model, ThermalMultiStart)
    @test run!(model) == RunStatus.SUCCESSFUL

    ######## Test with ThermalCompactDispatch ########
    template = get_thermal_standard_uc_template()
    c_sys5_uc = PSB.build_system(
        PSITestSystems,
        "c_sys5_pglib";
        add_single_time_series=true,
        force_build=true,
    )
    device_model = DeviceModel(PSY.ThermalStandard, PSI.ThermalCompactDispatch)
    set_device_model!(template, device_model)
    model = EmulationModel(template, c_sys5_uc; optimizer=HiGHS_optimizer)
    @test build!(model; executions=10, output_dir=mktempdir(cleanup=true)) ==
          BuildStatus.BUILT
end

@testset "Emulation Model initial_conditions test for Storage" begin
    ######## Test with BookKeeping ########
    template = get_thermal_dispatch_template_network()
    c_sys5_bat = PSB.build_system(
        PSITestSystems,
        "c_sys5_bat";
        add_single_time_series=true,
        force_build=true,
    )
    set_device_model!(template, GenericBattery, BookKeeping)
    model = EmulationModel(template, c_sys5_bat; optimizer=HiGHS_optimizer)
    @test build!(model; executions=10, output_dir=mktempdir(cleanup=true)) ==
          BuildStatus.BUILT
    ic_data = PSI.get_initial_condition(
        PSI.get_optimization_container(model),
        InitialEnergyLevel(),
        GenericBattery,
    )
    for ic in ic_data
        name = PSY.get_name(ic.component)
        e_var = JuMP.value(PSI.get_value(ic))
        @test PSY.get_initial_energy(ic.component) == e_var
    end
    @test run!(model) == RunStatus.SUCCESSFUL

    ######## Test with BatteryAncillaryServices ########
    template = get_thermal_dispatch_template_network()
    c_sys5_bat = PSB.build_system(
        PSITestSystems,
        "c_sys5_bat";
        add_single_time_series=true,
        force_build=true,
    )
    set_device_model!(template, GenericBattery, BatteryAncillaryServices)
    model = EmulationModel(template, c_sys5_bat; optimizer=HiGHS_optimizer)
    @test build!(model; executions=10, output_dir=mktempdir(cleanup=true)) ==
          BuildStatus.BUILT
    ic_data = PSI.get_initial_condition(
        PSI.get_optimization_container(model),
        InitialEnergyLevel(),
        GenericBattery,
    )
    for ic in ic_data
        name = PSY.get_name(ic.component)
        e_var = JuMP.value(PSI.get_value(ic))
        @test PSY.get_initial_energy(ic.component) == e_var
    end
    @test run!(model) == RunStatus.SUCCESSFUL

    ######## Test with EnergyTarget ########
    template = get_thermal_dispatch_template_network()
    c_sys5_bat = PSB.build_system(
        PSITestSystems,
        "c_sys5_bat_ems";
        add_single_time_series=true,
        force_build=true,
    )
    set_device_model!(template, BatteryEMS, EnergyTarget)
    model = EmulationModel(template, c_sys5_bat; optimizer=HiGHS_optimizer)
    @test build!(model; executions=10, output_dir=mktempdir(cleanup=true)) ==
          BuildStatus.BUILT
    ic_data = PSI.get_initial_condition(
        PSI.get_optimization_container(model),
        InitialEnergyLevel(),
        BatteryEMS,
    )
    for ic in ic_data
        name = PSY.get_name(ic.component)
        e_var = JuMP.value(PSI.get_value(ic))
        @test PSY.get_initial_energy(ic.component) == e_var
    end
    @test run!(model) == RunStatus.SUCCESSFUL
end

@testset "Emulation Model initial_conditions test for Hydro" begin
    ######## Test with HydroDispatchRunOfRiver ########
    template = get_thermal_dispatch_template_network()
    c_sys5_hyd = PSB.build_system(
        PSITestSystems,
        "c_sys5_hyd";
        add_single_time_series=true,
        force_build=true,
    )
    set_device_model!(template, HydroDispatch, HydroDispatchRunOfRiver)
    set_device_model!(template, HydroEnergyReservoir, HydroDispatchRunOfRiver)
    model = EmulationModel(template, c_sys5_hyd; optimizer=HiGHS_optimizer)
    @test build!(model; executions=10, output_dir=mktempdir(cleanup=true)) ==
          BuildStatus.BUILT
    initial_conditions_data =
        PSI.get_initial_conditions_data(PSI.get_optimization_container(model))
    @test !PSI.has_initial_condition_value(
        initial_conditions_data,
        ActivePowerVariable(),
        HydroEnergyReservoir,
    )
    @test run!(model) == RunStatus.SUCCESSFUL

    ######## Test with HydroCommitmentRunOfRiver ########
    template = get_thermal_dispatch_template_network()
    c_sys5_hyd = PSB.build_system(
        PSITestSystems,
        "c_sys5_hyd";
        add_single_time_series=true,
        force_build=true,
    )
    set_device_model!(template, HydroDispatch, HydroCommitmentRunOfRiver)
    set_device_model!(template, HydroEnergyReservoir, HydroCommitmentRunOfRiver)
    model = EmulationModel(template, c_sys5_hyd; optimizer=HiGHS_optimizer)

    @test build!(model; executions=10, output_dir=mktempdir(cleanup=true)) ==
          BuildStatus.BUILT
    initial_conditions_data =
        PSI.get_initial_conditions_data(PSI.get_optimization_container(model))
    @test PSI.has_initial_condition_value(
        initial_conditions_data,
        OnVariable(),
        HydroEnergyReservoir,
    )
    @test run!(model) == RunStatus.SUCCESSFUL

    ######## Test with HydroCommitmentRunOfRiver ########
    template = get_thermal_dispatch_template_network()
    c_sys5_hyd = PSB.build_system(
        PSITestSystems,
        "c_sys5_hyd";
        add_single_time_series=true,
        force_build=true,
    )
    set_device_model!(template, HydroEnergyReservoir, HydroDispatchReservoirBudget)
    model = EmulationModel(template, c_sys5_hyd; optimizer=HiGHS_optimizer)
    @test build!(model; executions=10, output_dir=mktempdir(cleanup=true)) ==
          BuildStatus.BUILT
    initial_conditions_data =
        PSI.get_initial_conditions_data(PSI.get_optimization_container(model))
    @test !PSI.has_initial_condition_value(
        initial_conditions_data,
        ActivePowerVariable(),
        HydroEnergyReservoir,
    )
    @test run!(model) == RunStatus.SUCCESSFUL

    ######## Test with HydroCommitmentRunOfRiver ########
    template = get_thermal_dispatch_template_network()
    c_sys5_hyd = PSB.build_system(
        PSITestSystems,
        "c_sys5_hyd";
        add_single_time_series=true,
        force_build=true,
    )
    set_device_model!(template, HydroEnergyReservoir, HydroCommitmentReservoirBudget)
    model = EmulationModel(template, c_sys5_hyd; optimizer=HiGHS_optimizer)
    @test build!(model; executions=10, output_dir=mktempdir(cleanup=true)) ==
          BuildStatus.BUILT
    initial_conditions_data =
        PSI.get_initial_conditions_data(PSI.get_optimization_container(model))
    @test PSI.has_initial_condition_value(
        initial_conditions_data,
        OnVariable(),
        HydroEnergyReservoir,
    )
    @test run!(model) == RunStatus.SUCCESSFUL

    ######## Test with HydroCommitmentRunOfRiver ########
    template = get_thermal_dispatch_template_network()
    c_sys5_hyd = PSB.build_system(
        PSITestSystems,
        "c_sys5_hyd_ems";
        add_single_time_series=true,
        force_build=true,
    )
    set_device_model!(template, HydroEnergyReservoir, HydroDispatchReservoirStorage)
    model = EmulationModel(template, c_sys5_hyd; optimizer=HiGHS_optimizer)
    @test build!(model; executions=10, output_dir=mktempdir(cleanup=true)) ==
          BuildStatus.BUILT
    initial_conditions_data =
        PSI.get_initial_conditions_data(PSI.get_optimization_container(model))
    @test !PSI.has_initial_condition_value(
        initial_conditions_data,
        ActivePowerVariable(),
        HydroEnergyReservoir,
    )
    ic_data = PSI.get_initial_condition(
        PSI.get_optimization_container(model),
        InitialEnergyLevel(),
        HydroEnergyReservoir,
    )
    for ic in ic_data
        name = PSY.get_name(ic.component)
        e_var = JuMP.value(PSI.get_value(ic))
        @test PSY.get_initial_storage(ic.component) == e_var
    end
    @test run!(model) == RunStatus.SUCCESSFUL

    ######## Test with HydroCommitmentRunOfRiver ########
    template = get_thermal_dispatch_template_network()
    c_sys5_hyd = PSB.build_system(
        PSITestSystems,
        "c_sys5_hyd_ems";
        add_single_time_series=true,
        force_build=true,
    )
    set_device_model!(template, HydroEnergyReservoir, HydroCommitmentReservoirStorage)
    model = EmulationModel(template, c_sys5_hyd; optimizer=HiGHS_optimizer)
    @test build!(model; executions=10, output_dir=mktempdir(cleanup=true)) ==
          BuildStatus.BUILT
    initial_conditions_data =
        PSI.get_initial_conditions_data(PSI.get_optimization_container(model))
    @test PSI.has_initial_condition_value(
        initial_conditions_data,
        OnVariable(),
        HydroEnergyReservoir,
    )
    ic_data = PSI.get_initial_condition(
        PSI.get_optimization_container(model),
        InitialEnergyLevel(),
        HydroEnergyReservoir,
    )
    for ic in ic_data
        name = PSY.get_name(ic.component)
        e_var = JuMP.value(PSI.get_value(ic))
        @test PSY.get_initial_storage(ic.component) == e_var
    end
    @test run!(model) == RunStatus.SUCCESSFUL
end

@testset "Emulation Model Results" begin
    template = get_thermal_dispatch_template_network()
    c_sys5 = PSB.build_system(
        PSITestSystems,
        "c_sys5_uc";
        add_single_time_series=true,
        force_build=true,
    )

    model = EmulationModel(template, c_sys5; optimizer=HiGHS_optimizer)
    executions = 10
    @test build!(model; executions=executions, output_dir=mktempdir(cleanup=true)) ==
          BuildStatus.BUILT
    @test run!(model) == RunStatus.SUCCESSFUL
    results = ProblemResults(model)
    @test list_aux_variable_names(results) == []
    @test list_aux_variable_keys(results) == []
    @test list_variable_names(results) == ["ActivePowerVariable__ThermalStandard"]
    @test list_variable_keys(results) ==
          [PSI.VariableKey(ActivePowerVariable, ThermalStandard)]
    @test list_dual_names(results) == []
    @test list_dual_keys(results) == []
    @test list_parameter_names(results) == ["ActivePowerTimeSeriesParameter__PowerLoad"]
    @test list_parameter_keys(results) ==
          [PSI.ParameterKey(ActivePowerTimeSeriesParameter, PowerLoad)]

    @test read_variable(results, "ActivePowerVariable__ThermalStandard") isa DataFrame
    @test read_variable(results, ActivePowerVariable, ThermalStandard) isa DataFrame
    @test read_variable(results, PSI.VariableKey(ActivePowerVariable, ThermalStandard)) isa
          DataFrame

    @test read_parameter(results, "ActivePowerTimeSeriesParameter__PowerLoad") isa DataFrame
    @test read_parameter(results, ActivePowerTimeSeriesParameter, PowerLoad) isa DataFrame
    @test read_parameter(
        results,
        PSI.ParameterKey(ActivePowerTimeSeriesParameter, PowerLoad),
    ) isa DataFrame

    @test read_optimizer_stats(model) isa DataFrame
    for n in names(read_optimizer_stats(model))
        stats_values = read_optimizer_stats(model)[!, n]
        if any(ismissing.(stats_values))
            @test ismissing.(stats_values) ==
                  ismissing.(read_optimizer_stats(results)[!, n])
        elseif any(isnan.(stats_values))
            @test isnan.(stats_values) == isnan.(read_optimizer_stats(results)[!, n])
        else
            @test stats_values == read_optimizer_stats(results)[!, n]
        end
    end

    for i in 1:executions
        @test get_objective_value(results, i) isa Float64
    end
end

@testset "Run EmulationModel with auto-build" begin
    for serialize in (true, false)
        template = get_thermal_dispatch_template_network()
        c_sys5 = PSB.build_system(
            PSITestSystems,
            "c_sys5_uc";
            add_single_time_series=true,
            force_build=true,
        )

        model = EmulationModel(template, c_sys5; optimizer=HiGHS_optimizer)
        @test_throws ErrorException run!(model, executions=10)
        @test run!(
            model,
            executions=10,
            output_dir=mktempdir(cleanup=true),
            serialize=serialize,
        ) == RunStatus.SUCCESSFUL
    end
end

@testset "Test serialization/deserialization of EmulationModel results" begin
    path = mktempdir(cleanup=true)
    template = get_thermal_dispatch_template_network()
    c_sys5 = PSB.build_system(
        PSITestSystems,
        "c_sys5_uc";
        add_single_time_series=true,
        force_build=true,
    )

    model = EmulationModel(template, c_sys5; optimizer=HiGHS_optimizer)
    executions = 10
    @test build!(model; executions=executions, output_dir=path) == BuildStatus.BUILT
    @test run!(model, export_problem_results=true) == RunStatus.SUCCESSFUL
    results1 = ProblemResults(model)
    var1_a = read_variable(results1, ActivePowerVariable, ThermalStandard)
    # Ensure that we can deserialize strings into keys.
    var1_b = read_variable(results1, "ActivePowerVariable__ThermalStandard")
    @test var1_a == var1_b

    # Results were automatically serialized here.
    results2 = ProblemResults(joinpath(PSI.get_output_dir(model)))
    var2 = read_variable(results2, ActivePowerVariable, ThermalStandard)
    @test var1_a == var2
    @test get_system(results2) !== nothing

    # Serialize to a new directory with the exported function.
    results_path = joinpath(path, "results")
    serialize_results(results1, results_path)
    @test isfile(joinpath(results_path, PSI._PROBLEM_RESULTS_FILENAME))
    results3 = ProblemResults(results_path)
    var3 = read_variable(results3, ActivePowerVariable, ThermalStandard)
    @test var1_a == var3
    @test get_system(results3) === nothing
    set_system!(results3, get_system(results1))
    @test get_system(results3) !== nothing

    exp_file =
        joinpath(path, "results", "variables", "ActivePowerVariable__ThermalStandard.csv")
    var4 = PSI.read_dataframe(exp_file)
    # Manually Multiply by the base power var1_a has natural units and export writes directly from the solver
    @test var1_a[:, propertynames(var1_a) .!= :DateTime] == var4 .* 100.0
end

@testset "Test deserialization and re-run of EmulationModel" begin
    path = mktempdir(cleanup=true)
    template = get_thermal_dispatch_template_network()
    c_sys5 = PSB.build_system(
        PSITestSystems,
        "c_sys5_uc";
        add_single_time_series=true,
        force_build=true,
    )

    model = EmulationModel(template, c_sys5; optimizer=HiGHS_optimizer)
    executions = 10
    @test build!(model; executions=executions, output_dir=path) == BuildStatus.BUILT
    @test run!(model) == RunStatus.SUCCESSFUL
    results = ProblemResults(model)
    var1 = read_variable(results, ActivePowerVariable, ThermalStandard)

    file_list = sort!(collect(readdir(path)))
    @test PSI._JUMP_MODEL_FILENAME in file_list
    @test PSI._SERIALIZED_MODEL_FILENAME in file_list
    path2 = joinpath(path, "tmp")
    model2 = EmulationModel(path, HiGHS_optimizer)
    build!(model2, output_dir=path2)
    @test run!(model2) == RunStatus.SUCCESSFUL
    results2 = ProblemResults(model2)
    var2 = read_variable(results, ActivePowerVariable, ThermalStandard)

    @test var1 == var2

    # Deserialize with different optimizer attributes.
    optimizer = JuMP.optimizer_with_attributes(HiGHS.Optimizer, "time_limit" => 110.0)
    @test_logs (:warn, r"Different optimizer attributes are set") match_mode = :any EmulationModel(
        path,
        optimizer,
    )

    # Deserialize with a different optimizer.
    @test_logs (:warn, r"Original solver was .* new solver is") match_mode = :any EmulationModel(
        path,
        GLPK_optimizer,
    )
end

@testset "Test serialization of InitialConditionsData" begin
    template = get_thermal_standard_uc_template()
    sys = PSB.build_system(
        PSITestSystems,
        "c_sys5_pglib";
        add_single_time_series=true,
        force_build=true,
    )
    optimizer = HiGHS_optimizer
    set_device_model!(template, ThermalMultiStart, ThermalMultiStartUnitCommitment)
    model = EmulationModel(template, sys; optimizer=HiGHS_optimizer)
    output_dir = mktempdir(cleanup=true)

    @test build!(model; executions=1, output_dir=output_dir) == BuildStatus.BUILT
    ic_file = PSI.get_initial_conditions_file(model)
    test_ic_serialization_outputs(model, ic_file_exists=true, message="make")
    @test run!(model) == RunStatus.SUCCESSFUL

    # Build again, use existing initial conditions.
    PSI.reset!(model)
    @test build!(model; executions=1, output_dir=output_dir) == PSI.BuildStatus.BUILT
    test_ic_serialization_outputs(model, ic_file_exists=true, message="make")
    @test run!(model) == RunStatus.SUCCESSFUL

    # Build again, use existing initial conditions.
    model = EmulationModel(
        template,
        sys;
        optimizer=optimizer,
        deserialize_initial_conditions=true,
    )
    @test build!(model; executions=1, output_dir=output_dir) == PSI.BuildStatus.BUILT
    test_ic_serialization_outputs(model, ic_file_exists=true, message="deserialize")
    @test run!(model) == RunStatus.SUCCESSFUL

    # Construct and build again with custom initial conditions file.
    initialization_file = joinpath(output_dir, ic_file * ".old")
    mv(ic_file, initialization_file)
    touch(ic_file)
    model = EmulationModel(
        template,
        sys;
        optimizer=optimizer,
        initialization_file=initialization_file,
        deserialize_initial_conditions=true,
    )
    @test build!(model; executions=1, output_dir=output_dir) == PSI.BuildStatus.BUILT
    test_ic_serialization_outputs(model, ic_file_exists=true, message="deserialize")

    # Construct and build again while skipping build of initial conditions.
    model = EmulationModel(template, sys; optimizer=optimizer, initialize_model=false)
    rm(ic_file)
    @test build!(model; executions=1, output_dir=output_dir) == PSI.BuildStatus.BUILT
    test_ic_serialization_outputs(model, ic_file_exists=false, message="skip")
    @test run!(model) == RunStatus.SUCCESSFUL
end


text:
export diric

"""
    kernel = diric(Ω::Real, n::Integer)

Dirichlet kernel, also known as periodic sinc function,
where `n` should be a positive integer.
Returns `sin(Ω * n/2) / (n * sin(Ω / 2))` typically,
but `±1` when `Ω` is a multiple of 2π.

In the usual case where 'n' is odd, the output is equivalent to
``1/n \\sum_{k=-(n-1)/2}^{(n-1)/2} e^{i k Ω}``,
which is the discrete-time Fourier transform (DTFT)
of a `n`-point moving average filter.

When `n` is odd or even, the function is 2π or 4π periodic, respectively.
The formula for general `n` is
`diric(Ω,n) = ``e^{-i (n-1) Ω/2}/n \\sum_{k=0}^{n-1} e^{i k Ω}``,
which is a real and symmetric function of `Ω`.

As of 2021-03-19, the Wikipedia definition has different factors.
The definition here is consistent with scipy and other software frameworks.

# Examples

```jldoctest
julia> round.(diric.((-2:0.5:2)*π, 5), digits=9)'
1×9 adjoint(::Vector{Float64}) with eltype Float64:
 1.0  -0.2  0.2  -0.2  1.0  -0.2  0.2  -0.2  1.0

julia> diric(0, 4)
1.0
```
"""
function diric(Ω::T, n::Integer) where T <: AbstractFloat
    n > 0 || throw(ArgumentError("n=$n not positive"))
    sign = one(T)
    if isodd(n)
        Ω = rem2pi(Ω, RoundNearest) # [-π,π)
    else # wrap to interval [-π,π) to improve precision near ±2π
        Ω = 2 * rem2pi(Ω/2, RoundNearest) # [-2π,2π)
        if Ω > π # [π,2π)
            sign = -one(T)
            Ω -= T(2π) # [-π,0)
        elseif Ω < -π # [-2π,-π)
            sign = -one(T)
            Ω += T(2π) # [0,π)
        end
    end

    denom = sin(Ω / 2)
    atol = eps(T)
    if abs(denom) ≤ atol # denom ≈ 0 ?
        return sign
    end

    return sign * sin(Ω * n/2) / (n * denom) # typical case
end

# handle non AbstractFloat types (e.g., Int, Rational, π)
diric(Ω::Real, n::Integer) = diric(float(Ω), n)


text:
# This function calculates the largest eigenvalue
# (absolute value wise) by power iteration.
function maxeig!(integrator, cache::StochasticDiffEqConstantCache)
  isfirst = integrator.iter == 1 || integrator.u_modified
  @unpack t, dt, uprev, u, p = integrator
  maxiter =  50
  safe =  1.2
  fsalfirst = integrator.f(uprev, p, t)
  # Initial guess for eigenvector `z`
  if isfirst
    fz = fsalfirst
    z = integrator.f(fz, p, t)
  else
    z = cache.zprev
  end
  # Perturbation
  u_norm = integrator.opts.internalnorm(uprev,t)
  z_norm = integrator.opts.internalnorm(z,t)
  pert   = eps(u_norm)
  sqrt_pert = sqrt(pert)
  is_u_zero = u_norm == zero(u_norm)
  is_z_zero = z_norm == zero(z_norm)
  # Normalize `z` such that z-u lie in a circle
  if ( !is_u_zero && !is_z_zero )
    dz_u = u_norm * sqrt_pert
    quot = dz_u/z_norm
    z = uprev + quot*z
  elseif !is_u_zero
    dz_u = u_norm * sqrt_pert
    z = uprev + uprev*dz_u
  elseif !is_z_zero
    dz_u = pert
    quot = dz_u/z_norm
    z *= quot
  else
    dz_u = pert
    z = dz_u .* (false .* z .+ one(eltype(z)))
  end # endif
  # Start power iteration
  integrator.eigen_est = 0
  for iter in 1:maxiter
    fz = integrator.f(z, p, t)
    fz = fz - fsalfirst

    Δ  = integrator.opts.internalnorm(fz,t)
    eig_prev = integrator.eigen_est
    integrator.eigen_est = Δ/dz_u * safe
    # Convergence
    if iter >= 2 && abs(eig_prev - integrator.eigen_est) < integrator.eigen_est*0.05
      # Store the eigenvector
      cache.zprev = z - uprev
      return true
    end

    # Next `z`
    if Δ != zero(Δ)
      quot = dz_u/Δ
      z = uprev + quot*fz
    else
        # An arbitrary change on `z`
        nind = length(z)
        if (nind != 1)
          ind = 1 + iter % nind
          # val = (uprev[ind] - (z[ind] - uprev[ind]))*one(eltype(z))*2
          vec(z) = vec(z) .* (1 .- 2 .*((1:length(z)) .== ind))
        else
          z = -z
        end
    end
  end
  return false
end

function maxeig!(integrator, cache::StochasticDiffEqMutableCache)
  isfirst = integrator.iter == 1 || integrator.u_modified
  @unpack t, dt, uprev, u, p = integrator
  fz, z, fsalfirst = cache.atmp, cache.tmp, cache.fsalfirst
  integrator.f(fsalfirst, uprev, p, t)
  ccache = cache.constantcache
  maxiter =  50
  safe = 1.2

  # Initial guess for eigenvector `z`
  if isfirst
    @.. fz = fsalfirst
    integrator.f(z, fz, p, t)
    # integrator.destats.nf += 1
  else
    @.. z = ccache.zprev
  end
  # Perturbation
  u_norm = integrator.opts.internalnorm(uprev,t)
  z_norm = integrator.opts.internalnorm(z,t)
  pert   = eps(u_norm)
  sqrt_pert = sqrt(pert)
  is_u_zero = u_norm == zero(u_norm)
  is_z_zero = z_norm == zero(z_norm)
  # Normalize `z` such that z-u lie in a circle
  if ( !is_u_zero && !is_z_zero )
    dz_u = u_norm * sqrt_pert
    quot = dz_u/z_norm
    @.. z = uprev + quot*z
  elseif !is_u_zero
    dz_u = u_norm * sqrt_pert
    @.. z = uprev + uprev*dz_u
  elseif !is_z_zero
    dz_u = pert
    quot = dz_u/z_norm
    @.. z *= quot
  else
    dz_u = pert
    @.. z = dz_u*(false*z + one(eltype(z)))
  end # endif
  # Start power iteration
  integrator.eigen_est = 0
  for iter in 1:maxiter
    integrator.f(fz, z, p, t)
    # integrator.destats.nf += 1
    @.. fz = fz - fsalfirst

    Δ  = integrator.opts.internalnorm(fz,t)
    eig_prev = integrator.eigen_est
    integrator.eigen_est = Δ/dz_u * safe
    # Convergence
    if iter >= 2 && abs(eig_prev - integrator.eigen_est) < integrator.eigen_est*0.05
      # Store the eigenvector
      @.. ccache.zprev = z - uprev
      return true
    end
    # Next `z`
    if Δ != zero(Δ)
      quot = dz_u/Δ
      @.. z = uprev + quot*fz
    else
      # An arbitrary change on `z`
       nind = length(uprev)
       if (nind != 1)
         ind = 1 + iter % nind
         # val = (uprev[ind] - (z[ind] - uprev[ind]))*one(eltype(z))
         vec(z) .= vec(z) .* (1 .- 2 .*((1:length(z)) .== ind))
       else
         z = -z
       end
    end
  end
  return false
end

function choose_deg!(integrator,cache::T) where T
  isconst = T <: StochasticDiffEqConstantCache
  isconst || ( cache = cache.constantcache )

  if integrator.alg isa SROCK1
    # binary search as stability domain is monotonically incrasing with number of stages
    mn_st, mx_st, mid_st = 3, 200, 3
    while (mx_st - mn_st > 1)
      mid_st = Int(floor((mn_st+mx_st)*0.5))
      if (mid_st^2*(131.97/197 - 0.45*mid_st/197) > cache.mdeg)
        mx_st = mid_st
      else
        mn_st = mid_st
      end
    end
    cache.mdeg = (cache.mdeg > mn_st^2*(131.97/197 - 0.45*mn_st/197)) ? mx_st : mn_st

    @inbounds for i in 1:size(cache.ms,1)
      if cache.ms[i] <= cache.mdeg
        cache.optimal_η = cache.mη[i]
      else
        break
      end
    end
  end

  if integrator.alg isa SROCK2
    start = 1
    @inbounds for i in 1:size(cache.ms,1)
      if cache.ms[i] >= cache.mdeg
        cache.deg_index = i;
        cache.mdeg = cache.ms[i]
        cache.start = start
        break
      else
        start += cache.ms[i]*2 - 1
      end
    end
  end

  if integrator.alg isa TangXiaoSROCK2
    start = 1
    @inbounds for i in 1:size(cache.ms,1)
      if cache.ms[i] >= cache.mdeg
        cache.deg_index = i;
        cache.mdeg = cache.ms[i]
        cache.start = start
        break
      else
        start += cache.ms[i]*2 - 1
      end
    end

    start = 1
    @inbounds for i in 1:5
      if integrator.alg.version_num == i
        cache.start_mcs = start
        break
      else
        start += cache.mn̂[i]
      end
    end
  end

  if integrator.alg isa SROCKEM
    @inbounds for i in 1:size(cache.ms,1)
      if cache.ms[i] <= cache.mdeg
        cache.optimal_η = cache.mη[i]
      else
        break
      end
    end
  end

  if integrator.alg isa KomBurSROCK2
    start = 1
    @inbounds for i in 1:size(cache.ms,1)
      if cache.ms[i] >= cache.mdeg
        cache.deg_index = i;
        cache.mdeg = cache.ms[i]
        cache.start = start
        break
      else
        start += cache.ms[i]*2 - 1
      end
    end
  end

  if integrator.alg isa SROCKC2
    start = 1
    @inbounds for i in 1:size(cache.ms,1)
      if cache.ms[i] >= cache.mdeg
        cache.deg_index = i;
        cache.mdeg = cache.ms[i]
        cache.start = start
        break
      else
        start += cache.ms[i]*2 - 1
      end
    end
  end

  return nothing
end


text:
@testset "geometry.embed2d" begin

@testset "embed2d" begin
    return
    mol = smilestomol("CCC1CCCC1C")
    zmatrix = coords2d(mol)
    display(zmatrix)
    println("")
    coords3d = cartesian(zmatrix)
    display(coords3d)
    println("")
    for (i, atom) in enumerate(atomvector(mol))
        atom.coords = tuple(coords3d[i, 1:2]...)
    end
    ASSETS_DIR = joinpath(dirname(@__FILE__), "..", "..", "assets")
    dest = open(joinpath(ASSETS_DIR, "image", "test.svg"), "w")
    write(dest, drawsvg(mol, 200, 200))
end

end # geometry.embed2d


text:
#=
   Ruin problem.

   Cf ruin_problem.jl

   Compared to ruin_problem.jl this version uses a recursive approach
   with concatenated arrays instead.

   The arrays always ends with a 0 (ruin).
   Note: We have to have an maxLen since otherwise the sequence goes to
   infinity...

   Summary Statistics
   parameters      mean       std   naive_se      mcse          ess      rhat   ess_per_sec
       Symbol   Float64   Float64    Float64   Float64      Float64   Float64       Float64

          len   11.0188   15.7735     0.1577    0.1505    9765.3342    1.0001    30233.2330
         ruin    0.8925    0.3098     0.0031    0.0028   10015.3143    0.9999    31007.1649

   Distributions of variable len (num:0)
   2.00000 =>    4999  (0.499900)
   4.00000 =>    1320  (0.132000)
   50.00000 =>    1097  (0.109700)
   6.00000 =>     603  (0.060300)
   8.00000 =>     393  (0.039300)
   10.00000 =>     265  (0.026500)
   12.00000 =>     204  (0.020400)
   14.00000 =>     148  (0.014800)
   16.00000 =>     126  (0.012600)
   18.00000 =>     126  (0.012600)
   20.00000 =>      99  (0.009900)
   22.00000 =>      87  (0.008700)
   24.00000 =>      65  (0.006500)
   30.00000 =>      58  (0.005800)
   26.00000 =>      56  (0.005600)
   28.00000 =>      51  (0.005100)
   36.00000 =>      48  (0.004800)
   32.00000 =>      42  (0.004200)
   38.00000 =>      41  (0.004100)
   34.00000 =>      39  (0.003900)
   40.00000 =>      30  (0.003000)
   46.00000 =>      29  (0.002900)
   42.00000 =>      28  (0.002800)
   44.00000 =>      27  (0.002700)
   48.00000 =>      19  (0.001900)

   Distributions of variable ruin (num:0)
   1.00000 =>    8925  (0.892500)
   0.00000 =>    1075  (0.107500)

.  Cf ~/webppl/ruin_problem2.wppl

=#
using Turing, StatsPlots
include("jl_utils.jl")

@model function ruin_problem2(maxLen=50,start=1,win=1,loose=1)
    head = 1
    tail = 2
    coins = [head, tail]

    function draw(arr)
        if length(arr) == 0
            return draw([start])
        elseif arr[end] == 0 || length(arr) >= maxLen
            return arr
        else
            lastVal = arr[end]
            c = rand(coins)
            if c == head
                return draw(push!(arr,lastVal+win))
            else
                if lastVal - loose <= 0
                    return push!(arr,0)
                else
                    return draw(push!(arr,lastVal-loose))
                end
            end
        end
    end

    a = draw([])
    len ~ Dirac(length(a))
    ruin ~ Dirac(a[end] == 0)

    # Use PG(5) when using these. Otherwise Prior() works and is faster.
    # true ~ Dirac(len == 6)
    # true ~ Dirac(a[end] ==8)

end

maxLen=50
start=1
win=1
loose=1
model = ruin_problem2(maxLen,start,win,loose)

chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns, :len)
show_var_dist_pct(chns, :ruin)

# chains_params = Turing.MCMCChains.get_sections(chns, :parameters)
# genq = generated_quantities(model, chains_params)
# show_var_dist_pct(genq)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
export TimeSolver

abstract type TimeSolver end

show(io::IO, s::TimeSolver) =
    try print(io,s.info)
    catch
        print(io,"Time solver: $(s.label)")
    end

Base.:(==)(a::TimeSolver, b::TimeSolver) = begin

    init = Random((N=10, L=10))

    Na = try size(a.U1,1)
    catch
        2
    end
    ma = Airy((μ=1, N=Na, L=10))
    Ua = ma.mapto(init)

    Nb = try size(b.U1,1)
    catch
        Complex.(rand(2,2))
    end
    mb = Airy((μ=1, N=Nb, L=10))
    Ub = mb.mapto(init)

    step!(a, ma, Ua, 1)
    step!(b, mb, Ub, 1)

    a.label == b.label && Ua == Ub

end


text:
WIDTH = 800
HEIGHT = 600
RATIO = WIDTH/(HEIGHT*1f0)
SIZE = WIDTH * HEIGHT
FOV = 60.0f0
CLIP_NEAR = 0.001f0
CLIP_FAR = 10000.0f0

"""
sets glfw window size + viewport
"""
function rezizeWindow(width,height)
  global WIDTH, HEIGHT, RATIO, SIZE
  WIDTH = width
  HEIGHT = height
  RATIO = WIDTH/(HEIGHT*1f0)
  SIZE = WIDTH * HEIGHT
  GLFW.SetWindowSize(window, WIDTH, HEIGHT)
  glViewport(0, 0, WIDTH, HEIGHT)
end

lastCursorPos = [0,0]
cursorPos_old = [0,0]
cursorPos = [0,0]
mouseMove = [0,0]

#MVP = eye(Mat4x4f)
#MVP *= transform([0.0f0,0.0f0,0.0f0],[0.0f0,0.0f0,0.0f0],Float32[])

keyFB = 0
keyLR = 0
keyUD = 0
speed = false

mouseKeyPressed = false
keyPressed = false
keyPressing = false
keyValue = 0

"""
camera object with holds position, rotation, scaling and various matrices like MVP
"""
type Camera
  moved::Bool

  position::AbstractArray
  rotation::AbstractArray
  scale::AbstractArray

  viewMat::AbstractArray
  projectionMat::AbstractArray
  modelMat::AbstractArray

  translateMat::AbstractArray
  rotationMat::AbstractArray
  scalingMat::AbstractArray

  MVP::AbstractArray

  Camera() = new(true, zerosVector3f,zerosVector3f,onesVector3f,
  eyeMat4x4f,eyeMat4x4f,eyeMat4x4f,eyeMat4x4f,eyeMat4x4f,eyeMat4x4f,eyeMat4x4f)
end

CAMERA = Camera()

"""
gets forward vector of camera direction
"""
forward(camera::Camera) = forward(camera.rotationMat)

"""
gets right vector of camera direction
"""
right(camera::Camera) = right(camera.rotationMat)

"""
gets up vector of camera direction
"""
up(camera::Camera) = up(camera.rotationMat)

"""
sets projection matrix
"""
setProjection(camera::Camera, m::AbstractArray) = (camera.projectionMat = m)

"""
sets view matrix
"""
setView(camera::Camera, m::AbstractArray) = (camera.viewMat = m)

VIEW_KEYS=false

"""
event which catches keyboard inputs.
here keys for wireframe, fullscreen and camera movement are defined
"""
function OnKey(window, key::Number, scancode::Number, action::Number, mods::Number)
  if key == 70 && action == 1 # f
    # ...
  elseif key == 71 && action == 1 # g
    #global wireframe
    #wireframe = !wireframe
    #println("Wireframe: ",wireframe)

  elseif key == 72 && action == 1 # h
    #global fullscreen=!fullscreen

  elseif key == 87 #w
    global keyFB = (action > 0)?1:0
  elseif key == 83 #s
    global keyFB = (action > 0)?-1:0
  elseif key == 65 #a
    global keyLR = (action > 0)?-1:0
  elseif key == 68 #d
    global keyLR = (action > 0)?1:0
  elseif key == 32 #space
    global keyUD = (action > 0)?1:0
  elseif key == 67 || key == 341 # c || lctrl
    global keyUD = (action > 0)?-1:0
  elseif key == 340 #lshift
    global speed = (action > 0)
  end

  if key == 75 && action == 1 #k
    global VIEW_KEYS = !VIEW_KEYS
    println("VIEW_KEYS = $VIEW_KEYS")
  end

  if VIEW_KEYS println("K: ", key, " A: ", action, " C: ", scancode, " M: ", mods) end

  if action == 1 global keyPressed = true end
  global keyPressing = action != 0
  global keyValue = key
  nothing
end

"""
event which catches mouse key inpits and hides/shows cursor when mouse button is pressed
"""
function OnMouseKey(window, key::Number, action::Number, mods::Number)
  global cursorPos_old, cursorPos
  global mouseKeyPressed = action != 0
  if action == 1
    cursorPos_old = cursorPos
    GLFW.SetInputMode(window, GLFW.CURSOR, GLFW.CURSOR_DISABLED)
  end
  if action == 0
    GLFW.SetInputMode(window, GLFW.CURSOR, GLFW.CURSOR_NORMAL)
  end
  nothing
end

"""
event which catches mouse position for camera rotation event
"""
function OnCursorPos(window, x::Number, y::Number)
  global CAMERA, lastCursorPos, cursorPos, mouseMove, mouseKeyPressed

  t = [x/800f0,y/600f0]
  cursorPos = t
  l = lastCursorPos; if l == [] l = t end
  n = l - t
  lastCursorPos = t
  mouseMove = n

  if mouseKeyPressed OnRotate(CAMERA) end
  nothing
end

"""
rotates camera
"""
function rotate(camera::Camera, rotation::AbstractArray)
  camera.rotation += rotation
  camera.moved=true
  #println("rotated")
end

"""
moves camera, adds vector to current position
"""
function move(camera::Camera, position::AbstractArray)
  camera.position += position
  camera.moved=true
  #println("moved")
end

"""
event which calculates cursor position shifts and calls rotate function
"""
function OnRotate(camera::Camera)
  global cursorPos, cursorPos_old
  mx = cursorPos[1] - cursorPos_old[1]
  my = cursorPos[2] - cursorPos_old[2]
  #println(mx)
  cursorPos_old = cursorPos
  rotate(camera, [-mx*2f0,my*2f0,0f0]) #[-mx*2,my*2,0f0] #Vec3f((-mx+0.5f0),(-0.5f0+my),0f0)
end

oldposition = CAMERA.position
shiftposition = [0.0f0,0,0]

"""
sets camera position
"""
function setPosition(camera::Camera, position::AbstractArray)
  camera.position = position
  global oldposition = position
  camera.translateMat = translation(camera.position)
  camera.moved = true
end

"""
event which updates positions shifts (left,right,up,down,forward,back)
key is (left,right,up,down,forward,back)
m is direction value with weight (positive, negative)
"""
function OnMove(camera::Camera, key::Symbol, m::Number)
  global shiftposition
  if key == :FORWARD  move(camera, forward(camera)*(m*0.05f0*(!speed?1f0:10f0)))
  elseif key == :RIGHT  move(camera, right(camera)*(-m*0.05f0*(!speed?1f0:10f0))) #+Vec3f(-right*0.02f0,-up*0.02f0,forward*0.02f0)
  elseif key == :UP  move(camera, up(camera)*(-m*0.05f0*(!speed?1f0:10f0)))
  end

  #dif = camera.position
  #step = trunc.(dif)
  #shiftposition = dif - step
  #step = abs.(step)

  #if step[1] >= 2 camera.position[1] = 0 end
  #if step[2] >= 2 camera.position[2] = 0 end
  #if step[3] >= 2 camera.position[3] = 0 end

  #  dist = sqrt(camera.position[1]^2+camera.position[2]^2+camera.position[3]^2)
  # if dist >= 1 camera.position = [0f0,0,0] end
end

"""
update function where camera translation is update only when camera was moved by input.
here cameras MVP Matrix is updated aswell
"""
function Update(camera::Camera)
  if OnTime(0.01)
    if keyFB != 0 OnMove(camera, :FORWARD, keyFB) end
    if keyLR != 0 OnMove(camera, :RIGHT, keyLR) end
    if keyUD != 0 OnMove(camera, :UP, keyUD) end
  end

  if camera.moved
    camera.translateMat = translation(camera.position)
    camera.rotationMat = computeRotation(camera.rotation)
    camera.viewMat = camera.scalingMat * camera.rotationMat * camera.translateMat
    camera.MVP = camera.modelMat * camera.projectionMat * camera.viewMat
  end
end

"""
event which is called by game loop and calls real update function
this event resets camera moved state
"""
function OnUpdate(camera::Camera)
  #if isFocus return end
  Update(camera)
  r = camera.moved
  camera.moved = false
  r
end

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
## ------------------------------------------------------------------
# Meta

get_sec(l::RBSymLink) = get_parent(l)
set_sec!(l::RBSymLink, sec::RBSection) = set_parent!(l, sec)

## ------------------------------------------------------------------
# Data

get_text(l::RBSymLink) = getproperty!(l, :text, "")
set_text!(l::RBSymLink, text::String) = setproperty!(l, :text, strip(text))

## ------------------------------------------------------------------
# show
function Base.show(io::IO, p::RBSymLink)
    print(io, "RBSymLink(\"", get_label(p), "\")")
    print(io, "\n - path: '", localrelpath(p), "'")
    text = get_text(p)
    !isempty(text) && print(io, "\n - text: \"", _preview(io, text), "\"")
    tags = join(get_tags(p), ", ")
    !isempty(tags) && print(io, "\n - tags: \"", _preview(io, tags), "\"")
    print(io, "\n", _preview(io, "-"^70))
end


text:
# ---
# title: 223. Rectangle Area
# id: problem223
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/rectangle-area/description/>
# hidden: true
# ---
#
# Find the total area covered by two **rectilinear** rectangles in a **2D**
# plane.
#
# Each rectangle is defined by its bottom left corner and top right corner as
# shown in the figure.
#
# ![Rectangle
# Area](https://assets.leetcode.com/uploads/2018/10/22/rectangle_area.png)
#
# **Example:**
#
#
#
#     Input: A = -3, B = 0, C = 3, D = 4, E = 0, F = -1, G = 9, H = 2
#     Output: 45
#
# **Note:**
#
# Assume that the total area is never beyond the maximum possible value of
# **int**.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module DiscriminantAnalysis
    using StatsBase, LinearAlgebra

    const LA = LinearAlgebra

    # Common functions across classifiers
    include("common.jl")
    include("whiten.jl")

    # Classifier-specific code
    """
        DiscriminantModel

    Abstract type representing a discriminant model
    """
    abstract type DiscriminantModel{T<:AbstractFloat} <: StatsBase.StatisticalModel end

    include("parameters.jl")
    for classifier in ["linear"] #, "quadratic"]
        include("classifiers/$(classifier).jl")
    end
    include("discriminants.jl")
end


text:
module mt4gen1
using StaticArrays
using MeshCore: nshapes
using MeshSteward: vtkwrite
using MeshSteward: T4blockx
using Test
function test()
    connectivity = T4blockx([0.0, 1.0, 3.0], [0.0, 1.0, 3.0], [0.0, 1.0, 3.0], :a)
    @test (nshapes(connectivity.right), nshapes(connectivity.left)) == (27, 48)
    vtkwrite("mt4gen1", connectivity)
    try rm("mt4gen1.vtu"); catch end
    true
end
end
using .mt4gen1
mt4gen1.test()

module mt4gen2
using StaticArrays
using MeshCore: nshapes
using MeshSteward: vtkwrite
using MeshSteward: T4block
using Test
function test()
    connectivity = T4block(1.0, 2.0, 3.0, 7, 9, 13, :a)
    # @show (nshapes(connectivity.right), nshapes(connectivity.left))
    @test (nshapes(connectivity.right), nshapes(connectivity.left)) == (1120, 4914)
    vtkwrite("mt4gen2", connectivity)
    try rm("mt4gen2.vtu"); catch end
    true
end
end
using .mt4gen2
mt4gen2.test()

module mt4gen3
using StaticArrays
using MeshCore: nshapes
using MeshSteward: vtkwrite
using MeshSteward: T4blockx, gradedspace
using Test
function test()
    connectivity = T4blockx([1.0, 2.0, 3.0], [1.0, 2.0, 3.0], gradedspace(0.0, 5.0, 7, 2), :a)
    # @show (nshapes(connectivity.right), nshapes(connectivity.left))
    @test (nshapes(connectivity.right), nshapes(connectivity.left)) == (63, 144)
    vtkwrite("mt4gen3", connectivity)
    try rm("mt4gen3.vtu"); catch end
    true
end
end
using .mt4gen3
mt4gen3.test()

module mt4gen4
using StaticArrays
using MeshCore: nshapes
using MeshSteward: vtkwrite
using MeshSteward: T4blockx, linearspace
using Test
function test()
    connectivity = T4blockx([1.0, 2.0, 3.0], [1.0, 2.0, 3.0], linearspace(0.0, 5.0, 7), :a)
    # @show (nshapes(connectivity.right), nshapes(connectivity.left))
    @test (nshapes(connectivity.right), nshapes(connectivity.left)) == (63, 144)
    vtkwrite("mt4gen4", connectivity)
    try rm("mt4gen4.vtu"); catch end
    true
end
end
using .mt4gen4
mt4gen4.test()

module mt4gen4b
using StaticArrays
using MeshCore: nshapes
using MeshSteward: vtkwrite
using MeshSteward: T4blockx, linearspace
using Test
function test()
    connectivity = T4blockx([1.0, 2.0, 3.0], [1.0, 2.0, 3.0], linearspace(0.0, 5.0, 7), :b)
    # @show (nshapes(connectivity.right), nshapes(connectivity.left))
    @test (nshapes(connectivity.right), nshapes(connectivity.left)) == (63, 144)
    vtkwrite("mt4gen4", connectivity)
    try rm("mt4gen4.vtu"); catch end
    true
end
end
using .mt4gen4b
mt4gen4b.test()

module mt4gen4ca
using StaticArrays
using MeshCore: nshapes
using MeshSteward: vtkwrite
using MeshSteward: T4blockx, linearspace
using Test
function test()
    connectivity = T4blockx([1.0, 2.0, 3.0], [1.0, 2.0, 3.0], linearspace(0.0, 5.0, 7), :ca)
    # @show (nshapes(connectivity.right), nshapes(connectivity.left))
    @test (nshapes(connectivity.right), nshapes(connectivity.left)) == (63, 120)
    vtkwrite("mt4gen4", connectivity)
    try rm("mt4gen4.vtu"); catch end
    true
end
end
using .mt4gen4ca
mt4gen4ca.test()

module mt4gen4cb
using StaticArrays
using MeshCore: nshapes
using MeshSteward: vtkwrite
using MeshSteward: T4blockx, linearspace
using Test
function test()
    connectivity = T4blockx([1.0, 2.0, 3.0], [1.0, 2.0, 3.0], linearspace(0.0, 5.0, 7), :cb)
    # @show (nshapes(connectivity.right), nshapes(connectivity.left))
    @test (nshapes(connectivity.right), nshapes(connectivity.left)) == (63, 120)
    vtkwrite("mt4gen4", connectivity)
    try rm("mt4gen4.vtu"); catch end
    true
end
end
using .mt4gen4cb
mt4gen4cb.test()


text:
function p53(data::Dict)

  # Setup basic dimensions of arrays

  # Parse & check FEdict data

  if :struc_el in keys(data)
    struc_el = data[:struc_el]
  else
    println("No fin_el type specified.")
    return
  end

  ndim = struc_el.ndim
  nst = struc_el.nst

  fin_el = struc_el.fin_el
  @assert typeof(fin_el) <: FiniteElement

  if typeof(fin_el) == Line
    (nels, nn) = mesh_size(fin_el, struc_el.nxe)
  elseif typeof(fin_el) == Triangle || typeof(fin_el) == Quadrilateral
    (nels, nn) = mesh_size(fin_el, struc_el.nxe, struc_el.nye)
  elseif typeof(fin_el) == Hexahedron
    (nels, nn) = mesh_size(fin_el, struc_el.nxe, struc_el.nye, struc_el.nze)
  else
    println("$(typeof(fin_el)) is not a known finite element.")
    return
  end

  nodof = fin_el.nodof           # Degrees of freedom per node
  ndof = fin_el.nod * nodof      # Degrees of freedom per fin_el

  # Update penalty if specified in FEdict

  penalty = 1e20
  if :penalty in keys(data)
    penalty = data[:penalty]
  end

  # Allocate all arrays

  # Start with arrays to be initialized from FEdict

  if :properties in keys(data)
    prop = zeros(size(data[:properties], 1), size(data[:properties], 2))
    for i in 1:size(data[:properties], 1)
      prop[i, :] = data[:properties][i, :]
    end
  else
    println("No :properties key found in FEdict")
  end

  nf = ones(Int, nodof, nn)
  if :support in keys(data)
    for i in 1:size(data[:support], 1)
      nf[:, data[:support][i][1]] = data[:support][i][2]
    end
  end

  x_coords = zeros(nn)
  if :x_coords in keys(data)
    x_coords = data[:x_coords]
  end

  y_coords = zeros(nn)
  if :y_coords in keys(data)
    y_coords = data[:y_coords]
  end

  z_coords = zeros(nn)
  if :z_coords in keys(data)
    z_coords = data[:z_coords]
  end

  etype = ones(Int, nels)
  if :etype in keys(data)
    etype = data[:etype]
  end

  # All other arrays

  points = zeros(struc_el.nip, ndim)
  g = zeros(Int, ndof)
  g_coord = zeros(ndim,nn)
  fun = zeros(fin_el.nod)
  coord = zeros(fin_el.nod, ndim)
  gamma = zeros(nels)
  jac = zeros(ndim, ndim)
  g_num = zeros(Int, fin_el.nod, nels)
  der = zeros(ndim, fin_el.nod)
  deriv = zeros(ndim, fin_el.nod)
  bee = zeros(nst,ndof)
  km = zeros(ndof, ndof)
  mm = zeros(ndof, ndof)
  gm = zeros(ndof, ndof)
  kg = zeros(ndof, ndof)
  eld = zeros(ndof)
  weights = zeros(struc_el.nip)
  g_g = zeros(Int, ndof, nels)
  num = zeros(Int, fin_el.nod)
  actions = zeros(ndof, nels)
  displacements = zeros(size(nf, 1), ndim)
  gc = ones(ndim)
  dee = zeros(nst,nst)
  sigma = zeros(nst)
  axial = zeros(nels)

  formnf!(nodof, nn, nf)
  neq = maximum(nf)
  println("There are $(neq) equations.\n")

  # Find global array sizes

  for iel in 1:nels
    hexahedron_xz!(iel, x_coords, y_coords, z_coords, coord, num)
    num_to_g!(num, nf, g)
    g_num[:, iel] = num
    g_coord[:, num] = coord'
    g_g[:, iel] = g
  end

  sample!(fin_el, points, weights)

  gsm = spzeros(neq, neq)
  for iel in 1:nels
    deemat!(dee, prop[etype[iel], 1], prop[etype[iel], 2])
    num = g_num[:, iel]
    coord = g_coord[:, num]'              # Transpose
    g = g_g[:, iel]
    km = zeros(ndof, ndof)
    for i in 1:struc_el.nip
      shape_der!(der, points, i)
      jac = der*coord
      detm = det(jac)
      jac = inv(jac)
      deriv = jac*der
      beemat!(bee, deriv)
      km += (bee')*dee*bee*detm*weights[i]
    end
    fsparm!(gsm, g, km)
  end

  loads = OffsetArray(zeros(neq + 1), 0:neq)
  if :loaded_nodes in keys(data)
    for i in 1:size(data[:loaded_nodes], 1)
      loads[nf[:, data[:loaded_nodes][i][1]]] = data[:loaded_nodes][i][2]
    end
  end

  fixed_freedoms = 0
  if :fixed_freedoms in keys(data)
    fixed_freedoms = size(data[:fixed_freedoms], 1)
  end
  no = zeros(Int, fixed_freedoms)
  node = zeros(Int, fixed_freedoms)
  sense = zeros(Int, fixed_freedoms)
  value = zeros(Float64, fixed_freedoms)
  if :fixed_freedoms in keys(data) && fixed_freedoms > 0
    for i in 1:fixed_freedoms
      no[i] = nf[data[:fixed_freedoms][i][2], data[:fixed_freedoms][i][1]]
      value[i] = data[:fixed_freedoms][i][3]
      gsm[no[i], no[i]] += penalty
      loads[no[i]] = gsm[no[i], no[i]] * value[i]
    end
  end

  cfgsm = cholesky(gsm)
  loads[1:neq] = cfgsm \ loads[1:neq]

  displacements = zeros(size(nf))
  for i in 1:size(displacements, 1)
    for j in 1:size(displacements, 2)
      if nf[i, j] > 0
        displacements[i,j] = loads[nf[i, j]]
      end
    end
  end
  displacements = displacements'

  dis_df = DataFrame(
    r_disp = displacements[:, 1],
    y_disp = displacements[:, 2],
    z_disp = displacements[:, 3]
  )

  struc_el.nip = 1
  points = zeros(struc_el.nip, ndim)
  weights = zeros(struc_el.nip)

  sample!(fin_el, points, weights)

  gc1 = Vector{Float64}()
  gc2 = Vector{Float64}()
  s1 = Vector{Float64}()
  s2 = Vector{Float64}()
  s3 = Vector{Float64}()
  t1 = Vector{Float64}()
  t2 = Vector{Float64}()
  t3 = Vector{Float64}()

  for iel in 1:nels
    deemat!(dee, prop[etype[iel], 1], prop[etype[iel], 2])
    num = g_num[:, iel]
    coord = g_coord[:, num]'
    g = g_g[:, iel]
    eld = loads[g]
    for i in 1:struc_el.nip
      shape_fun!(fun, points, i)
      shape_der!(der, points, i)
      gc = fun'*coord
      jac = inv(der*coord)
      deriv = jac*der
      beemat!(bee, deriv)
      sigma = dee*(bee*eld)
      gc1 = append!(gc1, gc[1])
      gc2 = append!(gc2, gc[2])
      s1 = append!(s1, sigma[1])
      s2 = append!(s2, sigma[2])
      s3 = append!(s3, sigma[3])
      t1 = append!(t1, sigma[4])
      t2 = append!(t2, sigma[5])
      t3 = append!(t3, sigma[6])
    end
  end

  sigma_df = DataFrame(
    r_coord = gc1,
    y_coord = gc2,
    sig_x = s1,
    sig_y = s2,
    sig_z = s3,
    tau_xy = t1,
    tau_yz = t2,
    tau_zx = t3
  )

  fem = PtFEM.jFEM(struc_el, fin_el, ndim, nels, nst, ndof, nn, nodof,
    neq, penalty, etype, g, g_g, g_num, nf, no,
    node, num, sense, actions, bee, coord, gamma, dee,
    der, deriv, displacements, eld, fun, gc, g_coord, jac,
    km, mm, kg, cfgsm, loads, points, prop, sigma, value,
    weights, x_coords, y_coords, z_coords, axial)

    (fem, dis_df, sigma_df)
  end

max_stars_repo_path: src/5 Elastic Solids/p53.jl

max_issues_repo_path: src/5 Elastic Solids/p53.jl

max_forks_repo_path: src/5 Elastic Solids/p53.jl


text:
#####
##### fallback rules
#####

#=
There's an idea at play here that's not made explicit by these default "no
primitive available" fallback implementations.

In some weird ideal sense, the fallback for e.g. `forward_rule` should actually
be "get the derivative via forward-mode AD". This is necessary to enable
mixed-mode rules, where e.g.  `forward_rule` is used within a `reverse_rule`
definition. For example, broadcasted functions may not themselves be
forward-mode *primitives*, but are often forward-mode *differentiable*.

ChainRules, by design, is decoupled from any specific AD implementation. How,
then, do we know which AD to fall back to when there isn't a primitive defined?

Well, if you're a greedy AD implementation, you can just overload `forward_rule`
and/or `reverse_rule` to use your AD directly. However, this won't place nice
with other AD packages doing the same thing, and thus could cause
load-order-dependent problems for downstream users.

It turns out, Cassette solves this problem nicely by allowing AD authors to
overload the fallbacks w.r.t. their own context. Example using ForwardDiff:

```
using ChainRules, ForwardDiff, Cassette

Cassette.@context MyChainRuleCtx

# ForwardDiff, itself, can call `my_forward_rule` instead of
# `forward_rule` to utilize the ForwardDiff-injected ChainRules
# infrastructure
my_forward_rule(args...) = Cassette.overdub(MyChainRuleCtx(), forward_rule, args...)

function Cassette.execute(::MyChainRuleCtx, ::typeof(forward_rule)
                          ::@sig(R → R), f, x)
    fx, df = forward_rule(sig, f, x)
    if isa(df, Nothing)
        fx, df = (f(x), ẋ -> ẋ * ForwardDiff.derivative(f, x))
    end
    return fx, df
end
```
=#

forward_rule(::Signature, ::Vararg{Any}) = (nothing, nothing)

reverse_rule(::Signature, ::Vararg{Any}) = (nothing, nothing)

# TODO: Should the default be to whitelist known holomorphic functions, or to
# blacklist known non-holomorphic functions? This implements the latter.
function forward_rule(signature::@sig(C → C), f, x)
    fx, df = forward_rule(Signature(RealScalar(), RealScalar()), f, x)
    return fx, ẋ -> (df(ẋ), false)
end

function reverse_rule(signature::@sig(R → R), f, x)
    fx, df = forward_rule(signature, f, x)
    return fx, (x̄, z̄) -> reverse_chain!(x̄, @thunk(df(z̄)))
end

#####
##### `@forward_rule`
#####

"""
    @forward_rule(signature, f(x, y, ...),
                  (df₁_dx, df₁_dy, ...),
                  (df₂_dx, df₂_dy, ...),
                  ...)

Define a method for `forward_rule` using the given signature, call expression,
and derivative expressions.

While this macro is convenient for avoiding boilerplate code when implementing
simple forward rules, note that more advanced rules will probably require
overloading `ChainRules.forward_rule` directly. For now, the macro only
supports real-domain rules.

Examples:

`@forward_rule(R → R, sin(x), cos(x))` expands to:

    function forward_rule(::@sig(R → R), ::typeof(sin), x)
        return sin(x), ẋ -> forward_chain(@thunk(cos(x)), ẋ)
    end

`@forward_rule(R⊕R → R, *(x, y), (y, x))` expands to:

    function forward_rule(::@sig(R → R), ::typeof(*), x, y)
        return *(x, y), (ẋ, ẏ) -> forward_chain(@thunk(y), ẋ, @thunk(x), ẏ)
    end

`@forward_rule(R⊕R → R, sincos(x), cos(x), -sin(x))` expands to:

    function forward_rule(::@sig(R → R⊕R), ::typeof(sincos), x)
        return sincos(x),
               (ẋ -> forward_chain(ẋ, @thunk(cos(x))),
                ẋ -> forward_chain(ẋ, @thunk(-sin(x))))
    end

Note that this last case is a good example of a primitive that is more
efficiently implemented with a manual `forward_rule` overload:

    function forward_rule(::@sig(R → R⊕R), ::typeof(sincos), x)
        sinx, cosx = sincos(x)
        return (sinx, cosx),
               (ẋ -> forward_chain(ẋ, @thunk(cosx)),
                ẋ -> forward_chain(ẋ, @thunk(-sinx)))
    end
"""
macro forward_rule(signature, call, derivs...)
    return generate_rule_definition(signature, call, nothing, derivs...)
end

"""
    @reverse_rule(signature, f(x, y, ...), (f̄₁, f̄₂, ...), df_dx, df_dy, ...)

Define a method for `reverse_rule` using the given signature, call expression,
and derivative expressions.

The third argument to the macro is a tuple of symbols naming the adjoints of the
outputs of `f(x, y, ...)` (or just a single symbol if there is only one output).

While this macro is convenient for avoiding boilerplate code when implementing
simple reverse rules, note that more advanced rules will probably require
overloading `ChainRules.reverse_rule` directly. For now, this macro only
supports real-domain rules.

Examples:

`@reverse_rule([R] → R, sum(x), ȳ, ȳ)` expands to:

    function reverse_rule(::@sig([R] → R), ::typeof(sum), x)
        return sum(x), (x̄, ȳ) -> reverse_chain!(x̄, @thunk(ȳ))
    end

`@reverse_rule([R]⊕[R] → [R], *(x, y), z̄, z̄ * y', x' * z̄)` expands to:

    function reverse_rule(::@sig([R]⊕[R] → R), ::typeof(*), x, y)
        return x * y,
               ((x̄, z̄) -> reverse_chain!(x̄, @thunk(z̄ * y')),
                (ȳ, z̄) -> reverse_chain!(ȳ, @thunk(x' * z̄)))
    end
"""
macro reverse_rule(signature, call, adjoint_names, derivs...)
    if Meta.isexpr(adjoint_names, :tuple)
        adjoint_names = convert(Vector{Symbol}, adjoint_names.args)
    elseif isa(adjoint_names, Symbol)
        adjoint_names = Symbol[adjoint_names]
    end
    return generate_rule_definition(signature, call, adjoint_names, derivs...)
end

# TODO: Expand this beyond real-domain rules by parsing signature ahead of time
function generate_rule_definition(signature, call,
                                  adjoint_names::Union{Nothing,Vector{Symbol}},
                                  derivs...)
    @assert Meta.isexpr(call, :call)
    call_function = esc(call.args[1])
    call_args = map(esc, call.args[2:end])
    seed_names = Any[Symbol(string(:seed_, i)) for i in 1:length(call_args)]
    chains = Any[]
    if isa(adjoint_names, Nothing) # we're doing forward mode
        for deriv in derivs
            thunkables = Meta.isexpr(deriv, :tuple) ? deriv.args : [deriv]
            thunks = [:(@thunk($(esc(t)))) for t in thunkables]
            @assert length(thunks) == length(call_args)
            chain_call = Expr(:call, :forward_chain)
            for i in 1:length(thunks)
                push!(chain_call.args, seed_names[i])
                push!(chain_call.args, thunks[i])
            end
            push!(chains, :($(Expr(:tuple, seed_names...)) -> $chain_call))
        end
    else # we're doing reverse mode
        @assert length(derivs) == length(call_args)
        @assert length(adjoint_names) > 0
        adjoint_names = map(esc, adjoint_names)
        for i in 1:length(derivs)
            seed_name = seed_names[i]
            chain_call = Expr(:call, :reverse_chain!)
            push!(chain_call.args, seed_name)
            push!(chain_call.args, :(@thunk($(esc(derivs[i])))))
            push!(chains, :($(Expr(:tuple, seed_name, adjoint_names...)) -> $chain_call))
        end
    end
    @assert length(chains) > 0
    chains = length(chains) > 1 ? Expr(:tuple, chains...) : chains[1]
    rule_function = isa(adjoint_names, Nothing) ? :forward_rule : :reverse_rule
    return quote
        @assert(@sig($signature) <: Signature{<:Tuple{Vararg{Union{RealTensor,RealScalar,Ignore}}},
                                              <:Tuple{Vararg{Union{RealTensor,RealScalar,Ignore}}}},
                "@forward_rule and @reverse_rule only support real-domain rules right now")
        function $ChainRules.$rule_function(::@sig($signature),
                                            ::typeof($call_function),
                                            $(call_args...))
            outputs = $(esc(call))
            return outputs, $chains
        end
    end
end

#####
##### forward rules
#####

@forward_rule(R → R, sin(x), cos(x))
@forward_rule(R → R, cos(x), -sin(x))
@forward_rule(R → R, log(x), inv(x))
@forward_rule(R⊕R → R, *(x, y), (y, x))

function forward_rule(::@sig(R⊕R → R), ::typeof(atan), y, x)
    h = hypot(y, x)
    return atan(y, x), (ẏ, ẋ) -> forward_chain(ẏ, @thunk(x / h), ẋ, @thunk(y / h))
end

function forward_rule(::@sig(R⊕R → R), ::typeof(hypot), x, y)
    h = hypot(x, y)
    return h, (ẋ, ẏ) -> forward_chain(ẋ, @thunk(x / h), ẏ, @thunk(y / h))
end

function forward_rule(::@sig(R → R⊕R), ::typeof(sincos), x)
    sinx, cosx = sincos(x)
    return (sinx, cosx),
           (ẋ -> forward_chain(ẋ, @thunk(cosx)),
            ẋ -> forward_chain(ẋ, @thunk(-sinx)))
end

forward_rule(::@sig(C → C), ::typeof(conj), x) = conj(x), ẋ -> (false, true)

#####
##### reverse rules
#####

@reverse_rule([R] → R, sum(x), ȳ, ȳ)
@reverse_rule([R]⊕[R] → [R], +(x, y), z̄, z̄, z̄)
@reverse_rule([R]⊕[R] → [R], *(x, y), z̄, z̄ * y', x' * z̄)

# TODO: This partial derivative extraction should be doable without the extra
# temporaries or preallocation utilized here, but AFAICT such an approach is
# hard to write without relying on inference hacks unless we have something
# akin to https://github.com/JuliaLang/julia/issues/22129
function reverse_rule(::@sig(_⊕[R] → [R]), ::typeof(broadcast), f, x)
    s = Signature(RealScalar(), RealScalar())
    f_rule = x -> begin
        y, d = forward_rule(s, f, x)
        y, d(one(x))
    end
    applied_f_rule = broadcast(f_rule, x)
    values = map(first, applied_f_rule)
    derivs = map(last, applied_f_rule)
    return values, (x̄, z̄) -> reverse_chain!(x̄, @thunk(broadcasted(*, z̄, derivs)))
end


text:
module deSolveDiffEq

using Reexport, RCall
@reexport using DiffEqBase

solver = Ref{Module}()

function __init__()
    solver[] = rimport("deSolve")
end

abstract type deSolveAlgoritm <: DiffEqBase.AbstractODEAlgorithm end

struct lsoda <: deSolveAlgoritm end
struct lsode <: deSolveAlgoritm end
struct lsodes <: deSolveAlgoritm end
struct lsodar <: deSolveAlgoritm end
struct vode <: deSolveAlgoritm end
struct daspk <: deSolveAlgoritm end
struct euler <: deSolveAlgoritm end
struct rk4 <: deSolveAlgoritm end
struct ode23 <: deSolveAlgoritm end
struct ode45 <: deSolveAlgoritm end
struct radau <: deSolveAlgoritm end
struct bdf <: deSolveAlgoritm end
struct bdf_d <: deSolveAlgoritm end
struct adams <: deSolveAlgoritm end
struct impAdams <: deSolveAlgoritm end
struct impAdams_d <: deSolveAlgoritm end
struct iteration <: deSolveAlgoritm end

function DiffEqBase.__solve(
    prob::DiffEqBase.AbstractODEProblem,
    alg::deSolveAlgoritm,timeseries=[],ts=[],ks=[];
    saveat=eltype(prob.tspan)[],timeseries_errors=true,
    reltol = 1e-3, abstol = 1e-6,
    maxiters = 100000,
    kwargs...)

    p = prob.p
    tspan = prob.tspan
    u0 = prob.u0

    if DiffEqBase.isinplace(prob)
        f = function (t,u,__p)
            du = similar(u)
            prob.f(du,u,p,t)
            R"list($du)" # Error message says a list return is required!
        end
    else
        f = function (t,u,__p)
            du = prob.f(u,p,t)
            R"list($du)" # Error message says a list return is required!
        end
    end

    _saveat = isempty(saveat) ? nothing : saveat
    if _saveat isa Array
        __saveat = _saveat
    elseif _saveat isa Number
        __saveat = Array(tspan[1]:_saveat:tspan[2])
    elseif _saveat isa Nothing
        __saveat = [tspan[1],tspan[2]]
    else
        __saveat = Array(_saveat)
    end

    out = rcopy(solver[].ode(times = __saveat, y = u0, func = f,
                             method = string(alg)[15:end-2],
                             parms=nothing, maxsteps = maxiters,
                             rtol = reltol, atol = abstol))

    ts = @view out[:,1]

    if typeof(u0) <: AbstractArray
        timeseries = Vector{typeof(u0)}(undef,length(ts))
        for i=1:length(ts)
            timeseries[i] = @view out[i,2:end]
        end
    else
        timeseries = out[i,end]
    end

    DiffEqBase.build_solution(prob,alg,ts,timeseries,
                              dense = false,
                              timeseries_errors = timeseries_errors)
end

end # module


text:
module Marker

using ..LIKWID: LibLikwid, PerfMon, capture_stderr!

"""
Initialize the Marker API. Must be called previous to all other functions.
"""
function init()
    LibLikwid.likwid_markerInit()
    Threads.@threads :static for i in 1:Threads.nthreads()
        LibLikwid.likwid_markerThreadInit()
    end
    return nothing
end

"""
Initialize the Marker API only on the main thread. `LIKWID.Marker.threadinit()` must be called manually.
"""
init_nothreads() = LibLikwid.likwid_markerInit()

"""
Add the current thread to the Marker API.
"""
threadinit() = LibLikwid.likwid_markerThreadInit()

"""
Register a region with name `regiontag` to the Marker API. On success, `true` is returned.

This is an optional function to reduce the overhead of region registration at `Marker.startregion`.
If you don't call `registerregion`, the registration is done at `startregion`.
"""
function registerregion(regiontag::AbstractString)
    ret = LibLikwid.likwid_markerRegisterRegion(regiontag)
    return ret == 0
end

"""
Start measurements under the name `regiontag`. On success, `true` is returned.
"""
function startregion(regiontag::AbstractString)
    ret = LibLikwid.likwid_markerStartRegion(regiontag)
    return ret == 0
end

"""
Stop measurements under the name `regiontag`. On success, `true` is returned.
"""
function stopregion(regiontag::AbstractString)
    ret = LibLikwid.likwid_markerStopRegion(regiontag)
    return ret == 0
end

"""
    getregion(regiontag::AbstractString, [num_events]) -> nevents, events, time, count
Get the intermediate results of the region identified by `regiontag`. On success, it returns
    * `nevents`: the number of events in the current group,
    * `events`: a list with all the aggregated event results,
    * `time`: the measurement time for the region and
    * `count`: the number of calls.
"""
function getregion(regiontag::AbstractString)
    current_group = PerfMon.get_id_of_active_group()
    return getregion(regiontag, PerfMon.get_number_of_events(current_group))
end

function getregion(regiontag::AbstractString, num_events::Integer)
    nevents = Ref(Int32(num_events))
    events = zeros(nevents[])
    time = Ref(0.0)
    count = Ref(Int32(0))
    LibLikwid.likwid_markerGetRegion(regiontag, nevents, events, time, count)
    return nevents[], events, time[], count[]
end

"""
Switch to the next event set in a round-robin fashion.
If you have set only one event set on the command line, this function performs no operation.
"""
nextgroup() = LibLikwid.likwid_markerNextGroup()

"""
Reset the values stored using the region name `regiontag`.
On success, `true` is returned.
"""
function resetregion(regiontag::AbstractString)
    ret = LibLikwid.likwid_markerResetRegion(regiontag)
    return ret == 0
end

"""
Close the connection to the LIKWID Marker API and write out measurement data to file.
This file will be evaluated by `likwid-perfctr`.
"""
close() = LibLikwid.likwid_markerClose()

"""
Checks whether the Marker API is active, i.e. julia has been started under `likwid-perfctr -C ... -g ... -m`.
"""
isactive() = !isnothing(get(ENV, "LIKWID_MODE", nothing))

"""
    region(f, regiontag::AbstractString)
Adds a LIKWID marker region around the execution of the given function `f` using [`Marker.startregion`](@ref),
[`Marker.stopregion`](@ref) under the hood.
Note that `LIKWID.Marker.init()` and `LIKWID.Marker.close()` must be called before and after, respectively.

# Examples
```julia
julia> using LIKWID

julia> Marker.init()

julia> region("sleeping...") do
           sleep(1)
       end
true

julia> region(()->rand(100), "create rand vec")
true

julia> Marker.close()

```
"""
function region(f, regiontag::AbstractString)
    Marker.startregion(regiontag)
    f()
    return Marker.stopregion(regiontag)
end

"""
Convenience macro for flanking code with [`Marker.startregion`](@ref) and [`Marker.stopregion`](@ref).

# Examples
```julia
julia> using LIKWID

julia> Marker.init()

julia> @region "sleeping..." sleep(1)
true

julia> @region "create rand vec" rand(100)
true

julia> Marker.close()

```
"""
macro region(regiontag, expr)
    q = quote
        LIKWID.Marker.startregion($regiontag)
        $(expr)
        LIKWID.Marker.stopregion($regiontag)
    end
    return esc(q)
end

end # module


text:
include("/opt/bcbi/predictmd/OfflineRegistry/startup.jl");

import Pkg;

offline_registry_directory = "/opt/bcbi/predictmd/OfflineRegistry/"
offline_registry_name = Pkg.TOML.parsefile(
    joinpath(splitpath(offline_registry_directory)..., "Registry.toml")
    )["name"];
offline_registry_uuid = Pkg.TOML.parsefile(
    joinpath(splitpath(offline_registry_directory)..., "Registry.toml")
    )["uuid"];
pushfirst!(Base.DEPOT_PATH,joinpath(splitpath(@__DIR__)..., "depot",),);
pushfirst!(
    Base.DEPOT_PATH,
    joinpath(
        homedir(),
        ".julia.isolated",
        offline_registry_name,
        offline_registry_uuid,
        ),
    );
filter!((x) -> !(Base.Filesystem.samefile(expanduser(x),expanduser("~/.julia")) || lowercase(strip(abspath(expanduser(x)))) == lowercase(strip(abspath(expanduser("~/.julia")))) || Base.Filesystem.samefile(x,"~/.julia") || lowercase(strip(abspath(x))) == lowercase(strip(abspath("~/.julia"))) || Base.Filesystem.samefile(expanduser(x),expanduser(joinpath(homedir(), ".julia"))) || lowercase(strip(abspath(expanduser(x)))) == lowercase(strip(abspath(expanduser(joinpath(homedir(), ".julia"))))) || Base.Filesystem.samefile(x,joinpath(homedir(), ".julia")) || lowercase(strip(abspath(x))) == lowercase(strip(abspath(joinpath(homedir(), ".julia"))))), Base.DEPOT_PATH);
unique!(Base.DEPOT_PATH);

predictmd = joinpath(
    homedir(),
    ".julia.isolated",
    offline_registry_name,
    offline_registry_uuid,
    "environments",
    "v$(VERSION.major).$(VERSION.minor)-PredictMD",
    );
predictmd_stable = joinpath(
    homedir(),
    ".julia.isolated",
    offline_registry_name,
    offline_registry_uuid,
    "environments",
    "v$(VERSION.major).$(VERSION.minor)-PredictMD-stable",
    );
predictmd_unstable = joinpath(
    homedir(),
    ".julia.isolated",
    offline_registry_name,
    offline_registry_uuid,
    "environments",
    "v$(VERSION.major).$(VERSION.minor)-PredictMD-unstable",
    );
predictmd_develop = joinpath(
    homedir(),
    ".julia.isolated",
    offline_registry_name,
    offline_registry_uuid,
    "environments",
    "v$(VERSION.major).$(VERSION.minor)-PredictMD-develop",
    );

predictmd_environments = String[
    predictmd,
    predictmd_stable,
    predictmd_unstable,
    predictmd_develop,
    ];

if !isfile(joinpath(predictmd, "Project.toml",)) ||
        !isfile(joinpath(predictmd, "Manifest.toml",))
    mkpath(predictmd);
    Pkg.activate(predictmd);
    Pkg.add(["PredictMD", "PredictMDExtra", "PredictMDFull",],);
end

if !isfile(joinpath(predictmd_stable, "Project.toml",)) ||
        !isfile(joinpath(predictmd_stable, "Manifest.toml",))
    mkpath(predictmd_stable);
    Pkg.activate(predictmd_stable);
    Pkg.add(["PredictMD", "PredictMDExtra", "PredictMDFull",],);
end

if !isfile(joinpath(predictmd_unstable, "Project.toml",)) ||
        !isfile(joinpath(predictmd_unstable, "Manifest.toml",))
    mkpath(predictmd_unstable);
    Pkg.activate(predictmd_unstable);
    Pkg.add(
        [
            Pkg.PackageSpec(name="PredictMD", rev="master",),
            Pkg.PackageSpec(name="PredictMDExtra", rev="master",),
            Pkg.PackageSpec(name="PredictMDFull", rev="master",),
            ],
        );
end

if !isfile(joinpath(predictmd_develop, "Project.toml",)) ||
        !isfile(joinpath(predictmd_develop, "Manifest.toml",))
    mkpath(predictmd_develop);
    Pkg.activate(predictmd_develop);
    Pkg.add(
        [
            Pkg.PackageSpec(name="PredictMD", rev="master",),
            Pkg.PackageSpec(name="PredictMDExtra", rev="master",),
            Pkg.PackageSpec(name="PredictMDFull", rev="master",),
            ],
        );
    Pkg.develop(
        ["PredictMD", "PredictMDExtra", "PredictMDFull",];
        shared = false,
        );
end

Pkg.activate(predictmd);


text:
# Local Geary test of Spatial Autocorrelation
struct LocalGeary <: AbstractLocalSpatialAutocorrelation
    n::Int
    C::Vector{Float64}
    p::Vector{Float64}
    Cperms::Matrix{Float64}
    Cpermsmean::Vector{Float64}
    Cpermsstd::Vector{Float64}
    z::Vector{Float64}
    q::Vector{Symbol}
    categories::Symbol
end

"""
    localgeary(x, W)
Compute the Local Geary test of spatial autocorrelation.

# Optional Arguments
- `permutations=9999`: number of permutations for the randomization test.
- `rng=default_rng()`: random number generator for the randomization test.
- `corrected=true`: divide the scaling factor by ``n-1`` instead of ``n``.
- `categories=:positivenegative`: assing observations to positive or negative spatial autocorrelation, or in combination with the `:moran` scatterplot.
"""
function localgeary(x::AbstractVector{T} where T<:Real, W::SpatialWeights; permutations::Int = 9999,
    corrected::Bool = true, categories::Symbol = :positivenegative,
    rng::AbstractRNG = default_rng())::LocalGeary

    (categories == :positivenegative) || (categories == :moran) || throw(ArgumentError("`categories` must be :positivenegative or :moran"))

    n = length(x)
    z = x .- mean(x)

    m2::Float64 = sum(z.^2)
    if corrected
        m2 = m2 ./ (n - 1)
    else
        m2 = m2 ./ n
    end

    # Auxiliar function to calculate Local Moran
    function localgeary_calc(zi::Number, wi::AbstractVector, zneighi::AbstractVector)::Float64
        return (1 ./ m2)  .* sum(wi .* (zi .- zneighi).^2 )
    end

    # Local Geary
    C = zeros(n)
    for i in 1:n
        zi = z[i]
        wi = weights(W, i)
        zneighi = z[neighbors(W, i)]
        C[i] = localgeary_calc(zi, wi, zneighi)
    end

    # Conditional randomizatoin
    Cperms = crand_local(permutations, z, W, localgeary_calc, rng)

    larger = sum(Cperms .>= repeat(C, 1, permutations), dims = 2)
    low = (permutations .- larger) .< larger
    larger[low] .= permutations .- larger[low]
    p = (larger .+ 1) ./ (permutations + 1)
    p = vec(p)

    Cpermsstd = vec(std(Cperms, dims = 2, corrected = false))
    Cpermsmean = vec(mean(Cperms, dims = 2))
    zval = (C .- Cpermsmean)  ./ Cpermsstd
    zval = vec(zval)

    # Classification
    Cmean = mean(C)
    Wz = slag(W, z)

    q = Array{Symbol}(undef, n)
    if categories == :positivenegative
        for i in 1:n
            if zval[i] < 0
                q[i] = :P
            else
                q[i] = :N
            end
        end
    elseif categories == :moran
        for i in 1:n
            if zval[i] < 0
                if (z[i] > 0) & (Wz[i] > 0)
                    q[i] = :HH
                elseif (z[i] < 0) & (Wz[i] < 0)
                    q[i] = :LL
                else
                    q[i] = :OP
                end
            else
                q[i] = :NE
            end
        end
    end

    return LocalGeary(n, C, p, Cperms, Cpermsmean, Cpermsstd, zval, q, categories)

end

score(x::LocalGeary) = x.C;

scoreperms(x::LocalGeary) = x.Cperms;

mean(x::LocalGeary) = x.Cpermsmean;

std(x::LocalGeary) = x.Cpermsstd;

zscore(x::LocalGeary) = x.z;

pvalue(x::LocalGeary) = x.p;

assignments(x::LocalGeary) = x.q;

testname(::LocalGeary) = "Local Geary";

function labelsorder(x::LocalGeary)
    if x.categories == :positivenegative
        return [:ns, :P, :N]
    elseif x.categories == :moran
        return [:ns; :HH; :LL; :OP; :NE];
    end
end

function labelsnames(x::LocalGeary)
    if x.categories == :positivenegative
        return ["Not Significant"; "Positive"; "Negative"]
    elseif x.categories == :moran
        return ["Not Significant"; "High-High"; "Low-Low"; "Other Positive"; "Negative"]
    end
end

function labelcolor(::LocalGeary, x::String)

    # :positivenegative categories
    if x == "P"
        catcolor = :red
        alpha = 1.0
    elseif x == "N"
        catcolor =:blue
        alpha = 1.0
    # :moran categories
    elseif x == "HH"
        catcolor = "#B2182B" # (178, 24, 43)
        alpha = 1.0
    elseif x == "LL"
        catcolor = "#EF8A62" # (239, 138, 98)
        alpha = 1.0
    elseif x == "OP"
        catcolor = "#FDDBC7" # (253, 219, 199)
        alpha = 1.0
    elseif x == "NE"
        catcolor = "#67ADC7" # (103, 173, 199)
        alpha = 1.0
    # not significant
    elseif x == "ns"
        catcolor = :lightgrey
        alpha = 0.4
    end

    return (catcolor, alpha)
end


text:
#=
    An implementation of the linear system solver in https://arxiv.org/abs/1605.02353 by Rasmus Kyng and Sushant Sachdeva.
    In addition to the setup in the paper, we also use a low stretch tree to approximate effective
    resistances on edges. To perform well cache wise, we implement a cache friendly list of
    linked lists - found in revampedLinkedListFloatStorage.jl

    The parameters used in the solver are described below.
=#

"""
Parameters for the sampling solver.
"""
type samplingParams{Tv,Ti}
    eps::Tv
    sampConst::Tv       # we create sampConst * log(n)^2 / eps^2 multiedges in the beginning
    beta::Tv            # we scale the tree up by beta
    startingSize::Ti    # the initial size of the linked list storage structure
    blockSize::Ti       # the size of each consecutive block of memory assigned to a certain element

    # debug parameters, we can get rid of them later
    verboseSS::Bool
    returnCN::Bool
    CNTol::Tv
end

defaultSamplingParams = samplingParams(0.5, 0.02, 1e3, 1000, 20,
                                false,false,1e-3)

"""
    An implementation of the linear system solver in https://arxiv.org/abs/1605.02353 by Rasmus Kyng and Sushant Sachdeva.
    In addition to the setup in the paper, we also use a low stretch tree to approximate effective
    resistances on edges. To perform well cache wise, we implement a cache friendly list of
    linked lists - found in revampedLinkedListFloatStorage.jl

    ~~~julia
    samplingLapSolver(a, tol, maxits, maxtime, params)
    ~~~
"""
function samplingSDDSolver{Tv,Ti}(SDDmat::SparseMatrixCSC{Tv,Ti};
                                tol::Tv=1e-6, maxits=1000, maxtime=Inf, verbose=false,
                                params::samplingParams{Tv,Ti}=defaultSamplingParams)

    # srand(1234)

    adjMat,diag = adj(SDDmat)

    a = extendMatrix(adjMat,diag)
    n = a.n

    F,gOp,_,_,ord,cn,cntime = buildSolver(a, params=params)

    if params.verboseSS
        println()
        println("Eps error: ", checkError(gOp))
        println("Condition number: ", cn)
        println()
    end

    la = lap(symPermuteCSC(a, ord))
    function f(b)
        #=
            We need to add an extra entry to b to make it match the size of a. The extra vertex in a is
            vertex n, thus, we will add the new entry in b on position n as well.
        =#
        auxb = copy(b)
        if norm(diag) != 0
            push!(auxb, -sum(auxb))
        end

        ret = pcg(la, auxb[ord], F, tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose)
        ret = ret[invperm(ord)]

        # We want to discard the nth element of ret (which corresponds to the first element in the permutation)
        ret = ret - ret[n]

        if norm(diag) != 0
            pop!(ret)
        end

        return ret
    end

    return f

end

"""
    An implementation of the linear system solver in https://arxiv.org/abs/1605.02353 by Rasmus Kyng and Sushant Sachdeva.
    In addition to the setup in the paper, we also use a low stretch tree to approximate effective
    resistances on edges. To perform well cache wise, we implement a cache friendly list of
    linked lists - found in revampedLinkedListFloatStorage.jl

    ~~~julia
    samplingLapSolver(a, tol, maxits, maxtime, params)
    ~~~
"""
function samplingLapSolver{Tv,Ti}(a::SparseMatrixCSC{Tv,Ti};
                                tol::Tv=1e-6, maxits=1000, maxtime=Inf, verbose=false,
                                params::samplingParams{Tv,Ti}=defaultSamplingParams)

    # srand(1234)

	n = a.n;

    F,gOp,_,_,ord,cn,cntime = buildSolver(a, params=params)

    if params.verboseSS
        println()
        println("Eps error: ", checkError(gOp))
        println("Condition number: ", cn)
        println()
    end

    la = lap(symPermuteCSC(a, ord))
    function f(b)
        ret = pcg(la, b[ord] - sum(b), F, tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose)
        return ret[invperm(ord)]
    end

    return f

end

# Add a new vertex to a with weights to the other vertices corresponding to diagonal surplus weight
function extendMatrix{Tv,Ti}(a::SparseMatrixCSC{Tv,Ti}, diag::Array{Tv,1})

    if norm(diag) == 0
        return a
    end

    n = a.n
    u,v,w = findnz(a)
    for i in 1:n
        if diag[i] > 0
            push!(u, i)
            push!(v, n + 1)
            push!(w, diag[i])

            push!(u, n + 1)
            push!(v, i)
            push!(w, diag[i])
        end
    end

    return sparse(u,v,w)

end

function buildSolver{Tv,Ti}(a::SparseMatrixCSC{Tv,Ti};
                            params::samplingParams{Tv,Ti}=defaultSamplingParams)

    # compute rho
    n = a.n;
    rho::Tv = params.sampConst * log(n) ^ 2 / params.eps ^ 2
    if params.verboseSS
	    println("rho = ", rho)
    end

    if params.verboseSS
	    tic()
    end
    #=
        Note: The n'th vertex will correspond to the vertex used to solve Laplacian systems
        with extra weight on the diagonal. Thus, we want to eliminate it last.
    =#

    tree = akpw(a);
    ord = reverse!(dfsOrder(tree, start = n));
    a = symPermuteCSC(a, ord)
    tree = symPermuteCSC(tree, ord)

    # Blow up the tree and compute the stretches
    a2 = copy(a)
    a = a + (params.beta - 1) * tree

    stretch = rho * compStretches(params.beta * tree, a)
    stretch.nzval = min(rho, stretch.nzval)

    if params.verboseSS
	    println("Initial number of multiedges = ", ceil(Int64,sum(stretch.nzval)), " . Nonzeros in a = ", nnz(a))
	    meanStretch = mean(stretch.nzval)
	    println("Average number of multiedges = ", mean(stretch.nzval))
	    maxStretch = maximum(stretch.nzval)
	    println("Maximum number of multiedges = ", maximum(stretch.nzval))
    end

  	if params.verboseSS
	    print("Time to build the tree and compute the stretch: ")
	    toc()
    end

    # Get u and d such that ut d u = -a (doesn't affect solver)
    Ut,d = samplingLDL(a, stretch, rho, params.startingSize, params.blockSize, params.verboseSS)
    U = Ut'

    if params.verboseSS
	    println("nnz in U matrix: ", length(U.data.nzval))
    end

    # Create the solver function
    f = function(b::Array{Float64,1})
        # center
        res = copy(b)
        res = res - sum(res) / n

        # forward solve
        res = Ut \ res

        # diag inverse
        for i in 1:(n - 1)
            res[i] = res[i] / d[i]
        end

        # backward solve
        res = U \ res

        # center
        res = res - sum(res) / n

        return res
    end

    # Create the error check function
    la = lap(a)
    g = function(b::Array{Float64,1})
        res = copy(b)
        res[n] = 0

        # diag sqrt inverse
        for i in 1:(n - 1)
            res[i] = res[i] * d[i]^(-1/2)
        end

        # backward solve #TODO?
        res = U \ res

        # apply lapl
        res = la * res

        # forward solve #TODO?
        res = Ut \ res

        # diag sqrt inverse
        for i in 1:(n - 1)
            res[i] = res[i] * d[i]^(-1/2)
        end

        # subtract identity, except we haven't zeroed out last coord
        res = res - b

        #zero out last coord
        res[n] = 0 #TODO?

        return res
    end
    gOp = SqLinOp(true,1.0,n,g)

    if params.returnCN
	    tic()
        cn = condNumber(lap(a2), U, d, tol=params.CNTol)
        print("computing the condition number takes: ")
	    cntime = toc()
    end

    if params.returnCN
        return f,gOp,U,d,ord,cn,cntime
    else
        return f,gOp,U,d,ord,(0.0,0.0),0.0
    end
end

# a is an adjacency matrix
function samplingLDL{Tv,Ti}(a::SparseMatrixCSC{Tv,Ti}, stretch::SparseMatrixCSC{Tv,Ti}, rho::Tv,
    startingSize::Ti, blockSize::Ti, verbose::Bool=false)

    n = a.n

    # some extra memory to be used later in the algorithm. this can be later pulled out of this function
    # into an external recipient, to be used on subsequent runs of the solver
    auxVal = zeros(Tv, n)                       # used to sum weights from multiedges
    auxMult = zeros(Tv, n)                      # used to count the number of multiedges

    wNeigh = zeros(Tv, n)
    multNeigh = zeros(Tv, n)
    indNeigh = zeros(Ti, n)

    ut = Array{Tuple{Tv,Ti},1}[[] for i in 1:n]  # the lower triangular u matrix part of u d u'
    d = zeros(Tv, n)                            # the d matrix part of u d u'

    # neigh[i] = the list of neighbors for vertex i with their corresponding weights
    # note neigh[i] only stores neighbors j such that j > i
    # neigh[i][1] is weight, [2] is number of multi-edges, [3] is neighboring vertex

    neigh = llsInit(a, startingSize = startingSize, blockSize = blockSize)

    # gather the info in a and put it into neigh and w
    for i in 1:length(a.colptr) - 1
        for j in a.colptr[i]:a.colptr[i + 1] - 1
            if a.rowval[j] > i
                llsAdd(neigh, i, (a.nzval[j], stretch.nzval[j], a.rowval[j]))
            end
        end
    end

    # Now, for every i, we will compute the i'th column in U
    for i in 1:(n-1)
        # We will get rid of duplicate edges
        # wSum - sum of weights of edges
        # multSum - sum of number of edges (including multiedges)
        # numPurged - the size in use of wNeigh, multNeigh and indNeigh
        # wNeigh - list of weights correspongind to each neighbors
        # multNeigh - list of number of multiedges to each neighbor
        # indNeigh - the indices of the neighboring vertices
        wSum, multSum, numPurged = llsPurge(neigh, i, auxVal, auxMult, wNeigh, multNeigh, indNeigh, rho = rho)

        # need to divide weights by the diagonal entry
        for j in 1:numPurged
            push!(ut[i], (-wNeigh[j] / wSum, indNeigh[j]))
        end
        push!(ut[i], (1, i)) #diag term

        d[i] = wSum

        multSum = ceil(Int64, multSum)
        wSamp = FastSampler(wNeigh[1:numPurged])
        multSamp = FastSampler(multNeigh[1:numPurged])

        jSamples = sampleMany(wSamp, multSum)
        kSamples = sampleMany(multSamp, multSum)

        # now propagate the clique to the neighbors of i
        for l in 1:multSum

            j = jSamples[l]
            k = kSamples[l]

            if j != k
                posj = indNeigh[j]
                posk = indNeigh[k]

                # swap so posj is smaller
                if posk < posj
                    j, k = k, j
                    posj, posk = posk, posj
                end

                wj = wNeigh[j]
                wk = wNeigh[k]

                sampScaling = wj * multNeigh[k] + wk * multNeigh[j]

                llsAdd(neigh, posj, (wj * wk / sampScaling, 1.0, posk))
            end
        end
    end

    # add the last diagonal term
    push!(ut[n], (1, n))
    d[n] = 0

    if verbose
	    println()
	    println("The total size of the linked list data structure should be at most ", ceil(Ti, sum(stretch) + rho * (n - 1)) + 20 * n)
	    println("The actual size is ", neigh.size * neigh.blockSize)
	    println()
    end

    return constructLowerTriangularMat(ut), d
end

# u is an array of arrays of tuples. to be useful, we need to convert it to a lowerTriangular matrix
function constructLowerTriangularMat{Tv,Ti}(u::Array{Array{Tuple{Tv,Ti},1},1})
    n = length(u)

    nnz = 0
    for i in 1:n
        nnz = nnz + length(u[i])
    end

    colptr = Array{Ti,1}(n + 1)
    rowval = Array{Ti,1}(nnz)
    nzval = Array{Tv,1}(nnz)

    colptr[1] = 1
    for i in 1:n
        colptr[i + 1] = colptr[i] + length(u[i])
    end
    index = copy(colptr)

    # We know that in u the values aren't necessarily ordered by row. So, we do a count sort-like algorith to keep linear time.
    helper = Array{Tuple{Tv,Ti},1}[[] for i in 1:n]
    for i in 1:n
        for j in 1:length(u[i])
            row = u[i][j][2]
            col = i
            val = u[i][j][1]
            push!(helper[row], (val, col))
        end
    end

    for i in 1:n
        for j in 1:length(helper[i])
            row = i
            col = helper[i][j][2]
            val = helper[i][j][1]

            rowval[index[col]] = row
            nzval[index[col]] = val
            index[col] = index[col] + 1
        end
    end

    return LowerTriangular(SparseMatrixCSC(n, n, colptr, rowval, nzval))
end

function checkError{Tv,Ti}(gOp::SqLinOp{Tv,Ti}; tol::Float64 = 0.0)
    return eigs(gOp;nev=1,which=:LM,tol=tol)[1][1]
end


text:
import Dates
using Compat

@compat abstract type RDate end
# Calendars
@compat abstract type CalendarManager end
@compat abstract type Calendar end
# Conventions
@compat abstract type HolidayRoundingConvention end
@compat abstract type InvalidDayConvention end
@compat abstract type MonthIncrementConvention end

"""
    NullCalendarManager() <: CalendarManager

The most primitive calendar manager, that will return an error for any request to
get a calendar. The default calendar manager that is available when applying rdates
using the + operator, without an explicit calendar manager.
"""
struct NullCalendarManager <: CalendarManager end

"""
    is_holiday(calendar::Calendar, date::Dates.Date)::Bool

Determine whether the date requested is a holiday for this calendar or not.
"""
is_holiday(x::Calendar, ::Dates.Date) = error("$(typeof(x)) does not support is_holiday")

"""
    calendar(calendarmgr::CalendarManager, names::Vector{String})::Calendar

Given a set of calendar names, request the calendar manager to retrieve the associated
calendar that supports the union of them.
"""
calendar(x::CalendarManager, names)::Calendar = error("$(typeof(x)) does not support calendar")

"""
    apply(rdate::RDate, date::Dates.Date, calendarmgr::CalendarManager)::Dates.Date

The application of an rdate to a specific date, given an explicit calendar manager.
"""
apply(rd::RDate, ::Dates.Date, ::CalendarManager)::Dates.Date = error("$(typeof(rd)) does not support date apply")

"""
    multiply_no_roll(rdate::RDate, count::Integer)::RDate

An internalised multiplication of an rdate which is generated without reapplication of a relative date multiple
times. This differs from multiply_roll when either calendars or month adjustments are involved.

For example "6m" + "6m" != "12m" for all dates, due to the fact that there are different days in each month
and invalid day conventions will kick in.
"""
multiply_no_roll(rd::RDate, count::Integer) = error("$(typeof(rd)) does not support no roll multiplication")

calendar_mgr = NullCalendarManager()
function with_cal_mgr(f::Function, cal_mgr::CalendarManager)
    global calendar_mgr
    current_cal_mgr = calendar_mgr
    calendar_mgr = cal_mgr
    try
        f()
    finally
        calendar_mgr = current_cal_mgr
    end
end

"""
    apply(rdate::RDate, date::Dates.Date)::Dates.Date

The application of an rdate to a specific date, without an explicit calendar manager.
This will use the globally specified calendar manager defined at the point of request.
"""
apply(rd::RDate, date::Dates.Date) = apply(rd, date, calendar_mgr)
Base.:+(rd::RDate, date::Dates.Date) = apply(rd, date)
Base.:+(date::Dates.Date, rd::RDate) = apply(rd, date)
Base.:-(date::Dates.Date, rd::RDate) = apply(-rd, date)
Base.:-(x::RDate)::RDate = error("$(typeof(x)) does not support negation")

"""
    apply(rounding::HolidayRoundingConvention, date::Dates.Date, calendar::Calendar)::Dates.Date

Apply the appropriate adjustment to the date if it falls on a holiday for the given calendar. Be aware
that there is no strict requirement that the resolved date will not be a holiday for the given calendar.
"""
apply(x::HolidayRoundingConvention, date::Dates.Date, calendar::Calendar) = error("$(typeof(x)) does not support apply")

"""
    apply(rounding::InvalidDayConvention, day::Integer, month::Integer, year::Integer)::Dates.Date

Given a day, month and year which do not generate a valid date, adjust them in some form to a valid date.
"""
adjust(x::InvalidDayConvention, day, month, year) = error("$(typeof(x)) does not support adjust")

"""
    apply(rounding::MonthIncrementConvention, from::Dates.Date, new_month::Integer, new_year::Integer, calendar_mgr::CalendarManager)::Tuple{Integer, Integer, Integer}

Given the initial day, month and year that we're moving from and a generated new month and new year, determine the new
day, month and year that should be used. Generally used to handle specific features around preservation of month ends
or days of week, if required.

Note that the generated (day, month, year) do not need to be able to produce a valid date. The invalid day convention
should be applied, if required, after this calculation.

May use a calendar manager if required as well
"""
adjust(x::MonthIncrementConvention, from::Dates.Date, new_month, new_year, calendar_mgr::CalendarManager) = error("$(typeof(x)) does not support adjust")