text:
using PCLCommon
using PCLIO
using PCLVisualization
using Libfreenect2
using Cxx

const w = 512
const h = 424
save_pcd = false

genfilename(ext=".pcd") =
    joinpath(dirname(@__FILE__), string(now(), "_", time_ns(), ext))

"""Get point cloud from undistored depth and reigstered color"""
function getPointCloudXYZRGB(registration::Registration, undistorted, registered)
    w = Libfreenect2.width(undistorted)
    h = Libfreenect2.height(undistorted)
    cloud = PointCloud{PointXYZRGB}(w, h)
    icxx"$(cloud.handle)->is_dense = false;"
    pointsptr = icxx"&$(cloud.handle)->points[0];"
    icxx"""
    for (size_t ri = 0; ri < $h; ++ri) {
        for (size_t ci = 0; ci < $w; ++ci) {
            auto p = $(pointsptr) + $w * ri + ci;
            $(registration.handle)->getPointXYZRGB($(undistorted.handle),
                $(registered.handle), ri, ci, p->x, p->y, p->z, p->rgb);
        }
    }
    """
    cloud
end

function ppcallback(event::cxxt"pcl::visualization::PointPickingEvent&",
        args::Ptr{Void})
    if Int(icxx"$(event).getPointIndex();") == -1
        return
    end
    p = PointXYZ()
    icxx"$(event).getPoint($p.x, $p.y, $p.z);"
    @show p
    return nothing::Void
end

f = Freenect2()
device = openDefaultDevice(f, OpenGLPacketPipeline())
listener = SyncMultiFrameListenerPtr()
setIrAndDepthFrameListener(device, listener)
setColorFrameListener(device, listener)

start(device)

# NOTE: must be called after start(device)
registration = Registration(getIrCameraParams(device),
    getColorCameraParams(device))
undistorted = FramePtr(w, h, 4, key=Libfreenect2.FRAME_DEPTH)
registered = FramePtr(w, h, 4, key=Libfreenect2.FRAME_COLOR)

info("Prepare PCL visualizer...")
global viewer = PCLVisualizer("pcl visualizer")
addCoordinateSystem(viewer, 1.2, 0,0,0)

registerPointPickingCallback(viewer, ppcallback)

global should_save = false
if !isdefined(:viewer_cb_defined)
cxx"""
void viewer_cb(const pcl::visualization::KeyboardEvent &event) {
    std::cout << "key event:" << event.getKeyCode() << std::endl;
    if (event.getKeyCode() == 's') {
        $:(global should_save = true; nothing);
    }
}
"""
const viewer_cb_defined = true
end

icxx"$(viewer.handle)->registerKeyboardCallback(viewer_cb);"

try
    while !wasStopped(viewer)
        frames = waitForNewFrame(listener)
        color = frames[FrameType.COLOR]
        ir = frames[FrameType.IR]
        depth = frames[FrameType.DEPTH]

        # Depth and color registration
        t = @elapsed begin
            apply(registration, color, depth, undistorted, registered)
            cloud = getPointCloudXYZRGB(registration, undistorted, registered)
        end
        println("Registration and getPointCloudXYZRGB: $(1/t) Hz")

        if should_save
            info("save pcd file...")
            savePCDFile(genfilename(), cloud; binary_mode=true)
            should_save = false
        end

        color_handler = PointCloudColorHandlerRGBField(cloud)
        if !updatePointCloud(viewer, cloud, color_handler, id="libfreenect2")
            addPointCloud(viewer, cloud, color_handler, id="libfreenect2")
        end
        try
            spinOnce(viewer, 1)
        catch e
            @show e
        end

        Libfreenect2.release(listener, frames)

        rand() > 0.95 && gc(false)
    end
finally
    stop(device)
    close(device)
    close(viewer)
    device=0;viewer=0;gc()
end


text:
struct GenericOpProblem <: PowerSimulationsOperationsProblem end

mutable struct OperationsProblem{M <: AbstractOperationsProblem}
    template::OperationsProblemTemplate
    sys::PSY.System
    psi_container::PSIContainer
end

"""
    OperationsProblem(::Type{M},
    template::OperationsProblemTemplate,
    sys::PSY.System,
    jump_model::Union{Nothing, JuMP.AbstractModel}=nothing;
    kwargs...) where {M<:AbstractOperationsProblem,
                      T<:PM.AbstractPowerFormulation}
This builds the optimization problem of type M with the specific system and template.
# Arguments
- `::Type{M} where M<:AbstractOperationsProblem`: The abstract operation model type
- `template::OperationsProblemTemplate`: The model reference made up of transmission, devices,
                                          branches, and services.
- `sys::PSY.System`: the system created using Power Systems
- `jump_model::Union{Nothing, JuMP.AbstractModel}`: Enables passing a custom JuMP model. Use with care
# Output
- `op_problem::OperationsProblem`: The operation model containing the model type, built JuMP model, Power
Systems system.
# Example
```julia
template = OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)
OpModel = OperationsProblem(TestOpProblem, template, system)
```
# Accepted Key Words
- `horizon::Int`: Manually specify the length of the forecast Horizon
- `initial_time::Dates.DateTime`: Initial Time for the model solve
- `use_forecast_data::Bool` : If true uses the data in the system forecasts. If false uses the data for current operating point in the system.
- `PTDF::PTDF`: Passes the PTDF matrix into the optimization model for StandardPTDFModel networks.
- `optimizer::JuMP.MOI.OptimizerWithAttributes`: The optimizer that will be used in the optimization model.
- `use_parameters::Bool`: True will substitute will implement formulations using ParameterJuMP parameters. Defatul is false.
- `warm_start::Bool` True will use the current operation point in the system to initialize variable values. False initializes all variables to zero. Default is true
- `balance_slack_variables::Bool` True will add slacks to the system balance constraints
- `services_slack_variables::Bool` True will add slacks to the services requirement constraints
"""
function OperationsProblem(
    ::Type{M},
    template::OperationsProblemTemplate,
    sys::PSY.System,
    jump_model::Union{Nothing, JuMP.AbstractModel} = nothing;
    kwargs...,
) where {M <: AbstractOperationsProblem}
    return OperationsProblem{M}(template, sys, jump_model; kwargs...)
end

function OperationsProblem{M}(
    template::OperationsProblemTemplate,
    sys::PSY.System,
    jump_model::Union{Nothing, JuMP.AbstractModel} = nothing;
    kwargs...,
) where {M <: AbstractOperationsProblem}
    check_kwargs(kwargs, OPERATIONS_ACCEPTED_KWARGS, "OperationsProblem")
    settings = PSISettings(sys; kwargs...)
    return OperationsProblem{M}(template, sys, jump_model, settings)
end

# The psi_container_init is called at the build! call in this constructor. This is meant to
# build an operation from a template.
function OperationsProblem{M}(
    template::OperationsProblemTemplate,
    sys::PSY.System,
    jump_model::Union{Nothing, JuMP.AbstractModel},
    settings::PSISettings,
) where {M <: AbstractOperationsProblem}
    op_problem =
        OperationsProblem{M}(template, sys, PSIContainer(sys, settings, jump_model))
    build!(op_problem)
    return op_problem
end

"""
    OperationsProblem(op_problem::Type{M},
                    ::Type{T},
                    sys::PSY.System,
                    jump_model::Union{Nothing, JuMP.AbstractModel}=nothing;
                    kwargs...) where {M<:AbstractOperationsProblem,
                                      T<:PM.AbstractPowerFormulation}
Return an unbuilt operation problem of type M with the specific system and network model T.
    This constructor doesn't build any device model; it is meant to built device models individually using [`construct_device!`](@ref)
# Arguments
- `::Type{M} where M<:AbstractOperationsProblem`: The abstract operation model type
- `::Type{T} where T<:AbstractPowerModel`: The abstract network formulation
- `sys::PSY.System`: the system created using Power Systems
- `jump_model::Union{Nothing, JuMP.AbstractModel}`: Enables passing a custom JuMP model. Use with care
# Output
- `op_problem::OperationsProblem`: The operation model containing the model type, unbuilt JuMP model, Power
Systems system.
# Example
```julia
OpModel = OperationsProblem(MyCustomOpProblem, DCPPowerModel, system)
model = DeviceModel(ThermalStandard, ThermalStandardUnitCommitment)
construct_device!(op_problem, :Thermal, model)
```
# Accepted Key Words
- `horizon::Int`: Manually specify the length of the forecast Horizon
- `initial_time::Dates.DateTime`: Initial Time for the model solve
- `use_forecast_data::Bool` : If true uses the data in the system forecasts. If false uses the data for current operating point in the system.
- `PTDF::PTDF`: Passes the PTDF matrix into the optimization model for StandardPTDFModel networks.
- `optimizer::JuMP.MOI.OptimizerWithAttributes`: The optimizer that will be used in the optimization model.
- `use_parameters::Bool`: True will substitute will implement formulations using ParameterJuMP parameters. Defatul is false.
- `warm_start::Bool` True will use the current operation point in the system to initialize variable values. False initializes all variables to zero. Default is true
- `balance_slack_variables::Bool` True will add slacks to the system balance constraints
- `services_slack_variables::Bool` True will add slacks to the services requirement constraints
"""
function OperationsProblem(
    ::Type{M},
    ::Type{T},
    sys::PSY.System,
    jump_model::Union{Nothing, JuMP.AbstractModel} = nothing;
    kwargs...,
) where {M <: AbstractOperationsProblem, T <: PM.AbstractPowerModel}
    return OperationsProblem{M}(T, sys, jump_model; kwargs...)
end

"""
    OperationsProblem(::Type{T},
                    sys::PSY.System,
                    jump_model::Union{Nothing, JuMP.AbstractModel}=nothing;
                    kwargs...) where {M<:AbstractOperationsProblem,
                                      T<:PM.AbstractPowerFormulation}
Return an unbuilt operation problem of type GenericOpProblem with the specific system and network model T.
    This constructor doesn't build any device model; it is meant to built device models individually using [`construct_device!`](@ref)
# Arguments
- `::Type{T} where T<:AbstractPowerModel`: The abstract network formulation
- `sys::PSY.System`: the system created using Power Systems
- `jump_model::Union{Nothing, JuMP.AbstractModel}`: Enables passing a custom JuMP model. Use with care
# Output
- `op_problem::OperationsProblem`: The operation model containing the model type, unbuilt JuMP model, Power
Systems system.
# Example
```julia
OpModel = OperationsProblem(DCPPowerModel, system)
model = DeviceModel(ThermalStandard, ThermalStandardUnitCommitment)
construct_device!(op_problem, :Thermal, model)
```
# Accepted Key Words
- `horizon::Int`: Manually specify the length of the forecast Horizon
- `initial_time::Dates.DateTime`: Initial Time for the model solve
- `use_forecast_data::Bool` : If true uses the data in the system forecasts. If false uses the data for current operating point in the system.
- `PTDF::PTDF`: Passes the PTDF matrix into the optimization model for StandardPTDFModel networks.
- `optimizer::JuMP.MOI.OptimizerWithAttributes`: The optimizer that will be used in the optimization model.
- `use_parameters::Bool`: True will substitute will implement formulations using ParameterJuMP parameters. Defatul is false.
- `warm_start::Bool` True will use the current operation point in the system to initialize variable values. False initializes all variables to zero. Default is true
- `balance_slack_variables::Bool` True will add slacks to the system balance constraints
- `services_slack_variables::Bool` True will add slacks to the services requirement constraints
"""
function OperationsProblem(
    ::Type{T},
    sys::PSY.System,
    jump_model::Union{Nothing, JuMP.AbstractModel} = nothing;
    kwargs...,
) where {T <: PM.AbstractPowerModel}
    return OperationsProblem{GenericOpProblem}(T, sys, jump_model; kwargs...)
end

# This constructor calls PSI container including the the PowerModels type in order to initialize
# the container and it is meant to construct operations problems using construct_device! function
# calls
function OperationsProblem{M}(
    ::Type{T},
    sys::PSY.System,
    jump_model::Union{Nothing, JuMP.AbstractModel} = nothing;
    kwargs...,
) where {M <: AbstractOperationsProblem, T <: PM.AbstractPowerModel}
    check_kwargs(kwargs, OPERATIONS_ACCEPTED_KWARGS, "OperationsProblem")
    settings = PSISettings(sys; kwargs...)
    return OperationsProblem{M}(
        OperationsProblemTemplate(T),
        sys,
        PSIContainer(T, sys, settings, jump_model),
    )

end

