text:
# Retrieve number of variables from equation instance
@inline nvariables(::AbstractEquations{NDIMS, NVARS}) where {NDIMS, NVARS} = NVARS

# TODO: Taal performance, 1:NVARS vs. Base.OneTo(NVARS) vs. SOneTo(NVARS)
@inline eachvariable(equations::AbstractEquations) = Base.OneTo(nvariables(equations))

"""
    get_name(equations::AbstractEquations)

Returns the canonical, human-readable name for the given system of equations.

# Examples
```jldoctest
julia> Trixi.get_name(CompressibleEulerEquations1D(1.4))
"CompressibleEulerEquations1D"
```
"""
get_name(equations::AbstractEquations) = equations |> typeof |> nameof |> string

# Add methods to show some information on systems of equations.
function Base.show(io::IO, equations::AbstractEquations)
  # Since this is not performance-critical, we can use `@nospecialize` to reduce latency.
  @nospecialize equations # reduce precompilation time

  print(io, get_name(equations), " with ")
  if nvariables(equations) == 1
    print(io, "one variable")
  else
    print(io, nvariables(equations), " variables")
  end
end

function Base.show(io::IO, ::MIME"text/plain", equations::AbstractEquations)
  # Since this is not performance-critical, we can use `@nospecialize` to reduce latency.
  @nospecialize equations # reduce precompilation time

  if get(io, :compact, false)
    show(io, equations)
  else
    summary_header(io, get_name(equations))
    summary_line(io, "#variables", nvariables(equations))
    for variable in eachvariable(equations)
      summary_line(increment_indent(io),
                   "variable " * string(variable),
                   varnames(cons2cons, equations)[variable])
    end
    summary_footer(io)
  end
end

@inline Base.ndims(::AbstractEquations{NDIMS}) where NDIMS = NDIMS

# equations act like scalars in broadcasting
Base.broadcastable(equations::AbstractEquations) = Ref(equations)

"""
    flux(u, orientation_or_normal, equations)

Given the conservative variables `u`, calculate the (physical) flux in Cartesian
direction `orientation::Integer` or in arbitrary direction `normal::AbstractVector`
for the corresponding set of governing `equations`.
`orientation` is `1`, `2`, and `3` for the x-, y-, and z-directions, respectively.
"""
function flux end

"""
    rotate_to_x(u, normal, equations)

Apply the rotation that maps `normal` onto the x-axis to the convservative variables `u`.
This is used by [`FluxRotated`](@ref) to calculate the numerical flux of rotationally
invariant equations in arbitrary normal directions.

See also: [`rotate_from_x`](@ref)
"""
function rotate_to_x end

"""
    rotate_from_x(u, normal, equations)

Apply the rotation that maps the x-axis onto `normal` to the convservative variables `u`.
This is used by [`FluxRotated`](@ref) to calculate the numerical flux of rotationally
invariant equations in arbitrary normal directions.

See also: [`rotate_to_x`](@ref)
"""
function rotate_from_x end

"""
    BoundaryConditionDirichlet(boundary_value_function)

Create a Dirichlet boundary condition that uses the function `boundary_value_function`
to specify the values at the boundary.
This can be used to create a boundary condition that specifies exact boundary values
by passing the exact solution of the equation.
The passed boundary value function will be called with the same arguments as an initial condition function is called, i.e., as
```julia
boundary_value_function(x, t, equations)
```
where `x` specifies the coordinates, `t` is the current time, and `equation` is the corresponding system of equations.

# Examples
```julia
julia> BoundaryConditionDirichlet(initial_condition_convergence_test)
```
"""
struct BoundaryConditionDirichlet{B}
  boundary_value_function::B
end

# Dirichlet-type boundary condition for use with TreeMesh or CurvedMesh
@inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner, orientation_or_normal,
                                                                  direction,
                                                                  x, t,
                                                                  surface_flux_function, equations)
  u_boundary = boundary_condition.boundary_value_function(x, t, equations)

  # Calculate boundary flux
  if direction in (2, 4, 6) # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation_or_normal, equations)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation_or_normal, equations)
  end

  return flux
end

# Dirichlet-type boundary condition for use with UnstructuredQuadMesh
# Note: For unstructured we lose the concept of an "absolute direction"
@inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner,
                                                                  normal_direction::AbstractVector,
                                                                  x, t,
                                                                  surface_flux_function, equations)
  # get the external value of the solution
  u_boundary = boundary_condition.boundary_value_function(x, t, equations)

  # Calculate boundary flux
  flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations)

  return flux
end

"""
    BoundaryConditionWall(boundary_value_function)

Create a generic wall type boundary condition that uses the function `boundary_value_function`
to specify the external solution values.
The boundary wall function is called with arguments for an internal solution state from inside an
element `u_inner`, an outward pointing `normal_direction` and a particular set of `equations`, e.g.,
```julia
boundary_value_function(u_inner, normal_direction, equations)
```
which will return an external solution state.

# Example
```julia
julia> BoundaryConditionWall(boundary_state_slip_wall)
```

!!! warning "Experimental code"
    This boundary condition can change any time and is currently only implemented for the
    [`CompressibleEulerEquations2D`](@ref) and [`AcousticPerturbationEquations2D`](@ref).
"""
struct BoundaryConditionWall{B}
  boundary_value_function::B
end

@inline function (boundary_condition::BoundaryConditionWall)(u_inner,
                                                             normal_direction::AbstractVector,
                                                             x, t,
                                                             surface_flux_function, equations)
  # get the external value of the solution
  u_boundary = boundary_condition.boundary_value_function(u_inner, normal_direction, equations)

  flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations)

  return flux
end

# set sensible default values that may be overwritten by specific equations
have_nonconservative_terms(::AbstractEquations) = Val(false)
have_constant_speed(::AbstractEquations) = Val(false)

default_analysis_errors(::AbstractEquations)     = (:l2_error, :linf_error)
default_analysis_integrals(::AbstractEquations)  = (entropy_timederivative,)

"""
    cons2cons(u, equations)

Return the conserved variables `u`. While this function is as trivial as `identity`,
it is also as useful.
"""
@inline cons2cons(u, ::AbstractEquations) = u
function cons2prim#=(u, ::AbstractEquations)=# end
@inline Base.first(u, ::AbstractEquations) = first(u)

"""
    cons2prim(u, equations)

Convert the conserved variables `u` to the primitive variables for a given set of
`equations`. The inverse conversion is performed by [`prim2cons`](@ref).
"""
function cons2prim end

"""
    prim2cons(u, equations)

Convert the conserved variables `u` to the primitive variables for a given set of
`equations`. The inverse conversion is performed by [`cons2prim`](@ref).
"""
function prim2cons end

"""
    entropy(u, equations)

Return the chosen entropy of the conserved variables `u` for a given set of
`equations`.
"""
function entropy end

"""
    cons2entropy(u, equations)

Convert the conserved variables `u` to the entropy variables for a given set of
`equations` with chosen standard [`entropy`](@ref). The inverse conversion is
performed by [`entropy2cons`](@ref).
"""
function cons2entropy end

"""
    entropy2cons(w, equations)

Convert the entropy variables `w` based on a standard [`entropy`](@ref) to the
conserved variables for a given set of `equations` . The inverse conversion is
performed by [`cons2entropy`](@ref).
"""
function entropy2cons end

# FIXME: Deprecations introduced in v0.3
@deprecate varnames_cons(equations) varnames(cons2cons, equations)
@deprecate varnames_prim(equations) varnames(cons2prim, equations)
@deprecate flux_upwind(u_ll, u_rr, orientation, equations) flux_godunov(u_ll, u_rr, orientation, equations)
@deprecate calcflux(u, orientation, equations) flux(u, orientation, equations)

####################################################################################################
# Include files with actual implementations for different systems of equations.

# Numerical flux formulations that are independent of the specific system of equations
include("numerical_fluxes.jl")

