text:
function filter_line_list(df::DataFrame, inst::IT ; λmin::Real = default_λmin, λmax::Real = default_λmax ) where { IT<:HARPS.AnyHARPS }
   df |> @filter(λmin <= _.lambda <= λmax) |>
    #    @filter( _.lambda < 6000.0 ) |>                       # Avoid tellurics at redder wavelengths
    #    @filter( _.lambda >6157 || _.lambda < 6155  ) |>   # Avoid "line" w/ large variability
    DataFrame
end

"""  read_telluric_ranges( filename )
Return DataFrame (keys lambda_lo and lambda_hi) with wavelength ranges to be avoided as containing tellurics
based on provided CSV file.  Assumes path is included in filename.
"""
function read_telluric_ranges(fn::String)
    @assert occursin(r"\.csv$",fn)
    #println("Trying to read from >",fn,"<.")
    @assert isfile(fn)
    df = CSV.read(fn, DataFrame)
    @assert size(df,1) >= 1
    @assert hasproperty(df,:lambda_lo)
    @assert hasproperty(df,:lambda_hi)
    return df
end

function make_λ_list_for_bad_columns(line_list::DataFrame, harps_data::DT )  where {
            T1<:Real, A1<:AbstractArray{T1}, T2<:Real, A2<:AbstractArray{T2}, T3<:Real, A3<:AbstractArray{T3},
             IT<:HARPS.AnyHARPS, ST<:Spectra2DBasic{T1,T2,T3,A1,A2,A3,IT}, DT<:AbstractArray{ST,1} }
    order_list = 1:size(first(harps_data).λ,2)
    df_bad_col_λs = DataFrame(:order=>Int[], :lambda_lo=>Float64[], :lambda_hi=>Float64[])
    for order in order_list
        λ_extrema = NaNMath.extrema(first(harps_data).λ[:,order])
        if isnan(first(λ_extrema)) ||  isnan(last(λ_extrema))   continue   end
        for bcr in bad_col_ranges(HARPS2D(), order)
            pixlo = first(bcr)
            pixhi = pixlo + 1
            if pixhi > size(first(harps_data).λ,1)
                pixhi -= 1
                pixlo -= 1
            end
            doppler_factors = map(obsid-> haskey(harps_data[obsid].metadata,:doppler_factor) ? harps_data[obsid].metadata[:doppler_factor] : 1 , 1:length(harps_data))
            Δλ_pixel = (first(harps_data).λ[pixhi,order] - first(harps_data).λ[pixlo,order]) * doppler_factors[1]
            (λ_lo, λ_hi) = mapreduce(obsid->extrema(harps_data[obsid].λ[bcr,order]) .* doppler_factors[obsid],
                               (a,b) -> (min(a[1],b[1]), max(a[2],b[2])),   1:length(harps_data) )
            λ_lo -= Δλ_pixel/2
            λ_hi += Δλ_pixel/2
            push!(df_bad_col_λs, Dict(:order=>order, :lambda_lo=>λ_lo, :lambda_hi=>λ_hi) )
        end
    end
    return df_bad_col_λs
end

#=
#old neid functions

function make_good_orders_pixels_df(neid_data::DT ; orders::A4 = orders_to_use_default(first(neid_data).inst) ) where {
                T1<:Real, A1<:AbstractArray{T1}, T2<:Real, A2<:AbstractArray{T2}, T3<:Real, A3<:AbstractArray{T3},
                 IT<:NEID.AnyNEID, ST<:Spectra2DBasic{T1,T2,T3,A1,A2,A3,IT}, DT<:AbstractArray{ST,1},
                 T4 <: Integer, A4<:AbstractArray{T4} }
    @warn "Deprecated use make_good_orders_pixels_df(inst) instead."
    df = DataFrame(order=Int[],pixels=UnitRange[])
    for ord in orders
        pixel_array = EchelleInstruments.calc_complement_index_ranges(get_pixel_range(NEID2D(),ord),EchelleInstruments.NEID.bad_col_ranges(NEID2D(),ord))
        order_array = repeat([ord],length(pixel_array))
        df_tmp = DataFrame(:order=>order_array,:pixels=>pixel_array)
        append!(df,df_tmp)
    end
    return df
end

function make_good_orders_pixels_df(inst::IT ; orders::A4 = orders_to_use_default(inst) ) where {
                 IT<:NEID.AnyNEID, T4 <: Integer, A4<:AbstractArray{T4} }
    df = DataFrame(order=Int[],pixels=UnitRange[])
    for ord in orders
        pixel_array = EchelleInstruments.calc_complement_index_ranges(get_pixel_range(NEID2D(),ord),EchelleInstruments.NEID.bad_col_ranges(NEID2D(),ord))
        order_array = repeat([ord],length(pixel_array))
        df_tmp = DataFrame(:order=>order_array,:pixels=>pixel_array)
        append!(df,df_tmp)
    end
    return df
end
=#

#=
function make_clean_line_list_from_bad_columns(line_list::DataFrame, neid_data::DT ;
                Δv_to_avoid_tellurics::Real = 0.0, v_center_to_avoid_tellurics::Real = 0.0
                )  where { T1<:Real, A1<:AbstractArray{T1}, T2<:Real, A2<:AbstractArray{T2}, T3<:Real, A3<:AbstractArray{T3}, IT<:NEID.AnyNEID, ST<:Spectra2DBasic{T1,T2,T3,A1,A2,A3,IT}, DT<:AbstractArray{ST,1} }

                if size(df_bad_col_λs,1) < 1 return line_list end
                line_list_to_search_for_tellurics = copy(line_list)
                line_list_to_search_for_tellurics.lambda_lo = line_list_to_search_for_tellurics.lambda./calc_doppler_factor(Δv_to_avoid_tellurics).*calc_doppler_factor(v_center_to_avoid_tellurics)
                line_list_to_search_for_tellurics.lambda_hi = line_list_to_search_for_tellurics.lambda.*calc_doppler_factor(Δv_to_avoid_tellurics).*calc_doppler_factor(v_center_to_avoid_tellurics)

                line_list_no_tellurics_df = line_list_to_search_for_tellurics |> # @filter(λmin <= _.lambda <= λmax) |>
                   @filter( !RvSpectMLBase.is_in_wavelength_range_list(_.lambda, order=_.order, list=df_bad_col_λs ) ) |>
                   DataFrame
end
=#

#old neid functions

"""
   `make_clean_line_list_from_tellurics(line_list, harps_data; Δv_to_avoid_tellurics )`
Returns a new line list that excludes lines with telluric contamination.
Inputs:
- `line_list`:  Dataframe containing field lambda
# - `harps_data`:  Array of spectra
- `Δv_to_avoid_tellurics`:  in m/s
- `v_center_to_avoid_tellurics` : in m/s
- `tellurics_filename`:  filename with wavelength ranges affected by tellurics
Outputs:
- `line_list_without_tellurics`:   DataFrame with fields: `lambda`, `weight`, `lambda_lo`, and `lambda_hi`.
Warning: Currently, reads in a list of wavelength ranges with tellurics and this file was created based on tellurics metadata in EXPRES data for HR 101501.
"""
function make_clean_line_list_from_tellurics(line_list::DataFrame, harps_data::DT ;
    Δv_to_avoid_tellurics::Real = 0.0, v_center_to_avoid_tellurics::Real = 0.0, tellurics_filename::String = joinpath(pkgdir(EchelleInstruments),"data/neid/telluric_ranges.csv")
               ) where { T1<:Real, A1<:AbstractArray{T1}, T2<:Real, A2<:AbstractArray{T2}, T3<:Real, A3<:AbstractArray{T3}, IT<:HARPS.AnyHARPS, ST<:Spectra2DBasic{T1,T2,T3,A1,A2,A3,IT}, DT<:AbstractArray{ST,1} }
   @assert 0.5*RvSpectMLBase.max_bc_earth_rotation <= Δv_to_avoid_tellurics <= 4*RvSpectMLBase.max_bc
   @error "This function does not yet work for HARPS data."
   if hasproperty(line_list,:lambda_lo) &&  hasproperty(line_list,:lambda_hi)
      line_list_to_search_for_tellurics = copy(line_list)
   else
      line_list_to_search_for_tellurics = RvSpectMLBase.add_line_boundaries_to_line_list(line_list, Δv_to_avoid_tellurics=Δv_to_avoid_tellurics, v_center_to_avoid_tellurics=v_center_to_avoid_tellurics )
   end
   #chunk_list_timeseries = RvSpectMLBase.make_chunk_list_timeseries_around_lines(harps_data,line_list_to_search_for_tellurics)

   if !hasproperty(line_list, :order)
       order_info = get_order_info(harps_data) # , orders_to_use=orders_to_use)
       #println("# order_info contains: ", names(order_info))
       #println("# line_list_df contains: ", names(line_list_df))
       line_list_to_search_for_tellurics = assign_lines_to_orders(line_list_to_search_for_tellurics, order_info, Δv_to_avoid_tellurics=Δv_to_avoid_tellurics, v_center=v_center_to_avoid_tellurics )
   end

   telluric_ranges = read_telluric_ranges(tellurics_filename)   # already includes range due to BCs
   # TODO: Check sign convention for v_center_to_avoid_tellurics
   # telluric_ranges.lambda_lo .*= calc_doppler_factor(v_center_to_avoid_tellurics) # / calc_doppler_factor(Δv_to_avoid_tellurics)
   # telluric_ranges.lambda_hi .*= calc_doppler_factor(v_center_to_avoid_tellurics) # * calc_doppler_factor(Δv_to_avoid_tellurics)

   #println("# make_clean_line_list_from_tellurics: λ_list_for_bad_cols_df /= calc_doppler_factor ")
   λ_list_for_bad_cols_df = make_λ_list_for_bad_columns(line_list_to_search_for_tellurics, harps_data)
   λ_list_for_bad_cols_df.lambda_lo .*= calc_doppler_factor(v_center_to_avoid_tellurics) / calc_doppler_factor(Δv_to_avoid_tellurics)
   λ_list_for_bad_cols_df.lambda_hi .*= calc_doppler_factor(v_center_to_avoid_tellurics) * calc_doppler_factor(Δv_to_avoid_tellurics)

   line_list_no_tellurics_df = line_list_to_search_for_tellurics |> # @filter(λmin <= _.lambda <= λmax) |>
      @filter( !RvSpectMLBase.is_in_wavelength_range_list(_.lambda, order=_.order, list=λ_list_for_bad_cols_df ) ) |>
      @filter( !RvSpectMLBase.is_in_wavelength_range_list(_.lambda, list=telluric_ranges ) ) |>
      DataFrame

end

#=
function add_line_boundaries_to_line_list(line_list::DataFrame; Δv_to_avoid_tellurics::Real = 0.0, v_center_to_avoid_tellurics::Real = 0.0 )
    @assert 0.5*RvSpectMLBase.max_bc_earth_rotation <= Δv_to_avoid_tellurics <= 4*RvSpectMLBase.max_bc
    line_list_to_search_for_tellurics = copy(line_list)
    line_list_to_search_for_tellurics.lambda_lo = line_list_to_search_for_tellurics.lambda./calc_doppler_factor(Δv_to_avoid_tellurics).*calc_doppler_factor(v_center_to_avoid_tellurics)
    line_list_to_search_for_tellurics.lambda_hi = line_list_to_search_for_tellurics.lambda.*calc_doppler_factor(Δv_to_avoid_tellurics).*calc_doppler_factor(v_center_to_avoid_tellurics)
    return line_list_to_search_for_tellurics
end
=#

function choose_obs_idx_for_init_guess(df::DataFrame, inst::IT )  where { IT<:HARPS.AnyHARPS }
   @assert size(df,1) >= 1
   return 1
   # TODO: Replace with something that estimates signal to noise
   #=
   @assert hasproperty(df,:snr_prelim)
   idx = findmax(df.snr_prelim)[2]
   @assert 1 <= idx <= size(df,1)
   return idx
   =#
end

max_stars_repo_name: RvSpectML/EchelleInstruments

max_issues_repo_name: RvSpectML/EchelleInstruments

max_forks_repo_name: RvSpectML/EchelleInstruments


text:
"""
    scale_data!(data; <keyword arguments>)

Scale lifetime and cost data.

See `_scale_time_data!`, `_scale_operational_cost_data!` and `_scale_investment_cost_data!`.

# Arguments
- `data`: a single-network data dictionary.
- `number_of_hours`: number of optimization periods (default: `dim_length(data, :hour)`).
- `year_scale_factor`: how many years a representative year should represent (default: `dim_meta(data, :year, "scale_factor")`).
- `number_of_years`: number of representative years (default: `dim_length(data, :year)`).
- `year_idx`: id of the representative year (default: `1`).
- `cost_scale_factor`: scale factor for all costs (default: `1.0`).
"""
function scale_data!(
        data::Dict{String,Any};
        number_of_hours::Int = haskey(data, "dim") ? dim_length(data, :hour) : 1,
        year_scale_factor::Int = haskey(data, "dim") ? dim_meta(data, :year, "scale_factor") : 1,
        number_of_years::Int = haskey(data, "dim") ? dim_length(data, :year) : 1,
        year_idx::Int = 1,
        cost_scale_factor::Float64 = 1.0
    )
    if _IM.ismultinetwork(data)
        Memento.error(_LOGGER, "`scale_data!` can only be applied to single-network data dictionaries.")
    end
    _scale_time_data!(data, year_scale_factor)
    _scale_operational_cost_data!(data, number_of_hours, year_scale_factor, cost_scale_factor)
    _scale_investment_cost_data!(data, number_of_years, year_idx, cost_scale_factor) # Must be called after `_scale_time_data!`
end

"""
    _scale_time_data!(data, year_scale_factor)

Scale lifetime data from years to periods of `year_scale_factor` years.

After applying this function, the step between consecutive years takes the value 1: in this
way it is easier to write the constraints that link variables belonging to different years.
"""
function _scale_time_data!(data, year_scale_factor)
    rescale = x -> x ÷ year_scale_factor
    for component in ("ne_branch", "branchdc_ne", "ne_storage", "convdc_ne", "load")
        for (key, val) in get(data, component, Dict{String,Any}())
            if !haskey(val, "lifetime")
                if component == "load"
                    continue # "lifetime" field is not used in OPF
                else
                    Memento.error(_LOGGER, "Missing `lifetime` key in `$component` $key.")
                end
            end
            if val["lifetime"] % year_scale_factor != 0
                Memento.error(_LOGGER, "Lifetime of $component $key ($(val["lifetime"])) must be a multiple of the year scale factor ($year_scale_factor).")
            end
            _PM._apply_func!(val, "lifetime", rescale)
        end
    end
