text:
Variable(vidx) = Variable(
    vidx,
    0,
    0,
    0,
    0,
    [],
    [],
    [],
    [],
    0,
    0,
    0,
    Vector{Vector{Tuple{Symbol,Int,Int,Int}}}(),
    false,
    false,
    false,
    false,
    0.0,
    Vector{Bool}()
)

MatchingInit() = MatchingInit(0, Int[], Int[], Int[], Int[], Int[], Int[], Bool[], Bool[])
SCCInit() = SCCInit(Int[], Int[], Int[], Bool[], Int[])

function ConstraintInternals(cidx::Int, fct, set, indices::Vector{Int})
    return ConstraintInternals(
        cidx,
        fct,
        set,
        indices,
        Int[],
        ImplementedConstraintFunctions(),
        false,
        false,
        false,
        Vector{BoundRhsVariable}(undef, 0),
    )
end

function ImplementedConstraintFunctions()
    return ImplementedConstraintFunctions([
        false for f in fieldnames(ImplementedConstraintFunctions)
    ]...)
end

function LinearConstraint(
    cidx,
    indices::Vector,
    coeffs::Vector{T},
    constant,
    set::MOI.AbstractScalarSet,
) where {T}
    @assert length(indices) == length(coeffs)
    len = length(coeffs)
    scalar_terms = Vector{MOI.ScalarAffineTerm{T}}(undef, len)
    for (i, idx, coeff) in zip(1:len, indices, coeffs)
        scalar_terms[i] = MOI.ScalarAffineTerm{T}(coeff, MOI.VariableIndex(idx))
    end
    saf = MOI.ScalarAffineFunction(scalar_terms, constant)
    return LinearConstraint(cidx, saf, set, indices)
end

function LinearConstraint(
    cidx::Int,
    fct::MOI.ScalarAffineFunction,
    set::MOI.AbstractScalarSet,
    indices::Vector{Int},
)
    # get common type for rhs and coeffs
    # use the first value (can be .upper, .lower, .value) and subtract left constant
    rhs = -fct.constant
    if isa(set, Union{MOI.EqualTo,CS.NotEqualTo})
        rhs += set.value
    elseif isa(set, MOI.LessThan)
        rhs += set.upper
    end
    coeffs = [t.coefficient for t in fct.terms]
    promote_T = promote_type(typeof(rhs), eltype(coeffs))
    if promote_T != eltype(coeffs)
        coeffs = convert.(promote_T, coeffs)
    end
    if promote_T != typeof(rhs)
        rhs = convert(promote_T, rhs)
    end
    maxs = zeros(promote_T, length(indices))
    mins = zeros(promote_T, length(indices))
    pre_maxs = zeros(promote_T, length(indices))
    pre_mins = zeros(promote_T, length(indices))
    # this can be changed later in `set_in_all_different!` but needs to be initialized with false
    in_all_different = false

    is_strict = !(isa(set, MOI.LessThan) || isa(set, MOI.EqualTo))
    is_equal = isa(set, MOI.EqualTo)
    internals = ConstraintInternals(cidx, fct, set, indices)
    # the rhs is filled in init_constraint
    lc = LinearConstraint(
        internals,
        in_all_different,
        is_strict,
        is_equal,
        false, # the rhs is not a strong rhs yet
        zero(promote_T),
        zero(promote_T),
        coeffs,
        mins,
        maxs,
        pre_mins,
        pre_maxs,
    )
    return lc
end

function ReifiedConstraint(std, act_std, inner_constraint, anti_constraint)
    return ReifiedConstraint(std, act_std, inner_constraint, anti_constraint, false, 0)
end

"""
    ConstraintSolverModel(T::DataType=Float64)

Create the constraint model object and specify the type of the solution
"""
function ConstraintSolverModel(::Type{T} = Float64) where {T<:Real}
    ConstraintSolverModel(
        nothing, # lp_model
        Vector{VariableRef}(), # lp_x
        Vector{Variable}(), # init_search_space
        Vector{Variable}(), # search_space
        Vector{Tuple{Int,Int}}(), # init_fixes
        Vector{Vector{Int}}(), # subscription
        Vector{Constraint}(), # constraints
        Vector{VarAndVal}(), # root_infeasible_vars
        Vector{Int}(), # bt_infeasible
        1, # c_backtrack_idx
        1, # c_step_nr
        Vector{BacktrackObj{T}}(), # backtrack_vec
        PriorityQueue{Int,Priority}(Base.Order.Reverse), # priority queue for `get_next_node`
        MOI.FEASIBILITY_SENSE, #
        NoObjective(), #
        Vector{Bool}(), # var_in_obj
        Val(:DFS),
        get_branch_strategy(),
        get_branch_split(),
        true,
        ActivityObj(),
        zero(T), # best_sol,
        zero(T), # best_bound
        Vector{Solution}(), # all solution objects
        CSInfo(0, false, 0, 0, 0, NumberConstraintTypes()), # info
        Dict{Symbol,Any}(), # input
        Vector{TreeLogNode{T}}(), # logs
        SolverOptions(), # options,
        -1.0, # solve start time
        -1.0, # solve time will be overwritten
    )
end

@deprecate init() ConstraintSolverModel()

function NumberConstraintTypes()
    return NumberConstraintTypes(zeros(Int, length(fieldnames(NumberConstraintTypes)))...)
end

function new_BacktrackObj(com::CS.CoM, parent_idx, vidx, lb, ub)
    parent = com.backtrack_vec[parent_idx]
    return BacktrackObj{parametric_type(com)}(
        length(com.backtrack_vec) + 1, # idx
        -1, # step_nr
        parent_idx,
        parent.depth + 1,
        :Open, # status
        true, # is feasible
        vidx,
        lb, # lb and ub only take effect if vidx != 0
        ub, # ub
        parent.best_bound,
        parent.solution,
        zeros(length(com.search_space)), # solution values of bound computation
    )
end

function BacktrackObj(com::CS.CoM)
    return BacktrackObj(
        1, # idx
        -1, # step_nr
        0, # parent_idx
        0, # depth
        :Closed, # status
        true, # is feasible until proven otherwise
        0, # vidx
        0, # lb and ub only take effect if vidx != 0
        0, # ub
        com.sense == MOI.MIN_SENSE ? typemax(com.best_bound) : typemin(com.best_bound),
        zeros(length(com.search_space)),
        zeros(length(com.search_space)),
    )
end

function BoolConstraintInternals()
    return BoolConstraintInternals(false, 0, false, 0)
end

function AndConstraint(internals::ConstraintInternals, lhs::Constraint, rhs::Constraint)
    return AndConstraint(internals, BoolConstraintInternals(), lhs, rhs)
end

function OrConstraint(internals::ConstraintInternals, lhs::Constraint, rhs::Constraint)
    return OrConstraint(internals, BoolConstraintInternals(), lhs, rhs)
end


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

module GenX

export configure_settings
export configure_solver
export load_inputs
export generate_model
export solve_model
export write_outputs
export cluster_inputs
export mga
export morris
export simple_operation

using JuMP # used for mathematical programming
using DataFrames # This package allows put together data into a matrix
using MathProgBase # for fix_integers
using CSV
using StatsBase
using LinearAlgebra
using YAML
using Dates
using Clustering
using Distances
using Combinatorics
using Documenter
using DiffEqSensitivity
using OrdinaryDiffEq
using QuasiMonteCarlo
using Random
using RecursiveArrayTools
using Statistics
using BenchmarkTools

# Uncomment if Gurobi or CPLEX active license and installations are there and the user intends to use either of them
using Gurobi
#using CPLEX
#using MOI
#using SCIP
# using Clp
# using Cbc
using GLPK

# Global scaling factor used when ParameterScale is on to shift values from MW to GW
# DO NOT CHANGE THIS (Unless you do so very carefully)
# To translate MW to GW, divide by ModelScalingFactor
# To translate $ to $M, multiply by ModelScalingFactor^2
# To translate $/MWh to $M/GWh, multiply by ModelScalingFactor
ModelScalingFactor = 1e+3

# Configure settings
include("configure_settings/configure_settings.jl")

# Configure optimizer instance
include("configure_solver/configure_gurobi.jl")
include("configure_solver/configure_scip.jl")
include("configure_solver/configure_glpk.jl")
include("configure_solver/configure_cplex.jl")
include("configure_solver/configure_clp.jl")
include("configure_solver/configure_cbc.jl")
include("configure_solver/configure_solver.jl")

# Load input data
include("load_inputs/load_generators_data.jl")
include("load_inputs/load_generators_variability.jl")
include("load_inputs/load_network_data.jl")
include("load_inputs/load_reserves.jl")
include("load_inputs/load_cap_reserve_margin.jl")
include("load_inputs/load_energy_share_requirement.jl")
include("load_inputs/load_co2_cap.jl")
include("load_inputs/load_period_map.jl")
include("load_inputs/load_minimum_capacity_requirement.jl")
include("load_inputs/load_load_data.jl")
include("load_inputs/load_fuels_data.jl")

include("load_inputs/load_inputs.jl")

include("time_domain_reduction/time_domain_reduction.jl")

# Core GenX Features
include("model/core/discharge/discharge.jl")
include("model/core/discharge/investment_discharge.jl")

include("model/core/non_served_energy.jl")
include("model/core/ucommit.jl")
include("model/core/emissions.jl")

include("model/core/reserves.jl")

include("model/core/transmission.jl")

include("model/resources/curtailable_variable_renewable/curtailable_variable_renewable.jl")

include("model/resources/flexible_demand/flexible_demand.jl")

include("model/resources/hydro/hydro_res.jl")
include("model/resources/hydro/hydro_inter_period_linkage.jl")

include("model/resources/must_run/must_run.jl")

include("model/resources/storage/storage.jl")
include("model/resources/storage/investment_energy.jl")
include("model/resources/storage/storage_all.jl")
include("model/resources/storage/long_duration_storage.jl")
include("model/resources/storage/investment_charge.jl")
include("model/resources/storage/storage_asymmetric.jl")
include("model/resources/storage/storage_symmetric.jl")

include("model/resources/thermal/thermal.jl")
include("model/resources/thermal/thermal_commit.jl")
include("model/resources/thermal/thermal_no_commit.jl")

include("model/policies/co2_cap.jl")
include("model/policies/energy_share_requirement.jl")
include("model/policies/cap_reserve_margin.jl")
include("model/policies/minimum_capacity_requirement.jl")

include("model/generate_model.jl")
include("model/solve_model.jl")

include("write_outputs/dftranspose.jl")
include("write_outputs/write_capacity.jl")
include("write_outputs/write_charge.jl")
include("write_outputs/write_charging_cost.jl")
include("write_outputs/write_costs.jl")
include("write_outputs/write_curtailment.jl")
include("write_outputs/write_emissions.jl")
include("write_outputs/write_energy_revenue.jl")
include("write_outputs/write_net_revenue.jl")
include("write_outputs/write_nse.jl")
include("write_outputs/write_power.jl")
include("write_outputs/write_power_balance.jl")
include("write_outputs/write_price.jl")
include("write_outputs/write_reliability.jl")
include("write_outputs/write_status.jl")
include("write_outputs/write_storage.jl")
include("write_outputs/write_storagedual.jl")
include("write_outputs/write_subsidy_revenue.jl")
include("write_outputs/write_time_weights.jl")
include("write_outputs/choose_output_dir.jl")

include("write_outputs/capacity_reserve_margin/write_capacity_value.jl")
include("write_outputs/capacity_reserve_margin/write_reserve_margin_revenue.jl")
include("write_outputs/capacity_reserve_margin/write_reserve_margin_w.jl")
include("write_outputs/capacity_reserve_margin/write_reserve_margin.jl")

