text:
abstract type Model end

struct SimpleModel{T} <: Model
    class::T
    args::Tuple
    kwargs::Tuple
end

Model(class, args...; kwargs...) = SimpleModel(class, tuple(args...), tuple(kwargs...))

struct FormulaModel{T} <: Model
    class::T
    formula::FormulaTerm
    args::Tuple
    kwargs::Tuple
end

Model(class, formula::FormulaTerm, args...; kwargs...) = FormulaModel(class, formula, tuple(args...), tuple(kwargs...))

abstract type BootstrapSampling end

abstract type ParametricSampling <: BootstrapSampling end
abstract type NonParametricSampling <: BootstrapSampling end

"""
Basic Sampling

```julia
BasicSampling(1000)
```

"""
struct BasicSampling <: BootstrapSampling
    nrun::Int
end

"""
Antithetic Sampling

```julia
AntitheticSampling(1000)
```

"""
struct AntitheticSampling <: BootstrapSampling
    nrun::Int
end

"""
Balanced Sampling

```julia
BalancedSampling(1000)
```

"""
struct BalancedSampling <: BootstrapSampling
    nrun::Int
end

"""
Residual Sampling

```julia
ResidualSampling(1000)
```

"""
struct ResidualSampling <: BootstrapSampling
    nrun::Int
end

"""
Wild Sampling

```julia
WildSampling(1000, rademacher)
WildSampling(1000, mammen)
```

"""
struct WildSampling <: BootstrapSampling
    nrun::Int
    noise::Function
end

"""
Exact Sampling

```julia
ExactSampling()
```

"""
mutable struct ExactSampling <: BootstrapSampling
    nrun::Int
end

ExactSampling() = ExactSampling(0)

"""
[Maximum Entropy Sampling](https://cran.r-project.org/web/packages/meboot/vignettes/meboot.pdf)

```julia
MaximumEntropySampling(100)
maximumEntropySampling(100, MaximumEntropyCache())
```

NOTE: Implementation based off [pymeboot](https://github.com/kirajcg/pymeboot) as the original
[R package](https://cran.r-project.org/web/packages/meboot/index.html) is GPL licensed.
"""
struct MaximumEntropySampling <: BootstrapSampling
    nrun::Int
    cache::MaximumEntropyCache
end

MaximumEntropySampling(nrun) = MaximumEntropySampling(nrun, MaximumEntropyCache())

abstract type BootstrapSample end

struct NonParametricBootstrapSample{T} <: BootstrapSample
    t0::Tuple
    t1::Tuple
    statistic::Function
    data::T
    sampling::BootstrapSampling
end

struct ParametricBootstrapSample{T,M} <: BootstrapSample
    t0::Tuple
    t1::Tuple
    statistic::Function
    data::T
    model::M
    sampling::BootstrapSampling
end

"""
Number of samples drawn from a bootstrap sampling

```julia
bs = BasicSampling(1000)
nrun(bs)

# output

1000
```

"""
nrun(bs::BootstrapSampling) = bs.nrun
nrun(bs::BootstrapSample) = nrun(sampling(bs))

"""
Return the statistic function of a `BootstrapSample`

```julia
bs = bootstrap(randn(20), mean, BasicSampling(100))
statistic(bs)
```

"""
statistic(bs::BootstrapSample) = bs.statistic

"""
Return the noise function of a wild bootstrap sampling
"""
noise(bs::WildSampling) = bs.noise

nvar(t::Tuple) = length(t)
nvar(bs::BootstrapSample) = nvar(original(bs))

tx(x) = tuple(x...)

zeros_tuple(t, n) = tuple([zeros(typeof(x), n) for x in t]...)

"""
bootstrap(statistic, data, BasicSampling())
"""
function bootstrap(statistic::Function, data, sampling::BasicSampling)
    t0 = tx(statistic(data))
    m = nrun(sampling)
    t1 = zeros_tuple(t0, m)
    data1 = copy(data)
    for i in 1:m
        draw!(data, data1)
        for (j, t) in enumerate(tx(statistic(data1)))
            t1[j][i] = t
        end
    end
    return NonParametricBootstrapSample(t0, t1, statistic, data, sampling)
end

"""
bootstrap(statistic, data, AntitheticSampling)
"""
function bootstrap(statistic::Function, data::AbstractVector, sampling::AntitheticSampling)
    t0 = tx(statistic(data))
    m = nrun(sampling)
    n = nobs(data)
    t1 = zeros_tuple(t0, m)
    idx = collect(1:n)
    idx1 = copy(idx)
    data1 = copy(data)
    data0 = copy(data)
    sort!(data0)
    for i in 1:m
        if isodd(i)
            sample!(idx, idx1)
        else
            idx1 = n .- idx1 .+ 1
        end
        data1 = pick(data0, idx1)
        for (j, t) in enumerate(tx(statistic(data1)))
            t1[j][i] = t
        end
    end
    return NonParametricBootstrapSample(t0, t1, statistic, data, sampling)
end

"""
bootstrap(statistic, data, sampling)
"""
function bootstrap(statistic::Function, data, sampling::BalancedSampling)
    n = nobs(data)
    m = nrun(sampling)
    t0 = tx(statistic(data))
    t1 = zeros_tuple(t0, m)
    idx = repeat([1:n;], m)
    ridx = zeros(Int, n, m)
    sample!(idx, ridx, replace = false)
    for i in 1:m
        for (j, t) in enumerate(tx(statistic(pick(data, ridx[:,i]))))
            t1[j][i] = t
        end
    end
    return NonParametricBootstrapSample(t0, t1, statistic, data, sampling)
end

struct ExactIterator{T}
    a::T
    k::Int64
end

exact(n::Int) = ExactIterator(1:n, n)

function iterate(itr::ExactIterator, s = ones(Int, itr.k))
    if length(s) > 0 && s[1] < 1
        return nothing
    end
    r = itr.a[s]
    for i = itr.k:-1:1
        if s[i] < lastindex(itr.a)
            s[i] = nextind(itr.a, s[i])
            for j = i:itr.k - 1
                s[j + 1] = s[j]
            end
            return (r, s)
        end
    end
    s[1] = 0
    return (r, s)
end

eltype(itr::ExactIterator) = typeof(itr.a)
eltype(itr::ExactIterator{UnitRange{T}}) where {T} = Array{T,1}
eltype(itr::ExactIterator{AbstractRange{T}}) where {T} = Array{T,1}

length(itr::ExactIterator) = binomial(length(itr.a) + itr.k - 1, itr.k)
size(itr::ExactIterator) = (itr.k,)

"""
bootstrap(statistic, data, sampling)
"""
function bootstrap(statistic::Function, data, sampling::ExactSampling)
    n = nobs(data)
    m = nrun_exact(n)
    t0 = tx(statistic(data))
    t1 = zeros_tuple(t0, m)
    for (i, idx) in enumerate(exact(n))
        for (j, t) in enumerate(tx(statistic(pick(data, idx))))
            t1[j][i] = t
        end
    end
    sampling.nrun = m
    return NonParametricBootstrapSample(t0, t1, statistic, data, sampling)
end

"""
bootstrap(statistic, data, MaximumEntropySampling)
"""
function bootstrap(statistic::Function, data, sampling::MaximumEntropySampling)
    init!(sampling.cache, data)

    t0 = tx(statistic(data))
    m = nrun(sampling)
    t1 = zeros_tuple(t0, m)
    data1 = copy(data)
    for i in 1:m
        draw!(sampling.cache, data, data1)
        for (j, t) in enumerate(tx(statistic(data1)))
            t1[j][i] = t
        end
    end
    return NonParametricBootstrapSample(t0, t1, statistic, data, sampling)
end

"""
bootstrap(statistic, data, model, sampling)
"""
function bootstrap(statistic::Function, data, model::SimpleModel, sampling::BootstrapSampling)
    f0 = fit(model.class, data, model.args...; model.kwargs...)
    t0 = tx(statistic(data))
    m = nrun(sampling)
    t1 = zeros_tuple(t0, m)
    data1 = copy(data)
    for i in 1:m
        draw!(f0, data1)
        for (j, t) in enumerate(tx(statistic(data1)))
            t1[j][i] = t
        end
    end
    return ParametricBootstrapSample(t0, t1, statistic, data, model, sampling)
end

struct ResidualTerm{S,T,R} <: StatsModels.AbstractTerm
    sampling::S
    t::T
    r0::R
    r1::R
    ResidualTerm(s::S, t::T, r::R) where {S,T,R} = new{S,T,R}(s, t, r, copy(r))
end

StatsModels.modelcols(t::ResidualTerm{S}, data) where {S} =
    throw(ArgumentError("ResidualTerm: don't know how to generate model columns for $S"))
StatsModels.modelcols(t::ResidualTerm{ResidualSampling}, data) =
    modelcols(t.t, data) + sample!(t.r0, t.r1)
StatsModels.modelcols(t::ResidualTerm{WildSampling}, data) =
    modelcols(t.t, data) + t.sampling.noise(t.r0)

# we only ever create one of these terms after the model has been fit once so we don't
# need to create a new schema for it every time...
StatsModels.needs_schema(t::ResidualTerm) = false
StatsModels.termvars(t::ResidualTerm) = StatsModels.termvars(t.t)

"""
bootstrap(statistic, data, model, formula, sampling)
bootstrap(statistic, data, model, formula, Wildsampling(nrun, noise))
"""
function bootstrap(statistic::Function, data::AbstractDataFrame, model::FormulaModel,
                   sampling::S) where {S <: BootstrapSampling}
    class = model.class
    formula = apply_schema(model.formula, schema(model.formula, data), model.class)

    args = model.args
    kwargs = model.kwargs

    f0 = fit(class, formula, data, args...; kwargs...)
    r0 = predict(f0) - response(f0)
    # replace LHS of formula with a special term that resamples from the residuals
    # whenever it's asked for model columns
    formula_resid = FormulaTerm(ResidualTerm(sampling, formula.lhs, r0), formula.rhs)

    m = nrun(sampling)
    t0 = tx(statistic(f0))
    t1 = zeros_tuple(t0, m)

    for i in 1:m
        f1 = fit(class, formula_resid, data, args...; kwargs...)
        for (j, t) in enumerate(tx(statistic(f1)))
            t1[j][i] = t
        end
    end
    return ParametricBootstrapSample(t0, t1, statistic, data, model, sampling)
end


text:
include(joinpath("caches", "link_transform.jl"))
include(joinpath("caches", "frame_point.jl"))

struct ControllerCache{n,S}
    link_caches::Vector{LinkFrameCache{n,S}}
    point_caches::Vector{FramePointCache{n,S}}
end
ControllerCache(link_caches::Vector{<:LinkFrameCache{n,S}}, point_caches) where {n,S} =
    ControllerCache{n,S}(link_caches, point_caches)

function setdirty!(cache::ControllerCache)
    map(setdirty!, cache.link_caches)
    map(setdirty!, cache.point_caches)
    nothing
end

ext: jl


text:
# ## Непосредственная реализация проверки
#
using Test;
using LaTeXStrings;
using Plots;

include("test/direct_check.jl")

# ### Визуализация

# Посмотрим на результат решения в сравнении с аналитическим решением

d = [missing, missing];
dd = [missing missing; missing missing];
make_gif(u, Xₙ, Tₘ, dd, dd, dd, d, d, u_model; name = "dicrect_check.gif")

# Найдем абсолютную погрешность численного решения

err = u .- u_model
heatmap(Xₙ, Tₘ, err', xlabel=L"X_n", ylabel=L"T_m", title="Absolute Error", size=(1200, 800))


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct KeyVaultMetaInfo <: SwaggerModel
    keyVaultResourceId::Any # spec type: Union{ Nothing, String } # spec name: keyVaultResourceId
    encryptionKeyName::Any # spec type: Union{ Nothing, String } # spec name: encryptionKeyName
    encryptionKeyVersion::Any # spec type: Union{ Nothing, String } # spec name: encryptionKeyVersion

    function KeyVaultMetaInfo(;keyVaultResourceId=nothing, encryptionKeyName=nothing, encryptionKeyVersion=nothing)
        o = new()
        validate_property(KeyVaultMetaInfo, Symbol("keyVaultResourceId"), keyVaultResourceId)
        setfield!(o, Symbol("keyVaultResourceId"), keyVaultResourceId)
        validate_property(KeyVaultMetaInfo, Symbol("encryptionKeyName"), encryptionKeyName)
        setfield!(o, Symbol("encryptionKeyName"), encryptionKeyName)
        validate_property(KeyVaultMetaInfo, Symbol("encryptionKeyVersion"), encryptionKeyVersion)
        setfield!(o, Symbol("encryptionKeyVersion"), encryptionKeyVersion)
        o
    end
end # type KeyVaultMetaInfo

const _property_map_KeyVaultMetaInfo = Dict{Symbol,Symbol}(Symbol("keyVaultResourceId")=>Symbol("keyVaultResourceId"), Symbol("encryptionKeyName")=>Symbol("encryptionKeyName"), Symbol("encryptionKeyVersion")=>Symbol("encryptionKeyVersion"))
const _property_types_KeyVaultMetaInfo = Dict{Symbol,String}(Symbol("keyVaultResourceId")=>"String", Symbol("encryptionKeyName")=>"String", Symbol("encryptionKeyVersion")=>"String")
Base.propertynames(::Type{ KeyVaultMetaInfo }) = collect(keys(_property_map_KeyVaultMetaInfo))
Swagger.property_type(::Type{ KeyVaultMetaInfo }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_KeyVaultMetaInfo[name]))}
Swagger.field_name(::Type{ KeyVaultMetaInfo }, property_name::Symbol) =  _property_map_KeyVaultMetaInfo[property_name]

function check_required(o::KeyVaultMetaInfo)
    (getproperty(o, Symbol("keyVaultResourceId")) === nothing) && (return false)
    (getproperty(o, Symbol("encryptionKeyName")) === nothing) && (return false)
    (getproperty(o, Symbol("encryptionKeyVersion")) === nothing) && (return false)
    true