end

"""
    _scale_operational_cost_data!(data, number_of_hours, year_scale_factor, cost_scale_factor)

Scale hourly costs to the planning horizon.

Scale hourly costs so that the sum of the costs over all optimization periods
(`number_of_hours` hours) represents the cost over the entire planning horizon
(`year_scale_factor` years). In this way it is possible to perform the optimization using a
reduced number of hours and still obtain a cost that approximates the cost that would be
obtained if 8760 hours were used for each year.
"""
function _scale_operational_cost_data!(data, number_of_hours, year_scale_factor, cost_scale_factor)
    rescale = x -> (8760*year_scale_factor / number_of_hours) * cost_scale_factor * x # scale hourly costs to the planning horizon
    for (g, gen) in data["gen"]
        _PM._apply_func!(gen, "cost", rescale)
    end
    for (l, load) in data["load"]
        _PM._apply_func!(load, "cost_shift", rescale) # Compensation for demand shifting
        _PM._apply_func!(load, "cost_curt", rescale)  # Compensation for load curtailment (i.e. involuntary demand reduction)
        _PM._apply_func!(load, "cost_red", rescale)   # Compensation for not consumed energy (i.e. voluntary demand reduction)
    end
    _PM._apply_func!(data, "co2_emission_cost", rescale)
end

"""
    _scale_investment_cost_data!(data, number_of_years, year_idx, cost_scale_factor)

Correct investment costs considering the residual value at the end of the planning horizon.

Linear depreciation is assumed.

This function _must_ be called after `_scale_time_data!`.
"""
function _scale_investment_cost_data!(data, number_of_years, year_idx, cost_scale_factor)
    # Assumption: the `lifetime` parameter of investment candidates has already been scaled
    # using `_scale_time_data!`.
    remaining_years = number_of_years - year_idx + 1
    for (b, branch) in get(data, "ne_branch", Dict{String,Any}())
        rescale = x -> min(remaining_years/branch["lifetime"], 1.0) * cost_scale_factor * x
        _PM._apply_func!(branch, "construction_cost", rescale)
        _PM._apply_func!(branch, "co2_cost", rescale)
    end
    for (b, branch) in get(data, "branchdc_ne", Dict{String,Any}())
        rescale = x -> min(remaining_years/branch["lifetime"], 1.0) * cost_scale_factor * x
        _PM._apply_func!(branch, "cost", rescale)
        _PM._apply_func!(branch, "co2_cost", rescale)
    end
    for (c, conv) in get(data, "convdc_ne", Dict{String,Any}())
        rescale = x -> min(remaining_years/conv["lifetime"], 1.0) * cost_scale_factor * x
        _PM._apply_func!(conv, "cost", rescale)
        _PM._apply_func!(conv, "co2_cost", rescale)
    end
    for (s, strg) in get(data, "ne_storage", Dict{String,Any}())
        rescale = x -> min(remaining_years/strg["lifetime"], 1.0) * cost_scale_factor * x
        _PM._apply_func!(strg, "eq_cost", rescale)
        _PM._apply_func!(strg, "inst_cost", rescale)
        _PM._apply_func!(strg, "co2_cost", rescale)
    end
    for (l, load) in data["load"]
        rescale = x -> min(remaining_years/load["lifetime"], 1.0) * cost_scale_factor * x
        _PM._apply_func!(load, "cost_inv", rescale)
        _PM._apply_func!(load, "co2_cost", rescale)
    end
end

ext: jl


text:
function clear_marketing(mdl::Model)

    Ψ = mdl.Ψ
    R,Rᴸ = mdl.Prices[ 1:2 ]
    T = mdl.DR
    ω, τ = mdl.AddPrices[ 1:2 ]
    grids = mdl.GridDef.Grids[ 1 : (length(mdl.GridDef.Grids)÷2) ]
    push!(grids,mdl.GridDef.Grids[end])
    grid_names = mdl.GridDef.VarNames[ 1 : (length(mdl.GridDef.Grids)÷2) ]
    grid_names = (grid_names...,mdl.GridDef.VarNames[end])
    (d₁, l₁, λᵢ) = (1, 2, 3)
    #     #Aggregate
    Lᴰ = 0.0
    L= 0.0 # Debt before defaults
    V = 0.0 #Claimed wages
    Vᴰ = 0.0 #Claimed deposits
    D = 0.0 # Deposits before defaults

    cube = CartesianIndices( tuple( length.(grids)...))

    for idx = eachindex(T)

        if (T[idx] == 0)

            Lᴰ += Ψ[idx] * grids[l₁][ cube[idx][l₁] ]
            Vᴰ  += Ψ[idx] * grids[d₁][ cube[idx][d₁] ]
            V  += Ψ[idx] * (grids[λᵢ][ cube[idx][λᵢ] ] == 1) * ( ( ( 1 - τ ) * ω ) - ϐ)

        else

            #Tutaj pomyślec czy nie powinno być  T[idx] - sprawdzić z równaniami poniżej
            L  += Ψ[idx] * grids[l₁][ cube[T[idx]][l₁] ]
            D  += Ψ[idx] * grids[d₁][ cube[T[idx]][d₁] ]

        end

    end

    V += Vᴰ

    B = D - (L - Lᴰ) - V

    #     # Solve et Coalgua
    function sys_of_eqs(Rates::Vector{T}) where T<:Real

        R, Rᴸ = Rates
        F = zeros(2)
        # Commercial Banks:
        F[1] = (L - Lᴰ) * Rᴸ  +  ( Lₛ * ( ρ - 1 )  +  B  -  ( D - Vᴰ )*ρ ) * R  + V
        # Central Bank:
        F[2] = (ρ * ( R/𝜋 ) ) * ( D - Vᴰ )  +  ((1 + R * (1-ρ) )^(-1))  -  (L - Lᴰ) * δᶠ * ( Rᴸ/𝜋 )  +  (L - (L - Lᴰ) * (1 - δᶠ) ) - (V - Vᴰ) - D - 1

        return(F)

    end

    function mx_of_J(Rates::Vector{T}) where T<:Real

        R, Rᴸ = Rates
        J = zeros(2,2)
        J[1,1] = ( Lₛ * ( ρ - 1 )  +  B  -  ( D - Vᴰ )*ρ )
        J[2,1] = (ρ /𝜋)  * D - ( R^(-2) )
        J[1,2] = (ρ * ( 1/𝜋 ) ) * ( D - Vᴰ )  -  (1 - ρ) / ((1 + R * (1-ρ) )^2)
        J[2,2] = (L - Lᴰ) * δᶠ * ( 1/𝜋 )
        return(J)

    end

    # Rates = nlsolve(sys_of_eqs, mx_of_J, [R,Rᴸ], method = :newton).zero
    Rates = nlsolve(sys_of_eqs, [R,Rᴸ], method = :newton).zero
    return(Rates)

end


text:
"""
Function transforms ModelInit or Model into Model with empty properties
"""
function init_model(mdl::ModelInit)

    VF = zeros(Float64, *(length.(mdl.GridDef.Grids[1:Int64(floor(length(mdl.GridDef.Grids)/2))])...), length(mdl.GridDef.Grids[end]) )
    DR = Int64.(VF)
    # TM = Array{Int64,1}()
    PM = VF
    hstry = Vector{NTuple{length(mdl.Prices)*2, Float64}}()

    mdl_new = Model(  mdl.Prices
    , mdl.GetAP(mdl)
    , mdl.GetAP
    , mdl.GridDef
    , mdl.λ
    , mdl.β
    # , mdl.statvar
    , mdl.HAutility
    , Utils()
    , mdl.HAconstraint
    , mdl.HAconstraint_Alt
    # , mdl.VFstart
    , VF
    , DR
    , zeros(Float64,length(mdl.GridDef.Grids[end]))
    , mdl.VF_afterprocess
    , VF
    # , mdl.Ψᵢₙᵢₜ
    , Vector{Float64}()
    , mdl.alternate_return
    , mdl.MCcondition
    , hstry
    , mdl.InitParams
    )

    return(mdl_new)

end
function init_model(mdl::Model)

    #TODO: make arrays with already properly defined dimmensions. TOdo after solver is moreless done.
    # utils = Array{Float64}()
    VF = zeros(Float64, *(length.(mdl.GridDef.Grids[1:Int64(floor(length(mdl.GridDef.Grids)/2))])...), length(mdl.GridDef.Grids[end]) )
    DR = Int64.(VF)
    # TM = Array{Int64,1}()
    PM = VF
    hstry = Vector{NTuple{length(mdl.Prices)*2, Float64}}()

    mdl_new = Model(  mdl.Prices
    , mdl.GetAP(mdl)
    , mdl.GetAP
    , mdl.GridDef
    , mdl.λ
    , mdl.β
    # , mdl.statvar
    , mdl.HAutility
    , Utils()
    , mdl.HAconstraint
    , mdl.HAconstraint_Alt
    # , mdl.VFstart
    , VF
    , DR
    , zeros(Float64,length(mdl.GridDef.Grids[end]))
    , mdl.VF_afterprocess
    , VF
    # , mdl.Ψᵢₙᵢₜ
    , Vector{Float64}()
    , mdl.alternate_return
    , mdl.MCcondition
    , hstry
    , mdl.InitParams
    )

    return(mdl_new)

end

"""
Function calculates convergeance objective
# Arguments
- `prices::Vector{Float64}`: Vector containing prices
- `mdl::Model`: Model struct.
- `method_cube::String`: Method to use while calculating Utils.
"""
function objective_function(prices::Vector{Float64}, mdl::Model, method_cube::String)

    #Initialize and start iteration
    printstyled( "(I) ",bold = true, color = 159)
    printstyled( string("?? Old: ", prices), color = 226 )
    mdl, new_prices = iterate_mdl( prices, mdl, method_cube )

    #Update message within REPL
    stl = length( string("?? Old: ", prices) )
    printstyled( string("\e[", stl,"C", " New: ", new_prices), color = 208)
    message = string(" Objective Function Value: ",sum(abs.(100 * (prices - new_prices))))
    message *="\n"

    #TODO: change message updates by adding save cursosr location/load cursosr location.
    #Will make future updates much more easier

    printstyled(message,color = 63, bold = true )
    return( (1 + Ψ_punish_degree * Ψ_punish + DR_punish_degree * DR_punish) * sum(abs.(100 * (prices - new_prices))))

end

"""
Single iteration of Equilibrium, for given prices.
# Arguments
- `prices::Vector{Float64}`: Vector containing prices
- `mdl::Model`: Model struct.
- `method_cube::String`: Method to use while calculating Utils.
"""
function iterate_mdl(prices::Vector{Float64}, mdl::Model, method_cube::String)

    mdl.Prices = prices
    global DR_punish = false
    global Ψ_punish = false

    # Initialize
    mdl = init_model(mdl)
    printstyled( "\e[1G(C) ", bold = true, color = 20)

    # Get Utils
    mdl.Utils = execute_cube(mdl, method_cube)
    printstyled( "\e[1G(D) ", bold = true, color = 9)

    # Get decission rules
    iterate_decission_rules(mdl)
    printstyled( "\e[1G(Ψ) ", bold = true, color = 196)

    # Get distribution Ψ
    Ψ_search(mdl)

    # Update prices
    new_prices = copy(mdl.MCcondition(mdl))
    printstyled( "\e[1G(U) ", bold = true, color = 46)
    push!( mdl.History, tuple(mdl.Prices..., new_prices...) )

    return(mdl, new_prices )

end