"""
    OperationsProblem(filename::AbstractString)

Construct an OperationsProblem from a serialized file.

# Arguments
- `filename::AbstractString`: path to serialized file
- `jump_model::Union{Nothing, JuMP.AbstractModel}` = nothing: The JuMP model does not get
   serialized. Callers should pass whatever they passed to the original problem.
- `optimizer::Union{Nothing,JuMP.MOI.OptimizerWithAttributes}` = nothing: The optimizer does
   not get serialized. Callers should pass whatever they passed to the original problem.
"""
function OperationsProblem(
    filename::AbstractString;
    jump_model::Union{Nothing, JuMP.AbstractModel} = nothing,
    optimizer::Union{Nothing, JuMP.MOI.OptimizerWithAttributes} = nothing,
)
    return deserialize(
        OperationsProblem,
        filename;
        jump_model = jump_model,
        optimizer = optimizer,
    )
end

get_transmission_ref(op_problem::OperationsProblem) = op_problem.template.transmission
get_devices_ref(op_problem::OperationsProblem) = op_problem.template.devices
get_branches_ref(op_problem::OperationsProblem) = op_problem.template.branches
get_services_ref(op_problem::OperationsProblem) = op_problem.template.services
get_system(op_problem::OperationsProblem) = op_problem.sys
get_psi_container(op_problem::OperationsProblem) = op_problem.psi_container
get_base_power(op_problem::OperationsProblem) = PSY.get_base_power(op_problem.sys)
get_jump_model(op_problem::OperationsProblem) = get_jump_model(op_problem.psi_container)

function reset!(op_problem::OperationsProblem)
    op_problem.psi_container =
        PSIContainer(op_problem.sys, op_problem.psi_container.settings, nothing)
    return
end

function set_transmission_model!(
    op_problem::OperationsProblem{M},
    transmission::Type{T},
) where {T <: PM.AbstractPowerModel, M <: AbstractOperationsProblem}
    op_problem.template.transmission = transmission
    reset!(op_problem)
    build!(op_problem)
    return
end

function set_devices_template!(
    op_problem::OperationsProblem{M},
    devices::Dict{Symbol, DeviceModel},
) where {M <: AbstractOperationsProblem}
    op_problem.template.devices = devices
    reset!(op_problem)
    build!(op_problem)
    return
end

function set_branches_template!(
    op_problem::OperationsProblem{M},
    branches::Dict{Symbol, DeviceModel},
) where {M <: AbstractOperationsProblem}
    op_problem.template.branches = branches
    reset!(op_problem)
    build!(op_problem)
    return
end

function set_services_template!(
    op_problem::OperationsProblem{M},
    services::Dict{Symbol, <:ServiceModel},
) where {M <: AbstractOperationsProblem}
    op_problem.template.services = services
    reset!(op_problem)
    build!(op_problem)
    return
end

function set_device_model!(
    op_problem::OperationsProblem{M},
    name::Symbol,
    device::DeviceModel{<:PSY.StaticInjection, <:AbstractDeviceFormulation},
) where {M <: AbstractOperationsProblem}
    if haskey(op_problem.template.devices, name)
        op_problem.template.devices[name] = device
        reset!(op_problem)
        build!(op_problem)
    else
        throw(IS.ConflictingInputsError("Device Model with name $(name) doesn't exist in the model"))
    end
    return
end

function set_branch_model!(
    op_problem::OperationsProblem{M},
    name::Symbol,
    branch::DeviceModel{<:PSY.Branch, <:AbstractDeviceFormulation},
) where {M <: AbstractOperationsProblem}
    if haskey(op_problem.template.branches, name)
        op_problem.template.branches[name] = branch
        reset!(op_problem)
        build!(op_problem)
    else
        throw(IS.ConflictingInputsError("Branch Model with name $(name) doesn't exist in the model"))
    end
    return
end

function set_services_model!(
    op_problem::OperationsProblem{M},
    name::Symbol,
    service::ServiceModel,
) where {M <: AbstractOperationsProblem}
    if haskey(op_problem.template.services, name)
        op_problem.template.services[name] = service
        reset!(op_problem)
        build!(op_problem)
    else
        throw(IS.ConflictingInputsError("Branch Model with name $(name) doesn't exist in the model"))
    end
    return
end

function set_model!(
    ::Type{D},
    op_problem::OperationsProblem,
    name::Symbol,
    device_model::DeviceModel,
) where {D <: PSY.Branch}
    op_problem.template.branches[name] = device_model
end

function set_model!(
    ::Type{D},
    op_problem::OperationsProblem,
    name::Symbol,
    device_model::DeviceModel,
) where {D <: PSY.StaticInjection}
    set_model!(op_problem.template, name, device_model)
end

function set_model!(
    ::Type{D},
    op_problem::OperationsProblem,
    name::Symbol,
    device_model::DeviceModel,
) where {D <: PSY.Service}
    op_problem.template.services[name] = device_model
end

function construct_device!(
    op_problem::OperationsProblem,
    name::Symbol,
    device_model::DeviceModel,
)
    set_model!(device_model.device_type, op_problem, name, device_model)

    construct_device!(
        op_problem.psi_container,
        get_system(op_problem),
        device_model,
        get_transmission_ref(op_problem),
    )
    JuMP.@objective(
        op_problem.psi_container.JuMPmodel,
        MOI.MIN_SENSE,
        op_problem.psi_container.cost_function
    )
    return
end

function construct_network!(op_problem::OperationsProblem)
    construct_network!(op_problem, op_problem.template.transmission)
    return
end

function construct_network!(
    op_problem::OperationsProblem,
    system_formulation::Type{T},
) where {T <: PM.AbstractPowerModel}
    construct_network!(op_problem.psi_container, get_system(op_problem), T)
    return
end

function get_initial_conditions(
    op_problem::OperationsProblem,
    ic::InitialConditionType,
    device::PSY.Device,
)
    psi_container = op_problem.psi_container
    key = ICKey(ic, device)

    return get_initial_conditions(psi_container, key)
end

function build!(op_problem::OperationsProblem{M}) where {M <: AbstractOperationsProblem}
    sys = get_system(op_problem)
    _build!(op_problem.psi_container, op_problem.template, sys)
    return
end

function check_problem_size(psi_container::PSIContainer)
    vars = JuMP.num_variables(psi_container.JuMPmodel)
    cons = 0
    for (exp, c_type) in JuMP.list_of_constraint_types(psi_container.JuMPmodel)
        cons += JuMP.num_constraints(psi_container.JuMPmodel, exp, c_type)
    end
    return "The current total number of variables is $(vars) and total number of constraints is $(cons)"
end

function _build!(
    psi_container::PSIContainer,
    template::OperationsProblemTemplate,
    sys::PSY.System,
)
    transmission = template.transmission
    # Order is required
    # The container is initialized here because this build! call for psi_container takes the
    # information from the template with cached PSISettings. It allows having the same build! call for operations problems
    # specified with template and simulation stage.
    psi_container_init!(psi_container, transmission, sys)
    construct_services!(psi_container, sys, template.services, template.devices)
    for device_model in values(template.devices)
        @debug "Building $(device_model.device_type) with $(device_model.formulation) formulation"
        construct_device!(psi_container, sys, device_model, transmission)
        @debug check_problem_size(psi_container)
    end
    @debug "Building $(transmission) network formulation"
    construct_network!(psi_container, sys, transmission)
    @debug check_problem_size(psi_container)

    for branch_model in values(template.branches)
        @debug "Building $(branch_model.device_type) with $(branch_model.formulation) formulation"
        construct_device!(psi_container, sys, branch_model, transmission)
        @debug check_problem_size(psi_container)
    end

    @debug "Building Objective"
    JuMP.@objective(psi_container.JuMPmodel, MOI.MIN_SENSE, psi_container.cost_function)
    @debug "Total operation count $(psi_container.JuMPmodel.operator_counter)"
    return
end

function get_variables_value(op_m::OperationsProblem)
    results_dict = Dict{Symbol, DataFrames.DataFrame}()
    for (k, v) in get_variables(op_m.psi_container)
        results_dict[k] = axis_array_to_dataframe(v)
    end
    return results_dict
end

function get_dual_values(op_m::OperationsProblem)
    return get_dual_values(op_m.psi_container)
end

"""
    solve!(op_problem::OperationsProblem; kwargs...)
This solves the operational model for a single instance and
outputs results of type OperationsProblemResult
# Arguments
- `op_problem::OperationModel = op_problem`: operation model
# Examples
```julia
results = solve!(OpModel)
```
# Accepted Key Words
- `save_path::String`: If a file path is provided the results
automatically get written to feather files
- `optimizer::MOI.OptimizerWithAttributes`: The optimizer that is used to solve the model
"""
function solve!(
    op_problem::OperationsProblem{T};
    kwargs...,
) where {T <: AbstractOperationsProblem}
    check_kwargs(kwargs, OPERATIONS_SOLVE_KWARGS, "Solve")
    timed_log = Dict{Symbol, Any}()
    save_path = get(kwargs, :save_path, nothing)

    if op_problem.psi_container.JuMPmodel.moi_backend.state == MOIU.NO_OPTIMIZER
        if !(:optimizer in keys(kwargs))
            error("No Optimizer has been defined, can't solve the operational problem")
        end
        JuMP.set_optimizer(op_problem.psi_container.JuMPmodel, kwargs[:optimizer])
        _,
        timed_log[:timed_solve_time],
        timed_log[:solve_bytes_alloc],
        timed_log[:sec_in_gc] = @timed JuMP.optimize!(op_problem.psi_container.JuMPmodel)
    else
        _,
        timed_log[:timed_solve_time],
        timed_log[:solve_bytes_alloc],
        timed_log[:sec_in_gc] = @timed JuMP.optimize!(op_problem.psi_container.JuMPmodel)
    end
    model_status = JuMP.primal_status(op_problem.psi_container.JuMPmodel)
    if model_status != MOI.FEASIBLE_POINT::MOI.ResultStatusCode
        error("The Operational Problem $(T) status is $(model_status)")
    end
    vars_result = get_variables_value(op_problem)
    param_values = get_parameters_value(get_psi_container(op_problem))
    optimizer_log = get_optimizer_log(op_problem)
    time_stamp = get_time_stamps(op_problem)
    time_stamp = shorten_time_stamp(time_stamp)
    base_power = PSY.get_base_power(op_problem.sys)
    dual_result = get_dual_values(op_problem)
    obj_value = Dict(
        :OBJECTIVE_FUNCTION => JuMP.objective_value(op_problem.psi_container.JuMPmodel),
    )
    base_power = get_base_power(op_problem)
    merge!(optimizer_log, timed_log)

    results = OperationsProblemResults(
        base_power,
        vars_result,
        obj_value,
        optimizer_log,
        time_stamp,
        dual_result,
        param_values,
    )

    !isnothing(save_path) && write_results(results, save_path)

    return results
end

# Function to create a dictionary for the optimizer log of the simulation
function get_optimizer_log(op_m::OperationsProblem)
    psi_container = op_m.psi_container
    optimizer_log = Dict{Symbol, Any}()
    optimizer_log[:obj_value] = JuMP.objective_value(psi_container.JuMPmodel)
    optimizer_log[:termination_status] = JuMP.termination_status(psi_container.JuMPmodel)
    optimizer_log[:primal_status] = JuMP.primal_status(psi_container.JuMPmodel)
    optimizer_log[:dual_status] = JuMP.dual_status(psi_container.JuMPmodel)
    optimizer_log[:solver] = JuMP.solver_name(psi_container.JuMPmodel)

    try
        optimizer_log[:solve_time] = MOI.get(psi_container.JuMPmodel, MOI.SolveTime())
    catch
        @warn("SolveTime() property not supported by $(optimizer_log[:solver])")
        optimizer_log[:solve_time] = "Not Supported by $(optimizer_log[:solver])"
    end
    return optimizer_log
end

# Function to create a dictionary for the time series of the simulation
function get_time_stamps(op_problem::OperationsProblem)
    initial_time = model_initial_time(get_psi_container(op_problem))
    interval = PSY.get_forecasts_resolution(op_problem.sys)
    horizon = get_horizon(get_settings(get_psi_container(op_problem)))
    range_time = collect(initial_time:interval:(initial_time + interval .* horizon))
    time_stamp = DataFrames.DataFrame(Range = range_time[:, 1])

    return time_stamp
end

function write_data(psi_container::PSIContainer, save_path::AbstractString; kwargs...)
    file_type = get(kwargs, :file_type, Feather)
    if file_type == Feather || file_type == CSV
        for (k, v) in get_variables(psi_container)
            file_path = joinpath(save_path, "$(k).$(lowercase("$file_type"))")
            variable = axis_array_to_dataframe(v)
            if isempty(variable)
                @debug "$(k) is empty, not writing $file_path"
            else
                file_type.write(file_path, variable)
            end
        end
    end
    return
end

function write_data(op_problem::OperationsProblem, save_path::String; kwargs...)
    write_data(op_problem.psi_container, save_path; kwargs...)
    return
end

