text:
using JuMP
using GLPK
using Plots  #; plotlyjs()
using Measures
using DataFrames
using NPZ
using CSV
using Dates
using DataFrames
using Statistics
using CartesianGeneticProgramming
using Cambrian
import Cambrian.mutate
using Plots, Measures
using StatsPlots
include("../utils.jl")
using Dierckx
using BSplineKit
using Distances

struct Month
    topo0
    topo1

    tp
    hs
    direction

    shoreline_change
    shoreline_change_class  # (landward, seaward)
end

function fit_poly_inf_norm(x,y,n)
	m = Model(GLPK.Optimizer)

    @variable(m, alpha[1:n+1])
	@variable(m, t >= 0.)
	@constraint(m, [i=1:length(y)], sum(alpha[j]*x[i]^(j-1) for j in 1:n+1) - y[i] <= t);
	@constraint(m, [i=1:length(y)], -t <= sum(alpha[j]*x[i]^(j-1) for j in 1:n+1) - y[i]);
	@objective(m, Min, t)

	optimize!(m)

    return JuMP.value.(alpha)
end

f2(x, args) = (args[2] * x) .+ args[1]

dataset_path = get_dataset_path(;dataset_name="series_dataset_path")
n_series = 6
series = []
count_months = 0
xs_csv_path = get_dataset_path(;dataset_name="xs_csv_path")
xs_csv = CSV.read(xs_csv_path, DataFrame)
xs = xs_csv.x

changes = []

for s in 1:n_series
    push!(series, [])
    series_dir = joinpath(dataset_path, string(s))

    files = readdir(series_dir)
    topo_files = filter(f-> !occursin("_waves.csv", f), files)
    topo_files = map(f->replace(f, "_topo.npy" => ""), topo_files)

    month_numbers = sort(map(f->parse(Int64, f), topo_files))
    n_months = length(month_numbers)
    global count_months += n_months

    println(s, " - ", n_months)

    for mi in 1:n_months
        i = month_numbers[mi]
        topo_path = string(dataset_path, '/', s, '/', i, "_topo.npy")
        if i == 312
            println("TEST")
            exit()
            continue
        end
        topo = npzread(topo_path)

        topo0 = topo[1, :]
        topo1 = topo[2, :]

        t0_index_last_positive = findlast(topo0.>0)
        t0_index_first_negative = findfirst(topo0.<0)

        t0_x_last_positive = xs[t0_index_last_positive]
        t0_x_first_negative = xs[t0_index_first_negative]
        t0_value_last_positive = topo0[t0_index_last_positive]
        t0_value_first_negative = topo0[t0_index_first_negative]

        t0_coeffs = fit_poly_inf_norm([t0_value_last_positive, t0_value_first_negative], [t0_x_last_positive, t0_x_first_negative], 1)
        t0_x_at_zero = f2([0], t0_coeffs)[1]

        t1_index_last_positive = findlast(topo1.>0)
        t1_index_first_negative = findfirst(topo1.<0)

        t1_x_last_positive = xs[t1_index_last_positive]
        t1_x_first_negative = xs[t1_index_first_negative]
        t1_value_last_positive = topo1[t1_index_last_positive]
        t1_value_first_negative = topo1[t1_index_first_negative]

        t1_coeffs = fit_poly_inf_norm([t1_value_last_positive, t1_value_first_negative], [t1_x_last_positive, t1_x_first_negative], 1)
        t1_x_at_zero = f2([0], t1_coeffs)[1]

        shoreline_change = t1_x_at_zero - t0_x_at_zero
        shoreline_change = norm_to_range(shoreline_change, -15.0, 15.0, -1.0, 1.0)
        if shoreline_change < 0
            class = (0, 1)  # (landward, seaward)
        else
            class = (1, 0)  # (landward, seaward)
        end

        waves_path = string(dataset_path, '/', s, '/', i, "_waves.csv")
        waves_csv = CSV.read(waves_path, DataFrame)

        tp_month = norm_to_range(waves_csv.tp, 0.0, 25.0, -1.0, 1.0)
        hs_month = norm_to_range(waves_csv.hs, 0.0, 10.0, -1.0, 1.0)
        direction_month = norm_to_range(waves_csv.direction, 0.0, 360.0, -1.0, 1.0)

        month = Month(topo0, topo1, tp_month, hs_month, direction_month, shoreline_change, class)
        push!(series[end], month)
    end
end

real_class = []
real_change = []
tp_min = []
tp_max = []
tp_avg = []
hs_min = []
hs_max = []
hs_avg = []
dir_min = []
dir_max = []
dir_avg = []

for si in 1:length(series)
    global correct
    global wrong
    months = series[si]
    for mi in 1:length(months)
        month = months[mi]
        push!(real_class, month.shoreline_change_class)
        push!(real_change, norm_to_range(month.shoreline_change, -1.0, 1.0, -15.0, 15.0))
        push!(tp_min, minimum(month.tp))
        push!(tp_max, maximum(month.tp))
        push!(tp_avg, mean(month.tp))
        push!(hs_min, minimum(month.hs))
        push!(hs_max, maximum(month.hs))
        push!(hs_avg, mean(month.hs))
        push!(dir_min, minimum(month.direction))
        push!(dir_max, maximum(month.direction))
        push!(dir_avg, mean(month.direction))
    end
end

println("Max change: ", maximum(real_change))
println("Min change: ", minimum(real_change))

# p_tp_min = scatter(tp_min, preds, xlabel="min(Tp)", ylabel="Pred (1=True, 0=False)")
# p_tp_max = scatter(tp_max, preds, xlabel="max(Tp)", ylabel="Pred (1=True, 0=False)")
# p_tp_avg = scatter(tp_avg, preds, xlabel="mean(Tp)", ylabel="Pred (1=True, 0=False)")
# tp_plot = plot(p_tp_min, p_tp_max, p_tp_avg, layout=(1, 3))

# p_hs_min = scatter(hs_min, preds, xlabel="min(Hs)", ylabel="Pred (1=True, 0=False)")
# p_hs_max = scatter(hs_max, preds, xlabel="max(Hs)", ylabel="Pred (1=True, 0=False)")
# p_hs_avg = scatter(hs_avg, preds, xlabel="mean(Hs)", ylabel="Pred (1=True, 0=False)")
# hs_plot = plot(p_hs_min, p_hs_max, p_hs_avg, layout=(1, 3))

# p_dir_min = scatter(dir_min, preds, xlabel="min(Dir)", ylabel="Pred (1=True, 0=False)")
# p_dir_max = scatter(dir_max, preds, xlabel="max(Dir)", ylabel="Pred (1=True, 0=False)")
# p_dir_avg = scatter(dir_avg, preds, xlabel="mean(Dir)", ylabel="Pred (1=True, 0=False)")
# dir_plot = plot(p_dir_min, p_dir_max, p_dir_avg, layout=(1, 3))

# display(plot(tp_plot, hs_plot, dir_plot, layout=(3, 1), size = (1600, 800)))

# display(plot(real_change, legend=false, ylabel="Shoreline change [m]", xlabel="Months", size=(1280, 300)))
# readline()
# exit()

################## dx correlations ##################
p_tp_min = scatter(tp_min, real_change, xlabel="min(Tp)", ylabel="Real change", title=string("Correlation: ", round(cor(tp_min, real_change), digits=2)))
p_tp_max = scatter(tp_max, real_change, xlabel="max(Tp)", ylabel="Real change", title=string("Correlation: ", round(cor(tp_max, real_change), digits=2)))
p_tp_avg = scatter(tp_avg, real_change, xlabel="mean(Tp)", ylabel="Real change", title=string("Correlation: ", round(cor(tp_avg, real_change), digits=2)))
tp_plot = plot(p_tp_min, p_tp_max, p_tp_avg, layout=(1, 3), bottom_margin=20mm)

p_hs_min = scatter(hs_min, real_change, xlabel="min(Hs)", ylabel="Real change", title=string("Correlation: ", round(cor(hs_min, real_change), digits=2)))
p_hs_max = scatter(hs_max, real_change, xlabel="max(Hs)", ylabel="Real change", title=string("Correlation: ", round(cor(hs_max, real_change), digits=2)))
p_hs_avg = scatter(hs_avg, real_change, xlabel="mean(Hs)", ylabel="Real change", title=string("Correlation: ", round(cor(hs_avg, real_change), digits=2)))
hs_plot = plot(p_hs_min, p_hs_max, p_hs_avg, layout=(1, 3), bottom_margin=20mm)

p_dir_min = scatter(dir_min, real_change, xlabel="min(dir)", ylabel="Real change", title=string("Correlation: ", round(cor(dir_min, real_change), digits=2)))
p_dir_max = scatter(dir_max, real_change, xlabel="max(dir)", ylabel="Real change", title=string("Correlation: ", round(cor(dir_max, real_change), digits=2)))
p_dir_avg = scatter(dir_avg, real_change, xlabel="mean(dir)", ylabel="Real change", title=string("Correlation: ", round(cor(dir_avg, real_change), digits=2)))
dir_plot = plot(p_dir_min, p_dir_max, p_dir_avg, layout=(1, 3))

dx_correlations_plot = plot(tp_plot, hs_plot, dir_plot, layout=(3, 1), size = (1600, 800), margin=10mm)

################## tp correlations ##################
p_hs_min = scatter(hs_min, tp_min, xlabel="min(Hs)", ylabel="min(Tp)", title=string("Correlation: ", round(cor(hs_min, tp_min), digits=2)))
p_hs_max = scatter(hs_max, tp_max, xlabel="max(Hs)", ylabel="max(Tp)", title=string("Correlation: ", round(cor(hs_max, tp_max), digits=2)))
p_hs_avg = scatter(hs_avg, tp_avg, xlabel="mean(Hs)", ylabel="mean(Tp)", title=string("Correlation: ", round(cor(hs_avg, tp_avg), digits=2)))
hs_plot = plot(p_hs_min, p_hs_max, p_hs_avg, layout=(1, 3), bottom_margin=20mm)

p_dir_min = scatter(dir_min, tp_min, xlabel="min(dir)", ylabel="min(Tp)", title=string("Correlation: ", round(cor(dir_min, tp_min), digits=2)))
p_dir_max = scatter(dir_max, tp_max, xlabel="max(dir)", ylabel="max(Tp)", title=string("Correlation: ", round(cor(dir_max, tp_max), digits=2)))
p_dir_avg = scatter(dir_avg, tp_avg, xlabel="mean(dir)", ylabel="mean(Tp)", title=string("Correlation: ", round(cor(dir_avg, tp_avg), digits=2)))
dir_plot = plot(p_dir_min, p_dir_max, p_dir_avg, layout=(1, 3))

tp_correlations_plot = plot(hs_plot, dir_plot, layout=(2, 1), size = (1600, 600), margin=10mm)
display(plot(tp_correlations_plot))
readline()

max_stars_repo_name: mahmoud-al-najar/cgp-morphcast

max_issues_repo_name: mahmoud-al-najar/cgp-morphcast

max_forks_repo_name: mahmoud-al-najar/cgp-morphcast


text:
export Math

module Math

const EPSILON = 0.0000001192092896

using ..Geometry

include("vector_2d.jl")

# Scratch vectors
v1 = Vector2D{Float64}()
v2 = Vector2D{Float64}()
v3 = Vector2D{Float64}()

include("transforms.jl")
include("interpolation.jl")

include("velocity.jl")

include("zoom_transform.jl")

end