"""
Function modifies Grid Definition based on `grid_power` variable. Grid power is understood as desired maximum number of positions on the grid. Two state variables with grid density 32 each have power 32`*`32=1024.
Decision rules (mdl.DR) are analysed and grid values that are never used are recylced and allocated to space where two or more neighbouring grid points are utilized. Spaces with more neighbouring points are prioritized when assigning points. Number of points is based on `grid_power`.

Returns __`GridDef`__ object

Requires total remake

# Arguments
- `mdl::Model`: Model with DR
- `grid_power::Int`: GridDef struct.
# Examples
```julia-repl
julia> update_grid(my_mdl,2048)
```
"""
function update_grid(mdl::Model,grid_power::Int=*(length.(mdl.GridDef.Grids[1:(end÷2)])...) )

    dr_temp = vec([1])

    while true

        check = length(dr_temp)
        dr_temp = append!(dr_temp, vec(mdl.DR[dr_temp,:]) )
        dr_temp = unique(dr_temp[dr_temp .!= 0])

        (check == length(dr_temp)) && (break)

    end

    ln_dim = length.(mdl.GridDef.Grids[1:(end÷2)])
    dr_ci = CartesianIndices(tuple(ln_dim...))[dr_temp]
    grid_n = length(ln_dim)
    new_grids = copy(mdl.GridDef.Grids)
    ind_in = (mdl.DR[mdl.Ψ .!= 0])
    ind_in = ind_in[ ind_in .!= 0 ]
    dr_to = CartesianIndices( tuple(ln_dim...) )[ind_in]
    ci_idx = Array{Vector{Int}}(undef,grid_n,4)

    for i in 1:grid_n

        i_ci = sort(unique(cat(collect.(Tuple.(dr_ci))...,dims=2)[i,:]))
        check_vec = (i_ci[2:end] .== (i_ci[1:end-1] .+ 1))
        idx_pos, seq_ln =  sequence(check_vec)
        seq_ln = seq_ln[check_vec[idx_pos] .== 1]
        idx_pos = idx_pos[check_vec[idx_pos] .== 1]

        # check if there are any sequences of data points AND if there is only one that conatain on last point and point before

        check_major = (!(isempty(seq_ln) | isempty(idx_pos))) & !( ( (i_ci[idx_pos[end]] + seq_ln[end]) == ln_dim[i]) & ( seq_ln[end] == 1 ) & (length(idx_pos) == 1))
        if  check_major

            #Check if last sequence icludes end point of grid

            if (i_ci[idx_pos[end]] + seq_ln[end]) == ln_dim[i]

                #check if last sequence is simply two points neighbouring each other

                if ( seq_ln[end] == 1 )

                    seq_ln = seq_ln[1:end-1]
                    idx_pos = idx_pos[1:end-1]

                else

                    seq_ln[end] -=1

                end

            end

            if (i_ci[idx_pos[1]]) == 1

                if ( seq_ln[end] == 1 )

                    seq_ln = seq_ln[2:end]
                    idx_pos = idx_pos[2:end]

                else

                    seq_ln[1] -=1
                    idx_pos[1] +=1

                end

            end

            idx_out = Vector{Int}()
            tmp = [range(x,length = y+1) for (x,y) = zip(i_ci[idx_pos], seq_ln )]
            tmp = collect.(tmp)
            append!.([idx_out],tmp)
            idx_out = idx_out[idx_out .∉  [[ 1,length(mdl.GridDef.Grids[2]) ]] ]
            (length(mdl.GridDef.Grids[i]) ∈ i_ci) ? (idx_leave = i_ci) : (idx_leave = [i_ci...,length(mdl.GridDef.Grids[i])])
            idx_leave = idx_leave[idx_leave .∉ [idx_out]]
            ci_idx[i,:] = [idx_leave, idx_out, i_ci[idx_pos], seq_ln]

        else

            ci_idx[i,:] = [ i_ci, [], [], [] ]

        end

    end

    # isdefined(ci_idx)

    if !all(length.(ci_idx[:,1]) .== ln_dim)

        tmp = sqrt( *(sum(length.(ci_idx[:,1:2]), dims = 2)...) \ grid_power )

        new_ln_dim = Int.(round.(sum(length.(ci_idx[:,1:2]) .* tmp, dims = 2)))

        for i in 1:grid_n

            to_i = sort(unique(cat(collect.(Tuple.(dr_to))...,dims=2)[i,:]))
            ci_i = ci_idx[i,:]
            new_grid_i = mdl.GridDef.Grids[i][ci_i[1]]

            if !isempty(ci_idx[i,2])

                points_left = new_ln_dim[i] - (length(ci_i[1]) + sum(ci_i[4] .+ 1))
                special = 0.0

                iter = 0
                bonus_Ψ  = [any(ci_i[3][m] .<=to_i .<= ci_i[3][m]+ci_i[4][m]) for m=1:length(ci_i[3]) ]
                maxiter = false

                let n_points

                    while (true)

                        # points_distro = (ci_i[4] .+ 1 .+special) .* (bonus_Ψ .+ 1))

                        points_distro = (ci_i[4] .+ 1 .+special) .* (bonus_Ψ)
                        n_points = round.(Int, ( (points_left * points_distro ./ sum(points_distro) ) ) )
                        ((sum(n_points) == points_left) | (maxiter)) && break
                        special -= .01 * min(round.(Int,points_left * (ci_i[4] .+ 1 .+special))... )/new_ln_dim[i]
                        iter += 1
                        maxiter = iter == GridDef_maxiter

                    end

                    for j = length.(ci_i[3])

                        k = ci_i[3][j]
                        l = ci_i[4][j]
                        range_i = mdl.GridDef.Grids[i][[k,k+l]]
                        midpoint = sum(range_i)/2
                        g_i = diff(range_i)/(l)
                        n_g_i = diff(range_i) / (n_points[j] + l - 1 * (n_points[j] != 0) )
                        n_grid_part = cumsum(repeat(n_g_i,  n_points[j] + l + 1))
                        n_grid_part .+= midpoint -(sum(n_grid_part) / length(n_grid_part))
                        append!(new_grid_i,n_grid_part)

                    end

                end

            end

            new_grids[[i,i+grid_n]] = repeat([sort(new_grid_i)],2)

        end

    end

    return(GridDef(mdl.GridDef.VarNames, new_grids))

end


text:
abstract internal_AbstractNLPEvaluator;

include("class_variables.jl")
include("class_nlp_cache.jl")
include("class_quadratic_program.jl")
#include("class_qcqp.jl")
include("class_nlp.jl")
include("class_nlo.jl")
include("class_settings.jl")

function validate_dimensions( nlp_eval::internal_AbstractNLPEvaluator, vars::class_variables )
	try
    @assert(n(vars) == n(nlp_eval))
    @assert(m(vars) == m(nlp_eval))
 	catch e
	  println("ERROR validate_dimensions")
		throw(e)
	end
end

function validate_dimensions( nlp_val::class_nlp_cache )
    try
        num_vars = nlp_val.n
        num_constraints = nlp_val.m

        @assert( size( nlp_val.val_hesslag_prod ) == (num_vars,  num_vars) )
        @assert( size( nlp_val.val_jac_a ) == ( num_constraints,  num_vars) )
        @assert( size( nlp_val.val_a ) == ( num_constraints ,))
        @assert( size( nlp_val.val_gradlag ) == ( num_vars,))
        @assert ( size( nlp_val.val_gradc ) == ( num_vars,) )
        @assert ( size( nlp_val.val_b ) == ( num_constraints ,) )
    catch e
        println("ERROR validate_dimensions")
        throw(e)
	  end
end

max_forks_repo_name: ohinder/Everest


text:
"""
harrison.jl

Harrison's tight-binding approach.

References:

* W.A. Harrison: Elementary Electronic Structure Revised Edition (2004).
* L. Shi and A. Papaconstantopoulos: Physical Review B, 70 (2004) 205101.
"""

const η_original = (-16.2, 8.75, -2.39)
const η_modified = (-21.22, 12.60, -2.29)
const η_dict = Dict(:original => η_original, :modified => η_modified)

const γ = 12 # coordination number

immutable WHTB <: TBPerturbation
  version::Symbol # :original or :modified parameters
  zd::Vector{Float64} # number of d-electrons
  rd::Vector{Float64} # d-state radius
  c::Vector{Float64} # composition
end

function WHTB(zd::Real, rd::Real; version=:original)
  WHTB(version, [zd], [rd], [1.0])
end

function WHTB(zd::Vector, rd::Vector, c::Vector; version=:original)
  WHTB(version, zd, rd, c)
end

ncomp(whtb::WHTB)::Int = length(whtb.zd)

function coreradius(whtb::WHTB)::Array{Float64,2}
  N = ncomp(whtb)
  zeros(N,N)
end

function cutoffradius(whtb::WHTB)::Array{Float64,2}
  @attach(whtb, rd)
  N::Int = ncomp(whtb)
  [ 10 * (rd[i] + rd[j]) / 2 for i in 1:N, j in 1:N ]
end

function transfercoefficient(whtb::WHTB)::Float64
  (η_ddσ, η_ddπ, η_ddδ) =
    get(η_dict, whtb.version, zeros(3))
  ( (η_ddσ^2 + 2η_ddπ^2 + 2η_ddδ^2) / 5 )^(1/2)
end

function transfermatrixelement(whtb::WHTB)::Array{Function,2}
  @attach(whtb, rd)
  N::Int = ncomp(whtb)
  A::Float64 = transfercoefficient(whtb)
  [ V(r::Float64)::Float64 = A * (rd[i]*rd[j])^(3/2) / r^5
    for i in 1:N, j in 1:N ]
end

function overlapcoefficient(whtb::WHTB)::Float64
  (η_ddσ, η_ddπ, η_ddδ) =
    get(η_dict, whtb.version, zeros(3))

  (σ_ddσ, σ_ddπ, σ_ddδ) = (5/π, -5/π, 5/2π)

  -2/5 * (σ_ddσ*η_ddσ + 2σ_ddπ*η_ddπ + 2σ_ddδ*η_ddδ)
end

function overlapmatrixelement(whtb::WHTB)::Array{Function,2}
  N::Int = ncomp(whtb)
  rd::Vector{Float64} = whtb.rd
  B::Float64 = overlapcoefficient(whtb)

  [ SV(r::Float64)::Float64 = B * (rd[i]*rd[j])^3 / r^8
    for i in 1:N, j in 1:N ]
end

function bandwidth(whtb::WHTB, ref::ReferenceSystem)::Float64
  @attach(whtb, zd, rd, c)

  N::Int = ncomp(ref)
  ρ::Float64 = totalnumberdensity(ref)
  g::Array{Function,2} = paircorrelation(ref)
  rmin::Array{Float64,2} = emptyradius(ref)
  rcut::Array{Float64,2} = cutoffradius(whtb)
  V::Array{Function,2} = transfermatrixelement(whtb)

  Wd²::Float64 = 0

  for i in 1:N, j in 1:N
    i > j && continue
    n = i == j ? 1 : 2

    gs = spline(g[i,j], rmin[i,j], rcut[i,j], 32)
    f(r) = V[i,j](r)^2 * gs(r) * r^2

    Wd² += 12 * 4π*ρ * n*c[i]*c[j] * ∫(f, rmin[i,j], rcut[i,j])
  end

  if Wd² < 0
    Wd = 0 # sometimes happens
  else
    Wd = √Wd²
  end
end

function pairpotential(whtb::WHTB, A::Float64, B::Float64, n::Int, m::Int)::Matrix{Function}
  @attach(whtb, zd, rd, c)
  N::Int = ncomp(whtb)
  ret = Matrix{Pairpotential}(N,N)
  for i in 1:N, j in 1:N
    i > j && continue

    z̄d = (zd[i] + zd[j]) / 2
    r̄d = √(rd[i]*rd[j])

    u(r) = -A * sqrt(12/γ) * z̄d * (1 - z̄d/10) * r̄d^3 / r^n +
           B * z̄d * r̄d^6 / r^m

    ret[i,j] = ret[j,i] = u
  end
  ret
end

function pairpotential(whtb::WHTB)::Matrix{Function}
  A = transfercoefficient(whtb)
  B = overlapcoefficient(whtb)
  pairpotential(whtb, A, B, 5, 8)
end

function pairpotential_attractive(whtb::WHTB)::Matrix{Function}
  A = transfercoefficient(whtb)
  B = 0.
  pairpotential(whtb, A, B, 5, 8)
end

function pairpotential_repulsive(whtb::WHTB)::Matrix{Function}
  A = 0.
  B = overlapcoefficient(whtb)
  pairpotential(whtb, A, B, 5, 8)
end

function pairpotential_derivative(whtb::WHTB)::Matrix{Function}
  @attach(whtb, zd, rd, c)
  N = ncomp(whtb)
  A = transfercoefficient(whtb)
  B = overlapcoefficient(whtb)
  pairpotential(whtb, -5A, -8B, 6, 9)
end

function entropy(whtb::WHTB, Wd::Float64, T::Float64)::Float64
  nd = 10 / Wd
  S_el = π^2 / 3 * kB^2 * T * nd / kB
end

function entropy(whtb::WHTB, ref::ReferenceSystem, T::Float64)::Float64
  entropy(whtb, bandwidth(whtb, ref), T)
end

function internal_band(whtb::WHTB, ref::ReferenceSystem, Wd::Float64)::Float64
  @attach(whtb, zd, rd, c)
  z̄d = sum(zd .* c)
  Eb = -1/2 * z̄d * (10 - z̄d)/10 * Wd
end

function internal_band(whtb::WHTB, ref::ReferenceSystem)::Float64
  internal_band(whtb, ref, bandwidth(whtb, ref))
end

function internal_rep(whtb::WHTB, ref::ReferenceSystem)::Float64
  N = ncomp(whtb)
  u_rep = pairpotential_repulsive(whtb)
  internal(ref, whtb, u_rep)
end

function internal(whtb::WHTB, ref::ReferenceSystem)::Float64
  U_band = internal_band(whtb, ref)
  U_rep = internal_rep(whtb, ref)
  U = U_band + U_rep
end


text:
module DoublePendulum

using DataDeps, CSV, DataFrames, MLUtils
using NeuralOperators, Flux
using CUDA, FluxTraining, BSON

function register_double_pendulum_chaotic()
    register(DataDep(
        "DoublePendulumChaotic",
        """
        Dataset was generated on the basis of 21 individual runs of a double pendulum.
        Each of the recorded sequences lasted around 40s and consisted of around 17500 frames.

        * `x_red`: Horizontal pixel coordinate of the red point (the central pivot to the first pendulum)
        * `y_red`: Vertical pixel coordinate of the red point (the central pivot to the first pendulum)
        * `x_green`: Horizontal pixel coordinate of the green point (the first pendulum)
        * `y_green`: Vertical pixel coordinate of the green point (the first pendulum)
        * `x_blue`: Horizontal pixel coordinate of the blue point (the second pendulum)
        * `y_blue`: Vertical pixel coordinate of the blue point (the second pendulum)

        Page: https://developer.ibm.com/exchanges/data/all/double-pendulum-chaotic/
        """,
        "https://dax-cdn.cdn.appdomain.cloud/dax-double-pendulum-chaotic/2.0.1/double-pendulum-chaotic.tar.gz",
        "4ca743b4b783094693d313ebedc2e8e53cf29821ee8b20abd99f8fb4c0866f8d",
        post_fetch_method=unpack
    ))
end

function get_data(; i=0, n=-1)
    data_path = joinpath(datadep"DoublePendulumChaotic", "original", "dpc_dataset_csv")
    df = CSV.read(
        joinpath(data_path, "$i.csv"),
        DataFrame,
        header=[:x_red, :y_red, :x_green, :y_green, :x_blue, :y_blue]
    )
    data = (n < 0) ? collect(Matrix(df)') : collect(Matrix(df)')[:, 1:n]

    return Float32.(data)
end

function preprocess(𝐱; Δt=1, nx=30, ny=30)
    # move red point to (0, 0)
    xs_red, ys_red = 𝐱[1, :], 𝐱[2, :]
    𝐱[3, :] -= xs_red; 𝐱[5, :] -= xs_red
    𝐱[4, :] -= ys_red; 𝐱[6, :] -= ys_red

    # needs only green and blue points
    𝐱 = reshape(𝐱[3:6, 1:Δt:end], 1, 4, :)
    # velocity of green and blue points
    ∇𝐱 = 𝐱[:, :, 2:end] - 𝐱[:, :, 1:(end-1)]
    # merge info of pos and velocity
    𝐱 = cat(𝐱[:, :, 1:(end-1)], ∇𝐱, dims=1)

    # with info of first nx steps to inference next ny steps
    n = size(𝐱)[end] - (nx + ny) + 1
    𝐱s = Array{Float32}(undef, size(𝐱)[1:2]..., nx, n)
    𝐲s = Array{Float32}(undef, size(𝐱)[1:2]..., ny, n)
    for i in 1:n
        𝐱s[:, :, :, i] .= 𝐱[:, :, i:(i+nx-1)]
        𝐲s[:, :, :, i] .= 𝐱[:, :, (i+nx):(i+nx+ny-1)]
    end

    return 𝐱s, 𝐲s