end

function validate_property(::Type{ KeyVaultMetaInfo }, name::Symbol, val)
end


text:
export VBasedPolicy

function default_value_action_mapping(env, value_learner; explorer = GreedyExplorer())
    A = legal_action_space(env)
    V = map(A) do a
        value_learner(child(env, a))
    end
    A[explorer(V)]
end

"""
    VBasedPolicy(;learner, mapping=default_value_action_mapping)

The `learner` must be a value learner. The `mapping` is a function which returns
an action given `env` and the `learner`. By default we iterate through all the
valid actions and select the best one which lead to the maximum state value.
"""
Base.@kwdef struct VBasedPolicy{L,M} <: AbstractPolicy
    learner::L
    mapping::M = default_value_action_mapping
end

(p::VBasedPolicy)(env::AbstractEnv) = p.mapping(env, p.learner)

function RLBase.update!(
    p::VBasedPolicy,
    t::AbstractTrajectory,
    e::AbstractEnv,
    s::AbstractStage,
)
    update!(p.learner, t, e, s)
end


text:
"""
    setDefaultAlgorithm(genj::GeneticJuliaStructure)

documentation
"""
function setDefaultAlgorithm(genj::GeneticJuliaStructure)
    genj._experimentInfo._algorithm = basicExperiment
    genj._experimentInfo._algorithmArgs = Array{Any}(undef, 0)
end # function

"""
    setDefaultStopCondition(genj::GeneticJuliaStructure)

documentation
"""
function setDefaultStopCondition(genj::GeneticJuliaStructure)
    setStopCondition(genj = genj, maxEvaluations = -1, maxIterations = 500, maxIterNotImproving = -1, maxTime = Inf)
end # function

"""
    setDefaultSelector()

documentation
"""
function setDefaultSelector(genj::GeneticJuliaStructure)
    aux = getDefaultSelector()
    setSelector(aux[1], aux[2]..., nSelected=genj._generator._popSize, genj=genj)
end # function

"""
    setDefaultCrossoverOp(genj::GeneticJuliaStructure)

documentation
"""
function setDefaultCrossoverOp(genj::GeneticJuliaStructure)
    aux = getDefaultCrossoverOp(genj._experimentInfo._individualType)
    setCrossoverOperator(aux[1], aux[2]..., genj=genj)
end # function

"""
    setDefaultMutationOp(genj::GeneticJuliaStructure)

documentation
"""
function setDefaultMutationOp(genj::GeneticJuliaStructure)
    aux = getDefaultMutationOp(genj._experimentInfo._individualType)
    setMutationOperator(aux[1], aux[2]..., genj=genj)
end # function

"""
    setDefaultReplacementOp(genj::GeneticJuliaStructure)

documentation
"""
function setDefaultReplacementOp(genj::GeneticJuliaStructure)
    aux = getDefaultReplacementOp()
    setReplacementOperator(aux[1], aux[2]..., genj=genj)
end # function

"""
    setDefaultSettings(; genj::GeneticJuliaStructure = GenJ)

documentation
"""
function setDefaultSettings(; genj::GeneticJuliaStructure = GenJ, i::Integer=0)

    if !isdefined(genj._experimentInfo, :_algorithm)
        setDefaultAlgorithm(genj)
    end

    if !isdefined(genj, :_stopCondition)
        setDefaultStopCondition(genj)
    end

    if isdefined(genj._experimentInfo, :_individualType)
        if isPredefinedType(genj._experimentInfo._individualType)
            if !isdefined(genj, :_crossoverOp)
                setDefaultCrossoverOp(genj)
            end

            if !isdefined(genj, :_mutationOp)
                setDefaultMutationOp(genj)
            end
        end
    else
        if i == 0
            error("individualType is mandatory and must be set")
        else
            error("individualType is mandatory and must be set in experiment $i")
        end
    end

    if !isdefined(genj, :_selector)
        setDefaultSelector(genj)
    else
        if !(typeof(genj._selector._nSelectedParents) <: Integer)
            popSize = genj._generator._popSize
            nSelectedParents = genj._selector._nSelectedParents
            nParents = genj._crossoverOp._nParents

            nSelectedParents = convert(Integer, round(nSelectedParents * popSize))
            remainder = nSelectedParents % nParents

            if remainder != 0
                nSelectedParents = nSelectedParents + nParents - remainder
            end
            if !genj._selector._samplingWithRep
                if nSelectedParents > popSize
                    nSelectedParents -= nParents
                end
            end

            genj._selector = SelectionOperator(genj._selector._method, nSelectedParents,
                                               genj._selector._needsComparison,
                                               genj._selector._samplingWithRep,
                                               genj._selector._varArgs)
        end
    end

    if !isdefined(genj, :_replacementOp)
        setDefaultReplacementOp(genj)
    end
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
##################################################
# Constraints
"""
List of constraints that the model will potentially have
"""
constr_list =
  [
   :startup_choose,
   :startup_restrict,
   :segprod_limit,
   :segprod_limita,
   :segprod_limitb,
   :prod_above_def,
   :prod_limit,
   :str_prod_limit,
   :binary_link,
   :switch_on_off,
   :ramp_up,
   :ramp_down,
   :str_ramp_up,
   :str_ramp_down,
   :startstop_limit,
   :startup_limit,
   :shutdown_limit,
   :min_uptime,
   :min_downtime,
   :power_balance,
   :net_injection_def,
   :min_reserve
  ]


text:
#!/usr/bin/env julia

# Lists all of the quarters of square degrees around the area
# imaged by SDSS. Output in this format can be piped to either infer-box.jl
# or estimate-box-runtime.jl, e.g.
#
#   ./list_quarters.jl | sort -R | xargs -n 4 estimate-box-runtime.jl

for ra in 0.01:.5:380, dec in -30.01:.5:90
    println("$(ra - 0.5) $ra $(dec - 0.5) $dec")
end


text:
## generation_tools.jl
export get_scan_pattern, ϕ_tuple, idx_to_dim

"""
    ϕ_tuple(t::NTuple)

Helper function to obtain the angle from a tuple of coordinates.
The azimuthal angle is calculated via the `atan`.
However, the order of the NTuple is (x,y).

Arguments:
+ `t`: an NTuple. Only the first two coordinates (x,y) are used and the azimuthal angle ϕ is returned
"""
ϕ_tuple(t::NTuple) = atan(t[2],t[1])

## conversion_tools.jl

"""
    idx_to_dim(idx_arr)  # this should be a view

Converts an N-dimensional array of NTuple to an N+1 dimensional array by orienting the (inner) tuple along the (outer) trailing+1 dimension.

Arguments:
+ `idx_arr`. The array of NTuple to convert

"""
function idx_to_dim(idx_arr::AbstractArray{T, N}) where {T, N}
    out_arr = Array{typeof(idx_arr[1][1])}(undef, size(idx_arr)..., length(idx_arr[begin]))

    # loop over array
    @inbounds for I in CartesianIndices(idx_arr)
        for (i, ind) in enumerate(axes(out_arr, N+1))
            out_arr[Tuple(I)..., ind] = idx_arr[I][i]
        end
    end
    return out_arr
end


text:
#=
Philosophy and convention:
A node is defined as the center of a grid cell. Grid cells may be those
in a primary grid or those in a dual grid.

The definitions "dual" and "primal" are made based on their conventional use
in holding fluid dynamic data. For example, pressure is held in primal nodes,
vorticity and streamfunction at dual nodes.

The definitions here are made with a "dual grid" in mind. That is, the grid is
defined by an integer number of dual cells in each direction. If a "primal grid"
is needed, then all of the defintions can be swapped (primal -> dual, dual -> primal).

Also, note that there might be dual cells that are "ghosts" (i.e. lie outside
the grid), but these are not distinguished in these basic definitions and operators.
=#

#=
Notes on data types:
All data types can be reduced to one of two different types in each direction.
(1) Lying on dual cell centers (C), numbered 1 through N
(2) Lying on dual cell edges (E), numbered 1 through N-1 (i.e. lying midway between C)

For example, Nodes{Dual} by definition lie on dual cell centers in both directions,
so they would be C x C. Nodes{Primal} are aligned with the corners of the dual
cells, so they lie along dual cell edges in both directions, E x E.

Nodes{Dual}    -- C x C
Nodes{Primal}  -- E x E
XEdges{Dual}   -- E x C
YEdges{Dual}   -- C x E
XEdges{Primal} -- C x E
YEdges{Primal} -- E x C

This is important when considering the interpolation and differentiation operations,
both of which must account for the relative indexing between E and C. The
key to remember is that E[1] is bounded by C[1] and C[2] in our indexing convention,
so that, for either interpolation or differentiation:
      E[i] <- C[i], C[i+1] for i in 1:N-1
and
      C[i] <- E[i-1], E[i] for i in 2:N-1
and
      C[i] <- C[i] for i in 1:N
and
      E[i] <- E[i] for i in 1:N-1
=#

module Fields

import Base: @propagate_inbounds, show, summary

#using Compat
using FFTW
using SpecialFunctions
using Statistics

using LinearAlgebra
using SparseArrays

import LinearAlgebra: mul!, ldiv!, cross, ×, dot, ⋅
import Base: parentindices
const GAMMA = MathConstants.γ

export Primal, Dual, ScalarGridData, VectorGridData, GridData,
       Edges, Nodes, XEdges, YEdges,
       EdgeGradient, NodePair,
       Points, ScalarData, VectorData, TensorData,
       celltype,
       diff!,interpolate!,
       curl, curl!, Curl, divergence, divergence!, Divergence,
       grad, grad!, Grad,
       laplacian, laplacian!, laplacian_symm!, plan_laplacian, plan_laplacian!,
       helmholtz, helmholtz!, plan_helmholtz, plan_helmholtz!,
       plan_intfact,plan_intfact!,Identity,
       product, product!, ∘,
       directional_derivative!, directional_derivative_conserve!, curl_cross!,
       convective_derivative!, convective_derivative_rot!,
       coordinates,
       DDF, GradDDF,
       Regularize, RegularizationMatrix, InterpolationMatrix,
       CircularConvolution

abstract type CellType end
abstract type Primal <: CellType end
abstract type Dual <: CellType end

abstract type GridData{NX,NY,T} <: AbstractMatrix{T} end

abstract type ScalarGridData{NX,NY,T} <: GridData{NX,NY,T} end

abstract type VectorGridData{NX,NY,T} <: GridData{NX,NY,T} end

macro wraparray(wrapper, field, N)
    S = eval(wrapper)
    @assert S <: AbstractArray "Wrapped type must be a subtype of AbstractArray"
    while supertype(S) <: AbstractArray
        S = supertype(S)
    end
    #T = supertype(eval(wrapper))
    #@assert T <: AbstractArray "Wrapped type must be a subtype of AbstractArray"
    #el_type, N = S.parameters

    quote
        Base.parent(A::$wrapper) = A.$field
        Base.size(A::$wrapper) = size(A.$field)
        parentindices(A::$wrapper) = parentindices(A.$field)

        if $N > 1
          function Base.show(io::IO, m::MIME"text/plain", A::$wrapper)
            println(io, "$(typeof(A)) data")
            println(io, "Printing in grid orientation (lower left is (1,1))")
            show(io,m, reverse(transpose(A.$field),dims=1))
          end
          #function Base.summary(io::IO, A::$wrapper)
          #  println(io, "$(typeof(A)) data")
          #  print(io, "Printing in grid orientation (lower left is (1,1))")
          #end
        end

        @propagate_inbounds Base.getindex(A::$wrapper, i::Int) = A.$field[i]
        @propagate_inbounds Base.setindex!(A::$wrapper, v, i::Int) = A.$field[i] = convert(eltype(A.$field), v)
        if $N > 1
          @propagate_inbounds Base.getindex(A::$wrapper, I::Vararg{Int, $N}) = A.$field[I...]
          @propagate_inbounds Base.setindex!(A::$wrapper, v, I::Vararg{Int, $N}) = A.$field[I...] = convert(eltype(A.$field), v)
        end
    end
end

celltype(w::GridData) = typeof(w).parameters[1]

function othertype end

macro othertype(celltype, k)
    esc(quote
        Fields.othertype(::$celltype) = $k
        Fields.othertype(::Type{$celltype}) = $k
    end)
end

@othertype Primal Dual
@othertype Dual Primal
@othertype CellType CellType

# This macro allows us to access scalar grid data via just the wrapper itself
@wraparray ScalarGridData data 2

include("fields/nodes.jl")
include("fields/edges.jl")
include("fields/collections.jl")

# (ctype,dnx,dny,shiftx,shifty)
scalarlist = ((:(Nodes{Primal,NX,NY,T}), 1,1,0.0,0.0),
              (:(Nodes{Dual,NX,NY,T}),   0,0,0.5,0.5),
              (:(XEdges{Primal,NX,NY,T}),0,1,0.5,0.0),
              (:(YEdges{Primal,NX,NY,T}),1,0,0.0,0.5),
              (:(XEdges{Dual,NX,NY,T}),  1,0,0.0,0.5),
              (:(YEdges{Dual,NX,NY,T}),  0,1,0.5,0.0))

# (ctype,dunx,duny,dvnx,dvny,shiftux,shiftuy,shiftvx,shiftvy)
vectorlist = ((:(Edges{Primal,NX,NY,T}),          0,1,1,0,0.5,0.0,0.0,0.5),
              (:(Edges{Dual,NX,NY,T}),            1,0,0,1,0.0,0.5,0.5,0.0),
              (:(NodePair{Dual,Dual,NX,NY,T}),    0,0,0,0,0.5,0.5,0.5,0.5),
              (:(NodePair{Primal,Dual,NX,NY,T}),  1,1,0,0,0.0,0.0,0.5,0.5),
              (:(NodePair{Dual,Primal,NX,NY,T}),  0,0,1,1,0.5,0.5,0.0,0.0),
              (:(NodePair{Primal,Primal,NX,NY,T}),1,1,1,1,0.0,0.0,0.0,0.0))