""" Exports the OpModel JuMP object in MathOptFormat"""
function export_op_model(op_problem::OperationsProblem, save_path::String)
    _write_psi_container(op_problem.psi_container, save_path)
    return
end

################ Functions to debug optimization models#####################################
""" "Each Tuple corresponds to (con_name, internal_index, moi_index)"""
function get_all_constraint_index(op_problem::OperationsProblem)
    con_index = Vector{Tuple{Symbol, Int, Int}}()
    for (key, value) in op_problem.psi_container.constraints
        for (idx, constraint) in enumerate(value)
            moi_index = JuMP.optimizer_index(constraint)
            push!(con_index, (key, idx, moi_index.value))
        end
    end
    return con_index
end

""" "Each Tuple corresponds to (con_name, internal_index, moi_index)"""
function get_all_var_index(op_problem::OperationsProblem)
    var_index = Vector{Tuple{Symbol, Int, Int}}()
    for (key, value) in op_problem.psi_container.variables
        for (idx, variable) in enumerate(value)
            moi_index = JuMP.optimizer_index(variable)
            push!(var_index, (key, idx, moi_index.value))
        end
    end
    return var_index
end

function get_con_index(op_problem::OperationsProblem, index::Int)
    for i in get_all_constraint_index(op_problem::OperationsProblem)
        if i[3] == index
            return op_problem.psi_container.constraints[i[1]].data[i[2]]
        end
    end

    @info "Index not found"
    return
end

function get_var_index(op_problem::OperationsProblem, index::Int)
    for i in get_all_var_index(op_problem::OperationsProblem)
        if i[3] == index
            return op_problem.psi_container.variables[i[1]].data[i[2]]
        end
    end
    @info "Index not found"
    return
end

function serialize(op_problem::OperationsProblem, filename::AbstractString)
    # A PowerSystem cannot be serialized in this format because of how it stores
    # time series data. Use its specialized serialization method instead.
    sys_filename = "$(basename(filename))-system-$(IS.get_uuid(op_problem.sys)).json"
    PSY.to_json(op_problem.sys, sys_filename)
    obj = OperationsProblemSerializationWrapper(
        op_problem.template,
        sys_filename,
        op_problem.psi_container.settings_copy,
        typeof(op_problem),
    )
    Serialization.serialize(filename, obj)
    @info "Serialized OperationsProblem to" filename
end

function deserialize(::Type{OperationsProblem}, filename::AbstractString; kwargs...)
    obj = Serialization.deserialize(filename)
    if !(obj isa OperationsProblemSerializationWrapper)
        throw(IS.DataFormatError("deserialized object has incorrect type $(typeof(obj))"))
    end

    if !ispath(obj.sys)
        throw(IS.DataFormatError("PowerSystems.System file $(obj.sys) does not exist"))
    end
    sys = PSY.System(obj.sys)

    return obj.op_problem_type(
        obj.template,
        sys,
        kwargs[:jump_model],
        restore_from_copy(obj.settings; optimizer = kwargs[:optimizer]),
    )
end

struct OperationsProblemSerializationWrapper
    template::OperationsProblemTemplate
    sys::String
    settings::PSISettings
    op_problem_type::DataType
end

ext: jl


text:
function mass_matrix_avr_entries!(
    mass_matrix,
    avr::T,
    global_index::Base.ImmutableDict{Symbol, Int64},
) where {T <: PSY.AVR}
    @debug "Using default mass matrix entries $T"
    return
end

function mass_matrix_avr_entries!(
    mass_matrix,
    avr::PSY.SEXS,
    global_index::Base.ImmutableDict{Symbol, Int64},
)
    mass_matrix[global_index[:Vf], global_index[:Vf]] = PSY.get_Te(avr)
    mass_matrix[global_index[:Vr], global_index[:Vr]] = PSY.get_Tb(avr)
    return
end

function mass_matrix_avr_entries!(
    mass_matrix,
    avr::PSY.EXST1,
    global_index::Base.ImmutableDict{Symbol, Int64},
)
    mass_matrix[global_index[:Vm], global_index[:Vm]] = PSY.get_Tr(avr)
    mass_matrix[global_index[:Vrll], global_index[:Vrll]] = PSY.get_Tb(avr)
    mass_matrix[global_index[:Vr], global_index[:Vr]] = PSY.get_Ta(avr)
    return
end

function mass_matrix_avr_entries!(
    mass_matrix,
    avr::PSY.EXAC1,
    global_index::Base.ImmutableDict{Symbol, Int64},
)
    mass_matrix[global_index[:Vm], global_index[:Vm]] = PSY.get_Tr(avr)
    mass_matrix[global_index[:Vr1], global_index[:Vr1]] = PSY.get_Tb(avr)
    mass_matrix[global_index[:Vr2], global_index[:Vr2]] = PSY.get_Ta(avr)
    return
end

function mdl_avr_ode!(
    ::AbstractArray{<:ACCEPTED_REAL_TYPES},
    ::AbstractArray{<:ACCEPTED_REAL_TYPES},
    inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{M, S, PSY.AVRFixed, TG, P}},
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}

    #Update Vf voltage on inner vars. In AVRFixed, Vf = V_ref
    inner_vars[Vf_var] = get_V_ref(dynamic_device)

    return
end

function mdl_avr_ode!(
    device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
    output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},
    inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{M, S, PSY.AVRSimple, TG, P}},
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}

    #Obtain references
    V_ref = get_V_ref(dynamic_device)

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(dynamic_device, PSY.AVRSimple)

    #Define inner states for component
    internal_states = @view device_states[local_ix]
    Vf = internal_states[1]

    #Define external states for device
    V_th = sqrt(inner_vars[VR_gen_var]^2 + inner_vars[VI_gen_var]^2)

    #Get Parameters
    Kv = PSY.get_Kv(PSY.get_avr(dynamic_device))

    #Compute ODEs
    output_ode[local_ix[1]] = Kv * (V_ref - V_th)

    #Update inner_vars
    inner_vars[Vf_var] = Vf

    return
end

function mdl_avr_ode!(
    device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
    output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},
    inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{M, S, PSY.AVRTypeI, TG, P}},
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}

    #Obtain references
    V0_ref = get_V_ref(dynamic_device)

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(dynamic_device, PSY.AVRTypeI)

    #Define inner states for component
    internal_states = @view device_states[local_ix]
    Vf = internal_states[1]
    Vr1 = internal_states[2]
    Vr2 = internal_states[3]
    Vm = internal_states[4]

    #Define external states for device
    V_th = sqrt(inner_vars[VR_gen_var]^2 + inner_vars[VI_gen_var]^2)
    Vs = inner_vars[V_pss_var]

    #Get parameters
    avr = PSY.get_avr(dynamic_device)
    Ka = PSY.get_Ka(avr)
    Ke = PSY.get_Ke(avr)
    Kf = PSY.get_Kf(avr)
    Ta = PSY.get_Ta(avr)
    Te = PSY.get_Te(avr)
    Tf = PSY.get_Tf(avr)
    Tr = PSY.get_Tr(avr)
    Ae = PSY.get_Ae(avr)
    Be = PSY.get_Be(avr)

    #Compute auxiliary parameters
    Se_Vf = Ae * exp(Be * abs(Vf)) #16.13
    V_ref = V0_ref + Vs

    # Compute block derivatives
    _, dVm_dt = low_pass(V_th, Vm, 1.0, Tr)
    y_hp, dVr2_dt = high_pass(Vf, Vr2, Kf, Tf)
    _, dVr1_dt = low_pass(V_ref - Vm - y_hp, Vr1, Ka, Ta)
    _, dVf_dt = low_pass_modified(Vr1, Vf, 1.0, Ke + Se_Vf, Te)

    #Compute 4 States AVR ODE:
    output_ode[local_ix[1]] = dVf_dt #16.12c
    output_ode[local_ix[2]] = dVr1_dt #16.12a
    output_ode[local_ix[3]] = dVr2_dt #16.12b
    output_ode[local_ix[4]] = dVm_dt #16.11

    #Update inner_vars
    inner_vars[Vf_var] = Vf

    return
end

function mdl_avr_ode!(
    device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
    output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},
    inner_vars,
    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{M, S, PSY.AVRTypeII, TG, P}},
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}

    #Obtain references
    V0_ref = get_V_ref(dynamic_device)

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(dynamic_device, PSY.AVRTypeII)

    #Define inner states for component
    internal_states = @view device_states[local_ix]
    Vf = internal_states[1]
    Vr1 = internal_states[2]
    Vr2 = internal_states[3]
    Vm = internal_states[4]

    #Define external states for device
    V_th = sqrt(inner_vars[VR_gen_var]^2 + inner_vars[VI_gen_var]^2)
    Vs = inner_vars[V_pss_var]

    #Get parameters
    avr = PSY.get_avr(dynamic_device)
    K0 = PSY.get_K0(avr)
    T1 = PSY.get_T1(avr)
    T2 = PSY.get_T2(avr)
    T3 = PSY.get_T3(avr)
    T4 = PSY.get_T4(avr)
    Te = PSY.get_Te(avr)
    Tr = PSY.get_Tr(avr)
    Ae = PSY.get_Ae(avr)
    Be = PSY.get_Be(avr)
    Va_min, Va_max = PSY.get_Va_lim(avr)

    #Compute auxiliary parameters
    Se_Vf = Ae * exp(Be * abs(Vf)) #16.13
    V_ref = V0_ref + Vs

    # Compute block derivatives
    _, dVm_dt = low_pass(V_th, Vm, 1.0, Tr)
    y_ll1, dVr1_dt = lead_lag(V_ref - Vm, Vr1, K0, T2, T1)
    y_ll2, dVr2_dt = lead_lag(y_ll1, K0 * Vr2, 1.0, K0 * T4, K0 * T3)
    Vr = clamp(y_ll2, Va_min, Va_max)
    _, dVf_dt = low_pass_modified(Vr, Vf, 1.0, 1.0 + Se_Vf, Te)

    #Compute 4 States AVR ODE:
    output_ode[local_ix[1]] = dVf_dt #16.18
    output_ode[local_ix[2]] = dVr1_dt #16.14
    output_ode[local_ix[3]] = dVr2_dt
    output_ode[local_ix[4]] = dVm_dt #16.11

    #Update inner_vars
    inner_vars[Vf_var] = Vf

    return
end

function mdl_avr_ode!(
    device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
    output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},
    inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{M, S, PSY.ESAC1A, TG, P}},
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}

    #Obtain references
    V_ref = get_V_ref(dynamic_device)

    #Obtain avr
    avr = PSY.get_avr(dynamic_device)

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(dynamic_device, typeof(avr))

    #Define inner states for component
    internal_states = @view device_states[local_ix]
    Vm = internal_states[1]
    Vr1 = internal_states[2]
    Vr2 = internal_states[3]
    Ve = internal_states[4]
    Vr3 = internal_states[5]

    #Define external states for device
    V_th = sqrt(inner_vars[VR_gen_var]^2 + inner_vars[VI_gen_var]^2)
    Vs = inner_vars[V_pss_var]
    Xad_Ifd = inner_vars[Xad_Ifd_var]

    #Get parameters
    Tr = PSY.get_Tr(avr)
    Ta = PSY.get_Ta(avr)
    Tb = PSY.get_Tb(avr)
    Tc = PSY.get_Tc(avr)
    inv_Tr = Tr < eps() ? 1.0 : 1.0 / Tr
    Ka = PSY.get_Ka(avr)
    Va_min, Va_max = PSY.get_Va_lim(avr) #Not used without UEL or OEL
    Te = PSY.get_Te(avr) # Te > 0
    Kf = PSY.get_Kf(avr)
    Tf = PSY.get_Tf(avr) # Te > 0
    Kc = PSY.get_Kc(avr)
    Kd = PSY.get_Kd(avr)
    Ke = PSY.get_Ke(avr)
    Vr_min, Vr_max = PSY.get_Vr_lim(avr)
    #Obtain saturation
    Se = saturation_function(avr, Ve)

    #Compute auxiliary parameters
    I_N = Kc * Xad_Ifd / Ve
    V_FE = Kd * Xad_Ifd + Ke * Ve + Se * Ve
    Vf = Ve * rectifier_function(I_N)

    # Compute blocks
    _, dVm_dt = low_pass(V_th, Vm, 1.0, 1.0 / inv_Tr)
    V_F, dVr3_dt = high_pass(V_FE, Vr3, Kf, Tf)
    V_in = V_ref + Vs - Vm - V_F
    y_ll, dVr1_dt = lead_lag(V_in, Vr1, 1.0, Tc, Tb)
    _, dVr2_dt = low_pass_nonwindup(y_ll, Vr2, Ka, Ta, Va_min, Va_max)

    #Set clamping for Vr2.
    V_R = clamp(Vr2, Vr_min, Vr_max)

    #Compute 4 States AVR ODE:
    output_ode[local_ix[1]] = dVm_dt
    output_ode[local_ix[2]] = dVr1_dt #dVr1/dt
    output_ode[local_ix[3]] = dVr2_dt
    output_ode[local_ix[4]] = (1.0 / Te) * (V_R - V_FE) #dVe/dt
    output_ode[local_ix[5]] = dVr3_dt

    #Update inner_vars
    inner_vars[Vf_var] = Vf

    return