include("write_outputs/energy_share_requirement/write_esr_prices.jl")
include("write_outputs/energy_share_requirement/write_esr_revenue.jl")

include("write_outputs/long_duration_storage/write_opwrap_lds_dstor.jl")
include("write_outputs/long_duration_storage/write_opwrap_lds_stor_init.jl")

include("write_outputs/reserves/write_reg.jl")
include("write_outputs/reserves/write_rsv.jl")

include("write_outputs/transmission/write_nw_expansion.jl")
include("write_outputs/transmission/write_transmission_flows.jl")
include("write_outputs/transmission/write_transmission_losses.jl")

include("write_outputs/ucommit/write_commit.jl")
include("write_outputs/ucommit/write_shutdown.jl")
include("write_outputs/ucommit/write_start.jl")

include("write_outputs/write_outputs.jl")

include("additional_tools/modeling_to_generate_alternatives.jl")
include("additional_tools/method_of_morris.jl")
end

ext: jl

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
__precompile__()

module RungeKuttaFehlberg

export rkf45_step, rkf45_step!, RKFBuffer, RKFConfig

"""
    RKFConfig(norm, tolerance, max_dt[, safety])

Container type for Fehlberg integration parameters (time adaptation).

* `norm::Function`: a metric for the difference between two estimates of `dx`.
  Must take both estimates as arguments and return a positive (non-zero) number.
  (Default is the `l1`-norm.)
* `tolerance::Float64`: target value for the error function, `dt` will be
  adapted to achieve this value
* `max_dt::Float64`: maximum allowed timestep.
* `safety::Float64`: parameter for the timestep estimator. Must be greater than
  zero and smaller than one. Small values lead to more conservative (smaller)
  timestep estimates which can improve convergence. (The default, 0.9, is almost
  always fine.)
"""
immutable RKFConfig
    norm::Function
    tolerance::Float64
    max_dt::Float64
    safety::Float64
    RKFConfig(norm, tol, max_dt, safety = 0.9) = new(norm, tol, max_dt, safety)
end

"""
    calulate_steps(f, x, t, dt)

Computes and returns the 4th- and 5th-order Runge-Kutta estimates for `dx` in the differential equation `dx / dt = f(x, t)` using Fehlberg's coefficients [1].

[1] E. Fehlberg, Computing **6**, 61 (1970).
"""
function calculate_steps(f::Function, x, t::Float64, dt::Float64)
            k1 = f(x, t)
            k2 = f(x .+ k1 * (dt * 1 / 4),
                   t + dt / 4)
            k3 = f(x .+ k1 * (dt * 3 / 32)
                     .+ k2 * (dt * 9 / 32),
                   t + 3 / 8 * dt)
            k4 = f(x .+ k1 * (dt * 1932 / 2197)
                     .- k2 * (dt * 7200 / 2197)
                     .+ k3 * (dt * 7296 / 2197),
                   t + 12 / 13 * dt)
            k5 = f(x .+ k1 * (dt * 439 / 216)
                     .- k2 * (dt * 8)
                     .+ k3 * (dt * 3680 / 513)
                     .- k4 * (dt * 845 / 4104),
                   t + dt)
            k6 = f(x .- k1 * (dt * 8 / 27)
                     .+ k2 * (dt * 2)
                     .- k3 * (dt * 3544 / 2565)
                     .+ k4 * (dt * 1859 / 4104)
                     .- k5 * (dt * 11 / 40),
                   t + dt / 2)
            step_rk4 = (k1 * (dt * 25 / 216)
                        .+ k3 * (dt * 1408 / 2565)
                        .+ k4 * (dt * 2197 / 4104)
                        .- k5 * (dt / 5))
            step_rk5 = (k1 * (dt * 16 / 135)
                        .+ k3 * (dt * 6656 / 12825)
                        .+ k4 * (dt * 28561 / 56430)
                        .- k5 * (dt * 9 / 50)
                        .+ k6 * (dt * 2 / 55))
            return step_rk4, step_rk5
end

immutable RKFBuffer{T}
    k1::T
    k2::T
    k3::T
    k4::T
    k5::T
    k6::T
    x_temp::T
end
RKFBuffer{T}(x::T) = RKFBuffer(similar(x), similar(x), similar(x), similar(x),
                               similar(x), similar(x), similar(x))

"""
    calulate_steps(f!, x, t, dt, b, dx4, dx5)

Like calculate_steps, but takes a mutating function `f!(x, t, dx)` which writes
the differential `dx` into its final argument.  Must not have any other side
effects.

`x` must be of a mutable type (it will not be mutated) and there must be a
`map!()` method for it.

`b` is a pre-allocated `RKFBuffer` object (it will be mutated).
The 4th and 5th order results will be written to `dx4` and  `dx5`, respectively.
"""
function calculate_steps!{T}(f!::Function, x::T, t::Float64, dt::Float64,
                             b::RKFBuffer{T}, dx4::T, dx5::T)
            f!(x, t, b.k1)
            map!((x, k1) -> x + k1 * 0.35dt, b.x_temp, x, b.k1)
            f!(b.x_temp, t + 0.25dt, b.k2)
            map!(b.x_temp, x, b.k1, b.k2) do x, k1, k2
                x + k1 * (dt * 3 / 32) + k2 * (dt * 9 / 32)
            end
            f!(b.x_temp, t + 0.375dt, b.k3)
            map!(b.x_temp, x, b.k1, b.k2, b.k3) do x, k1, k2, k3
                x + k1 * (dt * 1932 / 2197) - k2 * (dt * 7200 / 2197) +
                    k3 * (dt * 7296 / 2197)
            end
            f!(b.x_temp, t + 12 / 13 * dt, b.k4)
            map!(b.x_temp, x, b.k1, b.k2, b.k3, b.k4) do x, k1, k2, k3, k4
                x + k1 * (dt * 439 / 216) - k2 * (dt * 8) +
                    k3 * (dt * 3680 / 513) - k4 * (dt * 845 / 4104)
            end
            f!(b.x_temp, t + dt, b.k5)
            map!(b.x_temp, x, b.k1, b.k2, b.k3, b.k4, b.k5) do x, k1, k2, k3,
                                                               k4, k5
                x - k1 * (dt * 8 / 27) + k2 * (dt * 2) -
                    k3 * (dt * 3544 / 2565) + k4 * (dt * 1859 / 4104) -
                    k5 * (dt * 11 / 40)
            end
            f!(b.x_temp, t + 0.5dt, b.k6)
            map!(dx4, b.k1, b.k3, b.k4, b.k5) do k1, k3, k4, k5
                k1 * (dt * 25 / 216) + k3 * (dt * 1408 / 2565) +
                    k4 * (dt * 2197 / 4104) - k5 * 0.2dt
            end
            map!(dx5, b.k1, b.k3, b.k4, b.k5, b.k6) do k1, k3, k4, k5, k6
                k1 * (dt * 16 / 135) + k3 * (dt * 6656 / 12825) +
                    k4 * (dt * 28561 / 56430) - k5 * (dt * 9 / 50) +
                    k6 * (dt * 2 / 55)
            end
            return dx4, dx5
end

"""
Default metric for `rkf45_step()`, based on the `l1`-norm (aka 'Manhattan'
distance).
"""
function l1_metric{T}(x1::T, x2::T)
    return vecnorm(x1 .- x2, 1)
end

"""
    rkf45_step(f, x, t, tolerance, dt, config)

Computes the 5th-order Runge-Kutta estimate for `dx` in the differential
equation `dx / dt = f(x, t)` with an adaptive time step method (RKF45).
Returns the estimate for `dx`, the associated time step `dt` and a suggestion
for the next timestep `next_dt`, in that order.

# Arguments

* `f::Function`: the time derivative of `x`, must take two arguments: a state
  `x` and a time `t`. Should return a subtype of `typeof(x)`. Must not have any
  side effects.
* `x`: the state variable
* `t::Float64`: the time associated with `x`
* `dt::Float64`: estimate for an appropriate timestep, which should be the value
  for `dt` which was returned by the previous step. If `f` is well-behaved, `dt`
  will converge from any initial (positive) value after few steps.
* `config::RKFConfig`: container with integration parameters.
"""
function rkf45_step(f::Function,
                    x,
                    t::Float64,
                    dt::Float64,
                    config::RKFConfig)
    step_rk4, step_rk5 = calculate_steps(f, x, t, dt)
    err = config.norm(step_rk4, step_rk5)
    while err > config.tolerance
        dt *= config.safety * (config.tolerance / err)^(1 / 5)
        step_rk4, step_rk5 = calculate_steps(f, x, t, dt)
        err = config.norm(step_rk4, step_rk5)
    end
    next_dt = dt * config.safety * (config.tolerance / err)^(1 / 4)
    # Note that we need to catch cases where `err` approaches zero.
    return step_rk5, dt, min(next_dt, config.max_dt)
end

"""
    rkf45_step!(f!, x, t, dt, config, dx4, dx5, buffer)

Like `rkf45_step`, but instead of dynamical allocations uses `buffer`, `dx4` and
`dx5` to store data.
"""
function rkf45_step!{T}(f!::Function,
                        x::T,
                        t::Float64,
                        dt::Float64,
                        config::RKFConfig,
                        dx4::T,
                        dx5::T,
                        buffer::RKFBuffer{T})
    calculate_steps!(f!, x, t, dt, buffer, dx4, dx5)
    err = config.norm(dx4, dx5)
    while err > config.tolerance
        dt *= config.safety * (config.tolerance / err)^(1 / 5)
        calculate_steps!(f!, x, t, dt, buffer, dx4, dx5)
        err = config.norm(dx4, dx5)
    end
    next_dt = dt * config.safety * (config.tolerance / err)^(1 / 4)
    # Note that we need to catch cases where `err` approaches zero.
    return dx5, dt, min(next_dt, config.max_dt)
end

end # module


text:
module AbstractSphericalHarmonics

using ComputedFieldTypes
using LinearAlgebra
using StaticArrays

################################################################################

# FastSphericalHarmonics requires julia >= 1.7
# SSHT requires julia >= 1.6

# using FastSphericalHarmonics
using SSHT

export ash_grid_size, ash_nmodes
export ash_ntheta, ash_nphi, ash_thetas, ash_phis, ash_point_coord, ash_point_delta, ash_grid_as_phi_theta
export ash_mode_index, ash_mode_numbers
export ash_transform!, ash_transform, ash_evaluate!, ash_evaluate
export ash_eth!, ash_eth, ash_ethbar!, ash_ethbar

################################################################################

bitsign(i::Integer) = isodd(i) ? -one(i) : one(i)

################################################################################

include("sylm.jl")
include("tensor.jl")

end


text:
# ============================================================================ #
# DictSet
# ============================================================================ #
"""
Basic dict-like thing that holds multiple values for each key in a set.
"""
struct DictSet{K, V}
    data::Dict{K, Set{V}}
end

DictSet{K, V}() where {K, V} = DictSet{K, V}(Dict{K, Set{V}}())

function Base.getindex(ds::DictSet{K, V}, k::K) where {K, V}
    if haskey(ds.data, k)
        ds.data[k]
    else
        Set{V}()
    end
end

function Base.insert!(ds::DictSet{K, V}, k::K, v::V) where {K, V}
    if haskey(ds.data, k)
        push!(ds.data[k], v)
    else
        ds.data[k] = Set{V}(v)
    end
    v
end

function Base.delete!(ds::DictSet{K}, k::K) where K
    delete!(ds.data, k)
    ds
end

function Base.delete!(ds::DictSet{K, V}, k::K, v::V) where {K, V}
    if haskey(ds.data, k)
        s = ds.data[k]
        if length(s) > 1
            delete!(s, v)
        else
            delete!(ds.data, k)
        end
    end
    ds