end

function get_dataloader(; n_file=20, Δt=1, nx=30, ny=30, ratio=0.9, batchsize=100)
    𝐱s, 𝐲s = Array{Float32}(undef, 2, 4, nx, 0), Array{Float32}(undef, 2, 4, ny, 0)
    for i in 0:(n_file-1)
        𝐱s_i, 𝐲s_i = preprocess(get_data(i=i), Δt=Δt, nx=nx, ny=ny)
        𝐱s, 𝐲s = cat(𝐱s, 𝐱s_i, dims=4), cat(𝐲s, 𝐲s_i, dims=4)
    end

    data = shuffleobs((𝐱s, 𝐲s))
    data_train, data_test = splitobs(data, at=ratio)

    loader_train = Flux.DataLoader(data_train, batchsize=batchsize, shuffle=true)
    loader_test = Flux.DataLoader(data_test, batchsize=batchsize, shuffle=false)

    return loader_train, loader_test
end

__init__() = register_double_pendulum_chaotic()

function train(; cuda=true, Δt=1, η₀=1f-3, λ=1f-4, epochs=20)
    if cuda && CUDA.has_cuda()
        device = gpu
        CUDA.allowscalar(false)
        @info "Training on GPU"
    else
        device = cpu
        @info "Training on CPU"
    end

    model = FourierNeuralOperator(ch=(2, 64, 64, 64, 64, 64, 128, 2), modes=(4, 16), σ=gelu)
    data = get_dataloader(Δt=Δt)
    optimiser = Flux.Optimiser(WeightDecay(λ), Flux.ADAM(η₀))
    loss_func = l₂loss

    learner = Learner(
        model, data, optimiser, loss_func,
        ToDevice(device, device),
        Checkpointer(joinpath(@__DIR__, "../model/"))
    )

    fit!(learner, epochs)

    return learner
end

function get_model()
    model_path = joinpath(@__DIR__, "../model/")
    model_file = readdir(model_path)[end]

    return BSON.load(joinpath(model_path, model_file), @__MODULE__)[:model]
end

end # module


text:
using BenchmarkTools

fib(n) = n <= 1 ?  1 : fib(n - 2) + fib(n - 1)

suite = BenchmarkGroup()

suite["fib"] = BenchmarkGroup(["tag1", "tag2"])

suite["fib"][10] = @benchmarkable fib(10)
suite["fib"][20] = @benchmarkable fib(20)

tune!(suite)
results = run(suite, verbose = true)

BenchmarkTools.save("output.json", median(results))

max_stars_repo_name: mjochum/github-action-benchmark

max_issues_repo_name: mjochum/github-action-benchmark

max_forks_repo_name: mjochum/github-action-benchmark


text:
using DashBootstrapComponents, DashHtmlComponents

