text:
@doc raw"""
Splitting integrator for the solution of initial value problems
```math
\dot{q} (t) = v(t, q(t)) , \qquad q(t_{0}) = q_{0} ,
```
whose vector field ``v`` is given as a sum of vector fields
```math
v (t) = v_1 (t) + ... + v_r (t) .
```

`IntegratorSplitting` has two constructors:
```julia
IntegratorSplitting{DT,D}(solutions::Tuple, f::Vector{Int}, c::Vector, Δt)
IntegratorSplitting(equation::SODE, tableau::AbstractTableauSplitting, Δt)
```
In the first constructor, `DT` is the data type of the state vector and `D`
the dimension of the system. In the second constructor, this information
is extracted from the equation.
The tuple `solutions` contains functions implementing the flow (exact solution)
of the vector fields `v_i`. The vectors `f` and `c` define the actual splitting
method: `f` is a vector of indices of the flows in the split equation to be
solved and `c` is a vector of the same size `f` that contains the coefficients
for each splitting step, i.e., the resulting integrator has the form
```math
\varphi_{\tau} = \phi_{c[s] \tau}^{v_{f[s]}} \circ \dotsc \circ \phi_{c[2] \tau}^{v_{f[2]}} \circ \phi_{c[1] \tau}^{v_{f[1]}} .
```
In the second constructor, these vectors are constructed from the tableau and
the equation.

"""
struct IntegratorSplitting{DT, TT, D, S, QT <: Tuple} <: ODEIntegrator{DT,TT}
    q::QT
    f::NTuple{S,Int64}
    c::NTuple{S,TT}
    Δt::TT

    function IntegratorSplitting{DT,D}(solutions::solType, f::Vector{Int}, c::Vector{TT}, Δt::TT) where {DT, TT, D, solType <: Tuple}
        @assert length(f) == length(c)
        ft = Tuple(f)
        ct = Tuple(c)
        new{DT,TT,D,length(f),solType}(solutions, ft, ct, Δt)
    end
end

function IntegratorSplitting(equation::SODE{DT}, tableau::ST, Δt::TT) where {DT, TT, ST <: AbstractTableauSplitting{TT}}
    @assert hassolution(equation)
    IntegratorSplitting{DT, ndims(equation)}(get_solutions(equation), get_splitting_coefficients(length(equation.q), tableau)..., Δt)
end

@inline Base.ndims(::IntegratorSplitting{DT,TT,D}) where {DT,TT,D} = D
timestep(int::IntegratorSplitting) = int.Δt

function integrate_step!(int::IntegratorSplitting{DT,TT}, sol::AtomicSolutionODE{DT,TT}) where {DT,TT}
    local cᵢ::TT
    local tᵢ::TT

    # compute splitting steps
    for i in eachindex(int.f, int.c)
        if int.c[i] ≠ zero(TT)
            cᵢ = timestep(int) * int.c[i]
            tᵢ = sol.t̄ + cᵢ

            # reset atomic solution
            reset!(sol, cᵢ)

            # compute new solution
            int.q[int.f[i]](tᵢ, sol.q̄, sol.q, cᵢ)
        end
    end
end


text:
function add_forecast!(
                       component::T,
                       forecast::ForecastInternal,
                      ) where T <: InfrastructureSystemsType
    add_forecast!(_get_forecast_container(component), forecast)
    @debug "Added $forecast to $(typeof(component)) $(component.name) " *
           "num_forecasts=$(length(component._forecasts.data))."
end

"""
    remove_forecast_internal!(
                              ::Type{T},
                              component::InfrastructureSystemsType,
                              initial_time::Dates.DateTime,
                              label::AbstractString,
                             ) where T <: ForecastInternal

Removes the metadata for a forecast.
The caller must also remove the actual time series data.
"""
function remove_forecast_internal!(
                                   ::Type{T},
                                   component::InfrastructureSystemsType,
                                   initial_time::Dates.DateTime,
                                   label::AbstractString,
                                  ) where T <: ForecastInternal
    remove_forecast!(T, _get_forecast_container(component), initial_time, label)
    @debug "Removed forecast from $component:  $initial_time $label."
end

function clear_forecasts!(component::InfrastructureSystemsType)
    container = _get_forecast_container(component)
    if !isnothing(container)
        clear_forecasts!(container)
        @debug "Cleared forecasts in $component."
    end
end

"""
    get_forecast(
                 ::Type{T},
                 component::Component,
                 initial_time::Dates.DateTime,
                 label::AbstractString,
                ) where T <: Forecast

Return a forecast for the entire time series range stored for these parameters.
"""
function get_forecast(
                      ::Type{T},
                      component::InfrastructureSystemsType,
                      initial_time::Dates.DateTime,
                      label::AbstractString,
                     ) where T <: Forecast
    forecast_type = forecast_external_to_internal(T)
    forecast = get_forecast(forecast_type, component, initial_time, label)
    storage = _get_time_series_storage(component)
    ts = get_time_series(storage, get_time_series_uuid(forecast))
    return make_public_forecast(forecast, ts)
end

"""
    get_forecast(
                 ::Type{T},
                 component::InfrastructureSystemsType,
                 initial_time::Dates.DateTime,
                 label::AbstractString,
                 horizon::Int,
                ) where T <: Forecast

Return a forecast for a subset of the time series range stored for these parameters.
The range may span time series arrays as long as those timestamps are contiguous.
"""
function get_forecast(
                      ::Type{T},
                      component::InfrastructureSystemsType,
                      initial_time::Dates.DateTime,
                      label::AbstractString,
                      horizon::Int,
                     ) where T <: Forecast
    if !has_forecasts(component)
        throw(ArgumentError("no forecasts are stored in $component"))
    end

    first_forecast = iterate(iterate_forecasts(component))[1]
    resolution = get_resolution(first_forecast)
    sys_horizon = get_horizon(first_forecast)

    forecast = get_forecast(
        forecast_external_to_internal(T),
        component,
        initial_time,
        resolution,
        sys_horizon,
        label,
        horizon,
    )

    return forecast
end

function get_forecast(
                      ::Type{T},
                      component::InfrastructureSystemsType,
                      initial_time::Dates.DateTime,
                      label::AbstractString,
                     ) where T <: ForecastInternal
    return get_forecast(T, _get_forecast_container(component), initial_time, label)
end

function get_forecast(
                      ::Type{T},
                      component::InfrastructureSystemsType,
                      initial_time::Dates.DateTime,
                      sys_resolution::Dates.Period,
                      sys_horizon::Int,
                      label::AbstractString,
                      horizon::Int,
                     ) where T <: ForecastInternal
    forecast_type = forecast_internal_to_external(T)
    @debug "Requested forecast" get_name(component) forecast_type label initial_time horizon
    forecasts = Vector{forecast_type}()
    end_time = initial_time + sys_resolution * horizon
    initial_times = get_forecast_initial_times(T, _get_forecast_container(component), label)

    times_remaining = horizon
    found_start = false

    # This code concatenates ranges of contiguous forecasts.
    # Each initial_time represents one time series array that is stored.
    # Each array has a length equal to the system horizon.
    for it in initial_times
        len = 0
        if !found_start
            end_chunk = it + sys_resolution * sys_horizon
            if it <= initial_time && end_chunk > initial_time
                start_index = Int((initial_time - it) / sys_resolution) + 1
                found_start = true
            else
                # Keep looking for the start.
                continue
            end
            if end_chunk >= end_time
                end_index = sys_horizon - Int((end_chunk - end_time) / sys_resolution)
                len = end_index - start_index + 1
            else
                len = sys_horizon - start_index + 1
            end
        else
            start_index = 1
            len = times_remaining > sys_horizon ? sys_horizon : times_remaining
        end

        push!(forecasts, _make_forecast(T, component, start_index, len, it, label))
        times_remaining -= len
        if times_remaining == 0
            break
        end
    end

    if isempty(forecasts)
        throw(ArgumentError("did not find a forecast matching the requested parameters"))
    end

    @assert times_remaining == 0

    # Run the type-specificc constructor that concatenates forecasts.
    return forecast_type(forecasts)
end

function _make_forecast(
                        ::Type{T},
                        component::InfrastructureSystemsType,
                        start_index::Int,
                        len::Int,
                        initial_time::Dates.DateTime,
                        label::AbstractString
                       ) where T <: ForecastInternal
    forecast = get_forecast(T, _get_forecast_container(component), initial_time, label)
    ts = get_time_series(
        _get_time_series_storage(component),
        get_time_series_uuid(forecast);
        index=start_index,
        len=len,
    )
    return make_public_forecast(forecast, ts)
end

"""
    get_forecast_values(
                        ::Type{T},
                        component::InfrastructureSystemsType,
                        initial_time::Dates.DateTime,
                        label::AbstractString,
                       ) where T <: Forecast

Return a TimeSeries.TimeArray where the forecast data has been multiplied by the forecasted
component field.
"""
function get_forecast_values(
                             ::Type{T},
                             mod::Module,
                             component::InfrastructureSystemsType,
                             initial_time::Dates.DateTime,
                             label::AbstractString,
                            ) where T <: Forecast
    forecast = get_forecast(T, component, initial_time, label)
    return get_forecast_values(mod, component, forecast)
end

function get_forecast_values(
                             mod::Module,
                             component::InfrastructureSystemsType,
                             forecast::Forecast)
    scaling_factors = get_data(forecast)
    label = get_label(forecast)
    accessor_func = getfield(mod, Symbol(label))
    data = scaling_factors .* accessor_func(component)
    return data
end

function has_forecasts(component::InfrastructureSystemsType)
    container = _get_forecast_container(component)
    return !isnothing(container) && !isempty(container)
end

function get_forecast_initial_times(
                                    ::Type{T},
                                    component::InfrastructureSystemsType,
                                   ) where T <: Forecast
    if !has_forecasts(component)
        throw(ArgumentError("$(typeof(component)) does not have forecasts"))
    end
    return get_forecast_initial_times(forecast_external_to_internal(T),
                                      _get_forecast_container(component))
end

function get_forecast_initial_times(
                                    ::Type{T},
                                    component::InfrastructureSystemsType,
                                    label::AbstractString,
                                   ) where T <: Forecast
    if !has_forecasts(component)
        throw(ArgumentError("$(typeof(component)) does not have forecasts"))
    end
    return get_forecast_initial_times(forecast_external_to_internal(T),
                                      _get_forecast_container(component),
                                      label)
end

function get_forecast_initial_times!(
                                     initial_times::Set{Dates.DateTime},
                                     component::InfrastructureSystemsType,
                                    )
    if !has_forecasts(component)
        throw(ArgumentError("$(typeof(component)) does not have forecasts"))
    end

    get_forecast_initial_times!(initial_times, _get_forecast_container(component))
end

function get_forecast_initial_times(component::InfrastructureSystemsType)
    if !has_forecasts(component)
        throw(ArgumentError("$(typeof(component)) does not have forecasts"))
    end

    initial_times = Set{Dates.DateTime}()
    get_forecast_initial_times!(initial_times, component)

    return sort!(collect(initial_times))
end

"""
    generate_initial_times(
                           component::InfrastructureSystemsType,
                           interval::Dates.Period,
                           horizon::Int;
                           initial_time::Union{Nothing, Dates.DateTime}=nothing,
                          )

Generates all possible initial times for the stored forecasts. This should return the same
result regardless of whether the forecasts have been stored as one contiguous array or
chunks of contiguous arrays, such as one 365-day forecast vs 365 one-day forecasts.

Throws ArgumentError if there are no forecasts stored, interval is not a multiple of the
system's forecast resolution, or if the stored forecasts have overlapping timestamps.

# Arguments
- `component::InfrastructureSystemsType`: Component containing forecasts.
- `interval::Dates.Period`: Amount of time in between each initial time.
- `horizon::Int`: Length of each forecast array.
- `initial_time::Union{Nothing, Dates.DateTime}=nothing`: Start with this time. If nothing,
  use the first initial time.
"""
function generate_initial_times(
                                component::InfrastructureSystemsType,
                                interval::Dates.Period,
                                horizon::Int;
                                initial_time::Union{Nothing, Dates.DateTime}=nothing,
                               )
    # This throws if no forecasts.
    existing_initial_times = get_forecast_initial_times(component)

    first_forecast = iterate(iterate_forecasts(component))[1]
    resolution = Dates.Second(get_resolution(first_forecast))
    sys_horizon = get_horizon(first_forecast)

    first_initial_time, total_horizon = check_contiguous_forecasts(
        component, existing_initial_times, resolution, sys_horizon,
    )

    if isnothing(initial_time)
        initial_time = first_initial_time
    end

    interval = Dates.Second(interval)

    if interval % resolution != Dates.Second(0)
        throw(ArgumentError(
            "interval=$interval is not a multiple of resolution=$resolution"
        ))
    end

    last_initial_time = first_initial_time + total_horizon * resolution -
                        horizon * resolution
    initial_times = Vector{Dates.DateTime}()
    for it in range(initial_time, step=interval, stop=last_initial_time)
        push!(initial_times, it)
    end

    return initial_times