end

function mdl_avr_ode!(
    device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
    output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},
    inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{M, S, PSY.SEXS, TG, P}},
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}

    #Obtain references
    V0_ref = get_V_ref(dynamic_device)

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(dynamic_device, PSY.SEXS)

    #Define inner states for component
    internal_states = @view device_states[local_ix]
    Vf = internal_states[1]
    Vr = internal_states[2]

    #Define external states for device
    V_th = sqrt(inner_vars[VR_gen_var]^2 + inner_vars[VI_gen_var]^2)
    Vs = inner_vars[V_pss_var]

    #Get parameters
    avr = PSY.get_avr(dynamic_device)
    Ta_Tb = PSY.get_Ta_Tb(avr)
    Tb = PSY.get_Tb(avr)
    Ta = Tb * Ta_Tb
    Te = PSY.get_Te(avr)
    K = PSY.get_K(avr)
    V_min, V_max = PSY.get_V_lim(avr)

    #Compute auxiliary parameters
    V_in = V0_ref + Vs - V_th
    V_LL, dVr_dt = lead_lag_mass_matrix(V_in, Vr, 1.0, Ta, Tb)
    _, dVf_dt = low_pass_nonwindup_mass_matrix(V_LL, Vf, K, Te, V_min, V_max)

    #Compute 2 States AVR ODE:
    output_ode[local_ix[1]] = dVf_dt
    output_ode[local_ix[2]] = dVr_dt

    #Update inner_vars
    inner_vars[Vf_var] = Vf

    return
end

function mdl_avr_ode!(
    device_states::AbstractArray,
    output_ode::AbstractArray,
    inner_vars::AbstractArray,
    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{M, S, PSY.EXST1, TG, P}},
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}

    #Obtain references
    V0_ref = get_V_ref(dynamic_device)

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(dynamic_device, PSY.EXST1)

    #Define inner states for component
    internal_states = @view device_states[local_ix]
    Vm = internal_states[1]
    Vrll = internal_states[2]
    Vr = internal_states[3]
    Vfb = internal_states[4]

    #Define external states for device
    Vt = sqrt(inner_vars[VR_gen_var]^2 + inner_vars[VI_gen_var]^2) # machine's terminal voltage
    Vs = inner_vars[V_pss_var] # PSS output
    Ifd = inner_vars[Xad_Ifd_var] # machine's field current in exciter base

    #Get parameters
    avr = PSY.get_avr(dynamic_device)
    Tr = PSY.get_Tr(avr)
    Vi_min, Vi_max = PSY.get_Vi_lim(avr)
    Tc = PSY.get_Tc(avr)
    Tb = PSY.get_Tb(avr)
    Ka = PSY.get_Ka(avr)
    Ta = PSY.get_Ta(avr)
    Vr_min, Vr_max = PSY.get_Vr_lim(avr)
    Kc = PSY.get_Kc(avr)
    Kf = PSY.get_Kf(avr)
    Tf = PSY.get_Tf(avr)

    #Compute auxiliary parameters
    V_ref = V0_ref + Vs

    # Compute block derivatives
    _, dVm_dt = low_pass_mass_matrix(Vt, Vm, 1.0, Tr)
    y_hp, dVfb_dt = high_pass(Vr, Vfb, Kf, Tf)
    y_ll, dVrll_dt =
        lead_lag_mass_matrix(clamp(V_ref - Vm - y_hp, Vi_min, Vi_max), Vrll, 1.0, Tc, Tb)
    _, dVr_dt = low_pass_mass_matrix(y_ll, Vr, Ka, Ta)

    #Compute 4 States AVR ODE:
    output_ode[local_ix[1]] = dVm_dt
    output_ode[local_ix[2]] = dVrll_dt
    output_ode[local_ix[3]] = dVr_dt
    output_ode[local_ix[4]] = dVfb_dt

    #Update inner_vars
    Vf = clamp(Vr, Vt * Vr_min - Kc * Ifd, Vt * Vr_max - Kc * Ifd)
    inner_vars[Vf_var] = Vf
    return
end

function mdl_avr_ode!(
    device_states::AbstractArray,
    output_ode::AbstractArray,
    inner_vars::AbstractArray,
    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{M, S, PSY.EXAC1, TG, P}},
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}

    #Obtain references
    V_ref = get_V_ref(dynamic_device)

    #Obtain avr
    avr = PSY.get_avr(dynamic_device)

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(dynamic_device, typeof(avr))

    #Define inner states for component
    internal_states = @view device_states[local_ix]
    Vm = internal_states[1]
    Vr1 = internal_states[2]
    Vr2 = internal_states[3]
    Ve = internal_states[4]
    Vr3 = internal_states[5]

    #Define external states for device
    V_th = sqrt(inner_vars[VR_gen_var]^2 + inner_vars[VI_gen_var]^2) # machine's terminal voltage
    Vs = inner_vars[V_pss_var] # PSS output
    Xad_Ifd = inner_vars[Xad_Ifd_var] # machine's field current in exciter base

    #Get parameters
    Tr = PSY.get_Tr(avr)
    Tb = PSY.get_Tb(avr)
    Tc = PSY.get_Tc(avr)
    Ka = PSY.get_Ka(avr)
    Ta = PSY.get_Ta(avr)
    Vr_min, Vr_max = PSY.get_Vr_lim(avr)
    Te = PSY.get_Te(avr) # Te > 0
    Kf = PSY.get_Kf(avr)
    Tf = PSY.get_Tf(avr) # Tf > 0
    Kc = PSY.get_Kc(avr)
    Kd = PSY.get_Kd(avr)
    Ke = PSY.get_Ke(avr)

    #Obtain saturation
    Se = saturation_function(avr, Ve)

    #Compute auxiliary parameters
    I_N = Kc * Xad_Ifd / Ve
    V_FE = Kd * Xad_Ifd + Ke * Ve + Se * Ve
    Vf = Ve * rectifier_function(I_N)

    #Compute block derivatives
    _, dVm_dt = low_pass_mass_matrix(V_th, Vm, 1.0, Tr)
    V_F, dVr3_dt = high_pass(V_FE, Vr3, Kf, Tf)
    V_in = V_ref + Vs - Vm - V_F
    y_ll, dVr1_dt = lead_lag_mass_matrix(V_in, Vr1, 1.0, Tc, Tb)
    y_Vr, dVr2_dt = low_pass_nonwindup_mass_matrix(y_ll, Vr2, Ka, Ta, Vr_min, Vr_max)
    _, dVe_dt = integrator_nonwindup(y_Vr - V_FE, Ve, 1.0, Te, 0.0, Inf)

    #Compute 4 States AVR ODE:
    output_ode[local_ix[1]] = dVm_dt
    output_ode[local_ix[2]] = dVr1_dt
    output_ode[local_ix[3]] = dVr2_dt
    output_ode[local_ix[4]] = dVe_dt
    output_ode[local_ix[5]] = dVr3_dt

    #Update inner_vars
    inner_vars[Vf_var] = Vf
    return
end


text:
@with_kw struct TrackArtist
    id::UInt = 0
    name::String = ""
    resource_url::String = ""
    role::Union{String, Nothing} = nothing
    tracks::Union{String, Nothing} = nothing
    anv::Union{String, Nothing} = nothing
    join::Union{String, Nothing} = nothing
    thumbnail_url::Union{String, Nothing} = nothing
end

@with_kw struct Contributor
    username::String = ""
    resource_url::String = ""
end

@with_kw struct Rating
    average::Real = 0.0
    count::UInt = 0
end

@with_kw struct Community
    data_quality::String = ""
    want::UInt = 0
    rating::Union{Rating, Nothing} = nothing
    status::String = ""
    have::UInt = 0
    contributors::Union{Vector{Contributor}, Nothing} = nothing
    submitter::Union{Contributor, Nothing} = nothing
end

@with_kw struct Company
    id::UInt = 0
    name::String = ""
    resource_url::String = ""
    entity_type::Union{String, Nothing} = nothing
    entity_type_name::Union{String, Nothing} = nothing
    catno::Union{UInt, Nothing} = nothing
    thumbnail_url::Union{String, Nothing} = nothing
end

@with_kw struct Identifier
    value::String = ""
    type::String = ""
    description::Union{String, Nothing} = nothing
end

@with_kw struct Format
    name::String = ""
    descriptions::Vector{String} = String[]
    qty::String = ""
    text::Union{String, Nothing} = nothing
end

@with_kw struct Image
    uri::Union{String, Nothing} = nothing
    uri150::Union{String, Nothing} = nothing
    resource_url::Union{String, Nothing} = nothing
    height::UInt = 0
    width::UInt = 0
    type::String = ""
end

@with_kw struct Label
    id::UInt = 0
    name::String = ""
    resource_url::String = ""
    catno::String = ""
    entity_type::String = ""
    entity_type_name::String = ""
    thumbnail_url::Union{String, Nothing} = nothing
end

@with_kw struct Track
    position::Union{String, Nothing} = nothing
    title::String = ""
    type_::String = ""
    duration::Union{String, Nothing} = nothing
    artists::Union{Vector{TrackArtist}, Nothing} = nothing
    extraartists::Union{Vector{TrackArtist}, Nothing} = nothing
end

@with_kw struct Video
    duration::UInt = 0
    uri::String = ""
    embed::Bool = false
    title::String = ""
    description::Union{String, Nothing} = nothing
end

@with_kw struct Member
    id::UInt = 0
    name::String = ""
    resource_url::String = ""
    thumbnail_url::Union{String, Nothing} = nothing
    active::Union{Bool, Nothing} = nothing
end

@with_kw struct Artist
    id::UInt = 0
    name::String = ""
    resource_url::String = ""
    uri::String = ""
    releases_url::String = ""

    aliases::Union{Vector{Member}, Nothing} = nothing
    namevariations::Union{Vector{String}, Nothing} = nothing
    realname::Union{String, Nothing} = nothing
    urls::Union{Vector{String}, Nothing} = nothing
    data_quality::Union{String, Nothing} = nothing
    images::Union{Vector{Image}, Nothing} = nothing
    members::Union{Vector{Member}, Nothing} = nothing
    profile::Union{String, Nothing} = nothing
end

@with_kw struct Release
    id::UInt = 0
    title::String = ""
    resource_url::String = ""

    artists::Union{Vector{TrackArtist}, Nothing} = nothing
    artists_sort::Union{String, Nothing} = nothing
    data_quality::Union{String, Nothing} = nothing
    thumb::Union{String, Nothing} = nothing
    community::Union{Community, Nothing} = nothing
    companies::Union{Vector{Company}, Nothing} = nothing
    country::Union{String, Nothing} = nothing
    date_added::Union{ZonedDateTime, Nothing} = nothing
    date_changed::Union{ZonedDateTime, Nothing} = nothing
    estimated_weight::Union{UInt, Nothing} = nothing
    extraartists::Union{Vector{TrackArtist}, Nothing} = nothing
    format_quantity::Union{Int64, Nothing} = nothing
    formats::Union{Vector{Format}, Nothing} = nothing
    genres::Union{Vector{String}, Nothing} = nothing
    identifiers::Union{Vector{Identifier}, Nothing} = nothing
    images::Union{Vector{Image}, Nothing} = nothing
    labels::Union{Vector{Label}, Nothing} = nothing
    lowest_price::Union{Real, Nothing} = nothing
    master_id::Union{UInt, Nothing} = nothing
    master_url::Union{String, Nothing} = nothing
    notes::Union{String, Nothing} = nothing
    num_for_sale::Union{UInt, Nothing} = nothing
    released::Union{String, Nothing} = nothing
    released_formatted::Union{String, Nothing} = nothing
    series::Union{Vector{Label}, Nothing} = nothing
    status::Union{String, Nothing} = nothing
    styles::Union{Vector{String}, Nothing} = nothing
    tracklist::Union{Vector{Track}} = nothing
    uri::Union{String, Nothing} = nothing
    videos::Union{Vector{Video}, Nothing} = nothing
    year::Union{UInt, Nothing} = nothing
end


text:
using Gaston
# using PyCall
if VERSION >= v"0.7-"
    using DelimitedFiles
end

function loadcsv(inputcsv)
    contents = readdlm(inputcsv, ',', Float64, '\n'; header = true)
    return contents
end