badges = html_span([
    dbc_badge("Primary", pill=true, color="primary", className="mr-1"),
    dbc_badge("Secondary", pill=true, color="secondary", className="mr-1"),
    dbc_badge("Success", pill=true, color="success", className="mr-1"),
    dbc_badge("Warning", pill=true, color="warning", className="mr-1"),
    dbc_badge("Danger", pill=true, color="danger", className="mr-1"),
    dbc_badge("Info", pill=true, color="info", className="mr-1"),
    dbc_badge("Light", pill=true, color="light", className="mr-1"),
    dbc_badge("Dark", pill=true, color="dark"),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
export step!, dummystep

"""
    step!(model, agent_step!, n::Int = 1)
    step!(model, agent_step!, model_step!, n::Int = 1, agents_first::Bool=true)

Update agents `n` steps according to the stepping function `agent_step!`.
Agents will be activated as specified by the `model.scheduler`.
`model_step!` is triggered _after_ every scheduled agent has acted, unless
the argument `agents_first` is `false` (which then first calls `model_step!` and then
activates the agents).

    step!(model, agent_step!, model_step!, n::Function, agents_first::Bool=true)

In this version `n` is a function.
Then `step!` runs the model until `n(model, s)` returns `true`, where `s` is the
current amount of steps taken, starting from 0.
For this method of `step!`, `model_step!` must be provided always (use [`dummystep`](@ref)
if you have no model stepping dynamics).
"""
function step! end

"""
    dummystep(model)

Ignore the model dynamics. Use instead of `model_step!`.
"""
dummystep(model) = nothing
"""
    dummystep(agent, model)

Ignore the agent dynamics. Use instead of `agent_step!`.
"""
dummystep(agent, model) = nothing

until(ss, n::Int, model) = ss < n
until(ss, n, model) = !n(model, ss)

step!(model::ABM, agent_step!, n::Int=1, agents_first::Bool=true) = step!(model, agent_step!, dummystep, n, agents_first)

function step!(model::ABM, agent_step!, model_step!, n = 1, agents_first=true)
  s = 0
  while until(s, n, model)
    !agents_first && model_step!(model)
    activation_order = model.scheduler(model)
    for index in activation_order
      haskey(model.agents, index) || continue
      agent_step!(model.agents[index], model)
    end
    agents_first && model_step!(model)
    s += 1
  end
end


text:
using GPMLj, Test, Random, PyCall

Random.seed!(123)

@testset "Kernels" begin
    @testset "Matern12" begin
        kern = gpflow.kernels.Matern12(2)
        @test typeof(kern)<:gpflow.Kernel
        @test typeof(kern)<:gpflow.kernels.AbstractStationary
        @test typeof(kern)<:gpflow.kernels.Matern12
        @test typeof(kern.o)<:PyObject
        temp = kern.o
        @test instantiate!(kern) == temp
    end

    @testset "Matern32" begin
        kern = gpflow.kernels.Matern32(2)
        @test typeof(kern)<:gpflow.Kernel
        @test typeof(kern)<:gpflow.kernels.AbstractStationary
        @test typeof(kern)<:gpflow.kernels.Matern32
        @test typeof(kern.o)<:PyObject
        temp = kern.o
        @test instantiate!(kern) == temp
    end

    @testset "Matern52" begin
        kern = gpflow.kernels.Matern52(2)
        @test typeof(kern)<:gpflow.Kernel
        @test typeof(kern)<:gpflow.kernels.AbstractStationary
        @test typeof(kern)<:gpflow.kernels.Matern52
        @test typeof(kern.o)<:PyObject
        temp = kern.o
        @test instantiate!(kern) == temp
    end

    @testset "ArcCosine" begin
        kern = gpflow.kernels.ArcCosine(2)
        @test typeof(kern)<:gpflow.Kernel
        @test typeof(kern)<:gpflow.kernels.ArcCosine
        @test typeof(kern.o)<:PyObject
        temp = kern.o
        @test instantiate!(kern) == temp
    end

    @testset "Periodic" begin
        kern = gpflow.kernels.Periodic(2)
        @test typeof(kern)<:gpflow.Kernel
        @test typeof(kern)<:gpflow.kernels.Periodic
        @test typeof(kern.o)<:PyObject
        temp = kern.o
        @test instantiate!(kern) == temp
    end

    # TODO
    # @testset "Coregion" begin
    #     kern= gpflow.kernels.Coregion(2, 1)
    #     @test typeof(kern)<:gpflow.Kernel
    #     @test typeof(kern)<:gpflow.kernels.Coregion
    #     @test typeof(kern.o)<:PyObject
    #     temp = kern.o
    #     @test instantiate!(kern) == temp
    # end

    # Static Kernels
    @testset "Static" begin
        kern = gpflow.kernels.Static(2)
        @test typeof(kern)<:gpflow.Kernel
        @test typeof(kern)<:gpflow.kernels.AbstractStatic
        @test typeof(kern)<:gpflow.kernels.Static
        @test typeof(kern.o)<:PyObject
        temp = kern.o
        @test instantiate!(kern) == temp
    end

    @testset "Combination" begin
        kern1 = gpflow.kernels.Matern52(2)
        kern2 = gpflow.kernels.ArcCosine(2)
        kern3 = gpflow.kernels.Periodic(2)
        kernels = [kern1, kern2, kern3];
        kern = gpflow.kernels.Combination(kernels)
        @test typeof(kern)<:gpflow.Kernel
        @test typeof(kern)<:gpflow.kernels.AbstractCombination
        @test typeof(kern)<:gpflow.kernels.Combination
        @test typeof(kern.o)<:PyObject
        temp = kern.o
        @test instantiate!(kern) == temp
    end

    @testset "Product" begin
        kern1 = gpflow.kernels.Matern52(2)
        kern2 = gpflow.kernels.ArcCosine(2)
        kern3 = gpflow.kernels.Periodic(2)
        kernels = [kern1, kern2, kern3];
        kern = gpflow.kernels.Product(kernels)
        @test typeof(kern)<:gpflow.Kernel
        @test typeof(kern)<:gpflow.kernels.AbstractCombination
        @test typeof(kern)<:gpflow.kernels.Product
        @test typeof(kern.o)<:PyObject
        temp = kern.o
        @test instantiate!(kern) == temp
    end

    @testset "Sum" begin
        kern1 = gpflow.kernels.Matern52(2)
        kern2 = gpflow.kernels.ArcCosine(2)
        kern3 = gpflow.kernels.Periodic(2)
        kernels = [kern1, kern2, kern3];
        kern = gpflow.kernels.Sum(kernels)
        @test typeof(kern)<:gpflow.Kernel
        @test typeof(kern)<:gpflow.kernels.AbstractCombination
        @test typeof(kern)<:gpflow.kernels.Sum
        @test typeof(kern.o)<:PyObject
        temp = kern.o
        @test instantiate!(kern) == temp
    end

    @testset "Convolutional" begin
        basekern = gpflow.kernels.Matern52(4)
        kern = gpflow.kernels.Convolutional(basekern, [3, 3], [2, 2])
        @test typeof(kern)<:gpflow.Kernel
        @test typeof(kern)<:gpflow.kernels.Convolutional
        @test typeof(kern.o)<:PyObject
        temp = kern.o
        @test instantiate!(kern) == temp
    end

    @testset "WeightedConvolutional" begin
        basekern = gpflow.kernels.Matern52(4)
        kern = gpflow.kernels.Convolutional(basekern, [3, 3], [2, 2])
        @test typeof(kern)<:gpflow.Kernel
        @test typeof(kern)<:gpflow.kernels.Convolutional
        @test typeof(kern.o)<:PyObject
        temp = kern.o
        @test instantiate!(kern) == temp
    end

    @testset "Cosine" begin
        kern = gpflow.kernels.Cosine(2)
        @test typeof(kern)<:gpflow.Kernel
        @test typeof(kern)<:gpflow.kernels.AbstractStationary
        @test typeof(kern)<:gpflow.kernels.Cosine
        @test typeof(kern.o)<:PyObject
        temp = kern.o
        @test instantiate!(kern) == temp
    end

    @testset "Exponential" begin
        kern = gpflow.kernels.Exponential(2)
        @test typeof(kern)<:gpflow.Kernel
        @test typeof(kern)<:gpflow.kernels.AbstractStationary
        @test typeof(kern)<:gpflow.kernels.Exponential
        @test typeof(kern.o)<:PyObject
        temp = kern.o
        @test instantiate!(kern) == temp
    end

    @testset "RationalQuadratic" begin
        kern = gpflow.kernels.RationalQuadratic(2)
        @test typeof(kern)<:gpflow.Kernel
        @test typeof(kern)<:gpflow.kernels.AbstractStationary
        @test typeof(kern)<:gpflow.kernels.RationalQuadratic
        @test typeof(kern.o)<:PyObject
        temp = kern.o
        @test instantiate!(kern) == temp
    end

    @testset "SquaredExponential" begin
        kern = gpflow.kernels.SquaredExponential(2)
        @test typeof(kern)<:gpflow.Kernel
        @test typeof(kern)<:gpflow.kernels.AbstractStationary
        @test typeof(kern)<:gpflow.kernels.SquaredExponential
        @test typeof(kern.o)<:PyObject
        temp = kern.o
        @test instantiate!(kern) == temp
    end

    @testset "Stationary" begin
        kern = gpflow.kernels.Stationary(2)
        @test typeof(kern)<:gpflow.Kernel
        @test typeof(kern)<:gpflow.kernels.AbstractStationary
        @test typeof(kern)<:gpflow.kernels.Stationary
        @test typeof(kern.o)<:PyObject
        temp = kern.o
        @test instantiate!(kern) == temp
    end

    @testset "Linear" begin
        kern = gpflow.kernels.Linear(2)
        @test typeof(kern)<:gpflow.Kernel
        @test typeof(kern)<:gpflow.kernels.AbstractLinear
        @test typeof(kern)<:gpflow.kernels.Linear
        @test typeof(kern.o)<:PyObject
        temp = kern.o
        @test instantiate!(kern) == temp
    end

    @testset "Polynomial" begin
        kern = gpflow.kernels.Polynomial(2)
        @test typeof(kern)<:gpflow.Kernel
        @test typeof(kern)<:gpflow.kernels.AbstractLinear
        @test typeof(kern)<:gpflow.kernels.Polynomial
        @test typeof(kern.o)<:PyObject
        temp = kern.o
        @test instantiate!(kern) == temp
    end
end #module


text:
## --- Geochronology.jl

    @test eHf(0.2818792, 0.001009289, 1424.878) ≈ -0.9639408985107067

    (c, m, el, eu) = bin_bsr_eHf(1:4000, fill(0.2818792, 4000), fill(0.001009289,4000), 1:4000, 0, 4000, 4, fill(10,4000), fill(0.001,4000), fill(1e-6,4000), fill(10,4000), 500)
    @test c == 500:1000:3500
    isapprox(m, [-21.316, 0.724, 23.187, 45.97], atol=5)
    isapprox(el, [2.254, 2.216, 2.28, 2.319], atol=0.5)
    isapprox(eu, [2.245, 2.216, 2.273, 2.348], atol=0.5)

## --- GIS.jl

    # Calculate slope on a sphere
    A = [8 7 10 -1 -1; 5 2 2 -1 -1; 8 5 3 1 6; 8 10 3 9 5; -1 8 4 10 4]
    @test aveslope(A*1000, 1:4, 1:4, 1, Int64, minmatval=0, km_per_lat=111.1) ==
        [0 20 0 0 0; 19 23 8 0 0; 18 15 11 22 15; 4 0 15 35 16; 0 14 15 0 0]
    @test aveslope(A*1000, 1:4, 1:4, 1, Float16, minmatval=0, km_per_lat=111.1) ==
        Float16[0.0 19.72 0.0 0.0 0.0; 18.75 22.72 8.29 0.0 0.0; 18.45 15.375 11.2 22.45 14.63; 3.824 0.0 15.06 35.4 15.914; 0.0 13.59 15.39 0.0 0.0]
    @test maxslope(A*1000, 1:4, 1:4, 1, minmatval=0, km_per_lat=111.1) ==
        [0 45 72 0 0; 27 16 32 0 6; 38 36 22 14 45; 18 23 22 41 36; 0 32 38 45 0]
    @test maxslope(A*1000, 1:4, 1:4, 1, Float16, minmatval=0, km_per_lat=111.1) ==
        Float16[0.0 45.0 72.0 0.0 0.0; 27.02 15.914 31.5 0.0 6.367; 38.22 36.0 22.3 13.52 45.06; 18.05 22.56 22.3 40.5 36.1; 0.0 31.86 38.22 44.6 0.0]

    # Test ESRI arc/info ASCII grid import function
    grid = """
    ncols         4
    nrows         6
    xllcorner     0.0
    yllcorner     0.0
    cellsize      50.0
    NODATA_value  -9999
    -9999 -9999 5 2
    -9999 20 100 36
    3 8 35 10
    32 42 50 6
    88 75 27 9
    13 5 1 -9999
    """

    f = open("grid.asc","w")
    print(f, grid)
    close(f)

    (data, metadata) = importAAIGrid("grid.asc", Int64, undefval=-999)
    @test eltype(data) == Int64
    @test data == [-9999 -9999 5 2; -9999 20 100 36; 3 8 35 10; 32 42 50 6; 88 75 27 9; 13 5 1 -9999]
    @test metadata["nodata"] == -9999

    # Random lat-lon generation
    @test isa(randlatlon(), Tuple{Float64,Float64})
    lat,lon = randlatlon(100)
    @test maximum(lat) <= 90
    @test minimum(lat) >= -90
    @test maximum(lon) <= 180
    @test minimum(lon) >= -180

## ---


text:
using Documenter, DiffEqParamEstim

makedocs(
    sitename="DiffEqParamEstim.jl",
    authors="Chris Rackauckas et al.",
    clean=true,
    doctest=false,
    modules=[DiffEqParamEstim],

    format=Documenter.HTML(assets=["assets/favicon.ico"],
                           canonical="https://diffeqparamestim.sciml.ai/stable/"),

    pages=[
        "DiffEqParamEstim.jl: Parameter Estimation for Differential Equations" => "index.md",
        "Methods" => Any[
          "methods/recommended_methods.md",
          "methods/optimization_based_methods.md"
        ],
        "Examples" => Any[
          "examples/ODE_inference.md"
        ],
        "Tutorials" => Any[
          "tutorials/stochastic_evaluations.md"
        ],
    ]
)

deploydocs(
    repo="github.com/SciML/DiffEqParamEstim.jl";
    push_preview=true
)


text:
"""
This function creates the expected market data for each investor for each scenario using CEM.
"""
function create_expected_marketdata(investor_dir::String,
                                    market_names::Vector{Symbol},
                                    expected_portfolio::Vector{<: Project{<: BuildPhase}},
                                    zones::Vector{String},
                                    lines::Vector{ZonalLine},
                                    peak_load::Float64,
                                    rep_hour_weight::Vector{Float64},
                                    average_capital_cost_multiplier::Float64,
                                    scenario::Scenario,
                                    iteration_year::Int64,
                                    yearly_horizon::Int64,
                                    solver::JuMP.MOI.OptimizerWithAttributes,
                                    sys_results_dir::String,
                                    investor_name::String)

    system = create_cem_mkt_clr_problem(investor_dir,
                                        market_names,
                                        expected_portfolio,
                                        zones,
                                        lines,
                                        peak_load,
                                        rep_hour_weight,
                                        average_capital_cost_multiplier,
                                        scenario,
                                        iteration_year,
                                        yearly_horizon)

    capacity_price,
    energy_price,
    reserve_up_price,
    reserve_down_price,
    rec_price,
    capacity_factors,
    capacity_accepted_perc,
    new_options = cem(
                        system,
                        solver,
                        "C:/Users/manwar2/Documents/GitRepos/emt-tests/data/simulation_data/results/cem_results.txt"
                     )

    output_file = joinpath(investor_dir, "expected_market_data", "$(get_name(scenario))_year_$(iteration_year).jld2")

    FileIO.save(output_file,
                     "capacity_price", capacity_price,
                     "energy_price", energy_price,
                     "reserve_up_price", reserve_up_price,
                     "reserve_down_price", reserve_down_price,
                     "rec_price", rec_price,
                     "capacity_factors", capacity_factors,
                     "capacity_accepted_perc", capacity_accepted_perc,
                     "new_options", new_options
        )

    sys_results_file = joinpath(sys_results_dir, investor_name, "expected_market_data", "$(get_name(scenario))_year_$(iteration_year).jld2")

    FileIO.save(sys_results_file,
                     "capacity_price", capacity_price,
                     "energy_price", energy_price,
                     "reserve_up_price", reserve_up_price,
                     "reserve_down_price", reserve_down_price,
                     "rec_price", rec_price,
                     "capacity_factors", capacity_factors,
                     "capacity_accepted_perc", capacity_accepted_perc,
                     "new_options", new_options
        )

    return
end

"""
This function updates the maximum new options which can be built by the investor based on CEM predictions.
"""
function update_max_new_options!(max_new_options::Dict{String, Int64},
                                scenario_new_options::Dict{String, Int64},
                                option_projects::Vector{Project})
    for project in option_projects
        type = get_type(get_tech(project))
        zone = get_zone(get_tech(project))
        if scenario_new_options["option_$(type)_$(zone)"] > max_new_options[get_name(project)]
            max_new_options[get_name(project)] = scenario_new_options["option_$(type)_$(zone)"]
        end
    end
    return max_new_options
end


text:
@testset "997.find-the-town-judge.jl" begin
    @test find_judge(3, [[1,3],[2,3]]) == 3
    @test find_judge(3, [[1,3],[2,3],[3,1]]) == -1
    @test find_judge(3, [[1,2],[2,3]]) == -1
    @test find_judge(4, [[1,3],[1,4],[2,3],[2,4],[4,3]]) == 3
end


text:
#! /usr/bin/julia

# Rosetta Code, Constrained random points on a circle

const LO = 10
const HI = 15
const GOAL = 100

canvas = falses(2HI+1, 2HI+1)
i = 0

while i < GOAL
    x = rand(-HI:HI)
    y = rand(-HI:HI)
    LO^2-1 < x^2 + y^2 < HI^2+1 || continue
    i += 1
    canvas[x+HI+1, y+HI+1] = true
end

for i in 1:(2HI+1)
    println("    ", join(map(j -> j ? "\u25cf " : "  ", canvas[i,:])))
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
using Documenter
using GeCKOReader
using Bio

makedocs(
    format = :html,
    sitename = "GeCKOReader.jl",
    modules = [GeCKOReader],
    pages = [
        "Home" => "index.md"
    ]
)

# deploy off of julia 0.5 instead of nightly
deploydocs(
    repo="github.com/tlnagy/GeCKOReader.jl.git",
    julia  = "0.5",
    deps = nothing,
    make = nothing,
    target = "build"
)


text:
module test_julia_bytes

using Test

u8 = Vector{UInt8}([0x01, 0x02, 0x03, 0x04])
u32 = Vector{UInt32}([0x04030201])

@test reinterpret(UInt32, u8) == u32
@test reinterpret(UInt8, u32) == u8

end # module test_julia_bytes

module test_julia_bytes_b_str

using Test

@test b"hello" == @b_str("hello") == [0x68, 0x65, 0x6c, 0x6c, 0x6f]

end # module test_julia_bytes_b_str


text:
export LevenbergMarquardt_Result

@doc raw"""
```julia
struct LevenbergMarquardt_Result{T<:Real} <: Abstract_Solver_Result
   ...
end
```

This structure subtypes [`Abstract_Solver_Result`](@ref)
"""
struct LevenbergMarquardt_Result{T<:Real} <: Abstract_Solver_Result
    _converged::Bool
    _iter_count::Int
    _fobj::T
    _solution::Vector{T}

    # Do not use Base.@kwdef for Julia 1 compatibility reason
    function LevenbergMarquardt_Result(;
                                       _converged,
                                       _iter_count,
                                       _fobj,
                                       _solution::Vector{T}) where {T<:Real}

        new{T}(_converged, _iter_count, _fobj, _solution)
    end

end
converged(r::LevenbergMarquardt_Result) = r._converged
iteration_count(r::LevenbergMarquardt_Result) = r._iter_count
objective_value(r::LevenbergMarquardt_Result) = r._fobj
solution(r::LevenbergMarquardt_Result) = ReadOnlyArray(r._solution)


text:
using DataFrames, FITSIO, ProgressMeter, Dates

import DataFrames.DataFrame, Base.write

export fits2df

DataFrame(hdu::FITSIO.ImageHDU) = error("Can not convert an ImageHDU extension to a DataFrame")
function DataFrame(hdu::FITSIO.HDU)
    fld = join.(split.(FITSIO.colnames(hdu), '.'), '_')
    out = DataFrame()
		@showprogress 0.5 for i in 1:length(fld)
        tmp = read(hdu, (FITSIO.colnames(hdu))[i])
        if ndims(tmp) == 1
            out[!, Symbol(fld[i])] = tmp
        else
            app = DataFrame(Symbol(fld[i])=>Vector{typeof(tmp[1,1])}[])
            for j in 1:(size(tmp))[2]
                push!(app, [tmp[:,j]])
            end
            out = hcat(out, app)
        end
    end
    return out
end

function fits2df(filename::String, hdu=2)
    f = FITS(filename)
    out = DataFrame(f[hdu])
    close(f)
    return out
end

function write(f::FITSIO.FITS, dfr::DataFrame)
    df = deepcopy(dfr)
    data = Array{Any}(undef, 0)
    for name in names(df)
        tt = eltype(df[:, name])
        if isa(tt, Union)
            tt = nonmissingtype(tt)
            if tt == String
                @warn "Using empty string as missing value in field $name"
                df[ismissing.(df[:, name]), name] .= ""
            elseif tt <: AbstractFloat
                @warn "Using NaN as missing value in field $name"
                df[ismissing.(df[:, name]), name] .= NaN
            elseif tt <: Integer
                @warn "Using zero as missing value in field $name"
                df[ismissing.(df[:, name]), name] .= 0
            end
            disallowmissing!(df, name)
        end
        tt = eltype(df[:, name])
        @assert !isa(tt, Union)

        if tt == Symbol
            push!(data, string.(df[:, name]))
        elseif tt == String
            push!(data, string.(df[:, name]))
        elseif tt <: Dates.AbstractTime
            push!(data, string.(df[:, name]))
        else
            push!(data, df[:, name])
        end
    end
    write(f, string.(names(df)), data)
end


text:
module CorrNoise

using Random

export Flat128RNG, initflatrng128, GaussRNG, Oof2RNG, OofRNG, randoof2, randoof

const scalefactor = 1 / (1 + 0xFFFFFFFF)

function twiddle(var::UInt32)
    for i = 1:9
        var = xor(var, var << 13)
        var = xor(var, var >> 17)
        var = xor(var, var << 5)
    end

    var
end

"""
    Flat128RNG

State of the base-128 uniform random generator. Initialize this using the
function `initflatrng128`.
"""
mutable struct Flat128RNG <: AbstractRNG
    state::Array{UInt32}
end

function nextstate(state::Array{UInt32})
    tmp::UInt32 = xor(state[1], state[1] << 11)
    state[1] = state[2]
    state[2] = state[3]
    state[3] = state[4]
    state[4] = xor(xor(state[4], state[4] >> 19), xor(tmp, tmp >> 8))

    state
end

"""
    initflatrng128(xstart = 123456789, ystart = 362436069, zstart = 521288629, wstart = 88675123)

Initialize a flat random number generator with period 2^128. To draw random numbers,
use the `Base.rand` function as usual.

Example:
```
rng = initflatrng128()
print([rand(rng) for i in 1:4])
```
"""
function initflatrng128(xstart = 123456789, ystart = 362436069, zstart = 521288629, wstart = 88675123)
    rng = Flat128RNG(Array{UInt32}([xstart, ystart, zstart, wstart]))

    # Suffle the bits in the seeds
    for i = 1:4
        rng.state[i] = twiddle(rng.state[i])
    end

    # Burnin the RNG
    for i = 1:17
        rng.state = nextstate(rng.state)
    end

    rng
end

function Base.rand(rng::Flat128RNG)
    result = rng.state[4] * scalefactor
    nextstate(rng.state)
    result
end

# Gaussian distribution

mutable struct GaussRNG
    flatrng
    empty::Bool
    gset
end

"""
    GaussRNG(flatrng::AbstractRNG)
    GaussRNG(seed=0)

Initialize a Gaussian RNG. The parameter `flatrng` must be a uniform RNG.
If a `seed` is used, then a `MersenneTwister` RNG is used.
"""
GaussRNG(flatrng::AbstractRNG) = GaussRNG(flatrng, true, 0)
GaussRNG(seed=0) = GaussRNG(MersenneTwister(seed))

function Base.randn(state::GaussRNG)
    if state.empty
        local v1, v2, rsq

        while true
            v1 = 2 * rand(state.flatrng) - 1.0
            v2 = 2 * rand(state.flatrng) - 1.0
            rsq = v1^2 + v2^2
            if 0 < rsq < 1
                break
            end
        end

        fac = sqrt(-2log(rsq) / rsq)
        state.gset = v1 * fac
        state.empty = false

        v2 * fac
    else
        state.empty = true
        state.gset
    end
end

# 1/f^2 distribution

mutable struct Oof2RNG
    normrng
    c0::Float64
    c1::Float64
    d0::Float64
    x1::Float64
    y1::Float64
end

"""
    Oof2RNG(normrng, fmin, fknee, fsample)

Create a `Oof2RNG` RNG object. It requires a gaussian RNG generator in `normrng` (use `GaussRNG`),
the minimum frequency (longest period) in `fmin`, the knee frequency and the sampling frequency.
The measure unit of the three frequencies must be the same (e.g., Hz).

Use `randoof2` to draw samples from a `Oof2RNG` object, as in the following example:
```
rng = Oof2RNG(GaussRNG(), 1e-3, 1.0, 1e2)
print([randoof2(rng) for i in 1:4])
```
"""
function Oof2RNG(normrng, fmin::Number, fknee::Number, fsample::Number)
    w0 = π * fmin / fsample
    w1 = π * fknee / fsample

    Oof2RNG(normrng,
            (1.0 + w1) / (1.0 + w0),
            (w1 - 1.0) / (1.0 + w0),
            (1.0 - w0) / (1.0 + w0),
            0,
            0)
end

function oof2filter(rng::Oof2RNG, x2::Float64)
    y2 = rng.c0 * x2 + rng.c1 * rng.x1 + rng.d0 * rng.y1
    rng.x1 = x2
    rng.y1 = y2

    y2
end

"""
    randoof2(rng::Oof2RNG)

Draw a random sample from a 1/f^2 distribution.
"""
randoof2(rng::Oof2RNG) = oof2filter(rng, randn(rng.normrng))

# 1/f^alpha distribution

mutable struct OofRNG
    normrng
    slope
    fmin
    fknee
    fsample
    oof2states::Array{Oof2RNG}
end

wminmax(fmin::Number, fknee::Number) = (log10(2 * π * fmin), log10(2 * π * fknee))

function numofpoles(fmin::Number, fknee::Number, fsample::Number)
    (wmin, wmax) = wminmax(fmin, fknee)

    convert(Int32, floor(2(wmax - wmin) + log10(fsample)))
end

const OOF2STATESIZE = 5
oofstatesize(fmin::Number, fknee::Number, fsample::Number) = OOF2STATESIZE * numofpoles(fmin, fknee, fsample)

"""
    OofRNG(normrng, slope, fmin, fknee, fsample)

Create a `OofRNG` RNG object. It requires a gaussian RNG generator in `normrng`
(use `GaussRNG`), the slope α of the noise in `slope`, the minimum frequency
(longest period) in `fmin`, the knee frequency and the sampling frequency. The
measure unit of the three frequencies must be the same (e.g., Hz).

Use `randoof` to draw samples from a `OofRNG` object, as in the following example:
```
rng = OofRNG(GaussRNG(), -1.5, 1e-3, 1.0, 1e2)
print([randoof(rng) for i in 1:4])
```
"""
function OofRNG(normrng, slope::Number, fmin::Number, fknee::Number, fsample::Number)
    (wmin, wmax) = wminmax(fmin, fknee)
    a = -slope
    nproc = numofpoles(fmin, fknee, fsample)
    dp = (wmax - wmin) / nproc

    p = wmin + (1 - a / 2) * dp / 2
    z = p + a * dp / 2

    oof2states = Array{Oof2RNG}(undef, nproc)
    for i = 1:nproc
        oof2states[i] = Oof2RNG(normrng, 10^p / (2π), 10^z / (2π), fsample)

        p += dp
        z = p + a * dp / 2
    end

    OofRNG(normrng, slope, fmin, fknee, fsample, oof2states)
end

"""
    randoof(rng::OofRNG)

Draw a random sample from a 1/f^α distribution.
"""
function randoof(rng::OofRNG)
    x2 = randn(rng.normrng)
    for curstate in rng.oof2states
        x2 = oof2filter(curstate, x2)
    end

    x2
end

end # module


text:
include("preamble.jl")

include("plot_supfunc.jl")

include("plot_overapproximate.jl")

include("plot_polyhedron.jl")

include("plot_VTK.jl")

include("plot_TaylorModel.jl")


text:
using Plasm
using PyCall

p = PyCall.pyimport("pyplasm")
p_VIEW = p["VIEW"]
p_STRUCT = p["STRUCT"]
p_MKPOL = p["MKPOL"]

@testset "mini-FL Tests" begin
   @testset "comp" begin
		@test  Plasm.comp([acos,atan,tan,sin])(pi/4) * 4 - pi < 0.1e-14
		@test  Plasm.comp([acos,atan,tan,sin])(1)==(acos∘atan∘tan∘sin)(1)
		@test  Plasm.comp([atan,tan])(pi/4)*4 > 3.14159265358979
   end

   @testset "cons Tests" begin
		@test ((Plasm.cons([acos,atan,tan,sin])(pi/4) .<=
		[0.667458,0.665774,1.0,0.707107]))==[true,true,true,true]
		@test (Plasm.cons([sin,cos])(pi/2).>=[1.0,0.0])==[true,true]
		@test Plasm.cons([+,-])([2,3])==[[2,3],[-2,-3]]
		@test typeof(Plasm.cons([+,-])([2,3]))==Array{Array{Int64,1},1}
   end

   @testset "k Tests" begin
		@test Plasm.k(10)(100)==10
		@test Plasm.k(10)(sin)==10
		@test Plasm.k(sin)(cos)==sin
		@test Plasm.k(sin)(1000)==sin
		@test Plasm.k([1,2,3])(100)==[1,2,3]
		@test Plasm.k([1,2,3])(sin)==[1,2,3]
   end

   @testset "aa Tests" begin
		@test (Plasm.aa(sin)([0.,pi/2])==[0.0,1.0])
		@test typeof(Plasm.aa(sin)([0.,pi/6,pi/2]))==Array{Float64,1}
		@test Plasm.aa(sum)([[1,2],[3,4],[5,6],[7,8]])==[3,7,11,15]
   end
end

@testset "ascii codes Tests" begin
	@test Plasm.ascii_LAR[32] == ([0.0, 0.0], Array{Int64,1}[[1]])
	@test Plasm.ascii_LAR[33] == ([1.75 1.75 2.0 2.0 1.5 1.5; 1.75 5.5 0.5 1.0 1.0 0.5], Array{Int64,1}[[1, 2], [3, 4], [4, 5], [5, 6], [6, 3]])
	@test Plasm.ascii_LAR[34] == ([1.0 2.0 2.0 2.0 1.5 1.5 2.0 3.0 3.0 3.0 2.5 2.5; 4.0 5.0 5.5 6.0 6.0 5.5 4.0 5.0 5.5 6.0 6.0 5.5], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 3], [7, 8], [8, 9], [9, 10], [10, 11], [11, 12], [12, 9]])
	@test Plasm.ascii_LAR[35] == ([1.0 3.0 1.0 3.0 1.25 1.75 2.25 2.75; 2.5 2.5 3.5 3.5 1.75 4.0 1.75 4.0], Array{Int64,1}[[1, 2], [3, 4], [5, 6], [7, 8]])
	@test Plasm.ascii_LAR[36] == ([0.0 1.0 3.0 4.0 4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 2.0 2.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 4.0 5.0 6.0 6.0 5.0 -0.5 6.5], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10], [10, 11], [11, 12], [13, 14]])
	@test Plasm.ascii_LAR[37] == ([2.5 2.5 2.0 2.0 2.5 2.5 2.0 2.0 1.5 3.5; 0.0 0.5 0.5 0.0 5.5 6.0 6.0 5.5 5.5 11.5], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 1], [5, 6], [6, 7], [7, 8], [8, 5], [9, 10]])
	@test Plasm.ascii_LAR[38] == ([4.0 3.0 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 1.0 1.0 3.0 4.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 4.0 5.0 6.0 5.0 4.0 2.0 2.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10], [10, 11], [11, 12], [12, 13], [13, 14]])
	@test Plasm.ascii_LAR[39] == ([1.0 2.0 2.0 2.0 1.5 1.5; 4.0 5.0 5.5 6.0 6.0 5.5], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 3]])
	@test Plasm.ascii_LAR[40] == ([2.0 1.0 0.5 1.0 2.0; 0.0 1.0 3.0 5.0 6.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5]])
	@test Plasm.ascii_LAR[41] == ([2.0 3.0 3.5 3.0 2.0; 0.0 1.0 3.0 5.0 6.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5]])
	@test Plasm.ascii_LAR[42] == ([1.0 3.0 2.0 2.0 1.0 3.0 1.0 3.0; 3.0 3.0 2.0 4.0 2.0 4.0 4.0 2.0], Array{Int64,1}[[1, 2], [3, 4], [5, 6], [7, 8]])
	@test Plasm.ascii_LAR[43] == ([1.0 3.0 2.0 2.0; 3.0 3.0 2.0 4.0], Array{Int64,1}[[1, 2], [3, 4]])
	@test Plasm.ascii_LAR[44] == ([1.0 2.0 2.0 2.0 1.5 1.5; -1.0 0.0 0.5 1.0 1.0 0.5], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 3]])
	@test Plasm.ascii_LAR[45] == ([1.0 3.0; 3.0 3.0], Array{Int64,1}[[1, 2]])
	@test Plasm.ascii_LAR[46] == ([2.0 2.0 1.5 1.5 2.0; 0.0 0.5 0.5 0.0 0.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5]])
	@test Plasm.ascii_LAR[47] == ([1.0 3.0; 0.0 6.0], Array{Int64,1}[[1, 2]])
	@test Plasm.ascii_LAR[48] == ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0; 1.0 0.0 0.0 1.0 5.0 6.0 6.0 5.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 1], [4, 8]])
	@test Plasm.ascii_LAR[49] == ([0.0 2.0 2.0 0.0 4.0; 4.0 6.0 0.0 0.0 0.0], Array{Int64,1}[[1, 2], [2, 3], [4, 5]])
	@test Plasm.ascii_LAR[50] == ([0.0 0.0 1.0 3.0 4.0 4.0 0.0 4.0; 4.0 5.0 6.0 6.0 5.0 4.0 0.0 0.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8]])
	@test Plasm.ascii_LAR[51] == ([0.0 4.0 2.0 4.0 4.0 3.0 1.0 0.0 0.0; 6.0 6.0 4.0 2.0 1.0 0.0 0.0 1.0 2.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9]])
	@test Plasm.ascii_LAR[52] == ([4.0 0.0 4.0 4.0; 1.0 1.0 6.0 0.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4]])
	@test Plasm.ascii_LAR[53] == ([4.0 0.0 0.0 3.0 4.0 4.0 3.0 1.0 0.0 0.0; 6.0 6.0 4.0 4.0 3.0 1.0 0.0 0.0 1.0 2.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10]])
	@test Plasm.ascii_LAR[54] == ([4.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 3.0 1.0 0.0; 6.0 6.0 5.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 3.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10], [10, 11]])
	@test Plasm.ascii_LAR[55] == ([0.0 0.0 4.0 0.0; 5.0 6.0 6.0 0.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4]])
	@test Plasm.ascii_LAR[56] == ([1.0 3.0 4.0 4.0 3.0 1.0 0.0 0.0 4.0 4.0 3.0 1.0 0.0 0.0; 0.0 0.0 1.0 2.0 3.0 3.0 2.0 1.0 4.0 5.0 6.0 6.0 5.0 4.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 1], [6, 5], [5, 9], [9, 10], [10, 11], [11, 12], [12, 13], [13, 14], [14, 6]])
	@test Plasm.ascii_LAR[57] == ([0.0 3.0 4.0 4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0; 0.0 0.0 1.0 5.0 6.0 6.0 5.0 3.0 2.0 2.0 3.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10], [10, 11]])
	@test Plasm.ascii_LAR[58] == ([2.0 2.0 1.5 1.5 2.0 2.0 1.5 1.5; 1.0 1.5 1.5 1.0 3.0 3.5 3.5 3.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 1], [5, 6], [6, 7], [7, 8], [8, 5]])
	@test Plasm.ascii_LAR[59] == ([2.0 2.0 1.5 1.5 1.0 2.0 2.0 2.0 1.5 1.5; 3.0 3.5 3.5 3.0 -0.5 0.5 1.0 1.5 1.5 1.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 1], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10], [10, 7]])
	@test Plasm.ascii_LAR[60] == ([3.0 0.0 3.0; 6.0 3.0 0.0], Array{Int64,1}[[1, 2], [2, 3]])
	@test Plasm.ascii_LAR[61] == ([1.0 3.0 1.0 3.0; 2.5 2.5 3.5 3.5], Array{Int64,1}[[1, 2], [3, 4]])
	@test Plasm.ascii_LAR[62] == ([1.0 4.0 1.0; 6.0 3.0 0.0], Array{Int64,1}[[1, 2], [2, 3]])
	@test Plasm.ascii_LAR[63] == ([2.0 2.0 1.5 1.5 1.75 1.75 3.0 3.0 2.0 1.0 0.0 0.0; 0.0 0.5 0.5 0.0 1.0 2.75 4.0 5.0 6.0 6.0 5.0 4.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 1], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10], [10, 11], [11, 12]])
	@test Plasm.ascii_LAR[64] == ([4.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 2.0 1.0 2.0 3.0 2.0; 0.0 0.0 1.0 3.0 4.0 4.0 3.0 1.0 1.0 2.0 3.0 2.0 1.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10], [10, 11], [11, 12], [12, 13]])
	@test Plasm.ascii_LAR[65] == ([0.0 0.0 1.0 3.0 4.0 4.0 0.0 4.0; 0.0 5.0 6.0 6.0 5.0 0.0 2.0 2.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [7, 8]])
	@test Plasm.ascii_LAR[66] == ([0.0 0.0 3.0 4.0 4.0 3.0 4.0 4.0 3.0 0.0 0.0 3.0; 0.0 6.0 6.0 5.0 4.0 3.0 2.0 1.0 0.0 0.0 3.0 3.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10], [11, 12]])
	@test Plasm.ascii_LAR[67] == ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0; 1.0 0.0 0.0 1.0 5.0 6.0 6.0 5.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8]])
	@test Plasm.ascii_LAR[68] == ([0.0 0.0 3.0 4.0 4.0 3.0 0.0; 0.0 6.0 6.0 5.0 1.0 0.0 0.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7]])
	@test Plasm.ascii_LAR[69] == ([4.0 0.0 0.0 4.0 0.0 3.0; 0.0 0.0 6.0 6.0 3.0 3.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [5, 6]])
	@test Plasm.ascii_LAR[70] == ([0.0 0.0 4.0 0.0 3.0; 0.0 6.0 6.0 3.0 3.0], Array{Int64,1}[[1, 2], [2, 3], [4, 5]])
	@test Plasm.ascii_LAR[71] == ([2.0 4.0 4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0; 3.0 3.0 1.0 0.0 0.0 1.0 5.0 6.0 6.0 5.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10]])
	@test Plasm.ascii_LAR[72] == ([0.0 0.0 4.0 4.0 0.0 4.0; 0.0 6.0 6.0 0.0 3.0 3.0], Array{Int64,1}[[1, 2], [3, 4], [5, 6]])
	@test Plasm.ascii_LAR[73] == ([2.0 2.0 1.0 3.0 1.0 3.0; 0.0 6.0 0.0 0.0 6.0 6.0], Array{Int64,1}[[1, 2], [3, 4], [5, 6]])
	@test Plasm.ascii_LAR[74] == ([0.0 1.0 2.0 3.0 3.0 2.0 4.0; 1.0 0.0 0.0 1.0 6.0 6.0 6.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [6, 7]])
	@test Plasm.ascii_LAR[75] == ([4.0 0.0 4.0 0.0 0.0; 6.0 3.0 0.0 0.0 6.0], Array{Int64,1}[[1, 2], [2, 3], [4, 5]])
	@test Plasm.ascii_LAR[76] == ([4.0 0.0 0.0; 0.0 0.0 6.0], Array{Int64,1}[[1, 2], [2, 3]])
	@test Plasm.ascii_LAR[77] == ([0.0 0.0 2.0 4.0 4.0; 0.0 6.0 4.0 6.0 0.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5]])
	@test Plasm.ascii_LAR[78] == ([0.0 0.0 4.0 4.0 4.0; 0.0 6.0 2.0 0.0 6.0], Array{Int64,1}[[1, 2], [2, 3], [4, 5]])
	@test Plasm.ascii_LAR[79] == ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0; 1.0 0.0 0.0 1.0 5.0 6.0 6.0 5.0 1.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9]])
	@test Plasm.ascii_LAR[80] == ([0.0 0.0 3.0 4.0 4.0 3.0 0.0; 0.0 6.0 6.0 5.0 3.0 2.0 2.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7]])
	@test Plasm.ascii_LAR[81] == ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 3.0 4.0; 1.0 0.0 0.0 1.0 5.0 6.0 6.0 5.0 1.0 1.0 0.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [10, 11]])
	@test Plasm.ascii_LAR[82] == ([0.0 0.0 3.0 4.0 4.0 3.0 0.0 3.0 4.0; 0.0 6.0 6.0 5.0 3.0 2.0 2.0 2.0 0.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [8, 9]])
	@test Plasm.ascii_LAR[83] == ([0.0 1.0 3.0 4.0 4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 4.0 5.0 6.0 6.0 5.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10], [10, 11], [11, 12]])
	@test Plasm.ascii_LAR[84] == ([2.0 2.0 0.0 4.0; 0.0 6.0 6.0 6.0], Array{Int64,1}[[1, 2], [3, 4]])
	@test Plasm.ascii_LAR[85] == ([0.0 0.0 1.0 3.0 4.0 4.0; 6.0 1.0 0.0 0.0 1.0 6.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6]])
	@test Plasm.ascii_LAR[86] == ([0.0 2.0 4.0; 6.0 0.0 6.0], Array{Int64,1}[[1, 2], [2, 3]])
	@test Plasm.ascii_LAR[87] == ([0.0 0.0 1.0 2.0 3.0 4.0 4.0; 6.0 3.0 0.0 3.0 0.0 3.0 6.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7]])
	@test Plasm.ascii_LAR[88] == ([0.0 4.0 0.0 4.0; 0.0 6.0 6.0 0.0], Array{Int64,1}[[1, 2], [3, 4]])
	@test Plasm.ascii_LAR[89] == ([0.0 2.0 4.0 2.0 2.0; 6.0 2.0 6.0 2.0 0.0], Array{Int64,1}[[1, 2], [2, 3], [4, 5]])
	@test Plasm.ascii_LAR[90] == ([0.0 4.0 0.0 4.0; 6.0 6.0 0.0 0.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4]])
	@test Plasm.ascii_LAR[91] == ([2.0 1.0 1.0 2.0; 0.0 0.0 6.0 6.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4]])
	@test Plasm.ascii_LAR[92] == ([1.0 3.0; 6.0 0.0], Array{Int64,1}[[1, 2]])
	@test Plasm.ascii_LAR[93] == ([2.0 3.0 3.0 2.0; 0.0 0.0 6.0 6.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4]])
	@test Plasm.ascii_LAR[94] == ([1.0 2.0 3.0; 5.0 6.0 5.0], Array{Int64,1}[[1, 2], [2, 3]])
	@test Plasm.ascii_LAR[95] == ([1.0 4.0; 0.0 0.0], Array{Int64,1}[[1, 2]])
	@test Plasm.ascii_LAR[96] == ([2.0 2.0 3.0 2.5 2.5 2.0; 4.5 5.0 6.0 4.0 4.5 4.0], Array{Int64,1}[[1, 2], [2, 3], [4, 5], [5, 1], [1, 6], [6, 4]])
	@test Plasm.ascii_LAR[97] == ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 4.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 0.0 3.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 1], [9, 10]])
	@test Plasm.ascii_LAR[98] == ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 0.0 0.0 1.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 0.0 5.0 5.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 1], [9, 10], [10, 11]])
	@test Plasm.ascii_LAR[99] == ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8]])
	@test Plasm.ascii_LAR[100] == ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 4.0 3.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 0.0 5.0 5.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 1], [9, 10], [10, 11]])
	@test Plasm.ascii_LAR[101] == ([4.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 0.0; 0.0 0.0 1.0 2.0 3.0 3.0 2.0 1.0 1.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9]])
	@test Plasm.ascii_LAR[102] == ([4.0 4.0 3.0 2.0 1.0 1.0 0.0 2.0; 3.0 4.0 5.0 5.0 4.0 0.0 1.0 1.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [7, 8]])
	@test Plasm.ascii_LAR[103] == ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 3.0 1.0 0.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 0.0 -1.0 -1.0 0.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 1], [1, 9], [9, 10], [10, 11], [11, 12]])
	@test Plasm.ascii_LAR[104] == ([4.0 4.0 3.0 1.0 0.0 0.0 0.0 1.0; 0.0 2.0 3.0 3.0 2.0 0.0 5.0 5.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [6, 7], [7, 8]])
	@test Plasm.ascii_LAR[105] == ([1.0 3.0 1.0 3.0 2.0 2.0 2.25 2.25 1.75 1.75; 0.0 0.0 3.0 3.0 0.0 3.0 3.75 4.25 4.25 3.75], Array{Int64,1}[[1, 2], [3, 4], [5, 6], [7, 8], [8, 9], [9, 10], [10, 7]])
	@test Plasm.ascii_LAR[106] == ([1.0 3.0 2.0 2.0 1.0 0.0 2.25 2.25 1.75 1.75; 3.0 3.0 3.0 0.0 -1.0 0.0 3.75 4.25 4.25 3.75], Array{Int64,1}[[1, 2], [3, 4], [4, 5], [5, 6], [7, 8], [8, 9], [9, 10], [10, 7]])
	@test Plasm.ascii_LAR[107] == ([0.0 1.0 1.0 0.0 4.0 2.0 1.0 3.0 4.0; 0.0 0.0 3.0 3.0 0.0 0.0 1.0 3.0 3.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [5, 6], [6, 7], [7, 8], [8, 9]])
	@test Plasm.ascii_LAR[108] == ([2.0 2.0 1.0 1.0 3.0; 0.0 5.0 5.0 0.0 0.0], Array{Int64,1}[[1, 2], [2, 3], [4, 5]])
	@test Plasm.ascii_LAR[109] == ([4.0 4.0 2.0 2.0 2.0 0.0 0.0 0.0; 0.0 3.0 2.0 0.0 3.0 2.0 0.0 3.0], Array{Int64,1}[[1, 2], [2, 3], [4, 5], [5, 6], [7, 8]])
	@test Plasm.ascii_LAR[110] == ([3.0 3.0 1.0 1.0 1.0; 0.0 3.0 2.0 0.0 3.0], Array{Int64,1}[[1, 2], [2, 3], [4, 5]])
	@test Plasm.ascii_LAR[111] == ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 1.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9]])
	@test Plasm.ascii_LAR[112] == ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 0.0 0.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 3.0 -1.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 1], [9, 10]])
	@test Plasm.ascii_LAR[113] == ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 4.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 3.0 -1.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 1], [9, 10]])
	@test Plasm.ascii_LAR[114] == ([0.0 2.0 1.0 1.0 1.0 2.0 3.0 4.0; 0.0 0.0 0.0 3.0 2.0 3.0 3.0 2.0], Array{Int64,1}[[1, 2], [3, 4], [5, 6], [6, 7], [7, 8]])
	@test Plasm.ascii_LAR[115] == ([0.0 4.0 3.0 1.0 0.0 1.0 3.0 4.0; 0.0 0.0 1.0 1.0 2.0 3.0 3.0 2.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8]])
	@test Plasm.ascii_LAR[116] == ([1.0 3.0 2.0 2.0 2.0 3.0; 0.0 0.0 0.0 5.0 4.0 4.0], Array{Int64,1}[[1, 2], [3, 4], [5, 6]])
	@test Plasm.ascii_LAR[117] == ([0.0 1.0 1.0 2.0 3.0 4.0 4.0; 3.0 3.0 1.0 0.0 0.0 1.0 3.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7]])
	@test Plasm.ascii_LAR[118] == ([0.0 1.0 3.0 4.0; 3.0 0.0 3.0 3.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4]])
	@test Plasm.ascii_LAR[119] == ([0.0 0.0 1.0 2.0 3.0 4.0 4.0; 3.0 2.0 0.0 2.0 0.0 2.0 3.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7]])
	@test Plasm.ascii_LAR[120] == ([0.0 1.0 4.0 1.0 4.0; 3.0 3.0 0.0 0.0 3.0], Array{Int64,1}[[1, 2], [2, 3], [4, 5]])
	@test Plasm.ascii_LAR[121] == ([0.0 1.0 2.5 0.0 1.0 4.0; 3.0 3.0 1.5 0.0 0.0 3.0], Array{Int64,1}[[1, 2], [2, 3], [4, 5], [5, 6]])
	@test Plasm.ascii_LAR[122] == ([0.0 0.0 3.0 0.0 3.0 4.0; 2.0 3.0 3.0 0.0 0.0 1.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6]])
	@test Plasm.ascii_LAR[123] == ([2.5 2.0 2.0 1.5 2.0 2.0 2.5; 6.5 6.0 3.5 3.0 2.5 0.0 -0.5], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7]])
	@test Plasm.ascii_LAR[124] == ([2.0 2.0; 0.0 5.0], Array{Int64,1}[[1, 2]])
	@test Plasm.ascii_LAR[125] == ([1.5 2.0 2.0 2.5 2.0 2.0 1.5; 6.5 6.0 3.5 3.0 2.5 0.0 -0.5], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7]])
	@test Plasm.ascii_LAR[126] == ([1.0 1.75 2.75 3.5; 5.0 5.5 5.0 5.5], Array{Int64,1}[[1, 2], [2, 3], [3, 4]])