tensorlist = ((:(EdgeGradient{Dual,Primal,NX,NY,T}), 0,0,1,1,0.5,0.5,0.0,0.0),
              (:(EdgeGradient{Primal,Dual,NX,NY,T}), 1,1,0,0,0.0,0.0,0.5,0.5))

include("fields/basicoperations.jl")
include("fields/points.jl")

CollectedData = Union{EdgeGradient{R,S,NX,NY,T},NodePair{R,S,NX,NY,T}} where {R,S,NX,NY,T}

"""
    coordinates(w::GridData;[dx=1.0],[I0=(1,1)])

Return a tuple of the ranges of the physical coordinates in each direction for grid
data `w`. If `w` is of `Nodes` type, then it returns a tuple of the form
`xg,yg`. If `w` is of `Edges` or `NodePair` type, then it returns a tuple of
the form `xgu,ygu,xgv,ygv`.

The optional keyword argument `dx` sets the grid spacing; its default is `1.0`. The
optional keyword `I0` accepts a tuple of integers to set the index pair of the
primal nodes that coincide with the origin. The default is `(1,1)`.

# Example

```jldoctest
julia> w = Nodes(Dual,(12,22));

julia> xg, yg = coordinates(w,dx=0.1)
(-0.05:0.1:1.05, -0.05:0.1:2.0500000000000003)
```
"""
function coordinates end

for (ctype,dnx,dny,shiftx,shifty) in scalarlist
   @eval coordinates(w::$ctype;dx::Float64=1.0,I0::Tuple{Int,Int}=(1,1)) where {NX,NY,T} =
    dx.*((1-I0[1]-$shiftx):(NX-$dnx-I0[1]-$shiftx),
         (1-I0[2]-$shifty):(NY-$dny-I0[2]-$shifty))

end

for (ctype,dunx,duny,dvnx,dvny,shiftux,shiftuy,shiftvx,shiftvy) in vectorlist
   @eval coordinates(w::$ctype;dx::Float64=1.0,I0::Tuple{Int,Int}=(1,1)) where {NX,NY,T} =
    dx.*((1-I0[1]-$shiftux):(NX-$dunx-I0[1]-$shiftux),
         (1-I0[2]-$shiftuy):(NY-$duny-I0[2]-$shiftuy),
         (1-I0[1]-$shiftvx):(NX-$dvnx-I0[1]-$shiftvx),
         (1-I0[2]-$shiftvy):(NY-$dvny-I0[2]-$shiftvy))

end

include("fields/physicalgrid.jl")
include("fields/operators.jl")

end


text:
#=
julia custom startup file
to be placed at: ~/.julia/config/startup.jl

NOTE: this file is for julia versions v0.7 onward
=#

using REPL
function setup_repl(repl)
    repl.interface = REPL.setup_interface(repl)
    # default input prompt
    repl.interface.modes[1].prompt = ":> "
    # shell/terminal prompt
    repl.interface.modes[2].prompt = "\$> "
    # help prompt
    repl.interface.modes[3].prompt = "?> "
    #  # package manager prompt
    #  repl.interface.modes[6].prompt = "[> "
    return
end
atreplinit(setup_repl)
const global JULIA_HOME = expanduser("~/Applications/Julia/master/usr/bin")

max_stars_repo_name: dysonance/Preferences

max_issues_repo_name: dysonance/Preferences

max_forks_repo_name: dysonance/Preferences


text:
using MonteCarloIntegration
using HCubature
using Test

integrands = [
              x -> 1,
              x -> sum(x),
              x -> sum(sin.(x))
             ]

batch_f(f) = (pts) -> begin

    npts = size(pts, 1)
    ndims = size(pts, 2)
    fevals = zeros(npts)

    for i = 1:npts
        p = vec(pts[i,:])
        fevals[i] = f(p)
    end

    fevals
end

@testset "Standard Integrands" begin

    for f in integrands

        for dim = 1:5
            @info("Dimension = $dim")
            v, _ = vegas(f, zeros(dim), ones(dim))
            h, _ = hcubature(f, zeros(dim), ones(dim))
            @test isapprox(v, h, rtol = 1e-2)
        end
    end
end

@testset "Batched Integrands" begin

    for f in integrands

        f2 = batch_f(f)

        for dim = 1:5
            @info("Dimension = $dim")
            v, _ = vegas(f2, zeros(dim), ones(dim), batch = true)
            h, _ = hcubature(f, zeros(dim), ones(dim))
            @test isapprox(v, h, rtol = 1e-2)
        end
    end
end

#=@testset "Simple Monte Carlo Tests" begin
    for dim = 2:2
        @test isapprox(vegas(x -> 1, zeros(dim), ones(dim), nbins = 1)[1], 1., rtol = 1e-2)
        @test isapprox(vegas(x -> sum(x),zeros(dim), ones(dim), nbins = 1)[1], 1., rtol = 1e-2)
        @test isapprox(vegas(x -> sum(sin.(x)),zeros(dim), ones(dim), nbins = 1)[1], 0.919, rtol = 1e-2)
    end
end
@testset "Standard Integrands" begin
    for dim = 2:2
        @test isapprox(vegas(x -> 1, zeros(dim), ones(dim))[1], 1., rtol = 1e-2)
        @test isapprox(vegas(x -> sum(x), zeros(dim), ones(dim))[1], 1., rtol = 1e-2)
        @test isapprox(vegas(x -> sum(sin.(x)), zeros(dim), ones(dim))[1], 0.919, rtol = 1e-2)
    end
end=#


text:
using Pkg
if isfile("Project.toml") && isfile("Manifest.toml")
    Pkg.activate(".")
end

###################
# Terminal ########
###################

    usedayperiod = true
    startday = 201
    lastday = 269
    starthour = 1
    lasthour = 8760
	x_ag = 0.3
	x_re = 0.3
	x_ind = 0.4
	s_ag = 0.05
	# 0.025
	#0.0125
	s_re = 0.05
	# 0.025
	# 0.0125
	s_ind = 0.1
	# 0.05
	# 0.025
	workhours = 9:16
	trm = 0.25
    datapath = "data/"
	outputpath = "output/Output_TSOM/"

   if usedayperiod == true
        timeperiod = string((startday-1)*24+1)*":"*string((lastday*24))
   else
        timeperiod = string(starthour)*":"*string(lasthour)
   end

    timeperiod_str = replace(string(timeperiod), ":" => "-")

    println()
    # println("Case >" * case * "< selected.")
    println("Calculating for time period >" * timeperiod * "<.")
    println()
	isworkhour = Dict()
    for i in 0:24
        if i in workhours isworkhour[i] = 1
        else isworkhour[i] = 0
        end
    end

###################
# Run CM ##########
###################

    println("Loading required packages: ProgressMeter, CSV, DataFrames, Dates, JuMP, Gurobi.")
    using ProgressMeter
    using CSV
    using DataFrames
    using Dates
    using JuMP
    using Gurobi
    println("Loading of packages complete.")

    include("model/src/Sesam_ptdf_LoadShift.jl") #WHEN USING PTDF
    # include("model/src/Sesam.jl") WHEN USING NO PTDF
    using .Sesam

###################
# Data ############
###################
# Read nodes and load from .csv
    println("Reading data from csv...")
    prog = Progress(14)
	nodes_df = CSV.read(joinpath(datapath, "nodes.csv"))
		next!(prog)
	load_df = CSV.read(joinpath(datapath, "load.csv"))
		next!(prog)
	exchange_df = CSV.read(joinpath(datapath, "exchange.csv"))
		next!(prog)
	# Read lines from .csv
	lines_df = CSV.read(joinpath(datapath, "lines.csv"))
		next!(prog)
	# Read power plants from .csv
	powerplants_df = CSV.read(joinpath(datapath, "powerplants.csv"))
		next!(prog)
	avail_powerplants_df = CSV.read(joinpath(datapath, "avail_powerplants.csv"), header = 2, skipto = 3)
		next!(prog)
	mustrun_powerplants_df = CSV.read(joinpath(datapath, "mustrun_powerplants.csv"), header = 2, skipto = 3)
		next!(prog)
	# Read fuel costs from .csv
	fuelcost_df = CSV.read(joinpath(datapath, "fuelcost.csv"))
		next!(prog)
	# Read renewables and availabilities from .csv
	renewables_df = CSV.read(joinpath(datapath, "renewables.csv"))
		next!(prog)
	avail_solar_pv_df = CSV.read(joinpath(datapath, "avail_solar_pv.csv"))
		next!(prog)
	avail_wind_offshore_df = CSV.read(joinpath(datapath, "avail_wind_offshore.csv"))
		next!(prog)
	avail_wind_onshore_df = CSV.read(joinpath(datapath, "avail_wind_onshore.csv"))
		next!(prog)
	avail_other_res_df = CSV.read(joinpath(datapath, "avail_other_res.csv"))
		next!(prog)
	renewables_avail_dict = Dict(:SolarPV => avail_solar_pv_df,
		:WindOffshore => avail_wind_offshore_df,
		:WindOnshore => avail_wind_onshore_df,
		:other => avail_other_res_df)
		next!(prog)

    println("Reading of data complete.")
    println()

###################
# Model setup #####
###################

		println("Preparing model technologies, nodes, and grid...")
		prog = Progress(6)

	# Create sets
		# Create set of nodes and lines
		nodes = Nodes(nodes_df, load_df, exchange_df, 1000, 1)
			next!(prog)
		lines = Lines(lines_df, 1000, trm)
			next!(prog)
		# Create set of power plants	ED_L_S_opt, ED_p
		pp = PowerPlants(powerplants_df, avail_powerplants_df, mustrun_powerplants_df, fuelcost_df)
			next!(prog)
		# Create set of renewables
		res = Renewables(renewables_df, renewables_avail_dict)
			next!(prog)
		# Create a set of storages
		storages_df = CSV.read(joinpath(datapath, "storages.csv"))
			next!(prog)
		storages = Storages(storages_df)
			next!(prog)

	# Create Dict Zone to Gen
		N = nodes.id
		P = pp.unit
		R = res.unit
		S = storages.unit
		Node2PP = Dict((i, []) for i in nodes.id)
		for i in keys(pp.node)
			push!(Node2PP[pp.node[i]], i)
		end
		Node2Res = Dict((i, []) for i in nodes.id)
		for i in keys(res.node)
			push!(Node2Res[res.node[i]], i)
		end
		Node2Stor = Dict((i, []) for i in nodes.id)
		for i in keys(storages.node)
			push!(Node2Stor[storages.node[i]], i)
		end

	println("Preparation complete.")
	println()

		session_time = Dates.now()
		session = Dates.format(session_time, "yyyy-mm-dd__HH-MM-SS")*"_ptdf_LoadShift"

	# Create session folder for output

		mkdir(joinpath(outputpath, session*"__"*timeperiod_str))

	prog = Progress(length(dayslicer(timeperiod)))

	### Parameter outputs
		println("Exporting parameters.")
		println("Exporting merit order...")
		# Merit order
		output_ED_merit_order_pp = DataFrame(unit=Symbol[],fuel=String[],capacity=Float64[],mc=Float64[])

		try
			for g in pp.unit
					mc = pp.mc[g]
					fuel = pp.fuel[g]
					capacity = pp.capacity[g]
					push!(output_ED_merit_order_pp, [g fuel capacity mc])
			end

			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_merit_order_pp.csv"), output_ED_merit_order_pp)
		catch e
			println("The marginal costs for power plants are not available.")
		end

		output_ED_merit_order_res = DataFrame(unit=Symbol[],fuel=Symbol[],capacity=Float64[],mc=Float64[])

		try
			for r in res.unit
					mc = 0
					fuel = res.fuel[r]
					capacity = res.capacity[r]
					push!(output_ED_merit_order_res, [r fuel capacity mc])
			end

			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_merit_order_res.csv"), output_ED_merit_order_res)
		catch e
			println("The marginal costs for renewables are not available.")
		end

		output_ED_merit_order_storages = DataFrame(unit=Symbol[],fuel=String[],capacity=Float64[],mc=Float64[])

		try
			for s in storages.unit
					mc = 0
					fuel = storages.fuel[s]
					capacity = storages.power[s]
					push!(output_ED_merit_order_storages, [s fuel capacity mc])
			end

			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_merit_order_storages.csv"), output_ED_merit_order_storages)
		catch e
			println("The marginal costs for renewables are not available.")
		end

	println("Export complete.")
	println()

### Initial run of first dayslice
	println("Initial run. Calculating first dayslice: " * string(dayslicer(timeperiod)[1]))
	println()
	println()

	ED_mod, ED_P_opt, ED_P_R_opt, ED_P_R_opt_dist, ED_P_S_opt, ED_D_S_opt,
		ED_L_S_opt, ED_price = EconomicDispatch80(dayslicer(timeperiod)[1], nodes, pp, res, storages)

	CM_mod, CM_ΔP_up_opt, CM_ΔP_dn_opt, CM_ΔP_R_up_opt, CM_ΔP_R_dn_opt,
		CM_Θ_opt, CM_P_flow_opt, CM_price, CM_P_gen_lost_opt, CM_P_load_lost_opt, CM_P_S_up_opt, CM_Load_Shifted_Opt, CM_LoadShift_up_opt, CM_LoadShift_dn_opt = CongestionManagementPHS80_LoadShift(dayslicer(timeperiod)[1],
		nodes, lines, pp, res, storages, ED_price, ED_P_opt, ED_P_R_opt_dist, ED_P_S_opt, ED_D_S_opt, ED_L_S_opt, x_ag, x_re, x_ind, s_ag, s_re, s_ind, isworkhour)