function coldata(inputcsv, theset)
    contents = loadcsv(inputcsv)
    return contents[1][:, theset]
end

set(axis="loglog", plotstyle="linespoints", linewidth=2, pointsize = 2, xrange="[0.002:0.1]", xlabel = "Element size", ylabel = "Approximate error", grid="on", title = "")

# inputcsv = "All_EBC_2dir_T10_default_Stress.CSV"
# x = coldata(inputcsv, 1)
# y = coldata(inputcsv, 4)
# plot(x, abs.(y), legend = "T10", color = "black", marker = "x")

# inputcsv = "All_EBC_2dir_MST10_extrapmean_Stress.CSV"
# x = coldata(inputcsv, 1)
# y = coldata(inputcsv, 4)
# plot!(x, abs.(y), legend="MSOE", color = "black", marker = "edmd")

# inputcsv = "All_EBC_2dir_MST10_extraptrend_Stress.CSV"
# x = coldata(inputcsv, 1)
# y = coldata(inputcsv, 4)
# plot!(x, abs.(y), legend="TBE", color = "black", marker = "ecircle")

f = figure()
inputcsv = "All_EBC_2dir_MSH8_extraptrend_trapezoidal_Stress.CSV"
x = coldata(inputcsv, 1)
y = coldata(inputcsv, 4)
plot!(x, abs.(y), legend = "Trapezoidal", color = "black", marker = "x")

inputcsv = "All_EBC_2dir_MSH8_extraptrend_Stress.CSV"
x = coldata(inputcsv, 1)
y = coldata(inputcsv, 4)
plot!(x, abs.(y), legend = "Gauss", color = "black", marker = "ecircle")

figure(f)
# inputcsv = "All_EBC_2dir_MST10_extrapmean_Stress.CSV"
# x = coldata(inputcsv, 1)
# y = coldata(inputcsv, 4)
# plot!(x, abs.(y), legend="MSOE", color = "black", marker = "edmd")

# inputcsv = "All_EBC_2dir_MST10_extraptrend_Stress.CSV"
# x = coldata(inputcsv, 1)
# y = coldata(inputcsv, 4)
# plot!(x, abs.(y), legend="TBE", color = "black", marker = "ecircle")

# julia> include("All_EBC_2dir_examples.jl"); All_EBC_2dir_examples.allrun()
# WARNING: replacing module All_EBC_2dir_examples.
# #####################################################
# # All_EBC_2dir_MST10_conv
# Fiber-reinforced block: MST10

# Mesh: na, nb, nts = 3, 3, [3]
# count(fens) = 343
# Mesh: na, nb, nts = 6, 6, [6]
# count(fens) = 2197
# Mesh: na, nb, nts = 12, 12, [12]
# count(fens) = 15625
# Mesh: na, nb, nts = 24, 24, [24]
# count(fens) = 117649

# Stress RMS error
# Normalized Approximate Error = [0.0192602, 0.00861774, 0.00354827]
# Linear log-log fit: p = [1.22022, 0.214338]
# Wrote All_EBC_2dir_MST10_extrapmean_Stress.CSV

# Displacement RMS error
# Normalized Approximate Error = [0.00508317, 0.00143038, 0.000373055]
# Linear log-log fit: p = [1.88413, 1.13915]
# Wrote All_EBC_2dir_MST10_extrapmean_Displ.CSV
# Mesh: na, nb, nts = 3, 3, [3]
# count(fens) = 343
# Mesh: na, nb, nts = 6, 6, [6]
# count(fens) = 2197
# Mesh: na, nb, nts = 12, 12, [12]
# count(fens) = 15625
# Mesh: na, nb, nts = 24, 24, [24]
# count(fens) = 117649

# Stress RMS error
# Normalized Approximate Error = [0.0167626, 0.00613563, 0.00220595]
# Linear log-log fit: p = [1.46289, 0.889953]
# Wrote All_EBC_2dir_MST10_extraptrend_Stress.CSV

# Displacement RMS error
# Normalized Approximate Error = [0.00508317, 0.00143038, 0.000373055]
# Linear log-log fit: p = [1.88413, 1.13915]
# Wrote All_EBC_2dir_MST10_extraptrend_Displ.CSV
# Mesh: na, nb, nts = 3, 3, [3]
# count(fens) = 343
# Mesh: na, nb, nts = 6, 6, [6]
# count(fens) = 2197
# Mesh: na, nb, nts = 12, 12, [12]
# count(fens) = 15625
# Mesh: na, nb, nts = 24, 24, [24]
# count(fens) = 117649

# Stress RMS error
# Normalized Approximate Error = [0.017322, 0.00778243, 0.00321]
# Linear log-log fit: p = [1.21598, 0.0942623]
# Wrote All_EBC_2dir_MST10_default_Stress.CSV

# Displacement RMS error
# Normalized Approximate Error = [0.00508317, 0.00143038, 0.000373055]
# Linear log-log fit: p = [1.88413, 1.13915]
# Wrote All_EBC_2dir_MST10_default_Displ.CSV
# Done
# #####################################################
# # All_EBC_2dir_MST10_conv
# Fiber-reinforced block: MST10

# Mesh: na, nb, nts = 3, 3, [3]
# count(fens) = 64
# Mesh: na, nb, nts = 6, 6, [6]
# count(fens) = 343
# Mesh: na, nb, nts = 12, 12, [12]
# count(fens) = 2197
# Mesh: na, nb, nts = 24, 24, [24]
# count(fens) = 15625

# Stress RMS error
# Normalized Approximate Error = [0.0288266, 0.0128894, 0.00580383]
# Linear log-log fit: p = [1.15616, 0.384708]
# Wrote All_EBC_2dir_MST10_extrapmean_Stress.CSV

# Displacement RMS error
# Normalized Approximate Error = [0.00947204, 0.00264193, 0.000695543]
# Linear log-log fit: p = [1.88373, 1.75715]
# Wrote All_EBC_2dir_MST10_extrapmean_Displ.CSV
# Mesh: na, nb, nts = 3, 3, [3]
# count(fens) = 64
# Mesh: na, nb, nts = 6, 6, [6]
# count(fens) = 343
# Mesh: na, nb, nts = 12, 12, [12]
# count(fens) = 2197
# Mesh: na, nb, nts = 24, 24, [24]
# count(fens) = 15625

# Stress RMS error
# Normalized Approximate Error = [0.0253678, 0.0103545, 0.00442978]
# Linear log-log fit: p = [1.25885, 0.599482]
# Wrote All_EBC_2dir_MST10_extraptrend_Stress.CSV

# Displacement RMS error
# Normalized Approximate Error = [0.00947204, 0.00264193, 0.000695543]
# Linear log-log fit: p = [1.88373, 1.75715]
# Wrote All_EBC_2dir_MST10_extraptrend_Displ.CSV
# Mesh: na, nb, nts = 3, 3, [3]
# count(fens) = 64
# Mesh: na, nb, nts = 6, 6, [6]
# count(fens) = 343
# Mesh: na, nb, nts = 12, 12, [12]
# count(fens) = 2197
# Mesh: na, nb, nts = 24, 24, [24]
# count(fens) = 15625

# Stress RMS error
# Normalized Approximate Error = [0.0288266, 0.0128894, 0.00580383]
# Linear log-log fit: p = [1.15616, 0.384708]
# Wrote All_EBC_2dir_MST10_default_Stress.CSV

# Displacement RMS error
# Normalized Approximate Error = [0.00947204, 0.00264193, 0.000695543]
# Linear log-log fit: p = [1.88373, 1.75715]
# Wrote All_EBC_2dir_MST10_default_Displ.CSV
# Done
# #####################################################
# # All_EBC_2dir_T10_conv
# Fiber-reinforced block: T10

# Mesh: na, nb, nts = 3, 3, [3]
# count(fens) = 343
# Mesh: na, nb, nts = 6, 6, [6]
# count(fens) = 2197
# Mesh: na, nb, nts = 12, 12, [12]
# count(fens) = 15625
# Mesh: na, nb, nts = 24, 24, [24]
# count(fens) = 117649

# Stress RMS error
# Normalized Approximate Error = [0.0129395, 0.00559644, 0.00226381]
# Linear log-log fit: p = [1.25748, -0.0593887]
# Wrote All_EBC_2dir_T10_default_Stress.CSV

# Displacement RMS error
# Normalized Approximate Error = [0.00289237, 0.000501187, 8.09591e-5]
# Linear log-log fit: p = [2.57946, 2.93927]
# Wrote All_EBC_2dir_T10_default_Displ.CSV
# Done
# #####################################################
# # All_EBC_2dir_T10_conv
# Fiber-reinforced block: T10

# Mesh: na, nb, nts = 3, 3, [3]
# count(fens) = 64
# Mesh: na, nb, nts = 6, 6, [6]
# count(fens) = 343
# Mesh: na, nb, nts = 12, 12, [12]
# count(fens) = 2197
# Mesh: na, nb, nts = 24, 24, [24]
# count(fens) = 15625

# Stress RMS error
# Normalized Approximate Error = [0.023016, 0.0101339, 0.00453253]
# Linear log-log fit: p = [1.17213, 0.212446]
# Wrote All_EBC_2dir_T10_default_Stress.CSV

# Displacement RMS error
# Normalized Approximate Error = [0.00929339, 0.00262218, 0.00069389]
# Linear log-log fit: p = [1.87171, 1.69831]
# Wrote All_EBC_2dir_T10_default_Displ.CSV
# Done
# true

# julia> include("All_EBC_2dir_plots.jl")
# 1

# julia> include("All_EBC_2dir_plots.jl")
# 1

# julia> include("All_EBC_2dir_plots.jl")
# 1


text:
# Utility functions for SICR scripts to run set of models in tuple modeltypes on region

using Distributed
using DelimitedFiles
using Dates
using LinearAlgebra
using StatsBase
using Statistics
using Distributions
using DifferentialEquations
using DataFrames
using CSV

const txtstr = ".txt"

include("SICRPoisson.jl")

function runmcmc(total,datestrin,datestrout,temp=1.,N=1e7,pinstr= "pl",modeltype= ["rqfmt"])

    println(total,' ',temp)
    println(pinstr)

    pathin, pathout = makepaths(datestrin,datestrout)
    datafiles,regions = getdatafiles("data")
    println(datestrin,"->",datestrout)

    pstatstr = pathout * "pstats" * txtstr
    fpstat = open(pstatstr,"w")
    for j in eachindex(datafiles)
        datafile = datafiles[j]
        df = DataFrame(CSV.File("data/" * datafile))
        println(regions[j])
        for i in eachindex(modeltype)
            p,pfit,cols,data = mcmcargs(modeltype[i],regions[j],df,pathin,pinstr)
            # println(modeltype[i])
            outtuple = @spawn mcmc(data,cols,p,pfit,N,total,linearrd!,temp)

            model = modeltype[i]
            region = regions[j]
            writemeasures(outtuple,model,region,pstatstr)
            writecovariances(outtuple,model,region,pathout)
        end
    end
end

function runmcmcMP(total,datestrin,datestrout,temp=1.,N=1e7,pinstr= "pl",modeltype= ["rqfmt"])

    println(total,' ',temp)
    println(pinstr)

    pathin, pathout = makepaths(datestrin,datestrout)
    datafiles,regions = getdatafiles("data")
    println(datestrin,"->",datestrout)

    out = Array{Future,2}(undef,length(modeltype),length(datafiles))

    for j in eachindex(datafiles)
        datafile = datafiles[j]
        df = DataFrame(CSV.File("data/" * datafile))
        println(regions[j])
        for i in eachindex(modeltype)
            p,pfit,cols,data = mcmcargs(modeltype[i],regions[j],df,pathin,pinstr)
            # println(modeltype[i])
            out[i,j] = @spawn mcmc(data,cols,p,pfit,N,total,linearrd!,temp)
        end
    end

    pstatstr = pathout * "pstats" * txtstr
    fpstat = open(pstatstr,"w")
    for j in eachindex(datafiles), i in eachindex(modeltype)
        outtuple = fetch(out[i,j])
        model = modeltype[i]
        region = regions[j]
        writemeasures(outtuple,model,region,pstatstr)
        writecovariances(outtuple,model,region,pathout)
    end
end

# retries list of data files in folder datafolder
function getdatafiles(datafolder)
    datafiles = readdir("data")
    datafiles = datafiles[occursin.("covid_timeseries",datafiles)]
    regions = Array{String,1}(undef,0)
    for datafile in datafiles
        push!(regions,replace(replace(datafile,".csv"=>""),"covid_timeseries_"=>""))
    end
    return datafiles,regions
end

