text:
import CSV
import DataFrames
import JSON

# Kept for compatibility with legacy code.
function create_profile_data(number_of_periods, data, loadprofile = ones(length(data["load"]),number_of_periods), genprofile = ones(length(data["gen"]),number_of_periods))
    _FP.make_time_series(data, number_of_periods; loadprofile = permutedims(loadprofile), genprofile = permutedims(genprofile))
end

function create_profile_data_italy!(data)

    hours = _FP.dim_length(data, :hour)
    scenarios = _FP.dim_length(data, :scenario)

    genprofile = ones(length(data["gen"]), hours*scenarios)
    loadprofile = ones(length(data["load"]), hours*scenarios)

    monte_carlo = get(_FP.dim_meta(data, :scenario), "mc", false)

    for (s, scnr) in _FP.dim_prop(data, :scenario)
        pv_sicily, pv_south_central, wind_sicily = read_res_data(s; mc = monte_carlo)
        demand_center_north_pu, demand_north_pu, demand_center_south_pu, demand_south_pu, demand_sardinia_pu = read_demand_data(s; mc = monte_carlo)

        start_idx = (s-1) * hours
        if monte_carlo == false
            for h in 1 : hours
                h_idx = scnr["start"] + ((h-1) * 3600000)
                genprofile[3, start_idx + h] = pv_south_central["data"]["$h_idx"]["electricity"]
                genprofile[5, start_idx + h] = pv_sicily["data"]["$h_idx"]["electricity"]
                genprofile[6, start_idx + h] = wind_sicily["data"]["$h_idx"]["electricity"]
            end
        else
            genprofile[3, start_idx + 1 : start_idx + hours] = pv_south_central[1: hours]
            genprofile[5, start_idx + 1 : start_idx + hours] = pv_sicily[1: hours]
            genprofile[6, start_idx + 1 : start_idx + hours] = wind_sicily[1: hours]
        end
        loadprofile[:, start_idx + 1 : start_idx + hours] = [demand_center_north_pu'; demand_north_pu'; demand_center_south_pu'; demand_south_pu'; demand_sardinia_pu'][:, 1: hours]
        # loadprofile[:, start_idx + 1 : start_idx + number_of_hours] = repeat([demand_center_north_pu'; demand_north_pu'; demand_center_south_pu'; demand_south_pu'; demand_sardinia_pu'][:, 1],1,number_of_hours)
    end
    # Add bus locations to data dictionary
    data["bus"]["1"]["lat"] = 43.4894; data["bus"]["1"]["lon"] = 11.7946; # Italy central north
    data["bus"]["2"]["lat"] = 45.3411; data["bus"]["2"]["lon"] =  9.9489; # Italy north
    data["bus"]["3"]["lat"] = 41.8218; data["bus"]["3"]["lon"] = 13.8302; # Italy central south
    data["bus"]["4"]["lat"] = 40.5228; data["bus"]["4"]["lon"] = 16.2155; # Italy south
    data["bus"]["5"]["lat"] = 40.1717; data["bus"]["5"]["lon"] =  9.0738; # Sardinia
    data["bus"]["6"]["lat"] = 37.4844; data["bus"]["6"]["lon"] = 14.1568; # Sicily
    # Return info
    return data, loadprofile, genprofile
end

function create_profile_data_germany!(data)

    hours = _FP.dim_length(data, :hour)
    scenarios = _FP.dim_length(data, :scenario)

    genprofile = ones(length(data["gen"]), hours*scenarios)
    loadprofile = ones(length(data["load"]), hours*scenarios)

    monte_carlo = get(_FP.dim_meta(data, :scenario), "mc", false)

    for (s, scnr) in _FP.dim_prop(data, :scenario)
        wind_profile = read_res_data(s; mc = monte_carlo, country = "de")
        demand_profile = read_demand_data(s; mc = monte_carlo, country = "de")
        start_idx = (s-1) * hours
        if monte_carlo == false
            for h in 1 : hours
                h_idx = scnr["start"] + ((h-1) * 3600000)
                genprofile[2, start_idx + h] = wind_profile["2"]["data"]["$h_idx"]["electricity"]
                genprofile[4, start_idx + h] = wind_profile["5"]["data"]["$h_idx"]["electricity"]
                genprofile[20, start_idx + h] = wind_profile["67"]["data"]["$h_idx"]["electricity"]
                if length(data["gen"]) > 20
                    genprofile[21, start_idx + h] = wind_profile["23"]["data"]["$h_idx"]["electricity"]
                elseif length(data["gen"]) > 21
                    genprofile[22, start_idx + h] = wind_profile["54"]["data"]["$h_idx"]["electricity"]
                end
            end
        end
        loadprofile[:, start_idx + 1 : start_idx + hours] .= repeat(demand_profile[1:hours]', size(loadprofile, 1))
    end
    # Return info
    return data, loadprofile, genprofile
end

"Create load and generation profiles for CIGRE distribution network."
function create_profile_data_cigre(data, number_of_hours; start_period = 1, scale_load = 1.0, scale_gen = 1.0, file_profiles_pu = normpath(@__DIR__,"..","data","cigre_mv_eu","time_series","CIGRE_profiles_per_unit.csv"))

    ## Fixed parameters

    file_load_ind = normpath(@__DIR__,"..","data","cigre_mv_eu","CIGRE_industrial_loads.csv")
    file_load_res = normpath(@__DIR__,"..","data","cigre_mv_eu","CIGRE_residential_loads.csv")
    scale_unit    = 0.001 # scale factor from CSV power data to FlexPlan power base: here converts from kVA to MVA

    ## Import data

    load_ind    = CSV.read(file_load_ind, DataFrames.DataFrame)
    load_res    = CSV.read(file_load_res, DataFrames.DataFrame)
    profiles_pu = CSV.read(
        file_profiles_pu,
        DataFrames.DataFrame;
        skipto = start_period + 1, # +1 is for header line
        limit = number_of_hours,
        ntasks = 1 # To ensure exact row limit is applied
    )
    if DataFrames.nrow(profiles_pu) < number_of_hours
        Memento.error(_LOGGER, "insufficient number of rows in file \"$file_profiles_pu\" ($number_of_hours requested, $(DataFrames.nrow(profiles_pu)) found)")
    end
    DataFrames.select!(profiles_pu,
        :industrial_load  => :load_ind,
        :residential_load => :load_res,
        :photovoltaic     => :pv,
        :wind_turbine     => :wind,
        :fuel_cell        => :fuel_cell,
        :CHP_diesel       => :chp_diesel,
        :CHP_fuel_cell    => :chp_fuel_cell
    )
    profiles_pu = Dict(pairs(eachcol(profiles_pu)))

    ## Prepare output structure

    extradata = Dict{String,Any}()

    ## Loads

    # Compute active and reactive power base of industrial loads
    DataFrames.rename!(load_ind, [:bus, :s, :cosϕ])
    load_ind.p_ind = scale_load * scale_unit * load_ind.s .* load_ind.cosϕ
    load_ind.q_ind = scale_load * scale_unit * load_ind.s .* sin.(acos.(load_ind.cosϕ))
    DataFrames.select!(load_ind, :bus, :p_ind, :q_ind)

    # Compute active and reactive power base of residential loads
    DataFrames.rename!(load_res, [:bus, :s, :cosϕ])
    load_res.p_res = scale_load * scale_unit * load_res.s .* load_res.cosϕ
    load_res.q_res = scale_load * scale_unit * load_res.s .* sin.(acos.(load_res.cosϕ))
    DataFrames.select!(load_res, :bus, :p_res, :q_res)

    # Create a table of industrial and residential power bases, indexed by the load ids used by `data`
    load_base = coalesce.(DataFrames.outerjoin(load_ind, load_res; on=:bus), 0.0)
    load_base.bus = string.(load_base.bus)
    bus_load_lookup = Dict{String,String}()
    for (l, load) in data["load"]
        bus_load_lookup["$(load["load_bus"])"] = l
    end
    DataFrames.transform!(load_base, :bus => DataFrames.ByRow(b -> bus_load_lookup[b]) => :load_id)

    # Compute active and reactive power profiles of each load
    extradata["load"] = Dict{String,Any}()
    for l in eachrow(load_base)
        extradata["load"][l.load_id] = Dict{String,Any}()
        extradata["load"][l.load_id]["pd"] = l.p_ind .* profiles_pu[:load_ind] .+ l.p_res .* profiles_pu[:load_res]
        extradata["load"][l.load_id]["qd"] = l.q_ind .* profiles_pu[:load_ind] .+ l.q_res .* profiles_pu[:load_res]
    end

    ## Generators

    # Define a Dict for the technology of generators, indexed by the gen ids used by `data`
    gen_tech = Dict{String,Symbol}()
    gen_tech["1"]  = :pv
    gen_tech["2"]  = :pv
    gen_tech["3"]  = :pv
    gen_tech["4"]  = :fuel_cell
    gen_tech["5"]  = :pv
    gen_tech["6"]  = :wind
    gen_tech["7"]  = :pv
    gen_tech["8"]  = :pv
    gen_tech["9"]  = :chp_diesel
    gen_tech["10"] = :chp_fuel_cell
    gen_tech["11"] = :pv
    gen_tech["12"] = :fuel_cell
    gen_tech["13"] = :pv

    # Compute active and reactive power profiles of each generator
    extradata["gen"]  = Dict{String,Any}()
    for (g, gen) in data["gen"]
        if haskey(gen_tech, g)
            extradata["gen"][g] = Dict{String,Any}()
            extradata["gen"][g]["pmax"] = scale_gen * gen["pmax"] .* profiles_pu[gen_tech[g]]
            extradata["gen"][g]["pmin"] = scale_gen * gen["pmin"] .* profiles_pu[gen_tech[g]]
            extradata["gen"][g]["qmax"] = scale_gen * gen["qmax"] .* ones(number_of_hours)
            extradata["gen"][g]["qmin"] = scale_gen * gen["qmin"] .* ones(number_of_hours)
        end
    end

    return extradata
end

function read_demand_data(year; mc = false, country = "it")
    if country == "it"
        if mc == false
            if year > 2
                error("Only 2 scenarios are supported")
            end
            y = year + 2017
            # Read demand CSV files
            demand_north = CSV.read(normpath(@__DIR__,"../../test/data/case6/time_series/mc_false/demand_north_$y.csv"),DataFrames.DataFrame)[:,3]
            demand_center_north = CSV.read(normpath(@__DIR__,"../../test/data/case6/time_series/mc_false/demand_center_north_$y.csv"),DataFrames.DataFrame)[:,3]
            demand_center_south = CSV.read(normpath(@__DIR__,"../../test/data/case6/time_series/mc_false/demand_center_south_$y.csv"),DataFrames.DataFrame)[:,3]
            demand_south = CSV.read(normpath(@__DIR__,"../../test/data/case6/time_series/mc_false/demand_south_$y.csv"),DataFrames.DataFrame)[:,3]
            demand_sardinia = CSV.read(normpath(@__DIR__,"../../test/data/case6/time_series/mc_false/demand_sardinia_$y.csv"),DataFrames.DataFrame)[:,3]

            # Convert demand_profile to pu of maxximum
            demand_north_pu = demand_north ./ maximum(demand_north)
            demand_center_north_pu = demand_center_north ./ maximum(demand_center_north)
            demand_south_pu = demand_south ./ maximum(demand_south)
            demand_center_south_pu = demand_center_south ./ maximum(demand_center_south)
            demand_sardinia_pu = demand_sardinia ./ maximum(demand_sardinia)
        else
            if year > 35
                error("Only 35 scenarios are supported")
            end
            y = year - 1
            demand_north_pu = CSV.read(normpath(@__DIR__,"../../test/data/case6/time_series/mc_true/case_6_demand_$y.csv"),DataFrames.DataFrame)[:,3]
            demand_center_north_pu = CSV.read(normpath(@__DIR__,"../../test/data/case6/time_series/mc_true/case_6_demand_$y.csv"),DataFrames.DataFrame)[:,2]
            demand_center_south_pu = CSV.read(normpath(@__DIR__,"../../test/data/case6/time_series/mc_true/case_6_demand_$y.csv"),DataFrames.DataFrame)[:,4]
            demand_south_pu = CSV.read(normpath(@__DIR__,"../../test/data/case6/time_series/mc_true/case_6_demand_$y.csv"),DataFrames.DataFrame)[:,5]
            demand_sardinia_pu = CSV.read(normpath(@__DIR__,"../../test/data/case6/time_series/mc_true/case_6_demand_$y.csv"),DataFrames.DataFrame)[:,6]
        end

        return demand_north_pu, demand_center_north_pu, demand_center_south_pu, demand_south_pu, demand_sardinia_pu
    elseif country == "de"
        if year > 3
            error("Only 3 scenarios are supported")
        end
        y = year + 2016
        demand = CSV.read(normpath(@__DIR__,"../../test/data/case67/time_series/demand$y.csv"),DataFrames.DataFrame)[:,3]
        demand_pu = demand ./ maximum(demand)
        return demand_pu[1:4:end]
    end
end

function read_res_data(year; mc = false, country = "it")
    if country == "it"
        if mc == false
            if year > 2
                error("Only 2 scenarios are supported")
            end
            y = year + 2017
            pv_sicily = Dict()
            open(normpath(@__DIR__,"../../test/data/case6/time_series/mc_false/pv_sicily_$y.json")) do f
                dicttxt = read(f, String)  # file information to string
                pv_sicily = JSON.parse(dicttxt)  # parse and transform data
            end

            pv_south_central = Dict()
            open(normpath(@__DIR__,"../../test/data/case6/time_series/mc_false/pv_south_central_$y.json")) do f
                dicttxt = read(f, String)  # file information to string
                pv_south_central = JSON.parse(dicttxt)  # parse and transform data
            end

            wind_sicily = Dict()
            open(normpath(@__DIR__,"../../test/data/case6/time_series/mc_false/wind_sicily_$y.json")) do f
                dicttxt = read(f, String)  # file information to string
                wind_sicily = JSON.parse(dicttxt)  # parse and transform data
            end
        else
            if year > 35
                error("Only 35 scenarios are supported")
            end
            y = year - 1
            pv_sicily = CSV.read(normpath(@__DIR__,"../../test/data/case6/time_series/mc_true/case_6_PV_$y.csv"),DataFrames.DataFrame)[:,7]
            pv_south_central = CSV.read(normpath(@__DIR__,"../../test/data/case6/time_series/mc_true/case_6_PV_$y.csv"),DataFrames.DataFrame)[:,4]
            wind_sicily = CSV.read(normpath(@__DIR__,"../../test/data/case6/time_series/mc_true/case_6_wind_$y.csv"),DataFrames.DataFrame)[:,7]
        end
        return pv_sicily, pv_south_central, wind_sicily
    elseif country == "de"
        if year > 3
            error("Only 3 scenarios are supported")
        end
        y = year + 2016
        wind_profile = Dict{String, Any}()
        open(normpath(@__DIR__,"../../test/data/case67/time_series/wind_bus2_$y.json")) do f
            dicttxt = read(f, String)  # file information to string
            wind_profile["2"] = JSON.parse(dicttxt)  # parse and transform data
        end
        open(normpath(@__DIR__,"../../test/data/case67/time_series/wind_bus5_$y.json")) do f
            dicttxt = read(f, String)  # file information to string
            wind_profile["5"]  = JSON.parse(dicttxt)  # parse and transform data
        end
        open(normpath(@__DIR__,"../../test/data/case67/time_series/wind_bus23_$y.json")) do f
            dicttxt = read(f, String)  # file information to string
            wind_profile["23"]  = JSON.parse(dicttxt)  # parse and transform data
        end
        open(normpath(@__DIR__,"../../test/data/case67/time_series/wind_bus54_$y.json")) do f
            dicttxt = read(f, String)  # file information to string
            wind_profile["54"]  = JSON.parse(dicttxt)  # parse and transform data
        end
        open(normpath(@__DIR__,"../../test/data/case67/time_series/wind_bus67_$y.json")) do f
            dicttxt = read(f, String)  # file information to string
            wind_profile["67"]  = JSON.parse(dicttxt)  # parse and transform data
        end
        return wind_profile
    end
end


text:
# ================ Plotting 2D =================== #
using PyPlot
fig = figure()
ax = fig["add_subplot"](111, projection="3d")
stress_y =  200.0

function vm_upper(a, c)
    vals = f(a[1], a[2], c)
    vm(vals[1], vals[2], 200)
end
vm(a,b) = sqrt(a^2 - a*b + b^2) - stress_y
f(m,c) = [600*cos(c) 600*sin(c)].*m
x_vals = []
max_iter = 100
y_vals = []
for i=0:0.1:(2*pi+0.3)
    wf(x) = f(x, i)
    t = 0.01
    step = 2
    merkki = -1
    s11, s22 = wf(t)
    ii = 0
    while (abs(vm(s11, s22)) > 1e-7) && ii < max_iter
        val = vm(s11, s22)
        if sign(val) != merkki
            merkki *= -1
            step *= -0.5
        end
        t += step
        s11, s22 = wf(t)
        ii += 1
    end
    push!(x_vals, s11)
    push!(y_vals, s22)
end

#plot(x_vals, y_vals)
#PyPlot.title("von Mises yield surface in 2D")
#PyPlot.xlabel("Eig Stress 1")
#PyPlot.ylabel("Eig Stress 2")
#grid()
#show()

# ================ Plotting 3D =================== #
n(θ, ϕ) = [sin(θ)*cos(ϕ)
           sin(θ)*sin(ϕ)
           cos(θ)]

m(θ, ϕ, χ) = [-sin(ϕ)*cos(χ)-cos(θ)*cos(ϕ)*sin(χ)
               cos(ϕ)*cos(χ)-cos(θ)*sin(ϕ)*sin(χ)
               sin(θ)*sin(χ)]

w = [sqrt(2/3) * 200 * m(54.735 * pi / 180, 45 * pi/180, x) for x=0:0.15:(2*pi+0.1)]
base_vec =  [1 1 1] / sqrt(3)

for i=-5:5
    tt = [w[x] + vec(base_vec) + 50 * i for x=1:length(w)]
    x = map(x->tt[x][1], collect(1:length(w)))
    y = map(x->tt[x][2], collect(1:length(w)))
    z = map(x->tt[x][3], collect(1:length(w)))
    plot3D(x, y, z, color="blue")
end

tt = [w[x] + vec(base_vec) + 50 * -5 for x=1:length(w)]
x_start = map(x->tt[x][1], collect(1:length(w)))[1:5:end]
y_start = map(x->tt[x][2], collect(1:length(w)))[1:5:end]
z_start = map(x->tt[x][3], collect(1:length(w)))[1:5:end]

tt = [w[x] + vec(base_vec) + 50 * 5 for x=1:length(w)]
x_end = map(x->tt[x][1], collect(1:length(w)))[1:5:end]
y_end = map(x->tt[x][2], collect(1:length(w)))[1:5:end]
z_end = map(x->tt[x][3], collect(1:length(w)))[1:5:end]

for i=1:length(x_start)
    x = [x_start[i], x_end[i]]
    y = [y_start[i], y_end[i]]
    z = [z_start[i], z_end[i]]
    plot3D(x, y, z, color="blue")
end

info("Calculation finished")

# plot the surface
xx = zeros(10, 10)
yy = zeros(10, 10)

for i=1:10
    for j=1:10
        xx[i, j] = (i - 5) * 100
        yy[i, j] = (j - 5) * 100
    end
end

# calculate corresponding z
z = zeros(10, 10)
for i=1:10
    for j=1:10
        z[i, j] = 1
    end
end

# plot the surface
plot_surface(xx, yy, z, color="blue")
plot(x_vals, y_vals, zeros(length(y_vals)), color="yellow")
axis("equal")
PyPlot.title("von Mises yield surface in 3D")
PyPlot.xlabel("Eig Stress 1")
PyPlot.ylabel("Eig Stress 2")
PyPlot.zlabel("Eig Stress 3")
for i=1:360
    ax["view_init"](30, i)
    PyPlot.draw()
    savefig("$(lpad(i, 3, 0))_plot.png")
end

ext: jl


text:
using JLD
using Plots
using Images
using LsqFit

plotlyjs()

d = load("save.jld")
zs, κs, σs = d["zs"], d["ks"], d["ss"]
nzs = imfilter(zs, Kernel.gaussian(0.7))
plot(κs, σs, nzs', st=:surface, zlims=(0, 1), xlabel="κ" , ylabel="σ", zlabel="r_inf", title="R∞(σ, κ)")
display(current())

k = nzs[2:end, :] - nzs[1:end-1, :]

plot(κs, σs, k', zlims=(0,0.5), st=:surface, xlabel="κ" , ylabel="σ", zlabel="r_inf", title="R∞(σ, κ)")
display(current())

j = findall(k .> 0.09)

ys = map(j) do i; κs[i[1]]; end
xs = map(j) do i; σs[i[2]]; end

@. model(x, m) = m*x
fit = curve_fit(model, xs, ys, [0.])

plot(0:0.5:2.5, model(0:0.5:2.5, fit.param), yaxis="κ", xaxis="σ", xlims=(0,2.5), ylims=(0, 2.5))
plot!(xs, ys, st=:scatter, legend=:bottomright)

display(current())


text:
abstract type AbstractOperationModel end

struct DefaultOpModel<:AbstractOperationModel end

mutable struct ModelReference
    transmission::Type{<:PM.AbstractPowerModel}
    devices::Dict{Symbol, DeviceModel}
    branches::Dict{Symbol, DeviceModel}
    services::Dict{Symbol, ServiceModel}
end

"""
    ModelReference(::Type{T}) where {T<:PM.AbstractPowerFormulation}

Creates a model reference of the Power Formulation, devices, branches, and services.

# Arguments
-`model::Type{T} where {T<:PM.AbstractPowerFormulation} = CopperPlatePowerModel`:
-`devices::Dict{Symbol, DeviceModel} = devices`: device dictionary
-`branches::Dict{Symbol, BranchModel} = branches`: branch dictionary
-`services::Dict{Symbol, ServiceModel} = services`: service dictionary

# Example
```julia
model_ref= ModelReference(CopperPlatePowerModel, devices, branches, services)
```
"""
function ModelReference(::Type{T}) where {T<:PM.AbstractPowerModel}

    return  ModelReference(T,
                           Dict{Symbol, DeviceModel}(),
                           Dict{Symbol, DeviceModel}(),
                           Dict{Symbol, ServiceModel}())

end

mutable struct OperationModel{M<:AbstractOperationModel}
    model_ref::ModelReference
    sys::PSY.System
    canonical::CanonicalModel
end

"""
    OperationModel(::Type{M},
    model_ref::ModelReference,
    sys::PSY.System;
    optimizer::Union{Nothing, JuMP.OptimizerFactory}=nothing,
    kwargs...) where {M<:AbstractOperationModel,
                      T<:PM.AbstractPowerFormulation}

This builds the optimization model and populates the operation model

# Arguments
-`::Type{M} where {M<:AbstractOperationModel, T<:PM.AbstractPowerFormulation} = TestOptModel`:
The abstract operation model type
-`model_ref::ModelReference = model_ref`: The model reference made up of transmission, devices,
                                          branches, and services.
-`sys::PSY.System = c_sys5`: the system created using Power Systems

# Output
-`op_model::OperationModel`: The operation model contains the model type, model, Power
Systems system, and optimization model.

# Example
```julia
model_ref= ModelReference(CopperPlatePowerModel, devices, branches, services)
OpModel = OperationModel(TestOptModel, model_ref, c_sys5_re; PTDF = PTDF5, optimizer = GLPK_optimizer)
```

# Accepted Key Words
-`verbose::Bool = true`: verbose default is true
-`PTDF::PTDF = PTDF`: Passes the PTDF matrix into the optimization model
-`optimizer::union{Nothing,JuMP.OptimizerFactory} = GLPK_optimizer`: The optimizer gets passed
into the optimization model the default is nothing.
-`sequential_runs::Bool = false`: tells the model to do sequential runs
-`initial_conditions::DICKDA = DICKDA()`: default of Dict{ICKey, Array{InitialCondition}}
-`parameters::Bool = false`: enable JuMP parameters
-`forecast::Bool = true`: if true, forecast collects the time steps in Power Systems,
if false it runs for one time step
-`initial_time::Dates.DateTime = PSY.get_forecasts_initial_time(sys)`: initial time of forecast
"""
function OperationModel(::Type{M},
                        model_ref::ModelReference,
                        sys::PSY.System;
                        optimizer::Union{Nothing, JuMP.OptimizerFactory}=nothing,
                        kwargs...) where {M<:AbstractOperationModel}

    op_model = OperationModel{M}(model_ref,
                          sys,
                          CanonicalModel(model_ref.transmission, sys, optimizer; kwargs...))

    build_op_model!(op_model; kwargs...)

    return  op_model

end
"""
    OperationModel(op_model::Type{M},
                    ::Type{T},
                    sys::PSY.System;
                    kwargs...) where {M<:AbstractOperationModel,
                                    T<:PM.AbstractPowerFormulation}

This uses the Abstract Power Formulation to build the model reference and
the optimization model and populates the operation model struct.

# Arguments
-`op_model::Type{M} = where {M<:AbstractOperationModel`: Defines the type of the operation model
-`::Type{T} where T<:PM.AbstractPowerFormulation`: The power formulation used for model ref & optimization model
-`sys::PSY.System = c_sys5`: the system created in Power Systems

# Output
-`op_model::OperationModel`: The operation model contains the model type, model, Power
Systems system, and optimization model.

# Example
```julia
model_ref= ModelReference(CopperPlatePowerModel, devices, branches, services)
OpModel = OperationModel(TestOptModel, model_ref, c_sys5_re; PTDF = PTDF5, optimizer = GLPK_optimizer)
```

# Accepted Key Words
-`verbose::Bool = true`: verbose default is true
-`PTDF::PTDF = PTDF`: Passes the PTDF matrix into the optimization model
-`optimizer::union{Nothing,JuMP.OptimizerFactory} = GLPK_optimizer`: The optimizer gets passed
into the optimization model the default is nothing.
-`sequential_runs::Bool = false`: tells the model to do sequential runs
-`initial_conditions::DICKDA = DICKDA()`: default of Dict{ICKey, Array{InitialCondition}}
-`parameters::Bool = false`: enable JuMP parameters
-`forecast::Bool = true`: if true, forecast collects the time steps in Power Systems,
if false it runs for one time step
-`initial_time::Dates.DateTime = PSY.get_forecasts_initial_time(sys)`: initial time of forecast

"""
function OperationModel(::Type{M},
                        ::Type{T},
                        sys::PSY.System;
                        kwargs...) where {M<:AbstractOperationModel,
                                          T<:PM.AbstractPowerModel}

    optimizer = get(kwargs, :optimizer, nothing)

    return OperationModel{M}(ModelReference(T),
                          sys,
                          CanonicalModel(T, sys, optimizer; kwargs...))

end
"""
    OperationModel(::Type{T},
                    sys::PSY.System;
                    kwargs...) where {M<:AbstractOperationModel,
                                    T<:PM.AbstractPowerFormulation}

This uses the Abstract Power Formulation to build the model reference and
the optimization model and populates the operation model struct.

***Note:*** the abstract operation model is set to the default operation model

# Arguments
-`op_model::Type{M} = where {M<:AbstractOperationModel`: Defines the type of the operation model
-`::Type{T} where T<:PM.AbstractPowerFormulation`: The power formulation used for model ref & optimization model
-`sys::PSY.System = c_sys5`: the system created in Power Systems

# Output
-`op_model::OperationModel`: The operation model contains the model type, model, Power
Systems system, and optimization model.

# Example
```julia
model_ref= ModelReference(CopperPlatePowerModel, devices, branches, services)
OpModel = OperationModel(TestOptModel, model_ref, c_sys5_re; PTDF = PTDF5, optimizer = GLPK_optimizer)
```

# Accepted Key Words
-`verbose::Bool = true`: verbose default is true
-`PTDF::PTDF = PTDF`: Passes the PTDF matrix into the optimization model
-`optimizer::union{Nothing,JuMP.OptimizerFactory} = GLPK_optimizer`: The optimizer gets passed
into the optimization model the default is nothing.
-`sequential_runs::Bool = false`: tells the model to do sequential runs
-`initial_conditions::DICKDA = DICKDA()`: default of Dict{ICKey, Array{InitialCondition}}
-`parameters::Bool = false`: enable JuMP parameters
-`forecast::Bool = true`: if true, forecast collects the time steps in Power Systems,
if false it runs for one time step
-`initial_time::Dates.DateTime = PSY.get_forecasts_initial_time(sys)`: initial time of forecast

"""
function OperationModel(::Type{T},
                        sys::PSY.System;
                        kwargs...) where {T<:PM.AbstractPowerModel}

    return OperationModel(DefaultOpModel,
                         T,
                         sys; kwargs...)

end

get_transmission_ref(op_model::OperationModel) = op_model.model_ref.transmission
get_devices_ref(op_model::OperationModel) = op_model.model_ref.devices
get_branches_ref(op_model::OperationModel) = op_model.model_ref.branches
get_services_ref(op_model::OperationModel) = op_model.model_ref.services
get_system(op_model::OperationModel) = op_model.sys

function set_transmission_ref!(op_model::OperationModel{M},
                               transmission::Type{T}; kwargs...) where {T<:PM.AbstractPowerModel,
                                                                        M<:AbstractOperationModel}

    # Reset the canonical
    op_model.model_ref.transmission = transmission
    op_model.canonical = CanonicalModel(transmission,
                                        op_model.sys,
                                        op_model.canonical.optimizer_factory; kwargs...)

    build_op_model!(op_model; kwargs...)

    return
end

function set_devices_ref!(op_model::OperationModel{M},
                          devices::Dict{Symbol, DeviceModel}; kwargs...) where M<:AbstractOperationModel

    # Reset the canonical
    op_model.model_ref.devices = devices
    op_model.canonical = CanonicalModel(op_model.model_ref.transmission,
                                        op_model.sys,
                                        op_model.canonical.optimizer_factory; kwargs...)

    build_op_model!(op_model; kwargs...)

    return
end

function set_branches_ref!(op_model::OperationModel{M},
                           branches::Dict{Symbol, DeviceModel}; kwargs...) where M<:AbstractOperationModel

    # Reset the canonical
    op_model.model_ref.branches = branches
    op_model.canonical = CanonicalModel(op_model.model_ref.transmission,
                                        op_model.sys,
                                        op_model.canonical.optimizer_factory; kwargs...)

    build_op_model!(op_model; kwargs...)

    return
end

function set_services_ref!(op_model::OperationModel{M},
                           services::Dict{Symbol, DeviceModel}; kwargs...) where M<:AbstractOperationModel

    # Reset the canonical
    op_model.model_ref.services = services
    op_model.canonical = CanonicalModel(op_model.model_ref.transmission,
                                        op_model.sys,
                                        op_model.canonical.optimizer_factory; kwargs...)

    build_op_model!(op_model; kwargs...)

    return
end

function set_device_model!(op_model::OperationModel{M},
                           name::Symbol,
                           device::DeviceModel{D, B}; kwargs...) where {D<:PSY.Injection,
                                                                        B<:AbstractDeviceFormulation,
                                                                        M<:AbstractOperationModel}

    if haskey(op_model.model_ref.devices, name)
        op_model.model_ref.devices[name] = device
        op_model.canonical = CanonicalModel(op_model.model_ref.transmission,
                                            op_model.sys,
                                            op_model.canonical.optimizer_factory; kwargs...)
        build_op_model!(op_model; kwargs...)
    else
        error("Device Model with name $(name) doesn't exist in the model")
    end

    return

end

function set_branch_model!(op_model::OperationModel{M},
                           name::Symbol,
                           branch::DeviceModel{D, B}; kwargs...) where {D<:PSY.Branch,
                                                                        B<:AbstractDeviceFormulation,
                                                                        M<:AbstractOperationModel}

    if haskey(op_model.model_ref.branches, name)
        op_model.model_ref.branches[name] = branch
        op_model.canonical = CanonicalModel(op_model.model_ref.transmission,
                                            op_model.sys,
                                            op_model.canonical.optimizer_factory; kwargs...)
        build_op_model!(op_model; kwargs...)
    else
        error("Branch Model with name $(name) doesn't exist in the model")
    end

    return

end

function set_services_model!(op_model::OperationModel{M},
                             name::Symbol,
                             service::DeviceModel; kwargs...) where M<:AbstractOperationModel

    if haskey(op_model.model_ref.devices, name)
        op_model.model_ref.services[name] = service
        op_model.canonical = CanonicalModel(op_model.model_ref.transmission,
                                            op_model.sys,
                                            op_model.canonical.optimizer_factory; kwargs...)
        build_op_model!(op_model; kwargs...)
    else
        error("Branch Model with name $(name) doesn't exist in the model")
    end

    return

end

function construct_device!(op_model::OperationModel,
                           name::Symbol,
                           device_model::DeviceModel;
                           kwargs...)

    if haskey(op_model.model_ref.devices, name)
        error("Device with model name $(name) already exists in the Opertaion Model")
    end

    op_model.model_ref.devices[name] = device_model

    construct_device!(op_model.canonical,
                     get_system(op_model),
                      device_model,
                      get_transmission_ref(op_model);
                      kwargs...)

    JuMP.@objective(op_model.canonical.JuMPmodel,
                    MOI.MIN_SENSE,
                    op_model.canonical.cost_function)

    return

end

function construct_network!(op_model::OperationModel,
                            system_formulation::Type{T};
                            kwargs...) where {T<:PM.AbstractPowerModel}

    construct_network!(op_model.canonical, get_system(op_model), T; kwargs...)

    return
end

function get_initial_conditions(op_model::OperationModel)
    return op_model.canonical.initial_conditions
end

function get_initial_conditions(op_model::OperationModel,
                                ic::InitialConditionQuantity,
                                device::PSY.Device)

    canonical = op_model.canonical
    key = ICKey(ic, device)

    return get_ini_cond(canonical, key)

end


text:
"""
The `Utils` has common utils functions and types used in other modules.

```julia
using MIToS.Utils
```
"""
module Utils

using HTTP
using CodecZlib
using LightXML
using NamedArrays

export  # GeneralUtils.jl
        All,
        get_n_words,
        hascoordinates,
        select_element,
        matrix2list, list2matrix,
        check_pdbcode,
        getarray,
        # Read.jl
        FileFormat,
        lineiterator,
        check_file, isnotemptyfile,
        download_file,
        # Write.jl
        Commandline,
        # Imported from Base (and exported for docs)
        read,
        write

include("GeneralUtils.jl")
include("Read.jl")
include("Write.jl")

@deprecate deleteitems!(vector::Vector, items) filter!(x -> x ∉ items, vector)

# COMMAND LINE: Scripts Module
# ============================

module Scripts

export  parse_commandline,
        runscript, run_single_script,
        script,
        set_parallel,
        open_output, close_output,
        readorparse

using Pkg
using Distributed
using ArgParse, CodecZlib

include("Scripts.jl")

end # Scripts

end # Utils


text:
module Checkpoints

export checkpoint, @checkpoint, setdepot, getdepot, clear

# stdlib
using Serialization

const DEPOT = Ref{String}()

"""
    getdepot() -> String

Return the current checkpoint depot.
"""
getdepot() = DEPOT[]

"""
    setdepot(dir)

Set the checkpoint depot to `dir`. Create `dir` if it exists.
"""
function setdepot(dir::String)
    _make(dir)
    DEPOT[] = dir
    return dir
end

_make(dir) = !isdir(dir) && mkpath(dir)
_resolve(file) = joinpath(DEPOT[], file)

"""
    clear()

Clear all checkpoint in the current depot. Depots can be set with [`setdepot`](@ref) and
queried with [`getdepot`](@ref)
"""
function clear()
    rm(DEPOT[], recursive = true, force = true)
    _make(DEPOT[])
end

"""
    clear(checkpoint::String)

Remove `checkpoint` from the current depot.
"""
clear(file)  = rm(_resolve(file))

function checkpoint(f, file::String)
    path = _resolve(file)
    # Return cached object
    if ispath(path)
        return deserialize(path)
    else
        # Invoke function, save results
        object = f()
        serialize(path, object)
        return object
    end
end

"""
    @checkpoint f(args...; kwargs...) checkpoint

Execute f(args...:kwargs...) and save the results to the file `checkpoint` in the current depot. The
next time this expression is invoked, the results from the previous run will be returned.
"""
macro checkpoint(fn, path)
    return :(checkpoint(() -> $(esc(fn)), $(esc(path))))
end

end # module


text:
function lex(text::Vector{Char})
    nchars = length(text)

    if nchars < 1
        error("No text available for lexing")
    end

    index = 1
    c = text[index]

    state = WHITESPACE

    buffer = IOBuffer()

    tokens = Array(Token, 0)

    while index <= nchars
        if state == WHITESPACE
            # Continue to end of token
            while isspace(c) && index <= nchars
                write(buffer, c)
                index += 1
                if index <= nchars
                    c = text[index]
                end
            end

            # Store token
            contents = bytestring(buffer)
            truncate(buffer, 0)
            token = Token(WHITESPACE, contents)
            push!(tokens, token)

            # Transition into the next state
            if c == '@'
                # Check if directive
                write(buffer, '@')
                index += 1
                if index <= nchars
                    c = text[index]
                end
                if c == '@'
                    state = DIRECTIVE
                    write(buffer, '@')
                    index += 1
                    if index <= nchars
                        c = text[index]
                    end
                elseif c == ' '
                    error("Invalid @ encountered")
                else
                    state = TEXT
                    write(buffer, c)
                    index += 1
                    if index <= nchars
                        c = text[index]
                    end
                end
            else
                state = TEXT
                write(buffer, c)
                index += 1
                if index <= nchars
                    c = text[index]
                end
            end
        elseif state == DIRECTIVE
            # Continue to end of token
            while isalpha(c) && index <= nchars
                write(buffer, c)
                index += 1
                if index <= nchars
                    c = text[index]
                end
            end

            # Store token
            contents = bytestring(buffer)
            truncate(buffer, 0)
            token = Token(DIRECTIVE, contents)
            push!(tokens, token)

            # Transition into the next state
            if isspace(c)
                state = WHITESPACE
            else
                state = TEXT
            end
        elseif state == TEXT
            # Continue to end of token
            while !isspace(c) && index <= nchars
                write(buffer, c)
                index += 1
                if index <= nchars
                    c = text[index]
                end
            end

            # Store token
            contents = bytestring(buffer)
            truncate(buffer, 0)
            token = Token(TEXT, contents)
            push!(tokens, token)

            # Transition into the next state
            state = WHITESPACE
        end
    end

    token = Token(state, bytestring(buffer))
    push!(tokens, token)

    return tokens
end


text:
struct fdawarp
    f::Array{Float64}
    time::Vector{Float64}
    fn::Array{Float64}
    qn::Array{Float64}
    q0::Array{Float64}
    fmean::Vector{Float64}
    mqn::Vector{Float64}
    gam::Array{Float64}
    orig_var::Float64
    amp_var::Float64
    phase_var::Float64
    cost::Vector{Float64}
    lambda::Float64
    method::String
    omethod::String
    gamI::Vector{Float64}
    rsamps::Bool
    fs::Array{Float64}
    gams::Array{Float64}
    ft::Array{Float64}
    qs::Array{Float64}
end

struct vfpca
    q_pca::Array{Float64}
    f_pca::Array{Float64}
    latent::Vector{Float64}
    coef::Array{Float64}
    U::Array{Float64}
    id::Int64
    mqn::Vector{Float64}
    time::Vector{Float64}
end

struct hfpca
    gam_pca::Array{Float64}
    psi_pca::Array{Float64}
    latent::Vector{Float64}
    U::Array{Float64}
    coef::Array{Float64}
    vec::Array{Float64}
    gam_mu::Vector{Float64}
end

struct func
    x::Vector{Float64}
    y::Vector{Float64}
end

struct basis
    x::Vector{Float64}
    matrix::Array{Float64}
end

struct mcmc_results
    f2_warped::Vector
    gamma::Vector
    g_coef::Array{Float64}
    psi::Vector
    sigma1::Vector
    accept::Vector
    logl::Vector
    gamma_mat::Array{Float64}
    gamma_stats::Array{Float64}
    xdist::Vector
    ydist::Vector
end


text:
using Test, Pkg

ENV["JULIA_LOG_LEVEL"] = "WARN"

for submodule in ["Utilities/VariableTemplates",
                  "Utilities/ParametersType",
                  "Utilities/RootSolvers",
                  "Common/PlanetParameters",
                  "Common/MoistThermodynamics",
                  "Atmos/Parameterizations/SurfaceFluxes",
                  "Atmos/Parameterizations/TurbulenceConvection",
                  "Atmos/Parameterizations/Microphysics",
                  "Mesh",
                  "DGmethods",
                  "ODESolvers",
                  "Arrays",
                  "LinearSolvers",
                 ]

  println("Starting tests for $submodule")
  t = @elapsed include(joinpath(submodule,"runtests.jl"))
  println("Completed tests for $submodule, $(round(Int, t)) seconds elapsed")
end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

function prepareValueReference(md::fmi3ModelDescription, vr::fmi3ValueReferenceFormat)
    tvr = typeof(vr)
    if tvr == Array{fmi3ValueReference,1}
        return vr
    elseif tvr == fmi3ValueReference
        return [vr]
    elseif tvr == String
        return [fmi3String2ValueReference(md, vr)]
    elseif tvr == Array{String,1}
        return fmi3String2ValueReference(md, vr)
    elseif tvr == Int64
        return [fmi3ValueReference(vr)]
    elseif tvr == Array{Int64,1}
        return fmi3ValueReference.(vr)
    elseif tvr == Nothing
        return Array{fmi3ValueReference,1}()
    end

    @assert false "prepareValueReference(...): Unknown value reference structure `$tvr`."
end
function prepareValueReference(fmu::FMU3, vr::fmi3ValueReferenceFormat)
    prepareValueReference(fmu.modelDescription, vr)
end
function prepareValueReference(comp::FMU3Component, vr::fmi3ValueReferenceFormat)
    prepareValueReference(comp.fmu.modelDescription, vr)
end

"""
Returns an array of ValueReferences coresponding to the variable names.
"""
function fmi3String2ValueReference(md::fmi3ModelDescription, names::Array{String})
    vr = Array{fmi3ValueReference}(undef,0)
    for name in names
        reference = fmi3String2ValueReference(md, name)
        if reference === nothing
            @warn "Value reference for variable '$name' not found, skipping."
        else
            push!(vr, reference)
        end
    end
    vr
end

"""
Returns the ValueReference coresponding to the variable name.
"""
function fmi3String2ValueReference(md::fmi3ModelDescription, name::String)
    reference = nothing
    if haskey(md.stringValueReferences, name)
        reference = md.stringValueReferences[name]
    else
        @warn "No variable named '$name' found."
    end
    reference
end

function fmi3String2ValueReference(fmu::FMU3, name::Union{String, Array{String}})
    fmi3String2ValueReference(fmu.modelDescription, name)
end

"""
Returns an array of variable names matching a fmi3ValueReference.
"""
function fmi3ValueReference2String(md::fmi3ModelDescription, reference::fmi3ValueReference)
    [k for (k,v) in md.stringValueReferences if v == reference]
end
function fmi3ValueReference2String(md::fmi3ModelDescription, reference::Int64)
    fmi3ValueReference2String(md, fmi3ValueReference(reference))
end

function fmi3ValueReference2String(fmu::FMU3, reference::Union{fmi3ValueReference, Int64})
    fmi3ValueReference2String(fmu.modelDescription, reference)
end


text:
f = open(ENV["HOME"] * "/sandbox/org-wal.sh/events/juliacon2016.org")
lines = readlines(f)
close(f)

counter = 1
for line in lines
    m = match(r"^\*\* ", line)

    if m !== nothing
        dirname =lowercase(replace(strip("$line", ['*',':',',','.']), " ", "_"))
        println("$dirname")
    end
    counter += 1
end


text:
#---------------------------------------------------------
# # [Plot transect/cruise data](@id cruiseplots)
#---------------------------------------------------------

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/3_cruiseplot.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/3_cruiseplot.ipynb)

#md # !!! note
#md #     All the AIBECS tutorials and how-to guides are available as Jupyter notebooks.
#md #     You can view them with [nbviewer](https://nbviewer.jupyter.org/)or execute them online with [binder](https://mybinder.org/) by clicking on the badges above!
#md #     (Note that binder can be slow to launch and its memory caps can be a problem when running.)

# This guide follows the [basic plotting guide](@ref plots) and its goal is to plot data related to oceanographic observations from cruise expeditions.
# In this walkthrough, we will
# [1. Create a ficitious cruise track](@ref fictitiouscruise) with profile data
# [2. Plot transects](@ref transects)
# [3. Other plots](@ref other-plots)

# As in the [basic plotting guide](@ref plots), throughout this guide we will use the OCIM2 grid and we will create a `dummy` modelled tracer.

using AIBECS, Plots
grd, _ = OCIM2.load()
fdummy(lat, lon, depth) = @. cosd(lat) * sind(lon) + sqrt(depth) / 30
dummy = fdummy(latlondepthvecs(grd)...)

#-----------------------------------------------
# ## [A fictitious cruise](@id fictitiouscruise)
#-----------------------------------------------

# For cruise data we use the [OceanographyCruises.jl](https://github.com/briochemc/OceanographyCruises.jl) package.

#md # !!! note
#md #     The purpose of the [OceanographyCruises.jl](https://github.com/briochemc/OceanographyCruises.jl) package is to provide a Julia interface for handling discrete data from cruises.
#md #     One goal is to use it for, e.g., GEOTRACES data.
#md #     However here we merely use it to create fictitious cruise data.

#nb # > **Note:**
#nb # > The purpose of the [OceanographyCruises.jl](https://github.com/briochemc/OceanographyCruises.jl) package is to provide a Julia interface for handling discrete data from cruises.
#nb # > One goal is to use it for, e.g., GEOTRACES data.
#nb # > However here we merely use it to create fictitious cruise data.

# Let us create a station `Sydney` at (34°S, 152°E) and a `ALOHA` station at (22.75°N, 158°W).

using OceanographyCruises
Sydney = Station(name="Sydney", lat=-30, lon=156)
ALOHA = Station(name="ALOHA", lat=22.75, lon=-158)

# Now let's create a range of 10 stations from `Sydney` to `ALOHA`.

Nstations = 10
stations = range(Sydney, ALOHA, length=Nstations, westmostlon=0)

# (`westmostlon=0` ensures that the longitudes are in (0,360) to match the OCIM2 grid we use here.)

# We can now construct a fictitious cruise track

ct = CruiseTrack(name="CruisyMcCruiseFace", stations=stations)

# and check the station locations by overlaying a plot of the cruise's track over a surface map of the dummy tracer

surfacemap(dummy, grd, color=:grays)
plotcruisetrack!(ct, markercolor=:red)

# Let's create a transect of data that is almost equal to the dummy.

# First, a function for creating random depths

function randomdepths(n, max)
    depths = cumsum(rand(n+1))
    return max * view(depths,1:n) / maximum(depths)
end
Nobs = rand(1:20, Nstations) # number of obs per station/profile
depths = [randomdepths(Nobs[i], 4000) for i in 1:Nstations]
obs = [[fdummy(st.lat, st.lon, d) .+ 0.1randn() for d in depths[i]] for (i,st) in enumerate(stations)]
profiles = [DepthProfile(station=st, depths=depths[i], values=obs[i]) for (i,st) in enumerate(stations)]

t = Transect(tracer="dummy", cruise=ct.name, profiles=profiles)

#-----------------------------------------------
# ## [Transects](@id transects)
#-----------------------------------------------

# ### Transects

# We can plot the modelled `dummy` data along the `ct` cruise track in the zonal directiion (along longitudes) with

plottransect(dummy, grd, ct=ct)

# If we want the observations transect on top of it

plotscattertransect!(t)

# If you have the [GEOTRACES.jl](https://github.com/briochemc/GEOTRACES.jl) package installed and the GEOTRACES data installed at the right location, you can instead plot real data with something like
#
# ```julia
# using GEOTRACES
# t = GEOTRACES.transects("Fe")[1] # the first Fe transect is GA02
# plotscattertransect(t) # use a bang (!) to overlay on top of your model plot
# ```
#
# However, this cannot be showcased online because GEOTRACES decided its data should "not be distributed to third parties".

#----------------------------------------------------
# ## [Other plots](@id other-plots)
#----------------------------------------------------

# Plots the modelled ratio of two modelled tracers at a given station.

dummy1 = cosd.(latvec(grd)) + sqrt.(depthvec(grd)) / 30
dummy2 = cosd.(2latvec(grd)) + 0.5sqrt.(depthvec(grd)) / 30
ratioatstation(dummy1, dummy2, grd, ALOHA, xlabel="dummy 1", ylabel="dummy 2")

# This can be useful to compare stoichiometric ratios at different stations.

ratioatstation!(dummy1, dummy2, grd, Sydney, marker=:square)


text:
@testset "Polytopes" begin
  @testset "Segment" begin
    @test paramdim(Segment) == 1

    s1 = Segment(P2(0,0), P2(1,0))
    s2 = Segment(P2(0.5,0.0), P2(2,0))
    @test s1 ∩ s2 == Segment(P2(0.5,0.0), P2(1,0))
    @test s2 ∩ s1 == Segment(P2(0.5,0.0), P2(1,0))

    s1 = Segment(P2(0,0), P2(1,0))
    s2 = Segment(P2(0,0), P2(0,1))
    @test s1 ∩ s2 == P2(0,0)
    @test s2 ∩ s1 == P2(0,0)

    s1 = Segment(P2(0,0), P2(1,0))
    s2 = Segment(P2(0,0), P2(-1,0))
    @test s1 ∩ s2 == P2(0,0)
    @test s2 ∩ s1 == P2(0,0)

    s1 = Segment(P2(0,0), P2(0,1))
    s2 = Segment(P2(0,0), P2(0,-1))
    @test s1 ∩ s2 == P2(0,0)
    @test s2 ∩ s1 == P2(0,0)

    s1 = Segment(P2(1,1), P2(1,2))
    s2 = Segment(P2(1,1), P2(1,0))
    @test s1 ∩ s2 == P2(1,1)
    @test s2 ∩ s1 == P2(1,1)

    s1 = Segment(P2(1,1), P2(2,1))
    s2 = Segment(P2(1,0), P2(3,0))
    @test s1 ∩ s2 === nothing
    @test s2 ∩ s1 === nothing

    s1 = Segment(P2(0.181429364026879, 0.546811355144474),
                 P2(0.38282226144778, 0.107781953228536))
    s2 = Segment(P2(0.412498700935005, 0.212081819871479),
                 P2(0.395936725690311, 0.252041094122474))
    @test s1 ∩ s2 === nothing
    @test s2 ∩ s1 === nothing

    s1 = Segment(P2(1,2), P2(1,0))
    s2 = Segment(P2(1,0), P2(1,1))
    @test s1 ∩ s2 == Segment(P2(1,0), P2(1,1))

    s1 = Segment(P2(0,0), P2(2,0))
    s2 = Segment(P2(-2,0), P2(-1,0))
    s3 = Segment(P2(-1,0), P2(-2,0))
    @test s1 ∩ s2 === s2 ∩ s1 === nothing
    @test s1 ∩ s3 === s3 ∩ s1 === nothing

    s1 = Segment(P2(-1,0), P2(0,0))
    s2 = Segment(P2(0,0), P2(2,0))
    @test s1 ∩ s2 == s2 ∩ s1 == P2(0,0)

    s1 = Segment(P2(-1,0), P2(1,0))
    s2 = Segment(P2(0,0), P2(3,0))
    @test s1 ∩ s2 == s2 ∩ s1 == Segment(P2(0,0), P2(1,0))

    s1 = Segment(P2(0,0), P2(1,0))
    s2 = Segment(P2(0,0), P2(2,0))
    @test s1 ∩ s2 == s2 ∩ s1 == Segment(P2(0,0), P2(1,0))

    s1 = Segment(P2(0,0), P2(3,0))
    s2 = Segment(P2(1,0), P2(2,0))
    @test s1 ∩ s2 == s2 ∩ s1 == s2

    s1 = Segment(P2(0,0), P2(2,0))
    s2 = Segment(P2(1,0), P2(2,0))
    @test s1 ∩ s2 == s2 ∩ s1 == s2

    s1 = Segment(P2(0,0), P2(2,0))
    s2 = Segment(P2(1,0), P2(3,0))
    @test s1 ∩ s2 == s2 ∩ s1 == Segment(P2(1,0), P2(2,0))

    s1 = Segment(P2(0,0), P2(2,0))
    s2 = Segment(P2(2,0), P2(3,0))
    @test s1 ∩ s2 == s2 ∩ s1 == P2(2,0)

    s1 = Segment(P2(0,0), P2(2,0))
    s2 = Segment(P2(3,0), P2(4,0))
    @test s1 ∩ s2 === s2 ∩ s1 === nothing

    s1 = Segment(P2(2,1), P2(1,2))
    s2 = Segment(P2(1,0), P2(1,1))
    @test s1 ∩ s2 === s2 ∩ s1 === nothing

    s1 = Segment(P2(1.5,1.5), P2(3.0,1.5))
    s2 = Segment(P2(3.0,1.0), P2(2.0,2.0))
    @test s1 ∩ s2 == s2 ∩ s1 == P2(2.5,1.5)

    s1 = Segment(P2(0.94495744, 0.53224397), P2(0.94798386, 0.5344541))
    s2 = Segment(P2(0.94798386, 0.5344541), P2(0.9472896, 0.5340202))
    @test s1 ∩ s2 == s2 ∩ s1 == P2(0.94798386, 0.5344541)

    s1 = Segment(P2(0.,0.), P2(1., 1.))
    @test s1(T(0.)) == P2(0., 0.)
    @test s1(T(1.)) == P2(1., 1.)
    @test_throws DomainError(T(1.2), "s(t) is not defined for t outside [0, 1].") s1(T(1.2))
    @test_throws DomainError(T(-0.5), "s(t) is not defined for t outside [0, 1].") s1(T(-0.5))
  end

  @testset "Triangles" begin
    @test paramdim(Triangle) == 2

    t = Triangle(P2(0,0), P2(1,0), P2(0,1))
    @test signarea(t) == T(0.5)
    @test area(t) == T(0.5)

    t = Triangle(P2(0,0), P2(0,1), P2(1,0))
    @test signarea(t) == T(-0.5)
    @test area(t) == T(0.5)

    t = Triangle(P2(0,0), P2(1,0), P2(1,1))
    for p in P2[(0,0),(1,0),(1,1),(0.5,0.),(1.,0.5),(0.5,0.5)]
      @test p ∈ t
    end
    for p in P2[(-1,0),(0,-1),(0.5,1.)]
      @test p ∉ t
    end
  end

  @testset "Quadrangles" begin
    @test paramdim(Quadrangle) == 2

    q = Quadrangle(P2(0,0), P2(1,0), P2(1,1), P2(0,1))
    @test area(q) == T(1)

    q = Quadrangle(P2(0,0), P2(1,0), P2(1.5,1.0), P2(0.5,1.0))
    @test area(q) == T(1)

    q = Quadrangle(P2(0,0), P2(1,0), P2(1.5,1.0), P2(0.5,1.0))
    for p in P2[(0,0),(1,0),(1.5,1.0),(0.5,1.0),(0.5,0.5)]
      @test p ∈ q
    end
    for p in P2[(0,1),(1.5,0.0)]
      @test p ∉ q
    end
  end

  @testset "Chains" begin
    # constructors
    c1 = Chain(P2[(1,1),(2,2),(1,1)])
    c2 = Chain(P2(1,1),P2(2,2),P2(1,1))
    c3 = Chain(CircularVector(P2[(1,1),(2,2)]))
    c4 = Chain(T.((1,1.)),T.((2.,2.)),T.((1.,1.)))
    @test c2 isa Chain{2,T,Vector{P2}}
    @test c1 == c2 == c3 == c4

    # segments
    c = Chain(P2[(1,1),(2,2),(3,3)])
    @test collect(segments(c)) == [Segment(P2(1,1),P2(2,2)),Segment(P2(2,2),P2(3,3))]
    c = Chain(P2[(1,1),(2,2),(3,3),(1,1)])
    @test collect(segments(c)) == [Segment(P2(1,1),P2(2,2)),Segment(P2(2,2),P2(3,3)),Segment(P2(3,3),P2(1,1))]

    # unique vertices
    c = Chain(P2[(1,1),(2,2),(2,2),(3,3)])
    @test unique(c) == Chain(P2[(1,1),(2,2),(3,3)])
    @test c == Chain(P2[(1,1),(2,2),(2,2),(3,3)])
    unique!(c)
    @test c == Chain(P2[(1,1),(2,2),(3,3)])

    # closing/opening chains
    c = Chain(P2[(1,1),(2,2),(3,3)])
    close!(c)
    @test c == Chain(P2[(1,1),(2,2),(3,3),(1,1)])
    open!(c)
    @test c == Chain(P2[(1,1),(2,2),(3,3)])

    # reversing chains
    c = Chain(P2[(1,1),(2,2),(3,3)])
    reverse!(c)
    @test c == Chain(P2[(3,3),(2,2),(1,1)])
    c = Chain(P2[(1,1),(2,2),(3,3)])
    @test reverse(c) == Chain(P2[(3,3),(2,2),(1,1)])

    # angles and inner angles
    c = Chain(P2[(0,0),(1,0),(1,1),(0,1),(0,0)])
    @test angles(c) ≈ [-π/2, -π/2, -π/2, -π/2]
    c = Chain(P2[(0,0),(1,0),(1,1),(0,1)])
    @test angles(c) ≈ [-π/2, -π/2]
    c = Chain(P2[(0,0),(1,0),(1,1),(2,1),(2,2),(1,2),(0,0)])
    @test angles(c) ≈ [-atan(2), -π/2, +π/2, -π/2, -π/2, -(π-atan(2))]
    @test innerangles(c) ≈ [atan(2), π/2, 3π/2, π/2, π/2, π-atan(2)]

    # reconstruct chain from vertices
    c1 = Chain(P2[(0,0),(1,0),(1,1),(0,1),(0,0)])
    c2 = Chain(vertices(c1))
    @test c1 == c2

    # centroid
    c = Chain(P2[(0,0),(1,0),(1,1),(0,1),(0,0)])
    @test centroid(c) == P2(0.5, 0.5)

    # views
    c = Chain(P2[(0,0),(1,0),(1,1),(0,1),(0,0)])
    @test view(c, 1:3) == Chain(P2[(0,0),(1,0),(1,1)])
    @test view(c, 4:6) == Chain(P2[(0,1),(0,0),(1,0)])
  end

  @testset "Hexahedrons" begin
    @test paramdim(Hexahedron) == 3
  end

  @testset "Pyramids" begin
    @test paramdim(Pyramid) == 3
  end

  @testset "Tetrahedrons" begin
    @test paramdim(Tetrahedron) == 3
    t = Tetrahedron(P3[(0,0,0),(1,0,0),(0,1,0),(0,0,1)])
    @test measure(t) == T(1/6)
  end

  @testset "PolyAreas" begin
    @test paramdim(PolyArea) == 2

    # test accessor methods
    poly = PolyArea(P2[(1,2),(2,3),(1,2)])
    @test vertices(poly) == CircularVector(P2[(1,2),(2,3)])
    poly = PolyArea(P2[(1,2),(2,3),(1,2)], [P2[(1.1, 2.54),(1.4,1.5),(1.1,2.54)]])
    @test vertices(poly) == CircularVector(P2[(1,2),(2,3),(1.1,2.54),(1.4,1.5)])

    # COMMAND USED TO GENERATE TEST FILES (VARY --seed = 1, 2, ..., 5)
    # rpg --cluster 30 --algo 2opt --format line --seed 1 --output poly1
    fnames = ["poly$i.line" for i in 1:5]
    polys1 = [readpoly(T, joinpath(datadir, fname)) for fname in fnames]
    for poly in polys1
      @test nvertices(first(chains(poly))) == 30
      @test !hasholes(poly)
      @test issimple(poly)
      for algo in [WindingOrientation(), TriangleOrientation()]
        @test orientation(poly, algo) == [:CCW]
      end
      @test unique(poly) == poly
    end

    # COMMAND USED TO GENERATE TEST FILES (VARY --seed = 1, 2, ..., 5)
    # rpg --cluster 30 --algo 2opt --format line --seed 1 --output smooth1 --smooth 2
    fnames = ["smooth$i.line" for i in 1:5]
    polys2 = [readpoly(T, joinpath(datadir, fname)) for fname in fnames]
    for poly in polys2
      @test nvertices(first(chains(poly))) == 120
      @test !hasholes(poly)
      @test issimple(poly)
      for algo in [WindingOrientation(), TriangleOrientation()]
        @test orientation(poly, algo) == [:CCW]
      end
      @test unique(poly) == poly
    end

    # COMMAND USED TO GENERATE TEST FILES (VARY --seed = 1, 2, ..., 5)
    # rpg --cluster 30 --algo 2opt --format line --seed 1 --output hole1 --holes 2
    fnames = ["hole$i.line" for i in 1:5]
    polys3 = [readpoly(T, joinpath(datadir, fname)) for fname in fnames]
    for poly in polys3
      rings = chains(poly)
      @test nvertices(first(rings)) < 30
      @test all(nvertices.(rings[2:end]) .< 18)
      @test hasholes(poly)
      @test !issimple(poly)
      for algo in [WindingOrientation(), TriangleOrientation()]
        orients = orientation(poly, algo)
        @test orients[1] == :CCW
        @test all(orients[2:end] .== :CW)
      end
      @test unique(poly) == poly
    end

    # test bridges
    for poly in [polys1; polys2; polys3]
      b  = bridge(poly)
      nb = nvertices(b)
      np = nvertices.(chains(poly))
      @test nb ≥ sum(np)
      # triangle orientation always works even
      # in the presence of self-intersections
      @test orientation(b, TriangleOrientation()) == :CCW
      # winding orientation is only suitable
      # for simple polygonal chains
      if issimple(b)
        @test orientation(b, WindingOrientation()) == :CCW
      end
    end

    # bridges between holes
    outer = P2[(0,0),(1,0),(1,1),(0,1),(0,0)]
    hole1 = P2[(0.2,0.2),(0.4,0.2),(0.4,0.4),(0.2,0.4),(0.2,0.2)]
    hole2 = P2[(0.6,0.2),(0.8,0.2),(0.8,0.4),(0.6,0.4),(0.6,0.2)]
    poly  = PolyArea(outer, [hole1, hole2])
    @test vertices(poly) == P2[(0,0),(1,0),(1,1),(0,1),
                               (0.2,0.2),(0.2,0.4),(0.4,0.4),(0.4,0.2),
                               (0.6,0.2),(0.6,0.4),(0.8,0.4),(0.8,0.2)]
    chain = bridge(poly)
    target = T[
      0.0  0.2  0.2  0.4  0.4  0.6  0.6  0.8  0.8  0.6  0.4  0.2  0.0  1.0  1.0  0.0
      0.0  0.2  0.4  0.4  0.2  0.2  0.4  0.4  0.2  0.2  0.2  0.2  0.0  0.0  1.0  1.0
    ]
    @test vertices(chain) == Point.(eachcol(target))

    # test uniqueness
    points = P2[(1,1),(2,2),(2,2),(3,3),(1,1)]
    poly   = PolyArea(points)
    unique!(poly)
    @test first(chains(poly)) == Chain(points)

    # unique and bridges
    poly = PolyArea(P2[(0,0),(1,0),(1,0),(1,1),(1,2),(0,2),(0,1),(0,1),(0,0)])
    chain = poly |> unique |> bridge
    @test chain == Chain(P2[(0,0),(1,0),(1,1),(1,2),(0,2),(0,1),(0,0)])

    # centroid
    poly = PolyArea(P2[(0,0),(1,0),(1,1),(0,1),(0,0)])
    @test centroid(poly) == P2(0.5, 0.5)

    # point in polygonal area
    outer = P2[(0,0),(1,0),(1,1),(0,1),(0,0)]
    hole1 = P2[(0.2,0.2),(0.4,0.2),(0.4,0.4),(0.2,0.4),(0.2,0.2)]
    hole2 = P2[(0.6,0.2),(0.8,0.2),(0.8,0.4),(0.6,0.4),(0.6,0.2)]
    poly  = PolyArea(outer, [hole1, hole2])
    @test all(p ∈ poly for p in outer)
    @test P2(0.5,0.5) ∈ poly
    @test P2(0.2,0.6) ∈ poly
    @test P2(1.5,0.5) ∉ poly
    @test P2(-0.5,0.5) ∉ poly
    @test P2(0.25,0.25) ∉ poly
    @test P2(0.75,0.25) ∉ poly
    @test P2(0.75,0.75) ∈ poly
  end
end


text:
using Distributions
using SpecialFunctions
using JLD2
using FileIO
using Distributed

#### Analysis of double max

function build_double_opt_for_files(files,κ)
    num=length(files)
    δ= load(files[1],"delta")
    dsteps=length(δ)
    meanencouters_of_file=zeros(num)
    doubleopt=zeros(num,dsteps)
    aopt=zeros(num,dsteps)
    bopt=zeros(num,dsteps)
    dissopt=zeros(num,dsteps)
    usesopt=zeros(num,dsteps)
    for jj in 1:num
        δ= load(files[jj],"delta")
        β = load(files[jj],"beta")
        meanencounters2,repsOFproduction,asteps,bsteps,dsteps,amax= load(files[jj],"parameters")
        utilityopt1load = load(files[jj],"max_utility")
        utilityopt= load(files[jj],"max_utility")
        alphaopt= load(files[jj],"otimal_alpha")
        diss= load(files[jj],"dissmax")
        uses= load(files[jj],"meory_used_at_the_min")

        for ii in 1:dsteps
            for kk in 1:bsteps
                utilityopt[kk,ii] = utilityopt1load[kk,ii] -meanencounters2 *κ*β[kk]
            end
        end
        doubleopt[jj,:],aopt[jj,:],bopt[jj,:],dissopt[jj,:],usesopt[jj,:]=give_double_max_from_max_a_all(utilityopt./meanencounters2,alphaopt,β,diss./meanencounters2,uses)

        meanencouters_of_file[jj]=meanencounters2
    end
    return meanencouters_of_file,doubleopt,aopt,bopt,dissopt,usesopt,δ
end

function build_double_opt_for_files_amax(files,κ)
    num=length(files)
    δ= load(files[1],"delta")
    dsteps=length(δ)
    meanencouters_of_file=zeros(num)
    amaxes=zeros(num)
    doubleopt=zeros(num,dsteps)
    aopt=zeros(num,dsteps)
    bopt=zeros(num,dsteps)
    dissopt=zeros(num,dsteps)
    usesopt=zeros(num,dsteps)
    for jj in 1:num
        δ= load(files[jj],"delta")
        β = load(files[jj],"beta")
        meanencounters2,repsOFproduction,asteps,bsteps,dsteps,amax= load(files[jj],"parameters")
        utilityopt1load = load(files[jj],"max_utility")
        utilityopt= load(files[jj],"max_utility")
        alphaopt= load(files[jj],"otimal_alpha")
        diss= load(files[jj],"dissmax")
        uses= load(files[jj],"meory_used_at_the_min")

        for ii in 1:dsteps
            for kk in 1:bsteps
                utilityopt[kk,ii] = utilityopt1load[kk,ii] -meanencounters2 *κ*β[kk]
            end
        end
        doubleopt[jj,:],aopt[jj,:],bopt[jj,:],dissopt[jj,:],usesopt[jj,:]=give_double_max_from_max_a_all(utilityopt./meanencounters2,alphaopt,β,diss./meanencounters2,uses)
        amaxes[jj]=amax
        meanencouters_of_file[jj]=meanencounters2
    end
    return meanencouters_of_file,doubleopt,aopt,bopt,dissopt,usesopt,δ,amaxes
end

function give_double_max_from_max_a_all(utilityopt,alphaopt,β,diss,uses)
    ds=length(utilityopt[1,:])
    doubleopt=zeros(ds)
    aopt=zeros(ds)
    bopt=zeros(ds)
    dissopt=zeros(ds)
    usesopt=zeros(ds)
    for ii in 1:ds
        argbopt=argmax(utilityopt[:,ii])
        doubleopt[ii]=utilityopt[argbopt,ii]
        dissopt[ii]=diss[argbopt,ii]
        usesopt[ii]=uses[argbopt,ii]
        aopt[ii]=alphaopt[argbopt,ii]
        bopt[ii]=β[argbopt]
    end
    return     doubleopt,aopt,bopt,dissopt,usesopt
end

####File Analysis

function analysis_of_opt_for_kappa(filename,k_steps)
    κ=range(0,0.2,length=k_steps)
    utilityopt1load = load(filename,"max_utility")
    alphaopt2load= load(filename,"otimal_alpha")
    dissmax2load= load(filename,"dissmax")
    usesload= load(filename,"meory_used_at_the_min")
    δ= load(filename,"delta")
    β = load(filename,"beta")
    freeload=utilityopt1load.+dissmax2load
    meanencounters2,reps,asteps,bsteps,dsteps,amax= load(filename,"parameters")
    doubleopt=zeros(k_steps,dsteps)
    aopt=zeros(k_steps,dsteps)
    bopt=zeros(k_steps,dsteps)
    dissopt=zeros(k_steps,dsteps)
    usesopt=zeros(k_steps,dsteps)
    for kk in 1:k_steps
        utilityopt1 = zeros(bsteps,dsteps)
       for ii in 1:dsteps
            for jj in 1:bsteps
                 utilityopt1[jj,ii] = utilityopt1load[jj,ii] -meanencounters2 *κ[kk]*β[jj]
            end
        end
        doubleopt[kk,:],aopt[kk,:],bopt[kk,:],dissopt[kk,:],usesopt[kk,:]=give_double_max_from_max_a_all(utilityopt1,alphaopt2load,β,dissmax2load./meanencounters2,usesload)

    end
    return doubleopt,aopt,bopt,dissopt,usesopt, δ
end

function put_data_in_region_and_order(x,y,xlim,ylim)
    points=length(x)
    newpoints=0
    (xdown,xup)=xlim
    (ydown,yup)=ylim
    for ii in 1:points
        xi=x[ii]
        yi=y[ii]
        if (xi<xup && xi>xdown && yi<yup && yi>ydown)
            newpoints +=1
        end
    end
    xregion=zeros(newpoints)
    yregion=zeros(newpoints)
    jj=1
    for ii in 1:points
        xi=x[ii]
        yi=y[ii]
        if (xi<xup && xi>xdown && yi<yup && yi>ydown)
            xregion[jj]=xi *1
            yregion[jj]=yi*1
            jj +=1
        end
    end
     xorderd=zeros(newpoints)
    yorderd=zeros(newpoints)
    top=maximum(x)+10
    for jj in 1:newpoints
        argument= argmin(xregion)
        xorderd[jj]=xregion[argument]*1
        yorderd[jj]=yregion[argument]*1
        xregion[argument]=top*2
    end
    return xorderd,yorderd
end

function build_diff_meanencouters_for_files(files,meanencountersmax,reps,pos_delta,κ)
    num=length(files)
    utlity=zeros(num,meanencountersmax)
    dissopt=zeros(num,meanencountersmax)
    δ= load(files[1],"delta")
    meanencounters=range(1,length=meanencountersmax,step=1)
    meanencouters_of_file=zeros(num)
    uses=zeros(num,meanencountersmax)
    for jj in 1:num
        δ= load(files[jj],"delta")
        β = load(files[jj],"beta")
        meanencounters2,repsOFproduction,asteps,bsteps,dsteps,amax= load(files[jj],"parameters")
        utilityopt1load = load(files[jj],"max_utility")
        utilityopt= load(files[jj],"max_utility")
        alphaopt= load(files[jj],"otimal_alpha")
        for ii in 1:dsteps
            for kk in 1:bsteps
                utilityopt[kk,ii] = utilityopt1load[kk,ii] -meanencounters2 *κ*β[kk]
            end
        end
        doubleopt,aopt,bopt=give_double_max_from_max_a(utilityopt./meanencounters2,alphaopt,β)
        println("I found the opt its at $(bopt[pos_delta]), $(δ[pos_delta]) and $(aopt[pos_delta])")
        dummi0,dummi1,dummi2,dummi3,dummi4=test_dif_meanencouters(bopt[pos_delta],δ[pos_delta],aopt[pos_delta],meanencountersmax,reps)
        utlity[jj,:] = dummi0./dummi4 .- bopt[pos_delta] * κ
        dissopt[jj,:] = dummi3./dummi4
        uses[jj,:]=dummi1
        meanencouters_of_file[jj]=meanencounters2
    end
    return meanencouters_of_file,utlity,dissopt,meanencounters,uses
end

function test_dif_meanencouters(β,δ,α,meanencountersmax,reps)
    utilty=zeros(meanencountersmax)
    uses=zeros(meanencountersmax)
    free=zeros(meanencountersmax)
    diss=zeros(meanencountersmax)
    Threads.@threads for ii in 1:meanencountersmax
        utilty[ii],uses[ii],free[ii],diss[ii]=livesycle(β,δ,α,ii,reps)
    end
    meanencounters=range(1,length=meanencountersmax,step=1)
    return utilty,uses,free,diss,meanencounters
end

function builed_transistions_b_scaled(files,k_steps)
    κ=range(0,0.01/4,length=k_steps)
    κ=κ.*4/sqrt(π)
    transision_threshold=0.4 # persent of change which is a transition
     meanencouters_of_file,doubleopt1,aopt1,bopt1,dissopt1,usesopt1,δ= build_double_opt_for_files(files,0)
    dsteps=length(δ)
    encounter_steps=length(meanencouters_of_file)
    trans_points=zeros(encounter_steps,k_steps)
    for kk in 1:k_steps
        meanencouters_of_file,doubleopt,aopt,bopt,dissopt,usesopt,δ= build_double_opt_for_files(files,κ[kk])
        for ee in 1:encounter_steps
            ddiss=0.0
            transdiss=1
            for dd in 1:dsteps-1

                bdiss=(-bopt[ee,dd+1]+bopt[ee,dd])/bopt[ee,dd+1]
                if ((bdiss>transision_threshold) && bopt[ee,dd+1]>1)

                    transdiss=dd
                end
            end
            trans_points[ee,kk]=δ[transdiss]

        end
    end
    return meanencouters_of_file,δ,κ,trans_points
end

function give_double_max_from_max_a(utilityopt,alphaopt,β)
    ds=length(utilityopt[1,:])
    doubleopt=zeros(ds)
    aopt=zeros(ds)
    bopt=zeros(ds)
    for ii in 1:ds
        argbopt=argmax(utilityopt[:,ii])
        doubleopt[ii]=utilityopt[argbopt,ii]
        aopt[ii]=alphaopt[argbopt,ii]
        bopt[ii]=β[argbopt]
    end
    return     doubleopt,aopt,bopt
end

max_stars_repo_name: StatPhysBio/ImmuneMemoryDM

max_issues_repo_name: StatPhysBio/ImmuneMemoryDM

max_forks_repo_name: StatPhysBio/ImmuneMemoryDM


text:
#=
ASCIIStr type

Copyright 2017-2018 Gandalf Software, Inc., Scott P. Jones,
and other contributors to the Julia language
Licensed under MIT License, see LICENSE.md
Based in part on code for ASCIIString that used to be in Julia
=#

## overload methods for efficiency ##

function _string(coll)
    n = 0
    for str in coll
        n += ncodeunits(str)
    end
    buf, out = _allocate(UInt8, n)
    for str in coll
        @preserve str begin
            len = ncodeunits(str)
            unsafe_copyto!(out, pointer(str), len)
            out += len
        end
    end
    buf
end

string(c::MaybeSub{<:Str{ASCIICSE}}...) = length(c) == 1 ? c[1] : Str(ASCIICSE, _string(c))

## transcoding to ASCII ##

function convert(::Type{<:Str{ASCIICSE}}, str::AbstractString)
    isempty(str) && return empty_ascii
    len = length(str)
    buf, pnt = _allocate(UInt8, len)
    @inbounds for ch in str
        is_ascii(ch) || strerror(StrErrors.INVALID_ASCII, pnt - pointer(buf) + 1, ch)
        set_codeunit!(pnt, ch%UInt8)
        pnt += 1
    end
    Str(ASCIICSE, buf)
end

convert(::Type{<:Str{ASCIICSE}}, str::Str{<:SubSet_CSEs}) = strerror(StrErrors.INVALID_ASCII)

convert(::Type{<:Str{ASCIICSE}}, str::Str{<:Union{Text1CSE,BinaryCSE,LatinCSE,UTF8CSE}}) =
    (is_empty(str) ? empty_ascii
     : is_ascii(str) ? Str(ASCIICSE, str.data) : strerror(StrErrors.INVALID_ASCII))

function convert(::Type{<:Str{ASCIICSE}}, dat::Vector{UInt8})
    is_empty(dat) && return empty_ascii
    is_ascii(dat) || strerror(StrErrors.INVALID_ASCII)
    siz = sizeof(dat)
    buf = _allocate(siz)
    @preserve dat unsafe_copyto!(pointer(buf), pointer(dat), siz)
    Str(ASCIICSE, buf)
end

function convert(::Type{<:Str{ASCIICSE}}, str::SubString{T}) where
    {T<:Union{String, Str{<:Union{Binary_CSEs,LatinCSE,UTF8CSE}}}}
    is_empty(str) && return empty_ascii
    is_ascii(str) || strerror(StrErrors.INVALID_ASCII)
    @preserve str begin
        len = ncodeunits(str)
        buf, out = _allocate(UInt8, len)
        unsafe_copyto!(out, pointer(str), len)
        Str(ASCIICSE, buf)
    end
end

convert(::Type{<:Str{ASCIICSE}}, str::String) =
    (isempty(str) ? empty_ascii
     : is_ascii(str) ? Str(ASCIICSE, str) : strerror(StrErrors.INVALID_ASCII))

# TODO: allow transform as the first argument to replace?

ascii(str::MaybeSub{<:Str}) = is_ascii(str) ? convert(ASCIIStr, str).data : ascii_err()
ascii(str::MaybeSub{<:Str{<:SubSet_CSEs}}) = ascii_err()
ascii(str::Str{T}) where {T<:ASCIICSE} = str.data

# This should really use a keyword argument, and allow for the following possibilities:
# 1) use default substitution character \u{1a} for ASCII/Latin1 (SUB) or \u{fffd} for Unicode
# 2) use a different substitution character
# 3) use a substition string
# 4) use a substition function, which gets passed the location of the start of the invalid sequence
#    and returns a replacement string, and the new location after the invalid sequence.
#    Note: for that to work, would need to use the same function in the check_string function
# 5) throw an exception, with enough information that one can determine where in the input
#    the invalid character was, and what it was

# Todo: optimize this!
function _convert_ascii(a, invlen, invdat)
    len = length(a)
    cnt = 0
    @inbounds for i = 1:len ; cnt += (a[i] >= 0x80) ; end
    # Note: this doesn't make a copy if only ascii characters,
    # so that changes later to the vector will affect the string!
    # that is only safe if the vector came from something else immutable
    cnt == 0 && return Str(ASCIICSE, a)
    v = _allocate(len + cnt*invlen)
    out = 1
    @inbounds for i = 1:len
        if (ch = a[i]) <= 0x7f
            v[out] = ch
            out += 1
        else
            unsafe_copyto!(v, out, invdat, 1, invlen)
            out += invlen
        end
    end
    Str(ASCIICSE, v)
end

convert(::Type{<:Str{ASCIICSE}}, a::Vector{UInt8}, invalids_as::String) =
    _convert_ascii(a, sizeof(invalids_as), ascii(invalids_as))

convert(::Type{<:Str{ASCIICSE}}, a::Vector{UInt8}, invalids_as::Str{ASCIICSE}) =
    _convert_ascii(a, sizeof(invalids_as), invalids_as.data)

convert(::Type{<:Str{ASCIICSE}}, a::Vector{UInt8}, invalids_as::AbstractString) =
    convert(ASCIIStr, a, ascii(invalids_as))


text:
module Checker

"""
Get fitness by using control and test columns
"""
function get_fitness(controlCols::Array{Array{Int,1},1}, testCols::Array{Array{Int,1},1})::Int
    fitness = 0
    for (control, test) in zip(controlCols, testCols)
        zeroLen = length(control) - length(test)
        if zeroLen > 0
            push!(test, zeros(Int64, zeroLen)...)
        elseif zeroLen < 0
            push!(control, zeros(Int64, abs(zeroLen))...)
        end
        fitness -= sum(abs.(control - test))
    end
    return fitness
end
end


text:
using Cthulhu, OceanTurb

Δt = 1
N = 128
m_fe = Diffusion.Model(N=N, stepper=:ForwardEuler)

@descend OceanTurb.calc_explicit_rhs!(m_fe.timestepper.rhs, m_fe.timestepper.eqn, m_fe.solution, m_fe)


text:
using LinearAlgebra

include("./utils.jl")
include("./LossFunctions.jl")
include("./Dims.jl")
include("./Layers/Network.jl")
include("./Optimize/train.jl")

export NetWork, Dense, Conv, Flatten, Dropout, MaxPool, add_layer!, epochs


text:
include("behavior_cloning.jl")
include("CRR.jl")
include("FisherBRC.jl")
include("PLAS.jl")
include("common.jl")


text:
include("split.jl")
include("apply.jl")
include("combine.jl")
include("modify.jl")
include("meta.jl")
include("timearray.jl")
include("readwrite.jl")
include("timeseriesrc.jl")
include("deprecated.jl")


text:
# A functor for testing
struct Multiplier{T}
    x::T
end
(m::Multiplier)(y) = m.x * y
function ChainRulesCore.rrule(m::Multiplier, y)
    Multiplier_pullback(z̄) = Tangent{typeof(m)}(; x=y * z̄), m.x * z̄
    return m(y), Multiplier_pullback
end

# NoRules - has no rules defined
struct NoRules; end

@testset "test_helpers.jl" begin
    @testset "Multiplier functor test-helper" begin
        test_rrule(Multiplier(4.0), 3.0)
    end
end


text:
using Clang
using assimp_win_jll

function Base.islowercase(s::String)::Bool
    mapfoldl(Base.islowercase, (x, y) -> x && y, filter(isletter, s))
end

const ASSIMP_INCLUDE = joinpath(dirname(assimp_win_jll.assimp_path), "..", "include") |> normpath
const ASSIMP_HEADERS = String[
    joinpath(root, file)
    for (root, dirs, files) in walkdir(ASSIMP_INCLUDE) for file in files
    if endswith(file, ".h") && islowercase(file)
]
println.(ASSIMP_HEADERS)

ctx = DefaultContext()
STDLIB = String["-I", raw"C:\Program Files\mingw-w64\x86_64-8.1.0-win32-seh-rt_v6-rev0\mingw64\x86_64-w64-mingw32\include"]
parse_headers!(
    ctx, ASSIMP_HEADERS,
    args=STDLIB, includes=vcat(LLVM_INCLUDE, ASSIMP_INCLUDE),
)

# settings
ctx.libname = "assimp"
ctx.options["is_function_strictly_typed"] = false
ctx.options["is_struct_mutable"] = false

# write output
api_file = joinpath(@__DIR__, "libassimp_win_api.jl")
api_stream = open(api_file, "w")

for trans_unit in ctx.trans_units
    root_cursor = getcursor(trans_unit)
    push!(ctx.cursor_stack, root_cursor)
    header = spelling(root_cursor)
    @info "wrapping header: $header ..."
    # loop over all of the child cursors and wrap them, if appropriate.
    ctx.children = children(root_cursor)
    for (i, child) in enumerate(ctx.children)
        child_name = name(child)
        child_header = filename(child)
        ctx.children_index = i
        # choose which cursor to wrap
        startswith(child_name, "__") && continue  # skip compiler definitions
        child_name in keys(ctx.common_buffer) && continue  # already wrapped
        child_header != header && continue  # skip if cursor filename is not in the headers to be wrapped

        wrap!(ctx, child)
    end
    isempty(ctx.api_buffer) && continue

    @info "writing $(api_file)"
    println(api_stream, "# Julia wrapper for header: $(basename(header))")
    println(api_stream, "# Automatically generated using Clang.jl\n")
    print_buffer(api_stream, ctx.api_buffer)
    empty!(ctx.api_buffer)  # clean up api_buffer for the next header
end
close(api_stream)

# write "common" definitions: types, typealiases, etc.
common_file = joinpath(@__DIR__, "libassimp_win_common.jl")
open(common_file, "w") do f
    println(f, "# Automatically generated using Clang.jl\n")
    print_buffer(f, dump_to_buffer(ctx.common_buffer))
end