CM_mod
	next!(prog)

			LoadShift_df = DataFrame(time=Int64[], node=Int64[], OriginalLoad=Float64[], ShiftedLoad=Float64[], Shift=Float64[])
				for t in dayslicer(timeperiod)[1]
					for n in N
						originalload = nodes.load[n][t]
						shiftedload = CM_Load_Shifted_Opt[t, n]
						shift = CM_LoadShift_up_opt[t, n] - CM_LoadShift_dn_opt[t, n]
						push!(LoadShift_df, [t n originalload shiftedload shift])
					end
				end
				CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_LoadShift.csv"), LoadShift_df)

			C_V_RD_N = DataFrame(time=Int64[], node=Int64[], CostPP=Float64[], CostRes=Float64[], CostStor=Float64[], CostLostLoad=Float64[], CostLostGen=Float64[], CostTotal=Float64[],
				VolumePP=Float64[], VolumeRes=Float64[], VolumeStor=Float64[], VolumeLostLoad=Float64[], VolumeLostGen=Float64[], VolumeTotal=Float64[])
				for t in dayslicer(timeperiod)[1]
					for n in N
						if isempty(Node2PP[n])
							costpp = 0
						else
							costpp = sum((pp.mc[g][1]*CM_ΔP_up_opt[t, g] + (ED_price[t] - pp.mc[g][]) * CM_ΔP_dn_opt[t, g]) for g in Node2PP[n])
						end
						if isempty(Node2Res[n])
							costres = 0
						else
							costres = sum(((ED_price[t]) * CM_ΔP_R_dn_opt[t, res]) for res in Node2Res[n])
						end
						if isempty(Node2Stor[n])
							coststor = 0
						else
							coststor = sum(((ED_price[t]/storages.efficiency[stor]) * CM_P_S_up_opt[t, stor]) for stor in Node2Stor[n])
						end

						costlostload = 1000*CM_P_load_lost_opt[t, n]
						costlostgen = 1000*CM_P_gen_lost_opt[t, n] # ÜBERALL NULL!
						totalcost = costpp + costres + coststor + costlostload + costlostgen

						if isempty(Node2PP[n])
							volumepp = 0
						else
							volumepp = sum((CM_ΔP_up_opt[t, g] - CM_ΔP_dn_opt[t, g]) for g in Node2PP[n])
						end
						if isempty(Node2Res[n])
							volumeres = 0
						else
							volumeres = sum(-(CM_ΔP_R_dn_opt[t, res]) for res in Node2Res[n])
						end
						if isempty(Node2Stor[n])
							volumestor = 0
						else
							volumestor = sum((CM_P_S_up_opt[t, stor]) for stor in Node2Stor[n])
						end

						volumelostload = CM_P_load_lost_opt[t, n]
						volumelostgen = CM_P_gen_lost_opt[t, n]
						totalvolume = volumepp + volumeres + volumestor + volumelostload + volumelostgen

						push!(C_V_RD_N, [t n costpp costres coststor costlostload costlostgen totalcost volumepp volumeres volumestor volumelostload volumelostgen totalvolume])
					end
				end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "C_V_RD_N.csv"), C_V_RD_N)

	# Save ED Data
		# Dispatchable power plants
		output_ED_P_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])

		try
			for t in dayslicer(timeperiod)[1]
				for g in pp.unit
					val = ED_P_opt[t, g]
					fuel =  pp.fuel[g]
					push!(output_ED_P_opt, [t g fuel val pp.node[g]])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_P.csv"), output_ED_P_opt)
		catch e
			println("ED_P_opt does not exist.")
		end

		# RES generation
		output_ED_P_R_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=Symbol[],output=Float64[],node=Int64[])
		try
			for t in dayslicer(timeperiod)[1]
				for r in res.unit
					val = ED_P_R_opt[t, r]
					fuel =  res.fuel[r]
					push!(output_ED_P_R_opt, [t r fuel val res.node[r]])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_P_R.csv"), output_ED_P_R_opt)
		catch e
			println("ED_P_R_opt does not exist.")
		end

		# Maximum available RES infeed from input data
		output_MaxRES = DataFrame(time=Int64[],unit=Symbol[],fuel=Symbol[],output=Float64[],node=Int64[])
		try
			for t in dayslicer(timeperiod)[1]
				for r in res.unit
					output = res.infeed[r][t]
					fuel =  res.fuel[r]
					push!(output_MaxRES, [t r fuel output res.node[r]])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "MaxRES.csv"), output_MaxRES)
		catch e
			println("res.infeed does not exist.")
		end

		# Generation from storages
		output_ED_P_S_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
		try
			for t in dayslicer(timeperiod)[1]
				for s in storages.unit
					val = ED_P_S_opt[t, s]
					fuel =  storages.fuel[s]
					push!(output_ED_P_S_opt, [t s fuel val storages.node[s]])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_P_S.csv"), output_ED_P_S_opt)
		catch e
			println("ED_P_S_opt does not exist.")
		end

		# Electricity demand from storages
		output_ED_D_S_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
		try
			for t in dayslicer(timeperiod)[1]
				for s in storages.unit
					val = -1*ED_D_S_opt[t, s]
					fuel =  storages.fuel[s]
					push!(output_ED_D_S_opt, [t s fuel val storages.node[s]])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_D_S.csv"), output_ED_D_S_opt)
		catch e
			println("ED_D_S_opt does not exist.")
		end

		# Storage level
		output_ED_L_S_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],storage=Float64[],node=Int64[])
		try
			for t in dayslicer(timeperiod)[1]
				for s in storages.unit
					val = ED_L_S_opt[t, s]
					fuel =  storages.fuel[s]
					push!(output_ED_L_S_opt, [t s fuel val storages.node[s]])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_L_S.csv"), output_ED_L_S_opt)
		catch e
			println("ED_L_S_opt does not exist.")
		end

		# Load
		output_load = DataFrame(time=Int64[],load=Float64[])
		try
			for t in dayslicer(timeperiod)[1]
				load = nodes.systemload[t]
				push!(output_load, [t load])
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "load.csv"), output_load)
		catch e
			println("nodes.systemload does not exist.")
		end

		# Prices
		output_ED_price = DataFrame(time=Int64[],price=Float64[])
		try
			for t in dayslicer(timeperiod)[1]
				price =  ED_price[t]
				push!(output_ED_price, [t price])
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_price.csv"), output_ED_price)
		catch e
			println("ED_price does not exist.")
		end

		# Must-run
		output_pmin = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
		try
			for t in dayslicer(timeperiod)[1]
				for g in pp.unit
					val = pp.pmin[g][t]
					fuel =  pp.fuel[g]
					push!(output_pmin, [t g fuel val pp.node[g]])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "pmin.csv"), output_pmin)
		catch e
			println("pp.pmin does not exist.")
		end

		### Model parameters
		output_ED_info = DataFrame(timeperiod=String[],
										termination_status=String[],
										primal_status=String[],
										dual_status=String[],
										objective_value=Float64[],
										solve_time=Float64[])

		push!(output_ED_info, [string(dayslicer(timeperiod)[1]) string(termination_status(ED_mod)) string(primal_status(ED_mod)) string(dual_status(ED_mod)) objective_value(ED_mod) solve_time(ED_mod)])
		CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_info.csv"), output_ED_info)

		# Nodal exchange
		output_exchange = DataFrame(time=Int64[],node=Int64[],output=Float64[])
		try
			for t in dayslicer(timeperiod)[1]
				for j in nodes.id
					push!(output_exchange, [t j nodes.exchange[j][t]])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "exchange.csv"), output_exchange)
		catch e
			println("output_exchange does not exist.")
		end

###################
# CM ##############
###################
# Save CM Data
		# Dispatchable power plants
		output_CM_ΔP_up_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
		try
			for t in dayslicer(timeperiod)[1]
				for g in pp.unit
					val = CM_ΔP_up_opt[t, g]
					fuel =  pp.fuel[g]
					push!(output_CM_ΔP_up_opt, [t g fuel val pp.node[g]])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_up.csv"), output_CM_ΔP_up_opt)
		catch e
			println("CM_ΔP_up_opt does not exist.")
		end

		output_CM_ΔP_dn_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
		try
			for t in dayslicer(timeperiod)[1]
				for g in pp.unit
					val = -1*CM_ΔP_dn_opt[t, g]
					fuel =  pp.fuel[g]
					push!(output_CM_ΔP_dn_opt, [t g fuel val pp.node[g]])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_dn.csv"), output_CM_ΔP_dn_opt)
		catch e
			println("CM_ΔP_dn_opt does not exist.")
		end

		# RES generation
		output_CM_ΔP_R_up_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=Symbol[],output=Float64[],node=Int64[])
		try
			for t in dayslicer(timeperiod)[1]
				for r in res.unit
					val = CM_ΔP_R_up_opt[t, r]
					fuel =  res.fuel[r]
					push!(output_CM_ΔP_R_up_opt, [t r fuel val res.node[r]])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_R_up.csv"), output_CM_ΔP_R_up_opt)
		catch e
			println("CM_ΔP_R_up_opt does not exist.")
		end

		output_CM_ΔP_R_dn_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=Symbol[],output=Float64[],node=Int64[])
		try
			for t in dayslicer(timeperiod)[1]
				for r in res.unit
					val = -1*CM_ΔP_R_dn_opt[t, r]
					fuel =  res.fuel[r]
					push!(output_CM_ΔP_R_dn_opt, [t r fuel val res.node[r]])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_R_dn.csv"), output_CM_ΔP_R_dn_opt)
		catch e
			println("CM_ΔP_R_dn_opt does not exist.")
		end

		### Model parameters
		output_CM_info = DataFrame(timeperiod=String[],
										termination_status=String[],
										primal_status=String[],
										dual_status=String[],
										objective_value=Float64[],
										solve_time=Float64[])

		push!(output_CM_info, [string(dayslicer(timeperiod)[1]) string(termination_status(CM_mod)) string(primal_status(CM_mod)) string(dual_status(CM_mod)) objective_value(CM_mod) solve_time(CM_mod)])
		CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_info.csv"), output_CM_info)

		### Power flow results
		output_CM_P_flow_opt = DataFrame(time=Int64[],
										line=Int64[],
										from=Int64[],
										to=Int64[],
										pmax=Float64[],
										pflow=Float64[])

		try
			for t in dayslicer(timeperiod)[1]
				for l in lines.id
					from = parse(Int, string(lines.from[l]))
					to = parse(Int, string(lines.to[l]))
					pmax = lines.pmax[l]
					pflow = CM_P_flow_opt[t, l]
					push!(output_CM_P_flow_opt, [t l from to pmax pflow])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_flow.csv"), output_CM_P_flow_opt)
		catch e
			println("CM_P_flow_opt does not exist.")
		end

		# Lost generation
		output_CM_P_gen_lost = DataFrame(time=Int64[],node=Int64[],output=Float64[])
		try
			for t in dayslicer(timeperiod)[1]
				for j in nodes.id
					output = CM_P_gen_lost_opt[t, j]
					push!(output_CM_P_gen_lost, [t, j, output])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_gen_lost.csv"), output_CM_P_gen_lost)
		catch e
			println("CM_P_gen_lost does not exist.")
		end

		# Lost load
		output_CM_P_load_lost = DataFrame(time=Int64[],node=Int64[],output=Float64[])
		try
			for t in dayslicer(timeperiod)[1]
				for j in nodes.id
					output = CM_P_load_lost_opt[t, j]
					push!(output_CM_P_load_lost, [t, j, output])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_load_lost.csv"), output_CM_P_load_lost)
		catch e
			println("CM_P_load_lost does not exist.")
		end

		output_CM_price = DataFrame(time=Int64[],node=Int64[],price=Float64[])
		try
			for t in dayslicer(timeperiod)[1]
				for j in nodes.id
					price =  CM_price[j, t]
					push!(output_CM_price, [t j price])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_price.csv"), output_CM_price)
		catch e
			println("CM_price does not exist.")
		end

		# PHS output increase
		output_CM_P_S_up_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
		try
			for t in dayslicer(timeperiod)[1]
				for s in storages.unit
					val = CM_P_S_up_opt[t, s]
					fuel =  storages.fuel[s]
					push!(output_CM_P_S_up_opt, [t s fuel val storages.node[s]])
				end
			end
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_S_up.csv"), output_CM_P_S_up_opt)
		catch e
			println("CM_P_S_up_opt does not exist.")
		end