# returns prior for parameters for original r and u models
# param order: p =  beta,sigmac,sigmar,sigmad,q,f,mbase,mlocation,mrate,cmax,c50
function getpfit(modelname)
    if modelname == "r"
        return  [1;2;3;4;12]
    end
    if modelname == "rq"
        return [1;2;3;4;5;12]
    end
    if modelname == "rf"
        return [1;2;3;4;6;12]
    end
    if modelname == "rqf"
        return [1;2;3;4;5;6;12]
    end
    if modelname == "rm"
        return [1;2;3;4;7;8;9;12]
    end
    if modelname == "rqm"
        return [1;2;3;4;5;7;8;9;12]
    end
    if modelname == "rfm"
        return [1;2;3;4;6;7;8;9;12]
    end
    if modelname == "rqfm"
        return [1;2;3;4;5;6;7;8;9;12]
    end
    if modelname == "rmt"
        return [1;2;3;4;5;6;7;8;9;10;11;12]
    end
    if modelname == "rqmt" || modelname == "rfm"
        return [1;2;3;4;5;6;7;8;9;10;11;12]
    end
    if modelname == "rfmt"
        return [1;2;3;4;5;6;7;8;9;10;11;12]
    end
    if modelname == "rqfmt"
        return [1;2;3;4;5;6;7;8;9;10;11;12]
    end
    if modelname == "u"
        return [.25,.1,.1,.1]
    end
    if modelname == "uq" || modelname == "uf"
        return [.25,.1,.1,.1,.1]
    end
    if modelname == "um"
        return [.25,.1,.1,.1,5.,.1]
    end
    if modelname == "uqm" || modelname == "ufm"
        return [.25,.1,.1,.1,.1,5.,.1]
    end
    if modelname == "uqfm"
        return [.25,.1,.1,.1,.1,.1,5.,.1]
    end
end

# name = Dict(sicr! => "r", sicrq! => "rq", sicrf! => "rf", sicrqf! => "rqf", sicrm! => "rm", sicrqm! => "rqm", sicrfm! => "rfm", sicrqfm! => "rqfm",sicu! => "u", sicuq! => "uq", sicuf! => "uf", sicuq! => "uq", sicuqf! => "uqf", sicum! => "um", sicuqm! => "uqm", sicufm! => "ufm", sicuqfm! => "uqfm")

function makepaths(datestrin,datestrout)
    if ~isdir("Results")
        mkdir("Results")
    end
    pathin = "Results/$datestrin/"
    pathout = "Results/$datestrout/"
    if ~isdir(pathout)
        mkdir(pathout)
    end
    return pathin, pathout
end