# Linear scalar advection
abstract type AbstractLinearScalarAdvectionEquation{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
include("linear_scalar_advection_1d.jl")
include("linear_scalar_advection_2d.jl")
include("linear_scalar_advection_3d.jl")

# Inviscid Burgers
abstract type AbstractInviscidBurgersEquation{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
include("inviscid_burgers_1d.jl")

# CompressibleEulerEquations
abstract type AbstractCompressibleEulerEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
include("compressible_euler_1d.jl")
include("compressible_euler_2d.jl")
include("compressible_euler_3d.jl")

# CompressibleEulerMulticomponentEquations
abstract type AbstractCompressibleEulerMulticomponentEquations{NDIMS, NVARS, NCOMP} <: AbstractEquations{NDIMS, NVARS} end
include("compressible_euler_multicomponent_1d.jl")
include("compressible_euler_multicomponent_2d.jl")

# Retrieve number of components from equation instance for the multicomponent case
@inline ncomponents(::AbstractCompressibleEulerMulticomponentEquations{NDIMS, NVARS, NCOMP}) where {NDIMS, NVARS, NCOMP} = NCOMP
@inline eachcomponent(equations::AbstractCompressibleEulerMulticomponentEquations) = Base.OneTo(ncomponents(equations))

# Ideal MHD
abstract type AbstractIdealGlmMhdEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
include("ideal_glm_mhd_1d.jl")
include("ideal_glm_mhd_2d.jl")
include("ideal_glm_mhd_3d.jl")

# IdealGlmMhdMulticomponentEquations
abstract type AbstractIdealGlmMhdMulticomponentEquations{NDIMS, NVARS, NCOMP} <: AbstractEquations{NDIMS, NVARS} end
include("ideal_glm_mhd_multicomponent_1d.jl")
include("ideal_glm_mhd_multicomponent_2d.jl")

# Retrieve number of components from equation instance for the multicomponent case
@inline ncomponents(::AbstractIdealGlmMhdMulticomponentEquations{NDIMS, NVARS, NCOMP}) where {NDIMS, NVARS, NCOMP} = NCOMP
@inline eachcomponent(equations::AbstractIdealGlmMhdMulticomponentEquations) = Base.OneTo(ncomponents(equations))

# Diffusion equation: first order hyperbolic system
abstract type AbstractHyperbolicDiffusionEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
include("hyperbolic_diffusion_1d.jl")
include("hyperbolic_diffusion_2d.jl")
include("hyperbolic_diffusion_3d.jl")

# Lattice-Boltzmann equation (advection part only)
abstract type AbstractLatticeBoltzmannEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
include("lattice_boltzmann_2d.jl")
include("lattice_boltzmann_3d.jl")

# Acoustic perturbation equations
abstract type AbstractAcousticPerturbationEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
include("acoustic_perturbation_2d.jl")


text:
"""
    Robust()
Robust Scaler. This scaler scale data as:
``\\tilde{\\boldsymbol{x}} = \\frac{\\boldsymbol{x}-Q2}{Q3 - Q1}``

# Example
```jldoctest preprocessing
julia> scaler = Robust()
Robust(Float64[])

julia> fit!(scaler, x)
3×2 Matrix{Float64}:
 0.412913  0.739911
 0.602654  0.873014
 0.849116  0.905986

julia> transform!(scaler, x)
20×2 Matrix{Float64}:
 -0.0121192   0.041181
 -0.0260262  -0.28201
  0.649595    0.417263
  0.89357     0.764633
 -0.637774   -0.978423
  0.630675    0.0812893
  0.910919    0.3256
 -0.923097   -1.87636
 -0.0230992  -0.249283
  0.16723     0.0942497
 -0.261766   -0.742477
  0.258819   -0.041181
  0.547691    0.213832
 -0.367388   -0.392802
  0.616988    0.193438
 -0.970618   -2.10092
  0.283712    0.430314
  0.0121192   0.135449
 -1.02922    -2.64305
 -1.38159    -5.25675
```
"""
mutable struct Robust
    p::AbstractVecOrMat
    Robust() = new(Array{Float64}(undef, 0))
end

function fit!(scaler::Robust, x; dims=1)
    if dims == 1
        scaler.p = hcat([quantile(x[:, i], [0.25, 0.5, 0.75]) for i in 1 : size(x, 2)]...)
    elseif dims == 2
        scaler.p = hcat([quantile(x[i, :], [0.25, 0.5, 0.75]) for i in 1 : size(x, 1)]...)
    end
end

rs(x, q1, q2, q3) = @. (x-q2) / (q3-q1)

function transform!(scaler::Robust, x; dims=1)
    p = scaler.p
    check_size(x, p, dims)
    if dims == 1
        for i in 1 : size(x, 2)
            x[:, i] = rs(x[:, i], p[:, i]...)
        end
    elseif dims == 2
        for i in 1 : size(x, 1)
            x[i, :] = rs(x[i, :], p[:, i]...)
        end
    end
    return x
end

function fit_transform!(scaler::Robust, x; dims=1)
    fit!(scaler, x, dims=dims)
    transform!(scaler, x; dims=dims)
end

irs(x, q1, q2, q3) = @. x*(q3-q1)+q2

function inv_transform!(scaler::Robust, x; dims=1)
    p = scaler.p
    check_size(x, p, dims)
    if dims == 1
        for i in 1 : size(x, 2)
            x[:, i] = irs(x[:, i], p[:, i]...)
        end
    elseif dims == 2
        for i in 1 : size(x, 1)
            x[i, :] = irs(x[i, :], p[:, i]...)
        end
    end
    return x
end

ext: jl


text:
export Comm
export barrier, broadcastAll, gatherAll, sumAll, maxAll, minAll, scanSum
export myPid, numProc, createDistributor
export getComm

# methods (and docs) are currently based straight off Epetra_Comm
# tpetra's equivalent seemed to be a wrapper to other Trilinos packages

# following julia's convention, n processors are labled 1 through n
# count variables are removed, since that information is contained in the arrays

"""
The base type for types that represent communication in parallel computing.

All subtypes must have the following methods, with CommImpl standing in for the subtype:

    barrier(comm::CommImpl)
Each processor must wait until all processors have arrived

    broadcastAll(comm::CommImpl, myvals::AbstractArray{T}, Root::Integer)::Array{T} where T
Takes a list of input values from the root processor and sends to all
other processors.  The values are returned (including on the root process)

    gatherAll(comm::CommImpl, myVals::AbstractArray{T})::Array{T} where T
Takes a list of input values from all processors and returns an ordered
contiguous list of those values on each processor

    sumAll(comm::CommImpl, partialsums::AbstractArray{T})::Array{T} where T
Take a list of input values from all processors and returns the sum on each
processor.  The method +(::T, ::T)::T can be assumed to exist

    maxAll(comm::CommImpl, partialmaxes::AbstractArray{T})::Array{T} where T
Takes a list of input values from all processors and returns the max to all
processors.  The method <(::T, ::T)::Bool can be assumed to exist

    minAll(comm::CommImpl, partialmins::AbstractArray{T})::Array{T} where T
Takes a list of input values from all processors and returns the min to all
processors.  The method <(::T, ::T)::Bool can be assumed to exist

    scanSum(comm::CommImpl, myvals::AbstractArray{T})::Array{T} where T
Takes a list of input values from all processors, computes the scan sum and
returns it to all processors such that processor i contains the sum of
values from processor 1 up to and including processor i.  The method
`+(::T, ::T)::T` is used for the addition.

    myPid(comm::CommImpl{GID, PID, LID})::PID
Returns the process rank

    numProc(comm::CommImpl{GID, PID, LID})::PID
Returns the total number of processes

    createDistributor(comm::CommImpl{GID, PID, LID})::Distributor{GID, PID, LID}
Create a distributor object

"""
abstract type Comm{GID <: Integer, PID <:Integer, LID <: Integer}
end

function Base.show(io::IO, comm::Comm)
    print(io, String(Symbol(typeof(comm))), " with PID ", myPid(comm),
                " and ", numProc(comm), " processes")
end

"""
    broadcastAll(::Comm, myVals::T, root::Integer)::T
    broadcastAll(comm::Comm, myVals::AbstractArray{T, 1}, root::Integer)::Array{T, 1}

Takes a list of input values from the root processor and sends it to each
other processor.  The broadcasted values are then returned, including on
the root process.
"""
broadcastAll(comm::Comm, myVal, root) = broadcastAll(comm, [myVal], root)[1]

"""
    gatherAll(::Comm, myVal::T)::Array{T, 1}
    gatherAll(comm::Comm, myVals::AbstractArray{T, 1})::Array{T, 1}

Takes a list of input values from all processors and returns an ordered,
contiguous list of those values.
"""
gatherAll(comm::Comm, myVal) = gatherAll(comm, [myVal])

"""
    sumAll(::Comm, partialsum::T)::T
    sumAll(comm::Comm, partialsums::AbstractArray{T, 1})::Array{T, 1}

Takes a list of input values from all processors and returns the sum on each
processor.  The method `+(::T, ::T)::T` must exist.
"""
sumAll(comm::Comm, val) = sumAll(comm, [val])[1]

"""
    maxAll(::Comm, partialmax::T)::T
    maxAll(comm::Comm, partialmaxes::AbstractArray{T, 1})::Array{T, 1}

Takes a list of input values from all processors and returns the max to all
processors.  The method `<(::T, ::T)::Bool` must exist.
"""
maxAll(comm::Comm, val) = maxAll(comm, [val])[1]

"""
    minAll(::Comm, partialmin::T)::T
    minAll(comm::Comm, partialmins::AbstractArray{T, 1})::Array{T, 1}

Takes a list of input values from all processors and returns the min to all
processors.  The method `<(::T, ::T)::Bool` must exist.
"""
minAll(comm::Comm, val) = minAll(comm, [val])[1]

"""
    scanSum(::Comm, myval::T)::T
    scanSum(comm::Comm, myvals::AbstractArray{T, 1})::Array{T, 1}

Takes a list of input values from all processors, computes the scan sum and
returns it to all processors such that processor `i` contains the sum of
values from processor 1 up to, and including, processor `i`.  The method
+(::T, ::T)::T must exist
"""
scanSum(comm::Comm, val) = scanSum(comm, [val])[1]

"""
    getComm(obj)

Gets the Comm for the object, if applicable
"""
getComm(comm::Comm) = comm

#### documentation for required methods ####

"""
    barrier(::Comm)

Causes the process to pause until all processes have called barrier.  Used to synchronize the processes
"""
function barrier end

"""
    myPid(::Comm{GID, PID, LID})::PID

Returns the rank of the calling processor
"""
function myPid end

"""
    numProc(::Comm{GID, PID, LID})::PID

Returns the total number of processes
"""
function numProc end

"""
    createDistributor(comm::Comm{GID, PID, LID})::Distributor{GID, PID, LID}

Creates a distributor for the given Comm object
"""
function createDistributor end


text:
# Convolute two vectors of probabilities

# Input:
# p1:       Vector of probabilities P(X = 0), ..., P(X = n)
# p2:       Vector of probabilities P(Y = 0), ..., P(Y = n)

# Output:
# out:      Vector of probabilities P(X + Y = 0), ..., P(X + Y = n)

function convolution(p1::Vector{Float64},
                     p2::Vector{Float64})::Vector{Float64}
    #
    M = length(p1)
    out = zeros(M)
    @inbounds for i = 0:M-1
        for j = 0:M-1-i
            out[i+j+1] += p1[i+1]*p2[j+1]
        end
    end
    out
end

function convolution(p1::Vector{Float64},
                     p2::Vector{Float64},
                     x::Int64)::Vector{Float64}
    #
    out = zeros(x+1)
    @inbounds for i = 0:x
        for j = 0:x-i
            out[i+j+1] += p1[i+1]*p2[j+1]
        end
    end
    out
end

function convolutionRev(p1::Vector{Float64},
                     p2::Vector{Float64},
                     x::Int64)::Vector{Float64}
    #
    out = zeros(x+1)
    @inbounds for i = 0:x
        for j = 0:x-i
            out[x + 1 - i - j] += p1[i+1]*p2[j+1]
        end
    end
    out
end

function convolutionRev(p1::Vector{Float64},
                     p2::Matrix{Float64},
                     x::Int64)::Matrix{Float64}
    #
    ncol = size(p2)[2]
    out = zeros(x+1, ncol)
    @inbounds for i = 0:x
        for j = 0:x-i
            for col = 1:ncol
                out[x + 1 - i - j, col] += p1[i+1]*p2[j+1, col]
            end
        end
    end
    out
end

function convolutionRevScale(p1::Vector{Float64},
                             p2::Matrix{Float64},
                             x::Int64,
                             scale::Vector{Float64})::Matrix{Float64}
    ncol = size(p2)[2]
    out = zeros(ncol, x+1)
    @inbounds for i = 0:x
        for j = 0:x-i
            scalefactor = scale[x + 1 - i - j]
            p1Part = p1[i+1]
            for col = 1:ncol
                out[col, x + 1 - i - j] += p1Part*p2[j+1, col]*scalefactor
            end
        end
    end
    out
end


text:
using DSGE
using Base.Test

# Can we actually test? Require that the FRED_API_KEY ENV is populated.
if haskey(ENV, "FRED_API_KEY")
    m = Model990()
    m.testing=true

    # Specify exact vintage and dates
    m <= Setting(:data_vintage, "151127")
    m <= Setting(:date_presample_start, quartertodate("2015-Q2"))
    m <= Setting(:date_mainsample_end, quartertodate("2015-Q3"))

    df = load_data(m; try_disk=false, verbose=:none)

    # Silly comparisons
    delete!(df, :date)
    colmeans = colwise(mean, df)
    @test_approx_eq_eps colmeans[1][1] 0.4435 1e-4
    @test_approx_eq_eps colmeans[2][1] -51.2707 1e-4
    @test_approx_eq_eps colmeans[3][1] 0.1582 1e-4
    @test_approx_eq_eps colmeans[4][1] 0.4271 1e-4
    @test_approx_eq_eps colmeans[5][1] 0.3979 1e-4
    @test_approx_eq_eps colmeans[6][1] 0.0325 1e-4

    colstds  = colwise(std, df)
    @test_approx_eq_eps colstds[1][1] 0.3205 1e-4
    @test_approx_eq_eps colstds[2][1] 0.1015 1e-4
    @test_approx_eq_eps colstds[3][1] 0.3692 1e-4
    @test_approx_eq_eps colstds[4][1] 0.1367 1e-4
    @test_approx_eq_eps colstds[5][1] 0.0937 1e-4
    @test_approx_eq_eps colstds[6][1] 0.0000 1e-4
end


text:
const condadir = abspath(first(DEPOT_PATH), "conda")
const condadeps = joinpath(condadir, "deps.jl")

module DefaultDeps
    if isfile("deps.jl")
        include("deps.jl")
    elseif isfile(Main.condadeps)
        include(Main.condadeps)
    end
    if !isdefined(@__MODULE__, :MINICONDA_VERSION)
        const MINICONDA_VERSION = "3"
    end
    if !isdefined(@__MODULE__, :ROOTENV)
        const ROOTENV = joinpath(Main.condadir, MINICONDA_VERSION)
    end
end

ROOTENV = get(ENV, "CONDA_JL_HOME", DefaultDeps.ROOTENV)
MINICONDA_VERSION = get(ENV, "CONDA_JL_VERSION", DefaultDeps.MINICONDA_VERSION)

if isdir(ROOTENV) && MINICONDA_VERSION != DefaultDeps.MINICONDA_VERSION
    error("""Miniconda version changed, since last build.
However, a root enviroment already exists at $(ROOTENV).
Setting Miniconda version is not supported for existing root enviroments.
To leave Miniconda version as, it is unset the CONDA_JL_VERSION enviroment variable and rebuild.
To change Miniconda version, you must delete the root enviroment and rebuild.
WARNING: deleting the root enviroment will delete all the packages in it.
This will break many Julia packages that have used Conda to install their dependancies.
These will require rebuilding.
""")
end

deps = """
const ROOTENV = "$(escape_string(ROOTENV))"
const MINICONDA_VERSION = "$(escape_string(MINICONDA_VERSION))"
"""

mkpath(condadir)
mkpath(ROOTENV)

for depsfile in ("deps.jl", condadeps)
    if !isfile(depsfile) || read(depsfile, String) != deps
        write(depsfile, deps)
    end
end


text:
# A Julia-vs-Jax game, from here:
# https://twitter.com/cgarciae88/status/1254269041784561665
# https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213/19
# https://gist.github.com/cgarciae/a69fa609f8fcd0aacece92660b5c2315

# These versions are updated to use Float32, which is what Jax is using (by default).
# It helped a little to avoid Float64 constants, and to add @inbounds in a few places.

using Pkg; pkg"add LoopVectorization Einsum TensorCast https://github.com/mcabbott/Tullio.jl"
using LoopVectorization, Tullio, Einsum, TensorCast, Test, BenchmarkTools

a = -100 .+ 200 .* rand(Float32, 5000, 2);
b = -100 .+ 200 .* rand(Float32, 5000, 2);

const None = [CartesianIndex()]

function distances(data1, data2)
    data1 = deg2rad.(data1)
    data2 = deg2rad.(data2)
    lat1 = @view data1[:, 1]
    lng1 = @view data1[:, 2]
    lat2 = @view data2[:, 1]
    lng2 = @view data2[:, 2]
    diff_lat = @view(lat1[:, None]) .- @view(lat2[None, :])
    diff_lng = @view(lng1[:, None]) .- @view(lng2[None, :])
    data = (
        @. sin(diff_lat / 2)^2 +
        cos(@view(lat1[:, None])) * cos(@view(lat2[None,:])) * sin(diff_lng / 2)^2
    )
    data .= @. 2.0 * 6373.0 * atan(sqrt(abs(data)), sqrt(abs(1.0 - data)))
    return reshape(data, (size(data1, 1), size(data2, 1)))
end

res = distances(a, b);
@test eltype(res) == Float32

function distances_threaded(data1, data2)
    lat1 = [deg2rad(data1[i,1]) for i in 1:size(data1, 1)]
    lng1 = [deg2rad(data1[i,2]) for i in 1:size(data1, 1)]
    lat2 = [deg2rad(data2[i,1]) for i in 1:size(data2, 1)]
    lng2 = [deg2rad(data2[i,2]) for i in 1:size(data2, 1)]
    # data = Matrix{Float64}(undef, length(lat1), length(lat2))
    data = Matrix{eltype(data1)}(undef, length(lat1), length(lat2))
    @inbounds Threads.@threads for i in eachindex(lat2)
        lat, lng = lat2[i], lng2[i]
        data[:, i] .= @. sin((lat1 - lat) / 2)^2 + cos(lat1) * cos(lat) * sin((lng1 - lng) / 2)^2
    end
    Threads.@threads for i in eachindex(data)
        # data[i] = 2.0 * 6373.0 * atan(sqrt(abs(data[i])), sqrt(abs(1.0 - data[i])))
        @inbounds data[i] = 2 * 6373 * atan(sqrt(abs(data[i])), sqrt(abs(1 - data[i])))
    end
    return data
end

function distances_threaded_simd(data1, data2) # @baggepinnen
    lat1 = [deg2rad(data1[i,1]) for i in 1:size(data1, 1)]
    lng1 = [deg2rad(data1[i,2]) for i in 1:size(data1, 1)]
    lat2 = [deg2rad(data2[i,1]) for i in 1:size(data2, 1)]
    lng2 = [deg2rad(data2[i,2]) for i in 1:size(data2, 1)]
    # data = Matrix{Float64}(undef, length(lat1), length(lat2))
    data = similar(data1, length(lat1), length(lat2))
    cos_lat1 = @avx cos.(lat1)
    Threads.@threads for i in eachindex(lat2)
        # lat, lng = lat2[i], lng2[i]
        @inbounds lat, cos_lat, lng = lat2[i], cos(lat2[i]), lng2[i]
        # @avx data[:, i] .= @. sin((lat1 - lat) / 2)^2 + cos(lat1) * cos(lat) * sin((lng1 - lng) / 2)^2
        @avx data[:, i] .= @. sin((lat1 - lat) / 2)^2 + cos_lat1 * cos_lat * sin((lng1 - lng) / 2)^2
    end
    Threads.@threads for i in eachindex(data)
        # @avx data[i] = 2.0 * 6373.0 * atan(sqrt(abs(data[i])), sqrt(abs(1.0 - data[i])))
        @avx data[i] = 2 * 6373 * atan(sqrt(abs(data[i])), sqrt(abs(1 - data[i])))
    end
    return data
end

@test res ≈ distances_threaded(a, b)
@test eltype(distances_threaded(a, b)) == Float32
@test res ≈ distances_threaded_simd(a, b)
@test eltype(distances_threaded_simd(a, b)) == Float32

function distances_bcast(data1, data2) # @DNF
    data1 = deg2rad.(data1)
    data2 = deg2rad.(data2)
    lat1 = @view data1[:, 1]
    lng1 = @view data1[:, 2]
    lat2 = @view data2[:, 1]
    lng2 = @view data2[:, 2]
    data = sin.((lat1 .- lat2') ./ 2).^2 .+ cos.(lat1) .* cos.(lat2') .* sin.((lng1 .- lng2') ./ 2).^2
    @. data = 2 * 6373 * atan(sqrt(abs(data)), sqrt(abs(1 - data)))
    return data
end

function distances_bcast_simd(data1, data2)
    data1 = deg2rad.(data1)
    data2 = deg2rad.(data2)
    lat1 = @view data1[:, 1]
    lng1 = @view data1[:, 2]
    lat2 = @view data2[:, 1]
    lng2 = @view data2[:, 2]
    @avx data = sin.((lat1 .- lat2') ./ 2).^2 .+ cos.(lat1) .* cos.(lat2') .* sin.((lng1 .- lng2') ./ 2).^2
    @. data = 2 * 6373 * atan(sqrt(abs(data)), sqrt(abs(1 - data)))
    return data
end

@test res ≈ distances_bcast(a, b)
@test eltype(distances_bcast(a, b)) == Float32
@test res ≈ distances_bcast_simd(a, b)
@test eltype(distances_bcast_simd(a, b)) == Float32

function distances_einsum(data1deg, data2deg)
    data1 = deg2rad.(data1deg)
    data2 = deg2rad.(data2deg)

    @einsum cd1[n] := cos(data1[n,1])
    @einsum cd2[m] := cos(data2[m,1])

    @einsum data[n,m] := sin((data1[n,1] - data2[m,1])/2)^2 +
        cd1[n] * cd2[m] * sin((data1[n,2] - data2[m,2])/2)^2

    @einsum data[n,m] = 2 * 6373 * atan(sqrt(abs(data[n,m])), sqrt(abs(1 - data[n,m])))
end

function distances_vielsum(data1deg, data2deg)
    data1 = deg2rad.(data1deg)
    data2 = deg2rad.(data2deg)

    @vielsum cd1[n] := cos(data1[n,1])
    @vielsum cd2[m] := cos(data2[m,1])

    @vielsum data[n,m] := sin((data1[n,1] - data2[m,1])/2)^2 +
        cd1[n] * cd2[m] * sin((data1[n,2] - data2[m,2])/2)^2

    @vielsum data[n,m] = 2 * 6373 * atan(sqrt(abs(data[n,m])), sqrt(abs(1 - data[n,m])))
end

@test res ≈ distances_einsum(a, b)
@test eltype(distances_einsum(a, b)) == Float32
@test res ≈ distances_vielsum(a, b)
@test eltype(distances_vielsum(a, b)) == Float32

function distances_cast(data1deg, data2deg)
    data1 = deg2rad.(data1deg)
    data2 = deg2rad.(data2deg)

    @cast cd1[n] := cos(data1[n,1]) # pulling these out is worth 25%
    @cast cd2[m] := cos(data2[m,1])

    @cast data[n,m] := sin((data1[n,1] - data2[m,1])/2)^2 +
        cd1[n] * cd2[m] * sin((data1[n,2] - data2[m,2])/2)^2

    @cast data[n,m] = 2 * 6373 * atan(sqrt(abs(data[n,m])), sqrt(abs(1 - data[n,m])))
end

function distances_cast_avx(data1deg, data2deg)
    data1 = deg2rad.(data1deg)
    data2 = deg2rad.(data2deg)

    @cast cd1[n] := cos(data1[n,1])  avx
    @cast cd2[m] := cos(data2[m,1])  avx

    @cast data[n,m] := sin((data1[n,1] - data2[m,1])/2)^2 +
        cd1[n] * cd2[m] * sin((data1[n,2] - data2[m,2])/2)^2  avx

    @cast data[n,m] = 2 * 6373 * atan(sqrt(abs(data[n,m])), sqrt(abs(1 - data[n,m])))  avx
end

@test res ≈ distances_cast(a, b)
@test eltype(distances_cast(a, b)) == Float32
@test res ≈ distances_cast_avx(a, b)
@test eltype(distances_cast_avx(a, b)) == Float32

function distances_tullio(data1deg, data2deg)
    data1 = deg2rad.(data1deg)
    data2 = deg2rad.(data2deg)

    @tullio data[n,m] := sin((data1[n,1] - data2[m,1])/2)^2 +
        cos(data1[n,1]) * cos(data2[m,1]) * sin((data1[n,2] - data2[m,2])/2)^2

    @tullio data[n,m] = 2 * 6373 * atan(sqrt(abs(data[n,m])), sqrt(abs(1 - data[n,m])))
end

# function distances_tullio2(data1deg, data2deg)
#     data1 = deg2rad.(data1deg)
#     data2 = deg2rad.(data2deg)

#     @tullio cd1[n] := cos(data1[n,1]) # has no effect
#     @tullio cd2[m] := cos(data2[m,1])

#     @tullio data[n,m] := sin((data1[n,1] - data2[m,1])/2)^2 +
#         cd1[n] * cd2[m] * sin((data1[n,2] - data2[m,2])/2)^2

#     @tullio data[n,m] = 2 * 6373 * atan(sqrt(abs(data[n,m])), sqrt(abs(1 - data[n,m])))
# end

@test res ≈ distances_tullio(a, b)
@test eltype(distances_tullio(a, b)) == Float32

##### laptop (2 cores, 4 threads)

julia> a = -100 .+ 200 .* rand(Float32, 5000, 2);
julia> b = -100 .+ 200 .* rand(Float32, 5000, 2);

julia> @btime distances($a, $b);
  1.522 s (26 allocations: 286.18 MiB)

julia> @btime distances_threaded($a, $b);
  516.937 ms (64 allocations: 95.45 MiB)

julia> @btime distances_threaded_simd($a, $b);
  215.938 ms (66 allocations: 95.47 MiB)

julia> @btime distances_bcast($a, $b);
  1.352 s (10 allocations: 95.44 MiB)

julia> @btime distances_bcast_simd($a, $b);
  641.506 ms (43 allocations: 95.44 MiB)

julia> @btime distances_einsum($a, $b);
  983.168 ms (10 allocations: 95.48 MiB)

julia> @btime distances_vielsum($a, $b);
  389.831 ms (103 allocations: 95.49 MiB)

julia> @btime distances_cast($a, $b); # unlike distances_bcast, this pulls out cos(...)
  1.034 s (16 allocations: 95.48 MiB)

julia> @btime distances_cast_avx($a, $b); # and this applies more @avx than bcast_simd
  137.557 ms (43 allocations: 190.85 MiB)

julia> @btime distances_tullio($a, $b);
  51.442 ms (636 allocations: 95.47 MiB)

##### desktop (6 cores, 12 threads)

julia> a = -100 .+ 200 .* rand(Float32, 5000, 2);
julia> b = -100 .+ 200 .* rand(Float32, 5000, 2);

julia> @btime distances($a, $b);
  1.166 s (26 allocations: 286.18 MiB)

julia> @btime distances_threaded($a, $b);
  140.062 ms (144 allocations: 95.46 MiB)

julia> @btime distances_threaded_simd($a, $b);
  64.382 ms (147 allocations: 95.48 MiB)

julia> @btime distances_bcast($a, $b);
  1.033 s (10 allocations: 95.44 MiB)

julia> @btime distances_bcast_simd($a, $b);
  501.002 ms (43 allocations: 95.44 MiB)

julia> @btime distances_einsum($a, $b);
  756.749 ms (10 allocations: 95.48 MiB)

julia> @btime distances_vielsum($a, $b);
  108.200 ms (262 allocations: 95.51 MiB)

julia> @btime distances_cast($a, $b); # unlike distances_bcast, this pulls out cos(...)
  795.199 ms (16 allocations: 95.48 MiB)

julia> @btime distances_cast_avx($a, $b); # and this applies more @avx than bcast_simd
  112.824 ms (43 allocations: 190.85 MiB)

julia> @btime distances_tullio($a, $b);
  28.151 ms (788 allocations: 95.48 MiB)

julia> a = -100 .+ 200 .* rand(Float64, 5000, 2); ##### repeat everythinng in Float64
julia> b = -100 .+ 200 .* rand(Float64, 5000, 2);

julia> @btime distances($a, $b);
  1.308 s (26 allocations: 572.36 MiB)

julia> @btime distances_threaded($a, $b);
  146.374 ms (144 allocations: 190.90 MiB)

julia> @btime distances_threaded_simd($a, $b);
  92.097 ms (146 allocations: 190.94 MiB)

julia> @btime distances_bcast($a, $b);
  1.134 s (10 allocations: 190.89 MiB)

julia> @btime distances_bcast_simd($a, $b);
  728.564 ms (43 allocations: 190.89 MiB)

julia> @btime distances_einsum($a, $b);
  874.312 ms (10 allocations: 190.96 MiB)

julia> @btime distances_vielsum($a, $b);
  123.725 ms (262 allocations: 191.00 MiB)

julia> @btime distances_cast($a, $b);
  902.447 ms (16 allocations: 190.96 MiB)

julia> @btime distances_cast_avx($a, $b);
  431.136 ms (43 allocations: 381.70 MiB)

julia> @btime distances_tullio($a, $b);
  75.608 ms (786 allocations: 190.93 MiB)

##### GPU (an ancient one!)

julia> using CuArrays, KernelAbstractions # and then re-run defn. of distances_tullio

julia> CuArrays.allowscalar(false)

julia> ca = cu(a); cb = cu(b); # Float32

julia> cres = distances_bcast(ca, cb);

julia> @test cres ≈ distances_tullio(ca, cb)
Test Passed

julia> @test cres ≈ distances_cast(ca, cb)
Test Passed

julia> @btime CuArrays.@sync distances_bcast($ca, $cb);
  31.558 ms (420 allocations: 18.42 KiB)

julia> @btime CuArrays.@sync distances_cast($ca, $cb);
  29.728 ms (546 allocations: 22.48 KiB)

julia> @btime CuArrays.@sync distances_tullio($ca, $cb);
  187.258 ms (173551 allocations: 2.66 MiB)

##### Python
# From here, verbatim:
# https://gist.github.com/cgarciae/a69fa609f8fcd0aacece92660b5c2315

import typing as tp
from jax import numpy as jnp
import jax
import numpy as np
import time
@jax.jit
def distances_jax(data1, data2):
    # data1, data2 are the data arrays with 2 cols and they hold
    # lat., lng. values in those cols respectively
    np = jnp
    data1 = np.deg2rad(data1)
    data2 = np.deg2rad(data2)
    lat1 = data1[:, 0]
    lng1 = data1[:, 1]
    lat2 = data2[:, 0]
    lng2 = data2[:, 1]
    diff_lat = lat1[:, None] - lat2
    diff_lng = lng1[:, None] - lng2
    d = (
        np.sin(diff_lat / 2) ** 2
        + np.cos(lat1[:, None]) * np.cos(lat2) * np.sin(diff_lng / 2) ** 2
    )
    data = 2 * 6373 * np.arctan2(np.sqrt(np.abs(d)), np.sqrt(np.abs(1 - d)))
    return data.reshape(data1.shape[0], data2.shape[0])
def distances_np(data1, data2):
    # data1, data2 are the data arrays with 2 cols and they hold
    # lat., lng. values in those cols respectively
    data1 = np.deg2rad(data1)
    data2 = np.deg2rad(data2)
    lat1 = data1[:, 0]
    lng1 = data1[:, 1]
    lat2 = data2[:, 0]
    lng2 = data2[:, 1]
    diff_lat = lat1[:, None] - lat2
    diff_lng = lng1[:, None] - lng2
    d = (
        np.sin(diff_lat / 2) ** 2
        + np.cos(lat1[:, None]) * np.cos(lat2) * np.sin(diff_lng / 2) ** 2
    )
    data = 2 * 6373 * np.arctan2(np.sqrt(np.abs(d)), np.sqrt(np.abs(1 - d)))
    return data.reshape(data1.shape[0], data2.shape[0])
a = np.random.uniform(-100, 100, size=(5000, 2)).astype(np.float32)
b = np.random.uniform(-100, 100, size=(5000, 2)).astype(np.float32)
def dist_np_test():
    return distances_np(a, b)
# enforce eager evaluation
def dist_jax_test():
    return distances_jax(a, b).block_until_ready()

##### Times on the same laptop as above:

Python 3.7.5 (default, Nov  6 2019, 19:41:43)
Type 'copyright', 'credits' or 'license' for more information
IPython 7.9.0 -- An enhanced Interactive Python. Type '?' for help.

In [2]: dist_np_test()
Out[2]:
array([[ 4011.349 , 11679.735 ,  1918.837 , ...,  2963.0593, 13176.956 ,
        15359.288 ],
       ...,
       [10144.612 , 18684.783 ,  6158.844 , ..., 10165.801 , 13639.45  ,
         8931.506 ]], dtype=float32)

In [3]: %timeit dist_np_test()
988 ms ± 3.14 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [4]: %timeit dist_jax_test()
/Users/me/.pyenv/versions/3.7.5/lib/python3.7/site-packages/jax/lib/xla_bridge.py:114: UserWarning: No GPU/TPU found, falling back to CPU.
  warnings.warn('No GPU/TPU found, falling back to CPU.')
372 ms ± 14.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [5]: jax.__version__
Out[5]: '0.1.50'

In [6]: np.__version__
Out[6]: '1.17.3'

# Closest versions above are probably these:
# distances_einsum 983.168 ms -- simple loops, single-threaded
# distances_vielsum 389.831 ms -- ditto, multi-threaded

##### Times on the same desktop as above:

Python 3.6.5 (default, Jun 17 2018, 12:13:06)
Type 'copyright', 'credits' or 'license' for more information
IPython 7.13.0 -- An enhanced Interactive Python. Type '?' for help.

In [3]: %timeit dist_np_test()
822 ms ± 4.09 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [4]: %timeit dist_jax_test()
/Users/me/code/jax19/.direnv/python-3.6.5/lib/python3.6/site-packages/jax/lib/xla_bridge.py:116: UserWarning: No GPU/TPU found, falling back to CPU.
  warnings.warn('No GPU/TPU found, falling back to CPU.')
107 ms ± 1.14 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [5]: jax.__version__
Out[5]: '0.1.65'

# Again these nearly exactly match the following:
# distances_einsum 756.749 ms
# distances_vielsum 108.200 ms
# and distances_tullio is still faster by a factor of 3.


text:
"""
    SegmentMetadata(dataSource::DataSource; <keyword arguments>)

Returns a Dict with the query output.

A method with all arguments as keyword arguments is also provided.

# Arguments
- intervals::Vector{<:Interval} = nothing
- toInclude::Union{String, Vector{String}} = nothing
- merge::Bool = nothing
- analysisTypes::Vector{String} = nothing
- lenientAggregatorMerge::Bool = nothing
- virtualColumns::Vector{<:VirtualColumn} = nothing
- context::Dict = nothing
"""
mutable struct SegmentMetadata <: Query
    queryType::String
    dataSource::DataSource
    intervals
    toInclude
    merge
    analysisTypes
    lenientAggregatorMerge
    virtualColumns
    context
    function SegmentMetadata(
        dataSource; intervals=nothing, toInclude=nothing, merge=nothing,
        analysisTypes=nothing, lenientAggregatorMerge=nothing, virtualColumns=nothing, context=nothing
    )
        nothing_or_type(intervals, Vector{<:Interval})
        nothing_or_type(toInclude, Union{String, Vector{String}})
        nothing_or_type(merge, Bool)
        analysisTypes === nothing || (isa(analysisTypes, Vector{String}) &&
            all(at -> at ∈ ["cardinality", "minmax", "size", "interval", "timestampSpec", "queryGranularity", "aggregators", "rollup"], analysisTypes)) ||
            error("Invalid analysisTypes")
        nothing_or_type(lenientAggregatorMerge, Bool)
        nothing_or_type(virtualColumns, Vector{<:VirtualColumn})
        nothing_or_type(context, Dict)
        if isa(toInclude, String)
            toInclude ∈ ["all", "none"] || error("Invalid toInclude")
            toInclude = Dict("type" => toInclude)
        else
            toInclude = Dict("type" => "list", "columns" => toInclude)
        end
        new("segmentMetadata", dataSource, intervals, toInclude, merge, analysisTypes, lenientAggregatorMerge, virtualColumns, context)
    end
end
SegmentMetadata(
    ; dataSource, intervals=nothing, toInclude=nothing, merge=nothing,
    analysisTypes=nothing, lenientAggregatorMerge=nothing, virtualColumns=nothing, context=nothing
) = SegmentMetadata(
    dataSource; intervals, toInclude, merge,
    analysisTypes, lenientAggregatorMerge, virtualColumns, context
)

function execute(client::Client, query::SegmentMetadata; pretty=false)
    res = JSON.parse(execute_native_query(client, query; pretty))
    res[1]
end


text:
import Base.show
using CUDA: CuArray

# There is a data loader PR: https://github.com/FluxML/Flux.jl/pull/450
mutable struct DataSetIterator
    raw_x::AbstractArray
    raw_y::AbstractArray
    index::Int
    batch_size::Int
    nbatch::Int
end
function DataSetIterator(x,y,batch_size)
    DataSetIterator(x, y, 1, batch_size,
                    convert(Int, floor(size(x)[end]/batch_size)))
end

mutable struct CuDataSetIterator
    gpu_x::Array{CuArray}
    gpu_y::Array{CuArray}
    index::Int
    batch_size::Int
    nbatch::Int
end

function merge_dses(dses)
    # merge datasets
    # 1. merge raw_x and raw_y
    xdim = length(size(dses[1].raw_x))
    ydim = length(size(dses[1].raw_y))
    @show xdim, ydim
    # 4
    raw_x = cat([ds.raw_x for ds in dses]..., dims=xdim)
    # 3
    raw_y = cat([ds.raw_y for ds in dses]..., dims=ydim)
    # 2. assert batch_size
    batch_size = dses[1].batch_size
    # 4. construct new ds
    DataSetIterator(raw_x, raw_y, batch_size)
end

function Base.convert(::Type{CuDataSetIterator}, ds::DataSetIterator)
    CuDataSetIterator(ds.raw_x, ds.raw_y, ds.batch_size)
end

function CuDataSetIterator(ds::DataSetIterator)
    CuDataSetIterator(ds.raw_x, ds.raw_y, ds.batch_size)
end

function CuDataSetIterator(x,y,batch_size)
    xdims = collect(repeated(:, length(size(x)) - 1))
    ydims = collect(repeated(:, length(size(y)) - 1))

    N = size(x)[end]
    nbatch = Int(floor(N/batch_size))

    # reshuffle
    indices = shuffle(1:size(x)[end])
    x = x[xdims...,indices]
    y = y[ydims...,indices]

    gpu_x = map(1:nbatch) do i
        start = (i-1) * batch_size + 1
        stop = i * batch_size
        x[xdims..., start:stop]
    end
    gpu_y = map(1:nbatch) do i
        start = (i-1) * batch_size + 1
        stop = i * batch_size
        y[ydims..., start:stop]
    end
    # FIXME .|>
    gpu_x = gpu_x |> gpu
    gpu_y = gpu_y |> gpu
    CuDataSetIterator(gpu_x, gpu_y, 1, batch_size, nbatch)
end

Base.show(io::IO, x::DataSetIterator) = begin
    println(io, "DataSetIterator:")
    println(io, "  batch size: ", x.batch_size)
    println(io, "  number of batches: ", x.nbatch)
    println(io, "  x data shape: ", size(x.raw_x))
    println(io, "  y data shape: ", size(x.raw_y))
    print(io, "  current index: ", x.index)
end

Base.show(io::IO, x::CuDataSetIterator) = begin
    println(io, "CuDataSetIterator:")
    println(io, "  batch size: ", x.batch_size)
    println(io, "  number of batches: ", x.nbatch)
    println(io, "  x data shape: ", size(x.gpu_x[1]))
    println(io, "  y data shape: ", size(x.gpu_y[1]))
    print(io, "  current index: ", x.index)
end

function next_batch!(ds::DataSetIterator)
    if ds.index > ds.nbatch
        ds.index = 1
    end

    xdims = collect(repeated(:, length(size(ds.raw_x)) - 1))
    ydims = collect(repeated(:, length(size(ds.raw_y)) - 1))

    if ds.index == 1
        indices = shuffle(1:size(ds.raw_x)[end])
        ds.raw_x = ds.raw_x[xdims...,indices]
        ds.raw_y = ds.raw_y[ydims...,indices]
    end

    a = (ds.index - 1) * ds.batch_size + 1
    b = ds.index * ds.batch_size

    ds.index += 1
    # FIXME use view instead of copy
    # FIXME properly slice with channel-last format
    return ds.raw_x[xdims...,a:b], ds.raw_y[ydims...,a:b]
end

function next_batch!(ds::CuDataSetIterator)
    if ds.index > ds.nbatch
        ds.index = 1
    end
    ds.index += 1
    return ds.gpu_x[ds.index-1], ds.gpu_y[ds.index-1]
end

_next_batch_step = 0
function next_batch!(dses::Array{T, N}
                     where T<:Union{DataSetIterator, CuDataSetIterator}
                     where N)
    # keep an internal step?
    global _next_batch_step
    _next_batch_step += 1
    next_batch!(dses[_next_batch_step % length(dses) + 1])
end

"""
1. shuffle and divide *all* data into batches
2. move *all* data to GPU
3. each time next_batch! is called, a pointer to the GPU array is returned
"""
function move_to_gpu!(ds::DataSetIterator)
    xdims = collect(repeated(:, length(size(ds.raw_x)) - 1))
    ydims = collect(repeated(:, length(size(ds.raw_y)) - 1))

    N = size(ds.raw_x)[end]

    # reshuffle
    indices = shuffle(1:size(ds.raw_x)[end])
    ds.raw_x = ds.raw_x[xdims...,indices]
    ds.raw_y = ds.raw_y[ydims...,indices]

    gpu_x = map(1:Int(floor(N/ds.batch_size))) do i
        start = (i-1) * ds.batch_size + 1
        stop = i * ds.batch_size
        ds.raw_x[xdims..., start:stop]
    end
    gpu_y = map(1:Int(floor(N/ds.batch_size))) do i
        start = (i-1) * ds.batch_size + 1
        stop = i * ds.batch_size
        ds.raw_y[ydims..., start:stop]
    end
    ds.gpu_x = gpu_x |> gpu
    ds.gpu_y = gpu_y |> gpu
    ds.gpu_index = 1
    nothing
end

max_stars_repo_name: lihebi/DAG-EQ

max_issues_repo_name: lihebi/DAG-EQ

max_forks_repo_name: lihebi/DAG-EQ


text:
module Sdr

import ..Clarus

function lastprice(;params...)
  return Clarus.api_request("SDR","LastPrice";params...)
end

function lasttrade(;params...)
  return Clarus.api_request("SDR","LastTrade";params...)
end

function tradefilter(;params...)
  return Clarus.api_request("SDR","TradeFilter";params...)
end

function trades(;params...)
  return Clarus.api_request("SDR","Trades";params...)
end

function volume(;params...)
  return Clarus.api_request("SDR","Volume";params...)
end

function volumeadv(;params...)
  return Clarus.api_request("SDR","VolumeADV";params...)
end

end #MODULE-END


text:
"""This filter returns the `match` result if the level in the `Event` is the same or more specific than the configured level and the onMismatch value otherwise.
For example, if the `ThresholdFilter` is configured with `Level.ERROR` and the `Event` contains `Level.DEBUG` then the `mismatch` value will be returned since `ERROR` events are more specific than `DEBUG`.

Default values: `level` is `ERROR`, `match` is `NEUTRAL`, `mismatch` is `DENY`.
"""
type ThresholdFilter <: Filter
    level::Level.EventLevel
    match::FilterResult.ResultType
    mismatch::FilterResult.ResultType
    state::LifeCycle.State
    ThresholdFilter(lvl::Level.EventLevel, match::FilterResult.ResultType, mismatch::FilterResult.ResultType) =
        new(lvl, match, mismatch, LifeCycle.INITIALIZED)
end
ThresholdFilter(lvl::Level.EventLevel) = ThresholdFilter(lvl, FilterResult.NEUTRAL, FilterResult.DENY)
ThresholdFilter() = ThresholdFilter(Level.ERROR)

function ThresholdFilter(config::Dict{AbstractString,Any})
    lvl = if haskey(config, "level")
        get(getlevel(config["level"]), Level.ERROR)
    else
        Level.ERROR
    end
    rmatch = if haskey(config, "match")
        evaltype(config["match"], "FilterResult")
    else
        FilterResult.NEUTRAL
    end
    rmismatch = if haskey(config, "mismatch")
        evaltype(config["mismatch"], "FilterResult")
    else
        FilterResult.DENY
    end
    ThresholdFilter(lvl, rmatch, rmismatch)
end
ThresholdFilter(;kwargs...) = map(e->(string(e[1]),e[2]), kwargs) |> Dict{AbstractString,Any} |> ThresholdFilter

"Make the filter available for use."
function start(flt::ThresholdFilter)
    state!(flt, LifeCycle.STARTED)
end

"Disable the filter."
function stop(flt::ThresholdFilter)
    state!(flt, LifeCycle.STOPPED)
end

filter{E <: Event}(flt::ThresholdFilter, evnt::E) = filter(flt, level(evnt))
filter(flt::ThresholdFilter, level::Level.EventLevel, marker::MARKER, msg) = filter(flt, level)
filter(flt::ThresholdFilter, lvl::Level.EventLevel) = flt.level >= lvl ? flt.match : flt.mismatch


text:
using SubstitutionTilings
using SubstitutionTilings.Pinwheel

using Luxor

using Profile

K = Pinwheel.K

w = 800
h = 800
sc = 30
@png begin
    colors = ["#DD93FC", "#E7977A",]

    first_tile = wheel()
    tiling = Pinwheel.PinwheelElem(0,0,0,-i*20-20)*substitute(pinwheel(), [first_tile], 5)
    setline(1)

    for tile in tiling
        origin()
        draw(tile, sc, colors[Pinwheel.color(tile)], :fill)
        origin()
        draw(tile, sc, "black", :stroke)
    end
    #setline(2)
    #draw(first_tile, sc, "black", :stroke)
end w h "pinwheel-tiling.png"

@png begin
    w = 800
    h = 800
    sc = 100
    colors = ["#DD93FC", "#E7977A",]
    quadrants = Tiler(width, height, 2, 1, margin=5)

    for (pos, n) in quadrants
        first_tile = Pinwheel.PinwheelElem(0,0,0,-1//2*i-1//2)*wheel()
        tiling = substitute(pinwheel(), [first_tile], n-1)
        for tile in tiling
            origin()
            translate(pos)
            scale(sc)
            sethue(colors[Pinwheel.color(tile)])
            transform(embed_aff(tile[1]))
            draw(tile[2], :fill)
            origin()
            translate(pos)
            scale(sc)
            setline(0.5)
            sethue("black")
            transform(embed_aff(tile[1]))
            draw(tile[2], :stroke)
        end
    end
end width height "pinwheel-subst.png"

vertex_star = Dict([wheel(), wheel(0,0,1,K(0)), wheel(0,1,0,-i)])
initial_collar = collar_in(Dict(Pinwheel.PinwheelElem(0,0,0,-i-1)*substitute(pinwheel(), [first_tile], 2)), Pinwheel.PinwheelElem(0,0,0,K(0)))

vertex_star_2 = Dict([wheel(), wheel(0,0,1,K(0)), wheel(0,2,0,K(0)), wheel(0,2,1,K(0))])
@time frequency(pinwheel(), initial_collar, vertex_star,2)
@time frequency(pinwheel(), initial_collar, vertex_star_2,2) # 4x bigger than in Baake-Grimm because every tile gets counted

@png begin
    colors = ["#DD93FC", "#E7977A",]

    tiling = vertex_star
    setline(1)

    for tile in tiling
        origin()
        draw(tile, sc, colors[Pinwheel.color(tile)], :fill)
        origin()
        draw(tile, sc, "black", :stroke)
    end
    setline(2)
    #draw(first_tile, sc, "black", :stroke)
end w h "vertex_star_1.png"

@png begin
    colors = ["#DD93FC", "#E7977A",]

    tiling = vertex_star_2
    setline(1)

    for tile in tiling
        origin()
        draw(tile, sc, colors[Pinwheel.color(tile)], :fill)
        origin()
        draw(tile, sc, "black", :stroke)
    end
    setline(2)
    #draw(first_tile, sc, "black", :stroke)
end w h "vertex_star_2.png"


text:
using Test
using SymDiff

@test subs(:(1 + a + b), :a => 1) == :(1 + 1 + b)


text:
module Leaf
using PhysCon
using WaterVaporMod

using Parameters

export leaf_params, setLeafT!, BallBerry

# Scaling functions for Photosynthesis temperature response and inhibition
ft(tl, ha) = exp(ha/(physcon.Rgas*(physcon.tfrz+25)) * (1-(physcon.tfrz+25)/tl));
fth(tl, hd, se, fc) = fc / (1 + exp((-hd+se*tl)/(physcon.Rgas*tl)));
fth25(hd, se) = 1.0 + exp( (-hd + se * (physcon.tfrz+25.)) / (physcon.Rgas * (physcon.tfrz+25.)) );

# Structure with all parameter temperature dependencies of a Leaf (largely based on Bonan's Photosynthesis model but exported as struct)
@with_kw mutable struct leaf_params{TT<:Number}
    # Rate constants (arbitrary units here, only relative rates are important):
    Kf::TT = 0.05;                           # Rate constant for fluorescence (might try to fit this eventually)
    Kd::TT = 0.85;                           # Rate constant for thermal dissipation
    Kp::TT = 4.0;                            # Rate constant for photochemistry (all reaction centers open)
    Kn::TT = 0;                              # NPQ rate constant (initially zero)
    Kn_ss::TT = 0;                           # Kn in steady state
    maxPSII::TT = Kp/(Kp+Kf+Kd);             # Max PSII yield (Kn=0, all RC open)
    "Flexas et al derived Fluorescence model params"
    Knparams = [5.01, 1.93, 10.0]
    Fo::TT = 0;     # dark-adapted fluorescence yield Fo,0
    Fo′::TT = 0;     # light-adapted fluorescence yield in the dark Fo
    Fm::TT = 0;     # light-adapted fluorescence yield Fm
    Fm′::TT = 0;    # dark-adapted fluorescence yield Fm
    ϕs::TT = 0;     # steady-state (light-adapted) yield Ft (aka Fs)
    eta::TT = 0;
    qQ::TT = 0;     # photochemical quenching
    qE::TT = 0;     # non-photochemical quenching
    NPQ::TT = 0;
    "Dynamic Steady state boolean"
    dynamic_state::Bool = false

    effcon::TT = 1/5                       # [mol CO2/mol e-]  number of CO2 per electrons - typically 1/5 for C3 and 1/6 for C4, i.e. about 1/10 for both PSI and PSII!
    CO2_per_electron::TT = 1/5             # Similar to above but accounting for Photorespiration as well (which "steals" electrons)
    Ja::TT = 0                               # Actual electron transport rate (needed to compute fluorescence later on)

    C3::Bool = true                           # C3 (or C4) plant
    use_colim::Bool = true                    # Use co-limitation
    gstyp::Int = 1                            # stomatal conductance type (1=Ball Berry, 0 = Medlyn)
    o₂::TT = 0.209e3;                         # Standard O2 in mmol/mol
    T::TT = 298.15;                           # standard temperature (25C)
    esat::TT = 0.0
    desat::TT = 0.0
    #(esat, desat) = SatVap(T)
    # Curvature parameter:
    θ_j::TT = 0.90;

    # Conductances:
    # Ball Berry or Medlyn model
    g0::TT = 0.01;                          # Ball-Berry minimum leaf conductance, unstressed (mol H2O/m2/s)
    g1::TT = 9;                             # Ball-Berry slope of conductance-photosynthesis relationship, unstressed
    vpd_min = 100.0                            # Min VPD for Medlyn model (blows up otherwise?)
    # Add already here: Mesophyll conductance (Inf for now):
    gm::TT = Inf;                              # Mesophyll conductance (μmol/m2/s/Pa)
    gs::TT = 0.1;                              # Stomatal conductance (μmol/m2/s/Pa); just using some prior
    gs_ss::TT = 0.1;                           # Steady state Stomatal conductance (μmol/m2/s/Pa);

    #---------------------------------------------------------------------
    # kc, ko, cp at 25C: Bernacchi et al (2001) Plant, Cell Environment 24:253-259
    # Derive sco from cp with o2=0.209 mol/mol and re-calculate Γ to allow
    # variation in O2
    #---------------------------------------------------------------------
    kc_25::TT = 404.9;                         # Michaelis-Menten constant for CO2 at 25C μmol/mol
    ko_25::TT = 278.4;                         # Michaelis-Menten constant for O2 at 25C  mmol/mol
    Γ_25_::TT = 42.75;                         # CO2 compensation point at 25C μmol/mol
    sco::TT = 0.5 * 0.209 / (Γ_25_ * 1.e-06);  # Γ_25 (μmol/mol) -> (mol/mol)
    Γ_25::TT = 0.5 * o₂ / sco * 1000.;         # O2 is mmol/mol. Multiply by 1000 for μmol/mol

    #---------------------------------------------------------------------
    # Activation energy:
    # Bernacchi et al (2001) Plant, Cell Environment 24:253-259
    # Bernacchi et al (2003) Plant, Cell Environment 26:1419-1430
    # Acclimation from: Kattge and Knorr (2007) Plant, Cell Environment 30:1176-1190
    #---------------------------------------------------------------------
    kcha::TT    = 79430.;                      # Activation energy for kc (J/mol)
    koha::TT    = 36380.;                      # Activation energy for ko (J/mol)
    cpha::TT    = 37830.;                      # Activation energy for cp (J/mol)
    vcmaxha::TT = 65330.;                      # Activation energy for Vcmax (J/mol)
    jmaxha::TT  = 43540.;                      # Activation energy for Jmax (J/mol)
    rdha::TT    = 46390.;                      # Activation energy for Rd (J/mol)

    #---------------------------------------------------------------------
    # High temperature deactivation:
    # Leuning (2002) Plant, Cell Environment 25:1205-1210
    # The factor "c" scales the deactivation to a value of 1.0 at 25C
    # Acclimation from: Kattge and Knorr (2007) Plant, Cell Environment 30:1176-1190
    #---------------------------------------------------------------------
    vcmaxhd::TT = 150000.;                    # Deactivation energy for Vcmax (J/mol)
    jmaxhd::TT  = 150000.;                    # Deactivation energy for Jmax (J/mol)
    rdhd::TT    = 150000.;                    # Deactivation energy for Rd (J/mol)

    vcmaxse::TT = 490.;                       # Entropy term for Vcmax (J/mol/K)
    jmaxse::TT  = 490.;                       # Entropy term for Jmax (J/mol/K)
    rdse::TT    = 490.;                       # Entropy term for Rd (J/mol/K)

    vcmaxc::TT = fth25(vcmaxhd, vcmaxse);    # Scaling factor for high temperature inhibition (25 C = 1.0)
    jmaxc::TT  = fth25(jmaxhd, jmaxse);      # Scaling factor for high temperature inhibition (25 C = 1.0)
    rdc::TT    = fth25(rdhd, rdse);          # Scaling factor for high temperature inhibition (25 C = 1.0)

    # Initialize at standard T:
    kc::TT  = kc_25;
    ko::TT  = ko_25;
    Γstar::TT  = Γ_25;

    # Use some standard values first
    vcmax25::TT = 80.0;                        # Leaf maximum carboxylation rate at 25C for canopy layer (μmol/m2/s)
    jmax25::TT  = 1.97*vcmax25;                # C3 - maximum electron transport rate at 25C for canopy layer (μmol/m2/s)
    rd25::TT    = 0.7605;                      # Leaf respiration rate at 25C for canopy layer (μmol CO2/m2/s)
    vcmax::TT  = vcmax25;
    jmax::TT   = jmax25;
    rdleaf::TT = rd25;

    #
    Ci::TT = 0.0                               # CO2 concentration in mesophyll/internal

    # to be computed
    je::TT = 0.0                               # electron transport rate
end

# Set Leaf rates with vcmax, jmax and rd at 25C as well as actual T here:
# For some reason, this is slow and allocates a lot, can be improved!!
"Set Leaf rates with vcmax, jmax and rd at 25C as well as actual T here"
function setLeafT!(l::leaf_params,  T::Number)
    l.T = T;
    l.kc     = l.kc_25          * ft(T, l.kcha);
    l.ko     = l.ko_25          * ft(T, l.koha);
    l.Γstar  = l.Γ_25           * ft(T, l.cpha);
    l.vcmax   = l.vcmax25 * ft(T, l.vcmaxha) * fth(T, l.vcmaxhd, l.vcmaxse, l.vcmaxc);
    l.jmax    = l.jmax25  * ft(T, l.jmaxha)  * fth(T, l.jmaxhd, l.jmaxse, l.jmaxc);
    l.rdleaf  = l.rd25    * ft(T, l.rdha)    * fth(T, l.rdhd, l.rdse, l.rdc);
    (l.esat, l.desat) = SatVap(T);
    # l.kd = max(0.8738,  0.0301*(T-273.15)+ 0.0773); # Can implement that later.
end

# Ball-Berry stomatal conductance model:
function BallBerry(Cs, RH, A, l::leaf_params)
#  Cs  : CO2 at leaf surface
#  RH  : relative humidity
#  A   : Net assimilation in 'same units of CO2 as Cs'/m2/s
#  BallBerrySlope, BallBerry0,
#  minCi : minimum Ci as a fraction of Cs (in case RH is very low?)
#  Ci_input : will calculate gs if A is specified.

gs = l.BB_g1 * A * RH/Cs  + l.BB_g0;
end # function

end #Module

max_stars_repo_name: ruillercoaster/LSM-SPAM

max_issues_repo_name: ruillercoaster/LSM-SPAM

max_forks_repo_name: ruillercoaster/LSM-SPAM


text:
include("a.jl")
include("g.jl")

include("a_p.jl")
include("epsilon_p.jl")
include("kappa_p.jl")

include("kappa_tau.jl")

include("alpha_i_star.jl")
include("alpha_n_star.jl")
include("alpha_r_star.jl")

include("epsilon_b.jl")
include("r_cs.jl")
include("h_cs.jl")
include("b.jl")
include("c.jl")
include("d.jl")

include("perimeter.jl")
include("g_p.jl")

include("kappa_x.jl")
include("delta_x.jl")

include("cross_sectional_area.jl")
include("surface_area.jl")
include("volume.jl")


text:
cd(@__DIR__)
using Revise, DifferentialEvolutionMCMC, Random, Parameters
using SequentialSamplingModels, Distributions
Random.seed!(88484)

dist = LBA(ν=[3.0,2.0], A = .8, k = .2, τ = .3)
choice,rt = rand(dist, 100)
min_rt = minimum(rt)

function _prior_loglike(ν, A, k, τ, min_rt)
    LL = 0.0
    LL += sum(logpdf(Normal(1, 5), ν))
    LL += logpdf(Normal(.8, .2), A)
    LL += logpdf(Normal(.2, .1), k)
    LL += logpdf(Uniform(0, min_rt), τ)
    return LL
end

function sample_prior()
    ν = rand(Normal(1, 5), 2)
    A = rand(Normal(.8, .2))
    k = rand(Normal(.2, .1))
    τ = rand(Uniform(0, min_rt))
    return as_union([ν,A,k,τ])
end

wrapper(min_rt) = (ν, A, k, τ) -> _prior_loglike(ν, A, k, τ, min_rt)

prior_loglike = wrapper(min_rt)

bounds = ((0.0,Inf), (0.0,Inf), (0.0,Inf), (0.0,min_rt))
names = (:ν, :A, :k, :τ)

function loglike(data, ν, A, k, τ)
    dist = LBA(;ν, A, k, τ)
    return logpdf.(dist, data...) |> sum
end

model = DEModel(;
    sample_prior,
    prior_loglike,
    loglike,
    data = (choice,rt),
    names
)

de = DE(;sample_prior, sample_prior, bounds, burnin=1500, n_groups=3, Np=15)
n_iter = 3000
chain = sample(model, de, MCMCThreads(), n_iter, progress=true)


text:
import HierarchicalUtils: NodeType, children, InnerNode, LeafNode, printtree, noderepr

# for schema structures
NodeType(::Type{<:Union{Nothing, Entry}}) = LeafNode()  # because sometimes we have empty array extractor
NodeType(::Type{Pair{Symbol,Any}}) = InnerNode()  # because sometimes we have empty array extractor
NodeType(::Type{<:Union{ArrayEntry, DictEntry, MultiEntry}}) = InnerNode()

noderepr(n::Nothing) = "Nothing"
noderepr(n::Entry) = "[Scalar - $(join(sort(string.(types(n))), ","))], $(length(keys(n.counts))) unique values, updated = $(n.updated)"
noderepr(n::ArrayEntry) = "[" * (isnothing(n.items) ? "Empty " : "") * "List] (updated = $(n.updated))"
noderepr(n::DictEntry) = "[" * (isnothing(n.childs) ? "Empty " : "") * "Dict] (updated = $(n.updated))"
noderepr(n::MultiEntry) = "[" * (isempty(n.childs) ? "Empty " : "") * "MultiEntry] (updated = $(n.updated))"

children(n::ArrayEntry) = (n.items,)
# using vector of pairs because splatting to named tuple is not good for compiler
children(n::DictEntry) = collect(n.childs)
children(n::MultiEntry) = [Symbol(k) => v for (k,v) in enumerate(n.childs)]

# for extractor structures
# default extractor
NodeType(::Type{<:AbstractExtractor}) = LeafNode()
NodeType(::Type{<:Union{ExtractArray, ExtractDict, MultipleRepresentation, ExtractKeyAsField, AuxiliaryExtractor}}) = InnerNode()

noderepr(::T) where T <: AbstractExtractor = "$(Base.nameof(T))"
noderepr(n::ExtractArray) = "Array of"
noderepr(n::ExtractDict) = "Dict"
noderepr(n::ExtractCategorical) = "Categorical d = $(n.n)"
noderepr(n::ExtractScalar{T}) where {T} = "$(T)"
noderepr(n::ExtractString) = "String"
noderepr(n::ExtractVector) = "FeatureVector with $(n.n) items"
noderepr(n::MultipleRepresentation) = "MultiRepresentation"
noderepr(e::ExtractKeyAsField) = "KeyAsField"
noderepr(n::AuxiliaryExtractor) = "Auxiliary extractor with"

children(n::ExtractArray) = (n.item,)
children(n::Pair) = (n.second,)
children(n::MultipleRepresentation) = n.extractors
children(e::ExtractKeyAsField) = (e.key, e.item)
children(n::AuxiliaryExtractor) = (n.extractor,)
children(n::ExtractDict) = (; Dict(Symbol(k)=>v for (k,v) in n.dict)...)


text:
struct VectorDict{K,V} <: AbstractDict{K,V}
    data::Vector{Pair{K,V}}
end
VectorDict{K,V}() where {K,V} = VectorDict(Vector{Pair{K,V}}())

Base.keys(vd::VectorDict) = first.(vd.data)
Base.values(vd::VectorDict) = last.(vd.data)
Base.haskey(vd::VectorDict, key) = !isempty(get_all_matches(vd, key))

function Base.get(vd::VectorDict, key, default)
    matches = get_all_matches(vd, key)
    if !isempty(matches)
        first(matches)
    else
        default
    end
end

"""
    get_all_matches(vd::VectorDict, key)
Returns all values matching the key
"""
function get_all_matches(vd::VectorDict, key)
    (vv for (kk, vv) in vd.data if isequal(key, kk))
end

function Base.setindex!(vd::VectorDict, value, key)
    # insert it at the front so newest things are matched first
    # old ones are left behind,
    pushfirst!(vd.data, key=>value)
end

Base.iterate(vd::VectorDict) = iterate(vd.data)
Base.iterate(vd::VectorDict, state) = iterate(vd.data, state)


text:
function objective_min_cost_TNEP(pm::_PM.AbstractPowerModel)
    if pm.setting["FSprotection"] == true || pm.setting["NSprotection"] == true
        objective_min_cost_TNEP_FSNS(pm)
    elseif  pm.setting["Permanentloss"] == true
        objective_min_cost_TNEP_PL(pm)
    end
end

function objective_min_cost_TNEP_nocl(pm::_PM.AbstractPowerModel)
    if pm.setting["FSprotection"] == true || pm.setting["NSprotection"] == true
        objective_min_cost_TNEP_FSNS_nocl(pm)
    elseif  pm.setting["Permanentloss"] == true
        objective_min_cost_TNEP_PL_nocl(pm)
    end
end

function objective_min_cost_TNEP_FSNS_nocl(pm::_PM.AbstractPowerModel)
        base_nws = pm.setting["base_list"]
        cont_nws = pm.setting["Cont_list"]
        Total_sample = pm.setting["Total_sample"]
        curt_gen = pm.setting["curtailed_gen"]
        max_curt = pm.setting["max_curt"]
        year_base = pm.setting["year_base"]
        total_year = pm.setting["total_yr"]

        gen_cost = Dict()
        FFR_cost = Dict()
        FCR_cost = Dict()
        fail_prob = Dict()
        Inv_cost = _PM.var(pm)[:Inv_cost] = JuMP.@variable(pm.model, start = 0, lower_bound = 0)
        Gen_cost = _PM.var(pm)[:Gen_cost] = JuMP.@variable(pm.model, start = 0, lower_bound = 0)
        FFRReserves = _PM.var(pm)[:FFR_Reserves] = JuMP.@variable(pm.model, start = 0, lower_bound = 0)
        FCRReserves = _PM.var(pm)[:FCR_Reserves] = JuMP.@variable(pm.model, start = 0, lower_bound = 0)
        Cont = _PM.var(pm)[:Cont]     = JuMP.@variable(pm.model, start = 0, lower_bound = 0)
        Curt = _PM.var(pm)[:curt] = JuMP.@variable(pm.model, [i in 1:total_year], start = 0, lower_bound = 0)

        sol_component_value_mod_wonw(pm,   :Inv_cost, Inv_cost)
        sol_component_value_mod_wonw(pm,   :Gen_cost, Gen_cost)
        sol_component_value_mod_wonw(pm,   :Curt, Curt)
        sol_component_value_mod_wonw(pm,   :FFR_Reserves, FFRReserves)
        sol_component_value_mod_wonw(pm,   :FCR_Reserves, FCRReserves)
        sol_component_value_mod_wonw(pm,   :Cont, Cont)

        for (n, nw_ref) in _PM.nws(pm)
            for (r, reserves) in nw_ref[:reserves]
                FFR_cost[(n,r)] = reserves["Cf"]
                FCR_cost[(n,r)] = reserves["Cg"]
            end
            for (i,gen) in nw_ref[:gen]
                pg = _PM.var(pm, n, :pg, i)
                if length(gen["cost"]) == 1
                    gen_cost[(n,i)] = gen["cost"][1]
                elseif length(gen["cost"]) == 2
                    gen_cost[(n,i)] = gen["cost"][1]*pg + gen["cost"][2]
                elseif length(gen["cost"]) == 3
                    gen_cost[(n,i)] = gen["cost"][2]*pg + gen["cost"][3]
                else
                    gen_cost[(n,i)] = 0.0
                end
            end
        end

        FFR_list =  Dict()
        FCR_list = Dict()
        for (b,c,br) in cont_nws
            FFR_list[(b,c)] = FFR_cost[(c,2)]*100/10^6*_PM.var(pm, c, :Pff, 2)+
                              FCR_cost[(c,2)]*100/10^6*_PM.var(pm, c, :Pgg, 2)
            # a_ub, a_lb = _IM.variable_domain(_PM.var(pm, c, :Pff, 2))
            # display("a_ub:$a_ub a_lb:$a_lb")
        end

        Zff = _PM.var(pm)[:Zff] = Dict((b, c) => JuMP.@variable(pm.model, base_name = "zff[$(string(b)),$(string(c))]", binary = true, start = 0) for (b, c, br) in cont_nws)
        sol_component_value_mod_wonw(pm, :Zff, Zff)

        #maximum out of onshore converters
        FFRReserves_max = _PM.var(pm)[:FFRReserves_max] = JuMP.@variable(pm.model, [b in base_nws], base_name ="FFRReserves_max", start = 0)
        FCRReserves_max = _PM.var(pm)[:FCRReserves_max] = JuMP.@variable(pm.model, [b in base_nws], base_name ="FCRReserves_max", start = 0)
        for b in base_nws
            sol_component_value_mod(pm, b, :FFRReserves_max, FFRReserves_max[b])
            sol_component_value_mod(pm, b, :FCRReserves_max, FCRReserves_max[b])
        end

        base_cont= Dict()

        for b in base_nws
            # [item for item in a if item[2] == 4]
            base_cont = [tt for tt in cont_nws if tt[1] == b]

            ####_PM.var(pm)[:Zff] = JuMP.@variable(pm.model, [(z,w) in base_cont], base_name="zff", binary = true, start = 0 ) #just one variable, how to create multiple with an index
            for (z,w) in base_cont
            display(JuMP.@constraint(pm.model,  FFRReserves_max[b]  >= FFR_list[(z,w)] ))
            display(JuMP.@constraint(pm.model,  FFRReserves_max[b] <= FFR_list[(z,w)] + 100*(1-Zff[(z,w)])))
            end
            display(JuMP.@constraint(pm.model,  sum(Zff[(z,w)] for (z, w) in base_cont) == 1 ) )
        end

        Scale = 8760*5/Total_sample # 5 for no. of gap years between two time steps
        # multiperiod

		JuMP.@constraint(pm.model, Inv_cost == sum(sum(conv["cost"]*_PM.var(pm, 1, :conv_ne, i) for (i,conv) in _PM.nws(pm)[1][:convdc_ne]) for (n, nw_ref) in _PM.nws(pm)) +
		           sum(sum(branch["cost"]*_PM.var(pm, 1, :branchdc_ne, i) for (i,branch) in _PM.nws(pm)[1][:branchdc_ne]) for (n, nw_ref) in _PM.nws(pm)) )
		JuMP.@constraint(pm.model, Gen_cost == sum(sum(Scale/10^6*gen_cost[(b,i)] for (i,gen) in _PM.nws(pm)[b][:gen]) for b in base_nws) )
		JuMP.@constraint(pm.model, FFRReserves == Scale*sum(FFRReserves_max[b] for b in base_nws) )
		JuMP.@constraint(pm.model, FCRReserves == Scale*sum(FCRReserves_max[b] for b in base_nws) )

		# display(JuMP.@constraint(pm.model, FCRReserves ==  sum(FCR_cost[(c,2)]*100*Scale/10^6*_PM.var(pm, c, :Pgg, 2)  for (b,c,br) in cont_nws) ) )
		display(JuMP.@constraint(pm.model, Cont == sum(_PM.ref(pm, b, :branchdc_ne, br)["fail_prob"]*_PM.var(pm, b, :branchdc_ne, br)*
            	(Scale/10^6*gen_cost[(c,3)] - Scale/10^6*gen_cost[(b,3)]) for (b,c,br) in cont_nws) ) )

      curtailment = Dict()
      capacity = Dict()

      for b in base_nws
          for c in curt_gen
            curtailment[(b,c)] = JuMP.upper_bound(_PM.var(pm, b, :pg, c)) - _PM.var(pm, b, :pg, c)
            capacity[(b,c)] = JuMP.upper_bound(_PM.var(pm, b, :pg, c))
         end
     end
     for y = 1:total_year
            JuMP.@constraint(pm.model, Curt[y] == sum(sum(curtailment[(b,c)] for c in curt_gen) for b in year_base[y]) / sum(sum(capacity[(b,c)] for c in curt_gen) for b in year_base[y])  )
            JuMP.@constraint(pm.model, Curt[y] <= max_curt)
      end
        return  JuMP.@objective(pm.model, Min, Inv_cost + Gen_cost + FFRReserves + FCRReserves + Cont - Gen_cost)
end

function objective_min_cost_TNEP_FSNS(pm::_PM.AbstractPowerModel)
	display("objective_min_cost_TNEP_FSNS")
	base_nws = pm.setting["base_list"]
	cont_nws = pm.setting["Cont_list"]
	Total_sample = pm.setting["Total_sample"]
	curt_gen = pm.setting["curtailed_gen"]
	max_curt = pm.setting["max_curt"]
	year_base = pm.setting["year_base"]
	total_year = pm.setting["total_yr"]

	gen_cost = Dict()
	FFR_cost = Dict()
	FCR_cost = Dict()
	fail_prob = Dict()
        Inv_cost = _PM.var(pm)[:Inv_cost] = JuMP.@variable(pm.model, start = 0, lower_bound = 0)
        Gen_cost = _PM.var(pm)[:Gen_cost] = JuMP.@variable(pm.model, start = 0, lower_bound = 0)
        FFRReserves = _PM.var(pm)[:FFR_Reserves] = JuMP.@variable(pm.model, start = 0, lower_bound = 0)
        FCRReserves = _PM.var(pm)[:FCR_Reserves] = JuMP.@variable(pm.model, start = 0, lower_bound = 0)
        Curt = _PM.var(pm)[:curt] = JuMP.@variable(pm.model, start = 0)
        Cont = _PM.var(pm)[:Cont]     = JuMP.@variable(pm.model, start = 0, lower_bound = 0)
        weights = pm.setting["weights"]

        sol_component_value_mod_wonw(pm,   :Inv_cost, Inv_cost)
        sol_component_value_mod_wonw(pm,   :Gen_cost, Gen_cost)
		sol_component_value_mod_wonw(pm,   :FFR_Reserves, FFRReserves)
        sol_component_value_mod_wonw(pm,   :FCR_Reserves, FCRReserves)
        sol_component_value_mod_wonw(pm,   :Cont, Cont)
		sol_component_value_mod_wonw(pm,   :Curt, Curt)

        for (n, nw_ref) in _PM.nws(pm)
            for (r, reserves) in nw_ref[:reserves]
                FFR_cost[(n,r)] = reserves["Cf"]
                FCR_cost[(n,r)] = reserves["Cg"]
            end

             for (i,gen) in nw_ref[:gen]
                pg = _PM.var(pm, n, :pg, i)

                if length(gen["cost"]) == 1
                    gen_cost[(n,i)] = gen["cost"][1]
                elseif length(gen["cost"]) == 2
                    gen_cost[(n,i)] = gen["cost"][1]*pg + gen["cost"][2]
                elseif length(gen["cost"]) == 3
                    gen_cost[(n,i)] = gen["cost"][2]*pg + gen["cost"][3]
                else
                    gen_cost[(n,i)] = 0.0
                end
            end
        end

        FFR_list =  Dict()
        FCR_list = Dict()
        for (b,c,br) in cont_nws
            FFR_list[(b,c)] = FFR_cost[(c,2)]*100/10^6*_PM.var(pm, c, :Pff, 2)+
                              FCR_cost[(c,2)]*100/10^6*_PM.var(pm, c, :Pgg, 2)
            # a_ub, a_lb = _IM.variable_domain(_PM.var(pm, c, :Pff, 2))
            # display("a_ub:$a_ub a_lb:$a_lb")
        end
        Zff = _PM.var(pm)[:Zff] = Dict((b, c) => JuMP.@variable(pm.model, base_name = "zff[$(string(b)),$(string(c))]", binary = true, start = 0) for (b, c, br) in cont_nws)
        sol_component_value_mod_wonw(pm, :Zff, Zff)
        # _IM.sol_component_value(pm, nw, :reserves, :Zff, base_nws, Zff)
        FFRReserves_max = _PM.var(pm)[:FFRReserves_max] = JuMP.@variable(pm.model, [b in base_nws], base_name ="FFRReserves_max", start = 0)
        FCRReserves_max = _PM.var(pm)[:FCRReserves_max] = JuMP.@variable(pm.model, [b in base_nws], base_name ="FCRReserves_max", start = 0)
        for b in base_nws
            sol_component_value_mod(pm, b, :FFRReserves_max, FFRReserves_max[b])
            sol_component_value_mod(pm, b, :FCRReserves_max, FCRReserves_max[b])
        end
        base_cont= Dict()
        for b in base_nws
            # [item for item in a if item[2] == 4]
            base_cont = [tt for tt in cont_nws if tt[1] == b]

            ####_PM.var(pm)[:Zff] = JuMP.@variable(pm.model, [(z,w) in base_cont], base_name="zff", binary = true, start = 0 ) #just one variable, how to create multiple with an index
            for (z,w) in base_cont
            JuMP.@constraint(pm.model,  FFRReserves_max[b] >= FFR_list[(z,w)] )
            JuMP.@constraint(pm.model,  FFRReserves_max[b] <= FFR_list[(z,w)] + 100*(1-Zff[(z,w)]))
            end
            JuMP.@constraint(pm.model,  sum(Zff[(z,w)] for (z, w) in base_cont) == 1 )
        end

        Scale = 8760*10  # 5 for no. of gap years between two time steps

            # for (i,branch) in _PM.nws(pm)[1][:branchdc_ne]
            #     display(branch["cost"])
            #  end
            #
            #  for (i,conv) in _PM.nws(pm)[1][:convdc_ne]
            #      display(conv["cost"])
            #   end
            # display(gen_cost)

            JuMP.@constraint(pm.model, Inv_cost == sum(sum(conv["cost"]*_PM.var(pm, 1, :conv_ne, i) for (i,conv) in _PM.nws(pm)[1][:convdc_ne]) for (n, nw_ref) in _PM.nws(pm)) +
                       sum(sum(branch["cost"]*_PM.var(pm, 1, :branchdc_ne, i) for (i,branch) in _PM.nws(pm)[1][:branchdc_ne]) for (n, nw_ref) in _PM.nws(pm)) )
			JuMP.@constraint(pm.model, Gen_cost == sum(weights[b]*sum(Scale/10^6*gen_cost[(b,i)] for (i,gen) in _PM.nws(pm)[b][:gen]) for b in base_nws) )
			# display(JuMP.@constraint(pm.model, FFRReserves == sum(FFR_cost[(c,2)]*100*Scale/10^6* _PM.var(pm, c, :Pff, 2) for (b,c,br) in cont_nws) ))
            JuMP.@constraint(pm.model, FFRReserves == Scale*sum(weights[b]*FFRReserves_max[b] for b in base_nws) )
            JuMP.@constraint(pm.model, FCRReserves == Scale*sum(weights[b]*FCRReserves_max[b] for b in base_nws) )

			# display(JuMP.@constraint(pm.model, FCRReserves ==  sum(FCR_cost[(c,2)]*100*Scale/10^6*_PM.var(pm, c, :Pgg, 2)  for (b,c,br) in cont_nws) ) )
			# print(JuMP.@constraint(pm.model, Cont == sum( sum( (Scale/10^6) *gen_cost[(c,i)]*_PM.ref(pm, b, :branchdc_ne, br)["fail_prob"]*_PM.var(pm, 1, :branchdc_ne, br)
            #         for i in  curt_gen)
            #         for (b,c,br) in cont_nws) ) )

            # display(JuMP.@constraint(pm.model, Cont == sum(weights[b]*_PM.ref(pm, b, :branchdc_ne, br)["fail_prob"]*_PM.var(pm, b, :branchdc_ne, br)*
            # (Scale/10^6*gen_cost[(c,3)] - Scale/10^6*gen_cost[(b,3)]) for (b,c,br) in cont_nws) ) )
 			curtailment = Dict()
			capacity = Dict()
			#
			#  for b in base_nws
			# 	  for c in curt_gen
			# 		curtailment[(b,c)] = JuMP.upper_bound(_PM.var(pm, b, :pg, c)) - _PM.var(pm, b, :pg, c)
            #         # curtailment[(b,c)] = _PM.var(pm, b, :pg, c)
			# 		capacity[(b,c)] = JuMP.upper_bound(_PM.var(pm, b, :pg, c))
			# 	 end
			#  end
			 curtailment = Dict()
	         capacity = Dict()

	         for b in base_nws
	             for c in curt_gen
	               curtailment[(b,c)] = JuMP.upper_bound(_PM.var(pm, b, :pg, c)) - _PM.var(pm, b, :pg, c)
	               capacity[(b,c)] = JuMP.upper_bound(_PM.var(pm, b, :pg, c))
	            end
	        end
	        for y = 1:total_year
	               JuMP.@constraint(pm.model, Curt[y] == sum(sum(curtailment[(b,c)] for c in curt_gen) for b in year_base[y]) / sum(sum(capacity[(b,c)] for c in curt_gen) for b in year_base[y])  )
	               # JuMP.@constraint(pm.model, Curt[y] <= max_curt)
	        end
            return JuMP.@objective(pm.model, Min, Inv_cost + FFRReserves + FCRReserves + Gen_cost)
end

function objective_min_cost_TNEP_PL_nocl(pm::_PM.AbstractPowerModel)
	base_nws = pm.setting["base_list"]
	cont_nws = pm.setting["Cont_list"]
	Total_sample = pm.setting["Total_sample"]
	curt_gen = pm.setting["curtailed_gen"]
	max_curt = pm.setting["max_curt"]
	year_base = pm.setting["year_base"]
	total_year = pm.setting["total_yr"]
	weights = pm.setting["weights"]

	gen_cost = Dict()
	FFR_cost = Dict()
	FCR_cost = Dict()
	fail_prob = Dict()
	Inv_cost = _PM.var(pm)[:Inv_cost] = JuMP.@variable(pm.model, start = 0, lower_bound = 0)
	Gen_cost = _PM.var(pm)[:Gen_cost] = JuMP.@variable(pm.model, start = 0, lower_bound = 0)
	FFRReserves = _PM.var(pm)[:FFR_Reserves] = JuMP.@variable(pm.model, start = 0, lower_bound = 0)
	FCRReserves = _PM.var(pm)[:FCR_Reserves] = JuMP.@variable(pm.model, start = 0, lower_bound = 0)
	Cont = _PM.var(pm)[:Cont]     = JuMP.@variable(pm.model, start = 0, lower_bound = 0)
	Curt = _PM.var(pm)[:curt] = JuMP.@variable(pm.model, [i in 1:total_year], start = 0, lower_bound = 0)

	sol_component_value_mod_wonw(pm,   :Inv_cost, Inv_cost)
	sol_component_value_mod_wonw(pm,   :Gen_cost, Gen_cost)
	sol_component_value_mod_wonw(pm,   :Curt, Curt)
	sol_component_value_mod_wonw(pm,   :FFR_Reserves, FFRReserves)
	sol_component_value_mod_wonw(pm,   :FCR_Reserves, FCRReserves)
	sol_component_value_mod_wonw(pm,   :Cont, Cont)

	for (n, nw_ref) in _PM.nws(pm)
		for (r, reserves) in nw_ref[:reserves]
			FFR_cost[(n,r)] = reserves["Cf"]
			FCR_cost[(n,r)] = reserves["Cg"]
		end

		 for (i,gen) in nw_ref[:gen]
			pg = _PM.var(pm, n, :pg, i)

			if length(gen["cost"]) == 1
				gen_cost[(n,i)] = gen["cost"][1]
			elseif length(gen["cost"]) == 2
				gen_cost[(n,i)] = gen["cost"][1]*pg + gen["cost"][2]
			elseif length(gen["cost"]) == 3
				gen_cost[(n,i)] = gen["cost"][2]*pg + gen["cost"][3]
			else
				gen_cost[(n,i)] = 0.0
			end
		end
	end

	Scale = 8760/Total_sample*5 # 5 for no. of gap years between two time steps

	JuMP.@constraint(pm.model, Inv_cost == sum(sum(conv["cost"]*_PM.var(pm, 1, :conv_ne, i) for (i,conv) in _PM.nws(pm)[1][:convdc_ne]) for (n, nw_ref) in _PM.nws(pm)) +
			   sum(sum(branch["cost"]*_PM.var(pm, 1, :branchdc_ne, i) for (i,branch) in _PM.nws(pm)[1][:branchdc_ne]) for (n, nw_ref) in _PM.nws(pm)) )
	JuMP.@constraint(pm.model, Gen_cost == sum(sum(Scale/10^6*gen_cost[(b,i)] for (i,gen) in _PM.nws(pm)[b][:gen]) for b in base_nws) )
	JuMP.@constraint(pm.model, FFRReserves == sum(FFR_cost[(c,2)]*100*Scale/10^6* _PM.var(pm, c, :Pff, 2) for (b,c,br) in cont_nws) )
	JuMP.@constraint(pm.model, FCRReserves ==  sum(FCR_cost[(c,2)]*100*Scale/10^6*_PM.var(pm, c, :Pgg, 2)  for (b,c,br) in cont_nws) )
	JuMP.@constraint(pm.model, Cont == sum(_PM.ref(pm, b, :branchdc_ne, br)["fail_prob"]*_PM.var(pm, b, :branchdc_ne, br)*
	(Scale/10^6*gen_cost[(c,3)] - Scale/10^6*gen_cost[(b,3)]) for (b,c,br) in cont_nws) )
	curtailment = Dict()
	capacity = Dict()

	  for b in base_nws
	      for c in curt_gen
	        curtailment[(b,c)] = JuMP.upper_bound(_PM.var(pm, b, :pg, c)) - _PM.var(pm, b, :pg, c)
	        capacity[(b,c)] = JuMP.upper_bound(_PM.var(pm, b, :pg, c))
	     end
	 end
	 for y = 1:total_year
	        JuMP.@constraint(pm.model, Curt[y] == sum(sum(curtailment[(b,c)] for c in curt_gen) for b in year_base[y]) / sum(sum(capacity[(b,c)] for c in curt_gen) for b in year_base[y])  )
	        # JuMP.@constraint(pm.model, Curt[y] <= max_curt)
	  end

	return display( JuMP.@objective(pm.model, Min, Inv_cost + Gen_cost + FFRReserves + FCRReserves + Cont  ) )
end

function objective_min_cost_TNEP_PL(pm::_PM.AbstractPowerModel)
	base_nws = pm.setting["base_list"]
	cont_nws = pm.setting["Cont_list"]
	Total_sample = pm.setting["Total_sample"]
	curt_gen = pm.setting["curtailed_gen"]
	max_curt = pm.setting["max_curt"]
	year_base = pm.setting["year_base"]
	total_year = pm.setting["total_yr"]

	gen_cost = Dict()
	FFR_cost = Dict()
	FCR_cost = Dict()
	fail_prob = Dict()
        Inv_cost = _PM.var(pm)[:Inv_cost] = JuMP.@variable(pm.model, start = 0, lower_bound = 0)
        Gen_cost = _PM.var(pm)[:Gen_cost] = JuMP.@variable(pm.model, start = 0, lower_bound = 0)
        FFRReserves = _PM.var(pm)[:FFR_Reserves] = JuMP.@variable(pm.model, start = 0, lower_bound = 0)
        FCRReserves = _PM.var(pm)[:FCR_Reserves] = JuMP.@variable(pm.model, start = 0, lower_bound = 0)
        Curt = _PM.var(pm)[:curt] = JuMP.@variable(pm.model, start = 0)
        Cont = _PM.var(pm)[:Cont]     = JuMP.@variable(pm.model, start = 0, lower_bound = 0)
        weights = pm.setting["weights"]

        sol_component_value_mod_wonw(pm,   :Inv_cost, Inv_cost)
        sol_component_value_mod_wonw(pm,   :Gen_cost, Gen_cost)
		sol_component_value_mod_wonw(pm,   :FFR_Reserves, FFRReserves)
        sol_component_value_mod_wonw(pm,   :FCR_Reserves, FCRReserves)
        sol_component_value_mod_wonw(pm,   :Cont, Cont)
		sol_component_value_mod_wonw(pm,   :Curt, Curt)

        for (n, nw_ref) in _PM.nws(pm)
            for (r, reserves) in nw_ref[:reserves]
                FFR_cost[(n,r)] = reserves["Cf"]
                FCR_cost[(n,r)] = reserves["Cg"]
            end

             for (i,gen) in nw_ref[:gen]
                pg = _PM.var(pm, n, :pg, i)

                if length(gen["cost"]) == 1
                    gen_cost[(n,i)] = gen["cost"][1]
                elseif length(gen["cost"]) == 2
                    gen_cost[(n,i)] = gen["cost"][1]*pg + gen["cost"][2]
                elseif length(gen["cost"]) == 3
                    gen_cost[(n,i)] = gen["cost"][2]*pg + gen["cost"][3]
                else
                    gen_cost[(n,i)] = 0.0
                end
            end
        end

        FFR_list =  Dict()
        FCR_list = Dict()
        for (b,c,br) in cont_nws
            FFR_list[(b,c)] = FFR_cost[(c,2)]*100/10^6*_PM.var(pm, c, :Pff, 2)+
                              FCR_cost[(c,2)]*100/10^6*_PM.var(pm, c, :Pgg, 2)
            # a_ub, a_lb = _IM.variable_domain(_PM.var(pm, c, :Pff, 2))
            # display("a_ub:$a_ub a_lb:$a_lb")
        end
        Zff = _PM.var(pm)[:Zff] = Dict((b, c) => JuMP.@variable(pm.model, base_name = "zff[$(string(b)),$(string(c))]", binary = true, start = 0) for (b, c, br) in cont_nws)
        sol_component_value_mod_wonw(pm, :Zff, Zff)
        # _IM.sol_component_value(pm, nw, :reserves, :Zff, base_nws, Zff)
        FFRReserves_max = _PM.var(pm)[:FFRReserves_max] = JuMP.@variable(pm.model, [b in base_nws], base_name ="FFRReserves_max", start = 0)
        FCRReserves_max = _PM.var(pm)[:FCRReserves_max] = JuMP.@variable(pm.model, [b in base_nws], base_name ="FCRReserves_max", start = 0)
        for b in base_nws
            sol_component_value_mod(pm, b, :FFRReserves_max, FFRReserves_max[b])
            sol_component_value_mod(pm, b, :FCRReserves_max, FCRReserves_max[b])
        end
        base_cont= Dict()
        for b in base_nws
            # [item for item in a if item[2] == 4]
            base_cont = [tt for tt in cont_nws if tt[1] == b]

            ####_PM.var(pm)[:Zff] = JuMP.@variable(pm.model, [(z,w) in base_cont], base_name="zff", binary = true, start = 0 ) #just one variable, how to create multiple with an index
            for (z,w) in base_cont
            JuMP.@constraint(pm.model,  FFRReserves_max[b] >= FFR_list[(z,w)] )
            JuMP.@constraint(pm.model,  FFRReserves_max[b] <= FFR_list[(z,w)] + 100*(1-Zff[(z,w)]))
            end
            JuMP.@constraint(pm.model,  sum(Zff[(z,w)] for (z, w) in base_cont) == 1 )
        end

        Scale = 8760*10  # 5 for no. of gap years between two time steps

            # for (i,branch) in _PM.nws(pm)[1][:branchdc_ne]
            #     display(branch["cost"])
            #  end
            #
            #  for (i,conv) in _PM.nws(pm)[1][:convdc_ne]
            #      display(conv["cost"])
            #   end
            # display(gen_cost)

            JuMP.@constraint(pm.model, Inv_cost == sum(sum(conv["cost"]*_PM.var(pm, 1, :conv_ne, i) for (i,conv) in _PM.nws(pm)[1][:convdc_ne]) for (n, nw_ref) in _PM.nws(pm)) +
                       sum(sum(branch["cost"]*_PM.var(pm, 1, :branchdc_ne, i) for (i,branch) in _PM.nws(pm)[1][:branchdc_ne]) for (n, nw_ref) in _PM.nws(pm)) )
			JuMP.@constraint(pm.model, Gen_cost == sum(weights[b]*sum(Scale/10^6*gen_cost[(b,i)] for (i,gen) in _PM.nws(pm)[b][:gen]) for b in base_nws) )
			# display(JuMP.@constraint(pm.model, FFRReserves == sum(FFR_cost[(c,2)]*100*Scale/10^6* _PM.var(pm, c, :Pff, 2) for (b,c,br) in cont_nws) ))
            JuMP.@constraint(pm.model, FFRReserves == Scale*sum(weights[b]*FFRReserves_max[b] for b in base_nws) )
            JuMP.@constraint(pm.model, FCRReserves == Scale*sum(weights[b]*FCRReserves_max[b] for b in base_nws) )

			# display(JuMP.@constraint(pm.model, FCRReserves ==  sum(FCR_cost[(c,2)]*100*Scale/10^6*_PM.var(pm, c, :Pgg, 2)  for (b,c,br) in cont_nws) ) )
			# print(JuMP.@constraint(pm.model, Cont == sum( sum( (Scale/10^6) *gen_cost[(c,i)]*_PM.ref(pm, b, :branchdc_ne, br)["fail_prob"]*_PM.var(pm, 1, :branchdc_ne, br)
            #         for i in  curt_gen)
            #         for (b,c,br) in cont_nws) ) )

            # display(JuMP.@constraint(pm.model, Cont == sum(weights[b]*_PM.ref(pm, b, :branchdc_ne, br)["fail_prob"]*_PM.var(pm, b, :branchdc_ne, br)*
            # (Scale/10^6*gen_cost[(c,3)] - Scale/10^6*gen_cost[(b,3)]) for (b,c,br) in cont_nws) ) )``
 			curtailment = Dict()
			capacity = Dict()
			#
			#  for b in base_nws
			# 	  for c in curt_gen
			# 		curtailment[(b,c)] = JuMP.upper_bound(_PM.var(pm, b, :pg, c)) - _PM.var(pm, b, :pg, c)
            #         # curtailment[(b,c)] = _PM.var(pm, b, :pg, c)
			# 		capacity[(b,c)] = JuMP.upper_bound(_PM.var(pm, b, :pg, c))
			# 	 end
			#  end
			 curtailment = Dict()
	         capacity = Dict()

	         for b in base_nws
	             for c in curt_gen
	               curtailment[(b,c)] = JuMP.upper_bound(_PM.var(pm, b, :pg, c)) - _PM.var(pm, b, :pg, c)
	               capacity[(b,c)] = JuMP.upper_bound(_PM.var(pm, b, :pg, c))
	            end
	        end
	        for y = 1:total_year
	               JuMP.@constraint(pm.model, Curt[y] == sum(sum(curtailment[(b,c)] for c in curt_gen) for b in year_base[y]) / sum(sum(capacity[(b,c)] for c in curt_gen) for b in year_base[y])  )
	               # JuMP.@constraint(pm.model, Curt[y] <= max_curt)
	        end

	return display( JuMP.@objective(pm.model, Min, Inv_cost + Gen_cost + FFRReserves + FCRReserves) )
end

function objective_min_cost_TNEP_FSNS_rev1(pm::_PM.AbstractPowerModel)
	base_nws = pm.setting["base_list"]
	cont_nws = pm.setting["Cont_list"]
	Total_sample = pm.setting["Total_sample"]
	curt_gen = pm.setting["curtailed_gen"]
	max_curt = pm.setting["max_curt"]
	year_base = pm.setting["year_base"]
	total_year = pm.setting["total_yr"]

	gen_cost = Dict()
	FFR_cost = Dict()
	FCR_cost = Dict()
	fail_prob = Dict()
        Inv_cost = _PM.var(pm)[:Inv_cost] = JuMP.@variable(pm.model, start = 0, lower_bound = 0)
        Gen_cost = _PM.var(pm)[:Gen_cost] = JuMP.@variable(pm.model, start = 0, lower_bound = 0)
        FFRReserves = _PM.var(pm)[:FFR_Reserves] = JuMP.@variable(pm.model, start = 0, lower_bound = 0)
        FCRReserves = _PM.var(pm)[:FCR_Reserves] = JuMP.@variable(pm.model, start = 0, lower_bound = 0)
        Curt = _PM.var(pm)[:curt] = JuMP.@variable(pm.model, start = 0)
        Cont = _PM.var(pm)[:Cont]     = JuMP.@variable(pm.model, start = 0, lower_bound = 0)
        weights = pm.setting["weights"]

        sol_component_value_mod_wonw(pm,   :Inv_cost, Inv_cost)
        sol_component_value_mod_wonw(pm,   :Gen_cost, Gen_cost)
		sol_component_value_mod_wonw(pm,   :FFR_Reserves, FFRReserves)
        sol_component_value_mod_wonw(pm,   :FCR_Reserves, FCRReserves)
        sol_component_value_mod_wonw(pm,   :Cont, Cont)

        for (n, nw_ref) in _PM.nws(pm)
            for (r, reserves) in nw_ref[:reserves]
                FFR_cost[(n,r)] = reserves["Cf"]
                FCR_cost[(n,r)] = reserves["Cg"]
            end

             for (i,gen) in nw_ref[:gen]
                pg = _PM.var(pm, n, :pg, i)

                if length(gen["cost"]) == 1
                    gen_cost[(n,i)] = gen["cost"][1]
                elseif length(gen["cost"]) == 2
                    gen_cost[(n,i)] = gen["cost"][1]*pg + gen["cost"][2]
                elseif length(gen["cost"]) == 3
                    gen_cost[(n,i)] = gen["cost"][2]*pg + gen["cost"][3]
                else
                    gen_cost[(n,i)] = 0.0
                end
            end
        end

        FFR_list =  Dict()
        FCR_list = Dict()
        for (b,c,br) in cont_nws
            FFR_list[(b,c)] = FFR_cost[(c,2)]*100/10^6*_PM.var(pm, c, :Pff, 2)+
                              FCR_cost[(c,2)]*100/10^6*_PM.var(pm, c, :Pgg, 2)
            # a_ub, a_lb = _IM.variable_domain(_PM.var(pm, c, :Pff, 2))
            # display("a_ub:$a_ub a_lb:$a_lb")
        end
        Zff = _PM.var(pm)[:Zff] = Dict((b, c) => JuMP.@variable(pm.model, base_name = "zff[$(string(b)),$(string(c))]", binary = true, start = 0) for (b, c, br) in cont_nws)
        sol_component_value_mod_wonw(pm, :Zff, Zff)
        # _IM.sol_component_value(pm, nw, :reserves, :Zff, base_nws, Zff)
        FFRReserves_max = _PM.var(pm)[:FFRReserves_max] = JuMP.@variable(pm.model, [b in base_nws], base_name ="FFRReserves_max", start = 0)
        FCRReserves_max = _PM.var(pm)[:FCRReserves_max] = JuMP.@variable(pm.model, [b in base_nws], base_name ="FCRReserves_max", start = 0)
        for b in base_nws
            sol_component_value_mod(pm, b, :FFRReserves_max, FFRReserves_max[b])
            sol_component_value_mod(pm, b, :FCRReserves_max, FCRReserves_max[b])
        end
        base_cont= Dict()
        for b in base_nws
            # [item for item in a if item[2] == 4]
            base_cont = [tt for tt in cont_nws if tt[1] == b]

            ####_PM.var(pm)[:Zff] = JuMP.@variable(pm.model, [(z,w) in base_cont], base_name="zff", binary = true, start = 0 ) #just one variable, how to create multiple with an index
            for (z,w) in base_cont
            JuMP.@constraint(pm.model,  FFRReserves_max[b]  >= FFR_list[(z,w)] )
            JuMP.@constraint(pm.model,  FFRReserves_max[b] <= FFR_list[(z,w)] + 100*(1-Zff[(z,w)]))
            end
            JuMP.@constraint(pm.model,  sum(Zff[(z,w)] for (z, w) in base_cont) == 1 )
        end

        Scale = 8760*10  # 5 for no. of gap years between two time steps

            JuMP.@constraint(pm.model, Inv_cost == sum(sum(conv["cost"]*_PM.var(pm, 1, :conv_ne, i) for (i,conv) in _PM.nws(pm)[1][:convdc_ne]) for (n, nw_ref) in _PM.nws(pm)) +
                       sum(sum(branch["cost"]*_PM.var(pm, 1, :branchdc_ne, i) for (i,branch) in _PM.nws(pm)[1][:branchdc_ne]) for (n, nw_ref) in _PM.nws(pm)) )
			display(JuMP.@constraint(pm.model, Gen_cost == sum(weights[b]*sum(Scale/10^6*gen_cost[(b,i)] for i in curt_gen) for b in base_nws) ) )
			# display(JuMP.@constraint(pm.model, FFRReserves == sum(FFR_cost[(c,2)]*100*Scale/10^6* _PM.var(pm, c, :Pff, 2) for (b,c,br) in cont_nws) ))
            JuMP.@constraint(pm.model, FFRReserves == Scale*sum(weights[b]*FFRReserves_max[b] for b in base_nws) )
            JuMP.@constraint(pm.model, FCRReserves == Scale*sum(weights[b]*FCRReserves_max[b] for b in base_nws) )

            # JuMP.@constraint(pm.model, Cont == sum(weights[b]*_PM.ref(pm, b, :branchdc_ne, br)["fail_prob"]*_PM.var(pm, b, :branchdc_ne, br)*
            # (Scale/10^6*(gen_cost[(c,1)] + gen_cost[(c,2)]) - Scale/10^6*(gen_cost[(b,1)]+gen_cost[(b,2)]) ]) for (b,c,br) in cont_nws) )

			curtailment = Dict()
			capacity = Dict()

			 # for b in base_nws
				#   for c in curt_gen
				# 	curtailment[(b,c)] = JuMP.upper_bound(_PM.var(pm, b, :pg, c)) - _PM.var(pm, b, :pg, c)
             #        # curtailment[(b,c)] = _PM.var(pm, b, :pg, c)
				# 	capacity[(b,c)] = JuMP.upper_bound(_PM.var(pm, b, :pg, c))
				#  end
			 # end
			 # curtailment = Dict()
	         # capacity = Dict()

	        #  for b in base_nws
	        #      for c in curt_gen
	        #        curtailment[(b,c)] = JuMP.upper_bound(_PM.var(pm, b, :pg, c)) - _PM.var(pm, b, :pg, c)
	        #        capacity[(b,c)] = JuMP.upper_bound(_PM.var(pm, b, :pg, c))
	        #     end
	        # end
	        # for y = 1:total_year
	        #        JuMP.@constraint(pm.model, Curt[y] == sum(sum(curtailment[(b,c)] for c in curt_gen) for b in year_base[y]) / sum(sum(capacity[(b,c)] for c in curt_gen) for b in year_base[y])  )
	        #        JuMP.@constraint(pm.model, Curt[y] <= max_curt)
	        #  end
			 display("objective_min_cost_TNEP_FSNS_rev1")
            return JuMP.@objective(pm.model, Min,  FFRReserves + FCRReserves + Inv_cost  - Gen_cost)
end


text:
export DynamicObstacles

using Random

mutable struct DynamicObstacles{R} <: AbstractGridWorld
    world::GridWorldBase{Tuple{Empty,Wall,Obstacle,Goal}}
    agent_pos::CartesianIndex{2}
    agent::Agent
    num_obstacles::Int
    obstacle_pos::Array{CartesianIndex{2},1}
    obstacle_reward::Float64
    goal_reward::Float64
    reward::Float64
    rng::R
end

function DynamicObstacles(;n = 8, agent_start_pos = CartesianIndex(2,2), agent_start_dir = RIGHT, goal_pos = CartesianIndex(n-1, n-1), num_obstacles = n-3, rng = Random.GLOBAL_RNG)
    objects = (EMPTY, WALL, OBSTACLE, GOAL)
    w = GridWorldBase(objects, n, n)

    w[WALL, [1,n], 1:n] .= true
    w[WALL, 1:n, [1,n]] .= true

    obstacle_pos = Array{CartesianIndex{2},1}(undef,num_obstacles)

    obstacle_reward = -1.0
    goal_reward = 1.0
    reward = 0.0

    env = DynamicObstacles(w, agent_start_pos, Agent(dir=agent_start_dir), num_obstacles, obstacle_pos, obstacle_reward, goal_reward, reward, rng)

    reset!(env, agent_start_pos = agent_start_pos, agent_start_dir = agent_start_dir, goal_pos = goal_pos)

    return env
end

iscollision(env::DynamicObstacles) = env.world[OBSTACLE, env.agent_pos]

function (env::DynamicObstacles)(::MoveForward)
    env.reward = 0.0
    update_obstacles!(env)
    dir = get_dir(env.agent)
    dest = dir(env.agent_pos)
    if !env.world[WALL, dest]
        env.agent_pos = dest
    end
    if iscollision(env)
        env.reward = env.obstacle_reward
    end
    return env
end

function (env::DynamicObstacles)(action::Union{TurnRight, TurnLeft})
    env.reward = 0.0
    update_obstacles!(env)
    agent = get_agent(env)
    set_dir!(agent, action(get_dir(agent)))
    if iscollision(env)
        env.reward = env.obstacle_reward
    end
    return env
end

function valid_obstacle_dest(env::DynamicObstacles, pos::CartesianIndex{2})
    candidate_pos = [CartesianIndex(pos[1] + i, pos[2] + j) for i in [-1, 0, 1] for j in [-1, 0, 1]]
    filter(p -> env.world[EMPTY, p] || pos == p, candidate_pos)
end

function update_obstacles!(env::DynamicObstacles)
    for (i, pos) in enumerate(env.obstacle_pos)
        env.world[EMPTY, pos] = true
        env.world[OBSTACLE, pos] = false

        new_pos = rand(env.rng, valid_obstacle_dest(env, pos))
        env.obstacle_pos[i] = new_pos

        env.world[EMPTY, new_pos] = false
        env.world[OBSTACLE, new_pos] = true
    end
    return env
end

RLBase.get_terminal(env::DynamicObstacles) = iscollision(env) || env.world[GOAL, env.agent_pos]

function RLBase.reset!(env::DynamicObstacles; agent_start_pos = CartesianIndex(2, 2), agent_start_dir = RIGHT, goal_pos = CartesianIndex(size(env.world)[end] - 1, size(env.world)[end] - 1))

    n = size(env.world)[end]
    env.world[EMPTY, 2:n-1, 2:n-1] .= true
    env.world[GOAL, goal_pos] = true
    env.world[EMPTY, goal_pos] = false
    env.agent_pos = agent_start_pos
    agent = get_agent(env)
    set_dir!(agent, agent_start_dir)

    obstacles_placed = 0
    while obstacles_placed < env.num_obstacles
        pos = CartesianIndex(rand(env.rng, 2:n-1), rand(env.rng, 2:n-1))
        if (pos == env.agent_pos) || (env.world[OBSTACLE, pos] == true) || (pos == goal_pos)
            continue
        else
            env.world[OBSTACLE, pos] = true
            env.world[EMPTY, pos] = false
            obstacles_placed = obstacles_placed + 1
            env.obstacle_pos[obstacles_placed] = pos
        end
    end
    return env
end