end

"""
    are_forecasts_contiguous(component::InfrastructureSystemsType)
"""
function are_forecasts_contiguous(component::InfrastructureSystemsType)
    existing_initial_times = get_forecast_initial_times(component)
    first_initial_time = existing_initial_times[1]

    first_forecast = iterate(iterate_forecasts(component))[1]
    resolution = Dates.Second(get_resolution(first_forecast))
    horizon = get_horizon(first_forecast)
    total_horizon = horizon * length(existing_initial_times)

    return _are_forecasts_contiguous(existing_initial_times, resolution, horizon)
end

function _are_forecasts_contiguous(initial_times, resolution, horizon)
    if length(initial_times) == 1
        return true
    end

    for i in range(2, stop=length(initial_times))
        if initial_times[i] != initial_times[i - 1] + resolution * horizon
            return false
        end
    end

    return true
end

"""
Throws ArgumentError if the forecasts are not in consecutive order.
"""
function check_contiguous_forecasts(
                                    component::InfrastructureSystemsType,
                                    existing_initial_times,
                                    resolution::Dates.Period,
                                    horizon::Int
                                   )
    if !_are_forecasts_contiguous(existing_initial_times, resolution, horizon)
        throw(ArgumentError(
            "generate_initial_times is not allowed with overlapping timestamps"
        ))
    end

    first_initial_time = existing_initial_times[1]
    total_horizon = horizon * length(existing_initial_times)
    return first_initial_time, total_horizon
end

function get_forecast_keys(component::InfrastructureSystemsType)
    return keys(_get_forecast_container(component).data)
end

function get_forecast_labels(
                             ::Type{T},
                             component::InfrastructureSystemsType,
                             initial_time::Dates.DateTime,
                            ) where T <: Forecast
    return get_forecast_labels(forecast_external_to_internal(T),
                               _get_forecast_container(component),
                               initial_time)
end

function get_num_forecasts(component::InfrastructureSystemsType)
    container = _get_forecast_container(component)
    if isnothing(container)
        return 0
    end

    return length(container.data)
end

function get_time_series(component::InfrastructureSystemsType, forecast::Forecast)
    storage = _get_time_series_storage(component)
    return get_time_series(storage, get_time_series_uuid(forecast))
end

function get_time_series_uuids(component::InfrastructureSystemsType)
    container = _get_forecast_container(component)

    return [(get_time_series_uuid(container.data[key]), key.label)
             for key in get_forecast_keys(component)]
end

"""
This function must be called when a component is removed from a system.
"""
function prepare_for_removal!(component::InfrastructureSystemsType)
    # Forecasts can only be part of a component when that component is part of a system.
    clear_time_series!(component)
    set_time_series_storage!(component, nothing)
    clear_forecasts!(component)
    @debug "cleared all forecast data from" component
end

function iterate_forecasts(component::InfrastructureSystemsType)
    container = _get_forecast_container(component)
    forecast_keys = sort!(collect(keys(container.data)), by = x -> x.initial_time)

    Channel() do channel
        for key in forecast_keys
            put!(channel, container.data[key])
        end
    end
end

function clear_time_series!(component::InfrastructureSystemsType)
    storage = _get_time_series_storage(component)
    if !isnothing(storage)
        for (uuid, label) in get_time_series_uuids(component)
            remove_time_series!(storage, uuid, get_uuid(component), label)
        end
    end
end

function set_time_series_storage!(
                                  component::InfrastructureSystemsType,
                                  storage::Union{Nothing, TimeSeriesStorage},
                                 )
    container = _get_forecast_container(component)
    if !isnothing(container)
        set_time_series_storage!(container, storage)
    end
end

function validate_forecast_consistency(component::InfrastructureSystemsType)
    # Initial times for each label must be identical.
    initial_times = Dict{String, Vector{Dates.DateTime}}()
    for key in keys(_get_forecast_container(component).data)
        if !haskey(initial_times, key.label)
            initial_times[key.label] = Vector{Dates.DateTime}()
        end
        push!(initial_times[key.label], key.initial_time)
    end

    if isempty(initial_times)
        return true
    end

    base_its = nothing
    for (label, its) in initial_times
        sort!(its)
        if isnothing(base_its)
            base_its = its
        elseif its != base_its
            @error "initial times don't match" base_its, its
            return false
        end
    end

    return true
end

function _get_forecast_container(component::InfrastructureSystemsType)
    # TODO: The get__forecasts methods in PowerSystems need to be IS.get__forecasts.
    #container = get__forecasts(component)
    if :_forecasts in fieldnames(typeof(component))
        container = component._forecasts
    else
        container = nothing
    end

    return container
end

function _get_time_series_storage(component::InfrastructureSystemsType)
    container = _get_forecast_container(component)
    if isnothing(container)
        return nothing
    end

    return container.time_series_storage
end

ext: jl


text:
mutable struct KrylovCtrl{T,verb,min_iter,max_iter,abstol,reltol,bad_num,bad_val,kill_NaN}
    err0::T
    cnt::Int
    bad::Int
    exitval::String
end

function get_tiny(::Type{KrylovCtrl{T,verb,min_iter,max_iter,abstol,reltol,bad_num,bad_val,kill_NaN}}) where {T,verb,min_iter,max_iter,abstol,reltol,bad_num,bad_val,kill_NaN}
    if T === Float32
        tiny = Float32(1e-32)
    elseif T === Float64
        tiny = 1e-60
    end
    return tiny
end

@inline get_tiny(c::KrylovCtrl) = get_tiny(typeof(c))

for s in (:verb,:min_iter,:max_iter,:abstol,:reltol,:bad_num,:bad_val,:kill_NaN)
    f = Symbol(:get_,s)
    @eval begin
        ($f)(::Type{KrylovCtrl{T,verb,min_iter,max_iter,abstol,reltol,bad_num,bad_val,kill_NaN}}) where {T,verb,min_iter,max_iter,abstol,reltol,bad_num,bad_val,kill_NaN} = ($s)

        @inline ($f)(c::KrylovCtrl) = ($f)(typeof(c))
    end
end

function file_to_dict(i::AbstractString)
    open(i) do f
        d = Dict{Symbol,Any}()
        for l in eachline(f)
            a = split(l)
            if length(a) == 2
                d[Symbol(replace(a[1],":",""))] = a[2]
            elseif length(a) > 2
                d[Symbol(replace(a[1],":",""))] = (a[2:end]...)
            end
        end
        return d
    end
end

function KrylovCtrl{T}(infilename::AbstractString) where {T<:Real}
    d = file_to_dict(infilename)
    verb = parse(Int,d[:verbosity])
    min_iter = Int(parse(d[:min_max_iter][1]))
    max_iter = Int(parse(d[:min_max_iter][2]))
    abstol = parse(T,d[:abstol])
    reltol = parse(T,d[:reltol])
    badnum = parse(Int,d[:badtol][1])
    badval = parse(T,d[:badtol][2])
    kill_NaN = parse(Bool,d[:kill_NaN])

    return KrylovCtrl{T,verb,min_iter,max_iter,abstol,reltol,badnum,badval,kill_NaN}(0.0,0,0,"null")
end

KrylovCtrl(infile::AbstractString) = KrylovCtrl{Float64}(infile)
KrylovCtrl() = KrylovCtrl("PCG.set")

function init(ctrl::KrylovCtrl,ierr)
    ctrl.err0 = abs(ierr)
    ctrl.cnt = 0
    ctrl.bad = 0
    ctrl.exitval = "null"
    return nothing
end

function test(ctrl::KrylovCtrl,rho,beta=0.)
    ctrl.cnt += 1
    cnt = ctrl.cnt

    cnt <= get_min_iter(ctrl) && return true

    if beta > get_bad_val(ctrl)
        ctrl.bad += 1
    else
        ctrl.bad = 0
    end

    ar = abs(rho)

    if cnt > get_max_iter(ctrl)
        ctrl.exitval = "max"
        return false
    elseif ar < get_abstol(ctrl)
        ctrl.exitval = "abs"
        return false
    elseif ar < (ctrl.err0 * get_reltol(ctrl))
        ctrl.exitval = "rel"
        return false
    elseif ar < float(cnt)*get_tiny(ctrl)
        ctrl.exitval = "!!! round"
        return false
    elseif ctrl.bad >= get_bad_num(ctrl)
        ctrl.exitval = "!!! bad"
        return false
    end

    return true
end


text:
## Exercise 8-9
## A string slice can take a third index. The first specifies the start, the third the end and the second the “step size”; that is, the number of spaces between successive characters. A step size of 2 means every other character; 3 means every third, etc.

# julia> fruit = "banana"
# "banana"
# julia> fruit[1:2:6]
# "bnn"

## A step size of -1 goes through the word backwards, so the slice [end:-1:1] generates a reversed string.

## Use this idiom to write a one-line version of ispalindrome from Exercise 6-6.
println("Ans: ")

function isPalindrome(string::AbstractString)
    string == string[end:-1:1]
end

print("Please enter a word: ")
word = strip(readline())
if isPalindrome(word)
    println("The word is a palindrome.")
else
    println("The word is not a palindrome")
end

println("End.")


text:
using BesselK, SpecialFunctions, ForwardDiff, StaticArrays, FiniteDifferences

# for quick testing/debugging:
const _oo = @SVector ones(2)
const _zz = @SVector zeros(2)
const _pp = @SVector [1.1, 2.1, 1.3]

# Convenience tool that will remove allocs:
@inline pv(scale, range, v) = @SVector [scale, range, v]

# FD things to compare with:
const FDFAST(fn, x) = (fn(x+1e-6) - fn(x))/1e-6
const FDACC  = central_fdm(10,1)
const FD2    = central_fdm(2,1)

# Note that this function uses the specialized method for (x^v)*besselk(v,x).
# Which can actually be MORE accurate than using AMOS or whatever for very small
# x. Pretty interesting. Even for pretty large order (say, v \approx 5), it is
# pretty accurate. More importantly, it speeds up the AD a bit and makes THAT
# more accurate.
#
# This version of the function, which the most accuracy-oriented users should
# consider as a default, uses AMOS for normal evaluations and then switches to
# our implementation for AD derivatives.
function matern(x, y, params)
  (sg, rho, nu) = params
  dist = norm(x-y)
  iszero(dist) && return sg*sg
  arg = sqrt(2*nu)*dist/rho
  (sg*sg*(2^(1-nu))/gamma(nu))*BesselK.adbesselkxv(nu, arg)
end

# A pure BesselK.jl-backed version of the Matern covariance.
function matern_us(x, y, params)
  (sg, rho, nu) = params
  dist = norm(x-y)
  iszero(dist) && return sg*sg
  arg = sqrt(2*nu)*dist/rho
  (sg*sg*(2^(1-nu))/gamma(nu))*BesselK._besselkxv(nu, arg)
end

# A pure AMOS-backed version of the Matern covariance. Consequently, you'll have
# to use finite difference derivatives of this one for order derivatives.
@inline beskxv(v, x) = BesselK.besselk(v, x)*(x^v)
function matern_amos(x, y, params)
  (sg, rho, nu) = params
  dist = norm(x-y)
  iszero(dist) && return sg*sg
  arg = sqrt(2*nu)*dist/rho
  (sg*sg*(2^(1-nu))/gamma(nu))*beskxv(nu, arg)
end

# Easily computed analytically, so skipping this.
matern_d1(x, y, params)    = matern(x, y, (sqrt(2*params[1]), params[2], params[3]))
matern_d1_d1(x, y, params) = matern(x, y, (sqrt(2), params[2], params[3]))

# Can be computed analytically, so we'll go ahead and assume a dedicated FD user
# has grinded out the analytical derivatives.
function matern_d2(x, y, params)
  ForwardDiff.derivative(p2->matern(x, y, pv(params[1], p2, params[3])), params[2])
end
function matern_d1_d2(x, y, params)
  ForwardDiff.derivative(p2->matern_d1(x, y, pv(params[1], p2, params[3])), params[2])
end
function matern_d2_d2(x, y, params)
  ForwardDiff.derivative(p2->matern_d2(x, y, pv(params[1], p2, params[3])), params[2])
end

# First derivatives:
for (stem, fn) in ((:fdfast, FDFAST),
                   (:fd2,    FD2),
                   (:ad,     ForwardDiff.derivative))
  name3 = Symbol(:matern_, stem, :_d3)
  if stem == :ad
    @eval $name3(x,y,p) = $fn(_p->matern(x,y,pv(p[1],p[2],_p)), p[3])
  else
    @eval $name3(x,y,p) = $fn(_p->matern_amos(x,y,pv(p[1],p[2],_p)), p[3])
  end
end