###	Remaining periods: Initialising first hour with last hour of previous day
###	and append data
###################
### Calculate ED ##
###################
	for hours in dayslicer(timeperiod)[2:end]
		hours_extended = [hours[1]-1; hours]
		println()
		println()
		println("Calculating economic dispatch: " * string(hours))

	global ED_mod, ED_P_opt, ED_P_R_opt, ED_P_R_opt_dist, ED_P_S_opt, ED_D_S_opt,
		ED_L_S_opt, ED_price = EconomicDispatch80(hours, nodes, pp, res, storages)

	    ###################
    # Export ED ED ##############
	    ###################

	    # Dispatchable power plants
	    output_ED_P_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])

	    try
	    	for t in hours
	    		for g in pp.unit
	    			val = ED_P_opt[t, g]
	    			fuel =  pp.fuel[g]
	    			push!(output_ED_P_opt, [t g fuel val pp.node[g]])
	    		end
	    	end
	    	CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_P.csv"), output_ED_P_opt, append = true)
	    catch e
	    	println("ED_P_opt does not exist.")
	    end

	    # RES generation
	    output_ED_P_R_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=Symbol[],output=Float64[],node=Int64[])
	    try
	    	for t in hours
	    		for r in res.unit
	    			val = ED_P_R_opt[t, r]
	    			fuel =  res.fuel[r]
	    			push!(output_ED_P_R_opt, [t r fuel val res.node[r]])
	    		end
	    	end
	    	CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_P_R.csv"), output_ED_P_R_opt, append = true)
	    catch e
	    	println("ED_P_R_opt does not exist.")
	    end

	    # Maximum available RES infeed from input data
	    output_MaxRES = DataFrame(time=Int64[],unit=Symbol[],fuel=Symbol[],output=Float64[],node=Int64[])
	    try
	    	for t in hours
	    		for r in res.unit
	    			output = res.infeed[r][t]
	    			fuel =  res.fuel[r]
	    			push!(output_MaxRES, [t r fuel output res.node[r]])
	    		end
	    	end
	    	CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "MaxRES.csv"), output_MaxRES, append = true)
	    catch e
	    	println("res.infeed does not exist.")
	    end

	    # Generation from storages
	    output_ED_P_S_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
	    try
	    	for t in hours
	    		for s in storages.unit
	    			val = ED_P_S_opt[t, s]
	    			fuel =  storages.fuel[s]
	    			push!(output_ED_P_S_opt, [t s fuel val storages.node[s]])
	    		end
	    	end
	    	CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_P_S.csv"), output_ED_P_S_opt, append = true)
	    catch e
	    	println("ED_P_S_opt does not exist.")
	    end

	    # Electricity demand from storages
	    output_ED_D_S_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
	    try
	    	for t in hours
	    		for s in storages.unit
	    			val = -1*ED_D_S_opt[t, s]
	    			fuel =  storages.fuel[s]
	    			push!(output_ED_D_S_opt, [t s fuel val storages.node[s]])
	    		end
	    	end
	    	CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_D_S.csv"), output_ED_D_S_opt, append = true)
	    catch e
	    	println("ED_D_S_opt does not exist.")
	    end

	    # Storage level
	    output_ED_L_S_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],storage=Float64[],node=Int64[])
	    try
	    	for t in hours
	    		for s in storages.unit
	    			val = ED_L_S_opt[t, s]
	    			fuel =  storages.fuel[s]
	    			push!(output_ED_L_S_opt, [t s fuel val storages.node[s]])
	    		end
	    	end
	    	CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_L_S.csv"), output_ED_L_S_opt, append = true)
	    catch e
	    	println("ED_L_S_opt does not exist.")
	    end

	    # Load
	    output_load = DataFrame(time=Int64[],load=Float64[])
	    try
	    	for t in hours
	    		load = nodes.systemload[t]
	    		push!(output_load, [t load])
	    	end
	    	CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "load.csv"), output_load, append = true)
	    catch e
	    	println("nodes.systemload does not exist.")
	    end

	    # Prices
	    output_ED_price = DataFrame(time=Int64[],price=Float64[])
	    try
	    	for t in hours
	    		price =  ED_price[t]
	    		push!(output_ED_price, [t price])
	    	end
	    	CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_price.csv"), output_ED_price, append = true)
	    catch e
	    	println("ED_price does not exist.")
	    end

	    # Must-run
	    output_pmin = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
	    try
	    	for t in hours
	    		for g in pp.unit
	    			val = pp.pmin[g][t]
	    			fuel =  pp.fuel[g]
	    			push!(output_pmin, [t g fuel val pp.node[g]])
	    		end
	    	end
	    	CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "pmin.csv"), output_pmin, append = true)
	    catch e
	    	println("pp.pmin does not exist.")
	    end

		### Model parameters
		output_ED_info = DataFrame(timeperiod=String[],
										   termination_status=String[],
										   primal_status=String[],
										   dual_status=String[],
										   objective_value=Float64[],
										   solve_time=Float64[])

		push!(output_ED_info, [string(hours) string(termination_status(ED_mod)) string(primal_status(ED_mod)) string(dual_status(ED_mod)) objective_value(ED_mod) solve_time(ED_mod)])
		CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "ED_info.csv"), output_ED_info, append = true)

###################
### Calculate CM ##
###################
		println()
		println()
		println("Calculating congestion management: " * string(hours))

	global CM_mod, CM_ΔP_up_opt, CM_ΔP_dn_opt, CM_ΔP_R_up_opt, CM_ΔP_R_dn_opt,
		CM_Θ_opt, CM_P_flow_opt, CM_price, CM_P_gen_lost_opt, CM_P_load_lost_opt, CM_P_S_up_opt, CM_Load_Shifted_Opt, CM_LoadShift_up_opt, CM_LoadShift_dn_opt = CongestionManagementPHS80_LoadShift(hours, nodes, lines, pp, res, storages, ED_price,
		ED_P_opt,  ED_P_R_opt_dist, ED_P_S_opt, ED_D_S_opt, ED_L_S_opt, x_ag, x_re, x_ind, s_ag, s_re, s_ind, isworkhour)

			#	Save cost and volume of re-dispatch

		LoadShift_df = DataFrame(time=Int64[], node=Int64[], OriginalLoad=Float64[], ShiftedLoad=Float64[], Shift=Float64[])
		for t in hours
			for n in N
				originalload = nodes.load[n][t]
				shiftedload = CM_Load_Shifted_Opt[t, n]
				shift = CM_LoadShift_up_opt[t, n] - CM_LoadShift_dn_opt[t, n]
				push!(LoadShift_df, [t n originalload shiftedload shift])
			end
		end
		CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_LoadShift.csv"), LoadShift_df, append = true)

		C_V_RD_N = DataFrame(time=Int64[], node=Int64[], CostPP=Float64[], CostRes=Float64[], CostStor=Float64[], CostLostLoad=Float64[], CostLostGen=Float64[], CostTotal=Float64[],
			VolumePP=Float64[], VolumeRes=Float64[], VolumeStor=Float64[], VolumeLostLoad=Float64[], VolumeLostGen=Float64[], VolumeTotal=Float64[])
	    for t in hours
			for n in N
				if isempty(Node2PP[n])
					costpp = 0
				else
					costpp = sum((pp.mc[g][1]*CM_ΔP_up_opt[t, g] + (ED_price[t] - pp.mc[g][]) * CM_ΔP_dn_opt[t, g]) for g in Node2PP[n])
				end
				if isempty(Node2Res[n])
					costres = 0
				else
					costres = sum(((ED_price[t]) * CM_ΔP_R_dn_opt[t, res]) for res in Node2Res[n])
				end
				if isempty(Node2Stor[n])
					coststor = 0
				else
					coststor = sum(((ED_price[t]/storages.efficiency[stor]) * CM_P_S_up_opt[t, stor]) for stor in Node2Stor[n])
				end

				costlostload = 1000*CM_P_load_lost_opt[t, n]
				costlostgen = 1000*CM_P_gen_lost_opt[t, n] # ÜBERALL NULL!
				totalcost = costpp + costres + coststor + costlostload + costlostgen

				if isempty(Node2PP[n])
					volumepp = 0
				else
					volumepp = sum((CM_ΔP_up_opt[t, g] - CM_ΔP_dn_opt[t, g]) for g in Node2PP[n])
				end
				if isempty(Node2Res[n])
					volumeres = 0
				else
					volumeres = sum(-(CM_ΔP_R_dn_opt[t, res]) for res in Node2Res[n])
				end
				if isempty(Node2Stor[n])
					volumestor = 0
				else
					volumestor = sum((CM_P_S_up_opt[t, stor]) for stor in Node2Stor[n])
				end

				volumelostload = CM_P_load_lost_opt[t, n]
				volumelostgen = CM_P_gen_lost_opt[t, n]
				totalvolume = volumepp + volumeres + volumestor + volumelostload + volumelostgen

				push!(C_V_RD_N, [t n costpp costres coststor costlostload costlostgen totalcost volumepp volumeres volumestor volumelostload volumelostgen totalvolume])
			end
		end
		CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "C_V_RD_N.csv"), C_V_RD_N, append = true)

		# Dispatchable power plants
			output_CM_ΔP_up_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
			try
				for t in hours
					for g in pp.unit
						val = CM_ΔP_up_opt[t, g]
						fuel =  pp.fuel[g]
						push!(output_CM_ΔP_up_opt, [t g fuel val pp.node[g]])
					end
				end
				CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_up.csv"), output_CM_ΔP_up_opt, append = true)
			catch e
				println("CM_ΔP_up_opt does not exist.")
			end

			output_CM_ΔP_dn_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
			try
				for t in hours
					for g in pp.unit
						val = -1*CM_ΔP_dn_opt[t, g]
						fuel =  pp.fuel[g]
						push!(output_CM_ΔP_dn_opt, [t g fuel val pp.node[g]])
					end
				end
				CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_dn.csv"), output_CM_ΔP_dn_opt, append = true)
			catch e
				println("CM_ΔP_dn_opt does not exist.")
			end

			# RES generation
			output_CM_ΔP_R_up_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=Symbol[],output=Float64[],node=Int64[])
			try
				for t in hours
					for r in res.unit
						val = CM_ΔP_R_up_opt[t, r]
						fuel =  res.fuel[r]
						push!(output_CM_ΔP_R_up_opt, [t r fuel val res.node[r]])
					end
				end
				CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_R_up.csv"), output_CM_ΔP_R_up_opt, append = true)
			catch e
				println("CM_ΔP_R_up_opt does not exist.")
			end

			output_CM_ΔP_R_dn_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=Symbol[],output=Float64[],node=Int64[])
			try
				for t in hours
					for r in res.unit
						val = -1*CM_ΔP_R_dn_opt[t, r]
						fuel =  res.fuel[r]
						push!(output_CM_ΔP_R_dn_opt, [t r fuel val res.node[r]])
					end
				end
				CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_R_dn.csv"), output_CM_ΔP_R_dn_opt, append = true)
			catch e
				println("CM_ΔP_R_dn_opt does not exist.")
			end

			### Model parameters
			output_CM_info = DataFrame(timeperiod=String[],
											   termination_status=String[],
											   primal_status=String[],
											   dual_status=String[],
											   objective_value=Float64[],
											   solve_time=Float64[])

			push!(output_CM_info, [string(hours) string(termination_status(CM_mod)) string(primal_status(CM_mod)) string(dual_status(CM_mod)) objective_value(CM_mod) solve_time(CM_mod)])
			CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_info.csv"), output_CM_info, append = true)

			### Power flow results
			output_CM_P_flow_opt = DataFrame(time=Int64[],
											 line=Int64[],
											 from=Int64[],
											 to=Int64[],
											 pmax=Float64[],
											 pflow=Float64[])

			try
			   for t in hours
				   for l in lines.id
					   from = parse(Int, string(lines.from[l]))
					   to = parse(Int, string(lines.to[l]))
					   pmax = lines.pmax[l]
					   pflow = CM_P_flow_opt[t, l]
					   push!(output_CM_P_flow_opt, [t l from to pmax pflow])
				   end
			   end
			   CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_flow.csv"), output_CM_P_flow_opt, append = true)
		   catch e
			   println("CM_P_flow_opt does not exist.")
		   end

		   # Lost generation
		   output_CM_P_gen_lost = DataFrame(time=Int64[],node=Int64[],output=Float64[])
		   try
			   for t in hours
				   for j in nodes.id
					   output = CM_P_gen_lost_opt[t, j]
					   push!(output_CM_P_gen_lost, [t, j, output])
				   end
			   end
			   CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_gen_lost.csv"), output_CM_P_gen_lost, append = true)
		   catch e
			   println("CM_P_gen_lost does not exist.")
		   end

		   # Lost load
		   output_CM_P_load_lost = DataFrame(time=Int64[],node=Int64[],output=Float64[])
		   try
			   for t in hours
				   for j in nodes.id
					   output = CM_P_load_lost_opt[t, j]
					   push!(output_CM_P_load_lost, [t, j, output])
				   end
			   end
			   CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_load_lost.csv"), output_CM_P_load_lost, append = true)
		   catch e
			   println("CM_P_load_lost does not exist.")
		   end

		   # Noral Price
		   output_CM_price = DataFrame(time=Int64[],node=Int64[],price=Float64[])
	   	    try
	   	    	for t in hours
					for j in nodes.id
		   	    		price =  CM_price[j, t]
		   	    		push!(output_CM_price, [t j price])
					end
	   	    	end
	   	    	CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_price.csv"), output_CM_price, append = true)
	   	    catch e
	   	    	println("CM_price does not exist.")
	   	    end

			output_CM_P_S_up_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
			try
				for t in hours
					for s in storages.unit
						val = CM_P_S_up_opt[t, s]
						fuel =  storages.fuel[s]
						push!(output_CM_P_S_up_opt, [t s fuel val storages.node[s]])
					end
				end
				CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "CM_P_S_up.csv"), output_CM_P_S_up_opt, append = true)
			catch e
				println("CM_P_S_up_opt does not exist.")
			end

			# Nodal exchange
			output_exchange = DataFrame(time=Int64[],node=Int64[],output=Float64[])
			try
				for t in hours
					for j in nodes.id
						push!(output_exchange, [t j nodes.exchange[j][t]])
					end
				end
				CSV.write(joinpath(outputpath, session*"__"*timeperiod_str, "exchange.csv"), output_exchange, append = true)
			catch e
				println("output_exchange does not exist.")
			end

	next!(prog)

end

println()
println("Calculation complete.")

### ALTER CODE ENDE

#= ERSTELLT DICTIONARY ZONE ZU GENERATORS - ORIGINALCODE
N = nodes.id
P = pp.unit
R = res.unit
S = storages.unit
Node2PP = Dict((i, []) for i in nodes.id)
for i in keys(pp.node)
	push!(Node2PP[pp.node[i]], i)
end
Node2Res = Dict((i, []) for i in nodes.id)
for i in keys(res.node)
	push!(Node2Res[res.node[i]], i)
end
Node2Stor = Dict((i, []) for i in nodes.id)
for i in keys(storages.node)
	push!(Node2Stor[storages.node[i]], i)
end
Node2Res
Node2PP
Node2Stor
=#

#= TESTER
C_RD_N
C_test = copy(C_RD_N)
C_RD_N = nothing
C_RD_N = DataFrame(time=Int64[], node=Int64[], CostPP=Float64[], CostRes=Float64[], CostStor=Float64[], CostLostLoad=Float64[], CostLostGen=Float64[], Total=Float64[])
hours = dayslicer(1:24)[1]
for t in hours
	for n in N
		if isempty(Node2PP[n])
			costpp = 0
		else
			costpp = sum((pp.mc[g][1]*CM_ΔP_up_opt[t, g] + (ED_price[t] - pp.mc[g][]) * CM_ΔP_dn_opt[t, g]) for g in Node2PP[n])
		end
		if isempty(Node2Res[n])
			costres = 0
		else
			costres = sum(((ED_price[t]) * CM_ΔP_R_dn_opt[t, res]) for res in Node2Res[n])
		end
		if isempty(Node2Stor[n])
			coststor = 0
		else
			coststor = sum(((ED_price[t]/storages.efficiency[stor]) * CM_P_S_up_opt[t, stor]) for stor in Node2Stor[n])
		end

		costlostload = 1000*CM_P_load_lost_opt[t, n]
		costlostgen = 1000*CM_P_gen_lost_opt[t, n] # ÜBERALL NULL!
		totalcost = costpp + costres + coststor + costlostload + costlostgen

		push!(C_RD_N, [t n costpp costres coststor costlostload costlostgen totalcost])
	end