end

@testset "ascii codes Tests 2" begin
	@test length(Plasm.hpcs)==95
	@test Plasm.hpcs[1]==Plasm.ascii32
	@test Plasm.hpcs[10]==Plasm.ascii41
	@test Plasm.hpcs[12]==Plasm.ascii43
	@test Plasm.hpcs[20]==Plasm.ascii51
	@test Plasm.hpcs[30]==Plasm.ascii61
	@test Plasm.hpcs[40]==Plasm.ascii71
	@test Plasm.hpcs[50]==Plasm.ascii81
	@test Plasm.hpcs[60]==Plasm.ascii91
	@test Plasm.hpcs[70]==Plasm.ascii101
	@test Plasm.hpcs[80]==Plasm.ascii111
	@test Plasm.hpcs[95]==Plasm.ascii126
end

@testset "constants Tests" begin
	@test Plasm.textalignment == "centre"
	@test Plasm.textangle == pi/4
	@test Plasm.textwidth == 1.0
	@test Plasm.textheight == 1.0
	@test Plasm.textspacing == 0.25
	@test Plasm.fontwidth == 4.0
	@test Plasm.fontheight == 8.0
	@test Plasm.fontspacing == 1.0
end

@testset "various Tests" begin
	charlist = "Alberto"
	@test Plasm.charpols(charlist)[2][2]==[[1, 2],[2, 3],[4, 5]]
	@test typeof(Plasm.charpols(charlist))==
	Array{Tuple{Array{Float64,2},Array{Array{Int64,1},1}},1}
	@test typeof(Plasm.charpols(charlist)[2])==
	Tuple{Array{Float64,2},Array{Array{Int64,1},1}}
	@test typeof(Plasm.charpols(charlist)[2][2])==
	Array{Array{Int64,1},1}
	@test Plasm.distr([[1,2,3],0])==[[1,0],[2,0],[3,0]]
	@test Plasm.cat([[1,2],[3,4]])==[1,2,3,4]
	@test Plasm.translate(-1)(LinearAlgebraicRepresentation.cuboid([1,1]))==
	([-1.0 -1.0 0.0 0.0; 0.0 1.0 0.0 1.0], Array{Int64,1}[[1, 2, 3, 4]])