# Second derivatives:
for (stem, fn) in ((:fdfast, FDFAST),
                   (:ad,     ForwardDiff.derivative))
  # first names:
  name1 = Symbol(:matern_d1)
  name2 = Symbol(:matern_d2)
  name3 = Symbol(:matern_, stem, :_d3)
  # second names:
  name13 = Symbol(:matern_, stem, :_d1, :_d3)
  name22 = Symbol(:matern_, stem, :_d2, :_d2)
  name23 = Symbol(:matern_, stem, :_d2, :_d3)
  name33 = Symbol(:matern_, stem, :_d3, :_d3)
  # d1d3:
  @eval $name13(x,y,p) = $fn(_p->$name1(x,y,pv(p[1],p[2],_p)), p[3])
  # d2d3:
  @eval $name23(x,y,p) = $fn(_p->$name2(x,y,pv(p[1],p[2],_p)), p[3])
  # d3d3:
  @eval $name33(x,y,p) = $fn(_p->$name3(x,y,pv(p[1],p[2],_p)), p[3])
end

const AD_DFNS  = (matern_d1, matern_d2, matern_ad_d3)
const AD_D2FNS = ((matern_d1_d1, matern_d1_d2, matern_ad_d1_d3),
                  (matern_d2_d2, matern_ad_d2_d3),
                  (matern_ad_d3_d3,))

const FD_DFNS  = (matern_d1, matern_d2, matern_fdfast_d3)
const FD_D2FNS = ((matern_d1_d1, matern_d1_d2, matern_fdfast_d1_d3),
                  (matern_d2_d2, matern_fdfast_d2_d3),
                  (matern_fdfast_d3_d3,))


text:
# functions related to chi-square distribution

# R implementations
# For pdf and logpdf we use the Julia implementation
using .RFunctions:
    chisqcdf,
    chisqccdf,
    chisqlogcdf,
    chisqlogccdf,
    chisqinvcdf,
    chisqinvccdf,
    chisqinvlogcdf,
    chisqinvlogccdf

# Julia implementations
# promotion ensures that we do forward e.g. `chisqpdf(::Int, ::Float32)` to
# `gammapdf(::Float32, ::Int, ::Float32)` but not `gammapdf(::Float64, ::Int, ::Float32)`
chisqpdf(k::Real, x::Real) = chisqpdf(promote(k, x)...)
chisqpdf(k::T, x::T) where {T<:Real} = gammapdf(k / 2, 2, x)

chisqlogpdf(k::Real, x::Real) = chisqlogpdf(promote(k, x)...)
chisqlogpdf(k::T, x::T) where {T<:Real} = gammalogpdf(k / 2, 2, x)


text:
using QSM: V, W, F, Multigrid, cycle!
using QSM: BackwardSweep, ForwardSweep, SymmetricSweep, EarlyStopping, GaussSeidel
using QSM: restrict!, correct_prolong!
using QSM: residual!, tfill!

#####
##### TODO: Smoothing factor check
#####

@testset "Smoothers" begin
    include("smoothers.jl")
end

@testset "Transfer" begin
    include("transfer.jl")
end

@testset "Fixed point" for T in (Float32, Float64)
    # Cycling should not alter the exact solution.
    # Use exact solution as initial guess.
    for sz in [(128, 129, 130), (131, 132, 133)]
        A, u, d2u, d2uh = test_problem(T, sz)

        M = Multigrid(T, A,
            coarsesolver = GaussSeidel(SymmetricSweep(), iter=2500),
            presmoother = GaussSeidel(SymmetricSweep(), iter=1),
            postsmoother = GaussSeidel(SymmetricSweep(), iter=1),
        )

        for cycle in (V(), W(), F())
            x = M.workspace.x[1] = copy(u)
            M.workspace.b[1] = copy(d2uh)
            cycle!(M, cycle)
            @test x ≈ u
        end
    end
end

@testset "Two-level Cycle" for T in (Float32, Float64)
    # Test intergrid operators.
    # Postsmooth residuals should be lower than presmooth residuals.
    for sz in [(32, 33, 34), (35, 36, 37)]
        A, u, d2u, d2uh = test_problem(T, sz)

        M = Multigrid(T, A,
            coarsesolver = EarlyStopping(
                GaussSeidel(SymmetricSweep(), iter=10),
                sqrt(eps(T)),
                1_000_000
            ),
            presmoother = GaussSeidel(ForwardSweep(), iter=2),
            postsmoother = GaussSeidel(BackwardSweep(), iter=2),
            maxlevel = 2
        )

        nd0 = Inf
        nu0 = Inf
        M.workspace.x[1] = zeros(T, sz)
        M.workspace.b[1] = copy(d2uh)

        l = 1
        for i in 1:5
            A = M.grids[l].A
            Ac = M.grids[l+1].A

            presmooth! = M.grids[l].presmoother
            postsmooth! = M.grids[l].postsmoother

            x = M.workspace.x[l]
            b = M.workspace.b[l]
            r = M.workspace.r[l]

            xc = M.workspace.x[l+1]
            bc = M.workspace.b[l+1]

            presmooth!(x, A, b)
            residual!(r, b, A, x)
            restrict!(bc, Ac, r)

            xc = tfill!(xc, 0, Ac)
            M.coarsesolver(xc, Ac, bc)

            correct_prolong!(x, A, xc)
            postsmooth!(x, A, b)

            nd = L2norm(r)
            nu = L2norm(residual!(r, b, A, x))
            @test nd < nd0
            @test nu < nu0
            @test nu < nd
            nd0, nu0 = nd, nu
        end
    end
end


text:
using ProgressMeter
using Printf
using BSON: @save, @load
# using PyCall
# pygui(:qt)
# using PyPlot
# pygui(true)
using MATLAB
include("Dimer_Lagrangian_Descriptor.jl")

location="/media/brandon_behring/Extra_Space/BSON_FILES/"
file_name="BandF_Dimer_lagrangian_descriptors_data0_129250000000000_ 1000_2_2.250_4.5002019_12_01T12_53_22_675"

N=Int(1e5)
Energy=.129
t_end=5
Q_start=.2
Q_end=1.2

Q=range(Q_start, stop=Q_end, length=N)
mesh=[(q, 0.0) for q in Q]

LD=@showprogress map(Dimer_Lagrangian_Descriptor.Dimer_Lagrangian_Descriptor_Function, mesh, Energy*ones(N), t_end*ones(N))
dx=(Q_end-Q_start)/N
gradM=Dimer_Lagrangian_Descriptor.D_X(LD,dx)
gradMl=log10.(gradM)
# d2M=log10.(Dimer_Lagrangian_Descriptor.D_X(gradM,dx))
Q1=Q[4:N-3]
Q2=Q[7:N-6]

mat"figure();plot($Q, $LD,'k');hold on;"
mat"plot($Q1, $gradMl,'r')"
# mat"plot($Q2, $d2M,'b')"

# @load location*file_name*".bson" SAVE_DATA

# gradM=SAVE_DATA["gradM"]
# Energy=SAVE_DATA["Energy"]
# LD=SAVE_DATA["LD"]
# ArrP=SAVE_DATA["ArrP"]
# ArrQ=SAVE_DATA["ArrQ"]
# N=SAVE_DATA["N"]
# biggest=maximum(filter(!isnan,gradM))
# smallest=minimum(filter(!isnan,gradM))

# ArrP=ArrP[4:N-3]
# ArrQ=ArrQ[4:N-3]
# N=N-4
# # figure()
# # h_title = @sprintf("h= %.6f",Energy)
# # plot(ArrQ, LD[3,3:N-2], c=:blue)
# plot(ArrQ,gradM[1,:],c=:red)

# mat" imagesc($gradM) "
# grad=zeros(N,N)
# for i=1:N
#     for j=1:N
#        M=gradM[i,j]
#        if M<15.0
#             grad[i,j]=0.0
#        else
#             grad[i,j]=1.0
#        end
#     end
# end
# mat"figure"
# mat" imagesc($grad) "


text:
using Test
using Distributed
include(joinpath(Sys.BINDIR, "..", "share", "julia", "test", "testenv.jl"))
addprocs_with_testenv(3)
@test nworkers() == 3

@everywhere using Logging
for w in workers()
    # Do not log on worker nodes
    remotecall(()->global_logger(NullLogger()),w)
end

@everywhere using StochasticPrograms
using JuMP
using LinearAlgebra
using GLPKMathProgInterface

@everywhere import StochasticPrograms: probability, expected

include("../problems/problem_load.jl")
include("dfunctional_tests.jl")
include("dsolver_tests.jl")


text:
@testset "Heuristics" begin

# Load domains and problems
path = joinpath(dirname(pathof(ParallelPlanners)), "..", "domains", "gridworld")
gridworld = load_domain(joinpath(path, "domain.pddl"))
gw_problem = load_problem(joinpath(path, "problem-1.pddl"))
gw_state = init_state(gw_problem)

path = joinpath(dirname(pathof(ParallelPlanners)), "..", "domains", "doors-keys-gems")
doors_keys_gems = load_domain(joinpath(path, "domain.pddl"))
dkg_problem = load_problem(joinpath(path, "problem-1.pddl"))
dkg_state = init_state(dkg_problem)

path = joinpath(dirname(pathof(ParallelPlanners)), "..", "domains", "blocksworld")
blocksworld = load_domain(joinpath(path, "domain.pddl"))
bw_problem = load_problem(joinpath(path, "problem-0.pddl"))
bw_state = init_state(bw_problem)

@testset "Goal Count Heuristic" begin

goal_count = GoalCountHeuristic()
@test goal_count(gridworld, gw_state, gw_problem.goal) == 1
@test goal_count(doors_keys_gems, dkg_state, dkg_problem.goal) == 1
@test goal_count(blocksworld, bw_state, bw_problem.goal) == 2

end

@testset "Manhattan Heuristic" begin

manhattan = ManhattanHeuristic(@pddl("xpos", "ypos"))
@test manhattan(gridworld, gw_state, gw_problem.goal) == 2

end

@testset "HSP Heuristics" begin

@test HAdd()(blocksworld, bw_state, bw_problem.goal) == 4
@test HMax()(blocksworld, bw_state, bw_problem.goal) == 2

end

@testset "HSPr Heuristics" begin

bw_init = init_state(bw_problem)
bw_goal = goal_state(bw_problem)
h_add_r = precompute!(HAddR(), blocksworld, bw_init, bw_problem.goal)
h_max_r = precompute!(HMaxR(), blocksworld, bw_init, bw_problem.goal)

@test h_add_r(blocksworld, bw_goal, bw_problem.goal) == 4
@test h_max_r(blocksworld, bw_goal, bw_problem.goal) == 2

end

@testset "FF Heuristic" begin

ff = precompute!(FFHeuristic(), blocksworld, bw_state, bw_problem.goal)
@test ff(blocksworld, bw_state, bw_problem.goal) == 4

end

end


text:
"""
    A = JopDwt(T, n...[, wt=wavelet(WT.haar, WT.Lifting)])

`A` is the N-D (N=1,2 or 3) wavelet transform operator, operating on the domain
of type `T` and dimensions `n`.  The optional argument `wt` is the wavelet type.
The supported wavelet types are `WT.haar` and `WT.db` (Daubechies).

# Notes
* For 2 and 3 dimensional transforms the domain must be square.  In other-words, `size(dom,1)==size(dom,2)==size(dom,3)`.
* If your domain is rectangular, please consider using 1D wavelet transforms in combination wih the kronecker product (`JopKron`).
* The wavelet transform is provide by the Wavelets.jl package: `http://github.com/JuliaDSP/Wavelets.jl`
* You may try other wavelet classes supported by `Wavelets.jl`; however, these have yet to be tested for correctness with respect to the
transpose.

# Example
```julia
A = JopDwt(Float64, 64, 64; wt=wavelet(WT.db5))
d = A*rand(domain(A))
```
"""
function JopDwt(sp::JetSpace; wt=wavelet(WT.haar, WT.Lifting)) where {T,N}
    n = size(sp)
    @assert length(n) <= 3
    map(i->@assert(n[i]==n[1]), 2:length(n))
    @assert n[1] == nextpow(2, n[1])
    wt.name != "haar" && startswith(wt.name, "db") != true && error("JopDwt only supports Haar and Daubechies wavelet functions")
    JopLn(dom = sp, rng = sp, df! = JopDwt_df!, df′! = JopDwt_df′!, s = (wt=wt,))
end
JopDwt(::Type{T}, n::Vararg{N}; kwargs...) where {T,N} = JopDwt(JetSpace(T,n); kwargs...)
export JopDwt

function JopDwt_df!(d::AbstractArray, m::AbstractArray; wt, kwargs...)
    # dwt! does not work without copy for all wt
    # if ok with copy and gc
    copyto!(d, dwt(m, wt))
    #copyto!(d, m)
    #dwt!(d, op.wt, maxtransformlevels(m))