end
C_RD_N

C_RD_N == C_test

### hier drunter nur testzeilen

C_RD_N == C_test
C_asd = copy(C_test)
C_asd == C_test
C_test
CM_P_load_lost_opt[21, 272]
C_RD_N
N
C_RD_N[20*451+269, :CostPP]
C_RD_N[20*451+269, :CostLostLoad]
C_RD_N[20*451+269, :CostPP]+C_RD_N[20*451+269, :CostLostLoad]
costpp
pp.node[:BNA0015]

sum(CM_P_gen_lost_opt)
CM_P_load_lost_opt
	1000 * sum(P_load_lost[t, j] * Sbase for j in J, t in T) +
	1000 * sum(P_gen_lost[t, j] * Sbase for j in J, t in T)

coststor = sum(((ED_price[1]/storages.efficiency[stor]) * CM_P_S_up_opt[1, stor]) for stor in Node2Stor[331])
coststor = sum(((ED_price[1]/storages.efficiency[:BNA0467]) * CM_P_S_up_opt[1, :BNA0467]))

CM_P_S_up_opt
storages.node[:BNA0467]
storages.unit[12]
Node2Stor
sum((price[t]/storages.efficiency[s]) * P_S_up[t, s] * Sbase for s in S, t in T) +

costpp = sum((pp.mc[g][1]*CM_ΔP_up_opt[1, g] + (ED_price[1] - pp.mc[g][]) * CM_ΔP_dn_opt[1, g]) for g in Node2PP[46])
(pp.mc[:BNA0015][1]*CM_ΔP_up_opt[1, :BNA0015] + (ED_price[1] - pp.mc[:BNA0015][]) * CM_ΔP_dn_opt[1, :BNA0015])
CM_ΔP_up_opt[1, :BNA0015]* pp.mc[:BNA0015][1]
costpp

Node2PP[46]

pp.node[:BNA0015]
CM_ΔP_up_opt[1, :BNA0015] # ΔP_up von time und powerplant
CM_ΔP_dn_opt[1, :BNA0015] # ΔP_dn von time und powerplant

pp.mc[:BNA0015][] # marginalcost as single value, input powerplant
ED_price[1] # marketprice at hour 1

pp.mc
sum(powerplants.mc[g][1] * (ΔP_up[t, g]) * Sbase + (price[t] - powerplants.mc[g][1]) * ΔP_dn[t, g] * Sbase for g in G, t in T)
C_RD_N[:time=hours]

ED_price[1]

pp.node[:BNA0253]
keys(pp.node)
pp.unit
res.node
pp.node
storages.node

JuMP.objective_function(CM_mod)
JuMP.objective_value(CM_mod)
JuMP.all_variables(CM_mod)

a = JuMP.value.(CM_mod[:ΔP_up]) * 1000

=#


text:
export spatial_dvv

function spatial_dvv(InputDict::Dict)

    basefiname 			= InputDict["basefiname"]
    triangulationfiname = InputDict["triangulationfiname"]
    outputformat 		= InputDict["outputformat"]
    plotfig 			= InputDict["plotfig"]
    figname 			= InputDict["figname"]

    # load data
    t 	= jldopen(triangulationfiname)
    tris = t["tri"]

    # load existin dvv histories
    basefidir = basefiname[1:end-length(basename(basefiname))]
    dvvpath = ls(basefidir)
    dvvstalist = Array{String,1}(undef, 0)

    for path in dvvpath
        bn = split(basename(path)[length(basename(basefiname))+1:end], "-")
        try
            s1 = bn[1]
            s2 = bn[2]
            comp = bn[3][1:2]

            # TODO: 3 components cross-correlation.
            dvvstalist = vcat(dvvstalist, (s1, s2))
        catch
            #println(path)
            # remove this pair from list
            passid = findfirst(x -> x == path, dvvpath)
            dvvpath = dvvpath[setdiff(1:end, passid), :]
        end
    end

    # loop triangles

    loc_and_dvv = []
    # synchronize time to take mean

    stday = Inf
    etday = 0

    for tri in tris
        centerlonlat = tri["centerxy"]
        edgeinfo = tri["edgeinfo"]

        numofedge = 0
        dvvsum = []
        Tsum = []

		#println("tri $tri")

        for edge in edgeinfo
            stn1 = edge[1]
            stn2 = edge[2]

            #println((stn1, stn2))
            # search corresponding dvv history
            dvvid = findfirst(x -> x == (stn1, stn2), dvvstalist)
            dvvidrev = findfirst(x -> x == (stn2, stn1), dvvstalist)

            if !isnothing(dvvid)
                fipath = dvvpath[dvvid]
            elseif !isnothing(dvvidrev) # reverse
                fipath = dvvpath[dvvidrev]
            else
                #println("debug: no data on this edge.")
            end

            try
                fiedge = jldopen(fipath, "r")

                push!(Tsum, fiedge["T"])
                push!(dvvsum, fiedge["dvv"])

                close(fiedge)
            catch
                println("debug: load error edge loading.")
            end
        end

        avgdvv = []
        avgT = []

        println("start average")

        for tsumid = 1:length(Tsum)
            for ii = 1:length(Tsum[tsumid])
				# try to average at this time
                testtavg = Tsum[tsumid][ii]
                dvvtemp = []

                if isnothing(findfirst(x -> x == testtavg, avgT))
					# aveeraging this timestep

                    edgecount = 0
                    for jj = 1:length(Tsum)
                        oid = findfirst(x -> x == testtavg, Tsum[jj])

                        if !isnothing(oid) && !isnan(dvvsum[jj][oid])
							# avoid summing up NaN value
							# debuged: index is jj on dvvsum
                            push!(dvvtemp, dvvsum[jj][oid])
                            edgecount += 1
                        else

                        end
                    end

                    #println(dvvtemp)

                    if !isempty(dvvtemp) && edgecount > 0
                        push!(avgT, testtavg)
                        push!(avgdvv, mean(dvvtemp))
                    else
                        push!(avgT, testtavg)
                        push!(avgdvv, 0)
                    end
                end
            end
        end

        # if !any(iszero.(avgdvv))
        #     println("success")
        # else
        #     println(iszero.(avgdvv))
        #     println("avgdvv NaN")
        # end

        sortid = sortperm(avgT)
        avgT = avgT[sortid]
        avgdvv = avgdvv[sortid]
        avgT = convert(Array{Int64,1}, avgT)
        avgvv = convert(Array{Float64,1}, avgdvv)
        # println(sortid)

        # smoothing dv/v
        if InputDict["smoothing"]
            trace1 = PlotlyJS.scatter(x = timestamp.(avgT), y = avgdvv)
            SeisNoise.smooth!(avgdvv; half_win = 5)
            #println(avgdvv)
            trace2 = PlotlyJS.scatter(
                x = timestamp.(avgT),
                y = avgdvv,
                label = "after smooth"
            )
            p = PlotlyJS.plot([trace1, trace2])
            display(p)
            readline()
        end

        if !isempty(avgT)
            stday_test = avgT[1]
            etday_test = avgT[end]
            if stday_test < stday
                stday = stday_test
            end
            if etday_test > etday
                etday = etday_test
            end
        end

        push!(loc_and_dvv, (centerlonlat, edgeinfo, avgT, avgdvv))

    end

    dvvdata = Dict(
        "startday" => timestamp(stday),
        "endday" => timestamp(etday),
        "loc_and_dvv" => loc_and_dvv
    )

    # save data
    tstamp = timestamp(stday)

    if outputformat == "hdf5"

        tstamp = timestamp(stday)
        foname = basefidir * "/../spatialdvv_$(tstamp).h5"

        if ispath(foname)
            rm(foname)
        end

        h5open(foname, "w") do file

            gorigin = "spatiotemporal_dvv" #original group name

            write(file, gorigin * "/startday", dvvdata["startday"])
            write(file, gorigin * "/endday", dvvdata["endday"])

            loc_and_dvv = dvvdata["loc_and_dvv"]

            #create group on each element
            for i = 1:length(loc_and_dvv)
                centerlonlat = loc_and_dvv[i][1]
                edgeinfo = loc_and_dvv[i][2]
                T = loc_and_dvv[i][3]
                dvvsum = loc_and_dvv[i][4]

                groupname = join([gorigin, "ele$i"], "/")

                write(file, groupname * "/dvv", Array{Float64,1}(dvvsum))
                write(file, groupname * "/centerlon", centerlonlat[1])
                write(file, groupname * "/centerlat", centerlonlat[2])
                write(file, groupname * "/T", Array{Int64,1}(T))

                for edgeid = 1:length(edgeinfo)
                    edge = edgeinfo[edgeid]
                    sta1 = edge[1]
                    sta2 = edge[2]

                    edgegroupname = join([groupname, "edges$(edgeid)"], "/")

                    write(file, edgegroupname * "/sta1", sta1)
                    write(file, edgegroupname * "/sta1edge", edge[3])
                    write(file, edgegroupname * "/sta2", sta2)
                    write(file, edgegroupname * "/sta2edge", edge[4])

                end
            end
        end

    elseif outputformat == "jld2"

        foname = basefidir * "/../spatialdvv_$(tstamp).jld2"
        if ispath(foname)
            rm(foname)
        end

        outfile = jldopen(foname, "w")
        outfile["dvvdata"] = dvvdata
        close(outfile)

    else
        println("output format not exist. No output has been done.")
    end

    if plotfig
        fodir = "./fig"
        if !ispath(fodir)
            mkpath(fodir)
        end

        println("Plot figure is is working in progress. Please use Matlab or Python to plot.")
        #
        # mul = 400
        # height= abs(diff([studyarea[1], studyarea[3]])[1]) * mul
        # width = abs(diff([studyarea[2], studyarea[4]])[1]) * mul
        # layout = Layout(;width=width, height=height, title="Triangulation of station network",
        #                xaxis_range=[studyarea[4], studyarea[2]],
        #                yaxis_range=[studyarea[3], studyarea[1]],
        #                xaxis_title="lon",
        #                yaxis_title="lat",
        #                showlegend=false)
        #
        # p = PlotlyJS.plot(traceall, layout)
        # PlotlyJS.savefig(p, fodir*"/"*InputDict["figname"])
        #display(p)
        #println("type return...")
        #readline()
    end

    println("plot dvv has been successfully done.")
    return nothing

end


text:
using ModelingToolkit, StaticArrays, LinearAlgebra
using DiffEqBase, StochasticDiffEq, SparseArrays
using Test

# Define some variables
@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t

eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]

noiseeqs = [0.1*x,
            0.1*y,
            0.1*z]

de = SDESystem(eqs,noiseeqs,t,[x,y,z],[σ,ρ,β])
f = eval(generate_diffusion_function(de)[1])
@test f(ones(3),rand(3),nothing) == 0.1ones(3)

f = SDEFunction(de)
prob = SDEProblem(SDEFunction(de),f.g,[1.0,0.0,0.0],(0.0,100.0),(10.0,26.0,2.33))
sol = solve(prob,SRIW1())

noiseeqs_nd = [0.01*x 0.01*x*y 0.02*x*z
               σ      0.01*y   0.02*x*z
               ρ      β        0.01*z  ]
de = SDESystem(eqs,noiseeqs_nd,t,[x,y,z],[σ,ρ,β])
f = eval(generate_diffusion_function(de)[1])
@test f([1,2,3.0],[0.1,0.2,0.3],nothing) == [0.01*1   0.01*1*2   0.02*1*3
                                             0.1      0.01*2     0.02*1*3
                                             0.2      0.3        0.01*3  ]

f = eval(generate_diffusion_function(de)[2])
du = ones(3,3)
f(du,[1,2,3.0],[0.1,0.2,0.3],nothing)
@test du == [0.01*1   0.01*1*2   0.02*1*3
             0.1      0.01*2     0.02*1*3
             0.2      0.3        0.01*3  ]

f = SDEFunction(de)
prob = SDEProblem(SDEFunction(de),f.g,[1.0,0.0,0.0],(0.0,100.0),(10.0,26.0,2.33),
                  noise_rate_prototype = zeros(3,3))
sol = solve(prob,EM(),dt=0.001)

u0map = [
    x => 1.0,
    y => 0.0,
    z => 0.0
]

parammap = [
    σ => 10.0,
    β => 26.0,
    ρ => 2.33
]

prob = SDEProblem(de,u0map,(0.0,100.0),parammap)
@test size(prob.noise_rate_prototype) == (3,3)
@test prob.noise_rate_prototype isa Matrix
sol = solve(prob,EM(),dt=0.001)

prob = SDEProblem(de,u0map,(0.0,100.0),parammap,sparsenoise=true)
@test size(prob.noise_rate_prototype) == (3,3)
@test prob.noise_rate_prototype isa SparseMatrixCSC
sol = solve(prob,EM(),dt=0.001)


text:
using HttpServer

import Base.ismatch

"""
A `Route` represents an endpoint on the Fuji webserver. It currently has three
fields:

* `action`: The block of code that will be run when the endpoint is hit. The
    request and response are both passed to it
* `endpoint`: The endpoint that this route serves
* `methods`: The HTTP methods that this route should serve (e.g. GET, POST)
"""
immutable Route
    action::Function
    endpoint::String
    methods::Array{String,1}
end