ext: jl

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
@testset "Testing logic using strings" begin
    Base.one{T<:AbstractString}(::Type{T}) = "one "
    z = ["x " "y " "z "]

    want_1 = ["one " "x " "y " "z "]
    @test complete_polynomial(z, 1) == want_1

    want_2 = ["one " "x " "y " "z " "one x x " "one x y " "one x z " "one y y " "one y z " "one z z "]
    @test complete_polynomial(z, 2) == want_2
end

@testset "testing numerical routintes" begin
    z = rand(100, 10)
    the_ones = ones(size(z, 1))
    for d in 1:5
        n_comp = n_complete(size(z, 2), d)
        buff_d = Array(Float64, size(z, 1), n_comp)
        out_d = complete_polynomial(z, d)
        complete_polynomial!(z, d, buff_d)
        @test size(out_d, 1) == size(z, 1)
        @test size(out_d, 2) == n_comp
        @test out_d[:, 1] == the_ones
        @test out_d[:, 2:size(z, 2)+1] == z
        @test out_d == buff_d
    end

    z = rand(10, 3)
    out_2 = complete_polynomial(z, 2)
    # need to test columns size(z, 2) + 2:end
    @test out_2[:, 5] == z[:, 1].^2
    @test out_2[:, 6] == z[:, 1].*z[:, 2]
    @test out_2[:, 7] == z[:, 1].*z[:, 3]
    @test out_2[:, 8] == z[:, 2].*z[:, 2]
    @test out_2[:, 9] == z[:, 2].*z[:, 3]
    @test out_2[:, 10] == z[:, 3].*z[:, 3]
end


text:
using NamedTupleTools: namedtuple, merge

# -------------- Implementation of Product Basis

struct Product1pBasis{NB, TB <: Tuple, VALB} <: OneParticleBasis{Any}
   bases::TB
   indices::Vector{NTuple{NB, Int}}
   B_pool::VectorPool{VALB}
end

function Product1pBasis(bases)
   NB = length(bases)
   return Product1pBasis(bases, NTuple{NB, Int}[])
end

function Product1pBasis(bases, indices)
   # TODO: discuss whether to construct an optimal ordering, e.g.
   #       should the discrete bases come first once we implement the
   #       "strongzero" method?
   VT = _valtype(bases)
   Product1pBasis( tuple(bases...), indices, VectorPool{VT}() )
end

import Base.*
*(B1::OneParticleBasis, B2::OneParticleBasis) =
      Product1pBasis((B1, B2))
*(B1::Product1pBasis, B2::OneParticleBasis) =
      Product1pBasis((B1.bases..., B2))
*(B1::OneParticleBasis, B2::Product1pBasis) =
      Product1pBasis((B1, B2.bases...))
*(B1::Product1pBasis, B2::Product1pBasis) =
      Product1pBasis((B1.bases..., B2.bases...))

_numb(b::Product1pBasis{NB}) where {NB} = NB

Base.length(basis::Product1pBasis) = length(basis.indices)

# ------------------------- FIO CODES

==(B1::Product1pBasis, B2::Product1pBasis) =
      ( all(B1.bases .== B2.bases) &&
        B1.indices == B2.indices )

write_dict(B::Product1pBasis) =
      Dict("__id__" => "ACE_Product1pBasis",
            "bases" => write_dict.(B.bases),
          "indices" => B.indices )

function read_dict(::Val{:ACE_Product1pBasis}, D::Dict)
   bases = tuple( read_dict.(D["bases"])... )
   indices = [ tuple(v...) for v in D["indices"] ]
   return Product1pBasis(bases, indices)
end

# # -------- temporary hack for 1.6, should not be needed from 1.7 onwards

# function acquire_B!(basis::Product1pBasis, args...)
#    VT = valtype(basis, args...)
#    return acquire!(basis.B_pool, length(basis), VT)
# end

# function release_B!(basis::Product1pBasis, B)
#    return release!(basis.B_pool, B)
# end

# ------------------------------------

valtype(basis::Product1pBasis) = _valtype(basis.bases)

_valtype(bases) = promote_type(valtype.(bases)...)

valtype(basis::Product1pBasis, X::AbstractState) =
      promote_type(valtype.(basis.bases, Ref(X))...)

# valtype(basis::Product1pBasis, cfg::AbstractConfiguration) =
#       promote_type( valtype.(basis.bases, Ref(first(cfg)))... )

function valtype(basis::Product1pBasis, cfg::AbstractConfiguration)
   X = zero(eltype(cfg))
   return promote_type( valtype.(basis.bases, Ref(X))... )
end

gradtype(basis::Product1pBasis, cfg::Union{AbstractConfiguration, AbstractVector}) =
      gradtype(basis, zero(eltype(cfg)))

function gradtype(basis::Product1pBasis, X::AbstractState)
   VALT = valtype(basis, X)
   return dstate_type(VALT, X)
end

import Base.Cartesian: @nexprs

@generated function add_into_A!(A, basis::Product1pBasis{NB}, X) where {NB}
   quote
      @nexprs $NB i -> begin
         bas_i = basis.bases[i]
         B_i = acquire_B!(bas_i, X)
         evaluate!(B_i, bas_i, X)
      end
      for (iA, ϕ) in enumerate(basis.indices)
         t = one(eltype(A))
         @nexprs $NB i -> (t *= B_i[ϕ[i]])
         A[iA] += t
      end
      @nexprs $NB i -> release_B!(basis.bases[i], B_i)
      return nothing
   end
end

# this is a hack to resolve a method ambiguity.
add_into_A_dA!(A, dA, basis::Product1pBasis, X) =
               _add_into_A_dA!(A, dA, basis, X)
add_into_A_dA!(A, dA, basis::Product1pBasis, X, sym::Symbol) =
               _add_into_A_dA!(A, dA, basis, X, sym)

# args... could be a symbol to enable partial derivatives
# at the moment this is just passed into evaluate_ed!(...)
# to not evaluate 1p basis derivatives that aren't needed, but
# in the future a more complex construction could be envisioned that
# might considerably reduce the computational cost here...
@generated function _add_into_A_dA!(A, dA, basis::Product1pBasis{NB}, X,
                                   args...) where {NB}
   quote
      Base.Cartesian.@nexprs($NB, i -> begin   # for i = 1:NB
         bas_i = basis.bases[i]
         if !(bas_i isa Discrete1pBasis)
            # only evaluate basis gradients for a continuous basis
            B_i = acquire_B!(bas_i, X)
            dB_i = acquire_dB!(bas_i, X)
            Bt, dBt = evaluate_ed!(B_i, dB_i, bas_i, X, args...)
         else
            # we still need the basis values for the discrete basis though
            # TODO: maybe the d part should be a no-op and remove this
            # case distinction ...
            B_i = acquire_B!(bas_i, X)
            evaluate!(B_i, bas_i, X)
         end
      end)
      for (iA, ϕ) in enumerate(basis.indices)
         # evaluate A
         t = one(eltype(A))
         @nexprs $NB i -> (t *= B_i[ϕ[i]])
         A[iA] += t

         # evaluate dA
         # TODO: redo this with adjoints!!!!
         #     also reverse order of operations to make fewer multiplications!
         dA[iA] = zero(eltype(dA))
         Base.Cartesian.@nexprs($NB, a -> begin  # for a = 1:NB
            if !(basis.bases[a] isa Discrete1pBasis)
               dt = dB_a[ϕ[a]]
               Base.Cartesian.@nexprs($NB, b -> begin  # for b = 1:NB
                  if b != a
                     dt *= B_b[ϕ[b]]
                  end
               end)
               dA[iA] += dt
            end
         end)
      end
      Base.Cartesian.@nexprs($NB, i -> ( begin   # for i = 1:NB
         release_B!(bas_i, B_i)
         if !(basis.bases[i] isa Discrete1pBasis)
            release_dB!(bas_i, dB_i)
         end
      end))
      return nothing
   end
end

_symbols_prod(bases) = tuple(union( symbols.(bases)... )...)

symbols(basis::Product1pBasis) = _symbols_prod(basis.bases)

function indexrange(basis::Product1pBasis)
   allsyms = tuple(symbols(basis)...)
   rg = Dict{Symbol, Vector{Any}}([ sym => [] for sym in allsyms]...)
   for b in basis.bases
      rgb = indexrange(b)
      for sym in allsyms
         if haskey(rgb, sym)
            rg[sym] = union(rg[sym], rgb[sym])
         end
      end
   end
   # HACK: fix the m range based on the maximal l-range
   #       this needs to be suitably generalised if we have multiple
   #       (l, m) pairs, e.g. (l1, m1), (l2, m2)
   if haskey(rg, :m)
      maxl = maximum(rg[:l])
      rg[:m] = collect(-maxl:maxl)
   end

   # convert the range into a named tuple so that we remember the order!!
   return NamedTuple{allsyms}(ntuple(i -> rg[allsyms[i]], length(allsyms)))
end

isadmissible(b, basis::Product1pBasis) = all(isadmissible.(Ref(b), basis.bases))

function set_spec!(basis::Product1pBasis{NB}, spec) where {NB}
   empty!(basis.indices)
   for b in spec
      inds = ntuple(i -> get_index(basis.bases[i], b), NB)
      push!(basis.indices, inds)
   end
   return basis
end

get_spec(basis::Product1pBasis) = [ get_spec(basis, i) for i = 1:length(basis) ]

function get_spec(basis::Product1pBasis, i::Integer)
   inds = basis.indices[i]
   specs = get_spec.(basis.bases, inds)
   # TODO: here we should check that we are only merging compatible tuples,
   #       e.g. (n = 5, l = 2), (l = 2, m = -1) is ok
   #       but  (n = 5, l = 2), (l = 3, m = -1) is forbidden!
   return merge(specs...)
end

degree(b, basis::Product1pBasis) = sum( degree(b, B) for B in basis.bases )

degree(b::NamedTuple, basis::Product1pBasis, weight::Dict) =
      sum( degree(b, B, weight) for B in basis.bases )

# TODO: this looks like a horrible hack ...
function rand_radial(basis::Product1pBasis)
   for B in basis.bases
      if B isa ScalarACEBasis
         return rand_radial(B)
      end
   end
   return nothing
end

# --------------- AD codes

import ChainRules: rrule, NoTangent, ZeroTangent

_evaluate_bases(basis::Product1pBasis{NB}, X::AbstractState) where {NB} =
      ntuple(i -> evaluate(basis.bases[i], X), NB)

_evaluate_A(basis::Product1pBasis{NB}, BB) where {NB} =
      [ prod(BB[i][ϕ[i]] for i = 1:NB) for ϕ in basis.indices ]

evaluate(basis::Product1pBasis, X::AbstractState) =
      _evaluate_A(basis, _evaluate_bases(basis, X))

function _rrule_evaluate(basis::Product1pBasis{NB}, X::AbstractState,
                         w::AbstractVector{<: Number},
                         BB = _evaluate_bases(basis, X)) where {NB}
   VT = promote_type(valtype(basis, X), eltype(w))

   # dB = evaluate_d(basis, X)
   # return sum( (real(w) * real(db) + imag(w) * imag(db))
   #             for (w, db) in zip(w, dB) )

   # Compute the differentials for the individual sub-bases
   Wsub = ntuple(i -> zeros(VT, length(BB[i])), NB)
   for (ivv, vv) in enumerate(basis.indices)
      for t = 1:NB
         _A = one(VT)
         for s = 1:NB
            if s != t
               _A *= BB[s][vv[s]]
            end
         end
         Wsub[t][vv[t]] += w[ivv] * conj(_A)
      end
   end

   # now these can be propagated into the inner basis
   #  -> type instab to be fixed here
   g = sum( _rrule_evaluate(basis.bases[t], X, Wsub[t] )
            for t = 1:NB )
   return g