# write out parameter and R0 statistics
function writemeasures(outtuple,name,region,pstatstr)
    pstats,R0 = measures(outtuple[2],outtuple[5])
    f = open(pstatstr,"a")
    writedlm(f,[region * "_" * name])
    writedlm(f,[outtuple[4] R0[1] R0[2] R0[3] R0[4]'])
    for ptup in pstats
        writedlm(f, [ptup[1] ptup[2] ptup[3] ptup[4]'])
    end
    close(f)
end

# write out max likelihood params
function writemaxll(outtuple,name,region,pathout)
    mlstr = pathout * "maxll" * "_" * region * "_" * name * txtstr
    f=open(mlstr,"w")
    f=open(mlstr,"a")
    writedlm(f,[outtuple[4];outtuple[5]])
    close(f)
end

# write out all loglikelihood mcmc samples
function writellsamples(outtuple,name,region,pathout)
    llstr = pathout * "llsamples" * "_" * region * "_" * name * txtstr
    f=open(llstr,"w")
    f=open(llstr,"a")
    writedlm(f,outtuple[1])
    close(f)
end

# write out last parameter in run
function writelastparam(outtuple,name,region,pathout)
    pstr = pathout * "pl" * "_" * region * "_" * name * txtstr
    f=open(pstr,"w")
    f=open(pstr,"a")
    writedlm(f,[outtuple[4];outtuple[2][end,:]])
    close(f)
end

# write out all parameters parameter in run
function writeparamsamples(outtuple,name,region,pathout)
    pstr = pathout * "paramsamples" * "_" * region * "_" * name * txtstr
    f=open(pstr,"w")
    f=open(pstr,"a")
    for i in 1:size(outtuple[2],1)
        writedlm(f,[outtuple[2][i,:]])
    end
    close(f)
end

function writecovariances(outtuple,name,region,pathout)
    pstr = pathout * "covcor" * "_" * region * "_" * name * ".csv"
    f=open(pstr,"w")
    f=open(pstr,"a")
    writedlm(f,cov(outtuple[2]),',')
    writedlm(f,cor(outtuple[2]),',')
    close(f)
end

# output arguments used in mcmc function
function mcmcargs(name,region,df,pathin,pinstr)
    colsr = [5,6,7]
    colsu = [5,6]
    if occursin("u",name)
        data = sicuprep(df)
        cols = colsu
    else
        data = df
        cols = colsr
    end

    infile = pathin * pinstr * "_" * region * "_" * name * txtstr
    if isfile(infile) && ~isempty(read(infile))
        p = readdlm(infile)[2:end]
        pfit = getpfit(name)
    else
        pfit = getpfit(name)
        p = [.25,.1,.01,.01,.01,1.,.1,1.,1.,.1,10.,1.]
        println(name,": No prior")
    end
    return p,pfit,cols,data
end

max_stars_repo_name: cannin/covid-sicr-test

max_issues_repo_name: cannin/covid-sicr-test

max_forks_repo_name: cannin/covid-sicr-test


text:
using CALiPPSO
using DelimitedFiles, StatsBase, BenchmarkTools, HiGHS

const tol_optimality = CALiPPSO.default_tol_optimality

const highs_args = nothing
const highs_attributes = Dict("small_matrix_value"=>0.1*tol_optimality, "primal_feasibility_tolerance" => tol_optimality, "dual_feasibility_tolerance" => tol_optimality, "solver" => "ipm", "ipm_optimality_tolerance"=>tol_optimality,  "threads" => CALiPPSO.max_threads, "parallel"=>"on", "output_flag"=>false)
printstyled("\n______________________________________________________________________________________\n", color=:yellow)
precompile_main_function(HiGHS, highs_attributes, highs_args)

test_bench = @benchmarkable sin(1)
run(test_bench)

const N_LS = 1024 # size of configurations obtained after LS compression
const ds = [3, 4, 5] # dimensions to use
const L = 1.0 # size of each side of the system's volume

const io = IOContext(stdout)

# Testing CALiPPSO in systems of different dimensions.
# The initial conditions are configurations compressed with Lubachevsky--Stillinger up to a (reduced) pressure p>>10
for d in ds
    σ0 = readdlm("Centers-after-LS--N-$(N_LS)--d-$d.dat", skipstart=2)[1, 1] # initial diameter of the particles
    r0 = 0.5*σ0 # initial radius of particles
    Xs0 = PeriodicVectors(Matrix(readdlm("Centers-after-LS--N-$(N_LS)--d-$d.dat", skipstart=3)'), L) # initial particles' positions

    printstyled("\n\n++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n", color=:blue)
    printstyled("Using CALiPPSO to jam a system of N = ", N_LS, " in d = ", d, "\t Initial φ = ", packing_fraction(d, r0, N_LS, L), "\n", bold=:true, color=:blue)
    printstyled("++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n", color=:blue)

    benchmark_iLP = @benchmarkable produce_jammed_configuration!(Xs, r; verbose=false, solver=HiGHS, solver_attributes=highs_attributes, solver_args=highs_args) setup=(Xs = copy($Xs0); r=copy($r0) )
    res_benchmark = run(benchmark_iLP, seconds = 3600, samples=20)

    printstyled("\n\n - - - o o o - - - o o o - - - o o o - - - o o o - - - o o o - - - o o o - - - o o o - - - \n\n", color=:green)
    println("Results of benchmarking CALiPPSO with 20 samples:\n\n")
    show(io, MIME("text/plain"), res_benchmark)
    printstyled("\n\n - - - o o o - - - o o o - - - o o o - - - o o o - - - o o o - - - o o o - - - o o o - - - \n\n", color=:green)

    # io2 = open("Benchmarks--d-$d-after-LS--changing-solver.txt", "w")
    io2 = open("Benchmarks--d-$d-after-LS--base.txt", "w")
    show(io2, MIME("text/plain"), res_benchmark)
    close(io2)
end

max_stars_repo_name: rdhr/CALiPPSO

max_issues_repo_name: rdhr/CALiPPSO

max_forks_repo_name: rdhr/CALiPPSO


text:
function Euler6(N::Int64)

	println( div(N*(N+1)*(3*N+2)*(N-1),12))

end

@time Euler6(100)

max_stars_repo_name: gribeill/ProjectEuler

max_issues_repo_name: gribeill/ProjectEuler

max_forks_repo_name: gribeill/ProjectEuler


text:
const RectilinearPointLoad{dim, T, N, M} = Union{PointLoadCantilever{dim, T, N, M}, HalfMBB{dim, T, N, M}, LBeam{T, N, M}}

@params struct ElementMatrix{T, TM <: AbstractMatrix{T}} <: AbstractMatrix{T}
    matrix::TM
    mask
    meandiag::T
end
ElementMatrix(matrix, mask) = ElementMatrix(matrix, mask, sumdiag(matrix)/size(matrix, 1))

const StaticMatrices{m, T} = Union{StaticMatrix{m, m, T}, Symmetric{T, <:StaticMatrix{m, m, T}}}
@generated function sumdiag(K::StaticMatrices{m,T}) where {m,T}
    return reduce((ex1,ex2) -> :($ex1 + $ex2), [:(K[$j,$j]) for j in 1:m])
end

Base.size(m::ElementMatrix) = size(m.matrix)
Base.getindex(m::ElementMatrix, i...) = m.matrix[i...]

rawmatrix(m::ElementMatrix) = m.matrix
rawmatrix(m::Symmetric{T, <:ElementMatrix{T}}) where {T} = Symmetric(m.data.matrix)
@generated function bcmatrix(m::ElementMatrix{T, TM}) where {dim, T, TM <: StaticMatrix{dim, dim, T}}
    expr = Expr(:tuple)
    for j in 1:dim, i in 1:dim
        push!(expr.args, :(ifelse(m.mask[$i] && m.mask[$j], m.matrix[$i,$j], zero(T))))
    end
    return :($(Expr(:meta, :inline)); $TM($expr))
end
@generated function bcmatrix(m::Symmetric{T, <:ElementMatrix{T, TM}}) where {dim, T, TM <: StaticMatrix{dim, dim, T}}
    expr = Expr(:tuple)
    for j in 1:dim, i in 1:dim
        push!(expr.args, :(ifelse(m.data.mask[$i] && m.data.mask[$j], m.data.matrix[$i,$j], zero(T))))
    end
    return :($(Expr(:meta, :inline)); Symmetric($TM($expr)))
end

@params struct ElementFEAInfo{dim, T}
    Kes::AbstractVector{<:AbstractMatrix{T}}
    fes::AbstractVector{<:AbstractVector{T}}
    fixedload::AbstractVector{T}
    cellvolumes::AbstractVector{T}
    cellvalues::CellValues{dim, T}
    facevalues::FaceValues{<:Any, T}
    metadata::Metadata
    black::AbstractVector
    white::AbstractVector
    varind::AbstractVector{Int}
    cells
end
function ElementFEAInfo(sp, quad_order=2, ::Type{Val{mat_type}}=Val{:Static}) where {mat_type}
    Kes, weights, dloads, cellvalues, facevalues = make_Kes_and_fes(sp, quad_order, Val{mat_type})
    element_Kes = make_element_Kes(Kes, sp.ch.prescribed_dofs, sp.metadata.dof_cells)
    fixedload = Vector(make_cload(sp))
    assemble_f!(fixedload, sp, dloads)
    cellvolumes = get_cell_volumes(sp, cellvalues)
    cells = sp.ch.dh.grid.cells
    ElementFEAInfo(element_Kes, weights, fixedload, cellvolumes, cellvalues, facevalues, sp.metadata, sp.black, sp.white, sp.varind, cells)
end
function make_element_Kes(Kes::AbstractVector{TMorSymm}, bc_dofs, dof_cells) where {N, T, TM <: StaticMatrix{N, N, T}, TMorSymm <: Union{TM, Symmetric{T, TM}}}
    fill_matrix = zero(TM)
    fill_mask = ones(SVector{N, Bool})
    if TMorSymm <: Symmetric
        element_Kes = fill(Symmetric(ElementMatrix(fill_matrix, fill_mask)), length(Kes))
    else
        element_Kes = fill(ElementMatrix(fill_matrix, fill_mask), length(Kes))
    end
    for i in bc_dofs
        d_cells = dof_cells[i]
        for c in d_cells
            (cellid, localdof) = c
            if TMorSymm <: Symmetric
                Ke = element_Kes[cellid].data
            else
                Ke = element_Kes[cellid]
            end
            new_Ke = @set Ke.mask[localdof] = false
            element_Kes[cellid] = Symmetric(new_Ke)
        end
    end
    for e in 1:length(element_Kes)
        if eltype(element_Kes) <: Symmetric
            Ke = element_Kes[e].data
            matrix = Kes[e].data
            Ke = @set Ke.matrix = matrix
            element_Kes[e] = Symmetric(@set Ke.meandiag = sumdiag(Ke.matrix))
        else
            Ke = element_Kes[e]
            matrix = Kes[e]
            Ke = @set Ke.matrix = matrix
            element_Kes[e] = @set Ke.meandiag = sumdiag(Ke.matrix)
        end
    end
    element_Kes
end

function make_element_Kes(Kes::AbstractVector{TM}, bc_dofs, dof_cells) where {T, TM <: AbstractMatrix{T}, TMorSymm <: Union{TM, Symmetric{T, TM}}}
    N = size(Kes[1], 1)
    fill_matrix = zero(TM)
    fill_mask = ones(Bool, N)
    if TM <: Symmetric
        element_Kes = [Symmetric(deepcopy(ElementMatrix(fill_matrix, fill_mask))) for i in 1:length(Kes)]
    else
        element_Kes = [deepcopy(ElementMatrix(fill_matrix, fill_mask)) for i in 1:length(Kes)]
    end
    for i in bc_dofs
        d_cells = dof_cells[i]
        for c in d_cells
            (cellid, localdof) = c
            if TM <: Symmetric
                Ke = element_Kes[cellid].data
            else
                Ke = element_Kes[cellid]
            end
            Ke.mask[localdof] = false
        end
    end
    element_Kes
end

function get_cell_volumes(sp::StiffnessTopOptProblem{dim, T}, cellvalues) where {dim, T}
    dh = sp.ch.dh
    cellvolumes = zeros(T, getncells(dh.grid))
    for (i, cell) in enumerate(CellIterator(dh))
        reinit!(cellvalues, cell)
        cellvolumes[i] = sum(JuAFEM.getdetJdV(cellvalues, q_point) for q_point in 1:JuAFEM.getnquadpoints(cellvalues))
    end
    return cellvolumes
end

mutable struct GlobalFEAInfo{T, TK<:AbstractMatrix{T}, Tf<:AbstractVector{T}}
    K::TK
    f::Tf
end

GlobalFEAInfo(::Type{T}) where T = GlobalFEAInfo{T}()

GlobalFEAInfo() = GlobalFEAInfo{Float64}()

GlobalFEAInfo{T}() where T = GlobalFEAInfo{T, SparseMatrixCSC{T, Int}, Vector{T}}(sparse(zeros(T, 0, 0)), zeros(T, 0))

GlobalFEAInfo(sp::StiffnessTopOptProblem) = GlobalFEAInfo(make_empty_K(sp), make_empty_f(sp))

make_empty_K(sp::StiffnessTopOptProblem) = Symmetric(create_sparsity_pattern(sp.ch.dh))

make_empty_f(sp::StiffnessTopOptProblem{dim, T}) where {dim, T} = zeros(T, ndofs(sp.ch.dh))

function make_Kes_and_fes(problem, quad_order=2)
    make_Kes_and_fes(problem, quad_order, Val{:Static})
end

function make_Kes_and_fes(problem, ::Type{Val{mat_type}}) where mat_type
    make_Kes_and_fes(problem, 2, Val{mat_type})
end

function make_Kes_and_fes(problem, quad_order, ::Type{Val{mat_type}}) where {mat_type}
    T = floattype(problem)
    dim = getdim(problem)
    geom_order = getgeomorder(problem)
    dh = getdh(problem)
    E = getE(problem)
    ν = getν(problem)
    ρ = getdensity(problem)

    refshape = JuAFEM.getrefshape(dh.field_interpolations[1])

    λ = E*ν / ((1 + ν) * (1 - 2*ν))
    μ = E / (2*(1 + ν))
    δ(i,j) = i == j ? T(1) : T(0)
    g(i,j,k,l) = λ*δ(i,j)*δ(k,l) + μ*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k))
    C = SymmetricTensor{4, dim}(g)

    # Shape functions and quadrature rule
    interpolation_space = Lagrange{dim, refshape, geom_order}()
    quadrature_rule = QuadratureRule{dim, refshape}(quad_order)
    cellvalues = CellScalarValues(quadrature_rule, interpolation_space)
    facevalues = FaceScalarValues(QuadratureRule{dim-1, refshape}(quad_order), interpolation_space)

    # Calculate element stiffness matrices
    n_basefuncs = getnbasefunctions(cellvalues)

    Kes, weights = _make_Kes_and_weights(dh, Val{mat_type}, Val{n_basefuncs}, Val{dim*n_basefuncs}, C, ρ, quadrature_rule, cellvalues)
    dloads = _make_dloads(weights, problem, facevalues)

    return Kes, weights, dloads, cellvalues, facevalues
end

const g = [0., 9.81, 0.] # N/kg or m/s^2

function _make_Kes_and_weights(dh::DofHandler{dim, N, T}, ::Type{Val{mat_type}}, ::Type{Val{n_basefuncs}}, ::Type{Val{ndofs_per_cell}}, C, ρ, quadrature_rule, cellvalues) where {dim, N, T, mat_type, n_basefuncs, ndofs_per_cell}
    # Calculate element stiffness matrices
    Kesize = ndofs_per_cell
    nel = getncells(dh.grid)
    body_force = ρ .* g # Force per unit volume

    if !(T === BigFloat)
        if mat_type === :Static || mat_type === :SMatrix
            MatrixType = SMatrix{Kesize, Kesize, T, Kesize^2}
            VectorType = SVector{Kesize, T}
        elseif mat_type === :MMatrix
            MatrixType = MMatrix{Kesize, Kesize, T, Kesize^2}
            VectorType = MVector{Kesize, T}
        else
            MatrixType = Matrix{T}
            VectorType = Vector{T}
        end
    else
        if mat_type === :Static || mat_type === :SMatrix  || mat_type === :MMatrix
            MatrixType = SizedMatrix{Kesize, Kesize, T, Kesize^2}
            VectorType = SizedVector{Kesize, T}
        else
            MatrixType = Matrix{T}
            VectorType = Vector{T}
        end
    end

    if MatrixType <: StaticArray
        Kes = Symmetric{T, MatrixType}[]
        sizehint!(Kes, nel)
        weights = [zeros(VectorType) for i in 1:nel]

        Ke_e = zeros(T, dim, dim)
        fe = zeros(T, Kesize)
        Ke_0 = Matrix{T}(undef, Kesize, Kesize)
        celliterator = CellIterator(dh)
        for (k, cell) in enumerate(celliterator)
            Ke_0 .= 0
            reinit!(cellvalues, cell)
            fe = weights[k]
            for q_point in 1:getnquadpoints(cellvalues)
                dΩ = getdetJdV(cellvalues, q_point)
                for b in 1:n_basefuncs
                    ∇ϕb = shape_gradient(cellvalues, q_point, b)
                    ϕb = shape_value(cellvalues, q_point, b)
                    for d2 in 1:dim
                        fe = @set fe[(b-1)*dim + d2] += ϕb * body_force[d2] * dΩ
                        for a in 1:n_basefuncs
                            ∇ϕa = shape_gradient(cellvalues, q_point, a)
                            Ke_e .= dotdot(∇ϕa, C, ∇ϕb) * dΩ
                            for d1 in 1:dim
                                #if dim*(b-1) + d2 >= dim*(a-1) + d1
                                Ke_0[dim*(a-1) + d1, dim*(b-1) + d2] += Ke_e[d1,d2]
                                #end
                            end
                        end
                    end
                end
            end
            weights[k] = fe
            if MatrixType <: SizedMatrix # Work around because full constructor errors
                push!(Kes, Symmetric(SizedMatrix{Kesize,Kesize,T}(Ke_0)))
            else
                push!(Kes, Symmetric(MatrixType(Ke_0)))
            end
        end
    else
        Kes = let Kesize=Kesize, nel=nel
            [Symmetric(zeros(T, Kesize, Kesize), :U) for i = 1:nel]
        end
        weights = let Kesize=Kesize, nel=nel
            [zeros(T, Kesize) for i = 1:nel]
        end

        Ke_e = zeros(T, dim, dim)

        celliterator = CellIterator(dh)
        for (k, cell) in enumerate(celliterator)
            reinit!(cellvalues, cell)
            fe = weights[k]
            for q_point in 1:getnquadpoints(cellvalues)
                dΩ = getdetJdV(cellvalues, q_point)
                for b in 1:n_basefuncs
                    ∇ϕb = shape_gradient(cellvalues, q_point, b)
                    ϕb = shape_value(cellvalues, q_point, b)
                    for d2 in 1:dim
                        fe[(b-1)*dim + d2] += ϕb * body_force[d2] * dΩ
                        for a in 1:n_basefuncs
                            ∇ϕa = shape_gradient(cellvalues, q_point, a)
                            Ke_e .= dotdot(∇ϕa, C, ∇ϕb) * dΩ
                            for d1 in 1:dim
                                #if dim*(b-1) + d2 >= dim*(a-1) + d1
                                Kes[k].data[dim*(a-1) + d1, dim*(b-1) + d2] += Ke_e[d1,d2]
                                #end
                            end
                        end
                    end
                end
            end
        end
    end
    return Kes, weights
end

function _make_dloads(fes, problem, facevalues)
    dim = getdim(problem)
    N = nnodespercell(problem)
    T = floattype(problem)
    dloads = deepcopy(fes)
    for i in 1:length(dloads)
        if eltype(dloads) <: SArray
            dloads[i] = zero(eltype(dloads))
        else
            dloads[i] .= 0
        end
    end
    pressuredict = getpressuredict(problem)
    dh = getdh(problem)
    grid = dh.grid
    boundary_matrix = grid.boundary_matrix
    cell_coords = zeros(JuAFEM.Vec{dim, T}, N)
    n_basefuncs = getnbasefunctions(facevalues)
    for k in keys(pressuredict)
        t = -pressuredict[k] # traction = negative the pressure
        faceset = getfacesets(problem)[k]
        for (cellid, faceid) in faceset
            boundary_matrix[faceid, cellid] || throw("Face $((cellid, faceid)) not on boundary.")
            fe = dloads[cellid]
            getcoordinates!(cell_coords, grid, cellid)
            reinit!(facevalues, cell_coords, faceid)
            for q_point in 1:getnquadpoints(facevalues)
                dΓ = getdetJdV(facevalues, q_point) # Face area
                normal = getnormal(facevalues, q_point) # Nomral vector at quad point
                for i in 1:n_basefuncs
                    ϕ = shape_value(facevalues, q_point, i) # Shape function value
                    for d = 1:dim
                        if fe isa SArray
                            fe = @set fe[(i-1)*dim + d] += ϕ * t * normal[d] * dΓ
                        else
                            fe[(i-1)*dim + d] += ϕ * t * normal[d] * dΓ
                        end
                    end
                end
            end
            dloads[cellid] = fe
        end
    end

    return dloads
end

function make_cload(problem)
    T = floattype(problem)
    dim = getdim(problem)
    cloads = getcloaddict(problem)
    dh = getdh(problem)
    metadata = getmetadata(problem)
    node_dofs = metadata.node_dofs
    inds = Int[]
    vals = T[]
    for nodeidx in keys(cloads)
        for (dofidx, force) in enumerate(cloads[nodeidx])
            if force != 0
                dof = node_dofs[(nodeidx-1)*dim+dofidx]
                push!(inds, dof)
                push!(vals, force)
            end
        end
    end
    return sparsevec(inds, vals, ndofs(dh))
end


text:
maze = zeros(Int8, 8, 6)
#maze = [i+j*width  for j = 0:height for i = 0:width]

function generate_maze(width,height,arg...)
    maze = zeros(Int8,width,height)
    maze = [(i+j*width)%10 for i = 0:width-1,j = 0:height-1]
    return (m=maze,w=width,h=height,a=arg)
end

function print_maze(maze)
    for r = 1:size(maze,2)
        txt = ""
        for c in maze[:,r]
            txt *= "$c"
        end
        println(txt)
    end
end

function save_maze(maze,file)
    file = open(file,"w")
    for r = 1:size(maze,2)
        txt = ""
        for c in maze[:,r]
            txt *= "$c"
        end
        txt *= "\n"
        write(file,txt)
    end
    close(file)
end

res = generate_maze(8,6,"dungeon maze")
println(res)
print_maze(res.m)
save_maze(res.m,"maze.txt")

using Sockets
ip = "www.baidu.com"
connect(ip,80)
getaddrinfo(ip)
#print(maze[2:3,5:7])
cmd = `echo hello`
run(cmd)

max_stars_repo_name: fenixnix/Ascii

max_issues_repo_name: fenixnix/GoAscii

max_forks_repo_name: fenixnix/GoAscii


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -31.730360167469108 <= q <= 25.84895899492546)

                     add_NL_constraint(m, :(sigmoid(0.0048766473277672695 + 0.6370542772950833*sigmoid(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5])) + 0.7379502189935647*sigmoid(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5])) + -0.017781498731094914*sigmoid(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5])) + -0.30964053180650675*sigmoid(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5])) + -0.19693223852539354*sigmoid(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))) + sigmoid(0.6750746912592471 + -0.6689814871888631*sigmoid(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5])) + 0.08174772592079638*sigmoid(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5])) + 0.43040165086652493*sigmoid(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5])) + 0.6360324545158482*sigmoid(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5])) + 0.2590621611090005*sigmoid(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))) + sigmoid(-0.5603004507305789 + -0.6080760485894823*sigmoid(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5])) + -0.17827205367530619*sigmoid(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5])) + -0.3024929990799108*sigmoid(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5])) + 0.2876553916986979*sigmoid(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5])) + 0.1394466763243627*sigmoid(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))) + sigmoid(-0.3104079613918116 + -0.8967757209559282*sigmoid(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5])) + 0.9554659562201504*sigmoid(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5])) + 0.4284452595185697*sigmoid(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5])) + -0.8546417933790886*sigmoid(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5])) + -0.9550224035691546*sigmoid(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))) + sigmoid(0.11524778075403264 + -0.08842997013200282*sigmoid(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5])) + -0.2446780712327863*sigmoid(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5])) + 0.5457935761049009*sigmoid(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5])) + -0.1673283141057702*sigmoid(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5])) + 0.6343389483483679*sigmoid(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using CoupledODETools, DiffEqGPU, OrdinaryDiffEq, CuArrays

CuArrays.allowscalar(false)

function scan(net::Network, tspan::Tuple{Float32,Float32};
    alg = Tsit5(), saveat = .1f0, output_func = (sol,i) -> (sol, false))
    n = net()
    u0 = n.u0s
    p = n.space[1]
    parr = [[n.space[i]...] for i in 1:length(n.space)]

    prob = ODEProblem(eval(n.fs), n.u0s, tspan, p)
    prob_func = (prob,i,repeat) -> remake(prob,p=parr[i]);
    monteprob = EnsembleProblem(prob, prob_func = prob_func)
    sol = solve(monteprob,alg,EnsembleGPUArray(),trajectories=length(parr),saveat=saveat)
    resol = reshape(sol.u, size(n.space)...)
end

@Component loren begin
    Dx = σ*(y-x), .1
    Dy = x*(ρ -g*k -z) - y, 0
    Dz = x*y - β*z, 0
    ρout = k -> ~x
end begin
    σ = range(9f0, 11f0, length = 10) |> collect
    β = range(2f0,4f0, length = 10) |> collect
    ρ = 28#range(20f0,35f0, length = 3) |> collect
    g = 0
end
l = loren(:l)
net1 = Network([l])
net1()

sol = scan(net1, (0f0, 100f0))

l1 = loren(:l1, ρout = :l2)
l2 = loren(:l2, ρout = :l1)
net = Network([l1, l2])

sol = scan(net, (0f0, 100f0))

n =net()

n.space


text:
using Gtk, Plots, Test, VisualRegressionTests

# popup on local machine only
istravis = get(ENV, "CI", "false") == "true"

# load an image for testing
import FileIO, TestImages
img = TestImages.testimage("cameraman")

# save a temporary copy
reffn = tempname() * ".png"
FileIO.save(reffn, img)

# this is the test function, which saves an image to the given location
func = fn -> FileIO.save(fn, img)

@testset "VisualRegressionTests.jl" begin
  @testset "Basic" begin
    result = test_images(VisualTest(func, reffn))
    @test isa(result, VisualTestResult)
    @test result.status == EXACT_MATCH
    @test result.err == nothing
    @test result.diff == 0
  end

  @testset "Macros" begin
    @visualtest func "VisualTest.png" !istravis

    @plottest plot([1.,2.,3.]) "PlotTest.png" !istravis

    @plottest begin
      plot([1.,2.,3.])
      plot!([3.,2.,1.])
    end "MorePlotTest.png" !istravis

    plotit() = plot([1.,2.,3.])
    @plottest plotit "FuncPlotTest.png" !istravis
  end
end


text:
########1D functions for benchmarking###########
################################################
normal = (posterior=Normal(0,1),mode=0,mean=0,var=1,chi2=[9999],ks=[9999.],ahmi=[9999.])
exp1 = (posterior= Exponential(1),mode=0,mean=1,var=1,chi2=[9999],ks=[9999.],ahmi=[9999.])
cauchy = (posterior=Cauchy(0,1),mode=0,mean=0,var=Inf,chi2=[9999],ks=[9999.],ahmi=[9999.])
#cauchy_variance_for_edges = x -> 0.6365969042714238*x+-0.9834816540023621

funnel = (posterior=BAT.FunnelDistribution(1.2,0.,1),mode=0,mean=0,var=2.07350,chi2=[9999],ks=[9999.],ahmi=[9999.])
multi_cauchy = (posterior=BAT.MultimodalCauchy(µ=5.,σ=2.,n=1),mode=[-5,5],mean=0,var=Inf,chi2=[9999],ks=[9999.],ahmi=[9999.])

testfunctions_1D = Dict(
	"normal" 		=> normal,
	"exp"			=> exp1,
	"cauchy"		=> cauchy,
	#"multi_cauchy"	=> multi_cauchy,
	"funnel"		=> funnel
)

analytical_stats_name=["mode","mean","var"]	#could be taken into NamedTuple for easier addtions but would be needed to implmented into calcs anyway
sample_stats=[Vector{Float64}(undef,length(analytical_stats_name)) for i in 1:length(testfunctions_1D)]

run_stats_names = ["nsteps","nchains","Times"]
run_stats=[Vector{Float64}(undef,length(run_stats_names)) for i in 1:length(testfunctions_1D)]


text:
function readmathematica(filename)
    datastr = readline(filename)
    datastr = replace(datastr, "},{"=>"; ")
    datastr = replace(datastr, "}}"=>"]")
    datastr = replace(datastr, "{{"=>"[")
    datastr = replace(datastr, "*^"=>"e")
    datastr = replace(datastr, ","=>" ")
    dataMat = eval(Meta.parse(datastr))
    return dataMat[:,1], dataMat[:,2]
end

function getFilename(L,BETA,DISGAM,sigmas,nREALS,SPACING; index="",singlet=false,format="mathematica")

    gamString = rpad(DISGAM,4,"0")
    nRealsPrime = maximum(sigmas) == 0 ? 0 : nREALS
    initString = singlet ? "s" : "up"
    sigString = string("σJ",rpad(sigmas[1],4,"0"),"_","σγ",rpad(sigmas[2],4,"0"),"_","στ",rpad(sigmas[3],4,"0"),"_",lpad(nRealsPrime,5,"0"))

    return joinpath(pwd(),string("jdata",index),string(format,"_",L,"_",initString,"_β",rpad(BETA,4,"0"),"_γ",gamString,"_",sigString,"_",SPACING))
end

function saveFidelities(L,BETA,DISGAM,sigmas,nREALS,SPACING; index="",singlet=false,format="mathematica", data=nothing)

    filename = getFilename(L,BETA,DISGAM,sigmas,nREALS,SPACING; index=index, singlet=singlet,format="mathematica")
	mkpath(joinpath(pwd(),string("jdata",index)))

    if isfile(filename) && parse(Int,split(strip(read(`wc -c $filename`, String))," ")[1]) > 0
        println("File already found. Skipping.")
    else
        if isnothing(data)
            println("No file found. Calculating...")
            data = calculateFidelities(L,BETA,DISGAM,sigmas,nREALS,SPACING;singlet)
            println("Done.")
        end
        theseExponents = data[1]
        f = open(filename,"w")
        if format=="mathematica"
            # Mathematica plotting format
            print(f, "{")
            for i in 1:length(theseExponents)
                print(f,"{",data[1][i],",",replace(string(data[2][i]),"e"=>"*^"),"}")
                i == length(theseExponents) ? print(f,"}") : print(f,",")
            end
        elseif format == "julia"
            # Julia plotting format
            write(f, data[1], data[2])
        end
        close(f)
    end
    nothing
end

function averager(L,BETA,DISGAM,sigmas,nREALS,SPACING; singlet=false)
    #Currently only supports mathematica formats
    numExps = Int(3/SPACING + 1)
    avg = zeros(numExps)
    jSWAPs = zeros(numExps)

    n = 0
    for i in 1:100
        filename= getFilename(L, BETA, DISGAM, sigmas, nREALS, SPACING; index=i, singlet=singlet, format="mathematica")
        if isfile(filename)
            jSWAPs = readmathematica(filename)[1]
            avg += readmathematica(filename)[2]
            n += 1
        else
            continue
        end
    end

    thisData = jSWAPs, avg ./ n
    saveFidelities(L, BETA, DISGAM, sigmas, nREALS*n , SPACING; index="AVG", singlet=singlet, format="mathematica", data=thisData)
    return

end

function plotter!(L,BETA,DISGAM,sigmas,nREALS,SPACING;singlet=false,format="mathematica")

    n = length(range(0.0,3.0,step=SPACING))
    emptyarray = zeros(Float64,n,2)

    filename = getFilename(L,BETA,DISGAM,sigmas,nREALS,SPACING; singlet=singlet, format=format)

    if !isfile(filename)
        println("File is missing.")
        return nothing
    elseif format=="mathematica"
        plotdata = readmathematica(filename)
    else
        f = open(filename,"r")
        data = read!(f,emptyarray)
        close(f)
        plotdata = (data[:,1],data[:,2])
    end

    plot!(plotdata, scale=:log10, size=(500,500))
    return nothing
end


text:
function fill_bag!(cur_bag::AbstractBag)

  cur_bag_chars = Random.shuffle(Tetris.tetrominoes)

  for cur_char in cur_bag_chars
    make_piece!(cur_bag, cur_char)
  end

end


text:
#=
Frodo

Run in this directory with
    julia --project=.. frodo.jl Baggins \
        --auto-hide \
        --friends "Samwise Gangee" "Peregrin Took" "Meriadoc Brandybuck" "Fredegar Bolger"
=#

using ArgParse2

in_danger() = rand(Bool)

function julia_main()::Cint
    parser = ArgumentParser(description = "Welcome to Middle Earth",
                            epilog = "There is no real going back")

    add_argument!(parser, "surname", help = "Your surname")
    add_argument!(parser, "-s", "--ring-size", type = Int, help = "Ring size")
    add_argument!(parser,
        "--auto-hide",
        action = "store_true",
        default = false,
        help = "Turn invisible when needed")
    add_argument!(parser, "--friends", metavar="FRIEND", nargs = "+", required = true)

    args = parse_args(parser)

    println("Welcome, Frodo $(args.surname)!")
    if args.ring_size === nothing
        println("Let's get you fitted for a ring.  We'll need to measure your ring size.")
    else
        println("I see your ring size is $(args.ring_size)")
    end

    if in_danger() && args.auto_hide
        println("Orcs are near!  You turn invisible.")
    end

    if length(args.friends) > 0
        travel_companions = join(args.friends, ", ", " and ")
        println("You're traveling to Mordor with $travel_companions.")
    end

    return 0
end

julia_main()


text:
function call(::Types.Join, args::Vector, storage, conn::IO)
    ip_addr = string(first(getpeername(conn)))
    name = args[1]
    arr = Requests.get_channel(storage.sql_conn, name)

    if (length(arr) == 0)
        return (write(conn, "This channel doesnt exist\n"))
    end

    if (arr.protected == 1)
        if (length(args) != 2)
            return (write(conn, "This channel required a password\n"))
        end
        if (Bcrypt.CompareHashAndPassword(arr.password, String(args[2])) == false)
            return (write(conn, "Invalid password\n"))
        end
    end

    storage.active_clients[ip_addr].current_channel_id = string(arr.id)
    if (Network.channel_exist(storage, string(arr.id)) == false)
        storage.active_channels[string(arr.id)] = Types.Channel(string(arr.id),
        arr.name, arr.description, arr.protected, arr.password, string(arr.owner))
    end
end


text:
# load_usual_suspects.jl

@time using DataFrames, ExcelReaders
@time using Gadfly
@time using Makie, GeometryTypes, Colors

println("Done!\n")

max_stars_repo_name: anthony-wang/julia-scripts

max_issues_repo_name: anthony-wang/julia-scripts

max_forks_repo_name: anthony-wang/julia-scripts