"""
    ismatch(request, route)

Returns true if this route 'matches' this request, false if it doesn't. Note
that this is a request from HttpCommon, not from Fuji.
"""
function ismatch(request::Request, route::Route)
    if in(request.method, route.methods)
        named_param_regex = r":[A-z0-9\-_]+"
        splat_param_regex = r"\*"
        component_string = "([A-z0-9\\-_]+)"

        # Replace named parameters with regex string to see if it matches the request
        endpoint = replace(route.endpoint, named_param_regex, component_string)

        # Replace splat parameters with regex string to see if it matches the request
        endpoint = replace(endpoint, splat_param_regex, component_string)

        endpoint_regex = Regex(string("^", endpoint, "\$"))
        return ismatch(endpoint_regex, split(request.resource, "?")[1])
    end

    false
end


text:
# Data topmute
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: December 2017
#

export marineTopmute,opMarineTopmute

function marineTopmute(Dobs::joData, muteStart::Integer; mute=Array{Any}(3))
	# Data topmute for end-on spread marine streamer data
	Din = deepcopy(Dobs)

	# Design mute window
	j=1
	x0 = 1f0
	xend = length(Din.geometry.xloc[j])
	nt = Din.geometry.nt[j]
	nrec = length(Din.geometry.xloc[j])
	drec = abs(Din.geometry.xloc[j][1] - Din.geometry.xloc[j][2])
	offsetDirectWave = 1.5f0*Din.geometry.t[j]
	idxOffset = Int(round(offsetDirectWave/drec))
	dx = round(idxOffset - idxOffset/10f0)

	if j==1 && ~isassigned(mute)
		z0 = muteStart - Int(round(muteStart/10))
		slope = 0.95f0*(nt - z0)/dx
		mute[1] = x0
		mute[2] = z0
		mute[3] = slope
	elseif j==1 && isassigned(mute)
		x0 = mute[1]
		z0 = mute[2]
		slope = mute[3]
	end

	mask = ones(Float32,nt,nrec)
	mask[1:z0,:]=0f0

	# Linear mute
	if (nrec-x0 < dx)
		x = nrec
		zIntercept = Int(round(z0+slope*(x-x0)))
		zax = z0+1:1:zIntercept
	else
		x = x0+dx
		zax = z0+1:1:nt
	end
	if length(zax) > 1
		xax = Array{Int}(round.(linspace(x0,x,length(zax))))
	else
		xax = Int(round(x0))
	end
	for k=1:length(zax)
		mask[zax[k],xax[k]:end] = 0f0
	end

	for j=1:Din.nsrc
		Din.data[j] = Din.data[j].*mask
	end
	return Din
end

function opMarineTopmute(muteStart,geometry;params=Array{Any}(3))
# JOLI wrapper for the linear depth scaling function
	nsrc = length(geometry.xloc)
	N = 0
	for j=1:nsrc
		N += geometry.nt[j]*length(geometry.xloc[j])
	end
	D = joLinearFunctionFwdT(N,N,
							 v -> marineTopmute(v,muteStart;mute=params),
 							 w -> marineTopmute(w,muteStart;mute=params),
							 Float32,Float32,name="Data mute")
	return D
end

max_stars_repo_name: SINBADconsortium/opesciSLIM

max_issues_repo_name: SINBADconsortium/opesciSLIM

max_forks_repo_name: SINBADconsortium/opesciSLIM


text:
### Implement Laplace transform along a hyperbola contour ###

# 1.
# Weideman, J., and L. Trefethen.
# Parabolic and hyperbolic contours for computing the Bromwich integral.
# Mathematics of Computation 76.259 (2007): 1341-1356.

# 2.
# Liemert, André, and Alwin Kienle.
# Application of the Laplace transform in time-domain optical spectroscopy and imaging.
# Journal of biomedical optics 20.11 (2015): 110502.

### adaptive contour for each time point (can't precompute f(s) or sk)

# parameterized hyperbola contour
# parameters are given by ref [2] which are given in increased precision from calculations given in [1]
function s(θ, N, t)
    μ = 4.492075287 * N / t
    ϕ = 1.172104229
    return μ + im * μ * sinh(θ + im * ϕ)
end
# derivitive of hyperbola contour
function ds(θ, N, t)
    μ = 4.492075287 * N / t
    ϕ = 1.172104229
    return im * μ * cosh(θ + im * ϕ)
end

# compute the laplace transform along a hyperbola contour fixed time point
"""
    hyperbola(f::Function, t::AbstractFloat; N::Int = 16)
Evaluate the inverse Laplace transform of `f` at the point `t` by approximating the Bromwhich integral with a hyperbola contour.
A real valued time-domain signal is assumed with the integral being calculated by a computable series applying the midpoint rule.
The contour is defined by several parameters. μ controls the extreme nodes with a larger parameter moving the outlier nodes further into left half plane.
h is the uniform node spacing with N being the number of nodes on the contour. The parameters depend on both t and N such that the function must be computed N times for each t.
This is very inefficient especially if f is computational expensive, so this approach is recommended when only a few values of t are needed.
However, this method will provide more accurate results at the cost of computational time compared to fixed contour approaches.
The number of nodes, N, defaults to 16. The optimal N is case dependent and arbitrarily increasing N will not result in more accurate results. Smaller values of N will be faster.
Bromwhich contour approaches should only be applied for t > 0. We are aiming to optimize the convergence as N -> inf so our parameters are ~ N / t.
In other words we want the Bromwhich integral to converge rapidly, and this can be done by starting and ending our integration path in the left hand plane causing exp(st) to decay.
If t = 0, the exponential factor does not cause good convergence. If you want to evaluate for t = 0, you should try a small positive value like 1e-30.
# Example
```jldoctest
julia> InverseLaplace.hyperbola(s -> 1/(s + 1), 2.0)
0.13533528323665164
"""
function hyperbola(f::Function, t::AbstractFloat; N = 16)
    a =  zero(Complex{eltype(t)})
    N = convert(eltype(t), N)
    h = 1.081792140 / N
    for k in 0:N-1
        sk = s((k + 1/2) * h, N, t)
        dsk = ds((k + 1/2) * h, N, t)
        a += f(sk) * exp(sk * t) * dsk
    end
    return imag(a) * h / pi
end

# loop through an entire time array with multithreads
function hyperbola(f::Function, t::AbstractArray; N = 16)
    out = similar(t)
    Threads.@threads for ind in eachindex(t)
        out[ind] = hyperbola(f, t[ind], N = N)
    end
    return out
end

#### fixed integration path (can precompute) ####

# get the fixed integration components
function hyper_coef(N, t::AbstractArray; ϕ = 1.09)
    A = acosh(((π - 2 * ϕ) * t[end] / t[1] + 4 * ϕ - π) / ((4 * ϕ - π) * sin(ϕ)))
    μ = (4 * π * ϕ - π^2) * N / t[end] / A
    h = A / N
    return μ, h
end

# compute the fixed hyperbola contour
s_fixed(θ, μ; ϕ = 1.09) = μ + im * μ * sinh(θ + im * ϕ)
ds_fixed(θ, μ; ϕ = 1.09) = im * μ * cosh(θ + im * ϕ)

# compute the function values over the fixed contour nodes
function fixed_sk(f::Function, N, t::AbstractArray)
    μ, h = hyper_coef(N, t)
    a = zeros(Complex{eltype(h)}, Int(N))
    sk = similar(a)
    Threads.@threads for k in 0:Int(N)-1
        sk[k+1] = s_fixed((k + 1/2) * h, μ)
        dsk = ds_fixed((k + 1/2) * h, μ)
        a[k+1] = f(sk[k + 1]) * dsk
    end
    return a, sk, h
end

function hyper_fixed_points(a, sk, h, t::AbstractFloat)
    b = zero(eltype(a))
    for ind in eachindex(sk)
        b += a[ind] * exp(sk[ind] * t)
    end
    return imag(b) * h / π
end

"""
    hyper_fixed(f::Function, t::AbstractArray; N = 24)
Evaluate the inverse Laplace transform of `f` at an array of `t` by approximating the Bromwhich integral with a fixed hyperbola contour.
A real valued time-domain signal is assumed with the integral being calculated by a computable series applying the midpoint rule.
In contrast to `hyperbola` where the contour is dependent on 't', here we fix the contour for the entire array of `t` values.
This function should be considered when F(s) is complex and costly to compute where we want to minimize the number of evaulations of F(s).
This operates under the assumption that f(t) is needed for many values of `t` over some interval `t ∈ (t1:t2)`. The parameters are now time independent.
The angle value of the contour defaults to `ϕ = 1.09`. `ϕ` and the number of nodes `N` can be adjusted for accuracy.
This approach will in general be less accurate than `hyperbola` as we are using the same contour over an array of time values, instead of an optimized contour at each time.
Bromwhich contour approaches should only be applied for t > 0.
# Example
```jldoctest
julia> InverseLaplace.hyper_fixed(s -> 1/(s + 1), 2.0:3.0)
2-element Vector{Float64}:
 0.1353352832366098
 0.04978706836781721
"""
function hyper_fixed(f::Function, t::AbstractArray; N = 24)
    N = convert(eltype(t), N)
    a, sk, h = fixed_sk(f, N, t)
    out = zeros(eltype(h), length(t))
    Threads.@threads for ind in eachindex(t)
        out[ind] = hyper_fixed_points(a, sk, h, t[ind])
    end
    return out
end

#=

### Implement Laplace transform along improved Talbot contour ###

## adaptive contour for each time point (can't precompute f(ω) or sk)

function z(θ, N, t; α = 0.6407, c = 1.3580)
    B = c*(sin(α*pi))^2 / ( 2*α*c^2 * (sin(α*pi))^2 - pi*sin(2*α*pi)*(sinh(α*c))^2)
    σ = 2*α*c^2 * B
    μ = 2*(sinh(α*c))^2 * B
    ν = (sinh(2*α*c) - 2*α*c)*B

    # return both z(θ) and z'(θ)
    return N/t * (-σ + μ*θ*cot(α*θ) + ν*im*θ), N/t * (ν*im  + μ*(cot(α*θ) - α*θ*(csc(α*θ))^2))
end

# compute the laplace transform along a hyperbola contour fixed time point
function talbot_improved(f::Function, N, t::AbstractFloat)
    a =  zero(Complex{eltype(N)})
    h = 2*π/N
    for k in 1:N
        zk, dzk = z(-π + (k - 1/2)*h, N, t)
        a += f(zk)*exp(zk*t)*dzk
    end
    return a/N/im
end

function talbot_improved(f::Function, N, t::AbstractFloat)
    a =  zero(Complex{eltype(N)})
    h = 2*π/N
    for k in 1:N/2 # for real-valued f(t) need only to consider half of
        zk, dzk = z(-π + (k - 1/2)*h, N, t)
        a += f(zk)*exp(zk*t)*dzk
        println(-π + (k - 1/2)*h)
    end
    return 2*a/N/im
end

function talbot_improved(f::Function, N, t::AbstractFloat)
    a =  zero(Complex{eltype(N)})
    h = π/N
    for k in 1:N # for real-valued f(t) need only to consider half of evaluations
        zk, dzk = z(-π + (k - 1/2)*2*h, N, t)
        a += f(zk)*exp(zk*t)*dzk
    end
    return real(a/N/im)
end

=#


text:
using ZXCalculus, Graphs, Multigraphs, SparseArrays
using Documenter
using Test

@testset "scalar.jl" begin
    include("scalar.jl")
end

@testset "abstract_zx_diagram.jl" begin
    include("abstract_zx_diagram.jl")
end

@testset "zx_diagram.jl" begin
    include("zx_diagram.jl")
end

@testset "rules.jl" begin
    include("rules.jl")
end

@testset "zx_graph.jl" begin
    include("zx_graph.jl")
end

@testset "circuit_extraction.jl" begin
    include("circuit_extraction.jl")
end

@testset "phase_teleportation.jl" begin
    include("phase_teleportation.jl")
end

@testset "ir.jl" begin
    include("ir.jl")
end

@testset "phase.jl" begin
    include("phase.jl")
end

@testset "simplify.jl" begin
    include("simplify.jl")
end

@testset "ancilla_extraction.jl" begin
    include("ancilla_extraction.jl")
    include("challenge.jl")
end

doctest(ZXCalculus)


text:
function binary_op(op, W, @nospecialize(T))
    ty = LLVM_TYPES[T]
    if isone(W)
        V = T
    else
        ty = "<$W x $ty>"
        V = NTuple{W,VecElement{T}}
    end
    instrs = "%res = $op $ty %0, %1\nret $ty %res"
    call = :($LLVMCALL($instrs, $V, Tuple{$V,$V}, data(v1), data(v2)))
    W > 1 && (call = Expr(:call, :Vec, call))
    Expr(:block, Expr(:meta, :inline), call)
end

# Integer
for (op,f) ∈ [("add",:+),("sub",:-),("mul",:*),("shl",:<<)]
  ff = Symbol('v', op)
  fnsw = Symbol(ff,"_nsw")
  fnuw = Symbol(ff,"_nuw")
  fnw = Symbol(ff,"_nw")
  ff_fast = Symbol(ff, :_fast)
  @eval begin
    # @inline $ff(a,b) = $ff_fast(a,b)
    @inline $ff(a::T,b::T) where {T<:Union{FloatingTypes,IntegerTypesHW,AbstractSIMD}} = $ff_fast(a,b)
    # @generated $ff_fast(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:IntegerTypesHW} = binary_op($op * (T <: Signed ? " nsw" : " nuw"), W, T)
    @generated $ff_fast(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:IntegerTypesHW} = binary_op($op, W, T)
    @generated $fnsw(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:IntegerTypesHW} = binary_op($(op * " nsw"), W, T)
    @generated $fnuw(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:IntegerTypesHW} = binary_op($(op * " nuw"), W, T)
    @generated $fnw(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:IntegerTypesHW} = binary_op($(op * " nsw nuw"), W, T)
    # @generated Base.$f(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:IntegerTypesHW} = binary_op($op, W, T)
    @inline Base.$f(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:IntegerTypesHW} = $ff_fast(v1, v2)
    # @generated $ff(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:IntegerTypesHW} = binary_op($op, W, T)
    @inline $ff_fast(x, y) = $f(x,y)

    # @generated $ff_fast(v1::T, v2::T) where {T<:IntegerTypesHW} = binary_op($op * (T <: Signed ? " nsw" : " nuw"), 1, T)
    # @generated $ff_fast(v1::T, v2::T) where {T<:IntegerTypesHW} = binary_op($op, 1, T)
    @generated $fnsw(v1::T, v2::T) where {T<:IntegerTypesHW} = binary_op($(op * " nsw"), 1, T)
    @generated $fnuw(v1::T, v2::T) where {T<:IntegerTypesHW} = binary_op($(op * " nuw"), 1, T)
    @generated $fnw(v1::T, v2::T) where {T<:IntegerTypesHW} = binary_op($(op * " nsw nuw"), 1, T)
  end