end

function Base.delete!(ds::DictSet{K, V}, vals::AbstractSet{V}) where {K, V}
    for (k, v) in ds.data
        setdiff!(v, vals)
        if isempty(v)
            delete!(ds.data, k)
        end
    end
    ds
end

Base.keys(ds::DictSet) = keys(ds.data)
Base.values(ds::DictSet) = union(values(ds.data)...)
Base.haskey(ds::DictSet) = haskey(ds.data)
Base.isempty(ds::DictSet) = isempty(ds.data)

# ============================================================================ #
"""
  validatesimplex(σ)

Check if a simplex σ is valid. Throw an ArgumentError if not.
A simplex should contain Integers, be sorted and contain no duplicates.
"""
function validatesx(σ)
    if !(eltype(σ) <: Integer)
        throw(ArgumentError("Invalid simplex '$σ'!" *
                            "Labels should be Integers!"))
    end
    if !allunique(σ)
        throw(ArgumentError("Invalid simplex '$σ'!" *
                            "Vertices in simplex should be unique."))
    end
    if !issorted(σ)
        throw(ArgumentError("Invalid simplex '$σ'!" *
                            "The labels in simplex should be sorted."))
    end
    σ
end

"""
  isface(τ, σ)

Check if τ is a face of σ. Assumes the simplices are valid.
"""
function isface(τ, σ)
    i = 1; j = 1
    while i <= length(τ) && j <= length(σ)
        if τ[i] > σ[j]
            j += 1
        elseif τ[i] == σ[j]
            i += 1; j += 1
        else
            return false
        end
    end
    i == length(τ) + 1
end


text:
module Plates

using DocStringExtensions
using LinearAlgebra: rmul!

export Plate, bound_circulation, bound_circulation!,
       enforce_no_flow_through!, vorticity_flux, suction_parameters, unit_impulse, force,
       moment

using ..Points
using ..Blobs

using ..Elements
using ..RigidBodyMotions

import ..Elements: position, impulse, circulation
import ..Motions: induce_velocity, induce_velocity!, mutually_induce_velocity!, self_induce_velocity,
                  self_induce_velocity!, allocate_velocity, advect!

import ..Utils:@get, MappedVector

include("plates/chebyshev.jl")

"""
    Plate <: Elements.Element

An infinitely thin, flat plate, represented as a bound vortex sheet

# Constructors
- `Plate(N, L, c, α)`
"""
mutable struct Plate <: Element
    "chord length"
    L::Float64
    "centroid"
    c::ComplexF64
    "centroid velocity"
    α::Float64

    "total circulation"
    Γ::Float64

    "number of control points"
    N::Int
    "normalized positions (within [-1, 1]) of the control points"
    ss::Vector{Float64}
    "control point coordinates"
    zs::Vector{ComplexF64}
    "Chebyshev coefficients of the normal component of velocity induced along the plate by ambient vorticity"
    A::MappedVector{Float64, Vector{ComplexF64}, typeof(imag)}
    "Chebyshev coefficients of the velocity induced along the plate by ambient vorticity"
    C::Vector{ComplexF64}
    "zeroth Chebyshev coefficient associated with body motion"
    B₀::Float64
    "first Chebyshev coefficient associated with body motion"
    B₁::Float64

    "Preplanned discrete Chebyshev transform"
    dchebt!::Chebyshev.Transform{ComplexF64, true}
end
@kind Plate Singleton

#function deserialize(s::AbstractSerializer, t::Type{Plate})
#    fields = Tuple(deserialize(s) for i in 1:11)
#    deserialize(s)
#    dchebt! = Chebyshev.plan_transform!(fields[9])
#    Plate(fields..., dchebt!)
#end

function Plate(N, L, c, α)
    ss = Chebyshev.nodes(N)
    zs = c .+ 0.5L*ss*exp(im*α)

    C  = zeros(ComplexF64, N)
    A = MappedVector(imag, C, 1)

    dchebt! = Chebyshev.plan_transform!(C)

    Plate(L, c, α, 0.0, N, ss, zs, A, C, 0.0, 0.0, dchebt!)
end

Base.length(p::Plate) = p.N

circulation(p::Plate) = p.Γ
function impulse(p::Plate)
    @get p (c, B₀, α, Γ, L, A)
    -im*c*Γ - exp(im*α)*π*(0.5L)^2*im*(A[0] - 0.5A[2] - B₀)
end
function angularimpulse(p::Plate)
    @get p (c, B₀, α, Γ, L, A)
    -0.5*(c*conj(c)+L^2/8)*Γ - real(exp(im*α)*conj(c)*π*(0.5L)^2*(A[0] - 0.5A[2] - B₀)) -
          π*(0.25L)^3*(A[1]-A[3]-B₁)
end

function allocate_velocity(::Plate)
    @warn("Plate kinematics should be initialized manually.  This simply returns a stationary motion")
    RigidBodyMotion(0.0, 0.0)
end
function self_induce_velocity!(motion, ::Plate, t)
    motion.ċ, motion.c̈, motion.α̇ = motion.kin(t)
    motion
end

normal(z, α) = imag(exp(-im*α)*z)
tangent(z, α) = real(exp(-im*α)*z)

function induce_velocity(z::ComplexF64, p::Plate, t)
    @get p (α, L, c, B₀, B₁, Γ, A)

    z̃ = conj(2*(z - c)*exp(-im*α)/L)

    ρ = √(z̃ - 1)*√(z̃ + 1)
    J = z̃ - ρ

    w = (A[1] + 2Γ/(π*L))
    w += 2(A[0] - B₀)*J
    w -= B₁*J^2

    w /= ρ

    Jⁿ = J
    for n in 1:length(A)-1
        w -= 2A[n]*Jⁿ
        Jⁿ *= J
    end

      # note that this provides u+iv (JDE)
    0.5im*w*exp(im*α)
end

function induce_velocity!(ws::Vector, p::Plate, sources::T, t) where T <: Union{Tuple, AbstractArray}
    for source in sources
        induce_velocity!(ws, p, source, t)
    end
    ws
end
function induce_velocity(p::Plate, src, t)
    out = allocate_velocity(p.zs)
    induce_velocity!(out, p, src, t)
end

function induce_velocity!(ws::Vector, p::Plate, src, t)
    _singular_velocity!(ws, p, Elements.unwrap(src), t,
                        kind(Elements.unwrap_src(src)))
end

function _singular_velocity!(ws, p, src::Blob{T}, t, ::Type{Singleton}) where T
    induce_velocity!(ws, p.zs, Point{T}(src.z, src.S), t)
end

function _singular_velocity!(ws, p, src, t, ::Type{Singleton})
    induce_velocity!(ws, p.zs, src, t)
end

function _singular_velocity!(ws, p, src, t, ::Type{Group})
    for i in eachindex(src)
        induce_velocity!(ws, p, src[i], t)
    end
    ws
end

induce_velocity!(m::RigidBodyMotion, target::Plate, source, t) = m

function Elements.streamfunction(z::ComplexF64, p::Plate)
    @get p (N, L, c, α, Γ, A, B₀, B₁)
    z̃ = 2*(z - c)*exp(-im*α)/L
    J = z̃ - √(z̃ - 1)*√(z̃ + 1)

    ψ = -log(abs(J))*(2Γ/(π*L))
    ψ -= real(J)*2(A[0] - B₀)
    ψ -= real(0.5J^2)*(A[1] - B₁)

    for n in 2:N-1
        ψ -= A[n]*real(J^(n+1)/(n+1) - J^(n-1)/(n-1))
    end

    return -ψ*L/4
end

function Elements.complexpotential(z::ComplexF64, p::Plate)
    @get p (N, L, c, α, Γ, A, B₀, B₁)
    z̃ = 2*(z - c)*exp(-im*α)/L
    J = z̃ - √(z̃ - 1)*√(z̃ + 1)

    ψ = log(J)*(2Γ/(π*L))
    ψ += J*2(A[0] - B₀)
    ψ += 0.5J^2*(A[1] - B₁)

    for n in 2:N-1
        ψ += A[n]*(J^(n+1)/(n+1) - J^(n-1)/(n-1))
    end

    return im*ψ*L/4
end

function advect!(plate₊::Plate, plate₋::Plate, ṗ::RigidBodyMotion, Δt)
    if plate₊ != plate₋
        plate₊.L    = plate₋.L
        plate₊.Γ    = plate₋.Γ
        plate₊.B₀   = plate₋.B₀
        plate₊.B₁   = plate₋.B₁
        if plate₊.N != plate₋.N
            plate₊.N    = plate₋.N
            resize!(plate₊.ss, plate₊.N)
            resize!(plate₊.zs, plate₊.N)
            resize!(plate₊.C, plate₊.N)
        end
        plate₊.ss   .= plate₋.ss
        plate₊.zs   .= plate₋.zs
        plate₊.C    .= plate₋.C
        plate₊.dchebt! = plate₋.dchebt!
    end
    plate₊.c = plate₋.c + ṗ.ċ*Δt
    plate₊.α = plate₋.α + ṗ.α̇*Δt

    @get plate₊ (c, L, α)

    @. plate₊.zs = c + 0.5L*exp(im*α)*plate₊.ss

    plate₊
end

"""
    unit_impulse(src, plate::Plate)

Compute the impulse per unit circulation of `src` and its associated bound vortex sheet on `plate` (its image vortex)
`src` can be either a `ComplexF64` or a subtype of `Vortex.PointSource`.
"""
function unit_impulse(z::ComplexF64, plate::Plate)
    z̃ = 2(z - plate.c)*exp(-im*plate.α)/plate.L
    unit_impulse(z̃)
end
unit_impulse(z̃) = -im*(z̃ + real(√(z̃ - 1)*√(z̃ + 1) - z̃))
unit_impulse(src, plate::Plate) = unit_impulse(Elements.position(src), plate)

include("plates/boundary_conditions.jl")
include("plates/circulation.jl")
include("plates/force.jl")

@doc raw"""
    surface_pressure(plate, motion, te_sys, Γs₋, Δt)

Compute the pressure difference across the plate along Chebyshev nodes.

!!! note
    The pressure difference across the bound vortex sheet is given by:
    ```math
        [p]_-^+
      = -\rho \left[ \frac{1}{2}(\boldsymbol{v}^+ + \boldsymbol{v}^-)
                   - \boldsymbol{v}_b
             \right]
             \cdot ( \boldsymbol{\gamma} \times \boldsymbol{\hat{n}})
        +\rho \frac{\mathrm{d}\Gamma}{\mathrm{d}t}
    ```
    where ``\rho`` is the fluid density, ``\boldsymbol{v}^\pm`` is the
    velocity on either side of the plate, ``\boldsymbol{v}_b`` is the local
    velocity of the plate, ``\boldsymbol{\gamma}`` is the bound vortex
    sheet strength, and ``\Gamma`` is the integrated circulation.
    We will compute ``\frac{\mathrm{d}\Gamma}{\mathrm{d}t}`` using finite
    differences.  So we will need the circulation along the plate from a
    previous time-step in order to compute the current pressure
    distribution.  We assume that value of circulation at the trailing
    edge of the plate is equal the the net circulation of all the vorticity
    that has been shed from the trailing edge.

# Arguments

- `plate`: we assume that the `Plate` structure that is passed in
  already enforces the no-flow-through condition
- `motion`: the motion of the plate used to compute ``\boldsymbol{v}_b``
- `te_sys`: the system of vortex elements representing the vorticity
  shed from the trailing edge of the plate
- `Γs₋`: the circulation along the plate's Chebyshev nodes, this
  should be equivalent to calling
  `Vortex.circulation(te_sys) .+ Vortex.bound_circulation(plate)`
  from a previous time-step.
- `Δt`: time-step used to compute ``\frac{\mathrm{d}\Gamma}{\mathrm{d}t}``
  using finite differences

# Returns

- `Δp`: the pressure difference across the plate along Chebyshev nodes
- `Γs₊`: the circulation along the plate at the current time-step
  (this value is used in computing the current `Δp` and can be used as
  the `Γs₋` for computing pressure differences at the **next** time-step)
"""
function surface_pressure(plate, motion, ambient_sys, Γs₋, Δt)
    @get plate (C, ss, α)

    Δp = strength(plate) .* (Chebyshev.firstkind(real.(C), ss) .- tangent(motion.ċ, α))

    Γs₊ = circulation(ambient_sys) .+ bound_circulation(plate)
    Δp .+= (Γs₊ .- Γs₋)./Δt

    Δp, Γs₊