end

function JopDwt_df′!(m::AbstractArray, d::AbstractArray; wt, kwargs...)
    # idwt! does not work without copy for all wt
    # if ok with copy and gc
    copyto!(m, idwt(d, wt))
    #copyto!(m, d)
    #idwt!(m, op.wt, maxtransformlevels(d))
end


text:
using MolecularGraph:
    forward!, lgnot!, lghighand!, lgor!, lglowand!

@testset "smarts.logicaloperator" begin

@testset "logicaloperator" begin

    # Test function
    function parseabc(state)
        if read(state) == 'a'
            forward!(state)
            return :v => :a
        elseif read(state) == 'b'
            forward!(state)
            return :v => :b
        elseif read(state) == 'c'
            forward!(state)
            return :v => :c
        end
    end

    state = SmartsParser("!a", false)
    not1 = lgnot!(state, parseabc)
    @test not1 == (:not => (:v => :a))

    state = SmartsParser("a&b", false)
    and1 = lghighand!(state, parseabc)
    @test and1 == (:and => (:v => :a, :v => :b))

    state = SmartsParser("abc", false)
    and2 = lghighand!(state, parseabc)
    @test and2 == (:and => (:v => :a, :v => :b, :v => :c))

    state = SmartsParser("!ab", false)
    not2 = lghighand!(state, parseabc)
    @test not2 == (
        :and => (
            :not => (:v => :a),
            :v => :b
        )
    )

    state = SmartsParser("a,b", false)
    or1 = lgor!(state, parseabc)
    @test or1 == (:or => (:v => :a, :v => :b))

    state = SmartsParser("a;b", false)
    and3 = lglowand!(state, parseabc)
    @test and3 == (:and => (:v => :a, :v => :b))

    state = SmartsParser("", false)
    null = lglowand!(state, parseabc)
    @test null == nothing

    state = SmartsParser("!a!x", false)
    not3 = lglowand!(state, parseabc)
    @test not3 == (:not => (:v => :a))
    @test state.pos == 3

    state = SmartsParser("abcdef", false)
    and4 = lghighand!(state, parseabc)
    @test and4 == (:and => (:v => :a, :v => :b, :v => :c))
    @test state.pos == 4

    state = SmartsParser("a;b&c", false)
    and5 = lglowand!(state, parseabc)
    @test and5 == (
        :and => (
            :v => :a,
            :and => (:v => :b, :v => :c)
        )
    )

    state = SmartsParser("a,b&c", false)
    comp1 = lglowand!(state, parseabc)
    @test comp1 == (
        :or => (
            :v => :a,
            :and => (:v => :b, :v => :c)
        )
    )

    state = SmartsParser("a,b;c", false)
    comp2 = lglowand!(state, parseabc)
    @test comp2 == (
        :and => (
            :or => (:v => :a, :v => :b),
            :v => :c
        )
    )

    state = SmartsParser("ac;a!b,c;bx", false)
    comp3 = lglowand!(state, parseabc)
    @test comp3 == (
        :and => (
            :and => (:v => :a, :v => :c),
            :or => (
                :and => (
                    :v => :a,
                    :not => (:v => :b)
                ),
                :v => :c
            ),
            :v => :b
        )
    )
    @test state.pos == 11

end

end # logicalopelator


text:
# Tests convexity detection functions.
function gen_quad_expr(n::Int64 = 10)
    A = -5 .+ 10 .* rand(n, n)
    b = -5 .+ 10 .* rand(n)
    c = -5 + 10*rand()
    expr = :(x -> x .* $(A) .* x + sum(b .* x) + c)
    return expr
end

idxs = 1:25; # 25
# max_min = [26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38]
gms = Dict()
for idx in idxs
    try
        ex = sagemark_to_GlobalModel(idx, true);
        gms[idx] = ex
    catch
        println("Failed import of problem $(idx).")
    end
end

# Finding the bounded BBLs
bounded_bbls = []
bbl_idxs = Dict(idx => [] for idx in collect(idxs))
for idx in idxs
    gm = gms[idx]
    find_bounds!(gm)
    ct = 0
    for i = 1:length(gm.bbls)
        ct += 1
        try
            check_bounds(get_bounds(gm.bbls[i]))
            gm.bbls[i].name = gm.name * "_" * gm.bbls[i].name
            push!(bounded_bbls, gm.bbls[i])
            push!(bbl_idxs[idx], ct)
        catch ErrorException
        end
    end
end
bbl_idxs = Dict(key => val for (key, val) in bbl_idxs if !isempty(val))

@info "Found $(length(bounded_bbls)) bounded BBLs."

# Making sure to compute the convexity of the BBLs.
for bbl in bounded_bbls
    uniform_sample_and_eval!(bbl)
    update_vexity(bbl)
end
@info "Found $(sum(bbl.convex for bbl in bounded_bbls)) convex BBLs."
convex_bbls = [bbl for bbl in bounded_bbls if bbl.convex]

""" Checks the convexity of signomial constraints (sig >= 0). """
function check_cvx_con(alpha, c, lse::Bool = true)
    cvxity = true
    if lse # exponential form
        if length(c) == 1
        elseif length(c) >= 2
            if sum(c .> 0) != 1
                cvxity = false
            end
        end
    else # geometric form
        for i in 1:length(c)
            if c[i] >= 0
                if any(alpha[i,:] .> 0)
                    cvxity = false
                    break
                end
            else
                if any(alpha[i,:] .< 0)
                    cvxity = false
                    break
                end
            end
        end
    end
    return cvxity
end

""" Checks the convexity of signomial objectives. """
function check_cvx_obj(alpha, c, lse::Bool = true)
    cvxity = true
    if lse # exponential form
        if length(c) == 1
        elseif length(c) >= 2
            if any(c .< 0)
                cvxity = false
            end
        end
    else # geometric form
        for i in 1:length(c)
            if c[i] >= 0
                if any(alpha[i,:] .< 0)
                    cvxity = false
                    break
                end
            else
                if any(alpha[i,:] .> 0)
                    cvxity = false
                    break
                end
            end
        end
    end
    return cvxity
end

# Finding the actual convex ones analytically
lse = true
actual_cvx = Dict(idx => [] for idx in collect(idxs))

sagemarks = pyimport("sagebenchmarks.literature.solved");
for (idx, val) in bbl_idxs
    signomials, solver, run_fn = sagemarks.get_example(idx);
    f, greaters, equals = signomials;
    check_cvx_obj(f.alpha, f.c, lse) && push!(actual_cvx[idx], 1)
    for i = 1:length(greaters)
        # if to prune the bounds
        if !(sum(greaters[i].alpha .!= 0) == 1 && sum(greaters[i].alpha) == 1) && i+1 in val
            check_cvx_con(greaters[i].alpha, greaters[i].c, lse) &&
            push!(actual_cvx[idx], i+1)
        end
    end
end
actual_cvx = Dict(key => value for (key, value) in actual_cvx if !isempty(value))

# filenums = [2.15, 2.16, 2.17, 2.18, 3.13]
# filenames = ["problem" * string(filenum) * ".gms" for filenum in filenums]
# gms = Dict()
# for filename in filenames
#     try
#         gm = GAMS_to_GlobalModel(OCTHaGOn.GAMS_DIR, filename)
#         types = JuMP.list_of_constraint_types(gm.model)
#         if !isempty(types)
#             total_constraints = sum(length(all_constraints(gm.model, type[1], type[2])) for type in types)
#         end
#         gms[filename] = gm
#     catch
#         throw(OCTHaGOnException(filename * " has an import issue."))
#     end
# end

# gm = gms[filenames[1]]
# set_optimizer(gm, CPLEX.Optimizer)

# Functions to check
# g(x) = x^2/(x^2+1)

# Plotting
# using Plots
# bbl = gm.bbls[end]
# scatter(bbl.X[:,1], bbl.X[:,2], bbl.Y .>= 0)
# scatter(bbl.X[:,1], bbl.X[:,2], bbl.curvatures)

# pos_idxs = findall(x -> x >= 0, bbl.Y)
# neg_idxs = findall(x -> x < 0, bbl.Y)
# scatter(bbl.X[pos_idxs, 1], bbl.X[pos_idxs, 2], color="green")
# scatter!(bbl.X[neg_idxs, 1], bbl.X[neg_idxs, 2], color="red", legend=false)

# The implications of convexity finding.
# This allows us to exploit the properties of convex functions in the domains they are relevant, and come up with appropriate heuristic in the domains where they are invalid.


text:
export TigerProblemEnv

const REWARD_OF_LISTEN = -1
const REWARD_OF_TIGER = -100
const REWARD_OF_TREASURE = 10

"""
    TigerProblemEnv(;rng=Random>GLOBAL_RNG)

Here we use the [The Tiger
Proglem](https://cw.fel.cvut.cz/old/_media/courses/a4m36pah/pah-pomdp-tiger.pdf)
to demonstrate how to write a [POMDP](https://en.wikipedia.org/wiki/Partially_observable_Markov_decision_process) problem.
"""
Base.@kwdef mutable struct TigerProblemEnv <: AbstractEnv
    obs_prob::Float64 = 0.85
    action::Union{Nothing,Symbol} = nothing
    rng::AbstractRNG = Random.GLOBAL_RNG
    tiger_pos::Int = rand(rng, 1:2)
end

Random.seed!(env::TigerProblemEnv, s) = seed!(env.rng, s)

RLBase.action_space(::TigerProblemEnv) = Space((:listen, :open_left, :open_right))

(env::TigerProblemEnv)(action) = env.action = action

function RLBase.reward(env::TigerProblemEnv)
    if env.action == :listen
        REWARD_OF_LISTEN
    elseif (env.action == :open_left && env.tiger_pos == 1) ||
           (env.action == :open_right && env.tiger_pos == 2)
        REWARD_OF_TIGER
    else
        REWARD_OF_TREASURE
    end
end

RLBase.is_terminated(env::TigerProblemEnv) = !isnothing(env.action) && env.action != :listen

function RLBase.reset!(env::TigerProblemEnv)
    env.tiger_pos = rand(env.rng, 1:2)
    env.action = nothing
end

"""
The main difference compared to other environments is that, now we have two
kinds of *states*. The **observation** and the **internal state**. By default we
return the **observation**.
"""
RLBase.state(env::TigerProblemEnv) = state(env, Observation{Int}())

function RLBase.state(env::TigerProblemEnv, ::Observation)
    if isnothing(env.action)
        # game not started yet
        1
    elseif env.action == :listen
        if rand(env.rng) < env.obs_prob
            env.tiger_pos + 1 # shift by 1
        else
            2 - env.tiger_pos + 1 # shift by 1
        end
    else
        4 # terminal state, a dummy state
    end
end

RLBase.state(env::TigerProblemEnv, ::InternalState) = env.tiger_pos

RLBase.state_space(env::TigerProblemEnv) = state_space(env, Observation{Int}())
RLBase.state_space(env::TigerProblemEnv, ::Observation) = Space(1:4)
RLBase.state_space(env::TigerProblemEnv, ::InternalState) = Space(1:2)

RLBase.NumAgentStyle(::TigerProblemEnv) = SINGLE_AGENT
RLBase.DynamicStyle(::TigerProblemEnv) = SEQUENTIAL
RLBase.ActionStyle(::TigerProblemEnv) = MINIMAL_ACTION_SET
RLBase.InformationStyle(::TigerProblemEnv) = IMPERFECT_INFORMATION
RLBase.StateStyle(::TigerProblemEnv) = (Observation{Int}(), InternalState{Int}())
RLBase.RewardStyle(::TigerProblemEnv) = STEP_REWARD
RLBase.UtilityStyle(::TigerProblemEnv) = GENERAL_SUM
RLBase.ChanceStyle(::TigerProblemEnv) = STOCHASTIC


text:
function make_operation(@nospecialize(op), args)
    if op === (*)
        args = filter(!_isone, args)
        if isempty(args)
            return 1
        end
    elseif op === (+)
        args = filter(!_iszero, args)
        if isempty(args)
            return 0
        end
    end
    return op(args...)
end

function build_expr(head::Symbol, args)
    ex = Expr(head)
    append!(ex.args, args)
    ex
end

# used in parsing
isblock(x) = length(x) == 1 && x[1] isa Expr && x[1].head == :block
function flatten_expr!(x)
    isb = isblock(x)
    if isb
        x = MacroTools.striplines(x[1])
        filter!(z->z isa Symbol || z.head != :line, x.args)
        x = (x.args...,)
    end
    x
end

function detime_dvs(op::Term)
  if op.op isa Sym
      Sym{Real}(nameof(op.op))
  else
      Term{Real}(op.op,detime_dvs.(op.args))
  end
end
detime_dvs(op) = op

function retime_dvs(op::Sym,dvs,iv)
    Sym{FnType{Tuple{symtype(iv)}, Real}}(nameof(op))(iv)