end


text:
@safetestset "getpoints, getpointspec and getpointnum" begin
    using PoincareInvariants

    @test getpointspec(1, FirstFourierPlan) == 1
    @test getpointspec(50, FirstFourierPlan) == 50
    @test getpointspec(345, FirstFourierPlan) == 345

    @inferred Matrix{Float32} getpoints(x -> (x, 2x), Float32, 8, FirstFourierPlan)
    @inferred Vector{Float64} getpoints(x -> sin(x), Float64, 124, FirstFourierPlan)

    @test getpoints(identity, Float64, 4, FirstFourierPlan) ≈ [0, 0.25, 0.5, 0.75]
    @test getpoints(identity, Float64, 6, FirstFourierPlan) ≈ [n / 6 for n in 0:5]
    @test getpoints(cos, Float64, 124, FirstFourierPlan) ≈ [cos(n / 124) for n in 0:123]

    N = 6342
    testpnts = getpoints(identity, Float64, N, FirstFourierPlan)
    pnts3 = getpoints(Float32, N, FirstFourierPlan) do x
        x, x^2, sin(x)
    end
    @test pnts3 isa Matrix{Float32}
    @test size(pnts3) == (N, 3)

    @test pnts3[:, 1] ≈ testpnts
    @test pnts3[:, 2] ≈ testpnts .* testpnts
    @test pnts3[:, 3] ≈ sin.(testpnts)
