text:
using DashBootstrapComponents

alerts = html_div([
    dbc_alert(
        [
            html_i(className = "bi bi-info-circle-fill me-2"),
            "An example info alert with an icon",
        ],
        color = "info",
        className = "d-flex align-items-center",
    ),
    dbc_alert(
        [
            html_i(className = "bi bi-check-circle-fill me-2"),
            "An example success alert with an icon",
        ],
        color = "success",
        className = "d-flex align-items-center",
    ),
    dbc_alert(
        [
            html_i(className = "bi bi-exclamation-triangle-fill me-2"),
            "An example warning alert with an icon",
        ],
        color = "warning",
        className = "d-flex align-items-center",
    ),
    dbc_alert(
        [
            html_i(className = "bi bi-x-octagon-fill me-2"),
            "An example danger alert with an icon",
        ],
        color = "danger",
        className = "d-flex align-items-center",
    ),
])

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
module BetweenFlags

include("PerFlagFuncs.jl")
include("UtilityFuncs.jl")
include("FeaturedFuncs.jl")

export Flag
export FlagSet
export get_flat
export get_level
export get_level_flat
export remove_flat

end

ext: jl


text:
push!(LOAD_PATH, "/Users/czhang/Desktop/Projects/dw_/julialib/")
import DimmWitted
DimmWitted.set_libpath("/Users/czhang/Desktop/Projects/dw_/libdw_julia")

######################################
# The following piece of code creates a
# synthetic data set:
#    - Data type is Cdouble
#    - Modle type is Array{Cdouble}
#
nexp = 100000
nfeat = 100
examples = Array(Cdouble, nexp, nfeat+1)
for row = 1:nexp
	for col = 1:nfeat
		examples[row, col] = 1
	end
	if rand() > 0.8
		examples[row, nfeat+1] = 0
	else
		examples[row, nfeat+1] = 1
	end
end
model = Cdouble[0 for i = 1:nfeat]

######################################
# Define the loss function and gradient
# function for logistic regression
#
function loss(row::Array{Cdouble,1}, model::Array{Cdouble,1})
	@inbounds begin
		const label = row[length(row)]
		const nfeat = length(model)
		d = 0.0
		for i = 1:nfeat
			d = d + row[i]*model[i]
		end
	end
	return (-label * d + log(exp(d) + 1.0))
end

function grad(col::Array{Cdouble,1}, _colid::Cint, rows::Array{Array{Cdouble, 1}}, model::Array{Cdouble,1})
	colid = _colid + 1
	nfeat = length(model)
	nrows = length(rows)
	if colid > nfeat	# if it is the last column, then it is the labels, do nothing
		return 1.0
	end

	sum_term = 0.0
	pat_term = 0.0
	for ir = 1:length(rows)
		label = rows[ir][nfeat+1]
		d = 0.0
		for i = 1:nfeat
			d = d + rows[ir][i]*model[i]
		end
		sum_term = sum_term + label*rows[ir][colid]
		pat_term = pat_term + rows[ir][colid]*1.0/(1.0+exp(-d))
	end

	model[colid] = model[colid] - 0.00001* (-sum_term + pat_term)

	return 1.0
end

function avg(p_models::Ptr{Array{Cdouble,1}}, nrepl::Cint, irepl::Cint)
	models = pointer_to_array(p_models, convert(Int64, nrepl+1))
	julia_irepl = irepl + 1
	const nfeat = length(models[julia_irepl])
	for j = 1:nfeat
		d = 0.0
		for i = 1:nrepl
			d = d + models[i][j]
		end
		models[julia_irepl][j] = d/nrepl
	end
	return 1.0
end

######################################
# Create a DimmWitted object using data
# and model. You do not need to specify
# the type, they are infer'ed by the
# open() function, which is parametric.
#
dw = DimmWitted.open(examples, model,
                DimmWitted.MR_PERMACHINE,
                DimmWitted.DR_SHARDING,
                DimmWitted.AC_C2R)

######################################
# Register functions.
#
handle_loss = DimmWitted.register_row(dw, loss)
handle_grad = DimmWitted.register_c2r(dw, grad)
DimmWitted.register_model_avg(dw, handle_loss, avg, true)
DimmWitted.register_model_avg(dw, handle_grad, avg, true)

######################################
# Run 10 epoches.
#
for iepoch = 1:10
	rs = DimmWitted.exec(dw, handle_loss)
	println("LOSS: ", rs/nexp)
	println(sum(model))
	rs = DimmWitted.exec(dw, handle_grad)
end

max_stars_repo_name: GloriaGo/dw

max_issues_repo_name: GloriaGo/dw

max_forks_repo_name: GloriaGo/dw


text:
function Base.intersect(l::LineSegment{Point{2,T}},
                        g::RectilinearGrid{X,Y,0,T}) where {X,Y,T}
    # Intersect the bounding box
    # https://tavianator.com/2011/ray_box.html
    𝘂⁻¹ = 1 ./ (l[2] - l[1])
    X₁ = coordinates(l[1])
    r₁ = 𝘂⁻¹ ⊙ (Vec(xmin(g), ymin(g)) - X₁)
    r₂ = 𝘂⁻¹ ⊙ (Vec(xmax(g), ymax(g)) - X₁)
    rmin = maximum(min.(r₁, r₂))
    rmax = minimum(max.(r₁, r₂))
    # Valid line clipping
    # Assumes 0 ≤ rmin, rmax ≤ 1
    if rmin ≤ rmax
        pstart = l(rmin)
        pend = l(rmax)
        xlower, xupper = minmax(pstart[1], pend[1])
        ylower, yupper = minmax(pstart[2], pend[2])
        # Get the start and stop indices for the range of grid lines
        # that need to be tested
        # We keep the branches instead of using getsortedfirst, since
        # branched will be pruned at compile time.
        if X ≥ sorted_array_findfirst_threshold
            xlower_ind = searchsortedfirst(g.x, xlower)
            xupper_ind = searchsortedfirst(g.x, xupper)
        else
            xlower_ind = findsortedfirst(g.x, xlower)
            xupper_ind = findsortedfirst(g.x, xupper)
        end
        if Y ≥ sorted_array_findfirst_threshold
            ylower_ind = searchsortedfirst(g.y, ylower)
            yupper_ind = searchsortedfirst(g.y, yupper)
        else
            ylower_ind = findsortedfirst(g.y, ylower)
            yupper_ind = findsortedfirst(g.y, yupper)
        end

        xstart, xend = xlower_ind, xupper_ind
        xinc = 1
        if pstart[1] > pend[1]
            xend, xstart = xstart, xend
            xinc = -1
        end
        ixgen = xstart:xinc:xend

        ystart, yend = ylower_ind, yupper_ind
        yinc = 1
        if pstart[2] > pend[2]
            yend, ystart = ystart, yend
            yinc = -1
        end
        iygen = ystart:yinc:yend

        # Allocate for the max number of intersections and delete the end
        # if not needed
        nx = length(ixgen)
        ny = length(iygen)
        nmax = nx + ny
        intersections = Vector{Point{2,T}}(undef, nmax)
        r = 𝘂⁻¹[1]*(g.x[xstart] - X₁[1])
        s = 𝘂⁻¹[2]*(g.y[ystart] - X₁[2])
        i = 2
        j = 2
        k = 1
        # Add the intersection with smallest valid parametric coordinate
        # until we run out of x or y values to test
        while i ≤ nx && j ≤ ny
            if r < s
                if rmin ≤ r ≤ rmax
                    intersections[k] = l(r)
                    k += 1
                end
                r = 𝘂⁻¹[1]*(g.x[ixgen[i]] - X₁[1])
                i += 1
            else
                if rmin ≤ s ≤ rmax
                    intersections[k] = l(s)
                    k += 1
                end
                s = 𝘂⁻¹[2]*(g.y[iygen[j]] - X₁[2])
                j += 1
            end
        end

        # Set the remaining intersections
        while i ≤ nx
            if rmin ≤ r ≤ rmax
                intersections[k] = l(r)
                k += 1
            end
            r = 𝘂⁻¹[1]*(g.x[ixgen[i]] - X₁[1])
            i += 1
        end
        if rmin ≤ r ≤ rmax
            intersections[k] = l(r)
            k += 1
        end

        while j ≤ ny
            if rmin ≤ s ≤ rmax
                intersections[k] = l(s)
                k += 1
            end
            s = 𝘂⁻¹[2]*(g.y[iygen[j]] - X₁[2])
            j += 1
        end
        if rmin ≤ s ≤ rmax
            intersections[k] = l(s)
            k += 1
        end
        # Delete the unused end of the vector
        Base._deleteend!(intersections, nmax - k + 1)
        return intersections
    else
        return Point{2,T}[]
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
using DifferentialEquations
import DifferentialEquations.DiffEqBase: AbstractODESolution, HermiteInterpolation,AbstractDiffEqInterpolation

abstract type AbstractSimulation end
export AbstractSimulation

struct Simulation{Q<:AbstractODESolution,W<:AbstractDomain,L<:AbstractArray,G<:Function,G2<:AbstractArray} <: AbstractSimulation
    sol::Q
    domain::W
    names::L
    N::G
    Ns::G2
end

function Simulation(sol::Q,domain::W) where {Q<:AbstractODESolution,W<:AbstractDomain}
    names = getfield.(domain.phase.species,:name)
    Ns = sum(hcat(sol.interp.u...)[domain.indexes[1]:domain.indexes[2],:],dims=1)
    Nderivs = sum(hcat(sol.interp.du...)[domain.indexes[1]:domain.indexes[2],:],dims=1)
    N = HermiteInterpolation(sol.interp.t,Ns,Nderivs)
    F(t::T) where {T<:Real} = N(t,nothing,Val{0},sol.prob.p,:left)
    return Simulation(sol,domain,names,F,Ns)
end

export Simulation

length(p::T) where {T<:AbstractSimulation} = 1
export length

iterate(p::T) where {T<:AbstractSimulation} = p
export iterate

Broadcast.broadcastable(p::T) where {T<:AbstractSimulation} = Ref(p)
export broadcastable

spcindex(bsol::Z,name::Q) where {Z<:Simulation,Q<:AbstractString} = findfirst(isequal(name),getfield.(bsol.domain.phase.species,:name))
export spcindex

function molefractions(bsol::Q,name::W,t::E) where {Q<:AbstractSimulation, W<:String, E<:Real}
    @assert name in bsol.names
    ind = findfirst(isequal(name),bsol.names)
    return bsol.sol(t)[ind]/bsol.N(t)
end

function molefractions(bsol::Q, t::E) where {Q<:AbstractSimulation,E<:Real}
    return bsol.sol(t)[bsol.domain.indexes[1]:bsol.domain.indexes[2]]./bsol.N(t)
end

function molefractions(bsol::Q) where {Q<:AbstractSimulation}
    @views return hcat(bsol.sol.u...)[bsol.domain.indexes[1]:bsol.domain.indexes[2],:]./bsol.Ns
end

export molefractions

getT(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ConstantTPDomain,ConstantTVDomain},K<:Real,Q,G,L} = bsol.domain.T
getT(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ConstantVDomain,ParametrizedVDomain},K<:Real,Q,G,L} = bsol.sol(t)[bsol.domain.indexes[3]]
getT(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ParametrizedTConstantVDomain,ParametrizedTPDomain},K<:Real,Q,G,L} = bsol.domain.T(t)
export getT
getV(bsol::Simulation{Q,W,L,G}, t::K) where {W<:ConstantTPDomain,K<:Real,Q,G,L} = bsol.N(t)*getT(bsol,t)*R /bsol.domain.P
getV(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ConstantVDomain,ConstantTVDomain,ParametrizedTConstantVDomain},K<:Real,Q,G,L} = bsol.domain.V
getV(bsol::Simulation{Q,W,L,G}, t::K) where {W<:ParametrizedVDomain,K<:Real,Q,G,L} = bsol.domain.V(t)
getV(bsol::Simulation{Q,W,L,G}, t::K) where {W<:ParametrizedTPDomain,K<:Real,Q,G,L} = bsol.N(t)*getT(bsol,t)*R /bsol.domain.P(t)
export getV
getP(bsol::Simulation{Q,W,L,G}, t::K) where {W<:ConstantTPDomain,K<:Real,Q,G,L} = bsol.domain.P
getP(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ConstantTVDomain,ParametrizedTConstantVDomain},K<:Real,Q,G,L} = 1.0e6
getP(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ConstantVDomain,ParametrizedVDomain},K<:Real,Q,G,L} = bsol.N(t)*R*getT(bsol,t)/getV(bsol,t)
getP(bsol::Simulation{Q,W,L,G}, t::K) where {W<:ParametrizedTPDomain,K<:Real,Q,G,L} = bsol.domain.P(t)
export getP
getC(bsol::Simulation{Q,W,L,G}, t::K) where {W<:ConstantTPDomain,K<:Real,Q,G,L} = bsol.domain.P/(R*bsol.domain.T)
getC(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ConstantVDomain,ConstantTVDomain,ParametrizedTConstantVDomain},K<:Real,Q,G,L} = bsol.N(t)/bsol.domain.V
getC(bsol::Simulation{Q,W,L,G}, t::K) where {W<:ParametrizedVDomain,K<:Real,Q,G,L} = bsol.N(t)/bsol.domain.V(t)
getC(bsol::Simulation{Q,W,L,G}, t::K) where {W<:ParametrizedTPDomain,K<:Real,Q,G,L} = bsol.domain.P(t)/(R*bsol.domain.T(t))
export getC
"""
calculates the rates of production/loss at a given time point
this outputs a sparse matrix of  num reactions xnum species containing the production/loss
rate of that species associated with that reaction
"""
function rops(bsol::Q,t::X) where {Q<:Simulation,X<:Real}
    ropmat = zeros(length(bsol.domain.phase.reactions),length(bsol.domain.phase.species))
    cs,kfs,krevs = calcthermo(bsol.domain,bsol.sol(t),t)[[2,9,10]]
    @simd for i in 1:length(bsol.domain.phase.reactions)
        rxn = bsol.domain.phase.reactions[i]
        R = getrate(rxn,cs,kfs,krevs)
        for ind in rxn.productinds
            ropmat[i,ind] += R
        end
        for ind in rxn.reactantinds
            ropmat[i,ind] -= R
        end
    end
    return ropmat