end

function retime_dvs(op::Term, dvs, iv)
    similarterm(op, op.op, retime_dvs.(op.args,(dvs,),(iv,)))
end
retime_dvs(op,dvs,iv) = op

is_derivative(O::Term) = isa(O.op, Differential)
is_derivative(::Any) = false

"""
    get_variables(O) -> Vector{Union{Sym, Term}}

Returns the variables in the expression. Note that the returned variables are
not wrapped in the `Num` type.

# Examples
```julia
julia> @parameters t
(t,)

julia> @variables x y z(t)
(x, y, z(t))

julia> ex = x + y + sin(z)
(x + y) + sin(z(t))

julia> ModelingToolkit.get_variables(ex)
3-element Vector{Any}:
 x
 y
 z(t)
```
"""
get_variables(e::Num, varlist=nothing) = get_variables(value(e), varlist)
get_variables!(vars, e, varlist=nothing) = vars

is_singleton(e::Term) = e.op isa Sym
is_singleton(e::Sym) = true
is_singleton(e) = false

get_variables!(vars, e::Number, varlist=nothing) = vars

function get_variables!(vars, e::Symbolic, varlist=nothing)
    if is_singleton(e)
        if isnothing(varlist) || any(isequal(e), varlist)
            push!(vars, e)
        end
    else
        foreach(x -> get_variables!(vars, x, varlist), e.args)
    end
    return (vars isa AbstractVector) ? unique!(vars) : vars
end

function get_variables!(vars, e::Equation, varlist=nothing)
  get_variables!(vars, e.rhs, varlist)
end

get_variables(e, varlist=nothing) = get_variables!([], e, varlist)

modified_states!(mstates, e::Equation, statelist=nothing) = get_variables!(mstates, e.lhs, statelist)

# variable substitution
# Piracy but mild
"""
    substitute(expr, s::Pair)
    substitute(expr, s::Dict)
    substitute(expr, s::Vector)

Performs the substitution on `expr` according to rule(s) `s`.

# Examples
```julia
julia> @parameters t
(t,)

julia> @variables x y z(t)
(x, y, z(t))

julia> ex = x + y + sin(z)
(x + y) + sin(z(t))

julia> substitute(ex, Dict([x => z, sin(z) => z^2]))
(z(t) + y) + (z(t) ^ 2)
```
"""
substitute(expr::Num, s::Union{Pair, Vector, Dict}; kw...) = Num(substituter(s)(value(expr); kw...))
# TODO: move this to SymbolicUtils
substitute(expr::Term, s::Pair; kw...) = substituter([s[1] => s[2]])(expr; kw...)
substitute(expr::Term, s::Vector; kw...) = substituter(s)(expr; kw...)

substituter(pair::Pair) = substituter((pair,))
function substituter(pairs)
    dict = Dict(to_symbolic(k) => to_symbolic(v)  for (k, v) in pairs)
    (expr; kw...) -> SymbolicUtils.substitute(expr, dict; kw...)
end

macro showarr(x)
    n = string(x)
    quote
        y = $(esc(x))
        println($n, " = ", summary(y))
        Base.print_array(stdout, y)
        println()
        y
    end
end

@deprecate substitute_expr!(expr,s) substitute(expr,s)

function states_to_sym(states::Set)
    function _states_to_sym(O)
        if O isa Equation
            Expr(:(=), _states_to_sym(O.lhs), _states_to_sym(O.rhs))
        elseif O isa Term
            if isa(O.op, Sym)
                O in states && return tosymbol(O)
                # dependent variables
                return build_expr(:call, Any[O.op.name; _states_to_sym.(O.args)])
            else
                return build_expr(:call, Any[O.op; _states_to_sym.(O.args)])
            end
        elseif O isa Num
            return _states_to_sym(value(O))
        else
            return toexpr(O)
        end
    end
end
states_to_sym(states) = states_to_sym(Set(states))

"""
    toparam(s::Sym) -> Sym{<:Parameter}

Maps the variable to a paramter.
"""
toparam(s::Sym) = Sym{Parameter{symtype(s)}}(s.name)
toparam(s::Sym{<:Parameter}) = s

"""
    tovar(s::Sym) -> Sym{Real}
    tovar(s::Sym{<:Parameter}) -> Sym{Real}

Maps the variable to a state.
"""
tovar(s::Sym{<:Parameter}) = Sym{symtype(s)}(s.name)
tovar(s::Sym) = s

Base.Symbol(x::Union{Num,Symbolic}) = tosymbol(x)
tosymbol(x; kwargs...) = x
tosymbol(x::Sym; kwargs...) = nameof(x)
tosymbol(t::Num; kwargs...) = tosymbol(value(t); kwargs...)

"""
    tosymbol(x::Union{Num,Symbolic}; states=nothing, escape=true) -> Symbol

Convert `x` to a symbol. `states` are the states of a system, and `escape`
means if the target has escapes like `val"y⦗t⦘"`. If `escape` then it will only
output `y` instead of `y⦗t⦘`.

# Examples
```julia
julia> @parameters t; @variables z(t)
(z(t),)

julia> ModelingToolkit.tosymbol(z)
Symbol("z⦗t⦘")
```
"""
function tosymbol(t::Term; states=nothing, escape=true)
    if t.op isa Sym
        if states !== nothing && !(any(isequal(t), states))
            return nameof(t.op)
        end
        op = nameof(t.op)
        args = t.args
    elseif t.op isa Differential
        term = diff2term(t)
        op = Symbol(operation(term))
        args = arguments(term)
    else
        @goto err
    end

    return escape ? Symbol(op, "⦗", join(args, ", "), "⦘") : op
    @label err
    error("Cannot convert $t to a symbol")
end

"""
    makesym(x::Union{Num,Symbolic}, kwargs...) -> Sym

`makesym` takes the same arguments as [`tosymbol`](@ref), but it converts a
`Term` in the form of `x(t)` to a `Sym` in the form of `x⦗t⦘`.

# Examples
```julia
julia> @parameters t; @variables x(t)
(x(t),)

julia> ModelingToolkit.makesym(x)
x⦗t⦘
```
"""
makesym(t::Symbolic; kwargs...) = Sym{symtype(t)}(tosymbol(t; kwargs...))
makesym(t::Num; kwargs...) = makesym(value(t); kwargs...)

function lower_varname(var::Term, idv, order)
    order == 0 && return var
    name = string(nameof(var.op))
    underscore = 'ˍ'
    idx = findlast(underscore, name)
    append = string(idv)^order
    if idx === nothing
        newname = Symbol(name, underscore, append)
    else
        nidx = nextind(name, idx)
        newname = Symbol(name[1:idx], name[nidx:end], append)
    end
    return Sym{symtype(var.op)}(newname)(var.args[1])
end

function lower_varname(t::Term, iv)
    var, order = var_from_nested_derivative(t)
    lower_varname(var, iv, order)
end
lower_varname(t::Sym, iv) = t

function lower_mapnames(umap::AbstractArray{<:Pair}) where T
    [value(k) => value(v) for (k, v) in umap]
end
function lower_mapnames(umap::AbstractArray{<:Pair},name) where T
    [lower_varname(value(k), name) => value(v) for (k, v) in umap]
end
lower_mapnames(umap::AbstractArray{<:Number}) = umap # Ambiguity
lower_mapnames(umap::AbstractArray{<:Number},name) = umap

function flatten_differential(O::Term)
    @assert is_derivative(O) "invalid differential: $O"
    is_derivative(O.args[1]) || return (O.args[1], O.op.x, 1)
    (x, t, order) = flatten_differential(O.args[1])
    isequal(t, O.op.x) || throw(ArgumentError("non-matching differentials on lhs: $t, $(O.op.x)"))
    return (x, t, order + 1)
end

"""
    diff2term(x::Term) -> Term
    diff2term(x) -> x

Convert a differential variable to a `Term`. Note that it only takes a `Term`
not a `Num`.
```julia
julia> ModelingToolkit.diff2term(ModelingToolkit.value(D(D(x))))
xˍtt(t)
```
"""
function diff2term(O)
    isa(O, Term) || return O
    if is_derivative(O)
        (x, t, order) = flatten_differential(O)
        return lower_varname(x, t, order)
    end
    return Term{Real}(O.op, diff2term.(O.args))
end


text:
"""
    identify(indexmode, index)

Identify an `IndexMode` and its fields from the index content
and the existing `IndexMode`.

These methods guess which mode is most appropriate, encoding information
about the index in types that can later be used for dispatch. They also let
you fill in part of the information you need to specify, and guess the rest.

An example of the usefulness of identifying index traits up-front is if we
check that the index is ordered, we can be sure it remains ordered
unless it is indexed with a `Vector`. This means we will always use the
correct `searchsorted` method for it.
"""
function identify end
# No more identification required for some types
identify(mode::IndexMode, dimtype::Type, index) = mode
# AutoMode
identify(mode::AutoMode, dimtype::Type, index::AbstractArray{T}) where T = begin
    # A mixed type index is Categorical
    if isconcretetype(T)
        identify(Sampled(), dimtype, index)
    else
        identify(Categorical(), dimtype, index)
    end
end
identify(mode::AutoMode, dimtype::Type, index::AbstractArray{<:CategoricalEltypes}) =
    order(mode) isa AutoOrder ? Categorical(_orderof(index)) : Categorical(order(mode))
identify(mode::AutoMode, dimtype::Type, index::Val) =
    order(mode) isa AutoOrder ? Categorical(Unordered()) : Categorical(order(mode))
# Sampled
identify(mode::AbstractSampled, dimtype::Type, index) = begin
    o = identify(order(mode), dimtype, index)
    sp = identify(span(mode), dimtype, index)
    sa = identify(sampling(mode), sp, dimtype, index)
    rebuild(mode; order=o, span=sp, sampling=sa)
end
identify(order::Order, dimtype::Type, index) = order
identify(order::AutoOrder, dimtype::Type, index) = _orderof(index)
# Span
identify(span::AutoSpan, dimtype::Type, index::Union{AbstractArray,Val}) = Irregular()
identify(span::AutoSpan, dimtype::Type, index::AbstractRange) = Regular(step(index))
identify(span::Regular{AutoStep}, dimtype::Type, index::Union{AbstractArray,Val}) = _arraynosteperror()
identify(span::Regular, dimtype::Type, index::Union{AbstractArray,Val}) = span
identify(span::Regular{AutoStep}, dimtype::Type, index::AbstractRange) = Regular(step(index))
identify(span::Regular, dimtype::Type, index::AbstractRange) = begin
    step(span) isa Number && !(step(span) ≈ step(index)) && _steperror(index, span)
    span
end
identify(span::Irregular{AutoBounds}, dimtype, index) = Irregular(nothing, nothing)
identify(span::Irregular{<:Tuple}, dimtype, index) = span
identify(span::Explicit, dimtype, index) = span
# Sampling
identify(sampling::AutoSampling, span::Span, dimtype::Type, index) = Points()
identify(sampling::AutoSampling, span::Explicit, dimtype::Type, index) = Intervals(Center())
identify(sampling::Points, span::Span, dimtype::Type, index) = sampling
identify(sampling::Points, span::Explicit, dimtype::Type, index) = _explicitpoints_error()
identify(sampling::Intervals, span::Span, dimtype::Type, index) =
    rebuild(sampling, identify(locus(sampling), dimtype, index))
# Locus
identify(locus::AutoLocus, dimtype::Type, index) = Center()
identify(locus::Locus, dimtype::Type, index) = locus

_orderof(index::AbstractUnitRange) = Ordered()
_orderof(index::AbstractRange) = Ordered(index=_indexorder(index))
_orderof(index::Val) = _detectorder(unwrap(index))
_orderof(index::AbstractArray) = _detectorder(index)

function _detectorder(index)
    local sorted, indord
    # This is awful. But we don't know if we can
    # call `issorted` on the contents of `index`.
    # This may be resolved by: https://github.com/JuliaLang/julia/pull/37239
    try
        indord = _indexorder(index)
        sorted = issorted(index; rev=isrev(indord))
    catch
        sorted = false
    end
    sorted ? Ordered(index=indord) : Unordered()
end

_indexorder(index) = first(index) <= last(index) ? ForwardIndex() : ReverseIndex()

@noinline _explicitpoints_error() =
    throw(ArgumentError("Cannot use Explicit span with Points sampling"))
@noinline _steperror(index, span) =
    throw(ArgumentError("mode step $(step(span)) does not match index step $(step(index))"))
@noinline _arraynosteperror() =
    throw(ArgumentError("`Regular` must specify `step` size with an index other than `AbstractRange`"))


text:
#=
  Original Author: Anastasia Yendiki

  Copyright © 2022 The General Hospital Corporation (Boston, MA) "MGH"

  Terms and conditions for use, reproduction, distribution and contribution
  are found in the 'FreeSurfer Software License Agreement' contained
  in the file 'LICENSE' found in the FreeSurfer distribution, and here:

  https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense

  Reporting: freesurfer@nmr.mgh.harvard.edu