end
@inline vadd_fast(v1::T, v2::T) where {T<:IntegerTypesHW} = Base.add_int(v1,v2)
@inline vsub_fast(v1::T, v2::T) where {T<:IntegerTypesHW} = Base.sub_int(v1,v2)
@inline vmul_fast(v1::T, v2::T) where {T<:IntegerTypesHW} = Base.mul_int(v1,v2)
@inline vshl_fast(v1::T, v2::T) where {T<:IntegerTypesHW} = Base.shl_int(v1,v2)
for (op,f) ∈ [("div",:÷),("rem",:%)]
  ff = Symbol('v', op); #_ff = Symbol(:_, ff)
  sbf = Symbol('s', op, :_int)
  ubf = Symbol('u', op, :_int)
    @eval begin
      @generated $ff(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:Integer} = binary_op((T <: Signed ? 's' : 'u') * $op, W, T)
      @inline $ff(a::I, b::I) where {I<:SignedHW} = Base.$sbf(a, b)
      @inline $ff(a::I, b::I) where {I<:UnsignedHW} = Base.$ubf(a, b)
      @inline $ff(a::Int128, b::Int128) = Base.$sbf(a, b)
      @inline $ff(a::UInt128, b::UInt128) = Base.$ubf(a, b)
        # @generated $_ff(v1::T, v2::T) where {T<:Integer} = binary_op((T <: Signed ? 's' : 'u') * $op, 1, T)
        # @inline $ff(v1::T, v2::T) where {T<:IntegerTypesHW} = $_ff(v1, v2)
    end
end
# for (op,f) ∈ [("div",:÷),("rem",:%)]
#   @eval begin
#     @generated Base.$f(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:Integer} = binary_op((T <: Signed ? 's' : 'u') * $op, W, T)
#     @generated Base.$f(v1::T, v2::T) where {T<:IntegerTypesHW} = binary_op((T <: Signed ? 's' : 'u') * $op, 1, T)
#   end
# end
@inline vcld(x, y) = vadd(vdiv(vsub(x, one(x)), y), one(x))
@inline function vdivrem(x, y)
    d = vdiv(x, y)
    r = vsub(x, vmul(d, y))
    d, r
end
for (op,sub) ∈ [
    ("ashr",:SignedHW),
    ("lshr",:UnsignedHW),
    ("lshr",:IntegerTypesHW),
    ("and",:IntegerTypesHW),
    ("or",:IntegerTypesHW),
    ("xor",:IntegerTypesHW)
]
    ff = sub === :UnsignedHW ? :vashr : Symbol('v', op)
    @eval begin
        @generated $ff(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:$sub}  = binary_op($op, W, T)
        @generated $ff(v1::T, v2::T) where {T<:$sub}  = binary_op($op, 1, T)
    end
end

for (op,f) ∈ [("fadd",:vadd),("fsub",:vsub),("fmul",:vmul),("fdiv",:vfdiv)]#,("frem",:vrem)]
  ff = Symbol(f, :_fast)
  fop_fast = f === :vfdiv ? "fdiv fast" : op * ' ' * fast_flags(true)
  fop_contract = op * ' ' * fast_flags(false)
  @eval begin
    @generated  $f(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:Union{Float32,Float64}} = binary_op($fop_contract, W, T)
    @generated $ff(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:Union{Float32,Float64}} = binary_op($fop_fast, W, T)
    @inline $f(v1::Vec{W,Float16}, v2::Vec{W,Float16}) where {W} = $f(convert(Vec{W,Float32}, v1), convert(Vec{W,Float32}, v2))
    @inline $ff(v1::Vec{W,Float16}, v2::Vec{W,Float16}) where {W} = $ff(convert(Vec{W,Float32}, v1), convert(Vec{W,Float32}, v2))
  end
end
@inline vsub(a::T,b::T) where {T<:Union{Float32,Float64}} = Base.sub_float(a,b)
@inline vadd(a::T,b::T) where {T<:Union{Float32,Float64}} = Base.add_float(a,b)
@inline vmul(a::T,b::T) where {T<:Union{Float32,Float64}} = Base.mul_float(a,b)
@inline vsub_fast(a::T,b::T) where {T<:Union{Float32,Float64}} = Base.sub_float_fast(a,b)
@inline vadd_fast(a::T,b::T) where {T<:Union{Float32,Float64}} = Base.add_float_fast(a,b)
@inline vmul_fast(a::T,b::T) where {T<:Union{Float32,Float64}} = Base.mul_float_fast(a,b)

@inline vdiv(v1::AbstractSIMD{W,T}, v2::AbstractSIMD{W,T}) where {W,T<:FloatingTypes} = trunc(vfdiv_fast(v1, v2))
@inline vdiv_fast(v1::AbstractSIMD{W,T}, v2::AbstractSIMD{W,T}) where {W,T<:FloatingTypes} = trunc(vfdiv_fast(v1, v2))
@inline vdiv_fast(v1::T, v2::T) where {T<:FloatingTypes} = trunc(Base.FastMath.div_float_fast(v1, v2))
@inline vdiv_fast(v1::T, v2::T) where {T <: Number} = v1 ÷ v2
@inline vdiv(v1::T, v2::T) where {T <: Number} = v1 ÷ v2
@inline vdiv(v1::T, v2::T) where {T <: FloatingTypes} = vdiv_fast(v1, v2)
@inline vrem(a,b) = vfnmadd(vdiv_fast(a, b), b, a)
@inline vrem_fast(a,b) = vfnmadd(vdiv_fast(a, b), b, a)
# @inline vdiv_fast(v1::AbstractSIMD{W,T}, v2::AbstractSIMD{W,T}) where {W,T<:IntegerTypesHW} = trunc(T, vfloat_fast(v1) / vfloat_fast(v2))
@inline vdiv_fast(v1::AbstractSIMD{W,T}, v2::AbstractSIMD{W,T}) where {W,T<:IntegerTypesHW} = trunc(T, vfloat(v1) / vfloat(v2))
@inline function vdiv_fast(v1, v2)
    v3, v4 = promote_div(v1, v2)
    vdiv_fast(v3, v4)
end
@inline vdiv_fast(v1::VecUnroll{N,1,T,T}, s::T) where {N,T<:SignedHW} = VecUnroll(fmap(Base.sdiv_int, data(v1), s))
@inline vdiv_fast(v1::VecUnroll{N,1,T,T}, s::T) where {N,T<:UnsignedHW} = VecUnroll(fmap(Base.udiv_int, data(v1), s))
@inline vdiv_fast(v1::VecUnroll{N,1,T,T}, v2::VecUnroll{N,1,T,T}) where {N,T<:SignedHW} = VecUnroll(fmap(Base.sdiv_int, data(v1), data(v2)))
@inline vdiv_fast(v1::VecUnroll{N,1,T,T}, v2::VecUnroll{N,1,T,T}) where {N,T<:UnsignedHW} = VecUnroll(fmap(Base.udiv_int, data(v1), data(v2)))

@inline vfdiv(a::AbstractSIMDVector{W}, b::AbstractSIMDVector{W}) where {W} = vfdiv(vfloat(a), vfloat(b))
# @inline vfdiv_fast(a::AbstractSIMDVector{W}, b::AbstractSIMDVector{W}) where {W} = vfdiv_fast(vfloat_fast(a), vfloat_fast(b))
@inline vfdiv_fast(a::AbstractSIMDVector{W}, b::AbstractSIMDVector{W}) where {W} = vfdiv_fast(vfloat(a), vfloat(b))
@inline vfdiv(a, b) = a / b
@inline vfdiv_fast(a, b) = Base.FastMath.div_fast(a, b)

for (f,op) ∈ ((:vand,:and_int),(:vor,:or_int),(:vxor,:xor_int))
  @eval @inline $f(b1::Bool, b2::Bool) = Base.$op(b1, b2)
end

for f ∈ [:vadd,:vsub,:vmul]
  for s ∈ [Symbol(""),:_fast,:_nsw,:_nuw,:_nw]
    fs = Symbol(f,s)
    @eval begin
      @inline function $fs(a::Union{FloatingTypes,IntegerTypesHW,AbstractSIMD}, b::Union{FloatingTypes,IntegerTypesHW,AbstractSIMD})
        c, d = promote(a, b)
        $fs(c, d)
      end
    end
  end
end
# @inline vsub(a::T, b::T) where {T<:Base.BitInteger} = Base.sub_int(a, b)
for (vf,bf) ∈ [
  (:vadd,:add_int),(:vsub,:sub_int),(:vmul,:mul_int),
  (:vadd_fast,:add_int),(:vsub_fast,:sub_int),(:vmul_fast,:mul_int),
  (:vadd_nsw,:add_int),(:vsub_nsw,:sub_int),(:vmul_nsw,:mul_int),
  (:vadd_nuw,:add_int),(:vsub_nuw,:sub_int),(:vmul_nuw,:mul_int),
  (:vadd_nw,:add_int),(:vsub_nw,:sub_int),(:vmul_nw,:mul_int),
]
  @eval begin
    @inline $vf(a::Int128, b::Int128) = Base.$bf(a, b)
    @inline $vf(a::UInt128, b::UInt128) = Base.$bf(a, b)
  end
end
# @inline vrem(a::Float32, b::Float32) = Base.rem_float_fast(a, b)
# @inline vrem(a::Float64, b::Float64) = Base.rem_float_fast(a, b)

@inline function Base.FastMath.add_fast(a::AbstractSIMD, b::AbstractSIMD, c::AbstractSIMD)
  Base.FastMath.add_fast(Base.FastMath.add_fast(a,b),c)
end
@inline function Base.FastMath.add_fast(a::T, b::T, c::T) where {T<:AbstractSIMD}
  Base.FastMath.add_fast(Base.FastMath.add_fast(a,b),c)
end

@inline function Base.FastMath.add_fast(a::AbstractSIMD, b::AbstractSIMD, c::AbstractSIMD, d::AbstractSIMD)
  x = Base.FastMath.add_fast(a,b)
  y = Base.FastMath.add_fast(c,d)
  Base.FastMath.add_fast(x, y)
end
@inline function Base.FastMath.add_fast(a::T, b::T, c::T, d::T) where {T<:AbstractSIMD}
  x = Base.FastMath.add_fast(a,b)
  y = Base.FastMath.add_fast(c,d)
  Base.FastMath.add_fast(x, y)
end

@inline function Base.FastMath.add_fast(a::AbstractSIMD, b::AbstractSIMD, c::AbstractSIMD, d::AbstractSIMD, e::AbstractSIMD, f::Vararg{Number,K}) where {K}
  x = Base.FastMath.add_fast(a,b)
  y = Base.FastMath.add_fast(c,d)
  Base.FastMath.add_fast(Base.FastMath.add_fast(x, y), e, f...)
end
@inline function Base.FastMath.add_fast(a::T, b::T, c::T, d::T, e::T, f::Vararg{T,K}) where {T<:AbstractSIMD,K}
  x = Base.FastMath.add_fast(a,b)
  y = Base.FastMath.add_fast(c,d)
  Base.FastMath.add_fast(Base.FastMath.add_fast(x, y), e, f...)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineImageProperties <: SwaggerModel
    plan::Any # spec type: Union{ Nothing, PurchasePlan } # spec name: plan
    osDiskImage::Any # spec type: Union{ Nothing, OSDiskImage } # spec name: osDiskImage
    dataDiskImages::Any # spec type: Union{ Nothing, Vector{DataDiskImage} } # spec name: dataDiskImages

    function VirtualMachineImageProperties(;plan=nothing, osDiskImage=nothing, dataDiskImages=nothing)
        o = new()
        validate_property(VirtualMachineImageProperties, Symbol("plan"), plan)
        setfield!(o, Symbol("plan"), plan)
        validate_property(VirtualMachineImageProperties, Symbol("osDiskImage"), osDiskImage)
        setfield!(o, Symbol("osDiskImage"), osDiskImage)
        validate_property(VirtualMachineImageProperties, Symbol("dataDiskImages"), dataDiskImages)
        setfield!(o, Symbol("dataDiskImages"), dataDiskImages)
        o
    end
end # type VirtualMachineImageProperties

const _property_map_VirtualMachineImageProperties = Dict{Symbol,Symbol}(Symbol("plan")=>Symbol("plan"), Symbol("osDiskImage")=>Symbol("osDiskImage"), Symbol("dataDiskImages")=>Symbol("dataDiskImages"))
const _property_types_VirtualMachineImageProperties = Dict{Symbol,String}(Symbol("plan")=>"PurchasePlan", Symbol("osDiskImage")=>"OSDiskImage", Symbol("dataDiskImages")=>"Vector{DataDiskImage}")
Base.propertynames(::Type{ VirtualMachineImageProperties }) = collect(keys(_property_map_VirtualMachineImageProperties))
Swagger.property_type(::Type{ VirtualMachineImageProperties }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualMachineImageProperties[name]))}
Swagger.field_name(::Type{ VirtualMachineImageProperties }, property_name::Symbol) =  _property_map_VirtualMachineImageProperties[property_name]

function check_required(o::VirtualMachineImageProperties)
    true
end

function validate_property(::Type{ VirtualMachineImageProperties }, name::Symbol, val)
end