end

@safetestset "compute! with FirstFourierPlan" begin
    using PoincareInvariants

    T = Float64
    D = 2

    # rotating vector field
    θ(z, p, t) = [-z[2], z[1]]

    # circular path with radius R
    # ḟ = 2π .* R .* (-sinpi(2ϕ), cospi(2ϕ))
    # ∫ ḟ⋅θ dϕ = 2π * R
    R = 0.75
    f(ϕ) = R .* (cospi(2ϕ), sinpi(2ϕ))

    I = 2π * R^2

    N = 10
    plan = FirstFourierPlan{T, D}(θ, N)
    pinv = FirstPoincareInvariant{T, D}(θ, N, plan)
    pnts = getpoints(f, Float64, N, FirstFourierPlan)
    @test I ≈ compute!(pinv, pnts, 0.0, nothing) atol=10eps()
end


text:
#=
This is a very crude first stab at the Tables.jl interface
https://github.com/JuliaData/Tables.jl
=#

using Tables

Tables.istable(::Type{<:KeyedArray}) = true

Tables.rowaccess(::Type{<:KeyedArray}) = true

function Tables.rows(A::Union{KeyedArray, NdaKa})
    L = hasnames(A) ? (dimnames(A)..., :value) :  # should gensym() if :value in dimnames(A)
        (ntuple(d -> Symbol(:dim_,d), ndims(A))..., :value)
    R = keys_or_axes(A)
    nt(inds) = NamedTuple{L}((map(getindex, R, inds)..., A[inds...]))
    # (nt(inds) for inds in Iterators.product(axes(A)...)) # should flatten?
    (nt(inds) for inds in Vectorator(Iterators.product(axes(A)...)))
end

#=
rr = wrapdims(rand(2,3), 11:12, 21:23)
nn = wrapdims(rand(2,3), a=11:12, b=21:23)

Tables.rows(rr) |> collect |> vec
Tables.rows(nn) |> collect |> vec

Tables.Schema(nn) # define a struct? Now below...

# No error if Tables.rows's generator has size,
# it uses Vectorator mostly to give Tables.Schema something to find.
=#

Tables.columnaccess(::Type{<:KeyedArray{T,N,AT}}) where {T,N,AT} =
    IndexStyle(AT) === IndexLinear()

function Tables.columns(A::Union{KeyedArray, NdaKa})
    L = hasnames(A) ? (dimnames(A)..., :value) :
        (ntuple(d -> Symbol(:dim_,d), ndims(A))..., :value)
    R = keys_or_axes(A)
    G = ntuple(ndims(A)) do d
        vec([rs[d] for rs in Iterators.product(R...)])
        # _vec(rs[d] for rs in Iterators.product(R...))
    end
    C = (G..., vec(parent(A)))
    NamedTuple{L}(C)
end

function Tables.Schema(nt::NamedTuple) # 🏴☠
    L = keys(nt)
    T = map(v -> typeof(first(v)), values(nt))
    Tables.Schema(L,T)
end

#=
Ah, iterators aren't allowed for columns, must be indexable:
https://github.com/JuliaData/Tables.jl/issues/101
They could be something like this, but seems overkill?
https://github.com/MichielStock/Kronecker.jl
https://github.com/JuliaArrays/LazyArrays.jl#kronecker-products

Tables.columns(nn)
map(collect, Tables.columns(nn))

using DataFrames

DataFrame(rand(2,3))
DataFrame(nn) # doesn't see Tables

dd1 = DataFrame(Tables.rows(nn))
dd2 = DataFrame(Tables.columns(nn))
=#

"""
    Vectorator(iter)

Wrapper for iterators which ensures they do not have an n-dimensional size.
Tries to ensure that `collect(Vectorator(iter)) == vec(collect(iter))`.
"""
struct Vectorator{T}
    iter::T
end

_vec(iter) = (x for x in Vectorator(iter))

Base.iterate(x::Vectorator, s...) = iterate(x.iter, s...)

Base.length(x::Vectorator) = length(x.iter)

Base.IteratorSize(::Type{Vectorator{T}}) where {T} =
    Base.IteratorSize(T) isa Base.HasShape ? Base.HasLength() : IteratorSize(T)

Base.IteratorEltype(::Type{Vectorator{T}}) where {T} = Base.IteratorEltype(T)

Base.eltype(::Type{Vectorator{T}}) where {T} = eltype(T)

function Tables.Schema(rows::Base.Generator{<:Vectorator})
    row = first(rows)
    Tables.Schema(keys(row), map(typeof, values(row)))
end

# struct OneKron{T, AT} <: AbstractVector{T}
#     data::AT
#     inner::Int
#     outer::Int
# end

# Tables.materializer(A::KeyedArray) = wrapdims

# function wrapdims(tab)
#     sch = Tables.Schema(tab)
#     for r in Tables.rows(tab)

#     end
# end

"""
    AxisKeys.populate!(A, table, value; force=false)

Populate `A` with the contents of the `value` column in a provided `table`, matching the
[Tables.jl](https://github.com/JuliaData/Tables.jl) API. The `table` must contain columns
corresponding to the keys in `A` and implements `Tables.rows`. If the keys in `A` do not
uniquely identify rows in the `table` then an `ArgumentError` is throw. If `force` is true
then the duplicate (non-unique) entries will be overwritten.
"""
function populate!(A, table, value::Symbol; force=false)
    # Use a BitArray mask to detect duplicates and error instead of overwriting.
    mask = force ? falses() : falses(size(A))

    for r in Tables.rows(table)
        vals = Tuple(Tables.getcolumn(r, c) for c in dimnames(A))
        inds = map(findindex, vals, axiskeys(A))

        # Handle duplicate error checking if applicable
        if !force
            # Error if mask already set.
            mask[inds...] && throw(ArgumentError("Key $vals is not unique"))
            # Set mask, marking that we've set this index
            setindex!(mask, true, inds...)
        end

        # Insert our value into the data array
        setindex!(A, Tables.getcolumn(r, value), inds...)
    end

    return A
end

"""
    wrapdims(table, value, names...; default=undef, sort=false, force=false)

Construct `KeyedArray(NamedDimsArray(A,names),keys)` from a `table` matching
the [Tables.jl](https://github.com/JuliaData/Tables.jl) API.
(It must support both `Tables.columns` and `Tables.rows`.)

The contents of the array is taken from the column `value::Symbol` of the table.
Each symbol in `names` specifies a column whose unique entries
become the keys along a dimenension of the array.

If there is no row in the table matching a possible set of keys,
then this element of the array is undefined, unless you provide the `default` keyword.
If several rows share the same set of keys, then by default an `ArgumentError` is thrown.
Keyword `force=true` will instead cause these non-unique entries to be overwritten.

Setting `AxisKeys.nameouter() = false` will reverse the order of wrappers produced.
"""
function wrapdims(table, value::Symbol, names::Symbol...; kw...)
    if nameouter() == false
        _wrap_table(KeyedArray, identity, table, value, names...; kw...)
    else
        _wrap_table(NamedDimsArray, identity, table, value, names...; kw...)
    end
end

"""
    wrapdims(df, UniqueVector, :val, :x, :y)

Converts at Tables.jl table to a `KeyedArray` + `NamedDimsArray` pair,
using column `:val` for values, and columns `:x, :y` for names & keys.
Optional 2nd argument applies this type to all the key-vectors.
"""
function wrapdims(table, KT::Type, value::Symbol, names::Symbol...; kw...)
    if nameouter() == false
        _wrap_table(KeyedArray, KT, table, value, names...; kw...)
    else
        _wrap_table(NamedDimsArray, KT, table, value, names...; kw...)
    end
end

function _wrap_table(AT::Type, KT, table, value::Symbol, names::Symbol...; default=undef, sort::Bool=false, kwargs...)
    # get columns of the input table source
    cols = Tables.columns(table)

    # Extract key columns
    pairs = map(names) do k
        col = unique(Tables.getcolumn(cols, k))
        sort && Base.sort!(col)
        return k => KT(col)
    end

    # Extract data/value column
    vals = Tables.getcolumn(cols, value)

    # Initialize the KeyedArray
    sz = length.(last.(pairs))
    if default === undef
        data = similar(vals, sz)
    else
        data = similar(vals, Union{eltype(vals), typeof(default)}, sz)
        fill!(data, default)
    end
    A = AT(data; pairs...)

    populate!(A, table, value; kwargs...)
    return A
end


text:
abstract type AbstractTrainer end

mutable struct Trainer <: AbstractTrainer
    model
    optimizer
    loss_list
    eval_interval
    current_epoch
    function Trainer(model, optimizer)
        self = new()
        self.model = model
        self.optimizer = optimizer
        self.loss_list = []
        self.eval_interval = nothing
        self.current_epoch = 0
        return self
    end
end

function fit!(trainer::Trainer, x, t; max_epoch = 10, batch_size = 32, max_grad = nothing, eval_interval = 20)
    data_size = size(x, 1)
    max_iters = data_size ÷ batch_size
    trainer.eval_interval = eval_interval

    @time train!(trainer, data_size, batch_size, max_epoch, max_iters, x, t)
end

function train!(trainer::Trainer, data_size, batch_size, max_epoch, max_iters, train_x, train_t)
    total_loss = 0
    loss_count = 0

    for epoch = 1:max_epoch
        # シャッフルしたミニバッチの生成
        indices = rand(1:data_size, batch_size) # 1～data_sizeの範囲の一意な数値を要素に持つ要素数data_sizeのVectorを返す
        batch_x = selectdim(train_x, 1, indices)
        batch_t = selectdim(train_t, 1, indices)

        for iter = 1:max_iters
            total_loss, loss_count = iteration!(trainer, iter, max_iters, total_loss, loss_count, batch_x, batch_t)
        end

        trainer.current_epoch += 1
    end
end

function iteration!(trainer::Trainer, iter, max_iters, total_loss, loss_count, batch_x, batch_t)
    # 損失を計算
    loss = forward!(trainer.model, batch_x, batch_t) # BaseModel(TwoLayersNet)内に定義
    # 勾配を求める
    backward!(trainer.model) # BaseModel(TwoLayersNet)内に定義
    params, grads = remove_duplicate(trainer.model.params, trainer.model.grads)
    # パラメータを更新
    update!(trainer.optimizer, trainer.model.params, trainer.model.grads) # Optimizer内に定義
    # 損失履歴を追加
    push!(trainer.loss_list, loss)
    total_loss += loss
    loss_count += 1

    # 定期的に学習経過を出力（avg_lossはloss_listに追加しない）
    if (iter * trainer.current_epoch % trainer.eval_interval) == 0
        avg_loss = total_loss / loss_count
        println("| epoch $(trainer.current_epoch) | iter $(iter) / $(max_iters) | loss $(round(avg_loss, digits=2))")
        total_loss, loss_count = 0, 0
    end
    return total_loss, loss_count
end

function remove_duplicate(params, grads)
    #= パラメータ配列の重複する重みを一つに集約し，
    その重みに対応する勾配を加算する =#
    _params, _grads = params, grads

    # 重複を検出するまでループ
    for h = 1:lastindex(_params)
        for i = 1:(lastindex(_params[h])-1)
            for j = (i+1):lastindex(_params[h])
                # 重みを共有する場合
                if _params[h][i, :] == _params[h][j, :]
                    _grads[h][i, :] .+= _grads[h][j, :]
                    pop!(_params[h][j, :])
                    pop!(_grads[h][j, :])
                    return _params, _grads
                    # 転置行列として重みを共有する場合
                elseif (ndims(_params[h][i, :]) == 2) && (ndims(_params[h][j, :]) == 2) && (size(_params[h][i, :]') == size(_params[h][j, :])) && (_params[h][i, :]' == _params[h][j, :])
                    _grads[h][i, :] .+= _grads[h][j, :]'
                    pop!(_params[h][j, :])
                    pop!(_grads[h][j, :])
                    return _params, _grads
                end
            end
        end
    end
    return _params, _grads
end


text:
"""
    T = eigmat(M,x;ω=1.0,basis="Hermite",α=1.0)

Create a matrix of orthonormal mode functions for the chosen basis.
The matrix is in a form that allows the transformation
from the mode coefficients to the spatial grid `x`. In 1D:

```math
\\psi(x) = \\sum_j c_j\\phi_j(x) \\equiv Tc
```

where `c` is a column vector of coefficients in the basis.

`M`: number of modes in the basis.

`x`: spatial grid to which the coefficients are mapped.

`ω`: mode angular frequency, in units of the chosen reference frequency.

`α`: extra input for the `"Laguerre"` basis.

`basis`: set of eigenfunctions representing the c-field. At present `basis = "Hermite"` is implemented.

"""

function eigmat(M,x;ω=1.0,basis="Hermite",α=0.0)
        @assert M>=1
  if basis=="Hermite"
    M > 371 && error("Float64 Quadrature does not converge for M > 371.")
    ψ0 = exp.(-(√ω*x).^2/2)*(ω/π)^(1/4)
    ψ1 = sqrt(2)*exp.(-(√ω*x).^2/2).*(√ω*x)*(ω/π)^(1/4)
    n = collect(0:M-1)
    T = zeros(x*ones(1,M))
    T[:,1] = ψ0
    M >= 2 && (T[:,2] = ψ1)
    for m=1:M-2
      T[:,m+2]=sqrt(2/(n[m+2]))*(√ω*x).*T[:,m+1]-sqrt(n[m+1]/n[m+2])T[:,m]
    end
  elseif basis=="Laguerre"
    error(basis," ","basis not implemented")
    #M>26 && error("recursion unstable for M >=27")
    #alternate approach using ApproxFun.jl -same instability as
    #previous approach in logspace
#==
    T=zeros(x*ones(1,M))
    n=0:BigFloat(M)-1;
    Lnorm = exp.(lgamma.(n+α+1)-lgamma.(n+1))/ω

    for j = 0:M-1
      c=zeros(M)
      c[j+1]=1.
      Lna=Fun(Laguerre(α),c./sqrt(Lnorm));
      T[:,j+1]=Lna.(ω*x).*exp.(-ω*x/2).*(ω*x).^(α/2)
    end
==#
  else
    error(basis," ","basis not implemented")
  end
  return T
end