=#

using Statistics

export GQI, gqi_rec, gqi_write

"Container for outputs of a GQI fit"
struct GQI
  odf::MRI
  peak::Vector{MRI}
  qa::Vector{MRI}
end

"""
    GQIwork{T}

Pre-allocated workspace for GQI reconstruction computations

- `T::DataType`                : Data type for computations (default: `Float32`)
- `nvol::Int`                  : Number of volumes in DWI series
- `nvert::Int`                 : Number of ODF vertices (on the half sphere)
- `isort::Vector{Vector{Int}}` : Indices of ODF peak vertices (sorted) [`nvert`]
- `s::Vector{Vector{T}}`       : DWI signal [`nvol`]
- `o::Vector{Vector{T}}`       : ODF amplitudes [`nvert`]
- `odf_peak::Vector{Vector{T}}`: ODF amplitudes at local peaks [`nvert`]
- `faces::Matrix{Int}`         : ODF faces (on the half sphere) [nvert x 3]
- `A::Matrix{T}`               : System matrix [nvert x nvol]
"""
struct GQIwork{T}
  nvol::Int
  nvert::Int
  isort::Vector{Vector{Int}}
  s::Vector{Vector{T}}
  o::Vector{Vector{T}}
  odf_peak::Vector{Vector{T}}
  faces::Matrix{Int}
  A::Matrix{T}

  function GQIwork(bval::Vector{Float32}, bvec::Matrix{Float32}, odf_dirs::ODF=sphere_642, σ::Float32=Float32(1.25), T::DataType=Float32)

    # Number of volumes in DWI series
    nvol     = length(bval)

    # Number of ODF vertices (on the half sphere)
    nvert    = div(size(odf_dirs.vertices, 1), 2)

    # Indices of sorted ODF peaks
    isort    = [Vector{Int}(undef, nvert) for tid in 1:Threads.nthreads()]

    # DWI signal
    s        = [Vector{T}(undef, nvol) for tid in 1:Threads.nthreads()]

    # ODF amplitudes
    o        = [Vector{T}(undef, nvert) for tid in 1:Threads.nthreads()]

    # ODF amplitudes at local peaks
    odf_peak = [Vector{T}(undef, nvert) for tid in 1:Threads.nthreads()]

    # ODF faces (on the half sphere)
    faces = copy(odf_dirs.faces)
    faces[faces .> nvert] .-= nvert

    # System matrix
    A = Matrix{T}(undef, nvert, nvol)
    bq_vector = bvec .* (sqrt.(bval * T(0.01506)) * T(σ / π))
    A .= sinc.(odf_dirs.vertices[nvert+1:end, :] * bq_vector')

    new{T}(
      nvol,
      nvert,
      isort,
      s,
      o,
      odf_peak,
      faces,
      A
    )
  end
end

"""
    gqi_rec(dwi::MRI, mask::MRI, odf_dirs::ODF=sphere_642, σ::Float32=Float32(1.25))

Perform generalized q-space imaging (GQI) reconstruction of DWIs, and return a
`GQI` structure.

If you use this method, please cite:
Fang-Cheng Yeh, et al. (2010). Generalized q-sampling imaging. IEEE Transactions on Medical Imaging, 29(9), 1626–1635. https://doi.org/10.1109/TMI.2010.2045126

# Arguments
- `dwi::MRI`: A series of DWIs, stored in an `MRI` structure with valid `.bvec`
   and `.bval` fields
- `mask::MRI`: A brain mask volume, stored in an `MRI` structure
- `odf_dirs::ODF=sphere_642`: The vertices and faces of the ODF tessellation,
  stored in an `ODF` structure
- `σ::Float32=Float32(1.25)`: Diffusion sampling length factor

# Output
In the `GQI` structure:
- `.odf`: ODF amplitudes on the half sphere
- `.peak`: Orientation vectors of the 3 peak ODF amplitudes
- `.qa`: Quantitative anisotropy for each of the 3 peak orientations

"""
function gqi_rec(dwi::MRI, mask::MRI, odf_dirs::ODF=sphere_642, σ::Float32=Float32(1.25))

  if isempty(dwi.bval)
    error("Missing b-value table from input DWI structure")
  end

  if isempty(dwi.bvec)
    error("Missing gradient table from input DWI structure")
  end

  npeak = 3;

  W = GQIwork(dwi.bval, dwi.bvec, odf_dirs, σ)

  odf  = MRI(mask, W.nvert, Float32)
  peak = Vector{MRI}(undef, npeak)
  qa   = Vector{MRI}(undef, npeak)

  for ipeak in 1:npeak
    peak[ipeak] = MRI(mask, 3, Float32)
    qa[ipeak]   = MRI(mask, 1, Float32)
  end

  Threads.@threads for iz in 1:size(dwi.vol, 3)
    for iy in 1:size(dwi.vol, 2)
      for ix in 1:size(dwi.vol, 1)
        mask.vol[ix, iy, iz] == 0 && continue

        tid = Threads.threadid()

        W.s[tid] .= dwi.vol[ix, iy, iz, :]
        W.s[tid][W.s[tid] .< 0] .= 0

        maximum(W.s[tid]) == 0 && continue

        mul!(W.o[tid], W.A, W.s[tid])
        odf.vol[ix, iy, iz, :] = W.o[tid]

        odfmin = minimum(W.o[tid])

        nvalid = find_peaks!(W)

        n = min(nvalid, npeak)

        for ipeak in 1:n
          peak[ipeak].vol[ix, iy, iz, :] =
            odf_dirs.vertices[W.isort[tid][ipeak], :]

          qa[ipeak].vol[ix, iy, iz] =
            odf.vol[ix, iy, iz, W.isort[tid][ipeak]] - odfmin
        end
      end
    end
  end

  odfmax = maximum(mean(odf.vol, dims=4))

  for ipeak in 1:npeak
    qa[ipeak].vol /= odfmax
  end

  return GQI(odf, peak, qa)
end

"""
    find_peaks!(W::GQIwork)

Find the vertices whose amplitudes are local peaks and return them sorted
(assume that ODF amplitudes have been computed)
"""
function find_peaks!(W::GQIwork)

  tid = Threads.threadid()

  W.odf_peak[tid] .= W.o[tid]
  @views W.odf_peak[tid][W.faces[W.o[tid][W.faces[:,2]] .>=
                                 W.o[tid][W.faces[:,1]] .||
                                 W.o[tid][W.faces[:,3]] .>=
                                 W.o[tid][W.faces[:,1]], 1]] .= 0
  @views W.odf_peak[tid][W.faces[W.o[tid][W.faces[:,1]] .>=
                                 W.o[tid][W.faces[:,2]] .||
                                 W.o[tid][W.faces[:,3]] .>=
                                 W.o[tid][W.faces[:,2]], 2]] .= 0
  @views W.odf_peak[tid][W.faces[W.o[tid][W.faces[:,2]] .>=
                                 W.o[tid][W.faces[:,3]] .||
                                 W.o[tid][W.faces[:,1]] .>=
                                 W.o[tid][W.faces[:,3]], 3]] .= 0

  sortperm!(W.isort[tid], W.odf_peak[tid], rev=true)

  return length(findall(W.odf_peak[tid] .> 0))
end

"""
    gqi_write(gqi::GQI, basename::String)

Write the volumes from a `GQI` structure that was created by `gqi_rec()`
to files whose names start with the specified base name.
"""
function gqi_write(gqi::GQI, basename::String)

  for var in fieldnames(GQI)
    vartype = fieldtype(GQI, var)

    if vartype == MRI
      fname = basename * "_" * string(var) * ".nii.gz"
      mri_write(getfield(gqi, var), fname)
    elseif vartype == Vector{MRI}
      for ivol = 1:length(getfield(gqi, var))
        fname = basename * "_" * string(var) * string(ivol) * ".nii.gz"
        mri_write(getfield(gqi, var)[ivol], fname)
      end
    end
  end
end

max_stars_repo_name: freesurfer/jupysurf

max_issues_repo_name: freesurfer/jupysurf

max_forks_repo_name: freesurfer/jupysurf


text:
"""
    brim(N::NT, L::Dict{E,Int64}) where {NT<:AbstractEcologicalNetwork,E}

Uses BRIM to optimize the modularity of an ecological network. The `L` argument
is a dictionary mapping every species in the network to its module. This
function returns a tuple of the network and its module assignment.

#### References

- Barber, M.J., 2007. Modularity and community detection in bipartite networks.
  Phys. Rev. E 76, 066102. https://doi.org/10.1103/PhysRevE.76.066102

- Newman, M.E., 2006. Modularity and community structure in networks.
  Proceedings of the national academy of sciences 103, 8577–8582.

- Thébault, E., 2013. Identifying compartments in presence–absence matrices and
  bipartite networks: insights into modularity measures. Journal of Biogeography
  40, 759–768. https://doi.org/10.1111/jbi.12015
"""
function brim(N::NT, L::Dict{E,Int64}) where {NT<:AbstractEcologicalNetwork,E}
  for s in species(N)
    @assert haskey(L, s)
  end
  EcologicalNetworks.tidy_modules!(L)

  old_Q = Q(N, L)
  new_Q = old_Q+0.00001

  m = links(N)

  c = length(unique(collect(values(L))))

  R = zeros(Int64, (richness(N; dims=1),c))
  for si in eachindex(species(N; dims=1))
    R[si,L[species(N; dims=1)[si]]] = 1
  end

  T = zeros(Int64, (richness(N; dims=2),c))
  for sj in eachindex(species(N; dims=2))
    T[sj,L[species(N; dims=2)[sj]]] = 1
  end

  d_in, d_out = degree(N; dims=2), degree(N; dims=1)
  s_d_in = [d_in[s] for s in species(N; dims=2)]
  s_d_out = [d_out[s] for s in species(N; dims=1)]
  P = kron(s_d_out, s_d_in')./m
  B = N.edges .- P

  while old_Q < new_Q

    t_tilde = B*T
    R = map(Int64, t_tilde .== maximum(t_tilde; dims=2))
    r_tilde = B'*R
    T = map(Int64, r_tilde .== maximum(r_tilde; dims=2))
    S = vcat(R, T)
    L = Dict(zip(species(N), vec(mapslices(r -> StatsBase.sample(findall(x->x!=0, r)), S; dims=2))))
    EcologicalNetworks.tidy_modules!(L)
    old_Q = new_Q
    new_Q = Q(N,L)

  end

  EcologicalNetworks.tidy_modules!(L)
  return (N, L)

end


text:
using JLD2, Test

struct SingleFieldWrapper{T}
    x::T
end
Base.:(==)(a::SingleFieldWrapper, b::SingleFieldWrapper) = a.x == b.x

struct MultiFieldWrapper{T}
    x::T
    y::Int
end
Base.:(==)(a::MultiFieldWrapper, b::MultiFieldWrapper) = (a.x == b.x && a.y == b.y)

struct UntypedWrapper
    x
end
Base.:(==)(a::UntypedWrapper, b::UntypedWrapper) = a.x == b.x

# This is a type with a custom serialization, where the original type has data
# but the custom serialization is empty
struct CSA
    x::Ptr{Cvoid}
end
a = CSA(Ptr{Cvoid}(0))

struct CSASerialization end
JLD2.writeas(::Type{CSA}) = CSASerialization
function JLD2.wconvert(::Type{CSASerialization}, x::CSA)
    global converted = true
    CSASerialization()
end
function JLD2.rconvert(::Type{CSA}, x::CSASerialization)
    global converted = true
    CSA(Ptr{Cvoid}(0))
end

# This is a type with a custom serialization, where the original type has no
# data but the custom serialization does
struct CSB end
b = CSB()

struct CSBSerialization
    x::Int
end
JLD2.writeas(::Type{CSB}) = CSBSerialization
function JLD2.wconvert(::Type{CSBSerialization}, x::CSB)
    global converted = true
    CSBSerialization(9018620829326368991)
end
function JLD2.rconvert(::Type{CSB}, x::CSBSerialization)
    global converted = true
    x.x == 9018620829326368991 ? CSB() : error("invalid deserialized data")
end

# This is a type where the custom serialized data can be stored inline when it
# is a field of another type, but the original data could not
mutable struct CSC
    x::Vector{Int}
end
Base.:(==)(a::CSC, b::CSC) = a.x == b.x
c = CSC(rand(Int, 2))

struct CSCSerialization
    a::Int
    b::Int
end
JLD2.writeas(::Type{CSC}) = CSCSerialization
function JLD2.wconvert(::Type{CSCSerialization}, x::CSC)
    global converted = true
    CSCSerialization(x.x[1], x.x[2])
end
function JLD2.rconvert(::Type{CSC}, x::CSCSerialization)
    global converted = true
    CSC([x.a, x.b])
end

# This is a type where the original data could be stored inline when it is a
# field of another type, but the custom serialized data cannot
struct CSD
    a::Int
    b::Int
end
d = CSD(rand(Int), rand(Int))

struct CSDSerialization
    x::Vector{Int}
end
JLD2.writeas(::Type{CSD}) = CSDSerialization
function JLD2.wconvert(::Type{CSDSerialization}, x::CSD)
    global converted = true
    CSDSerialization([x.a, x.b])
end
function JLD2.rconvert(::Type{CSD}, x::CSDSerialization)
    global converted = true
    CSD(x.x[1], x.x[2])
end

# This is a type that gets written as an array
struct CSE
    a::Int
    b::Int
    c::Int
end
e = CSE(rand(Int), rand(Int), rand(Int))

JLD2.writeas(::Type{CSE}) = Vector{Int}
function JLD2.wconvert(::Type{Vector{Int}}, x::CSE)
    global converted = true
    [x.a, x.b, x.c]
end
function JLD2.rconvert(::Type{CSE}, x::Vector{Int})
    global converted = true
    CSE(x[1], x[2], x[3])
end

# This is a type that is custom-serialized as an Int
struct CSF
    x::Int
end
f = CSF(rand(Int))

JLD2.writeas(::Type{CSF}) = Int
function Base.convert(::Type{Int}, x::CSF)
    global converted = true
    x.x
end
function Base.convert(::Type{CSF}, x::Int)
    global converted = true
    CSF(x)
end

# This is a type that is custom-serialized as a String
struct CSG
    x::Int
end
g = CSG(rand(Int))

JLD2.writeas(::Type{CSG}) = String
function Base.convert(::Type{String}, x::CSG)
    global converted = true
    convert(String, string(x.x))
end
function Base.convert(::Type{CSG}, x::String)
    global converted = true
    CSG(parse(Int, x))
end

# This is a type that is custom-serialized as a DataType
struct CSH
    T::DataType
    N::Int
end
h = CSH(Int, 2)

JLD2.writeas(::Type{CSH}) = DataType
function Base.convert(::Type{DataType}, x::CSH)
    global converted = true
    Array{x.T, x.N}
end
function Base.convert(::Type{CSH}, x::Type{Array{T,N}}) where {T,N}
    global converted = true
    CSH(T,N)
end

# This is a type that is custom-serialized as a Union
struct CSK
    T::DataType
    S::DataType
end
k = CSK(Int, Float64)
Base.:(==)(x::CSK, y::CSK) = (x.T == y.T && x.S == y.S) ||
                             (x.T == y.S && x.S == y.T)

JLD2.writeas(::Type{CSK}) = Union
function Base.convert(::Type{Union}, x::CSK)
    global converted = true
    Union{x.T,x.S}
end
function Base.convert(::Type{CSK}, x::Union)
    global converted = true
    CSK(x.types...)
end

function write_tests(file, prefix, obj)
    write(file, prefix, obj)
    write(file, "$(prefix)_singlefieldwrapper", SingleFieldWrapper(obj))
    write(file, "$(prefix)_multifieldwrapper", MultiFieldWrapper(obj, 5935250212119237787))
    write(file, "$(prefix)_untypedwrapper", UntypedWrapper(obj))
    write(file, "$(prefix)_arr", [obj])
    write(file, "$(prefix)_empty_arr", typeof(obj)[])
end

function read_tests(file, prefix, obj)
    global converted = false
    @test read(file, prefix) == obj
    @test converted
    @test read(file, "$(prefix)_singlefieldwrapper") == SingleFieldWrapper(obj)
    @test read(file, "$(prefix)_multifieldwrapper") == MultiFieldWrapper(obj, 5935250212119237787)
    @test read(file, "$(prefix)_untypedwrapper") == UntypedWrapper(obj)
    arr = read(file, "$(prefix)_arr")
    @test typeof(arr) == Vector{typeof(obj)} && length(arr) == 1 && arr[1] == obj
    empty_arr = read(file, "$(prefix)_empty_arr")
    @test typeof(empty_arr) == Vector{typeof(obj)} && length(empty_arr) == 0
end

fn = joinpath(tempdir(),"test.jld")
println(fn)
file = jldopen(fn, "w")
write_tests(file, "a", a)
write_tests(file, "b", b)
write_tests(file, "c", c)
write_tests(file, "d", d)
write_tests(file, "e", e)
write_tests(file, "f", f)
write_tests(file, "g", g)
write_tests(file, "h", h)
write_tests(file, "k", k)
close(file)

file = jldopen(fn, "r")
read_tests(file, "a", a)
read_tests(file, "b", b)
read_tests(file, "c", c)
read_tests(file, "d", d)
read_tests(file, "e", e)
read_tests(file, "f", f)
read_tests(file, "g", g)
read_tests(file, "h", h)
read_tests(file, "k", k)
close(file)


text:
# # Intorduction

# As we know, triangular meshes can be stored in a computer in multiple different ways, each having strength and weaknesses in a particular case at hand. But it is not always clear which data structure would be most suitable for a specific task. Thus it is wise to write a data structure generic code which is the precise purpose of this package for closed oriented closed surfaces.

# The most straightforward representation of triangular mesh topology is in array `Array{Faces{3,Int},1}` containing a list of triangular faces which are defined by their vertices. That as name stands allows quick iteration over faces and also edges. However, often in a numerical code one wants not only to iterate over faces or vertices but also in case of singularity subtraction, integration and local property estimation like in normal vector and curvature calculations to know what are neighbouring vertices surrounding a given vertex while keeping track of the orientation of the normals. Also, one wishes to modify the topology itself by collapsing, flipping and splitting edges. And that is why different data structures are needed for different problems.

# Fortunately, it is possible to abstract mesh topology queries through iterators:
# ```@docs
# Faces
# Edges
# ```
# and circulators:
# ```@docs
# VertexRing
# EdgeRing
# ```

# ## API

# The package implements multiple kinds of data structures. The simplest one is `PlainDS` one which stores a list of faces and is just an alias to `Array{Faces{3,Int},1}`. As an example of how that works, let's define the data structure.

using GeometryTypes
using SurfaceTopology

faces = Face{3,Int64}[
    [1, 12, 6], [1, 6, 2], [1, 2, 8], [1, 8, 11], [1, 11, 12], [2, 6, 10], [6, 12, 5],
    [12, 11, 3], [11, 8, 7], [8, 2, 9], [4, 10, 5], [4, 5, 3], [4, 3, 7], [4, 7, 9],
    [4, 9, 10], [5, 10, 6], [3, 5, 12], [7, 3, 11], [9, 7, 8], [10, 9, 2]
]

# We can use the data structure `PlainDS` for the queries. The iterators, for example.
collect(Faces(faces))
# and
collect(Edges(faces))
# giving us desirable output.

# We can also ask what neighbouring vertices and edges for a particular vertex by using circulators. For this simple data structure that requires us to do a full lookup on the face list, which is nicely abstracted away:
collect(VertexRing(3,faces))
# and
collect(EdgeRing(3,faces))
# In practice, one should use `EdgeRing` over `VertexRing` since, in the latter one, vertices are not ordered and thus can not be used for example to deduce the direction of the normal vector.

# ## Data structures

# The same API works for all other data structures. There is a data structure `CachedDS` built on top of `PlainDS` stores closest vertices (vertex ring). Then there is a data structure `FaceDS` which with `PlainDS` also stores neighbouring faces which have a common edge. And then there is the most commonly used data structure in numerics `HalfEdgeDS` (implemented as `EdgeDS`).

# The most straightforward extension of `PlainDS` is just a plain caching of neighbouring vertices for each vertex which are stored in `CacheDS` also with the list of faces.
# ```@docs
# CachedDS
# CachedDS(::SurfaceTopology.PlainDS)
# ```
# which can be initialised from `PlainDS`
cachedtopology = CachedDS(faces)
# And the same API can be used for querries:
collect(VertexRing(3,cachedtopology))

# A more advanced data structure is a face based data structure `FaceDS` which additionally for each face stores three neighbouring face indices.
# ```@docs
# FaceDS
# FaceDS(::SurfaceTopology.PlainDS)
# ```
# which again can be initialised from `PlainDS`
facedstopology = FaceDS(faces)
# and what would one expect
collect(VertexRing(3,facedstopology))
# works.

# All previous ones were some forms of face-based data structures. More common (by my own impression) the numerical world uses edge-based data structures. This package implements half-edge data structure `EdgeDS` which stores a list of edges by three numbers - base vertex index, next edge index and twin edge index.
# ```@docs
# EdgeDS
# ```
# To initialise this datastructure one executes:
edgedstopology = EdgeDS(faces)
#
collect(VertexRing(3,edgedstopology))

# ## Wishlist

# At the moment the package is able only to answer queries, but it would be desirable also to be able to do topological surgery operations. For completeness, those would include.

#   + `edgeflip(topology,edge)`
#   + `edgesplit(topology,edge)`
#   + `edgecollapse(topology,edge)`

# And with them even a method for `defragmenting` the topology (actually trivial if we generalize constructors as in `CachedDS`). Unfortunately, at the moment, I am not working with anything geometry related thus the development of that on my own will be slow. I hope that the clarity and simplicity of this package could serve someone as a first step, and so eventually, topological operations would be implemented out of necessity.


text:
using CSV, DataFrames, PyCall, PyPlot

PyDict(matplotlib["rcParams"])["axes.formatter.use_locale"] = true

function sosread(csvraw)
	csvform = tempname()
	if !(ispath(csvform))
		run(pipeline(`./soscan.awk`, stdin=csvraw, stdout=csvform))
	end
	df = CSV.read(csvform, delim=';')
	df[:Kohort] = df[:År] .- (df[:Ålder] .+ 2)
	return df
end

function agelab(age)
	if age == 85
		return "85\u2013"
	else
		return "$(age)\u2013$(age+4)"
	end
end

function caplot(frame, sex, age)
	subframe = frame[frame[:Ålder].==age, :]
	salab = "$sex $(agelab(age))"
	plot(subframe[:Kohort], log.(subframe[sex]./10^5), "-*", label=salab)
end

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
import Base: copy!, splice!, append!, push!, setindex!, start, next, done
import Base: getindex, map, length, eltype, endof, ndims, size, resize!

resize!(A::GPUArray, newdims::Int...) = resize!(A, newdims)
function resize!(A::GPUArray{T, 1}, newdims::NTuple{1, Int}) where T
    newdims == size(A) && return A
    newlength = newdims[1]
    real_length = length(buffer(A))
    if real_length >= newlength # underlying buffer already big enough
        A.size = newdims
        return A
    else
        B = similar(A, newdims)
        if length(A) > 0
            max_len = min(length(A), newlength) #might also shrink
            copy!(B, (1:max_len,), A, (1:max_len,))
        end
        A.size = newdims
        A.buffer = buffer(B)
        return A
    end
end

function reshape!(A::GPUArray{T, NDim}, newdims::NTuple{NDim, Int}) where {T, NDim}
    size(A) == newdims && return A
    if prod(newdims) == length(A)
        throw(DimensionMismatch("new dimensions $newdims must be consistent with array length $(length(A))"))
    end
    GPUArray(buffer(A), context(A), newdims)
end

"""
This updates an array, even if dimensions and sizes don't match.
Will resize accordingly!
"""
function update!(A::GPUArray{T, N}, value::Array{T, N}) where {T, N}
    size(A) != size(value) && resize!(A, size(value))
    copy!(A, value)
    return
end

function grow_dimensions(
        real_length::Int, _size::Int, additonal_size::Int,
        growfactor::Real = 1.5
    )
    new_dim = round(Int, real_length * growfactor)
    return max(new_dim, additonal_size + _size)
end

push!(v::GPUVector{T}, x) where {T} = push!(v, convert(T, x))
push!(v::GPUVector{T}, x::T) where {T} = append!(v, [x])
push!(v::GPUVector{T}, x::T...) where {T} = append!(v, [x...])

function append!(v::GPUVector{T}, value) where T
    x = array_convert(Vector{T}, value)
    lv, lx = length(v), length(x)
    real_length = length(buffer(v))
    if real_length < (lv + lx)
        resize!(v, grow_dimensions(real_length, lv, lx))
    end
    v.size = (lv + lx,)
    v[(lv + 1) : (lv + lx)] = value
    v
end

function grow_at(v::GPUVector, index::Int, amount::Int)
    resize!(v, length(v) + amount)
    copy!(v, index, v, index + amount, amount)
end

splice!(v::GPUVector{T}, index::Int, x::T) where {T} = (v[index] = x)
function splice!(v::GPUVector{T}, index::Int, x::Vector = T[]) where T
    splice!(v, index:index, map(T, x))
end

function splice!(v::GPUVector{T}, index::UnitRange, x::Vector=T[]) where T
    lenv = length(v)
    elements_to_grow = length(x) - length(index) # -1
    buffer = similar(buffer(v), length(v) + elements_to_grow)
    copy!(v.buffer, 1, buffer, 1, first(index) - 1) # copy first half
    copy!(v.buffer, last(index) + 1, buffer, first(index) + length(x), lenv - last(index)) # shift second half
    v.buffer = buffer
    v.real_length = length(buffer)
    v.size = (v.real_length,)
    copy!(x, 1, buffer, first(index), length(x)) # copy contents of insertion vector
    return
end


text:
using TestPackage2
using Test

@testset "TestPackage2.jl" begin
    @test TestPackage2.answer_to_life_universe() == 42
end


text:
"""
    run_opt(ts_data::ClustData,opt_data::OptDataCEP,opt_config::Dict{String,Any},optimizer::DataType)
Organizing the actual setup and run of the CEP-Problem. This function shouldn't be called by a user, but from within the other `run_opt`-functions
Required elements are:
- `ts_data`: The time-series data.
- `opt_data`: In this case the OptDataCEP that contains information on costs, nodes, techs and for transmission also on lines.
- `opt_config`: This includes all the settings for the design optimization problem formulation.
- `optimizer`: The used optimizer, which could e.g. be Clp: `using Clp` `optimizer=Clp.Optimizer` or Gurobi: `using Gurobi` `optimizer=Gurobi.Optimizer`.
"""
function run_opt(ts_data::ClustData,
                    opt_data::OptDataCEP,
                    opt_config::Dict{String,Any},
                    optimizer::DataType
                    )
  #Check the consistency of the data provided
  check_opt_data_cep(opt_data)
  #Setup the basic elements
  cep=setup_opt_cep_basic(ts_data, opt_data, opt_config, optimizer, opt_config["optimizer_config"])
  #Setup the basic variables
  setup_opt_cep_basic_variables!(cep, ts_data, opt_data)
  #If lost load costs aren't Inf, setup lost load (LL, SLACK)
  if opt_config["lost_load_cost"]["el"]!=Inf
    setup_opt_cep_lost_load!(cep, ts_data, opt_data)
  end
  #If lost emission costs aren't Inf, setup lost emission (LE)
  if opt_config["lost_emission_cost"]["CO2"]!=Inf
    setup_opt_cep_lost_emission!(cep, ts_data, opt_data)
  end
  #If storage and seasonalstorage, setup seasonal storage (continous)
  if opt_config["storage_in"] && opt_config["storage_out"] && opt_config["storage_e"] && opt_config["seasonalstorage"]
    setup_opt_cep_storage!(cep, ts_data, opt_data)
    setup_opt_cep_seasonalstorage!(cep, ts_data, opt_data)
  #Else if storage, but no seasonalstorage, setup intra-day storage (same level in first and last time step for each period)
  elseif opt_config["storage_in"] && opt_config["storage_out"] && opt_config["storage_e"] && !(opt_config["seasonalstorage"])
    setup_opt_cep_storage!(cep, ts_data, opt_data)
    setup_opt_cep_simplestorage!(cep, ts_data, opt_data)
  end
  #If transmission, setup TRANS and FLOW
  if opt_config["transmission"]
      setup_opt_cep_transmission!(cep, ts_data, opt_data)
  end
  #Setup the electricity generation
  setup_opt_cep_generation_el!(cep, ts_data, opt_data)
  #If co2-limit isn't Inf, limit the total co2 output of the energy system
  if opt_config["co2_limit"]!=Inf
    setup_opt_cep_co2_limit!(cep, ts_data, opt_data; co2_limit=opt_config["co2_limit"],  lost_emission_cost=opt_config["lost_emission_cost"])
  end
  # Setup the energy balences to match the demand
  setup_opt_cep_demand!(cep, ts_data, opt_data; lost_load_cost=opt_config["lost_load_cost"])
  #If fixed_design_variables are provided, fix the installed capacities to them
  if "fixed_design_variables" in keys(opt_config)
    setup_opt_cep_fix_design_variables!(cep, ts_data, opt_data, opt_config["fixed_design_variables"])
  end
  #If existing_infrastructure, add existing infrastructure
  if opt_config["existing_infrastructure"]
      setup_opt_cep_existing_infrastructure!(cep, ts_data, opt_data)
  end
  #If limit_infrastructure, limit the infrastructure expansion
  if opt_config["limit_infrastructure"]
      setup_opt_cep_limit_infrastructure!(cep, ts_data, opt_data)
  end
  #Setup the objective
  setup_opt_cep_objective!(cep, ts_data, opt_data; lost_load_cost=opt_config["lost_load_cost"], lost_emission_cost=opt_config["lost_emission_cost"])
  # solve and return the CEP
  return solve_opt_cep(cep, ts_data, opt_data, opt_config)
end

"""
     run_opt(ts_data::ClustData,opt_data::OptDataCEP,opt_config::Dict{String,Any},fixed_design_variables::Dict{String,Any},optimizer::DataTyple;lost_el_load_cost::Number=Inf,lost_CO2_emission_cost::Number)
This problem runs the operational optimization problem only, with fixed design variables.
provide the fixed design variables and the `opt_config` of the previous step (design run or another opterational run)
Required elements are:
- `ts_data`: The time-series data, which should be be the original time-series data for this operational run. The `keys(ts_data.data)` need to match the `[time_series_name]-[node]`
- `opt_data`: In this case the OptDataCEP that contains information on costs, nodes, techs and for transmission also on lines. - Should be the same as in the design run.
- `opt_config`: This includes all the previous settings for the design optimization problem formulation and ensures that the configuration is the same.
- `fixed_design_variables`: All the design variables that are determined by the previous design run.
- `optimizer`: The used optimizer, which could e.g. be Clp: `using Clp` `optimizer=Clp.Optimizer` or Gurobi: `using Gurobi` `optimizer=Gurobi.Optimizer`.
What you can change in the `opt_config`:
- `lost_el_load_cost`: Number indicating the lost load price/MWh (should be greater than 1e6), give a number lower than Inf for SLACK and LL (Lost Load - a variable for unmet demand by the installed capacities). No SLACK and LL can lead to infeasibilities.
- `lost_CO2_emission_cost`: Number indicating the emission price/kg-CO2 (should be greater than 1e6), give Inf for no LE (Lost Emissions - a variable for emissions that will exceed the limit in order to provide the demand with the installed capacities). No LE can lead to a higher SLACK and LL, as the optimization is not allowed to break any emission limit then.
"""
function run_opt(ts_data::ClustData,
                    opt_data::OptDataCEP,
                    opt_config::Dict{String,Any},
                    fixed_design_variables::Dict{String,Any},
                    optimizer::DataType;
                    lost_el_load_cost::Number=Inf,
                    lost_CO2_emission_cost::Number=Inf)
  # Create dictionary for lost_load_cost of the single elements
  lost_load_cost=Dict{String,Number}("el"=>lost_el_load_cost)
  # Create dictionary for lost_emission_cost of the single elements
  lost_emission_cost=Dict{String,Number}("CO2"=>lost_CO2_emission_cost)
  # Add the fixed_design_variables and new setting for slack costs to the existing config
  set_opt_config_cep!(opt_config;fixed_design_variables=fixed_design_variables, lost_load_cost=lost_load_cost, lost_emission_cost=lost_emission_cost)

  return run_opt(ts_data,opt_data,opt_config,optimizer)
end

"""
     run_opt(ts_data::ClustData,opt_data::OptDataCEP,optimizer::DataTyple;co2_limit::Number=Inf,lost_el_load_cost::Number=Inf,lost_CO2_emission_cost::Number=Inf,existing_infrastructure::Bool=false,limit_infrastructure::Bool=false,storage::String="none",transmission::Bool=false,descriptor::String="",print_flag::Bool=true,optimizer_config::Dict{Symbol,Any}=Dict{Symbol,Any}(),round_sigdigits::Int=9)
Wrapper function for type of optimization problem for the CEP-Problem (NOTE: identifier is the type of `opt_data` - in this case OptDataCEP - so identification as CEP problem).
Required elements are:
- `ts_data`: The time-series data, which could either be the original input data or some aggregated time-series data. The `keys(ts_data.data)` need to match the `[time_series_name]-[node]`
- `opt_data`: The OptDataCEP that contains information on costs, nodes, techs and for transmission also on lines.
- `optimizer`: The used optimizer, which could e.g. be Clp: `using Clp` `optimizer=Clp.Optimizer` or Gurobi: `using Gurobi` `optimizer=Gurobi.Optimizer`.
Options to tweak the model are:
- `co2_limit`: A number limiting the kg.-CO2-eq./MWh (normally in a range from 5-1250 kg-CO2-eq/MWh), give Inf or no kw if unlimited
- `lost_el_load_cost`: Number indicating the lost load price/MWh (should be greater than 1e6), give Inf for no SLACK and LL (Lost Load - a variable for unmet demand by the installed capacities)
- `lost_CO2_emission_cost`: Number indicating the emission price/kg-CO2 (should be greater than 1e6), give Inf for no LE (Lost Emissions - a variable for emissions that will exceed the limit in order to provide the demand with the installed capacities)
- `existing_infrastructure`: true or false to include or exclude existing infrastructure to the model
- `storage`: String "none" for no storage or "simple" to include simple (only intra-day storage) or "seasonal" to include seasonal storage (inter-day)
Optional elements are:
- `descriptor`: String with the name of this paricular model like "kmeans-10-co2-500"
- `print_flag`: Bool to decide if a summary of the Optimization result should be printed.
- `optimizer_config`: Each Symbol and the corresponding value in the Dictionary is passed on to the `with_optimizer` function in addition to the `optimizer`. For Gurobi an example Dictionary could look like `Dict{Symbol,Any}(:Method => 2, :OutputFlag => 0, :Threads => 2)` more information can be found in the optimizer specific documentation.
- `round_sigdigits`: Can be used to round the values of the result to a certain number of `sigdigits`.
"""
function run_opt(ts_data::ClustData,
                 opt_data::OptDataCEP,
                 optimizer::DataType;
                 co2_limit::Number=Inf,
                 lost_el_load_cost::Number=Inf,
                 lost_CO2_emission_cost::Number=Inf,
                 existing_infrastructure::Bool=false,
                 limit_infrastructure::Bool=false,
                 storage::String="none",
                 transmission::Bool=false,
                 descriptor::String="",
                 print_flag::Bool=true,
                 optimizer_config::Dict{Symbol,Any}=Dict{Symbol,Any}(),
                 round_sigdigits::Int=9)
   # Activated seasonal or simple storage corresponds with storage
   if storage=="seasonal"
       storage=true
       seasonalstorage=true
   elseif storage=="simple"
       storage=true
       seasonalstorage=false
   elseif storage =="none"
       storage=false
       seasonalstorage=false
  else
      storage=false
      seasonalstorage=false
      @warn("String indicating storage not identified as 'none', 'seasonal' or 'simple' → no storage")
   end
  # Create dictionary for lost_load_cost of the single elements
  lost_load_cost=Dict{String,Number}("el"=>lost_el_load_cost)
  # Create dictionary for lost_emission_cost of the single elements
  lost_emission_cost=Dict{String,Number}("CO2"=>lost_CO2_emission_cost)

  #Setup the opt_config file based on the data input and
  opt_config=set_opt_config_cep(opt_data; descriptor=descriptor, co2_limit=co2_limit, lost_load_cost=lost_load_cost, lost_emission_cost=lost_emission_cost, existing_infrastructure=existing_infrastructure, limit_infrastructure=limit_infrastructure, storage_e=storage, storage_in=storage, storage_out=storage, seasonalstorage=seasonalstorage, transmission=transmission, print_flag=print_flag, optimizer_config=optimizer_config, round_sigdigits=round_sigdigits)
  #Run the optimization problem
  run_opt(ts_data, opt_data, opt_config, optimizer)
end # run_opt


text:
function f()
    steps1 = 0
    steps2 = 0
    while !isborder(robot, Nord)
        steps1 += 1
        move!(robot, Nord)
    end
    while !isborder(robot, West)
        steps2 += 1
        move!(robot, West)
        putmarker!(robot)
    end
    steps2 = 11 - steps2
    while !isborder(robot, Sud)
        move!(robot, Sud)
        putmarker!(robot)
    end
    while !isborder(robot, Ost)
        move!(robot, Ost)
        putmarker!(robot)
    end
    while !isborder(robot, Nord)
        move!(robot, Nord)
        putmarker!(robot)
    end
    while steps2 > 0
        move!(robot, West)
        putmarker!(robot)
        steps2 -= 1
    end
    while steps1 > 0
        move!(robot, Sud)
        steps1 -= 1
    end
end

f()

max_stars_repo_name: TheXlebushek/mirea-progs

max_issues_repo_name: TheXlebushek/mirea-progs

max_forks_repo_name: TheXlebushek/mirea-progs