end

"""
    edges(plate)

Return the coordinates of the leading and trailing edges

# Example

```jldoctest
julia> p = Plate(128, 1.0, 0, π/4)
Plate: N = 128, L = 1.0, c = 0.0 + 0.0im, α = 45.0ᵒ
       LESP = 0.0, TESP = 0.0

julia> Plates.edges(p)
(0.3535533905932738 + 0.35355339059327373im, -0.3535533905932738 - 0.35355339059327373im)
```
"""
edges(plate) = plate.zs[end], plate.zs[1]

include("plates/pressure.jl")

function Base.show(io::IO, p::Plate)
    lesp, tesp = suction_parameters(p)
    println(io, "Plate: N = $(p.N), L = $(p.L), c = $(p.c), α = $(round(rad2deg(p.α); digits=2))ᵒ")
    print(io, "       LESP = $(round(lesp; digits=2)), TESP = $(round(tesp; digits=2))")
end

end


text:
using Statistics
using JLD2
using Printf
using Oceananigans
using Oceananigans.Units

using Oceananigans.MultiRegion
using Oceananigans.MultiRegion: multi_region_object_from_array
using Oceananigans.Fields: interpolate, Field
using Oceananigans.Architectures: arch_array
using Oceananigans.Coriolis: HydrostaticSphericalCoriolis
using Oceananigans.BoundaryConditions
using Oceananigans.Grids: boundary_node, inactive_node, peripheral_node
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBoundary
using CUDA: @allowscalar, device!
using Oceananigans.Operators
using Oceananigans.Operators: Δzᵃᵃᶜ
using Oceananigans: prognostic_fields

@inline function visualize(field, lev, dims)
    (dims == 1) && (idx = (lev, :, :))
    (dims == 2) && (idx = (:, lev, :))
    (dims == 3) && (idx = (:, :, lev))

    r = deepcopy(Array(interior(field)))[idx...]
    r[ r.==0 ] .= NaN
    return r
end

device!(2)

#####
##### Grid
#####

arch = GPU()
reference_density = 1029

latitude = (-75, 75)

# 0.25 degree resolution
Nx = 1440
Ny = 600

const Nyears  = 2
const Nmonths = 12
const thirty_days = 30days

output_prefix = "near_global_shallow_water_$(Nx)_$(Ny)"
pickup_file   = false

#####
##### Load forcing files and inital conditions from ECCO version 4
##### https://ecco.jpl.nasa.gov/drive/files
##### Bathymetry is interpolated from ETOPO1 https://www.ngdc.noaa.gov/mgg/global/
#####

using DataDeps

path = "https://github.com/CliMA/OceananigansArtifacts.jl/raw/ss/new_hydrostatic_data_after_cleared_bugs/quarter_degree_near_global_input_data/"

datanames = ["tau_x-1440x600-latitude-75",
             "tau_y-1440x600-latitude-75"]

dh = DataDep("quarter_degree_near_global_lat_lon",
    "Forcing data for global latitude longitude simulation",
    [path * data * ".jld2" for data in datanames]
)

DataDeps.register(dh)

datadep"quarter_degree_near_global_lat_lon"

files = [:file_tau_x, :file_tau_y]
for (data, file) in zip(datanames, files)
    datadep_path = @datadep_str "quarter_degree_near_global_lat_lon/" * data * ".jld2"
    @eval $file = jldopen($datadep_path)
end

τˣ = zeros(Nx, Ny, Nmonths)
τʸ = zeros(Nx, Ny, Nmonths)

# Files contain 1 year (1992) of 12 monthly averages
τˣ = file_tau_x["field"] ./ reference_density
τʸ = file_tau_y["field"] ./ reference_density
τˣ = arch_array(arch, τˣ)
τʸ = arch_array(arch, τʸ)

smoothed_bathymetry = jldopen("smooth-bathymetry.jld2")
bat = smoothed_bathymetry["bathymetry"]

bat[ bat .> 0 ] .= 1000.0
bat .-= 10.0

boundary = Int.(bat .> 0)
bat = -bat

# A spherical domain
@show underlying_grid = LatitudeLongitudeGrid(arch,
                                              size = (Nx, Ny),
                                              longitude = (-180, 180),
                                              latitude = latitude,
                                              halo = (4, 4),
                                              topology = (Periodic, Bounded, Flat),
                                              precompute_metrics = true)

grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBoundary(boundary))

#####
##### Boundary conditions / time-dependent fluxes
#####

@inline current_time_index(time, tot_months)     = mod(unsafe_trunc(Int32, time / thirty_days), tot_months) + 1
@inline next_time_index(time, tot_months)        = mod(unsafe_trunc(Int32, time / thirty_days) + 1, tot_months) + 1
@inline cyclic_interpolate(u₁::Number, u₂, time) = u₁ + mod(time / thirty_days, 1) * (u₂ - u₁)

using Oceananigans.Operators: ℑxᶠᵃᵃ, ℑyᵃᶠᵃ

@inline function boundary_stress_u(i, j, k, grid, clock, fields, p)
    time = clock.time
    n₁ = current_time_index(time, Nmonths)
    n₂ = next_time_index(time, Nmonths)

    @inbounds begin
        τ₁ = p.τ[i, j, n₁]
        τ₂ = p.τ[i, j, n₂]
    end

    h_int = ℑxᶠᵃᵃ(i, j, k, grid, fields.h)
    if h_int > 0
        return (cyclic_interpolate(τ₁, τ₂, time) - p.μ * fields.u[i, j, k]) / h_int
    else
        return 0.0
    end
end

@inline function boundary_stress_v(i, j, k, grid, clock, fields, p)
    time = clock.time
    n₁ = current_time_index(time, Nmonths)
    n₂ = next_time_index(time, Nmonths)

    @inbounds begin
        τ₁ = p.τ[i, j, n₁]
        τ₂ = p.τ[i, j, n₂]
    end

    h_int =  ℑyᵃᶠᵃ(i, j, k, grid, fields.h)
    if h_int > 0
        return (cyclic_interpolate(τ₁, τ₂, time) - p.μ * fields.v[i, j, k]) / h_int
    else
        return 0.0
    end
end

# Linear bottom drag:
μ = 0.001 # ms⁻¹

Fu = Forcing(boundary_stress_u, discrete_form = true, parameters = (μ = μ, τ = τˣ))
Fv = Forcing(boundary_stress_v, discrete_form = true, parameters = (μ = μ, τ = τʸ))

using Oceananigans.Models.ShallowWaterModels: VectorInvariantFormulation
using Oceananigans.Advection: VelocityStencil, VorticityStencil
using Oceananigans.TurbulenceClosures: HorizontalDivergenceFormulation

νh = 0e+1

using Oceananigans.Operators: Δx, Δy
using Oceananigans.TurbulenceClosures

@inline νhb(i, j, k, grid, lx, ly, lz) = (1 / (1 / Δx(i, j, k, grid, lx, ly, lz)^2 + 1 / Δy(i, j, k, grid, lx, ly, lz)^2 ))^2 / 5days

horizontal_diffusivity = HorizontalScalarDiffusivity(ν=νh)
biharmonic_viscosity   = HorizontalScalarBiharmonicDiffusivity(ν=νhb, discrete_form=true)

model = ShallowWaterModel(grid = grid,
			              gravitational_acceleration = 9.8055,
                          momentum_advection = WENO5(vector_invariant = VorticityStencil()),
                          mass_advection = WENO5(bounds = (0, Inf)),
                          bathymetry = bat,
                          coriolis = HydrostaticSphericalCoriolis(),
                          forcing = (u=Fu, v=Fv),
			              formulation = VectorInvariantFormulation())

#####
##### Initial condition:
#####

h_init = bat
set!(model, h=h_init)
fill_halo_regions!(model.solution.h)

@info "model initialized"

#####
##### Simulation setup
#####

Δt = 20seconds

simulation = Simulation(model, Δt = Δt, stop_time = Nyears*years)

start_time = [time_ns()]

function progress(sim)
    wall_time = (time_ns() - start_time[1]) * 1e-9

    u = sim.model.solution.u
    h = sim.model.solution.h

    @info @sprintf("Time: % 12s, iteration: %d, max(|u|): %.2e ms⁻¹, min(h): %.2e ms⁻¹, wall time: %s",
                    prettytime(sim.model.clock.time),
                    sim.model.clock.iteration, maximum(abs, u), minimum(h),
                    prettytime(wall_time))

    start_time[1] = time_ns()

    return nothing
end

simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))

u, v, h = model.solution

ζ_op = KernelFunctionOperation{Face, Face, Center}(ζ₃ᶠᶠᶜ, grid; computed_dependencies=(u, v));
ζ = Field(ζ_op)
compute!(ζ)

save_interval = 1days

simulation.output_writers[:surface_fields] = JLD2OutputWriter(model, (; u, v, h, ζ),
                                                            schedule = TimeInterval(save_interval),
                                                            filename = output_prefix * "_surface",
                                                            overwrite_existing = true)

# Let's go!
@info "Running with Δt = $(prettytime(simulation.Δt))"
run!(simulation)


text:
function NewtonianAttraction(arg0::jdouble)
    return NewtonianAttraction((jdouble,), arg0)
end

function acceleration(obj::NewtonianAttraction, arg0::FieldSpacecraftState, arg1::Vector{RealFieldElement})
    return jcall(obj, "acceleration", FieldVector3D, (FieldSpacecraftState, Vector{RealFieldElement}), arg0, arg1)
end

function acceleration(obj::NewtonianAttraction, arg0::SpacecraftState, arg1::Vector{jdouble})
    return jcall(obj, "acceleration", Vector3D, (SpacecraftState, Vector{jdouble}), arg0, arg1)
end

function add_contribution(obj::NewtonianAttraction, arg0::FieldSpacecraftState, arg1::FieldTimeDerivativesEquations)
    return jcall(obj, "addContribution", void, (FieldSpacecraftState, FieldTimeDerivativesEquations), arg0, arg1)
end

function add_contribution(obj::NewtonianAttraction, arg0::SpacecraftState, arg1::TimeDerivativesEquations)
    return jcall(obj, "addContribution", void, (SpacecraftState, TimeDerivativesEquations), arg0, arg1)
end

function depends_on_position_only(obj::NewtonianAttraction)
    return jcall(obj, "dependsOnPositionOnly", jboolean, ())
end

function get_events_detectors(obj::NewtonianAttraction)
    return jcall(obj, "getEventsDetectors", Stream, ())
end

function get_field_events_detectors(obj::NewtonianAttraction, arg0::Field)
    return jcall(obj, "getFieldEventsDetectors", Stream, (Field,), arg0)
end

function get_mu(obj::NewtonianAttraction)
    return jcall(obj, "getMu", jdouble, ())
end

function get_mu(obj::NewtonianAttraction, arg0::Field)
    return jcall(obj, "getMu", RealFieldElement, (Field,), arg0)