end

function rrule(::typeof(evaluate), basis::Product1pBasis, X::AbstractState)
   BB = _evaluate_bases(basis, X)
   A = _evaluate_A(basis, BB)
   return A,
      w -> (NoTangent(), NoTangent(), _rrule_evaluate(basis, X, w, BB))
end

#    function _rrule_evaluate(basis::Scal1pBasis, X::AbstractState,
#       w::AbstractVector{<: Number})
# x = _val(X, basis)
# a = _rrule_evaluate(basis.P, x, w)
# TDX = ACE.dstate_type(a, X)
# return TDX( NamedTuple{(_varsym(basis),)}( (a,) ) )
# end


text:
module JacobiAlgorithm

    struct EquilibriumProblem{T<:Real}
        nbz::Int64
        e_z::Function
        pmin::T
        pmax::T
        q_z::Vector{T}
        p0_z::Vector{T}
        function EquilibriumProblem(nbz::Int64,e_z::Function,pmin::T,pmax::T,q_z::Vector{T},p0_z::Vector{T})  where {T<:Real}
            this = new{T}(nbz,e_z,pmin,pmax,q_z,p0_z) #,NaN,[],[],[],NaN)
        end
    end

    function bisection(f::Function, a::T, b::T;
        tol::T=1e-5, maxiter::Int64=100) where T<:Real
        fa = f(a)
        fa*f(b) <= 0 || error("No real root in [a,b]")
        i = 0
        local c
        while b-a > tol
                i += 1
                i != maxiter || error("Max iteration exceeded")
                c = (a+b)/2
                fc = f(c)
            if fc == 0
                break
            elseif fa*fc > 0
                a = c  # Root is in the right half of [a,b].
                fa = fc
            else
                b = c  # Root is in the left half of [a,b].
            end
        end
        return c
    end

    function coordinate_update(problem::EquilibriumProblem{T}, z::Int64, price_z::Vector{T}) where T<:Real
        old_price = price_z[z]
        f(p) = problem.e_z(setindex!(price_z,p,z))[z] - problem.q_z[z]
        new_price = bisection(f, problem.pmin, problem.pmax,tol=1e-5,maxiter=1000)
        price_z[z] = old_price
        return new_price, new_price - old_price
    end

    function fJ_z(problem::EquilibriumProblem{T}, p_z::Vector{T}) where T<:Real
        nbz = size(p_z,1)
        pp_z = zeros(T, nbz)
        diffp_z = zeros(T, nbz)
        for z = 1:nbz
            pp_z[z], diffp_z[z] = coordinate_update(problem,z,p_z)
        end
        return pp_z, diffp_z
    end

    function fJ_threaded_z(problem::EquilibriumProblem{T}, p_z::Vector{T}) where T<:Real
        nbz = size(p_z,1)
        pp_z = zeros(T, nbz)
        diffp_z = zeros(T, nbz)
        Threads.@threads for z = 1:nbz
            p_z_safe = deepcopy(p_z)
            pp_z[z], diffp_z[z] = coordinate_update(problem,z,p_z_safe)
        end
        return pp_z, diffp_z
    end

    function solve(problem::EquilibriumProblem{T};
        maxit::Int64 = 10_000,
        method::Symbol = :jacobian,
        valtol::T = 1e-5,
        steptol::T = 1e-9,
        output::Int64 = 0) where T<:Real

        code = 0
        start_time = time()
        delta_z = zeros(T,problem.nbz)
        diffp_z = zeros(T,problem.nbz)
        iterations = 0
        p_z = deepcopy(problem.p0_z)
        if method != :jacobian
            error("Method not implemented")
        end
        for i = 1:maxit
            p_z, diffp_z = fJ_z(problem, p_z)
            delta_z = problem.e_z(p_z) .- problem.q_z'
            if output>1
                println("p = $p_z")
            end
            if maximum(abs.(delta_z)) < valtol
                code = 0
                iterations = i
                break
            end
            if maximum(abs.(diffp_z)) < steptol
                code = 1
                iterations = i
                break
            end
            code = 2
        end
        comp_time = time() - start_time
        if output > 0
            println("$(String(method)) method converged in $iterations iterations.")
            println("Value of p = $(p_z)")
            println("Value of p' - p = $(diffp_z)")
            println("Value of e(p)-q = $(delta_z)")
            println("Time elapsed: $comp_time")
            println("Code: $code")
        end
        return code
    end

    function solve_threaded(problem::EquilibriumProblem{T};
        maxit::Int64 = 10_000,
        method::Symbol = :jacobian,
        valtol::T = 1e-5,
        steptol::T = 1e-9,
        output::Int64 = 0) where T<:Real

        code = 0
        start_time = time()
        delta_z = zeros(T,problem.nbz)
        diffp_z = zeros(T,problem.nbz)
        iterations = 0
        p_z = deepcopy(problem.p0_z)
        if method != :jacobian
            error("Method not implemented")
        end
        for i = 1:maxit
            p_z, diffp_z = fJ_threaded_z(problem, p_z)
            delta_z = problem.e_z(p_z) .- problem.q_z'
            if output>1
                println("p = $p_z")
            end
            if maximum(abs.(delta_z)) < valtol
                code = 0
                iterations = i
                break
            end
            if maximum(abs.(diffp_z)) < steptol
                code = 1
                iterations = i
                break
            end
            code = 2
        end
        comp_time = time() - start_time
        if output > 0
            println("$(String(method)) method converged in $iterations iterations.")
            println("Value of p = $(p_z)")
            println("Value of p' - p = $(diffp_z)")
            println("Value of e(p)-q = $(delta_z)")
            println("Time elapsed: $comp_time")
            println("Code: $code")
        end
        return code
    end
end


text:
export M4, M5, M6

const σ4 = (2/3, 10/7/π, 1/π)
const σ5 = (1/24, 96/1199/π, 1/20/π)
const σ6 = (1/120, 7/478/π, 1/120/π)

const M4 = σ4 .* @piecewise_polynomial begin
    0.0 => p"0.0"
    1.0 => p"0.25*(2-q)^3 - (1 - q)^3"
    2.0 => p"0.25*(2-q)^3"
    _   => p"0.0"
end

const M5 = σ5 .* @piecewise_polynomial begin
    0.0 => p"0.0"
    0.5 => p"(2.5 - q)^4 - 5*(1.5 - q)^4 + 10*(0.5 - q)^4"
    1.5 => p"(2.5 - q)^4 - 5*(1.5 - q)^4"
    2.5 => p"(2.5 - q)^4"
    _   => p"0.0"
end

const M6 = σ6 .* @piecewise_polynomial begin
    0.0 => p"0.0"
    1.0 => p"(3 - q)^5 - 6*(2 - q)^5 + 15*(1 - q)^5"
    2.0 => p"(3 - q)^5 - 6*(2 - q)^5"
    3.0 => p"(3 - q)^5"
    _   => p"0.0"
end


text:
using BinaryBuilder