end

"""
calculates the rates of production/loss at a given time point for a single species
this outputs a sparse vector of length num reactions containing the production/loss
rate associated with that reaction for the given species
"""
function rops(bsol::Y,name::X,t::Z) where {Y<:Simulation, X<:AbstractString, Z<:Real}
    rop = spzeros(length(bsol.domain.phase.reactions))
    cs,kfs,krevs = calcthermo(bsol.domain,bsol.sol(t),t)[[2,9,10]]
    ind = findfirst(isequal(name),getfield.(bsol.domain.phase.species,:name))
    @assert !isa(ind,Nothing) "species $name not in species array"
    for (i,rxn) in enumerate(bsol.domain.phase.reactions)
        c = 0
        R = getrate(rxn,cs,kfs,krevs)
        c -= count(isequal(ind),rxn.reactantinds)
        c += count(isequal(ind),rxn.productinds)
        if c != 0
            rop[i] = c*R
        end
    end
    return rop
end
export rops

function getconcentrationsensitivity(bsol::Simulation{Q,W,L,G}, numerator::String, denominator::String, t::K) where {W<:Union{ConstantVDomain,ConstantTVDomain},K<:Real,Q,G,L}
    @assert numerator in bsol.names
    @assert denominator in bsol.names
    indnum = findfirst(isequal(numerator),bsol.names)
    inddeno = findfirst(isequal(denominator),bsol.names)
    Nvars = bsol.domain.indexes[end]-bsol.domain.indexes[1]+1
    Nrxns = length(bsol.domain.phase.reactions)
    arr = bsol.sol(t)
    s = arr[Nvars+Nrxns*Nvars+(inddeno-1)*Nvars+indnum]
    val = s/arr[indnum] #constant volume
    if t == 0
        return 0.0
    else
        return val
    end
end

function getconcentrationsensitivity(bsol::Simulation{Q,W,L,G}, numerator::String, denominator::String, t::K) where {W<:ConstantTPDomain,K<:Real,Q,G,L}
    @assert numerator in bsol.names
    @assert denominator in bsol.names
    indnum = findfirst(isequal(numerator),bsol.names)
    inddeno = findfirst(isequal(denominator),bsol.names)
    Nvars = bsol.domain.indexes[end]-bsol.domain.indexes[1]+1
    Nrxns = length(bsol.domain.phase.reactions)
    arr = bsol.sol(t)
    s = arr[Nvars+Nrxns*Nvars+(inddeno-1)*Nvars+indnum]
    V = getV(bsol,t)
    c = arr[indnum]/V
    val =  (s-c*sum(arr[Nvars+Nrxns*Nvars+(inddeno-1)*Nvars+1:Nvars+Nrxns*Nvars+inddeno*Nvars])*R*bsol.domain.T/bsol.domain.P)/(c*V) #known T and P
    if t == 0
        return 0.0
    else
        return val
    end
end

function getconcentrationsensitivity(bsol::Simulation{Q,W,L,G}, numerator::String, denominator::Z, t::K) where {W<:Union{ConstantVDomain,ConstantTVDomain},K<:Real,Z<:Integer,Q,G,L}
    @assert numerator in bsol.names
    indnum = findfirst(isequal(numerator),bsol.names)
    inddeno = denominator
    Nvars = bsol.domain.indexes[end]-bsol.domain.indexes[1]+1
    Nrxns = length(bsol.domain.phase.reactions)
    arr = bsol.sol(t)
    s = arr[Nvars+(inddeno-1)*Nvars+indnum]
    T = getT(bsol,t)
    P = getP(bsol,t)
    C = getC(bsol,t)
    k = bsol.domain.phase.reactions[inddeno].kinetics(T=T,P=P,C=C)
    val = s*k/arr[indnum] #constant volume
    if t == 0
        return 0.0
    else
        return val
    end
end

function getconcentrationsensitivity(bsol::Simulation{Q,W,L,G}, numerator::String, denominator::Z, t::K) where {W<:ConstantTPDomain,K<:Real,Z<:Integer,Q,G,L}
    @assert numerator in bsol.names
    indnum = findfirst(isequal(numerator),bsol.names)
    inddeno = denominator
    Nvars = bsol.domain.indexes[end]-bsol.domain.indexes[1]+1
    Nrxns = length(bsol.domain.phase.reactions)
    arr = bsol.sol(t)
    s = arr[Nvars+(inddeno-1)*Nvars+indnum]
    V = getV(bsol,t)
    T = getT(bsol,t)
    P = getP(bsol,t)
    C = getC(bsol,t)
    c = arr[indnum]/V
    k = bsol.domain.phase.reactions[inddeno].kinetics(T=T,P=P,C=C)
    val = k*(s-c*sum(arr[Nvars+(inddeno-1)*Nvars+1:Nvars+inddeno*Nvars])*R*bsol.domain.T/bsol.domain.P)/(c*V) #known T and P
    if t == 0
        return 0.0
    else
        return val
    end
end

export getconcentrationsensitivity

"""
calculate the rates of all reactions at time t
"""
function rates(bsol::Q,t::X) where {Q<:Simulation,X<:Real}
    cs,kfs,krevs = calcthermo(bsol.domain,bsol.sol(t),t)[[2,9,10]]
    return [getrate(rxn,cs,kfs,krevs) for rxn in bsol.domain.phase.reactions]
end

"""
calculate the rates of all reactions at given times ts
defaults to using bsol.sol.t if ts is not supplied
"""
function rates(bsol::Q;ts::X=Array{Float64,1}()) where {Q<:Simulation,X<:AbstractArray}
    if length(ts) == 0
        ts = bsol.sol.t
    end
    return hcat([rates(bsol,t) for t in ts]...)
end

export rates


text:
module EL

import ..ED
using ..SimLib
using ..SimLib: Maybe
using LinearAlgebra
using SharedArrays: sdata

export EigenstateLocality, ELDataDescriptor, ELData, load_el

### Descriptor

struct ELDataDescriptor{T} <: ED.EDDerivedDataDescriptor
    operator::Maybe{T}
    operatorname::String
    derivedfrom::ED.EDDataDescriptor
end

ELDataDescriptor(operatorname::String, args...; kwargs...) = ELDataDescriptor(nothing, operatorname, args...; kwargs...)
ELDataDescriptor(operator, operatorname::String, args...; kwargs...) = ELDataDescriptor(operator, operatorname, EDDataDescriptor(args...; kwargs...))

### Data obj

struct ELData{T, N} <: SimLib.AbstractSimpleData
    descriptor::ELDataDescriptor{T}
    data::Array{Float64, N}
end

ED._default_folder(::ELDataDescriptor) = "locality"
ED._filename_addition(eldd::ELDataDescriptor) = "_" * eldd.operatorname

"""
    load_el(operatorname, edd)
    load_el(operatorname, model, diagtype[, location])
"""
load_el(args...; kwargs...) = load(ELDataDescriptor(args...; kwargs...))

## functions

eigenstatelocality(eigen, operator) = eigenstatelocality!(Vector{Float64}(undef, length(eigen.values)-1), eigen, operator)

eigenstatelocality!(res, eigen, operator) = eigenstatelocality!(res, eigen.values, eigen.vectors, operator)

function eigenstatelocality!(res, evals, evecs, operator)
    energies = [eva + real(dot(eve, operator, eve)) for (eva, eve) in zip(evals, eachcol(evecs))]
    order = sortperm(energies)
    for (i, (index, nextindex)) in enumerate(zip(order[1:end-1], order[2:end]))
        res[i] = log(abs(dot(view(evecs, :, index), operator, view(evecs, :, nextindex))) / (energies[nextindex] - energies[index]))
    end
    res
end

### Task

mutable struct ELTask{T} <: ED.EDTask
    operatorname::String
    operator::T
    data
end

EigenstateLocality(operatorname, operator) = ELTask(operatorname, operator, nothing)

function ED.initialize!(task::ELTask, arrayconstructor, spectral_size)
    task.data = arrayconstructor(Float64, spectral_size-1)
end

function ED.compute_task!(task::ELTask, evals, evecs, inds...)
    n = min(size(task.data,1), length(evals)-1)
    eigenstatelocality!(view(task.data, 1:n, inds...), view(evals, 1:n+1), evecs, task.operator)
end

function ED.failed_task!(task::ELTask, inds...)
    task.data[:, inds...] .= NaN64
end

function ED.assemble(task::ELTask, edd)
    ELData(ELDataDescriptor(task.operator, task.operatorname, edd), sdata(task.data))
end

Base.summary(task::ELTask) = string(typeof(task)) * "($(task.operatorname))"

end #module


text:
function import_distributions(t=Distribution, count=0)
    for dist in subtypes(t)
        name = dist.name.name
        obj = Distributions.(name)
        if isgeneric(obj)
            # Concrete.
            eval(parse("import Distributions: $name; export $name"))
            eval(:($name(s::Store, k::Function, address, args...) = k(s, $name(args...))))
            count += 1
        else
            # Abstract.
            count += import_distributions(obj, count)
        end
    end
    count
end

_dist_count = import_distributions()

const Dir = Dirichlet

const _flip_map = (false,true)

@pp function flip(args...)
    _flip_map[~Bernoulli(args...)+1]
end

@pp function randominteger(k)
    ~Categorical(k)
end

# Gradient of the score (logpdf) w.r.t. the parameters.
function gradient(erp::Bernoulli, x)
    g = (x-erp.p) / (erp.p*(1-erp.p))
    @assert isfinite(g)
    [g]
end

immutable Discrete
    x::Vector
    p::Vector{Float64}
    dict::Dict # Potentially un-normalized.
    z::Float64 # Normalizing constant.
    function Discrete(dict)
        x = collect(keys(dict))
        p = collect(values(dict))
        z = sum(p)
        z == 0 && error("Distribution has zero probability mass.")
        z != 1.0 && (p /= z)
        new(x,p,dict,z)
    end
end

rand(erp::Discrete) = erp.x[rand(erp.p)]
support(erp::Discrete) = erp.x
score(erp::Discrete, x) = log(erp.dict[x]/erp.z)

# @pp
Discrete(s::Store, k::Function, address, x, p) = k(s, Discrete(Dict(x,p)))

# This is experimental and maybe removed. It repeatedly calls a thunk
# and wraps the result in an ERP. This seems convinient as existing
# function can be used to plot the results. This has different
# semantics to pmcmc(1,n,thunk) which simulates n executions of the
# program in which thunk is called once, rather than been a single
# execution of a program which calls the thunk n times.

# @pp
function Discrete(s::Store, k::Function, address, comp::Function, n)
    partial(repeat, s)(address, comp, n) do store::Store, samples
        counts = Dict{Any,Int64}()
        for s in samples
            counts[s] = get(counts,s,0) + 1
        end
        k(store, Discrete(counts))
    end
end

# Based on Base.showdict
function show(io::IO, erp::Discrete)
    t = sort(collect(zip(erp.p, erp.x)), rev=true)

    xs = Array(String, length(t))
    keylen = 0
    for (i, (_,x)) in enumerate(t)
        xs[i] = sprint(show, x)
        keylen = max(length(xs[i]), keylen)
    end

    for (i, (p,_)) in enumerate(t)
        print(io, rpad(xs[i], keylen, " "))
        print(io, " | ")
        print(io, p)
        i < length(t) && print(io, "\n")
    end
end

const ERP = Union(Distribution, Discrete)

function hellingerdistance(p::ERP, q::Discrete)
    psupp = support(p)
    qsupp = support(q)
    @assert issubset(qsupp, psupp)
    acc = 0.
    for x in psupp
        px = exp(score(p,x))
        qx = x in qsupp ? exp(score(q,x)) : 0.
        acc += (sqrt(px)-sqrt(qx))^2
    end
    sqrt(acc/2.)
end

function kl(p::ERP, q::ERP)
    psupp = support(p)
    qsupp = support(q)
    @assert length(psupp) == length(qsupp)
    for x in psupp; @assert x in qsupp; end
    sum([exp(score(p,x)) * (score(p,x) - score(q,x)) for x in psupp])
end

kl(s::Store, k::Function, address, p::ERP, q::ERP) = k(s, kl(p,q))


text:
module LarSurf
    using Logging
    using Distributed
	using SparseArrays
	using JLD2, FileIO
	# using Pkg
    # @everywhere using LinearAlgebraicRepresentation
    # Lar = LinearAlgebraicRepresentation
    @info "Loaded M on $(myid())"

    export version
    # export convert
    hello(who::String) = "Hello, $who"

    # @everywhere LarSurf.ahoj("ku")
    # include("ahoj.jl")
    # @everywhere include("ahoj.jl")

	include("lar.jl")
    include("convert.jl")
    include("print_function.jl")
    include("arr_fcn.jl")
	include("triangulation.jl")
    include("read.jl")
    include("plasm.jl") ;
    include("surface.jl") ;
    include("representation.jl") ;
    include("import3d.jl") ;
    include("block.jl") ;
    include("sampledata.jl") ;
    # include("io3d.jl") ;
    include("boundary_operator.jl") ;
    include("inner_boundary.jl") ;
    include("datasets.jl") ;
    include("surface_extraction.jl") ;
    # include("experiment_support.jl") ;
    include("brick_surf_extraction.jl") ;
    include("surface_extraction_parallel_low_comunication.jl") ;
    # include("surface_extraction_parallel.jl") ;
	include("surface_extraction_per_voxel.jl")
	include("smoothing.jl")
	# if "ExSup" in keys(Pkg.installed())
		include("larsurf_experiments.jl")
	# end
	# include("dicom_support.jl")

    function version()
        return "0.0.2"
    end

end


text:
using Pkg; Pkg.activate(@__DIR__)
using MuJoCo
using JLD2
mj_activate("/home/taylor/.mujoco/bin/mjkey.txt") # set location to MuJoCo key path

using LyceumMuJoCo, LyceumMuJoCoViz
using FiniteDiff

using IterativeLQR
using LinearAlgebra
using Random

# ## load MuJoCo model
path = joinpath(@__DIR__, "../../../env/raiberthopper/deps/hopper.xml")
path = joinpath(@__DIR__, "../../../environments/rexhopper/deps/rexhopper.xml")

include("mujoco_model.jl")
hopper = MuJoCoModel(path)
sim = LyceumMuJoCo.MJSim(hopper.m, hopper.d)

# @show hopper.m.body_mass[2] = 1.0
# @show hopper.m.body_mass[3] = 0.1
# @show hopper.m.body_inertia[:, 2] .= 0.004
# @show hopper.m.body_inertia[:, 3] .= 0.0001

@show hopper.m.body_mass[:]
@show hopper.m.body_inertia[:, :]

# @show length(hopper.d.ctrl)

# LyceumMuJoCoViz.visualize(sim) #, trajectories=[states]) # ctrl + LEFT (access trajectory mode)

T = 100
states = Array(undef, statespace(sim), T)
sim.d.qpos[:] = [0.0; 0.0; 0.55; 0.0; 0.0; 0.0; 0.5]
sim.d.qvel[:] = zeros(7)
sim.d.ctrl
for t = 1:T
    sim.d.ctrl[:] = [0.0; 0.0; 1.0 * hopper.m.body_mass[2] * 9.81 / 200; 0.0; 0.0; 0.0; 0.0]
    # sim.d.ctrl .= 1.0 * randn(3)#[0.0; 0.0; 1.0 * hopper.m.body_mass[2] * 9.81]#[5.0; 10.0; 0.0]
    LyceumMuJoCo.step!(sim)
    states[:, t] .= getstate(sim)
end
visualize(sim, trajectories=[states])
# 2.0 * hopper.m.body_mass[2] * 9.81

# ## horizon
T = 101
Tm = convert(Int, floor((T - 1) / 2))

# ## hopper
nx = hopper.nx
nu = hopper.nu

# ## model
dyn = Dynamics(
    (y, x, u, w) -> f!(y, hopper, x, u),
    (dx, x, u, w) -> fx!(dx, hopper, x, u),
    (du, x, u, w) -> fu!(du, hopper, x, u),
    nx, nx, nu)

model = [dyn for t = 1:T-1]

# ## initial conditions
x1 = [0.0; 0.0; 0.55; 0.0; 0.0; 0.0; 0.5; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0]
xM = [0.5; 0.5; 0.55 + 0.5; 0.0; 0.0; 0.0; 0.5; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0]
xT = [0.5; 0.5; 0.55 + 0.0; 0.0; 0.0; 0.0; 0.5; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0]

# # rollout
ū = [[0.0; 0.0; 1.0 * hopper.m.body_mass[2] * 9.81 / 200; 0.0; 0.0; 0.0; 0.0] for t = 1:T-1]
w = [zeros(0) for t = 1:T-1]
x̄ = rollout(model, x1, ū)

# ## objective
obj1 = (x, u, w) -> (transpose(x - xM) * Diagonal([1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 0.001; 0.001; 0.001; 0.001; 0.001; 0.001; 0.001]) * (x - xM) + transpose(u) * Diagonal([1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0]) * u)
obj2 = (x, u, w) -> (transpose(x - xT) * Diagonal([1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 0.001; 0.001; 0.001; 0.001; 0.001; 0.001; 0.001]) * (x - xT) + transpose(u) * Diagonal([1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0]) * u)
objT = (x, u, w) -> (transpose(x - xT) * Diagonal([1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 0.001; 0.001; 0.001; 0.001; 0.001; 0.001; 0.001]) * (x - xT))

ct1 = IterativeLQR.Cost(obj1, nx, nu)
ct2 = IterativeLQR.Cost(obj2, nx, nu)
cT = IterativeLQR.Cost(objT, nx, 0)
obj = [[ct1 for t = 1:Tm]..., [ct2 for t = 1:Tm]..., cT]

# ## constraints
god_control(x, u, w) = u[4:7]
goal(x, u, w) = 1.0 * (x - xT)[collect(1:7)]

jcont = Constraint(god_control, nx, nu)
conT = Constraint(goal, nx, 0)
cons = [[cont for t = 1:T-1]..., conT]

# ## problem
prob = solver(model, obj, cons,
    opts=Options(linesearch=:armijo,
        α_min=1.0e-5,
        obj_tol=1.0e-3,
        grad_tol=1.0e-3,
        con_tol=0.005,
        max_iter=100,
        max_al_iter=5,
        ρ_init=1.0,
        ρ_scale=10.0,
        verbose=true))
initialize_controls!(prob, ū)
initialize_states!(prob, x̄)

# ## solve
@time solve!(prob)

# ## solution
x_sol, u_sol = get_trajectory(prob)
@show IterativeLQR.eval_obj(prob.m_data.obj.costs, prob.m_data.x, prob.m_data.u, prob.m_data.w)
@show prob.s_data.iter[1]
@show norm(goal(prob.m_data.x[T], zeros(0), zeros(0)), Inf)

path = joinpath(@__DIR__, "hopper_solution.jld2")
# @save path x_sol u_sol

# ##
# using Plot
# plot(hcat(x_sol...)')
# plot(hcat(u_sol..., u_sol[end])', linetype=:steppost)

# ## visualize
x_vis = [[x̄[1] for t = 1:10]..., x̄..., [x̄[end] for t = 1:10]...]
x_vis = [[x_sol[1] for t = 1:100]..., x_sol..., [x_sol[end] for t = 1:100]...]
states = Array(undef, statespace(sim), length(x_vis))
for t = 1:length(x_vis)
    sim.d.qpos[:] .= x_vis[t][1:7]
    sim.d.qvel[:] .= x_vis[t][7 .+ (1:7)]
    # sim.d.ctrl .=
    # LyceumMuJoCo.step!(sim)
    states[:, t] .= getstate(sim)
end
visualize(sim, trajectories=[states])

using Plots
plot(hcat(u_sol...)[1:3, :]')
plot(hcat(u_sol...)[4:7, :]')


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    RegularGridDiscretizer(dims)

Discretize spatial region into regular grid of given
dimensions `dims`.
"""
struct RegularGridDiscretizer{N} <: AbstractDiscretizer
  dims::Dims{N}
end

discretize(r::RectangleRegion{T,N}, d::RegularGridDiscretizer{N}) where {N,T} =
  RegularGrid(lowerleft(r), upperright(r), dims=d.dims)


text:
using Gadfly, RDatasets

set_default_plot_size(6inch, 3inch)

plot(dataset("lattice", "singer"), x="VoicePart", y="Height", Geom.beeswarm)


text:
############################################################
## joAbstractLinearOperatorInplace - outer constructors
############################################################

############################################################
## joMatrixInplace - outer constructors

"""
joMatrixInplace outer constructor

    joMatrixInplace(array::AbstractMatrix;
             DDT::DataType=eltype(array),
             RDT::DataType=promote_type(eltype(array),DDT),
             name::String="joMatrixInplace")

Look up argument names in help to joMatrixInplace type.

# Example
- joMatrixInplace(rand(4,3)) # implicit domain and range
- joMatrixInplace(rand(4,3);DDT=Float32) # implicit range
- joMatrixInplace(rand(4,3);DDT=Float32,RDT=Float64)
- joMatrixInplace(rand(4,3);name="my matrix") # adding name

# Notes
- if DDT:<Real for complex matrix then imaginary part will be neglected for transpose/adjoint operator
- if RDT:<Real for complex matrix then imaginary part will be neglected for forward/conjugate operator

"""
function joMatrixInplace(array::AbstractMatrix{EDT};
    DDT::DataType=EDT,RDT::DataType=promote_type(EDT,DDT),name::String="joMatrixInplace") where {EDT}

        (typeof(array)<:DArray || typeof(array)<:SharedArray) && @warn "Creating joMatrixInplace from non-local array like $(typeof(array)) is likely going to have adverse impact on JOLI's health. Please, avoid it."

        farray=factorize(array)
        return joMatrixInplace{DDT,RDT}(name,size(array,1),size(array,2),
            (y1,x1)->mul!(y1,array,x1),
            (y2,x2)->mul!(y2,transpose(array),x2),
            (y3,x3)->mul!(y3,adjoint(array),x3),
            (y4,x4)->mul!(y4,conj(array),x4),
            (y5,x5)->ldiv!(y5,farray,x5),
            (y6,x6)->ldiv!(y6,transpose(farray),x6),
            (y7,x7)->ldiv!(y7,adjoint(farray),x7),
            (y8,x8)->ldiv!(y8,factorize(conj(array)),x8),
            )
end

############################################################
## joLooseMatrixInplace - outer constructors

"""
joLooseMatrixInplace outer constructor

    joLooseMatrixInplace(array::AbstractMatrix;
             DDT::DataType=eltype(array),
             RDT::DataType=promote_type(eltype(array),DDT),
             name::String="joLooseMatrixInplace")

Look up argument names in help to joLooseMatrixInplace type.

# Example
- joLooseMatrixInplace(rand(4,3)) # implicit domain and range
- joLooseMatrixInplace(rand(4,3);DDT=Float32) # implicit range
- joLooseMatrixInplace(rand(4,3);DDT=Float32,RDT=Float64)
- joLooseMatrixInplace(rand(4,3);name="my matrix") # adding name

# Notes
- if DDT:<Real for complex matrix then imaginary part will be neglected for transpose/adjoint operator
- if RDT:<Real for complex matrix then imaginary part will be neglected for forward/conjugate operator

"""
function joLooseMatrixInplace(array::AbstractMatrix{EDT};
    DDT::DataType=EDT,RDT::DataType=promote_type(EDT,DDT),name::String="joLooseMatrixInplace") where {EDT}

        (typeof(array)<:DArray || typeof(array)<:SharedArray) && @warn "Creating joLooseMatrixInplace from non-local array like $(typeof(array)) is likely going to have adverse impact on JOLI's health. Please, avoid it."

        farray=factorize(array)
        return joLooseMatrixInplace{DDT,RDT}(name,size(array,1),size(array,2),
            (y1,x1)->mul!(y1,array,x1),
            (y2,x2)->mul!(y2,transpose(array),x2),
            (y3,x3)->mul!(y3,adjoint(array),x3),
            (y4,x4)->mul!(y4,conj(array),x4),
            (y5,x5)->ldiv!(y5,farray,x5),
            (y6,x6)->ldiv!(y6,transpose(farray),x6),
            (y7,x7)->ldiv!(y7,adjoint(farray),x7),
            (y8,x8)->ldiv!(y8,factorize(conj(array)),x8),
            )
end

############################################################
## joLinearFunctionInplace - outer constructors

export joLinearFunctionInplaceAll, joLinearFunctionInplace_T, joLinearFunctionInplace_A,
       joLinearFunctionInplaceFwd, joLinearFunctionInplaceFwd_T, joLinearFunctionInplaceFwd_A

"""
joLinearFunctionInplace outer constructor

    joLinearFunctionInplaceAll(m::Integer,n::Integer,
        fop::Function,fop_T::Function,fop_A::Function,fop_C::Function,
        iop::Function,iop_T::Function,iop_A::Function,iop_C::Function,
        DDT::DataType,RDT::DataType=DDT;
        name::String="joLinearFunctionInplaceAll")

Look up argument names in help to joLinearFunctionInplace type.

# Notes
- the developer is responsible for ensuring that used functions provide correct DDT & RDT

"""
joLinearFunctionInplaceAll(m::Integer,n::Integer,
    fop::Function,fop_T::Function,fop_A::Function,fop_C::Function,
    iop::Function,iop_T::Function,iop_A::Function,iop_C::Function,
    DDT::DataType,RDT::DataType=DDT;
    name::String="joLinearFunctionInplaceAll") =
        joLinearFunctionInplace{DDT,RDT}(name,m,n,
            fop,fop_T,fop_A,fop_C,
            iop,iop_T,iop_A,iop_C
            )
"""
joLinearFunctionInplace outer constructor

    joLinearFunctionInplace_T(m::Integer,n::Integer,
        fop::Function,fop_T::Function, iop::Function,iop_T::Function,
        DDT::DataType,RDT::DataType=DDT;
        name::String="joLinearFunctionInplace_T")

Look up argument names in help to joLinearFunctionInplace type.

# Notes
- the developer is responsible for ensuring that used functions provide correct DDT & RDT

"""
joLinearFunctionInplace_T(m::Integer,n::Integer,
    fop::Function,fop_T::Function, iop::Function,iop_T::Function,
    DDT::DataType,RDT::DataType=DDT;
    name::String="joLinearFunctionInplace_T") =
        joLinearFunctionInplace{DDT,RDT}(name,m,n,
            fop, fop_T, @joNF, @joNF,
            iop, iop_T, @joNF, @joNF
            )
"""
joLinearFunctionInplace outer constructor

    joLinearFunctionInplace_A(m::Integer,n::Integer,
        fop::Function,fop_A::Function, iop::Function,iop_A::Function,
        DDT::DataType,RDT::DataType=DDT;
        name::String="joLinearFunctionInplace_A")

Look up argument names in help to joLinearFunctionInplace type.

# Notes
- the developer is responsible for ensuring that used functions provide correct DDT & RDT

"""
joLinearFunctionInplace_A(m::Integer,n::Integer,
    fop::Function,fop_A::Function, iop::Function,iop_A::Function,
    DDT::DataType,RDT::DataType=DDT;
    name::String="joLinearFunctionInplace_A") =
        joLinearFunctionInplace{DDT,RDT}(name,m,n,
            fop, @joNF, fop_A, @joNF,
            iop, @joNF, iop_A, @joNF
            )
"""
joLinearFunctionInplace outer constructor

    joLinearFunctionInplaceFwd(m::Integer,n::Integer,
        fop::Function,fop_T::Function,fop_A::Function,fop_C::Function,
        DDT::DataType,RDT::DataType=DDT;
        name::String="joLinearFunctionInplaceAll")

Look up argument names in help to joLinearFunctionInplace type.

# Notes
- the developer is responsible for ensuring that used functions provide correct DDT & RDT

"""
joLinearFunctionInplaceFwd(m::Integer,n::Integer,
    fop::Function,fop_T::Function,fop_A::Function,fop_C::Function,
    DDT::DataType,RDT::DataType=DDT;
    name::String="joLinearFunctionInplaceAll") =
        joLinearFunctionInplace{DDT,RDT}(name,m,n,
            fop,fop_T,fop_A,fop_C,
            @joNF, @joNF, @joNF, @joNF
            )
"""
joLinearFunctionInplace outer constructor

    joLinearFunctionInplaceFwd_T(m::Integer,n::Integer,
        fop::Function,fop_T::Function,
        DDT::DataType,RDT::DataType=DDT;
        name::String="joLinearFunctionInplaceFwd_T")

Look up argument names in help to joLinearFunctionInplace type.

# Notes
- the developer is responsible for ensuring that used functions provide correct DDT & RDT

"""
joLinearFunctionInplaceFwd_T(m::Integer,n::Integer,
    fop::Function,fop_T::Function,
    DDT::DataType,RDT::DataType=DDT;
    name::String="joLinearFunctionInplaceFwd_T") =
        joLinearFunctionInplace{DDT,RDT}(name,m,n,
            fop, fop_T, @joNF, @joNF,
            @joNF, @joNF, @joNF, @joNF
            )
"""
joLinearFunctionInplace outer constructor

    joLinearFunctionInplaceFwd_A(m::Integer,n::Integer,
        fop::Function,fop_A::Function,
        DDT::DataType,RDT::DataType=DDT;
        name::String="joLinearFunctionInplaceFwd_A")

Look up argument names in help to joLinearFunctionInplace type.

# Notes
- the developer is responsible for ensuring that used functions provide correct DDT & RDT

"""
joLinearFunctionInplaceFwd_A(m::Integer,n::Integer,
    fop::Function,fop_A::Function,
    DDT::DataType,RDT::DataType=DDT;
    name::String="joLinearFunctionInplaceFwd_A") =
        joLinearFunctionInplace{DDT,RDT}(name,m,n,
            fop, @joNF, fop_A, @joNF,
            @joNF, @joNF, @joNF, @joNF
            )

############################################################
## joLooseLinearFunctionInplace - outer constructors

export joLooseLinearFunctionInplaceAll, joLooseLinearFunctionInplace_T, joLooseLinearFunctionInplace_A,
       joLooseLinearFunctionInplaceFwd, joLooseLinearFunctionInplaceFwd_T, joLooseLinearFunctionInplaceFwd_A

"""
joLooseLinearFunctionInplace outer constructor

    joLooseLinearFunctionInplaceAll(m::Integer,n::Integer,
        fop::Function,fop_T::Function,fop_A::Function,fop_C::Function,
        iop::Function,iop_T::Function,iop_A::Function,iop_C::Function,
        DDT::DataType,RDT::DataType=DDT;
        name::String="joLooseLinearFunctionInplaceAll")

Look up argument names in help to joLooseLinearFunctionInplace type.

# Notes
- the developer is responsible for ensuring that used functions provide correct DDT & RDT

"""
joLooseLinearFunctionInplaceAll(m::Integer,n::Integer,
    fop::Function,fop_T::Function,fop_A::Function,fop_C::Function,
    iop::Function,iop_T::Function,iop_A::Function,iop_C::Function,
    DDT::DataType,RDT::DataType=DDT;
    name::String="joLooseLinearFunctionInplaceAll") =
        joLooseLinearFunctionInplace{DDT,RDT}(name,m,n,
            fop,fop_T,fop_A,fop_C,
            iop,iop_T,iop_A,iop_C
            )
"""
joLooseLinearFunctionInplace outer constructor

    joLooseLinearFunctionInplace_T(m::Integer,n::Integer,
        fop::Function,fop_T::Function, iop::Function,iop_T::Function,
        DDT::DataType,RDT::DataType=DDT;
        name::String="joLooseLinearFunctionInplace_T")

Look up argument names in help to joLooseLinearFunctionInplace type.

# Notes
- the developer is responsible for ensuring that used functions provide correct DDT & RDT

"""
joLooseLinearFunctionInplace_T(m::Integer,n::Integer,
    fop::Function,fop_T::Function, iop::Function,iop_T::Function,
    DDT::DataType,RDT::DataType=DDT;
    name::String="joLooseLinearFunctionInplace_T") =
        joLooseLinearFunctionInplace{DDT,RDT}(name,m,n,
            fop, fop_T, @joNF, @joNF,
            iop, iop_T, @joNF, @joNF
            )
"""
joLooseLinearFunctionInplace outer constructor

    joLooseLinearFunctionInplace_A(m::Integer,n::Integer,
        fop::Function,fop_A::Function, iop::Function,iop_A::Function,
        DDT::DataType,RDT::DataType=DDT;
        name::String="joLooseLinearFunctionInplace_A")

Look up argument names in help to joLooseLinearFunctionInplace type.

# Notes
- the developer is responsible for ensuring that used functions provide correct DDT & RDT

"""
joLooseLinearFunctionInplace_A(m::Integer,n::Integer,
    fop::Function,fop_A::Function, iop::Function,iop_A::Function,
    DDT::DataType,RDT::DataType=DDT;
    name::String="joLooseLinearFunctionInplace_A") =
        joLooseLinearFunctionInplace{DDT,RDT}(name,m,n,
            fop, @joNF, fop_A, @joNF,
            iop, @joNF, iop_A, @joNF
            )
"""
joLooseLinearFunctionInplace outer constructor

    joLooseLinearFunctionInplaceFwd(m::Integer,n::Integer,
        fop::Function,fop_T::Function,fop_A::Function,fop_C::Function,
        DDT::DataType,RDT::DataType=DDT;
        name::String="joLooseLinearFunctionInplaceAll")

Look up argument names in help to joLooseLinearFunctionInplace type.

# Notes
- the developer is responsible for ensuring that used functions provide correct DDT & RDT

"""
joLooseLinearFunctionInplaceFwd(m::Integer,n::Integer,
    fop::Function,fop_T::Function,fop_A::Function,fop_C::Function,
    DDT::DataType,RDT::DataType=DDT;
    name::String="joLooseLinearFunctionInplaceAll") =
        joLooseLinearFunctionInplace{DDT,RDT}(name,m,n,
            fop,fop_T,fop_A,fop_C,
            @joNF, @joNF, @joNF, @joNF
            )
"""
joLooseLinearFunctionInplace outer constructor

    joLooseLinearFunctionInplaceFwd_T(m::Integer,n::Integer,
        fop::Function,fop_T::Function,
        DDT::DataType,RDT::DataType=DDT;
        name::String="joLooseLinearFunctionInplaceFwd_T")

Look up argument names in help to joLooseLinearFunctionInplace type.

# Notes
- the developer is responsible for ensuring that used functions provide correct DDT & RDT

"""
joLooseLinearFunctionInplaceFwd_T(m::Integer,n::Integer,
    fop::Function,fop_T::Function,
    DDT::DataType,RDT::DataType=DDT;
    name::String="joLooseLinearFunctionInplaceFwd_T") =
        joLooseLinearFunctionInplace{DDT,RDT}(name,m,n,
            fop, fop_T, @joNF, @joNF,
            @joNF, @joNF, @joNF, @joNF
            )
"""
joLooseLinearFunctionInplace outer constructor

    joLooseLinearFunctionInplaceFwd_A(m::Integer,n::Integer,
        fop::Function,fop_A::Function,
        DDT::DataType,RDT::DataType=DDT;
        name::String="joLooseLinearFunctionInplaceFwd_A")

Look up argument names in help to joLooseLinearFunctionInplace type.

# Notes
- the developer is responsible for ensuring that used functions provide correct DDT & RDT

"""
joLooseLinearFunctionInplaceFwd_A(m::Integer,n::Integer,
    fop::Function,fop_A::Function,
    DDT::DataType,RDT::DataType=DDT;
    name::String="joLooseLinearFunctionInplaceFwd_A") =
        joLooseLinearFunctionInplace{DDT,RDT}(name,m,n,
            fop, @joNF, fop_A, @joNF,
            @joNF, @joNF, @joNF, @joNF
            )


text:
# Functions to setup the CapacityExpansion optimization model
"""
    setup_opt_basic(ts_data::ClustData,opt_data::CEPData,config::Dict{String,Any},optimizer::DataType,optimizer_config::Dict{Symbol,Any})
setting up the basic core elements for a CEP-model
- a JuMP Model is setup and the optimizer is configured. The optimizer itself is passed on as a `optimizer`. It's configuration with `optimizer_config` - Each Symbol and the corresponding value in the Dictionary is passed on to the `with_optimizer` function in addition to the optimizer. For Gurobi an example Dictionary could look like `Dict{Symbol,Any}(:Method => 2, :OutputFlag => 0, :Threads => 2)`
- the sets are setup
"""
function setup_opt_basic(ts_data::ClustData,
                            opt_data::OptDataCEP,
                            config::Dict{String,Any},
                            optimizer::DataType,
                            optimizer_config::Dict{Symbol,Any})
   ## MODEL CEP ##
   # Initialize model
   model =  JuMP.Model(with_optimizer(optimizer;optimizer_config...))
   # Initialize info
   info=[config["descriptor"]]
   # Setup set
   set=setup_opt_set(ts_data, opt_data, config)
   # Setup Model CEP
   return OptModelCEP(model,info,set)
 end

"""
    setup_opt_basic_variables!(cep::OptModelCEP,ts_data::ClustData,opt_data::OptDataCEP, scale::Dict{Symbol,Int})
Adding basic variables COST, CAP and GEN based on set
"""
function setup_opt_basic_variables!(cep::OptModelCEP,
                                  ts_data::ClustData,
                                  opt_data::OptDataCEP)
  ## DATA ##
  set=cep.set

  ## VARIABLES ##
  # Cost
  push!(cep.info,"Variable COST[account, impact, tech] in $(set["impact"]["all"].*" "...)") #Note that variable COST is scaled only within the model with the value scale[:COST]: Real-COST [`EUR` or `USD`] = scale[:COST] ⋅ COST (numeric variable within model)
  @variable(cep.model, COST[account=set["account"]["all"],impact=set["impact"]["all"],tech=set["tech"]["all"]])
  # Capacity
  push!(cep.info,"Variable CAP[tech_n, infrastruct, nodes] ≥ 0 in MW") #Note that variable CAP is scaled only within the model with the value scale[:CAP]: Real-CAP ['MW'] = scale[:CAP] ⋅ CAP (numeric variable within model)
  @variable(cep.model, CAP[tech=set["tech"]["node"],infrastruct=set["infrastruct"]["all"] ,node=set["nodes"]["all"]]>=0)
  # Generation #
  push!(cep.info,"Variable GEN[tech_power, carrier, t, k, node] in MW") #Note that variable is scaled only within the model
  @variable(cep.model, GEN[tech=set["tech"]["all"], carrier=set["carrier"][tech], t=set["time_T_period"]["all"], k=set["time_K"]["all"], node=set["nodes"]["all"]])
  #end
  return cep
end

"""
     setup_opt_lost_load!(cep::OptModelCEP, ts_data::ClustData, opt_data::OptDataCEP) set::Dict)
Adding variable SLACK, LL (LostLoad - if demand cannot be met with installed capacity -> Lost Load can be "purchased" to meet demand)
"""
function setup_opt_lost_load!(cep::OptModelCEP,
                                  ts_data::ClustData,
                                  opt_data::OptDataCEP,
                                  scale::Dict{Symbol,Int})
  ## DATA ##
  set=cep.set
  #ts_weights: k - weight of each period:
  ts_weights=ts_data.weights
  #ts_deltas:  t x k - Δt of each segment x period
  ts_deltas=ts_data.delta_t

  ## LOST LOAD ##
  # Slack variable #
  push!(cep.info,"Variable SLACK[carrier, t, k, node] ≥ 0 in MW") #Note that variable is scaled only within the model
  @variable(cep.model, SLACK[carrier=set["carrier"]["lost_load"], t=set["time_T_period"]["all"], k=set["time_K"]["all"], node=set["nodes"]["all"]] >=0)
  # Lost Load variable #
  push!(cep.info,"Variable LL[carrier, node] ≥ 0 in MWh") #Note that variable is scaled only within the model
  @variable(cep.model, LL[carrier=set["carrier"]["lost_load"], node=set["nodes"]["all"]] >=0)
  # Calculation of Lost Load
  ### Scaling: Scaling is applied to all variables based on the parameters provided in the Dictionary scale. Each variable is multiplied with the scaling parameter 'scale[:VARNAME]' for numerical speedup of the code. We typically devide the equations with the scaling parameter of the left side variable: In this example we divided the entire equation with the scaling parameter of :SLACK which is 'scale[:SLACK]'
  push!(cep.info,"LL[carrier, node] = Σ SLACK[carrier, t, k, node] ⋅ ts_weights[k] ⋅ Δt[t,k] ∀ carrier, node")
  @constraint(cep.model, [carrier=set["carrier"]["lost_load"], node=set["nodes"]["all"]], cep.model[:LL][carrier, node]==sum(cep.model[:SLACK][carrier, t, k, node]*ts_weights[k]*ts_deltas[t,k] for t=set["time_T_period"]["all"], k=set["time_K"]["all"])*scale[:SLACK]/scale[:LL])
  return cep
end

"""
     setup_opt_lost_emission!(cep::OptModelCEP, ts_data::ClustData, opt_data::OptDataCEP)
Adding variable LE (LostEmission - if demand cannot be met without breaking Emission-constraint -> Lost Emission can be "purchased" to meet demand with "dirty" production)
"""
function setup_opt_lost_emission!(cep::OptModelCEP,
                              ts_data::ClustData,
                              opt_data::OptDataCEP)
  ## DATA ##
  set=cep.set

  ## LOST EMISSION ##
  # Lost Emission variable #
  push!(cep.info,"Variable LE[impact_{lost_emission}] ≥ 0 in kg")
  @variable(cep.model, LE[impact=set["impact"]["lost_emission"]] >=0)
  return cep
end

"""
     setup_opt_fix_design_variables!(cep::OptModelCEP,ts_data::ClustData, opt_data::OptDataCEP,scale::Dict{Symbol,Int},fixed_design_variables::Dict{String,Any})
Fixing variables CAP based on first stage vars
"""
function setup_opt_fix_design_variables!(cep::OptModelCEP,
                                  ts_data::ClustData,
                                  opt_data::OptDataCEP,
                                  scale::Dict{Symbol,Int},
                                  fixed_design_variables::Dict{String,Any})
  ## DATA ##
  set=cep.set
  cap=fixed_design_variables["CAP"]
  ## VARIABLES ##
  # Line based
  if haskey(set["tech"],"line")
    trans=fixed_design_variables["TRANS"]
    push!(cep.info,"TRANS[tech, 'new', line] = existing infrastructure ∀ tech_trans, line")
    @constraint(cep.model, [line=set["lines"]["all"], tech=set["tech"]["line"]], cep.model[:TRANS][tech,"new",line]==trans[tech, "new", line]/scale[:TRANS])
  end
  # Node based
  push!(cep.info,"CAP[tech, 'new', node] = existing infrastructure ∀ tech_n, node")
  @constraint(cep.model, [node=set["nodes"]["all"], tech=set["tech"]["node"]], cep.model[:CAP][tech,"new",node]==cap[tech, "new", node]/scale[:CAP])
  return cep
end

"""
     setup_opt_cost_var!(cep::OptModelCEP,ts_data::ClustData,opt_data::OptDataCEP,tech_group::String;dir_main_carrier::Symbol=:input,sign_generation::Number=1)
add variable and fixed Costs for the technology defined by `tech_group`
"""
function setup_opt_cost_var!(cep::OptModelCEP,
                            ts_data::ClustData,
                            opt_data::OptDataCEP,
                            scale::Dict{Symbol,Int},
                            tech_group::String;
                            dir_main_carrier::Symbol=:input,
                            sign_generation::Number=1)
    ## DATA ##
    set=cep.set
    #`costs::OptVariable`: costs[tech,node,year,account,impact] - annulized costs [USD in USD/MW_el, CO2 in kg-CO₂-eq./MW_el]`
    costs = opt_data.costs
    #`techs::OptVariable`: techs[tech][tech_group] - OptDataCEPTech
    techs = opt_data.techs
    #ts_weights: k - weight of each period:
    ts_weights=ts_data.weights
    #ts_deltas:  t x k - Δt of each segment x period
    ts_deltas=ts_data.delta_t

    ## COST ##
    # Calculate Variable Costs
    push!(cep.info,"COST['var',impact,tech] = ($sign_generation) ⋅ Σ_{t,k,node}GEN[carrier_$dir_main_carrier,t,k,node]⋅ ts_weights[k] ⋅ ts_deltas[t,k]⋅ var_costs[tech,impact] ∀ impact, tech_$tech_group")
    @constraint(cep.model, [impact=set["impact"]["all"], tech=set["tech"][tech_group]], cep.model[:COST]["var",impact,tech]==sign_generation*sum(cep.model[:GEN][tech,getfield(techs[tech],dir_main_carrier)["carrier"],t,k,node]*ts_weights[k]*ts_deltas[t,k]*costs[tech,node,set["year"]["all"][1],"var",impact] for node=set["nodes"]["all"], t=set["time_T_period"]["all"], k=set["time_K"]["all"])*scale[:GEN]/scale[:COST])
    return cep
  end

"""
     setup_opt_cost_cap!(cep::OptModelCEP,ts_data::ClustData,opt_data::OptDataCEP,tech_group::String)
add variable and fixed Costs for the technology defined by `tech_group`
"""
function setup_opt_cost_cap!(cep::OptModelCEP,
                            ts_data::ClustData,
                            opt_data::OptDataCEP,
                            scale::Dict{Symbol,Int},
                            tech_group::String)
    ## DATA ##
    set=cep.set
    #`costs::OptVariable`: costs[tech,node,year,account,impact] - annulized costs [USD in USD/MW_el, CO2 in kg-CO₂-eq./MW_el]`
    costs = opt_data.costs
    #ts_weights: k - weight of each period:
    ts_weights=ts_data.weights
    #ts_deltas:  t x k - Δt of each segment x period
    ts_deltas=ts_data.delta_t

    # Calculate Fixed Costs
    push!(cep.info,"COST['cap_fix',impact,tech] = Σ_{t,k}(ts_weights ⋅ ts_deltas[t,k])/8760h ⋅ Σ_{node}CAP[tech,'new',node] ⋅ cap_costs[tech,impact] ∀ impact, tech_$tech_group")
    @constraint(cep.model, [impact=set["impact"]["all"], tech=set["tech"][tech_group]], cep.model[:COST]["cap_fix",impact,tech]==sum(ts_weights[k]*ts_deltas[t,k] for t=set["time_T_period"]["all"], k=set["time_K"]["all"])/8760* sum(cep.model[:CAP][tech,"new",node]*costs[tech,node,set["year"]["all"][1],"cap_fix",impact] for node=set["nodes"]["all"])*scale[:CAP]/scale[:COST])
    return cep
end

"""
     setup_opt_cost_trans!(cep::OptModelCEP,ts_data::ClustData,opt_data::OptDataCEP,tech_group::String)
add variable and fixed Costs for the technology defined by `tech_group`
"""
function setup_opt_cost_trans!(cep::OptModelCEP,
                            ts_data::ClustData,
                            opt_data::OptDataCEP,
                            scale::Dict{Symbol,Int},
                            tech_group::String)
    ## DATA ##
    set=cep.set
    #`costs::OptVariable`: costs[tech,node,year,account,impact] - annulized costs [USD in USD/MW_el, CO2 in kg-CO₂-eq./MW_el]`
    costs = opt_data.costs
    #`lines::OptVarible`: lines[tech,line] - OptDataCEPLine
    lines = opt_data.lines
    #ts_weights: k - weight of each period:
    ts_weights=ts_data.weights
    #ts_deltas:  t x k - Δt of each segment x period
    ts_deltas=ts_data.delta_t

    # Calculate Fixed Costs
    push!(cep.info,"COST['cap-fix',impact,tech] = Σ_{t,k}(ts_weights ⋅ ts_deltas[t,k])/8760h ⋅ Σ_{node}(TRANS[tech,'new',line] ⋅ length[line]) ⋅ (cap_costs[tech,impact]+fix_costs[tech,impact]) ∀ impact, tech_$tech_group")
    @constraint(cep.model, [impact=set["impact"]["all"], tech=set["tech"][tech_group]], cep.model[:COST]["cap_fix",impact,tech] == sum(ts_weights[k]*ts_deltas[t,k] for t=set["time_T_period"]["all"], k=set["time_K"]["all"])/8760* sum(cep.model[:TRANS][tech,"new",line]*lines[tech,line].length *(costs[tech,lines[tech,line].node_start,set["year"]["all"][1],"cap_fix",impact]) for line=set["lines"]["all"])*scale[:TRANS]/scale[:COST])
    return cep
end

"""
     setup_opt_demand!(cep::OptModelCEP,ts_data::ClustData,opt_data::OptDataCEP,scale::Dict{Symbol,Int})
add variable and fixed Costs and limit generation to installed capacity (and limiting time_series, if dependency in techs defined) for fossil and renewable power plants
"""
function setup_opt_demand!(cep::OptModelCEP,
                            ts_data::ClustData,
                            opt_data::OptDataCEP,
                            scale::Dict{Symbol,Int})
    ## DATA ##
    set=cep.set
    #`costs::OptVariable`: costs[tech,node,year,account,impact] - annulized costs [USD in USD/MW_el, CO2 in kg-CO₂-eq./MW_el]`
    costs = opt_data.costs
    #`techs::OptVariable`: techs[tech][tech_group] - OptDataCEPTech
    techs = opt_data.techs
    #`nodes::OptVariable`: nodes[tech,node] - OptDataCEPNode
    nodes = opt_data.nodes
    #ts          Dict( tech-node ): t x k
    ts=ts_data.data
    #ts_weights: k - weight of each period:
    ts_weights=ts_data.weights
    #ts_deltas:  t x k - Δt of each segment x period
    ts_deltas=ts_data.delta_t

    ## DEMAND COST ##
    #Variable Costs
    setup_opt_cost_var!(cep,ts_data,opt_data,scale,"demand")
    #Fixed Costs based on installed capacity
    setup_opt_cost_cap!(cep,ts_data,opt_data,scale,"demand")

    ## DEMAND GENERATION ##
    # Fix the demand
    push!(cep.info,"GEN[tech, carrier t, k, node] = (-1) ⋅ Σ_{infrastruct} CAP[tech,infrastruct,node] * ts[tech-node,t,k] ∀ node, tech_demand, t, k")
    @constraint(cep.model, [tech=set["tech"]["demand"], node=set["nodes"]["all"], t=set["time_T_period"]["all"], k=set["time_K"]["all"]], cep.model[:GEN][tech,techs[tech].input["carrier"],t,k,node] == (-1) * sum(cep.model[:CAP][tech,infrastruct,node] for infrastruct=set["infrastruct"]["all"])*ts[techs[tech].output["timeseries"]*"-"*node][t,k]*scale[:CAP]/scale[:GEN])
    return cep
end

"""
     setup_opt_dispatchable_generation!(cep::OptModelCEP,ts_data::ClustData,opt_data::OptDataCEP,scale::Dict{Symbol,Int})
add variable and fixed Costs and limit generation to installed capacity (and limiting time_series, if dependency in techs defined) for fossil and renewable power plants
"""
function setup_opt_dispatchable_generation!(cep::OptModelCEP,
                            ts_data::ClustData,
                            opt_data::OptDataCEP,
                            scale::Dict{Symbol,Int})
    ## DATA ##
    set=cep.set
    #`costs::OptVariable`: costs[tech,node,year,account,impact] - annulized costs [USD in USD/MW_el, CO2 in kg-CO₂-eq./MW_el]`
    costs = opt_data.costs
    #`techs::OptVariable`: techs[tech][tech_group] - OptDataCEPTech
    techs = opt_data.techs
    #`nodes::OptVariable`: nodes[tech,node] - OptDataCEPNode
    nodes = opt_data.nodes
    #ts          Dict( tech-node ): t x k
    ts=ts_data.data
    #ts_weights: k - weight of each period:
    ts_weights=ts_data.weights
    #ts_deltas:  t x k - Δt of each segment x period
    ts_deltas=ts_data.delta_t

    ## DISPACHABLE COST ##
    #Variable Costs
    setup_opt_cost_var!(cep,ts_data,opt_data,scale,"dispatchable_generation";dir_main_carrier=:output)
    #Fixed Costs based on installed capacity
    setup_opt_cost_cap!(cep,ts_data,opt_data,scale,"dispatchable_generation")

    ## DISPATCHABLE GENERATION ##
    # Limit the generation of dispathables to the infrastructing capacity of dispachable power plants
    push!(cep.info,"0 ≤ GEN[tech, carrier t, k, node] ≤ Σ_{infrastruct} CAP[tech,infrastruct,node] ∀ node, tech_dispatchable_generation, t, k")
    # Limit the generation of dispathables to the infrastructing capacity of dispachable power plants
    @constraint(cep.model, [tech=set["tech"]["dispatchable_generation"], node=set["nodes"]["all"], t=set["time_T_period"]["all"], k=set["time_K"]["all"]], 0 <=cep.model[:GEN][tech, techs[tech].output["carrier"], t, k, node])
    @constraint(cep.model, [tech=set["tech"]["dispatchable_generation"], node=set["nodes"]["all"], t=set["time_T_period"]["all"], k=set["time_K"]["all"]],     cep.model[:GEN][tech, techs[tech].output["carrier"], t, k, node] <=sum(cep.model[:CAP][tech,infrastruct,node] for infrastruct=set["infrastruct"]["all"])*scale[:CAP]/scale[:GEN])
    return cep
end

"""
     setup_opt_non_dispatchable_generation!(cep::OptModelCEP,ts_data::ClustData,opt_data::OptDataCEP,scale::Dict{Symbol,Int})
add variable and fixed Costs and limit generation to installed capacity (and limiting time_series, if dependency in techs defined) for fossil and renewable power plants
"""
function setup_opt_non_dispatchable_generation!(cep::OptModelCEP,
                            ts_data::ClustData,
                            opt_data::OptDataCEP,
                            scale::Dict{Symbol,Int})
    ## DATA ##
    set=cep.set
    #`costs::OptVariable`: costs[tech,node,year,account,impact] - annulized costs [USD in USD/MW_el, CO2 in kg-CO₂-eq./MW_el]`
    costs = opt_data.costs
    #`techs::OptVariable`: techs[tech][tech_group] - OptDataCEPTech
    techs = opt_data.techs
    #`nodes::OptVariable`: nodes[tech,node] - OptDataCEPNode
    nodes = opt_data.nodes
    #ts          Dict( tech-node ): t x k
    ts=ts_data.data

    ## GENERATION COST ##
    #Variable Costs
    setup_opt_cost_var!(cep,ts_data,opt_data,scale,"non_dispatchable_generation";dir_main_carrier=:output)
    #Fixed Costs based on installed capacity
    setup_opt_cost_cap!(cep,ts_data,opt_data,scale,"non_dispatchable_generation")

    ## GENERATION ELECTRICITY ##
    # Limit the generation of dispathables to the infrastructing capacity of non-dispachable power plants
    push!(cep.info,"0 ≤ GEN[tech, carrier t, k, node] ≤ Σ_{infrastruct}CAP[tech,infrastruct,node]*ts[tech-node,t,k] ∀ node, tech_generation{non_dispatchable}, t, k")
    # Limit the generation of non-dispathable generation to the infrastructing capacity of non-dispachable power plants
    @constraint(cep.model, [tech=set["tech"]["non_dispatchable_generation"], node=set["nodes"]["all"], t=set["time_T_period"]["all"], k=set["time_K"]["all"]], 0 <=cep.model[:GEN][tech, techs[tech].output["carrier"], t, k, node])
    @constraint(cep.model, [tech=set["tech"]["non_dispatchable_generation"], node=set["nodes"]["all"], t=set["time_T_period"]["all"], k=set["time_K"]["all"]],  cep.model[:GEN][tech, techs[tech].output["carrier"], t,k,node] <= sum(cep.model[:CAP][tech,infrastruct,node] for infrastruct=set["infrastruct"]["all"])*ts[techs[tech].input["timeseries"]*"-"*node][t,k]*scale[:CAP]/scale[:GEN])
    return cep
end

"""
     setup_opt_conversion!(cep::OptModelCEP,ts_data::ClustData,opt_data::OptDataCEP,scale::Dict{Symbol,Int})

A conversion technology converts the input carrier to an output carrier with a certain efficiency. The costs and capacities are scaled with the input carriers unit
"""
function setup_opt_conversion!(cep::OptModelCEP,
                            ts_data::ClustData,
                            opt_data::OptDataCEP,
                            scale::Dict{Symbol,Int})
    ## DATA ##
    set=cep.set
    #`costs::OptVariable`: costs[tech,node,year,account,impact] - annulized costs [USD in USD/MW_el, CO2 in kg-CO₂-eq./MW_el]`
    costs = opt_data.costs
    #`techs::OptVariable`: techs[tech][tech_group] - OptDataCEPTech
    techs = opt_data.techs
    #ts_weights: k - weight of each period:
    ts_weights=ts_data.weights
    #ts_deltas:  t x k - Δt of each segment x period
    ts_deltas=ts_data.delta_t

    ## CONVERSION COST ##
    #Variable Costs
    setup_opt_cost_var!(cep,ts_data,opt_data,scale,"conversion";sign_generation=(-1))
    #Fixed Costs based on installed capacity
    setup_opt_cost_cap!(cep,ts_data,opt_data,scale,"conversion")

    ##CONVERSION GEN ##
    #Calculate the input generation
    push!(cep.info,"0 ≥ GEN[carrier_input, tech, t, k, node] ≥ (-1) ⋅ Σ_{infrastruct} CAP[tech,infrastruct,node] ∀ node, tech_conversion, t, k")
    @constraint(cep.model, [tech=set["tech"]["conversion"], carrier=[techs[tech].input["carrier"]], node=set["nodes"]["all"], t=set["time_T_period"]["all"], k=set["time_K"]["all"]], 0 >= cep.model[:GEN][tech,carrier,t,k,node])
    @constraint(cep.model, [tech=set["tech"]["conversion"], carrier=[techs[tech].input["carrier"]], node=set["nodes"]["all"], t=set["time_T_period"]["all"], k=set["time_K"]["all"]], cep.model[:GEN][tech,carrier,t,k,node]>= (-1)*sum(cep.model[:CAP][tech,infrastruct,node] for infrastruct=set["infrastruct"]["all"])*scale[:CAP]/scale[:GEN])
    #Calculate the output generation
    push!(cep.info,"GEN[carrier_output,tech, t, k, node] = (-1) ⋅ η[carrier_output, tech] ⋅ GEN[carrier_input, tech, t, k, node] ∀ node, tech_storage_in, t, k")
    @constraint(cep.model, [tech=set["tech"]["conversion"], carrier_in=[techs[tech].input["carrier"]], carrier_out=[techs[tech].output["carrier"]], node=set["nodes"]["all"], t=set["time_T_period"]["all"], k=set["time_K"]["all"]], cep.model[:GEN][tech,carrier_out,t,k,node] == (-1) * techs[tech].constraints["efficiency"] * cep.model[:GEN][tech,carrier_in,t,k,node])
    #TODO add multiple outputs and inputs
    return cep
end

"""
    setup_opt_intertech_cap!(cep::OptModelCEP,
                                      ts_data::ClustData,
                                      opt_data::OptDataCEP,
                                      scale::Dict{Symbol,Int})
Setup constraints for capacities of technologies between different technologies
"""
function setup_opt_intertech_cap!(cep::OptModelCEP,
                                      ts_data::ClustData,
                                      opt_data::OptDataCEP,
                                      scale::Dict{Symbol,Int})
    ## DATA ##
    set=cep.set
    #`techs::OptVariable`: techs[tech][tech_group] - OptDataCEPTech
    techs = opt_data.techs

    #Check for specific constraint that bounds the installed capacity of one technology to another
    for tech in set["tech"]["all"]
      if haskey(techs[tech].constraints,"cap_eq")
          # With 'cap_eq_multiply', one can define the relationship between two technologies. This is useful for setting fixed storage hours: For a system
          # with 4 hours of storage, set this value for 'bat_e' to '4.0' (assuming hourly delta t).
          if haskey(techs[tech].constraints,"cap_eq_multiply")
            cap_eq_multiply =   techs[tech].constraints["cap_eq_multiply"]  #parse(Float64,techs[tech].constraints["cap_eq_multiply"])
          else
            cap_eq_multiply = 1.0
          end
          push!(cep.info,"CAP[$tech, 'new', node] = CAP[$(techs[tech].constraints["cap_eq"]), 'new', node] *$cap_eq_multiply ∀ node, tech_{EUR-Cap-Cost out/in==0}")
          @constraint(cep.model, [node=set["nodes"]["all"]], cep.model[:CAP][tech,"new",node] == cap_eq_multiply * cep.model[:CAP][techs[tech].constraints["cap_eq"],"new",node])
      end
    end
    return cep
end

"""
     setup_opt_storage!(cep::OptModelCEP,ts_data::ClustData,opt_data::OptDataCEP, scale::Dict{Symbol,Int})
add variables INTRASTORGEN and INTRASTOR, variable and fixed Costs, limit generation to installed power-capacity, connect simple-storage levels (within period) with generation
basis for either simplestorage or seasonalstorage
"""
function setup_opt_storage!(cep::OptModelCEP,
                            ts_data::ClustData,
                            opt_data::OptDataCEP,
                            scale::Dict{Symbol,Int})
    ## DATA ##
    set=cep.set
    #`costs::OptVariable`: costs[tech,node,year,account,impact] - annulized costs [USD in USD/MW_el, CO2 in kg-CO₂-eq./MW_el]`
    costs = opt_data.costs
    #`techs::OptVariable`: techs[tech][tech_group] - OptDataCEPTech
    techs = opt_data.techs
    #ts_weights: k - weight of each period:
    ts_weights=ts_data.weights
    #ts_deltas:  t x k - Δt of each segment x period
    ts_deltas=ts_data.delta_t

    ## VARIABLE ##
    # Storage has additional element 0 for storage at hour 0 of day
    push!(cep.info,"Variable INTRASTOR[tech_storage, carrier, t, k, node] ≥ 0 in MWh") #Note that variable is scaled only within the model
    @variable(cep.model, INTRASTOR[tech=set["tech"]["storage"], carrier=set["carrier"][tech],  t=set["time_T_point"]["all"], k=set["time_K"]["all"], node=set["nodes"]["all"]] >=0)

    ## STORAGE COST ##
    #Variable Costs
    # Calculate Variable Costs
    push!(cep.info,"COST['var',impact,tech] = 0 ∀ impact, tech_storage")
    @constraint(cep.model, [impact=set["impact"]["all"], tech=set["tech"]["storage"]], cep.model[:COST]["var",impact,tech]==0)
    #Fixed Costs based on installed capacity
    setup_opt_cost_cap!(cep,ts_data,opt_data,scale,"storage")

    ## STORAGE LEVEL ##
    # Connect the previous storage level and the integral of the flows with the new storage level
    push!(cep.info,"INTRASTOR[carrier,tech, t, k, node] = INTRASTOR[carrier,tech, t-1, k, node] η[tech]^(ts_deltas[t,k]/732h) + ts_deltas[t,k] ⋅ (-1) ⋅ (GEN[tech, carrier t, k, node] ∀ carrier(tech), node, tech_storage, t, k")
    @constraint(cep.model, [node=set["nodes"]["all"], tech=set["tech"]["storage"], t in set["time_T_period"]["all"], k=set["time_K"]["all"]], cep.model[:INTRASTOR][tech, techs[tech].input["carrier"],t,k,node]==cep.model[:INTRASTOR][tech, techs[tech].input["carrier"],t-1,k,node]*(techs[tech].constraints["efficiency"])^(ts_deltas[t,k]/732) - ts_deltas[t,k] * (cep.model[:GEN][tech, techs[tech].input["carrier"], t,k,node])*scale[:GEN]/scale[:INTRASTOR])

    return cep
end

"""
     setup_opt_simplestorage!(cep::OptModelCEP,ts_data::ClustData,opt_data::OptDataCEP, scale::Dict{Symbol,Int})
Adding only intra-day storage:
Looping constraint for each period (same start and end level for all periods) and limit storage to installed energy-capacity
"""
function setup_opt_simplestorage!(cep::OptModelCEP,
                            ts_data::ClustData,
                            opt_data::OptDataCEP,
                            scale::Dict{Symbol,Int})
    ## DATA ##
    set=cep.set
    #`techs::OptVariable`: techs[tech][tech_group] - OptDataCEPTech
    techs = opt_data.techs

    ## INTRASTORAGE ##
    # Limit the storage of the energy part of the battery to its installed power
    push!(cep.info,"INTRASTOR[carrier,tech, t, k, node] ≤ Σ_{infrastruct} CAP[tech,infrastruct,node] ∀ node, tech_storage, t, k")
    @constraint(cep.model, [node=set["nodes"]["all"], tech=set["tech"]["storage"], t=set["time_T_period"]["all"], k=set["time_K"]["all"]], cep.model[:INTRASTOR][tech, techs[tech].input["carrier"], t,k,node]<=sum(cep.model[:CAP][tech, infrastruct, node] for infrastruct=set["infrastruct"]["all"])*scale[:CAP]/scale[:INTRASTOR])
    # Set storage level at beginning and end of day equal
    push!(cep.info,"INTRASTOR[carrier,tech, '0', k, node] = INTRASTOR[carrier,tech, 't[end]', k, node] ∀ node, tech_storage, k")
    @constraint(cep.model, [node=set["nodes"]["all"], tech=set["tech"]["storage"], k=set["time_K"]["all"]], cep.model[:INTRASTOR][tech, techs[tech].input["carrier"], 0, k, node]== cep.model[:INTRASTOR][tech,techs[tech].input["carrier"],set["time_T_point"]["all"][end],k,node])
    # Set the storage level at the beginning of each representative day to the same
    push!(cep.info,"INTRASTOR[carrier,tech, '0', k, node] = INTRASTOR[carrier,tech, '0', k, node] ∀ node, tech_storage, k")
    @constraint(cep.model, [node=set["nodes"]["all"], tech=set["tech"]["storage"], k=set["time_K"]["all"]], cep.model[:INTRASTOR][tech, techs[tech].input["carrier"], 0, k, node]== cep.model[:INTRASTOR][tech, techs[tech].input["carrier"], 0, 1, node])
    return cep
end

"""
     setup_opt_seasonalstorage!(cep::OptModelCEP,ts_data::ClustData,opt_data::OptDataCEP, scale::Dict{Symbol,Int})
Adding inter-day storage:
add variable INTERSTOR, calculate seasonal-storage-level and limit total storage to installed energy-capacity
"""
function setup_opt_seasonalstorage!(cep::OptModelCEP,
                            ts_data::ClustData,
                            opt_data::OptDataCEP,
                            scale::Dict{Symbol,Int})
    ## DATA ##
    set=cep.set
    #`techs::OptVariable`: techs[tech][tech_group] - OptDataCEPTech
    techs = opt_data.techs
    #K identification numbers
    k_ids=ts_data.k_ids

    ## VARIABLE ##
    # Storage
    push!(cep.info,"Variable INTERSTOR[carrier, tech, i, node] ≥ 0 in MWh") #Note that variable is scaled only within the model
    @variable(cep.model, INTERSTOR[tech=set["tech"]["storage"], carrier=set["carrier"][tech], i=set["time_I_point"]["all"], node=set["nodes"]["all"]]>=0)

    ## INTERSTORAGE ##
    # Set storage level at the beginning of the year equal to the end of the year
    push!(cep.info,"INTERSTOR[carrier,tech, '0', node] = INTERSTOR[carrier,tech, 'end', node] ∀ node, tech_storage, t, k")
    @constraint(cep.model, [node=set["nodes"]["all"], tech=set["tech"]["storage"]], cep.model[:INTERSTOR][tech,techs[tech].input["carrier"],0,node]== cep.model[:INTERSTOR][tech, techs[tech].input["carrier"], set["time_I_point"]["all"][end],node])
    # Connect the previous seasonal-storage level and the daily difference of the corresponding simple-storage with the new seasonal-storage level
    push!(cep.info,"INTERSTOR[carrier,tech, i+1, node] = INTERSTOR[carrier,tech, i, node] + INTRASTOR[carrier,tech, 'k[i]', 't[end]', node] - INTRASTOR[carrier,tech, 'k[i]', '0', node] ∀ node, tech_storage, i")
    # Limit the total storage (seasonal and simple) to be greater than zero and less than total storage cap
    push!(cep.info,"0 ≤ INTERSTOR[carrier,tech, i, node] + INTRASTOR[carrier,tech, t, k[i], node] ≤ Σ_{infrastruct} CAP[tech,infrastruct,node] ∀ node, tech_storage, i, t")
    push!(cep.info,"0 ≤ INTERSTOR[carrier,tech, i, node] + INTRASTOR[carrier,tech, t, k[i], node] ≤ Σ_{infrastruct} CAP[tech,infrastruct,node] ∀ node, tech_storage, i, t")
    for i in set["time_I_period"]["all"]
        @constraint(cep.model, [node=set["nodes"]["all"], tech=set["tech"]["storage"]], cep.model[:INTERSTOR][tech,techs[tech].input["carrier"],i,node] == cep.model[:INTERSTOR][tech,techs[tech].input["carrier"],i-1,node] + (cep.model[:INTRASTOR][tech,techs[tech].input["carrier"],set["time_T_period"]["all"][end],k_ids[i],node] - cep.model[:INTRASTOR][tech,techs[tech].input["carrier"],0,k_ids[i],node])*scale[:INTRASTOR]/scale[:INTERSTOR])
        @constraint(cep.model, [node=set["nodes"]["all"], tech=set["tech"]["storage"], t=set["time_T_point"]["all"]], 0 <= cep.model[:INTERSTOR][tech,techs[tech].input["carrier"],i,node]+cep.model[:INTRASTOR][tech,techs[tech].input["carrier"],t,k_ids[i],node]*scale[:INTRASTOR]/scale[:INTERSTOR])
        @constraint(cep.model, [node=set["nodes"]["all"], tech=set["tech"]["storage"], t=set["time_T_point"]["all"]], cep.model[:INTERSTOR][tech,techs[tech].input["carrier"],i,node]+cep.model[:INTRASTOR][tech,techs[tech].input["carrier"],t,k_ids[i],node]*scale[:INTRASTOR]/scale[:INTERSTOR] <= sum(cep.model[:CAP][tech,infrastruct,node] for infrastruct=set["infrastruct"]["all"])*scale[:CAP]/scale[:INTERSTOR])
    end
    return cep
end

"""
     setup_opt_transmission!(cep::OptModelCEP,ts_data::ClustData,opt_data::OptDataCEP, scale::Dict{Symbol,Int})
Setup variable FLOW and TRANS, calculate fixed and variable COSTs, set CAP-trans to zero, limit FLOW with TRANS, calculate GEN-trans for each node
"""
function setup_opt_transmission!(cep::OptModelCEP,
                            ts_data::ClustData,
                            opt_data::OptDataCEP,
                            scale::Dict{Symbol,Int})
    ## DATA ##
    set=cep.set
    #`costs::OptVariable`: costs[tech,node,year,account,impact] - annulized costs [USD in USD/MW_el, CO2 in kg-CO₂-eq./MW_el]`
    costs = opt_data.costs
    #`lines::OptVarible`: lines[tech,line] - OptDataCEPLine
    lines = opt_data.lines
    #ts_weights: k - weight of each period:
    ts_weights=ts_data.weights
    #ts_deltas:  t x k - Δt of each segment x period
    ts_deltas=ts_data.delta_t
    ## VARIABLE ##
    # Add varibale FLOW
    push!(cep.info,"Variable FLOW[tech_transmission, carrier, dir, t, k, line] ≥ 0 in MW") #Note that variable is scaled only within the model
    @variable(cep.model, FLOW[tech=set["tech"]["transmission"], carrier=set["carrier"][tech], dir=set["dir_transmission"]["all"], t=set["time_T_period"]["all"], k=set["time_K"]["all"], line=set["lines"]["all"]] >= 0)
    # Add variable TRANS
    push!(cep.info,"Variable TRANS[tech_transmission,  infrastruct, lines] ≥ 0 in MW") #Note that variable is scaled only within the model
    @variable(cep.model, TRANS[tech=set["tech"]["transmission"], infrastruct=set["infrastruct"]["all"], line=set["lines"]["all"]] >= 0)

    ## TRANSMISSION COST ##
    #Variable Costs
    push!(cep.info,"COST['var',impact,tech] = 0 ∀ impact, tech_storage")
    @constraint(cep.model, [impact=set["impact"]["all"], tech=set["tech"]["transmission"]], cep.model[:COST]["var",impact,tech]==0)
    #Fixed Costs based on installed transmission
    setup_opt_cost_trans!(cep,ts_data,opt_data,scale,"transmission")

    ## TRANSMISSION TRANS ##
    # Limit the flow per line to the existing infrastructure
    push!(cep.info,"| FLOW[carrier, dir, tech, t, k, line] | ≤ Σ_{infrastruct}TRANS[tech,infrastruct,line] ∀ line, tech_transmission, t, k")
    @constraint(cep.model, [tech=set["tech"]["transmission"], carrier=set["carrier"][tech], line=set["lines"]["all"], dir=set["dir_transmission"]["all"],  t=set["time_T_period"]["all"], k=set["time_K"]["all"]], cep.model[:FLOW][tech, carrier,dir, t, k, line] <= sum(cep.model[:TRANS][tech,infrastruct,line] for infrastruct=set["infrastruct"]["all"])*scale[:TRANS]/scale[:FLOW])
    # Calculate the sum of the flows for each node
    push!(cep.info,"GEN[tech, carrier t, k, node] = Σ_{line-end(node)} FLOW[tech, carrier, 'uniform', t, k, line] - Σ_{line_pos} FLOW[tech, carrier, 'opposite', t, k, line] / (η[tech]⋅length[line]) + Σ_{line-start(node)} Σ_{line_pos} FLOW[tech, carrier, 'opposite', t, k, line] - FLOW[tech, carrier, 'uniform', t, k, line] / (η[tech]⋅length[line])∀ tech_transmission, t, k")
    for node in set["nodes"]["all"]
      @constraint(cep.model, [tech=set["tech"]["transmission"], carrier=set["carrier"][tech], t=set["time_T_period"]["all"], k=set["time_K"]["all"]], cep.model[:GEN][tech, carrier, t, k, node] == (sum(cep.model[:FLOW][tech, carrier,"uniform", t, k, line_end] - cep.model[:FLOW][tech, carrier,"opposite", t, k, line_end]/lines[tech,line_end].eff for line_end=set["lines"]["all"][getfield.(lines[tech,:], :node_end).==node]) + sum(cep.model[:FLOW][tech, carrier,"opposite", t, k, line_start] - cep.model[:FLOW][tech, carrier,"uniform", t, k, line_start]/lines[tech,line_start].eff for line_start=set["lines"]["all"][getfield.(lines[tech,:], :node_start).==node]))*scale[:FLOW]/scale[:GEN])
    end
    return cep
end

"""
    setup_opt_energy_balance_transmission!(cep::OptModelCEP,ts_data::ClustData,opt_data::OptDataCEP,scale::Dict{Symbol,Int})
Add energy-balance which shall be matched by the generation (GEN) taking the transmission into account
"""
function setup_opt_energy_balance_transmission!(cep::OptModelCEP,
                            ts_data::ClustData,
                            opt_data::OptDataCEP,
                            scale::Dict{Symbol,Int})
  ## DATA ##
  set=cep.set

  ## ENERGY BALANCE ##
  # Force the demand and slack to match the generation either with transmission
  push!(cep.info,"Σ_{tech}GEN[tech, carriert,k,node] + SLACK[carrier,t,k,node] = 0 ∀ node,t,k")
  @constraint(cep.model, [carrier=intersect(set["carrier"]["all"],set["carrier"]["lost_load"]),node=set["nodes"]["all"], t=set["time_T_period"]["all"], k=set["time_K"]["all"]], sum(cep.model[:GEN][tech,carrier,t,k,node] for tech=set["tech"][carrier]) + (cep.model[:SLACK][carrier,t,k,node]*scale[:SLACK])/scale[:GEN] ==0)
  # Force the demand without slack to match the generation either with transmission
  push!(cep.info,"Σ_{tech}GEN[tech, carriert,k,node] = 0 ∀ node,t,k")
  @constraint(cep.model, [carrier = setdiff(set["carrier"]["all"], set["carrier"]["lost_load"]), node=set["nodes"]["all"], t=set["time_T_period"]["all"], k=set["time_K"]["all"]], sum(cep.model[:GEN][tech,carrier,t,k,node] for tech=set["tech"][carrier]) == 0)
  return cep
end

"""
    setup_opt_energy_balance_copperplate!(cep::OptModelCEP,ts_data::ClustData,opt_data::OptDataCEP,scale::Dict{Symbol,Int})
Add energy-balance using a copperplate assumption that the generation just has to be matched across all nodes without taking transmission into account
"""
function setup_opt_energy_balance_copperplate!(cep::OptModelCEP,
                            ts_data::ClustData,
                            opt_data::OptDataCEP,
                            scale::Dict{Symbol,Int})
  ## DATA ##
  set=cep.set

  ## ENERGY BALANCE ##
  push!(cep.info,"Σ_{tech,node}GEN[tech, carriert,k,node] + SLACK[carrier,t,k,node] = 0 ∀ t,k")
  @constraint(cep.model, [carrier=intersect(set["carrier"]["all"],set["carrier"]["lost_load"]),t=set["time_T_period"]["all"], k=set["time_K"]["all"]], sum(cep.model[:GEN][tech,carrier,t,k,node] for node=set["nodes"]["all"], tech=set["tech"][carrier]) + sum(cep.model[:SLACK][carrier,t,k,node]*scale[:SLACK] for node=set["nodes"]["all"])/scale[:GEN] == 0)
  # or on copperplate
  push!(cep.info,"Σ_{tech,node}GEN[tech, carriert,k,node] = 0 ∀ t,k")
  @constraint(cep.model, [carrier=setdiff(set["carrier"]["all"],set["carrier"]["lost_load"]),t=set["time_T_period"]["all"], k=set["time_K"]["all"]], sum(cep.model[:GEN][tech,carrier,t,k,node] for node=set["nodes"]["all"], tech=set["tech"][carrier]) == 0)
  return cep
end

"""
    setup_opt_limit_emission!(cep::OptModelCEP, ts_data::ClustData, opt_data::OptDataCEP, scale::Dict{Symbol,Int}; limit::Dict{String,Dict} = Dict{String,Dict}("CO2"=>Dict{String,Number}("electricity"=>Inf)), lost_emission_cost::Dict{String,Number} = Dict{String,Number}("CO2"=>Inf))
Add emission limits constraints
"""
function setup_opt_limit_emission!(cep::OptModelCEP,
                            ts_data::ClustData,
                            opt_data::OptDataCEP,
                            scale::Dict{Symbol,Int};
                            limit_emission::Dict{String,Dict}=Dict{String,Dict}(),
                            lost_emission_cost::Dict{String,Number}=Dict{String,Number}())
  ## DATA ##
  set=cep.set
  #ts          Dict( tech-node ): t x k
  ts=ts_data.data
  #ts_weights: k - weight of each period:
  ts_weights=ts_data.weights
  #ts_deltas:  t x k - Δt of each segment x period
  ts_deltas=ts_data.delta_t

  ## EMISSIONS ##
  # Limit the Emissions with limit[impact] and allow lost emissions for all impacts that are both within the tech-group `limit` and `lost_emission`
  push!(cep.info,"ΣCOST_{account,tech}[account,impact_limit&lost_emission,tech] ≤ LE[impact] + ∑limit_emission[impactpercarrier] Σ_{node,t,k} ts[demand_carrier-node,t,k] ⋅ ts_weights[k] ⋅ ts_deltas[t,k]")
  @constraint(cep.model, [impact=intersect(set["impact"]["limit"],set["impact"]["lost_emission"])], sum(cep.model[:COST][account,impact,tech] for account=set["account"]["all"], tech=set["tech"]["all"])<= (cep.model[:LE][impact]*scale[:LE] +  sum(sum(limit_emission[impact][carrier]*ts["demand_"*carrier*"-"*node][t,k] for carrier=keys(limit_emission[impact]))*ts_deltas[t,k]*ts_weights[k] for t=set["time_T_period"]["all"], k=set["time_K"]["all"], node=set["nodes"]["all"]))/scale[:COST])
  # Limit the Emissions with limit[impact] and do NOT allow lost emissions for all impacts that are only within the tech-group `limit` and not `lost_emission`
  # Total demand can also be determined with the function get_total_demand() edit both in case of changes of e.g. ts_deltas
  push!(cep.info,"ΣCOST_{account,tech}[account,impact_limit|lost_emission,tech] ≤ limit_emission[impactpercarrier] ⋅ Σ_{node,t,k} ts[demand_carrier-node,t,k] ⋅ ts_weights[k] ⋅ ts_deltas[t,k]")
  @constraint(cep.model, [impact=setdiff(set["impact"]["limit"],set["impact"]["lost_emission"])], sum(cep.model[:COST][account,impact,tech] for account=set["account"]["all"], tech=set["tech"]["all"])<= sum(sum(limit_emission[impact][carrier]*ts["demand_"*carrier*"-"*node][t,k] for carrier=keys(limit_emission[impact]))*ts_deltas[t,k]*ts_weights[k] for t=set["time_T_period"]["all"], k=set["time_K"]["all"], node=set["nodes"]["all"])/scale[:COST])
  return cep
end

"""
     setup_opt_existing_infrastructure!(cep::OptModelCEP,ts_data::ClustData,opt_data::OptDataCEP, scale::Dict{Symbol,Int})
fixing existing infrastructure to CAP[tech, 'ex', node] and  TRANS[tech, 'ex', line] to the value defined in `nodes.csv` for all tech within the set set["tech"]["exist_inf"] as defined in the dictionary config['infrastructure']['existing']
"""
function setup_opt_existing_infrastructure!(cep::OptModelCEP,
                            ts_data::ClustData,
                            opt_data::OptDataCEP,
                            scale::Dict{Symbol,Int})
  ## DATA ##
  set=cep.set
  #`nodes::OptVariable`: nodes[tech,node] - OptDataCEPNode
  nodes = opt_data.nodes
  #`lines::OptVarible`: lines[tech,line] - OptDataCEPLine
  lines = opt_data.lines

  ## ASSIGN VALUES ##
  # Assign the existing capacity from the nodes table
  push!(cep.info,"CAP[tech, 'ex', node] = existing infrastructure ∀ node, tech ∈ tech_group_ex")
  push!(cep.info,"CAP[tech, 'ex', node] = 0  ∀ node, tech ∉ tech_group_ex")
  @constraint(cep.model, [node=set["nodes"]["all"], tech=intersect(set["tech"]["node"],set["tech"]["exist_inf"])], cep.model[:CAP][tech,"ex",node]==nodes[tech,node].power_ex/scale[:CAP])
  @constraint(cep.model, [node=set["nodes"]["all"], tech=setdiff(set["tech"]["node"],set["tech"]["exist_inf"])], cep.model[:CAP][tech,"ex",node]==0)
  if haskey(set["tech"],"transmission")
    push!(cep.info,"TRANS[tech, 'ex', line] = existing infrastructure ∀ tech ∈ tech_group_ex, line")
    push!(cep.info,"TRANS[tech, 'ex', line] = 0 ∀ tech ∉ tech_group_ex, line")
    @constraint(cep.model, [line=set["lines"]["all"], tech=intersect(set["tech"]["line"],set["tech"]["exist_inf"])], cep.model[:TRANS][tech,"ex",line]==lines[tech,line].power_ex/scale[:TRANS])
    @constraint(cep.model, [line=set["lines"]["all"], tech=setdiff(set["tech"]["line"],set["tech"]["exist_inf"])], cep.model[:TRANS][tech,"ex",line]==0)
  end
  return cep
end
"""
     setup_opt_limit_infrastructure!(cep::OptModelCEP,ts_data::ClustData,opt_data::OptDataCEP, scale::Dict{Symbol,Int})
limit infrastructure setup of CAP[tech, sum(infrastuct), node] and  TRANS[tech, infrastruct, line] to be smaller than the limit defined in `nodes.csv` for all tech within the set set["tech"]["limit"] as defined in the dictionary config['infrastructure']['limit']
"""
function setup_opt_limit_infrastructure!(cep::OptModelCEP,
                            ts_data::ClustData,
                            opt_data::OptDataCEP,
                            scale::Dict{Symbol,Int})
  ## DATA ##
  set=cep.set
  #`nodes::OptVariable`: nodes[tech,node] - OptDataCEPNode
  nodes = opt_data.nodes
  #`lines::OptVarible`: lines[tech,line] - OptDataCEPLine
  lines = opt_data.lines

  ## ASSIGN VALUES ##
  # Limit the capacity for each tech at each node with the limit provided in nodes table in column infrastruct
  push!(cep.info,"∑_{infrastuct} CAP[tech, infrastruct, node] <= limit infrastructure ∀ tech_n, node")
  @constraint(cep.model, [node=set["nodes"]["all"], tech=intersect(set["tech"]["node"],set["tech"]["limit"])], sum(cep.model[:CAP][tech,infrastruct,node] for infrastruct=set["infrastruct"]["all"]) <= nodes[tech,node].power_lim/scale[:CAP])
  if haskey(set["tech"],"transmission")
    push!(cep.info,"∑_{infrastuct} TRANS[tech, infrastruct, line] <= limit infrastructure ∀ tech_trans, line")
    @constraint(cep.model, [line=set["lines"]["all"], tech=intersect(set["tech"]["line"],set["tech"]["limit"])], sum(cep.model[:TRANS][tech,infrastruct,line] for infrastruct=set["infrastruct"]["all"]) <= lines[tech,line].power_lim/scale[:TRANS])
  end
  return cep
end

"""
     setup_opt_objective!(cep::OptModelCEP,ts_data::ClustData,opt_data::OptDataCEP, scale::Dict{Symbol,Int})
Calculate total system costs and set as objective
"""
function setup_opt_objective!(cep::OptModelCEP,
                            ts_data::ClustData,
                            opt_data::OptDataCEP,
                            scale::Dict{Symbol,Int};
                            lost_load_cost::Dict{String,Number}=Dict{String,Number}(),
                            lost_emission_cost::Dict{String,Number}=Dict{String,Number}())
  ## DATA ##
  set=cep.set

  ## OBJECTIVE ##
  # Minimize the total €-Costs s.t. the Constraints introduced above
  push!(cep.info,"min Σ_{account,tech}COST[account,'$(first(set["impact"]["mon"]))',tech] + Σ_{node,carrier_ll} LL[carrier,node] ⋅ lost_load_cost[carrier]) +  Σ_{impact_le} LE[impact] ⋅ lost_emission_cost[impact] st. above")
  @objective(cep.model, Min,  sum(cep.model[:COST][account,first(set["impact"]["mon"]),tech] for account=set["account"]["all"], tech=set["tech"]["all"]) + sum(cep.model[:LL][carrier,node]*lost_load_cost[carrier] for node=set["nodes"]["all"], carrier=set["carrier"]["lost_load"])*scale[:LL]/scale[:COST] + sum(cep.model[:LE][impact]*lost_emission_cost[impact] for impact=set["impact"]["lost_emission"])*scale[:LE]/scale[:COST])
  return cep
end

"""
     solve_opt_cep(cep::OptModelCEP,ts_data::ClustData,opt_data::OptDataCEP,config::Dict{String,Any})
solving the cep model and writing it's results and `limit` into an OptResult-Struct
"""
function solve_opt_cep(cep::OptModelCEP,
                            ts_data::ClustData,
                            opt_data::OptDataCEP,
                            config::Dict{String,Any})
  optimize!(cep.model)
  status=Symbol(termination_status(cep.model))
  scale=config["scale"]
  objective=objective_value(cep.model)*scale[:COST]
  total_demand=get_total_demand(cep,ts_data)
  variables=Dict{String,Any}()
  # cv - Cost variable, dv - design variable, which is used to fix variables in a dispatch model, ov - operational variable
  variables["COST"]=OptVariable(cep,:COST,"cv",scale)
  variables["CAP"]=OptVariable(cep,:CAP,"dv",scale)
  variables["GEN"]=OptVariable(cep,:GEN,"ov",scale)
  lost_load=0
  lost_emission=0
  if !isempty(config["lost_load_cost"])
    variables["SLACK"]=OptVariable(cep,:SLACK,"sv",scale)
    variables["LL"]=OptVariable(cep,:LL,"sv",scale)
    lost_load=sum(variables["LL"].data)
  end
  if !isempty(config["lost_emission_cost"])
    variables["LE"]=OptVariable(cep,:LE,"sv",scale)
    lost_emission=sum(variables["LE"].data)
  end
  if config["storage"]
    variables["INTRASTOR"]=OptVariable(cep,:INTRASTOR,"ov",scale)
    if config["seasonalstorage"]
      variables["INTERSTOR"]=OptVariable(cep,:INTERSTOR,"ov",scale)
    end
  end
  if config["transmission"]
    variables["TRANS"]=OptVariable(cep,:TRANS,"dv",scale)
    variables["FLOW"]=OptVariable(cep,:FLOW,"ov",scale)
  end
  get_met_cap_limit(cep, opt_data, variables)
  currency=variables["COST"].axes[2][1]
  if lost_load==0 && lost_emission==0
    config["print_flag"] && @info("Solved Scenario $(config["descriptor"]): "*String(status)*" min COST: $(round(objective,sigdigits=4)) [$currency] ⇨ $(round(objective/total_demand,sigdigits=4)) [$currency per MWh] s.t. $(text_limit_emission(config["limit_emission"]))")
  else
    cost=variables["COST"]
    config["print_flag"] && @info("Solved Scenario $(config["descriptor"]): "*String(status)*" min COST: $(round(sum(cost[:,axes(cost,"impact")[1],:]),sigdigits=4)) [$currency] ⇨ $(round(sum(cost[:,axes(cost,"impact")[1],:])/total_demand,sigdigits=4)) [$currency per MWh] with LL: $lost_load [MWh] s.t. $(text_limit_emission(config["limit_emission"])) + $(round(lost_emission/total_demand,sigdigits=4)) (violation) [kg-CO₂-eq. per MWh]")
  end
  info=Dict{String,Any}("total_demand"=>total_demand,"model"=>cep.info,)
  return OptResult(status,objective,variables,cep.set,config,info)
end


text:
# makes a good dir in default temp directory
let
  @pytest function(tempdir_fixture)
    @test isdir(tempdir_fixture)
    @test dirname(tempdir_fixture) == Base.tempdir()
  end
end

# deletes the tmpdir
let
  safe = Dict()
  @pytest function(tempdir_fixture)
    safe["tempdir_fixture"] = tempdir_fixture
  end
  @test !ispath(safe["tempdir_fixture"])
end


text:
using CoordinateSystems
using Base.Test
using ImmutableArrays

# write your own tests here
@test CoordinateSystems.normalize(Vector3(2.0,0.0,0.0)) == Vector3(1.0,0.0,0.0)

v0 = Vector3(0.0,0.0,0.0)
vx = Vector3(1.0,0.0,0.0)
vy = Vector3(0.0,1.0,0.0)
vz = Vector3(0.0,0.0,1.0)
vn = Vector3(1.0,2.0,3.0)

normal_cs = CoordinateSystem(v0,vz,vx)
@test to_global(normal_cs,vn) == vn
@test from_global(normal_cs,vn) == vn

flipxy_cs = CoordinateSystem(v0,vz,vy)
@test to_global(flipxy_cs,vn) == Vector3(-2.0,1.0,3.0)
@test from_global(flipxy_cs,vn) == Vector3(2.0,-1.0,3.0)

offsetx_cs = CoordinateSystem(vx,vx+vz,vx+vx)
@test to_global(offsetx_cs,vn) == Vector3(2.0,2.0,3.0)
@test from_global(offsetx_cs,vn) == Vector3(0.0,2.0,3.0)


text:
"""
    test(ref)

Determine 'ref' and print it to the console
"""

function test(ref)
	println("This is a function that prints $(ref)")
end


text:
"""
Sensitivity vector with QR
"""

"""
stsen(Q::Matrix{Tv}, pis::Vector{Tv}, b::Vector{Tv})

Get a sensitivity vector for stationary vector of CTMC.

Parameters:
- Q: CTMC Kernal
- pis: A stationary vector, i.e., pis * Q = 0
- b: A vector. In the case of first derivative, b = pis * Qdash where Qdash is the first derivate of Q. It may change when we want to the high-order derivatives.
Return value:
A tuple of
- x: sensitivity vector
"""

function stsen(Q::Matrix{Tv}, pis::Vector{Tv}, b::Vector{Tv})::Vector{Tv} where Tv
    m, n = size(Q)
    @assert m == n
    @assert ctmcstcheck(Q, pis)
    qm, rm = qr(Q')
    rm[m,n] = Tv(1)
    xx = rm \ (- qm' * b)
    xx - sum(xx) * pis
end