end

function get_parameters_drivers(obj::NewtonianAttraction)
    return jcall(obj, "getParametersDrivers", Vector{ParameterDriver}, ())
end


text: using MKtest, CSV, DataFrames, RCall, Plots, StatsBase


text:
using MLLibrary
using RunTests
using Base.Test
module RegressionTest

end

RunTests.set_show_progress(false)

println("Running Tests")
println("...")

tests = ["unit_tests"]

for t in tests
  run_tests(t)
end

println("...")
println("Tests Done")

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
import Flux
import StatsBase
using Random: randn!

export CircularSquaredExponentialKernel

"""
    CircularKernel

An abstract covariance kernel defined over a circle, or more generally a
``d``-dimensional torus.
"""
abstract type CircularKernel <: EuclideanKernel end

function hyperprior_logpdf(self::CircularKernel)
  self.hyperprior.log_variance(self.log_variance) .+
    self.hyperprior.log_length_scales(self.log_length_scales)
end

Flux.trainable(k::CircularKernel) = (k.log_length_scales,k.log_variance)

"""
    CircularSquaredExponentialKernel

A squared exponential kernel on the circle, defined on the circle or, if
``d > 1``, the ``d``-dimensional torus. The expression is given by

``k(x,x') = \\sum_{n \\in 2\\pi\\mathbb{Z}^d} \\sigma^2\\exp\\left(
           -\\left\\lVert\\frac{x - x' + n}{\\kappa} \\right\\rVert^2 \\right)``

parameterized by log-variance ``\\ln(\\sigma^2)`` (trainable by default) and
log-length-scales ``\\ln(\\kappa)`` (trainable by default) applied
element-wise to each dimension. The sum is computed to a specified truncation
level.
"""
mutable struct CircularSquaredExponentialKernel{V<:AbstractVector,H<:Hyperprior} <: CircularKernel
  dims              :: Tuple{Integer, Integer}
  log_variance      :: V
  log_length_scales :: V
  truncation_level  :: Int
  reference_length_scale :: Rational{Int}
  hyperprior        :: NamedTuple{(:log_variance, :log_length_scales), Tuple{H,H}}
end

Flux.@functor CircularSquaredExponentialKernel

"""
    CircularSquaredExponentialKernel(dim::Int)

Creates a squared exponential kernel defined on the circle, or more generally a
torus of dimension `dim`.
"""
function CircularSquaredExponentialKernel(dim::Int)
  if dim!=1 error("Not implemented") end
  dims = (dim, 1)
  log_variance = [0.0]
  log_length_scales = zeros(dim)
  truncation_level = 10
  reference_length_scale = 1//1
  hyperprior = (log_variance = NormalHyperprior([0.0],[1.0]), log_length_scales = NormalHyperprior(zeros(dim),ones(dim)))
  CircularSquaredExponentialKernel(dims, log_variance, log_length_scales, truncation_level, reference_length_scale, hyperprior)
end

"""
    (k::CircularSquaredExponentialKernel)(x1::AbstractMatrix,
                                          x2::AbstractMatrix)

Computes the kernel matrix for the given circular squared exponential kernel.
"""
function (k::CircularSquaredExponentialKernel)(x1::AbstractMatrix, x2::AbstractMatrix)
  (_,m) = size(x1)
  (_,n) = size(x2)
  (d,_) = k.dims
  Fl = eltype(x1)

  loop = Fl(2*pi) .* (-k.truncation_level:1:k.truncation_level) # HACK: only correct for d=1!
  dist = (reshape(pairwise_column_difference(x1,x2), (1,d,m,n)) .+ loop) ./ reshape(exp.(k.log_length_scales), (1,d,1,1))
  sq_dist = dropdims(sum(dist.^2; dims=2); dims=2)
  kernel = exp.(k.log_variance) .* dropdims(sum(exp.(.-sq_dist); dims=1); dims=1)

  m == n ? (kernel + kernel')./2 : kernel # symmetrize to account for roundoff error
end

function (k::LeftGradientKernel{<:CircularSquaredExponentialKernel})(x1::AbstractMatrix, x2::AbstractMatrix)
  (_,m) = size(x1)
  (_,n) = size(x2)
  (d,_) = k.parent.dims
  Fl = eltype(x1)

  loop = Fl(2*pi) .* (-k.parent.truncation_level:1:k.parent.truncation_level) # HACK: only correct for d=1!
  dist = (reshape(pairwise_column_difference(x1,x2), (1,d,m,n)) .+ loop) ./ reshape(exp.(k.parent.log_length_scales), (1,d,1,1))
  sq_dist = dropdims(sum(dist.^2; dims=2); dims=2)
  kernel = exp.(k.parent.log_variance) .* exp.(.-sq_dist)

  dist_sc = dist ./ reshape(exp.(k.parent.log_length_scales), (1,d,1,1))
  @reduce out[D,M,N] := sum(L) -2 * kernel[L,M,N] * dist_sc[L,D,M,N]
  reshape(out, (d*m,n))
end

function (k::RightGradientKernel{<:CircularSquaredExponentialKernel})(x1::AbstractMatrix, x2::AbstractMatrix)
  (_,m) = size(x1)
  (_,n) = size(x2)
  (d,_) = k.parent.dims
  Fl = eltype(x1)

  loop = Fl(2*pi) .* (-k.parent.truncation_level:1:k.parent.truncation_level) # HACK: only correct for d=1!
  dist = (reshape(pairwise_column_difference(x1,x2), (1,d,m,n)) .+ loop) ./ reshape(exp.(k.parent.log_length_scales), (1,d,1,1))
  sq_dist = dropdims(sum(dist.^2; dims=2); dims=2)
  kernel = exp.(k.parent.log_variance) .* exp.(.-sq_dist)

  dist_sc = dist ./ reshape(exp.(k.parent.log_length_scales), (1,d,1,1))
  @reduce out[M,D,N] := sum(L) 2 * kernel[L,M,N] * dist_sc[L,D,M,N]
  reshape(out, (m,d*n))
end

function (k::GradientKernel{<:CircularSquaredExponentialKernel})(x1::AbstractMatrix, x2::AbstractMatrix)
  (_,m) = size(x1)
  (_,n) = size(x2)
  (d,_) = k.parent.dims
  Fl = eltype(x1)

  loop = Fl(2*pi) .* (-k.parent.truncation_level:1:k.parent.truncation_level) # HACK: only correct for d=1!
  dist = (reshape(pairwise_column_difference(x1,x2), (1,d,m,n)) .+ loop) ./ reshape(exp.(k.parent.log_length_scales), (1,d,1,1))
  sq_dist = dropdims(sum(dist.^2; dims=2); dims=2)
  kernel = exp.(k.parent.log_variance) .* exp.(.-sq_dist)

  dist_sc = dist ./ reshape(exp.(k.parent.log_length_scales), (1,d,1,1))
  scale_matrix = Diagonal(exp.(k.parent.log_length_scales))
  @reduce out[D1,M,D2,N] := sum(L) (2*scale_matrix[D1,D2] - 4*dist_sc[L,D1,M,N] * dist_sc[L,D2,M,N]) * kernel[L,M,N]
  reshape(out, (d*m,d*n))
end

"""
    spectral_distribution(k::CircularSquaredExponentialKernel,
                          num_samples::Integer)

Draws `n` samples from the spectral distribution of a standard squared
exponential kernel on the circle.
"""
function spectral_distribution(k::CircularSquaredExponentialKernel; num_samples::Integer)
  (id,_) = k.dims
  support = (-k.truncation_level):(k.truncation_level) |> collect
  measure = [exp(-1//4 * k.reference_length_scale^2 * n^2) for n in support] |> StatsBase.Weights
  frequency = StatsBase.sample(support, measure, num_samples*id)
  reshape(frequency, (id,1,num_samples))
end

function spectral_weights(k::CircularSquaredExponentialKernel, frequency::AbstractArray{<:Any,3})
  Fl = eltype(frequency)
  loop = -k.truncation_level:1:k.truncation_level # HACK: only correct for d=1!
  norm_const_ls = sum(exp.(-1//4 .* reshape(loop, (1,:)).^2 .* exp.(2 .* k.log_length_scales)); dims=(1,2))
  norm_const_std = sum(exp.(-1//4 .* loop.^2 .* k.reference_length_scale.^2); dims=1)
  ratio = dropdims(sum(exp.(-1//8 .* (exp.(2 .* k.log_length_scales) .- k.reference_length_scale.^2) .* frequency.^2); dims=(1,2)); dims=(1,2))
  (exp.(k.log_variance ./ 2) .* ratio .* sqrt.(norm_const_std ./ norm_const_ls), 1)
end


text:
using OrdinaryDiffEq, Test

jac_called = Ref(false)
tgrad_called = Ref(false)

function Lotka(du,u,p,t)
  du[1] = u[1] - u[1] * u[2] # REPL[7], line 3:
  du[2] = -3 * u[2] + 1 * u[1] * u[2]
  nothing
end

function Lotka_jac(J,u,p,t)
  jac_called.x = true
  J[1,1] = 1.0 - u[2]
  J[1,2] = -u[1]
  J[2,1] = 1 * u[2]
  J[2,2] = -3 + u[1]
  nothing
end

function Lotka_tgrad(grad,u,p,t)
  tgrad_called.x = true
  grad[1] = 1 * 0
  grad[2] = 1 * 0
  nothing
end

Lotka_f = ODEFunction(Lotka,jac=Lotka_jac,tgrad=Lotka_tgrad)
prob = ODEProblem(Lotka_f,ones(2),(0.0,10.0))

good_sol = solve(prob,Rosenbrock23())

@test jac_called[]
@test tgrad_called[]

prob2 = ODEProblem(Lotka,ones(2),(0.0,10.0))

sol = solve(prob2,Rosenbrock23(autodiff=true))
@test ≈(good_sol[end],sol[end],rtol=1e-2)

sol = solve(prob2,Rosenbrock23(autodiff=true,chunk_size=1))
@test ≈(good_sol[end],sol[end],rtol=1e-2)

sol = solve(prob2,Rosenbrock23(autodiff=false))
@test ≈(good_sol[end],sol[end],rtol=1e-2)

invW_called = Ref(false)

function Lotka_invW(iW,u,p,γ,t)
  invW_called.x = true
  iW[1,1] = (1 - (1 * u[1] * u[2] * γ ^ 2) / (((1 - (-3 + u[1]) * γ) + (1 * u[1] * u[2] * γ ^ 2) / (1 - (1 - 1 * u[2]) * γ)) * (1 - (1 - 1 * u[2]) * γ))) / (1 - (1 - 1 * u[2]) * γ)
  iW[1,2] = (-(1) * u[1] * γ) / (((1 - (-3 + u[1]) * γ) + (1 * u[1] * u[2] * γ ^ 2) / (1 - (1 - 1 * u[2]) * γ)) * (1 - (1 - 1 * u[2]) * γ))
  iW[2,1] = (u[2] * γ) / (((1 - (-3 + u[1]) * γ) + (1 * u[1] * u[2] * γ ^ 2) / (1 - (1 - 1 * u[2]) * γ)) * (1 - (1 - 1 * u[2]) * γ))
  iW[2,2] = ((1 - (-3 + u[1]) * γ) + (1 * u[1] * u[2] * γ ^ 2) / (1 - (1 - 1 * u[2]) * γ)) ^ -1
  nothing
end

Lotka_f2 = ODEFunction(Lotka,jac=Lotka_jac,tgrad=Lotka_tgrad,invW = Lotka_invW)
prob3 = ODEProblem(Lotka_f2,ones(2),(0.0,10.0))
inv_sol = solve(prob3,Rosenbrock23())

@test invW_called[]
@test ≈(good_sol[end],inv_sol[end],rtol=1e-2)


text:
 innertexelRes = 100
 PETminValue =200.0
 PETmaxValue= 8000.0
 PETValueRange= 8000.0-200
 texelRes=  clamp(innertexelRes, PETminValue,PETmaxValue );
 normalized = (texelRes/Float32(PETValueRange))*6;
 indexx = Int64(floor(normalized)) #// so we normalize floor  in order to get index of color from color list
 colorFloats = [0.0,0.0,1.0,1.0,1.0,1.0]  ;
 normalizedColorPercent= normalized-Float32(indexx) ;
xx=  clamp(colorFloats[indexx]*normalizedColorPercent  +  colorFloats[indexx-1]*(1- normalizedColorPercent) ,0.0,1.0 )

# ######### Defining Helper Functions and imports

# #I use Simple ITK as most robust
# using MedEye3d, Conda,PyCall,Pkg

# Conda.pip_interop(true)
# Conda.pip("install", "SimpleITK")
# Conda.pip("install", "h5py")
# sitk = pyimport("SimpleITK")
# np= pyimport("numpy")

# import MedEye3d
# import MedEye3d.ForDisplayStructs
# import MedEye3d.ForDisplayStructs.TextureSpec
# using ColorTypes
# import MedEye3d.SegmentationDisplay

# import MedEye3d.DataStructs.ThreeDimRawDat
# import MedEye3d.DataStructs.DataToScrollDims
# import MedEye3d.DataStructs.FullScrollableDat
# import MedEye3d.ForDisplayStructs.KeyboardStruct
# import MedEye3d.ForDisplayStructs.MouseStruct
# import MedEye3d.ForDisplayStructs.ActorWithOpenGlObjects
# import MedEye3d.OpenGLDisplayUtils

# """
# given directory (dirString) to file/files it will return the simple ITK image for futher processing
# isMHD when true - data in form of folder with dicom files
# isMHD  when true - we deal with MHD data
# """
# function getImageFromDirectory(dirString,isMHD::Bool, isDicomList::Bool)
#     #simpleITK object used to read from disk
#     reader = sitk.ImageSeriesReader()
#     if(isDicomList)# data in form of folder with dicom files
#         dicom_names = reader.GetGDCMSeriesFileNames(dirString)
#         reader.SetFileNames(dicom_names)
#         return reader.Execute()
#     elseif(isMHD) #mhd file
#         return sitk.ReadImage(dirString)
#     end
# end#getPixelDataAndSpacing

# """
# becouse Julia arrays is column wise contiguus in memory and open GL expects row wise we need to rotate and flip images
# pixels - 3 dimensional array of pixel data
# """
# function permuteAndReverse(pixels)
#     sizz=size(pixels)

# 1#
# 2# all upside down
# 3
# 21# on 3 dim correct
# 12
# 13
# 23
# 32# liver on the left
# 32
# 23
# 123
# 321# 3 upside down
# 312
# 213# 1,2 good ; 3 upside down
# 231

# 321 ; 2#3 upside down

# 321 ; 1#3 ok 1 upside down
# 2;12 # 3 ok 1 and 2 upside down

# 2; 13 # 3 upside down
# 2; 3# 3 ok 1 and 2 upside down
# 2;2# 1,2, ok 3 upside down
# 2;23# liver on the left evrywhere

# 32;1 # all upside down
# 32;21 # 3 left to right
# 32;213 # 3 left to right
# 32;13# 3 is ok 1 and 2 upside down
# 32;1 # 3 upside down

# 32; 3#- upside down
# 21 ; 1 #- upside down

# # for i in 1:sizz[3]
# #     pixels[:,:,i] =  reverse(pixels[:,:,i])
# #     end#

# #  for i in 1:sizz[2]
# #     pixels[:,i,:] =  reverse(pixels[:,i,:])
# #     end#

# #     for i in 1:sizz[1]
# #         pixels[i,:,:] =  reverse(pixels[i,:,:])
# #     end#

#     pixels=  permutedims(pixels, (3,2,1))
#     sizz=size(pixels)

#     for i in 1:sizz[1]
#         for j in 1:sizz[3]
#             pixels[i,:,j] =  reverse(pixels[i,:,j])
#         end#
#     end#

#     # for i in 1:sizz[1]
#     #     pixels[i,:,:] =  reverse(pixels[i,:,:])
#     # end#

#     # for i in 1:sizz[1]
#     #     pixels[i,:,:] =  reverse(pixels[i,:,:])
#     # end#

#     # for i in 1:sizz[2]
#     # pixels[:,i,:] =  reverse(pixels[:,i,:])
#     # end#

#     # for i in 1:sizz[1]
#     #     pixels[i,:,:] =  reverse(pixels[i,:,:])
#     # end#

#     return pixels
#   end#permuteAndReverse

# """
# given simple ITK image it reads associated pixel data - and transforms it by permuteAndReverse functions
# it will also return voxel spacing associated with the image
# """
# function getPixelsAndSpacing(image)
#     pixelsArr = np.array(sitk.GetArrayViewFromImage(image))# we need numpy in order for pycall to automatically change it into julia array
#     spacings = image.GetSpacing()
#     return ( permuteAndReverse(pixelsArr), spacings  )
# end#getPixelsAndSpacing

# # directories of PET/CT Data - from not published (yet) dataset - single example available from https://wwsi365-my.sharepoint.com/:f:/g/personal/s9956jm_ms_wwsi_edu_pl/Eq3cL7Md5bhPvnUlFLAMKZAB3nsbl6Q18fG96iVajvnNqA?e=bzX68X
# dirOfExample ="C:\\GitHub\\JuliaMedPipe\\data\\PETphd\\slicerExp\\all17\\bad17NL-bad17NL\\20150518-PET^1_PET_CT_WholeBody_140_70_Recon (Adult)\\4-CT AC WB  1.5  B30f"
# dirOfExamplePET ="C:\\GitHub\\JuliaMedPipe\\data\\PETphd\\slicerExp\\all17\\bad17NL-bad17NL\\20150518-PET^1_PET_CT_WholeBody_140_70_Recon (Adult)\\3-PET WB"

# # in most cases dimensions of PET and CT data arrays will be diffrent in order to make possible to display them we need to resample and make dimensions equal
# imagePET= getImageFromDirectory(dirOfExamplePET,false,true)
# ctImage= getImageFromDirectory(dirOfExample,false,true)
# pet_image_resampled = sitk.Resample(imagePET, ctImage)
# purePetPixels, PurePetSpacing = getPixelsAndSpacing(imagePET)
# ctPixels, ctSpacing = getPixelsAndSpacing(ctImage)
# # In my case PET data holds 64 bit floats what is unsupported by Opengl
# petPixels, petSpacing =getPixelsAndSpacing(pet_image_resampled)
# petPixels = Float32.(petPixels)

# # we need to pass some metadata about image array size and voxel dimensions to enable proper display
# datToScrollDimsB= MedEye3d.ForDisplayStructs.DataToScrollDims(imageSize=  size(ctPixels) ,voxelSize=PurePetSpacing, dimensionToScroll = 3 );
# # example of texture specification used - we need to describe all arrays we want to display
# textureSpecificationsPETCT = [
#   TextureSpec{Float32}(
#       name = "PET",
#       isNuclearMask=true,
#       # we point out that we will supply multiple colors
#       isContinuusMask=true,
#       #by the number 1 we will reference this data by for example making it visible or not
#       numb= Int32(1),
#       colorSet = [RGB(0.0,0.0,0.0),RGB(1.0,1.0,0.0),RGB(1.0,0.5,0.0),RGB(1.0,0.0,0.0)]
#       #display cutoff all values below 200 will be set 2000 and above 8000 to 8000 but only in display - source array will not be modified
#       ,minAndMaxValue= Float32.([200,8000])
#      ),
#   TextureSpec{UInt8}(
#       name = "manualModif",
#       numb= Int32(2),
#       color = RGB(0.0,1.0,0.0)
#       ,minAndMaxValue= UInt8.([0,1])
#       ,isEditable = true
#      ),

#      TextureSpec{Int16}(
#       name= "CTIm",
#       numb= Int32(3),
#       isMainImage = true,
#       minAndMaxValue= Int16.([0,100]))
# ];
# # We need also to specify how big part of the screen should be occupied by the main image and how much by text fractionOfMainIm= Float32(0.8);
# fractionOfMainIm= Float32(0.8);
# """
# If we want to display some text we need to pass it as a vector of SimpleLineTextStructs - utility function to achieve this is
# textLinesFromStrings() where we pass resies of strings, if we want we can also edit those structs to controll size of text and space to next line look into SimpleLineTextStruct doc
# mainLines - will be displayed over all slices
# supplLines - will be displayed over each slice where is defined - below just dummy data
# """
# import MedEye3d.DisplayWords.textLinesFromStrings

# mainLines= textLinesFromStrings(["main Line1", "main Line 2"]);
# supplLines=map(x->  textLinesFromStrings(["sub  Line 1 in $(x)", "sub  Line 2 in $(x)"]), 1:size(ctPixels)[3] );

# """
# If we want to pass 3 dimensional array of scrollable data we need to supply it via vector ThreeDimRawDat's struct
# utility function to make creation of those easier is getThreeDims which creates series of ThreeDimRawDat from list of tuples where
#     first entry is String and second entry is 3 dimensional array with data
#     strings needs to be the same as we  defined in texture specifications at the bagining
#     data arrays needs to be o the same size and be of the same type we specified in texture specification
# """
# import MedEye3d.StructsManag
# tupleVect = [("PET",petPixels) ,("CTIm",ctPixels),("manualModif",zeros(UInt8,size(petPixels)) ) ]
# slicesDat= MedEye3d.StructsManag.getThreeDims(tupleVect )
# """
# Holds data necessary to display scrollable data
# """
# mainScrollDat = FullScrollableDat(dataToScrollDims =datToScrollDimsB
#                                  ,dimensionToScroll=1 # what is the dimension of plane we will look into at the beginning for example transverse, coronal ...
#                                  ,dataToScroll= slicesDat
#                                  ,mainTextToDisp= mainLines
#                                  ,sliceTextToDisp=supplLines );

# tupleVect = [("PET",petPixels) ,("CTIm",ctPixels),("manualModif",zeros(UInt8,size(petPixels)) ) ]
# slicesDat=  MedEye3d.StructsManag.getThreeDims(tupleVect )
# """
# Holds data necessary to display scrollable data
# """
# mainScrollDat = FullScrollableDat(dataToScrollDims =datToScrollDimsB
#                                  ,dimensionToScroll=1 # what is the dimension of plane we will look into at the beginning for example transverse, coronal ...
#                                  ,dataToScroll= slicesDat
#                                  ,mainTextToDisp= mainLines
#                                  ,sliceTextToDisp=supplLines );

# """
# This function prepares all for display; 1000 in the end is responsible for setting window width for more look into SegmentationDisplay.coordinateDisplay
# """
# SegmentationDisplay.coordinateDisplay(textureSpecificationsPETCT ,fractionOfMainIm ,datToScrollDimsB ,1200);
# """
# as all is ready we can finally display image
# """
# Main.SegmentationDisplay.passDataForScrolling(mainScrollDat);
# Main.SegmentationDisplay.mainActor.actor.currentDisplayedSlice = 800

# #########Interactions

# """
# Next all Interactions are done either by mouse or by keyboard shortcuts
# left click and drag - will mark active texture (look below - set with alt ...)
#     if it is set to be modifiable in the texture specifications, to the set value and size (by tab...)
# right click and drag - sets remembered position - when we will change plane of crossection
#      for example from tranverse to coonal this point will be also visible on new plane
# all keyboard shortcuts will be activated on RELEASE of keys or by pressing enter while still pressing
# shift + number - make mask associated with given number visible
# ctrl + number -  make mask associated with given number invisible
# alt + number -  make mask associated with given number active for mouse interaction
# tab + number - sets the number that will be  used as an input to masks modified by mouse
#     when tab plus (and then no number) will be pressed it will increase stroke width
#     when tab minus (and then no number) will be pressed it will increase stroke width
# shift + numberA + "-"(minus sign) +numberB  - display diffrence between masks associated with numberA and numberB - also it makes automaticall mask A and B invisible
# ctrl + numberA + "-"(minus sign) +numberB  - stops displaying diffrence between masks associated with numberA and numberB - also it makes automaticall mask A and B visible
# space + 1 or 2 or 3 - change the plane of view (transverse, coronal, sagittal)
# ctrl + z - undo last action
# tab +/- increase or decrease stroke width
# F1 - will display wide window for bone Int32(1000),Int32(-1000)
# F2 - will display window for soft tissues Int32(400),Int32(-200)
# F3 - will display wide window for lung viewing  Int32(0),Int32(-1000)
# KEY_F4,  KEY_F5 -
#     sets minimum (F4) and maximum (KEY_F5) value for display (with combination of + and minus signs - to increase or decrease given treshold) -
#         in case of continuus colors it will clamp values - so all above max will be equaled to max ; and min if smallert than min
#         in case of main CT mask - it will controll min shown white and max shown black
#         in case of maks with single color associated we will step data so if data is outside the rande it will return 0 - so will not affect display
# """

# # ######### Benchmark PET/CT
# # #For transparency I include Below code used to benchark  PET/CT data

# # window = Main.SegmentationDisplay.mainActor.actor.mainForDisplayObjects.window
# # syncActor = Main.SegmentationDisplay.mainActor

# # using GLFW,DataTypesBasic, ModernGL,Setfield
# # using BenchmarkTools

# # BenchmarkTools.DEFAULT_PARAMETERS.samples = 100
# # BenchmarkTools.DEFAULT_PARAMETERS.seconds =5000
# # BenchmarkTools.DEFAULT_PARAMETERS.gcsample = true

# # function toBenchmarkScroll(toSc)
# #     MedEye3d.ReactToScroll.reactToScroll(toSc ,syncActor, false)
# # end

# # function toBenchmarkPaint(carts)
# #     MedEye3d.ReactOnMouseClickAndDrag.reactToMouseDrag(MouseStruct(true,false, carts),syncActor )
# # end

# # function toBenchmarkPlaneTranslation(toScroll)
# #     MedEye3d.ReactOnKeyboard.processKeysInfo(Option(toScroll),syncActor,KeyboardStruct(),false    )
# #     OpenGLDisplayUtils.basicRender(syncActor.actor.mainForDisplayObjects.window)
# #     glFinish()
# # end

# # function prepareRAndomCart(randInt)
# #     return [CartesianIndex(12+randInt,13+randInt),CartesianIndex(12+randInt,15+randInt),CartesianIndex(12+randInt,18+randInt),CartesianIndex(2+randInt,10+randInt),CartesianIndex(2+randInt,14+randInt)]
# # end

# # #we want some integers but not 0
# # sc = @benchmarkable toBenchmarkScroll(y) setup=(y = filter(it->it!=0, rand(-5:5,20))[1]  )

# # paint =  @benchmarkable toBenchmarkPaint(y) setup=(y =  prepareRAndomCart(rand(1:40,1)[1]  ))

# # translations =  @benchmarkable toBenchmarkPlaneTranslation(y) setup=(y = setproperties(syncActor.actor.onScrollData.dataToScrollDims,  (dimensionToScroll=rand(1:3,2)[1])) )

# # using BenchmarkPlots, StatsPlots
# # # Define a parent BenchmarkGroup to contain our suite

# # scrollingPETCT = run(sc)
# # mouseInteractionPETCT = run(paint)
# # translationsPETCT = run(translations)

# # plot(scrollingPETCT)

# # #### PURE CT image exaple , MHD files  taken from https://sliver07.grand-challenge.org/

# # exampleLabel = "/home/jakub/JuliaProjects/Probabilistic-medical-segmentation/data/liverPrimData/training-labels/label/liver-seg002.mhd"
# # exampleCTscann = "/home/jakub/JuliaProjects/Probabilistic-medical-segmentation/data/liverPrimData/training-scans/scan/liver-orig002.mhd"

# # # loading data
# # imagePureCT= getImageFromDirectory(exampleCTscann,true,false)
# # imageMask= getImageFromDirectory(exampleLabel,true,false)

# # ctPixelsPure, ctSpacingPure = getPixelsAndSpacing(imagePureCT)
# # maskPixels, maskSpacing =getPixelsAndSpacing(imageMask)

# # # we need to pass some metadata about image array size and voxel dimensions to enable proper display
# # datToScrollDimsB= MedEye3d.ForDisplayStructs.DataToScrollDims(imageSize=  size(ctPixelsPure) ,voxelSize=ctSpacingPure, dimensionToScroll = 3 );
# # # example of texture specification used - we need to describe all arrays we want to display
# # listOfTexturesSpec = [
# #     TextureSpec{UInt8}(
# #         name = "goldStandardLiver",
# #         numb= Int32(1),
# #         color = RGB(1.0,0.0,0.0)
# #         ,minAndMaxValue= Int8.([0,1])
# #        ),
# #     TextureSpec{UInt8}(
# #         name = "manualModif",
# #         numb= Int32(2),
# #         color = RGB(0.0,1.0,0.0)
# #         ,minAndMaxValue= UInt8.([0,1])
# #         ,isEditable = true
# #        ),

# #        TextureSpec{Int16}(
# #         name= "CTIm",
# #         numb= Int32(3),
# #         isMainImage = true,
# #         minAndMaxValue= Int16.([0,100]))
# #  ];
# # # We need also to specify how big part of the screen should be occupied by the main image and how much by text fractionOfMainIm= Float32(0.8);
# # fractionOfMainIm= Float32(0.8);
# # """
# # If we want to display some text we need to pass it as a vector of SimpleLineTextStructs
# # """
# # import MedEye3d.DisplayWords.textLinesFromStrings

# # mainLines= textLinesFromStrings(["main Line1", "main Line 2"]);
# # supplLines=map(x->  textLinesFromStrings(["sub  Line 1 in $(x)", "sub  Line 2 in $(x)"]), 1:size(ctPixelsPure)[3] );

# # """
# # If we want to pass 3 dimensional array of scrollable data"""
# # import MedEye3d.StructsManag.getThreeDims

# # tupleVect = [("goldStandardLiver",maskPixels) ,("CTIm",ctPixelsPure),("manualModif",zeros(UInt8,size(ctPixelsPure)) ) ]
# # slicesDat= MedEye3d.StructsManag.getThreeDims(tupleVect )

# # """
# # Holds data necessary to display scrollable data
# # """
# # mainScrollDat = FullScrollableDat(dataToScrollDims =datToScrollDimsB
# #                                  ,dimensionToScroll=1 # what is the dimension of plane we will look into at the beginning for example transverse, coronal ...
# #                                  ,dataToScroll= slicesDat
# #                                  ,mainTextToDisp= mainLines
# #                                  ,sliceTextToDisp=supplLines );

# # """
# # This function prepares all for display; 1000 in the end is responsible for setting window width for more look into SegmentationDisplay.coordinateDisplay
# # """
# # SegmentationDisplay.coordinateDisplay(listOfTexturesSpec ,fractionOfMainIm ,datToScrollDimsB ,1000);
# # """
# # as all is ready we can finally display image
# # """
# # Main.SegmentationDisplay.passDataForScrolling(mainScrollDat);

# # ###########next part of benchmark

# # #we want some integers but not 0
# # scPureCt = @benchmarkable toBenchmarkScroll(y) setup=(y = filter(it->it!=0, rand(-5:5,20))[1]  )

# # paintPureCt =  @benchmarkable toBenchmarkPaint(y) setup=(y =  prepareRAndomCart(rand(1:40,1)[1]  ))

# # translationsPureCt =  @benchmarkable toBenchmarkPlaneTranslation(y) setup=(y = setproperties(syncActor.actor.onScrollData.dataToScrollDims,  (dimensionToScroll=rand(1:3,2)[1])) )

# # using BenchmarkPlots, StatsPlots

# # scB = @benchmarkable toBenchmarkScroll(y) setup=(y = filter(it->it!=0, rand(-5:5,20))[1]  )

# # paintB =  @benchmarkable toBenchmarkPaint(y) setup=(y =  prepareRAndomCart(rand(1:40,1)[1]  ))

# # translationsB =  @benchmarkable toBenchmarkPlaneTranslation(y) setup=(y = setproperties(syncActor.actor.onScrollData.dataToScrollDims,  (dimensionToScroll=rand(1:3,2)[1])) )

# # using BenchmarkPlots, StatsPlots
# # # Define a parent BenchmarkGroup to contain our suite

# # scrollingPureCT = run(scB)
# # mouseInteractionPureCT = run(paintB)
# # translationsPureCT = run(translationsB)

# #  using GLFW
# #  GLFW.PollEvents()


text:
# Dantzig-Wolfe Reformulation and Column Generation
# Functions for Master Problem
# Edward J. Xu
# 2019.5.25
########################################################################################################################

function setGAPModel(w, p, c)
    (numSub, numXPerSub) = size(w)
    # mod = Model(solver = CplexSolver(CPX_PARAM_TILIM = 1))
    modAll = Model(solver = GurobiSolver())
    @variable(modAll, x[1:numSub, 1:numXPerSub], Bin)
    @objective(modAll, Min, sum(p[i,j] * x[i,j] for i = 1:numSub for j = 1:numXPerSub))
    # each job must be served
    @constraint(modAll, [j = 1:numXPerSub], sum(x[i,j] for i = 1:numSub) == 1)
    @constraint(modAll, [i = 1:numSub], sum(w[i,j] * x[i,j] for j = 1:numXPerSub) <= c[i])
    ##
    # define blocks (Each block becomes a sub-problem)
    blocks = []
    for i = 1:numSub
        # constraint (numXPerSub + i) goes into block i
        push!(blocks, [numXPerSub + i])
    end
    # When using cplex we need to start solving the model before we can
    # export the constraint matrix and so on.
    # solve(mod)
    JuMP.build(modAll)
    return (modAll, numXPerSub, blocks)
end

function readData(strFileName::String, series::Int64)
    # function for reading GAP files from the OR-library
    open(strFileName) do f
        allFile = read(f, String)
        allNumbers = map(x -> parse(Int64,x), split(allFile))
        # return allNumbers
        numSeries = allNumbers[1]
        println("File contains $numSeries instances.")
        if series > numSeries
            throw("readData(...): specified file does not contain $(series) example")
        end
        index = 2
        for inst = 1:numSeries
            numSub = allNumbers[index]
            numXPerSub = allNumbers[index + 1]
            println("numSub = $(numSub), numXPerSub = $(numXPerSub)")
            w = zeros(numSub, numXPerSub)
            p = zeros(numSub, numXPerSub)
            c = zeros(numSub)
            index += 2
            println("start index when reading p: $index")
            for i = 1:numSub
                for j = 1:numXPerSub
                    p[i,j] = allNumbers[index]
                    index += 1
                end
            end
            println("start index when reading w: $index")
            for i = 1:numSub
                for j = 1:numXPerSub
                    w[i,j] = allNumbers[index]
                    index += 1
                end
            end
            println("start index when reading c: $index")
            for i=1:numSub
                c[i] = allNumbers[index]
                index += 1
            end
            if inst == series
                return (w, p, c)
            end
        end
    end
end

function getMat(mod)
    JuMP.build(mod)
    mpb = MathProgBase
    modInternal = internalmodel(mod)
    vecBoundLowerCons = mpb.getconstrLB(modInternal)
    vecBoundUpperCons = mpb.getconstrUB(modInternal)
    vecSense = Vector{Sense}(undef, length(vecBoundLowerCons))
    vec_b = Vector{Float64}(undef, length(vecBoundLowerCons))
    for i = 1:length(vecBoundLowerCons)
        if vecBoundLowerCons[i] == vecBoundUpperCons[i]
            vecSense[i] = eq
            vec_b[i] = vecBoundLowerCons[i]
        elseif vecBoundLowerCons[i] == -Inf
            vecSense[i] = leq
            vec_b[i] = vecBoundUpperCons[i]
        elseif vecBoundUpperCons[i] == Inf
            vecSense[i] = geq
            vec_b[i] = vecBoundLowerCons[i]
        else
            throw("Error: One of the constraints in the input model is a \"range\" constraint. This is not supported.")
        end
    end
    mat_a = mpb.getconstrmatrix(modInternal)
    vec_c = mpb.getobj(modInternal)
    vecBoundLowerVar = mpb.getvarLB(modInternal)
    vecBoundUpperVar = mpb.getvarUB(modInternal)
    vecTypeVar = mpb.getvartype(modInternal)
    return (mat_a, vec_b, vec_c, vecSense)
    # return (mat_a, vec_b, vec_c, vecBoundLowerVar, vecBoundUpperVar, vecTypeVar, vecSense)
end

function addVarsFromRow(Arow, varBitSet)
    for idx=1:length(Arow.nzval)
        if Arow.nzval[idx] != 0
            push!(varBitSet, Arow.nzind[idx])
        end
    end
end

function getIndex(mat_a, blocks)
    # find out which variables that belong to which block:
    numSub = length(blocks)
    if numSub < 1
        throw("error. We expect at least one block")
    end
    (numCons, numVar) = size(mat_a)
    varsInSub = Vector{BitSet}(undef, numSub)
    for k = 1:numSub
        varsInSub[k] = BitSet()
        for row in blocks[k]
            addVarsFromRow(mat_a[row,:], varsInSub[k])
        end
    end
    # check that the blocks are not sharing any variables
    for k1 = 1:numSub
        for k2 = k1 + 1:numSub
            if !isempty(intersect(varsInSub[k1],varsInSub[k2]))
                throw("Blocks $k1 and $k2 share variables. This is not supported.")
            end
        end
    end
    # Check that all variables belong to a block
    setUnion = deepcopy(varsInSub[1])
    for k = 2:numSub
        union!(setUnion, varsInSub[k])
    end
    for j = 1:numVar
        if !(j in setUnion)
            throw("variable $j does not belong to any block. This is not supported")
        end
    end
    # Find the rows that are left in master problem:
    indexMas = BitSet(1:numCons)
    for block in blocks
        setdiff!(indexMas, block)
    end
    ## IndexSub
    indexSub = Vector{Vector{Int64}}(undef, numSub)
    for k = 1:numSub
        indexSub[k] = collect(varsInSub[k])
    end
    return (indexMas, indexSub)
end

function setGeneralAssignProbDate(strFileName, series)
    println("#### 1.1/5,  Read Data #########################################################") ########################
    (w, p, c) = readData(strFileName, series)
    (modAll, numXPerSub, blocks) = setGAPModel(w, p, c)
    println("#### 1.2/5,  Get Index #########################################################") ########################
    (mat_a, vec_b, vec_c, vecSenseAll) = getMat(modAll)
    (indexMas, indexSub) = getIndex(mat_a, blocks)
    return (mat_a, vec_b, vec_c, vecSenseAll, indexMas, blocks, indexSub, numXPerSub, numXPerSub)
end


text:
export create,
    retrieve,
    update,
    delete

"""
    create(c::Client, ::Type{T}, args...; kwargs...)

Create, add, send, etc.

# Examples
```julia-repl
# Send a message.
julia> create(c, Message, channel; content="foo")

# Create a new channel.
julia> create(c, DiscordChannel, guild; name="bar")

# Ban a member.
julia> create(c, Ban, guild, member; reason="baz")
```
"""
function create end

"""
    retrieve(c::Client, ::Type{T}, args...; kwargs...)

Retrieve, get, list, etc.

# Examples
```julia-repl
# Get the client user.
julia> retrieve(c, User)

# Get a guild's channels.
julia> retrieve(c, DiscordChannel, guild)

# Get an Invite to a guild  by code.
julia> retrieve(c, Invite, "abcdef")
```
"""
function retrieve end

"""
    update(c::Client, x::T, args...; kwargs...)

Update, edit, modify, etc.

# Examples
```julia-repl
# Edit a message.
julia> update(c, message; content="foo2")

# Modify a webhook.
julia> update(c, webhook; name="bar2")

# Update a role.
julia> update(c, role, guild; permissions=8)
```
"""
function update end

"""
    delete(c::Client, x::T, args...)

Delete, remove, discard, etc.

# Examples
```julia-repl
# Kick a member.
julia> delete(c, member)

# Unban a member.
julia> delete(c, ban, guild)

# Delete all reactions from a message.
# This is the only update/delete method which takes a type parameter.
delete(c, Reaction, message)
```
"""
function delete end

include(joinpath("crud", "audit_log.jl"))
include(joinpath("crud", "ban.jl"))
include(joinpath("crud", "channel.jl"))
include(joinpath("crud", "emoji.jl"))
include(joinpath("crud", "guild_embed.jl"))
include(joinpath("crud", "guild.jl"))
include(joinpath("crud", "integration.jl"))
include(joinpath("crud", "invite.jl"))
include(joinpath("crud", "member.jl"))
include(joinpath("crud", "message.jl"))
include(joinpath("crud", "overwrite.jl"))
include(joinpath("crud", "reaction.jl"))
include(joinpath("crud", "role.jl"))
include(joinpath("crud", "user.jl"))
include(joinpath("crud", "voice_region.jl"))
include(joinpath("crud", "webhook.jl"))

#=
Functions not covered here:
- trigger_typing_indicator
- group_dm_add_recipient
- group_dm_remove_recipient
- modify_guild_channel_positions
- modify_current_user_nick
- modify_guild_role_positions
- get_guild_prune_count
- begin_guild_prune
- sync_guild_integration (combine with update?)
- get_guild_vanity_url
- get_guild_widget_image
- leave_guild (combine with delete?)
- execute_webhook
- execute_slack_compatible_webhook
- execute_github_compatible_webhook
=#


text:
export Uninformative

struct Uninformative end

@node Uninformative Deterministic [ out ]


text:
StructTypes.subtypes(::Type{AbstractBenchmark}) = (benchmark_group=BenchmarkGroup,)
const BENCHMARKEXT_VERSION = v"0.1.0"

function snake_case(camelsymbol::Symbol)
    camelstring = String(camelsymbol)

    wordpat = r"
    ^[a-z]+ |                  #match initial lower case part
    [A-Z][a-z]+ |              #match Words Like This
    \d*([A-Z](?=[A-Z]|$))+ |   #match ABBREV 30MW
    \d+                        #match 1234 (numbers without units)
    "x

    smartlower(word) = any(islowercase, word) ? lowercase(word) : word
    words = [smartlower(m.match) for m in eachmatch(wordpat, camelstring)]

    Symbol(join(words, "_"))
end

# TODO: Add any new types as they're added
const SUPPORTED_TYPES = Dict{Symbol, Symbol}(Base.typename(x).name => snake_case(Base.typename(x).name) for x in [Parameters, BenchmarkGroup, TagFilter, Trial,
    TrialEstimate, TrialJudgement, TrialRatio])
# n.b. Benchmark type not included here, since it is gensym'd

# We can make a macro, which generate this structure when new types appear.
# Or we can just add new fields manually
mutable struct Serializator
    julia::String
    benchmark_version::String
    parameters::Dict{String, Parameters}
    benchmark_group::Dict{String, BenchmarkGroup}
    tag_filter::Dict{String, TagFilter}
    trial::Dict{String, Trial}
    trial_estimate::Dict{String, TrialEstimate}
    trial_judgement::Dict{String, TrialJudgement}
    trial_ratio::Dict{String, TrialRatio}
end
StructTypes.StructType(::Type{Serializator}) = StructTypes.Mutable()
StructTypes.omitempties(::Type{Serializator}) = true
Serializator(julia, version) = Serializator(julia, version, Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), Dict())
Serializator() = Serializator("", "")

function badext(filename)
    noext, ext = splitext(filename)
    msg = "Only JSON serialization is supported, please provide file name in the form \"$noext.json\"."
    throw(ArgumentError(msg))
end

"""
    save(file::Union{String, IO}, args...)

Save `args` to a json file. Arguments `args` can be either `BenchmarkExt` objects or pairs of the form `name => arg`.
"""
function save(filename::AbstractString, args...)
    endswith(filename, ".json") || badext(filename)
    open(filename, "w") do io
        save(io, args...)
    end
end

function save(io::IO, args...)
    isempty(args) && throw(ArgumentError("Nothing to save"))
    ser = Serializator(string(VERSION), string(BENCHMARKEXT_VERSION))
    store = false
    for arg in args
        name = typeof(arg).name.name
        if !haskey(SUPPORTED_TYPES, name)
            arg isa Pair || throw(ArgumentError("Only BenchmarkExt types can be serialized, encountered \"$name\"."))
            name2 = typeof(arg[2]).name.name
            haskey(SUPPORTED_TYPES, name2) || throw(ArgumentError("Only BenchmarkExt types can be serialized, encountered \"$name2\"."))
            getfield(ser, SUPPORTED_TYPES[name2])[string(arg[1])] = arg[2]
        else
            field = getfield(ser, SUPPORTED_TYPES[name])
            field[string(SUPPORTED_TYPES[name])*"_"*string(length(field) + 1)] = arg
        end
        store = true
    end
    isempty(store) && error("Nothing to save")
    JSON3.write(io, ser)

    return nothing
end

"""
    load(file::Union{String, IO}, pattern = "", group = true)

Load benchmark values from json `file`. If `pattern` is not `nothing` then results will be filtered according to `pattern`. If `group` is `true`, then result is given as `Vector{Any}` where all values from all types combined together, otherwise it is returned as `BenchmarkExt.Serializator`.
"""
function load(filename::AbstractString, pattern = nothing, group = true)
    endswith(filename, ".json") || badext(filename)
    open(filename, "r") do f
        load(f, pattern, group)
    end
end

function load(io::IO, pattern = nothing, group = true)
    parsed = JSON3.read(io, Serializator)
    if isempty(parsed.julia) || isempty(parsed.benchmark_version)
        error("Unexpected JSON format")
    end
    if pattern !== nothing
        for name in values(SUPPORTED_TYPES)
            field = getfield(parsed, name)
            filter!(x -> occursin(pattern, x[1]), field)
        end
    end

    group || return parsed

    res = []
    for name in fieldnames(Serializator)
        name in values(SUPPORTED_TYPES) || continue
        append!(res, values(getfield(parsed, name)))
    end

    return res
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

insts = OrderedDict()
insts["minimal"] = [
    ((2, false),),
    ((2, true),),
    ((4, false),),
    ((4, true),),
    ]
insts["fast"] = [
    ((50, false),),
    ((50, true),),
    ((1000, false),),
    ((250, true),),
    ]
insts["various"] = [
    ((1000, false),),
    ((500, true),),
    ((4000, false),),
    ((1000, true),),
    ((8000, false),),
    ((2000, true),),
    ((10000, false),),
    ]
return (DiscreteMaxLikelihood, insts)