# Collection of sources required to build WCSLIB
sources = [
    "ftp://ftp.atnf.csiro.au/pub/software/wcslib/wcslib-5.18.tar.bz2" =>
    "b693fbf14f2553507bc0c72bca531f23c59885be8f7d3c3cb889a5349129509a",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd wcslib-5.18
# configure.ac uses the AC_CANONICAL_BUILD macro when it really should be using AC_CANONICAL_TARGET
chmod +w configure.ac
patch configure.ac <<EOF
23,24c23,24
< AC_CANONICAL_BUILD
< ARCH="\${build_cpu}-\$build_os"
---
> AC_CANONICAL_TARGET
> ARCH="\${target_cpu}-\$target_os"
241c241
<   case "\$build_os" in
---
>   case "\$target_os" in
249c249
<     case "\$build_cpu" in
---
>     case "\$target_cpu" in
EOF
autoconf
# mingw defines a wcsset function that clashes with a definition in libwcs, setting -DNO_OLDNAMES deletes the mingw one
if [[ "${target}" == *mingw* ]]; then
    ./configure --prefix=$prefix --host=$target --disable-fortran --without-cfitsio --without-pgplot --disable-utils CFLAGS=-DNO_OLDNAMES
else
    ./configure --prefix=$prefix --host=$target --disable-fortran --without-cfitsio --without-pgplot --disable-utils
fi
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    BinaryProvider.Linux(:i686, :glibc, :blank_abi),
    BinaryProvider.Linux(:x86_64, :glibc, :blank_abi),
    BinaryProvider.Linux(:armv7l, :glibc, :eabihf),
    BinaryProvider.Linux(:powerpc64le, :glibc, :blank_abi),
    BinaryProvider.MacOS(:x86_64, :blank_libc, :blank_abi),
    BinaryProvider.Windows(:i686, :blank_libc, :blank_abi),
    BinaryProvider.Windows(:x86_64, :blank_libc, :blank_abi)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libwcs", :libwcs)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "WCSLIB", sources, script, platforms, products, dependencies)

max_stars_repo_name: mweastwood/WCSLIBBuilder

max_issues_repo_name: mweastwood/WCSLIBBuilder

max_forks_repo_name: mweastwood/WCSLIBBuilder


text:
# vim: set wrap
#: by David Thierry 2021
#: Set all the eqns to hslice - 1
using JuMP
using SCIP
using DataFrames
using CSV

#: Data frames section
#: Load parameters

df_gas = DataFrame(CSV.File("../reg/gas_coeffs.csv"))
df_steam_full_power = DataFrame(CSV.File("../reg/steam_coeffs.csv"))
df_steam_full_steam = DataFrame(CSV.File("../reg/steam_coeffs_v3.csv"))

#: Load Prices
df_pow_c = DataFrame(CSV.File("../resources/FLECCSPriceSeriesData.csv"))
df_ng_c = DataFrame(CSV.File("../resources/natural_gas_price.csv"))

#: Assign parameters

#: Gas parameters
bPowGasTeLoad = df_gas[1, 2]
aPowGasTeLoad = df_gas[1, 3]

bFuelEload = df_gas[2, 2]
aFuelEload = df_gas[2, 3]
lbcoToTonneco = 0.4535924 / 1000
bEmissFactEload = df_gas[3, 2] * lbcoToTonneco
aEmissFactEload = df_gas[3, 3] * lbcoToTonneco

bPowHpEload = df_gas[4, 2] / 1000  #: To scale the kW to MW
aPowHpEload = df_gas[4, 3] / 1000

bPowIpEload = df_gas[5, 2] / 1000
aPowIpEload = df_gas[5, 3] / 1000

bAuxRateGas = df_gas[6, 2] / 1000
aAuxRateGas = df_gas[6, 3] / 1000

#: Steam params
bCcRebDutyEload = df_steam_full_steam[2, 2]
aCcRebDutyEload = df_steam_full_steam[2, 3]

#: Full power gives you the min steam
bDacSteaBaseEload = df_steam_full_power[3, 2]
aDacSteaBaseEload = df_steam_full_power[3, 3]

bSideSteaEload = df_steam_full_steam[3, 2] - df_steam_full_power[3, 2]
aSideSteaEload = df_steam_full_steam[3, 3] - df_steam_full_power[3, 3]

bAuxRateStea = df_steam_full_power[4, 2] / 1000
aAuxRateStea = df_steam_full_power[4, 3] / 1000

aLpSteaToPow = 78.60233832  # MMBtu to kwh

kwhToMmbtu = 3412.1416416 / 1e+06
#aSteaUseRateDacAir = 1944 * kwhToMmbtu
#aSteaUseRateDacFlue = 1944 * kwhToMmbtu
# 7 GJ/tonneCO
aSteaUseRateDacAir = 5 * (7 * 1e+06 / 3600) * kwhToMmbtu
aSteaUseRateDacFlue = 5 * (7 * 1e+06 / 3600) * kwhToMmbtu

aPowUseRateDacAir = 500 / 1000
aPowUseRateDacFlue = 250 / 1000
# 1 mmol/gSorb #
# per gCo/gSorb
gCogSorbRatio = 1e-03 * 44.0095

aSorbCo2CapFlue = 1 * gCogSorbRatio
aSorbCo2CapAir = gCogSorbRatio
#aSorbAmountFreshFlue = 176. * 10  # Tonne sorb
#aSorbAmountFreshAir = 176. * 10  # Tonne sorb

aSorbAmountFreshFlue = 176. * 10   # Tonne sorb (Max. heat basis)
aSorbAmountFreshAir = 3162.18 - 176. * 10  # Tonne sorb (Max. heat basis)

aCapRatePcc = 0.97
# 2.4 MJ/kg (1,050 Btu/lb) CO2 page 379/
#aSteaUseRatePcc = aSteaUseRateDacFlue * 0.2
#aSteaUseRatePcc = 2.4 * 1000 * 1000 / 3600 * kwhToMmbtu
#println(aSteaUseRatePcc)
# aPowUseRatePcc = 0.173514487  # MWh/tonneCoi2 (old)
aSteaUseRatePcc = 2.69 + 0.0218 + 0.00127 # MMBTU/tonne CO2 (trimeric)
println(aSteaUseRatePcc)
aPowUseRatePcc = 0.047 # MWh/tonne CO2 (trimeric)

#: Horizon Lenght
tHorz = 24 * 10 - 1

#: Slices per hour
hSlice = 4  # the number of slices of a given hour
# There are tHorz - 1 slices
# Each slice has hSlice points, but only states have the 0th

# If there's several slices in an hour we kind of need to divide the
# hourly-based quantities :(

# Correction for low load
slopeFactor =  (0.01/ 0.1) / ((.95 - .75) / (0.9 - 0.6))

# Dictionaries with the disjunction parameters

## Slope
aPowGt = Dict(0 => 1, 1 => aPowGasTeLoad * slopeFactor, 2 => aPowGasTeLoad)
aFuel = Dict(0 => 1, 1 => aFuelEload, 2 => aFuelEload)
aEmis = Dict(0 => 1, 1 => aEmissFactEload, 2 => aEmissFactEload)
aAuxGt = Dict(0 => 1,
    1 => aAuxRateGas * slopeFactor,
    2 => aAuxRateGas)
aPowHp = Dict(0 => 1, 1 => aPowHpEload * slopeFactor,
    2 => aPowHpEload)
aPowIp = Dict(0 => 1, 1 => aPowIpEload * slopeFactor,
    2 => aPowIpEload)
aCcReb = Dict(
    0 => 1,
    1 => aCcRebDutyEload * slopeFactor,
    2 => aCcRebDutyEload)
aDacSb = Dict(0 => 1, 1 => aDacSteaBaseEload * slopeFactor,
    2 => aDacSteaBaseEload)
aSideS = Dict(0 => 1, 1 => aSideSteaEload * slopeFactor, 2 => aSideSteaEload)
aAuxSt = Dict(0 => 1, 1 => aAuxRateStea * slopeFactor, 2 => aAuxRateStea)

ldics = ["aPowGt", "aFuel", "aEmis", "aAuxGt", "aPowHp", "aPowIp", "aCcReb", "aDacSb", "aSideS", "aAuxSt"]
println("Slopes")
for v in ldics
    s = Symbol(v)
    println(v, eval(s))
end

## b
bPowGt = Dict(0 => 0,
    1 => 60 * (-aPowGt[1] + aPowGt[2]) + bPowGasTeLoad,
    2 => bPowGasTeLoad)
bFuel = Dict(0 => 0, 1 => bFuelEload, 2 => bFuelEload)
bEmis = Dict(0 => 0, 1 => bEmissFactEload, 2 => bEmissFactEload)
bAuxGt = Dict(0 => 0,
    1 => 60 * (-aAuxGt[1] + aAuxGt[2]) + bAuxRateGas,
    2 => bAuxRateGas)
bPowHp = Dict(0 => 0,
    1 => 60 * (-aPowHp[1] + aPowHp[2]) + bPowHpEload,
    2 => bPowHpEload)
bPowIp = Dict(0 => 0,
    1 => 60 * (-aPowIp[1] + aPowIp[2]) + bPowIpEload,
    2 => bPowIpEload
    )
bCcReb = Dict(0 => 0,
    1 => 60 * (-aCcReb[1] + aCcReb[2]) + bCcRebDutyEload,
    2 => bCcRebDutyEload)
bDacSb = Dict(0 => 0,
    1 => 60 * (-aDacSb[1] + aDacSb[2]) + bDacSteaBaseEload,
    2 => bDacSteaBaseEload)
bSideS = Dict(0 => 0,
    1 => 60 * (-aSideS[1] + aSideS[2]) + bSideSteaEload,
    2 => bSideSteaEload)
bAuxSt = Dict(0 => 0,
    1 => 60 * (-aAuxSt[1] + aAuxSt[2]) + bAuxRateStea,
    2 => bAuxRateStea)

ldics = ["bPowGt", "bFuel", "bEmis", "bAuxGt", "bPowHp", "bPowIp", "bCcReb", "bDacSb", "bSideS", "bAuxSt"]
println("InterC")
for v in ldics
    s = Symbol(v)
    println(v, eval(s))
end

# USD/MWh
pow_price =(df_pow_c[!, "MiNg_150_ERCOT"])  # USD/MWh

# pow_price =(df_pow_c[!, "MiNg_150_PJM-W"])  # USD/MWh
#: Natural gas price
# 0.056 lb/cuft STP
#std_w_ng1000cuft = 0.056 * 1000
#cNgPerLbUsd = (3.5 / 1000) / 0.056

# Cost of natural gas
cNgPerMmbtu = 3.5

m = Model()

# aPowUseRateComp = 0.279751187  # MWh/tonneCo2
aPowUseRateComp = 0.076 # MWh/tonneCo2 (Trimeric)

# Other costs
cCostInvCombTurb = 1e+02
cCostInvSteaTurb = 1e+02
cCostInvTransInter = 1e+02
cCostInvPcc = 1e+02
cCostInvDac = 1e+03
cCostInvComp = 1e+01

# Cost parameters.
cEmissionPrice = 1.5e+02 # USD/tonne CO2
cCo2TranspPrice = 1e+00
pCo2Credit = 1e+00

#vCapCombTurb = 3.
vCapSteaTurb = 2.
vCapTransInter = 5.
vCapPcc = 20.
vCapComp = 1000.
# Capital Cost DAC USD/tCo2/yr
cCostInvDacUsdtCo2yr = 750
cCostFixedDacUsdtCo2yr = 25
cCostVariableDacUsdtCo2yr = 12

nMod = 2

# 0 -> off
# 1 -> warm-up
# 2 -> on

# 0 for off 2 for other ones.
extrPoint = Dict(0 => 0, 1 => 0:1, 2 => 0:1)
@variable(m, 0 <= lambda[0:tHorz, i=0:nMod, extrPoint[i]] <= 1)
# On/Off
@variable(m, y[0:tHorz, 0:nMod], Bin, start = 0)  # On and off
# Transition
@variable(m, z[0:tHorz, j1 = 0:nMod, j2 = 0:nMod], Bin, start = 0)

for i in 0:tHorz
    set_start_value(y[i, 1], 1)
end

@variable(m, xLoad[0:tHorz, 0:nMod])

@variable(m, 0 <= xeLoad[0:tHorz] <= 100)

@variable(m, 0 <= xPowGasTur[0:tHorz])
@variable(m, 0 <= xPowGross[0:tHorz])
@variable(m, 0 <= xPowOut[0:tHorz])

@variable(m, 0 <= xAuxPowGasT[0:tHorz])

# Steam Turbine
@variable(m, 0 <= xPowHp[0:tHorz])
@variable(m, 0 <= xPowIp[0:tHorz])
@variable(m, 0 <= xPowLp[0:tHorz])

@variable(m, 0 <= xFuel[0:tHorz])
@variable(m, 0 <= xCo2Fuel[0:tHorz])
@variable(m, 0 <= xDacSteaDuty[0:tHorz])

@variable(m, 0 <= xCcRebDuty[0:tHorz])
@variable(m, 0 <= xDacSteaBaseDuty[0:tHorz])

@variable(m, 0 <= xAllocSteam[0:tHorz])
@variable(m, 0 <= xSteaPowLp[0:tHorz])
@variable(m, 0 <= xSideSteaDac[0:tHorz])

#
@variable(m, 0 <= xPowSteaTur[0:tHorz])
@variable(m, 0 <= xAuxPowSteaT[0:tHorz])

# Disagretated variables

@variable(m, 0 <= xPowGt[0:tHorz, 0:nMod])
@variable(m, 0 <= xFuelNu[0:tHorz, 0:nMod])
@variable(m, 0 <= xEmis[0:tHorz, 0:nMod])
@variable(m, 0 <= xAuxGt[0:tHorz, 0:nMod])
@variable(m, 0 <= xPowHpNu[0:tHorz, 0:nMod])
@variable(m, 0 <= xPowIpNu[0:tHorz, 0:nMod])
@variable(m, 0 <= xPccReb[0:tHorz, 0:nMod])
@variable(m, 0 <= xDacSb[0:tHorz, 0:nMod])
@variable(m, 0 <= xAllocS[0:tHorz, 0:nMod])
@variable(m, 0 <= xAuxSt[0:tHorz, 0:nMod])

#@constraint(m, powgtDisEq[i = 0:tHorz, k = 0:nMod],
#    xPowGt[i, k] = aPowGt[k] * xLoad[i, k] + bPowGt[k])
#@constraint(m, powGasTur[i = 0:tHorz],
#            xPowGasTur[i] == sum(xPowGt[i, k] for k in 0:nMod)
#           )

# Pcc
#@variable(m, 0 <= xCo2CapPcc[0:tHorz - 1] <= vCapPcc)
@variable(m, 0 <= xCo2CapPcc[0:tHorz])
@variable(m, 0 <= xSteaUsePcc[0:tHorz])
@variable(m, 0 <= xPowUsePcc[0:tHorz])
@variable(m, 0 <= xCo2PccOut[0:tHorz])
#@variable(m, 0 <= vCo2PccVent[0:tHorz - 1] <= 0.1)
@variable(m, 0 <= vCo2PccVent[0:tHorz])
#vCo2PccVent = 0.0
@variable(m, 0 <= xCo2DacFlueIn[0:tHorz])
@variable(m, 0 <= xPccSteaSlack[0:tHorz])

# Dac-Flue
@variable(m, 0 <= xA0Flue[0:tHorz, 0:hSlice-1]) # Kind of input

@variable(m, 0 <= xA1Flue[0:tHorz, 0:hSlice]) # State

@variable(m, 0 <= xR0Flue[0:tHorz, 0:hSlice-1])  # Kind of input

@variable(m, 0 <= xR1Flue[0:tHorz, 0:hSlice])  # State
@variable(m, 0 <= xFflue[0:tHorz, 0:hSlice])  # State
@variable(m, 0 <= xSflue[0:tHorz, 0:hSlice])  # State

@variable(m, 0 <= xCo2CapDacFlue[0:tHorz, 0:hSlice-1])
@variable(m, 0 <= xSteaUseDacFlue[0:tHorz, 0:hSlice-1])
@variable(m, 0 <= xPowUseDacFlue[0:tHorz, 0:hSlice-1])
@variable(m, 0 <= xCo2DacVentFlue[0:tHorz])

# Dac-Air
@variable(m, 0 <= xA0Air[0:tHorz, 0:hSlice-1]) # Kind of input

@variable(m, 0 <= xA1Air[0:tHorz, 0:hSlice]) # State
@variable(m, 0 <= xA2Air[0:tHorz, 0:hSlice]) # State

@variable(m, 0 <= xR0Air[0:tHorz, 0:hSlice-1])  # Kind of input

@variable(m, 0 <= xR1Air[0:tHorz, 0:hSlice])  # State
@variable(m, 0 <= xFair[0:tHorz, 0:hSlice])  # State
@variable(m, 0 <= xSair[0:tHorz, 0:hSlice])  # State

@variable(m, 0 <= xCo2CapDacAir[0:tHorz, 0:hSlice-1])
@variable(m, 0 <= xSteaUseDacAir[0:tHorz, 0:hSlice-1])
@variable(m, 0 <= xPowUseDacAir[0:tHorz, 0:hSlice-1])

@variable(m, 0 <= xDacSteaSlack[0:tHorz])
# DAC hourly capacity
#

# CO2 compression
@variable(m, 0 <= xCo2Comp[0:tHorz])
@variable(m, 0 <= xPowUseComp[0:tHorz])
#@variable(m, 0 <= vCapComp)
@variable(m, xCo2Vent[0:tHorz])  # This used to be only positive.

@variable(m, 0 <= xAuxPow[0:tHorz])

# Constraints
# Op mode
# Down times
# Up times

# Disjunction 0 (off)

extreme_d_0 = [0]
extreme_d_1 = [20.0, 50.0]
extreme_d_2 = [60.0, 100.0]

ep_list = [extreme_d_0, extreme_d_1, extreme_d_2]

# Convex conbination
@constraint(m, cConvxEq[i = 0:tHorz, mod = 0:nMod],
    sum(lambda[i, mod, k] * ep_list[mod + 1][k + 1] for k in extrPoint[mod]) ==
    xLoad[i, mod]  # There is only a single extreme
    )

# jth mode
# Lambda constraint
@constraint(m, xLambdaEq[i = 0:tHorz, mod = 0:nMod],
    sum(lambda[i, mod, k] for k in extrPoint[mod]) == y[i, mod]
    )

# Big M
@constraint(m, bmconEq[i = 0:tHorz, mod = 0:nMod],
    xLoad[i, mod] <= maximum(ep_list[mod + 1])* y[i, mod]
    )

# Overall Disjunction

@constraint(m, gasLeq[i = 0:tHorz],
    xeLoad[i] == sum(xLoad[i, mod] for mod in 0:nMod)
    )

@constraint(m, sumyEq[i = 0:tHorz],
    sum(y[i, mod] for mod in 0:nMod) == 1
    )

# Switch

@constraint(m, switchCon[i = 1:tHorz, mod = 0:nMod],
    sum(z[i, k, mod] for k in 0:nMod) -
    sum(z[i, mod, k] for k in 0:nMod)
    == y[i, mod] - y[i - 1, mod])

# Forbidden
# from off to full
@constraint(m, forbConOffFullEq[i = 0:tHorz],
    z[i, 0, 2] == 0)
# from full to warm-up
@constraint(m, forbConEq2[i = 0:tHorz],
    z[i, 2, 1] == 0)
# from warm-up to off
#@constraint(m, forbConwarmpuoff[i = 0:tHorz],
#    z[i, 1, 0] == 0)

# Ramping Constraints
rup = [0, 15, 30]
rdown = [0, 15, 30]

@constraints(m, begin
    rampup1eq[i=1:tHorz], xLoad[i, 1] - xLoad[i-1, 1] <= rup[2] * y[i, 1]
    rampup2eq[i=1:tHorz], xLoad[i, 2] - xLoad[i-1, 2] <= rup[3] * y[i, 2]
end)

#@constraints(m, begin
    #rampdo1eq[i=1:tHorz], xLoad[i - 1, 1] - xLoad[i, 1] <= rdown[2] * y[i, 1]
    #rampdo2eq[i=1:tHorz], xLoad[i - 1, 2] - xLoad[i, 2] <= rdown[3] * y[i, 2]
#end)

# Minimum stay
KminOff = [[0, 36, 0], [0, 0, 0], [48, 0, 0]]

# off to on(1)
@constraint(m, minstay01[i = 1:tHorz],
    y[i, 1] >= sum(z[i - k, 0, 1] for k in 0:(24-1) if (i-k) >= 0)
    )
# on(2) to off
@constraint(m, minstay20[i = 1:tHorz],
    y[i, 0] >= sum(z[i - k, 2, 0] for k in 0:(48-1) if (i-k) >= 0)
    )

# Disagregated variables constraints
@constraint(m, powgtdeq[i = 0:tHorz, mod = 0:nMod],
    xPowGt[i, mod] == aPowGt[mod] * xLoad[i, mod] + bPowGt[mod] * y[i, mod]
    )
@constraint(m, fueldeq[i = 0:tHorz, mod = 0:nMod],
    xFuelNu[i, mod] == aFuel[mod] * xLoad[i, mod] + bFuel[mod] * y[i, mod]
    )
@constraint(m, emissdeq[i = 0:tHorz, mod = 0:nMod],
    xEmis[i, mod] == aEmis[mod] * xLoad[i, mod] + bEmis[mod] * y[i, mod]
    )
@constraint(m, auxgtdeq[i = 0:tHorz, mod = 0:nMod],
    xAuxGt[i, mod] == aAuxGt[mod] * xLoad[i, mod] + bAuxGt[mod] * y[i, mod]
    )
@constraint(m, powhpdeq[i = 0:tHorz, mod = 0:nMod],
    xPowHpNu[i, mod] == aPowHp[mod] * xLoad[i, mod] + bPowHp[mod] * y[i, mod]
    )
@constraint(m, powipdeq[i = 0:tHorz, mod = 0:nMod],
    xPowIpNu[i, mod] == aPowIp[mod] * xLoad[i, mod] + bPowIp[mod] * y[i, mod]
    )
@constraint(m, pccrebdeq[i = 0:tHorz, mod = 0:nMod],
    xPccReb[i, mod] == aCcReb[mod] * xLoad[i, mod] + bCcReb[mod] * y[i, mod]
    )
@constraint(m, dacsbtdeq[i = 0:tHorz, mod = 0:nMod],
    xDacSb[i, mod] == aDacSb[mod] * xLoad[i, mod] + bDacSb[mod] * y[i, mod]
    )
@constraint(m, allocstdeq[i = 0:tHorz, mod = 0:nMod],
    xAllocS[i, mod] == aSideS[mod] * xLoad[i, mod] + bSideS[mod] * y[i, mod]
    )
@constraint(m, auxstdeq[i = 0:tHorz, mod = 0:nMod],
    xAuxSt[i, mod] == aAuxSt[mod] * xLoad[i, mod] + bAuxSt[mod] * y[i, mod]
    )

# Gas Turbine
@constraint(m, powGasTur[i = 0:tHorz],
            xPowGasTur[i] == sum(xPowGt[i, mod] for mod in 0:nMod)
           )
#
@constraint(m, fuelEq[i = 0:tHorz],
            xFuel[i] == sum(xFuelNu[i, mod] for mod in 0:nMod)
           )
#
@constraint(m, co2FuelEq[i = 0:tHorz],
            xCo2Fuel[i] == sum(xEmis[i, mod] for mod in 0:nMod)
           )

@constraint(m, auxPowGasT[i = 0:tHorz],
            xAuxPowGasT[i] == sum(xAuxGt[i, mod] for mod in 0:nMod)
           )
#
# Steam
#
@constraint(m, powHpEq[i = 0:tHorz],
            xPowHp[i] == sum(xPowHpNu[i, mod] for mod in 0:nMod)
           )
#
@constraint(m, powIpEq[i = 0:tHorz],
            xPowIp[i] == sum(xPowIpNu[i, mod] for mod in 0:nMod)
           )

#
@constraint(m, powLpEq[i = 0:tHorz],
            xPowLp[i] == xSteaPowLp[i] * aLpSteaToPow / 1000
           )
#
@constraint(m, powerSteaEq[i = 0:tHorz],
            xPowSteaTur[i] ==
            xPowHp[i] + xPowIp[i] + xPowLp[i]
           )

@constraint(m, ccRebDutyEq[i = 0:tHorz],
            xCcRebDuty[i] == sum(xPccReb[i, mod] for mod in 0:nMod)
           )

@constraint(m, dacSteaDutyEq[i = 0:tHorz],
            xDacSteaBaseDuty[i] == sum(xDacSb[i, mod] for mod in 0:nMod)
           )

@constraint(m, sideSteaEloadEq[i = 0:tHorz],
            xAllocSteam[i] == sum(xAllocS[i, mod] for mod in 0:nMod)
           )

@constraint(m, sideSteaRatioEq[i = 0:tHorz],
            xAllocSteam[i] == xSideSteaDac[i] + xSteaPowLp[i]
           )

@constraint(m, availSteaDacEq[i = 0:tHorz],
            xDacSteaDuty[i] == xDacSteaBaseDuty[i] + xSideSteaDac[i]
           )

@constraint(m, auxPowSteaTEq[i = 0:tHorz],
            xAuxPowSteaT[i] == sum(xAuxSt[i, mod] for mod in 0:nMod)
           )

# PCC
#
#@constraint(m, co2CapPccEq[i = 0:tHorz - 1],
#xCo2CapPcc[i] == aCo2PccCapRate * xCo2Fuel[i])
@constraint(m, co2CapPccEq[i = 0:tHorz],
            xCo2CapPcc[i] == aCapRatePcc * xCo2Fuel[i])
#
@constraint(m, co2PccOutEq[i = 0:tHorz],
            xCo2PccOut[i] == xCo2Fuel[i] - xCo2CapPcc[i])
#
@constraint(m, co2DacFlueInEq[i = 0:tHorz],
            xCo2DacFlueIn[i] == xCo2PccOut[i] - vCo2PccVent[i])
#
# @constraint(m, co2CapPccIn[i = 0:tHorz - 1], xCo2CapPcc[i] <= vCapPcc)
# Dav: Sometimes there is not enough steam, so we have to relax this constraint
@constraint(m, steamUsePccEq[i = 0:tHorz],
            xSteaUsePcc[i] <= aSteaUseRatePcc * xCo2CapPcc[i])
#
@constraint(m, powerUsePccEq[i = 0:tHorz],
            xPowUsePcc[i] == aPowUseRatePcc * xCo2CapPcc[i]
           )

@constraint(m, pccSteaSlack[i = 0:tHorz],
            xPccSteaSlack[i] == xCcRebDuty[i] - xSteaUsePcc[i])

# DAC-Flue
# Flue gas takes 15 minutes to saturation?
#: "State equation"
@constraint(m, a1dFlueEq[i = 0:tHorz, j=1:hSlice],
            xA1Flue[i, j] == xA0Flue[i, j-1]
           )

#: "State equation"
@constraint(m, aRdFlueEq[i = 0:tHorz, j=1:hSlice],
            xR1Flue[i, j] == xR0Flue[i, j-1]
           )
#: "State equation"
@constraint(m, storeFflueeq[i = 0:tHorz, j = 1:hSlice],
            xFflue[i, j] == xFflue[i, j-1] - xA0Flue[i, j-1] + xR1Flue[i, j-1]
           )
#: "State equation"
@constraint(m, storeSflueeq[i = 0:tHorz, j = 1:hSlice],
            xSflue[i, j] == xSflue[i, j-1] - xR0Flue[i, j-1] + xA1Flue[i, j-1]
           )
# Initial conditions
@constraint(m, icXa1FlueEq, xA1Flue[0, 0] == 0.)
@constraint(m, icAR1FlueEq, xR1Flue[0, 0] == 0.)
@constraint(m, capDacFlueEq, xFflue[0, 0] == aSorbAmountFreshFlue)
@constraint(m, icSsFlueEq, xSflue[0, 0] == 0.)
# End-point constraints we need to get rid of them and then put them back
@constraint(m, endXa1FlueEq, xA1Flue[tHorz, hSlice] == 0.)
@constraint(m, endAR1FlueEq, xR1Flue[tHorz, hSlice] == 0.)
@constraint(m, endDacFlueEq, xFflue[tHorz, hSlice] == aSorbAmountFreshFlue)
@constraint(m, endSsFlueEq, xSflue[tHorz, hSlice] == 0.)
#
#These dac related variables must start at 0 and end at 1-hslice
@constraint(m, co2CapDacFlueEq[i = 0:tHorz, j = 0:hSlice-1],
            xCo2CapDacFlue[i, j] ==
            #aSorbCo2CapFlue * xR1Flue[i, j]
            aSorbCo2CapFlue * xA1Flue[i, j]
           )
#
@constraint(m, steamUseDacFlueEq[i = 0:tHorz, j = 0:hSlice-1],
            xSteaUseDacFlue[i, j] ==
            #aSteaUseRateDacFlue * xCo2CapDacFlue[i, j]
            aSteaUseRateDacFlue * aSorbCo2CapFlue * xR1Flue[i, j]
           )
#
@constraint(m, powUseDacFlueEq[i = 0:tHorz, j = 0:hSlice-1],
            xPowUseDacFlue[i, j] == aPowUseRateDacFlue * xCo2CapDacFlue[i, j]
           )

# DAC-Air
# Bluntly assume we can just take CO2 from air in pure form.
# "State equation"
@constraint(m, a1dAirEq[i = 0:tHorz, j = 1:hSlice],
            xA1Air[i, j] == xA0Air[i, j - 1]
           )
# "State equation"
@constraint(m, a2dAirEq[i = 0:tHorz, j = 1:hSlice],
            xA2Air[i, j] == xA1Air[i, j - 1]
           )
# "State equation"
@constraint(m, aRdAirEq[i = 0:tHorz, j = 1:hSlice],
            xR1Air[i, j] == xR0Air[i, j - 1]
           )
# "State equation"
@constraint(m, storeFairEq[i = 0:tHorz, j = 1:hSlice],
            xFair[i, j] == xFair[i, j-1] - xA0Air[i, j-1] + xR1Air[i, j-1]
           )
# "State equation"
@constraint(m, storeSaireq[i = 0:tHorz, j = 1:hSlice],
            xSair[i, j] == xSair[i, j-1] - xR0Air[i, j-1] + xA2Air[i, j-1]
           )

# Initial conditions - Air
#@constraint(m, capDacAirEq, xFair[0] == xSorbFreshAir)
@constraint(m, capDacAirEq, xFair[0, 0] == aSorbAmountFreshAir)
@constraint(m, icA1AirEq, xA1Air[0, 0] == 0.)
@constraint(m, icA2AirEq, xA2Air[0, 0] == 0.)
@constraint(m, icAR1AirEq, xR1Air[0, 0] == 0.)
@constraint(m, icSsAirEq, xSair[0, 0] == 0.)
# End-point conditions - Air

@constraint(m, endDacAirEq, xFair[tHorz, hSlice] == aSorbAmountFreshAir)
@constraint(m, endA1AirEq, xA1Air[tHorz, hSlice] == 0.)
@constraint(m, endA2AirEq, xA2Air[tHorz, hSlice] == 0.)
@constraint(m, endAR1AirEq, xR1Air[tHorz, hSlice] == 0.)
@constraint(m, endSsAirEq, xSair[tHorz, hSlice] == 0.)

#@constraint(m, endDacAirEq, xFair[0] == xSorbFreshAir)
#@constraint(m, endDacAirEq, xFair[0, 0] == aSorbAmountFreshAir)
#@constraint(m, endA1AirEq, xA1Air[0, 0] == 0.)
#@constraint(m, endA2AirEq, xA2Air[0, 0] == 0.)
#@constraint(m, endAR1AirEq, xR1Air[0, 0] == 0.)
#@constraint(m, endSsAirEq, xSair[0, 0] == 0.)

#
@constraint(m, co2CapDacAirEq[i = 0:tHorz, j=0:hSlice-1],
            xCo2CapDacAir[i, j] ==
            #aSorbCo2CapAir * xR1Air[i, j]
            (aSorbCo2CapAir * xA1Air[i, j]/2 + aSorbCo2CapAir * xA2Air[i, j]/2)
           )
#
@constraint(m, steamUseDacAirEq[i = 0:tHorz, j=0:hSlice-1],
            xSteaUseDacAir[i, j] ==
            #aSteaUseRateDacAir * xCo2CapDacAir[i, j]
            aSteaUseRateDacAir * aSorbCo2CapAir * xR1Air[i, j]
           )
#
@constraint(m, powUseDacAirEq[i = 0:tHorz, j=0:hSlice-1],
            xPowUseDacAir[i, j] == aPowUseRateDacAir * xCo2CapDacAir[i, j]
           )
# We need the integral over the whole time for DAC
@constraint(m, dacSteaSlackEq[i = 0:tHorz],
            xDacSteaSlack[i] == xDacSteaDuty[i]
            - sum(xSteaUseDacFlue[i, j] for j in 0:hSlice - 1)
            - sum(xSteaUseDacAir[i, j] for j in 0:hSlice - 1)
           )

# Equal to the amount vented, at least in flue mode.
# Integral for DAC
@constraint(m, co2DacFlueVentEq[i = 0:tHorz],
            xCo2DacVentFlue[i] == xCo2DacFlueIn[i]
            - sum(xCo2CapDacFlue[i, j] for j in  0:hSlice-1)
           )

# Co2 Compression
#
# We need the integral over the whole time for DAC
@constraint(m, co2CompEq[i = 0:tHorz],
            xCo2Comp[i] == xCo2CapPcc[i]
            + sum(xCo2CapDacFlue[i, j] for j in 0:hSlice-1)
            + sum(xCo2CapDacAir[i, j] for j in 0:hSlice-1)
           )
#
@constraint(m, powUseCompEq[i = 0:tHorz],
            xPowUseComp[i] == aPowUseRateComp * xCo2Comp[i]
           )
#
@constraint(m, co2VentEq[i = 0:tHorz],
            xCo2Vent[i] == vCo2PccVent[i]
            + xCo2DacVentFlue[i] - sum(xCo2CapDacAir[i, j]
                for j in 0:hSlice-1)
           )

## Overall

#
@constraint(m, powGrossEq[i = 0:tHorz],
            xPowGross[i] == xPowGasTur[i] + xPowSteaTur[i]
           )
@constraint(m, auxPowEq[i = 0:tHorz],
            xAuxPow[i] == xAuxPowGasT[i] + xAuxPowSteaT[i])

@constraint(m, powOutEq[i = 0:tHorz],
            xPowOut[i] == xPowGross[i]
            - xPowUsePcc[i]
            - sum(xPowUseDacFlue[i, j] for j in 0:hSlice-1)
            - sum(xPowUseDacAir[i, j] for j in 0:hSlice-1)
            - xPowUseComp[i]
            - xAuxPow[i]
           )

# Continuity of states
@constraint(m,
            contxfflue[i = 1:tHorz], xFflue[i, 0] == xFflue[i - 1, hSlice])
@constraint(m,
            conta1flue[i = 1:tHorz], xA1Flue[i, 0] == xA1Flue[i - 1, hSlice])

@constraint(m,
            contcxsflue[i = 1:tHorz], xSflue[i, 0] == xSflue[i - 1, hSlice])
@constraint(m,
            contr1flue[i = 1:tHorz], xR1Flue[i, 0] == xR1Flue[i - 1, hSlice])

@constraint(m,
            contxfair[i = 1:tHorz], xFair[i, 0] == xFair[i - 1, hSlice])
@constraint(m,
            conta1air[i = 1:tHorz], xA1Air[i, 0] == xA1Air[i - 1, hSlice])
@constraint(m,
            conta2air[i = 1:tHorz], xA2Air[i, 0] == xA2Air[i - 1, hSlice])
@constraint(m,
            contcxsair[i = 1:tHorz], xSair[i, 0] == xSair[i - 1, hSlice])
@constraint(m,
            contr1air[i = 1:tHorz], xR1Air[i, 0] == xR1Air[i - 1, hSlice])

shutdowncost = cNgPerMmbtu * (bFuelEload + 60 * aFuelEload)
shutdownfromwarm = cNgPerMmbtu * (bFuelEload + 60 * aFuelEload) * 24
startupcost = cNgPerMmbtu * (bFuelEload + 60 * aFuelEload) / 2
# Objective function expression
@expression(m, eObjfExpr,
    sum(cNgPerMmbtu * xFuel[i]
        + cEmissionPrice * xCo2Vent[i]
        + cCo2TranspPrice * xCo2Comp[i]
        - pow_price[i + 1] * xPowOut[i]
        for i in 0:tHorz)
        +
    sum(shutdowncost * z[i, 2, 0] + startupcost * z[i, 0, 1]
        + shutdownfromwarm * z[i, 1, 0]
        for i in 0:tHorz)
    )

@objective(m, Min, eObjfExpr)

print("The number of variables\t")
println(num_variables(m))
print("The number of constraints\n")

n = 0
for i in list_of_constraint_types(m)
    global n
    println(num_constraints(m, i[1], i[2]))
    n += num_constraints(m, i[1], i[2])
end

println()

# Set optimizer options
set_optimizer(m, SCIP.Optimizer)
#set_optimizer_attribute(m, "LogLevel", 3)
#set_optimizer_attribute(m, "PresolveType", 1)

optimize!(m)
println(termination_status(m))

#f = open("model.lp", "w")
#print(f, m)
#close(f)

#write_to_file(m, "lp_mk0.mps")
#write_to_file(m, "lp_mk10.lp", format=MOI.FileFormats.FORMAT_LP)

#format::MOI.FileFormats.FileFormat = MOI.FileFormats.FORMAT_AUTOMATIC

# Co2 Data Frame
df_co = DataFrame(Symbol("Co2Fuel") => Float64[], # Pairs.
                  Symbol("Co2CapPcc") => Float64[],
                  Symbol("Co2PccOut") => Float64[],
                  Symbol("vCo2PccVent") => Float64[],
                  Symbol("Co2DacFlueIn") => Float64[],
                  Symbol("Co2CapDacFlue") => Float64[],
                  Symbol("Co2CapDacAir") => Float64[],
                  Symbol("Co2DacVentFlue") => Float64[],
                  Symbol("Co2Vent") => Float64[],
                 )

# Co2 / hSlice
for i in 0:tHorz
    co2fuel = value(xCo2Fuel[i])
    co2pcc = value(xCo2CapPcc[i])
    co2pccout = value(xCo2PccOut[i])
    co2pccvent = value(vCo2PccVent[i])
    co2dacfluein = value(xCo2DacFlueIn[i])
    co2dacflue = sum(value(xCo2CapDacFlue[i, j]) for j in 0:hSlice-1)
    co2dacair = sum(value(xCo2CapDacAir[i, j]) for j in 0:hSlice-1)
    co2dacventflue = value(xCo2DacVentFlue[i])
    co2vent = value(xCo2Vent[i])
    push!(df_co,(
        value(co2fuel),
        value(co2pcc),
        value(co2pccout),
        value(co2pccvent),
        value(co2dacfluein),
        value(co2dacflue),
        value(co2dacair),
        value(co2dacventflue),
        value(co2vent),
        ))
end

# Power Data Frame.
df_pow = DataFrame(
                  Symbol("PowGasTur") => Float64[], # Pairs.
                  Symbol("PowSteaTurb") => Float64[],
                  Symbol("PowHp") => Float64[],
                  Symbol("PowIp") => Float64[],
                  Symbol("PowLp") => Float64[],
                  Symbol("PowUsePcc") => Float64[],
                  Symbol("PowUseDacFlue") => Float64[],
                  Symbol("PowUseDacAir") => Float64[],
                  Symbol("PowUseComp") => Float64[],
                  Symbol("AuxPowGasT") => Float64[],
                  Symbol("AuxPowSteaT") => Float64[],
                  Symbol("PowGross") => Float64[],
                  Symbol("PowOut") => Float64[],
                  Symbol("xeLoad") => Float64[],
                 )

# Pow / hSlice
for i in 0:tHorz
    ygastelecload = value(xeLoad[i])
    powgastur = value(xPowGasTur[i])
    powsteatur = value(xPowSteaTur[i])
    powhp = value(xPowHp[i])
    powip = value(xPowIp[i])
    powlp = value(xPowLp[i])
    powusepcc = value(xPowUsePcc[i])
    powusedacflue = sum(value(xPowUseDacFlue[i, j]) for j in 0:hSlice-1)
    powusedacair = sum(value(xPowUseDacAir[i, j]) for j in 0:hSlice-1)
    powusecomp = value(xPowUseComp[i])
    auxpowgast = value(xAuxPowGasT[i])
    auxpowsteat = value(xAuxPowSteaT[i])
    powgross = value(xPowGross[i])
    powout = value(xPowOut[i])
    push!(df_pow, (
        powgastur,
        powsteatur,
        powhp,
        powip,
        powlp,
        powusepcc,
        powusedacflue,
        powusedacair,
        powusecomp,
        auxpowgast,
        auxpowsteat,
        powgross,
        powout,
        ygastelecload
        ))
end

# Steam DataFrame
df_steam = DataFrame(
                     Symbol("CcRebDuty") => Float64[],
                     Symbol("SteaUsePcc") => Float64[],
                     Symbol("PccSteaSlack") => Float64[],
                     Symbol("DacSteaDuty") => Float64[],
                     Symbol("SteaUseDacFlue") => Float64[],
                     Symbol("SteaUseDacAir") => Float64[],
                     Symbol("DacSteaSlack") => Float64[],
                     Symbol("SideStea") => Float64[],
                     Symbol("DacSteaBaseDuty") => Float64[],
                     Symbol("SideSteaDac") => Float64[],
                     Symbol("Fuel") => Float64[]
                    )

# Steam / hSlice
for i in 0:tHorz
    ccrebduty = value(xCcRebDuty[i])
    steausepcc = value(xSteaUsePcc[i])
    pccsteaslack = value(xPccSteaSlack[i])
    dacsteaduty = value(xDacSteaDuty[i])
    steausedacflue = sum(value(xSteaUseDacFlue[i, j]) for j in 0:hSlice-1)
    steausedacair = sum(value(xSteaUseDacAir[i, j]) for j in 0:hSlice-1)
    dacsteaslack = value(xDacSteaSlack[i])
    sidesteam = value(xAllocSteam[i])
    dacsteabaseduty = value(xDacSteaBaseDuty[i])
    sidesteadac = value(xSideSteaDac[i])
    xfuel = value(xFuel[i])
    push!(df_steam,
        (
            ccrebduty,
            steausepcc,
            pccsteaslack,
            dacsteaduty,
            steausedacflue,
            steausedacair,
            dacsteaslack,
            sidesteam,
            dacsteabaseduty,
            sidesteadac,
            xfuel
            ))
end

# DAC-flue DataFrame
df_dac_flue = DataFrame(
    :time => Float64[],
    :xFflue => Float64[],
    :xSflue => Float64[],
    :xA0Flue => Float64[],
    :xA1Flue => Float64[],
    :xR0Flue => Float64[],
    :xR1Flue => Float64[]
    )

sliceFact = 1/hSlice
for i in 0:tHorz
    for j in 0:hSlice-1
        currtime = i + j * sliceFact
        push!(df_dac_flue,(
            currtime,
            value(xFflue[i, j]),
            value(xSflue[i, j]),
            value(xA0Flue[i, j]),
            value(xA1Flue[i, j]),
            value(xR0Flue[i, j]),
            value(xR1Flue[i, j]))
        )
    end
end

# DAC-air DataFrame
df_dac_air = DataFrame(
    :time => Float64[],
    :xFair => Float64[],
    :xSair => Float64[],
    :xA0Air => Float64[],
    :xA1Air => Float64[],
    :xA2Air => Float64[],
    :xR0Air => Float64[],
    :xR1Air => Float64[])

for i in 0:tHorz
    for j in 0:hSlice-1
        currtime = i + j * sliceFact
        push!(df_dac_air,(
            currtime,
            value(xFair[i, j]),
            value(xSair[i, j]),
            value(xA0Air[i, j]),
            value(xA1Air[i, j]),
            value(xA2Air[i, j]),
            value(xR0Air[i, j]),
            value(xR1Air[i, j])))
    end
end

df_pow_price = DataFrame(
                         price = Float64[]
                        )
for i in 0:tHorz
    push!(df_pow_price, (pow_price[i + 1],))
end

# Cost DataFrame
df_cost = DataFrame(
                    cNG = Float64[],
                    cCo2Em = Float64[],
                    cTransp = Float64[],
                    PowSales = Float64[]
                   )

for i in 0:tHorz
    cng = cNgPerMmbtu * value(xFuel[i])
    cco = cEmissionPrice * value(xCo2Vent[i])
    ctr = cCo2TranspPrice * value(xCo2Comp[i])
    cpow = pow_price[i + 1] * value(xPowOut[i])
    push!(df_cost, (cng, cco, ctr, cpow))
end

df_time_slice = DataFrame(:time => Float64[])

for i in 0:tHorz
    for j in 0:hSlice-1
        currtime = i + j * sliceFact
        push!(df_time_slice, (currtime, ))
    end
end

df_binary = DataFrame(
    :y0 => Float64[],
    :y1 => Float64[],
    :y2 => Float64[],
    :z00 => Float64[],
    :z01 => Float64[],
    :z02 => Float64[],
    :z10 => Float64[],
    :z11 => Float64[],
    :z12 => Float64[],
    :z20 => Float64[],
    :z21 => Float64[],
    :z22 => Float64[],
    )

for i in 0:tHorz
    push!(df_binary,
        (
            value(y[i, 0]),
            value(y[i, 1]),
            value(y[i, 2]),
            value(z[i, 0, 0]),
            value(z[i, 0, 1]),
            value(z[i, 0, 2]),
            value(z[i, 1, 0]),
            value(z[i, 1, 1]),
            value(z[i, 1, 2]),
            value(z[i, 2, 0]),
            value(z[i, 2, 1]),
            value(z[i, 2, 2]),
        ))
end

df_pr = DataFrame(
    :load0 => Float64[],
    :load1 => Float64[],
    :load2 => Float64[],
    :lambda0 => Float64[],
    :lambda10 => Float64[],
    :lambda11 => Float64[],
    :lambda20 => Float64[],
    :lambda21 => Float64[]
    )

for i in 0:tHorz
    push!(df_pr, (
        value(xLoad[i, 0]),
        value(xLoad[i, 1]),
        value(xLoad[i, 2]),
        value(lambda[i, 0, 0]),
        value(lambda[i, 1, 0]),
        value(lambda[i, 1, 1]),
        value(lambda[i, 2, 0]),
        value(lambda[i, 2, 1]),
        ))
end

# Write CSV
CSV.write("df_co.csv", df_co)
CSV.write("df_pow.csv", df_pow)
CSV.write("df_steam.csv", df_steam)
CSV.write("df_dac_flue.csv", df_dac_flue)
CSV.write("df_dac_air.csv", df_dac_air)
CSV.write("df_pow_price.csv", df_pow_price)
CSV.write("df_cost.csv", df_cost)
CSV.write("df_time_slice.csv", df_time_slice)
CSV.write("df_binary.csv", df_binary)
CSV.write("df_pr.csv", df_pr)

df_diss = DataFrame(
    "load0" => Float64[],
    "load1" => Float64[],
    "load2" => Float64[],
    "powgt0" => Float64[],
    "powgt1" => Float64[],
    "powgt2" => Float64[],
    "fuel0" => Float64[],
    "fuel1" => Float64[],
    "fuel2" => Float64[],
    "emiss0" => Float64[],
    "emiss1" => Float64[],
    "emiss2" => Float64[],
    "auxgt0" => Float64[],
    "auxgt1" => Float64[],
    "auxgt2" => Float64[],
    "powhp0" => Float64[],
    "powhp1" => Float64[],
    "powhp2" => Float64[],
    "powip0" => Float64[],
    "powip1" => Float64[],
    "powip2" => Float64[],
    "pccreb0" => Float64[],
    "pccreb1" => Float64[],
    "pccreb2" => Float64[],
    "dacsb0" => Float64[],
    "dacsb1" => Float64[],
    "dacsb2" => Float64[],
    "allocst0" => Float64[],
    "allocst1" => Float64[],
    "allocst2" => Float64[],
    "auxst0" => Float64[],
    "auxst1" => Float64[],
    "auxst2" => Float64[])

dvars = ["xLoad", "xPowGt", "xFuelNu", "xEmis", "xAuxGt", "xPowHpNu", "xPowIpNu", "xPccReb", "xDacSb", "xAllocS", "xAuxSt"]
for i in 0:tHorz
    l0 = []
    for v in dvars
        s = Symbol(v)
        l = [value(m[s][i, mod]) for mod in 0:nMod]
        l0 = vcat(l0, l)
    end
    push!(df_diss, l0)
end

CSV.write("df_diss.csv", df_diss)

max_stars_repo_name: dthierry/princetonDacLti

max_issues_repo_name: dthierry/princetonDacLti

max_forks_repo_name: dthierry/princetonDacLti


text:
function getlastdeath(barcode::PersistenceBarcode{T}, dims) where T
    lastdeath = zero(T)
    for i in dims
        isempty(barcode[i]) && continue
        lastdeath = max(lastdeath,
                        maximum(map(barcode[i]) do p
                                d = death(p)
                                isfinite(d) ? d : typemin(T)
                                end))
    end
    lastdeath
end

function checkdims(barcode::PersistenceBarcode, dims)
    dims == nothing           && return 0:dim(barcode)
    any(dims .< 0)            && throw(ArgumentError("Invalid `dims`!"))
    any(dims .> dim(barcode)) && throw(ArgumentError("Invalid `dims`!"))
    dims
end

@recipe function f(barcode::PersistenceBarcode{T};
                   dims = nothing, infinity = nothing) where T
    #TODO: Clean this part up?
    _dims = checkdims(barcode, dims)
    lastdeath = getlastdeath(barcode, _dims)
    if infinity == nothing
        infty = round(lastdeath, RoundUp) +
            (lastdeath > 1 ? length(digits(round(Int, lastdeath))) : 0)
    else
        infty = infinity
    end

    padding = lastdeath * 0.1
    xlim --> [0, infty + padding]
    h = 1
    for dim in _dims
        @series begin
            seriestype := :path
            label := "dim = $dim"
            linewidth --> 3
            marker --> 1

            xs = T[]; ys = T[]
            for (bth, dth) in barcode[dim]
                dth = isfinite(dth) ? dth : infty
                append!(xs, [bth, dth, NaN])
                append!(ys, [h,   h,   NaN])
                h += 1
            end
            xs, ys
        end
    end
    ylim --> [0, h+1]

    @series begin
        seriestype := :path
        label := "infinity"
        color := :grey
        [infty, infty], [0, h+1]
    end
end

@userplot PersistenceDiagram
@recipe function f(pd::PersistenceDiagram; dims=nothing, infinity=nothing)
    if length(pd.args) ≠ 1 || !(typeof(first(pd.args)) <: PersistenceBarcode)
        throw(ArgumentError("barcode is expecting a single PersistenceBarcode" *
                            " argument. Got: $(typeof(pd.args))"))
    end
    #TODO: Clean this part up?
    barcode = pd.args[1]
    _dims = checkdims(barcode, dims)
    lastdeath = getlastdeath(barcode, _dims)
    if infinity == nothing
        infty = round(lastdeath, RoundUp) +
            (lastdeath > 1 ? length(digits(round(Int, lastdeath))) : 0)
    else
        infty = infinity
    end
    padding = lastdeath * 0.1

    xlims --> (-padding, lastdeath)
    ylims --> (-padding, infty + padding)

    xlabel := "birth"
    ylabel := "death"

    # Births and deaths.
    for dim in _dims
        @series begin
            seriestype := :scatter
            label := "dim = $dim"

            xs = map(birth, barcode[dim])
            ys = map(death, barcode[dim])
            map!(y -> isfinite(y) ? y : infty, ys, ys)
            xs, ys
        end
    end
    # x = y line
    @series begin
        seriestype := :path
        label := ""
        color := :black
        [-padding, infty], [-padding, infty]
    end
    # infinity
    @series begin
        seriestype := :path
        label := "infinity"
        color := :grey
        [-padding, infty], [infty, infty]
    end
end


text:
"""
Transformers visit each node of the tree, and run the appropriate
method on it according to the node's data.

The methods (provided by the user) are dispatched according to the
type and the node's grammar rule. The returned value replaces the old
one in the structure.

They work bottom-up (or depth-first), starting with the leaves and
ending at the root of the tree.  Transformers can be used to
implement map - reduce patterns. Because nodes are reduced from
leaf to root, at any point the callbacks may assume the children
have already been transformed (if applicable).

`Transformer` can do anything `Visitor` can do, but because it
reconstructs the tree, it is slightly less efficient.

NOTE: A transformer without methods essentially performs a non-memoized deepcopy.
"""
abstract type Transformer end

"""
    visit_tokens(t::Transformer)

Should the transformer visit tokens in addition to rules.  Setting
this to false is an order of magnitude faster. For consistency with
Lark Defaults to `true`.  (For processing ignored tokens, use the
``lexer_callbacks`` options)
"""
visit_tokens(t::Transformer) = true

#== Python decorator

The _apply_decorator Transformer classmethod finds all
members of the subclass and submits them to the supplied
decorator function (with arguments as necessary), replacing
the previous value.

The main decorator is _visitor_args_func_dec, which
takes the function/method to be decorated together with
options.  This uses the utility smart_decorator function,
which calls the supplied "stupid" decorator in the
correct way depending on whether the argument is a
function, method, builtin function, partial application
or callable, by using the "with_self" argument to
the internal create_decorator function.

_visitor_args_func_dec seems to simply remove "self"
from the argument list if the with_self argument is
supplied, and adds the inline/meta/whole_tree arguments
to the function attributes. It appears that only
"inline" is used.

So: if a Transformer subclass is decorated with @inline_args,
every method not beginning with "_" and belonging to the
subclass has attribute 'inline' added to it as a True/False
value. Then, when creating a callback during parse tree construction,
if inline is 'False', then the user-written function is
splatted over the children. If 'inline' is True, then the
method is assumed to be written assuming that a list of children
is one of the arguments and no splatting is done.

As an example, EBNF_TO_BNF has inline declared as True. This
means that any of the internal methods are written assuming a
list of children is supplied.

We can duplicate this with macros, if we want.
==#

"""
    Transformer_InPlace

Non-recursive `Transformer`. Changes the tree in-place instead of returning new instances
"""
abstract type Transformer_InPlace <: Transformer end

"""
    Transformer_InPlaceRecursive

Recursive. Changes the tree in-place instead of returning new instances.
"""
abstract type Transformer_InPlaceRecursive <: Transformer end

"""
    transformer_func(t::Transformer,::Val{production},meta,args)

`transformer_func` transforms `production` with children `args`.
Methods of this function are defined using the @rule and @inline_rule
macros and chosen by Julia method dispatch. `meta` is only used by the
default `transformer_func`.
"""
function transformer_func end

"""
    transformer_func(t::Transformer,::Val{N},meta,children) where N

Default `transformer_func`
"""
transformer_func(t::Transformer,::Val{N},meta::Meta,children) where N = begin
    Tree(String(N),children,meta)
end

transformer_func(::Nothing,::Val{N},::Meta,children) where N = Tree(String(N),children)

"""
A token transformer function, works like the tree transformer
"""
function token_func end

# default
token_func(t::Transformer,::Val{N},tok) where N = begin
    return tok
end

# Note that `new_children` replaces `tr.children` as the idea is that they
# have been processed if they are provided as an argument
_call_userfunc(t::Transformer,tr::Tree; new_children = nothing) = begin
    children = new_children
    if new_children == nothing
        children = tr.children
    end
    # println("Processing $(tr.data)")
    return transformer_func(t,Val{Symbol(tr.data)}(),tr._meta,children)
end

_call_userfunc_token(t::Transformer,token) = begin
    return token_func(t,Val{Symbol(token)}(),token)
end

# No need for a channel here as we just collect the results anyway.
# We annotate the return type as otherwise Julia will assume the type
# is the type of the first returned value.
_transform_children(t::Transformer,children)::Array{Any} = begin
    map(children) do c
        try
            if c isa Tree
                _transform_tree(t,c)
            elseif visit_tokens(t) && c isa Token
                _call_userfunc_token(t,c)
            else
                c
            end
        catch e
            # For production, drop the error...
            # println("Warning: discarding $e")
            rethrow(e)
        end
    end
end

_transform_tree(t::Transformer,tree) = begin
    children = _transform_children(t,tree.children)
    #println("Children were $(tree.children), now $children")
    return _call_userfunc(t,tree,new_children=children)
end

"""
    transform(tr::Transformer,tree)

Transform parse tree `tree` according to rules defined for type `tr`.
Typically `tr` is a subtype of `Transformer`.
"""
transform(tr::Transformer,tree) = _transform_tree(tr,tree)

Base.:*(t1::Transformer,t2::Transformer) = TransformerChain([t1,t2])

__default__(tr::Transformer,data,children,meta) = begin
    #println("Warning: calling default transformer for $data and $tr")
    Tree(data,children,meta)
end

# Weird introspective decorator magic skipped

struct TransformerChain
    transformers::Array{Transformer}
end

TransformerChain(transf...) = TransformerChain(transf)

Base.:*(tc::TransformerChain,tc2::TransformerChain) = begin
    new_transformers = tc.transformers
    append!(new_transformers,tc2.transformers)   #TODO: Check we have a nice flat list
    TransformerChain(new_transformers)
end

Base.:*(tc::TransformerChain,t::Transformer) = begin
    new_transformers = tc.transformers
    push!(new_transformers,t)
    TransformerChain(new_transformers)
end

transform(tc::TransformerChain,tree) = begin
    for t in tc.transformers
        tree = transform(t,tree)
    end
    return tree
end

_transform_tree(tip::Transformer_InPlace,tree) = begin
    _call_userfunc(tip,tree)
end

_transform_tree(tipr::Transformer_InPlaceRecursive,tree) = begin
    tree.children = _transform_children(tipr,tree.children)
    return _call_userfunc(tipr,tree)
end

transform(tip::Transformer_InPlace, tree) = begin
    for subtree in tree
        subtree.children = _transform_children(tip,subtree.children)
    end
    return _transform_tree(tip,tree)
end

struct Transformer_NonRecursive <: Transformer end

visit_tokens(t::Transformer_NonRecursive) = false

transform(tnr::Transformer_NonRecursive, tree) = begin
    rev_postfix = []
    q = Union{Token,Tree}[tree]
    while !isempty(q)
        t = pop!(q)
        push!(rev_postfix,t)
        if t isa Tree
            append!(q,t.children)
        end
    end
    stack = []
    for x in reverse(rev_postfix)
        if x isa Tree
            size = length(x.children)
            if size > 0
                args = stack[end-size+1:end]
                stack = stack[1:end-size]
            else
                args = []
            end
            push!(stack,_call_userfunc(tnr,x,new_children=args))
        else
            push!(stack,x)
        end
    end
    t = stack[]
    return t
end

#== Visitors ==#

abstract type VisitorBase end

"""
    Visitor

User-defined methods called with an instance of a subtype of Visitor will visit
each node of the parse tree in top-down order without altering it or returning
a transformed tree.
"""
abstract type Visitor <: VisitorBase end

"""
    When called by `visit`, subtypes of `Visitor_Recursive` visit nodes of the parse tree
    in bottom-up order.
"""
abstract type Visitor_Recursive <: VisitorBase end

_call_userfunc(v::VisitorBase,tree) = transformer_func(v,Val{Symbol(tree.data)}(),Meta(),tree)

transformer_func(v::VisitorBase,::Val,tree) = tree

"""
    visit(v::Visitor,tree)

Visit each node of `tree`, calling methods defined for `v` on each node.
"""
visit(v::Visitor,tree) = begin
    for subtree in tree
        _call_userfunc(v,subtree)
    end
    return tree
end

visit_topdown(v::Visitor,tree) = begin
    for subtree in iter_subtrees_topdown(tree)
        _call_userfunc(v,subtree)
    end
    return tree
end

visit(v::Visitor_Recursive,tree) = begin
    for child in tree.children
        if child isa Tree
            visit(v,child)
        end
    end
    _call_userfunc(v,tree)
    return tree
end

visit_topdown(v::Visitor_Recursive,tree) = begin
    _call_userfunc(tree)
    for child in tree.children
        if child isa Tree
            visit_topdown(v,child)
        end
    end
    return tree
end

"""
    Subtypes of the `Interpreter` type do not automatically visit children
    of the parse tree node when called by `visit`.
"""
abstract type Interpreter end

"""
    transformer_func(i::Interpreter,::Val,meta,tree)

Default interpreter visitor does not automatically visit children
"""
transformer_func(i::Interpreter,::Val,meta,tree) = visit_children(i,tree)

visit(inter::Interpreter,tree) = begin
    transformer_func(inter,Val{Symbol(tree.data)}(),Meta(),tree)
end

visit_children(inter::Interpreter,tree) = begin
    map(tree.children) do child
        if typeof(child) == Tree
            visit(inter,child)
        else
            child
        end
    end
end

abstract type CollapseAmbiguities <: Transformer end

@rule _ambig(c::CollapseAmbiguities,options) = sum(options,[])

__default__(c::CollapseAmbiguities,data,children_lists,meta) = begin
    [Tree(data,children,meta) for children in combine_alternatives(children_lists)]
end

__default_token__(c::CollapseAmbiguities,t) = [t]


text:
include("Mortgage.jl")
mutable struct carMortgage <: Mortgage
    #fields
end

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia