text:
using ApproxFun, SpecialFunctions, LinearAlgebra, Test

@testset "Readme" begin
    @testset "Calculus and algebra" begin
        x = Fun(identity,0..10)
        f = sin(x^2)
        g = cos(x)

        @test ≈(f(.1),sin(.1^2);atol=1000eps())

        h = f + g^2
        r = roots(h)
        rp = roots(differentiate(h))

        @test norm(h.(r))<1000eps()
        @test norm(h'.(rp))<100000eps()
    end

    @testset "Differentiation and Integration" begin
        f = Fun(x->exp(x),-1..1)
        @test norm(f-f')<1000eps()

        g = cumsum(f)
        g = g + f(-1)
        @test norm(f-g)<100eps()

        x = Fun(identity)
        f = exp(x)
        g = f/sqrt(1-x^2)

        space(f),domain(f)
        space(g),domain(g)

        f = Fun(x->cospi(5x))
        g = abs(f)
        space(f)
        space(g)

        x = Fun()
        f = erf(x)
        g = besselj(3,exp(f))
        h = airyai(10asin(f)+2g)
    end

    @testset "ODE" begin
        a,b = -1000,200
        x = Fun(identity, a..b)
        d = domain(x)
        D = Derivative(d)
        B = Dirichlet(d)
        L = D^2 - x
        u = [B;L] \ [[airyai(a),airyai(b)],0]

        @test ≈(u(0.),airyai(0.);atol=10000eps())

        ##) Nonlinear BVPs
        x=Fun()
        u0=0.0x

        N=u->[u(-1.)-1.,u(1.)+0.5,0.001u''+6*(1-x^2)*u'+u^2-1.]
        u=newton(N,u0)

        @test norm(N(u)[end]) ≤ 1000eps()
    end

    @testset "Periodic" begin
        f = Fun(cos,Fourier(-π..π))
        @test norm(differentiate(f) + Fun(sin,Fourier(-π..π))) < 100eps()

        s = Chebyshev(-π..π)
        a = Fun(t-> 1+sin(cos(2t)),s)
        L = Derivative() + a
        f = Fun(t->exp(sin(10t)),s)
        B = periodic(s,0)
        uChebyshev = [B;L]\[0.,f]

        s = Fourier(-π..π)
        a = Fun(t-> 1+sin(cos(2t)),s)
        L = Derivative() + a
        f = Fun(t->exp(sin(10t)),s)
        uFourier = L\f

        @test uChebyshev(0.) ≈ uFourier(0.)
    end

    @testset "Sampling" begin
        ## Sampling
        f = abs(Fun(sin,-5..5))
        x = ApproxFun.sample(f,10)
    end

    @testset "PDE" begin
        d = ChebyshevInterval()^2                            # Defines a rectangle

        # @time u = \([Dirichlet(d);Laplacian(d)+100I],
        #                     [ones(∂(d));0.];tolerance=1E-10)      # First four entries of rhs are
        #

        QR = qr([Dirichlet(d);Laplacian()+100I])
                @time ApproxFun.resizedata!(QR,:,4000)
                @time u = \(QR,
                                [ones(∂(d));0.];tolerance=1E-7)

        @test u(0.1,1.) ≈ 1.0
        @test ≈(u(0.1,0.2),-0.02768276827514463;atol=1E-8)
    end

    @testset "BigFloat" begin
        setprecision(1000) do
            d=BigFloat(0)..BigFloat(1)
            D=Derivative(d)
            u=[ldirichlet();D-I]\[1;0]
            @test u(1) ≈ exp(BigFloat(1))
        end
    end
end


text:
#! /usr/bin/julia

# Rosetta Code, Playing cards

type DeckDesign{T<:Integer,U<:String}
    rlen::T
    slen::T
    ranks::Array{U,1}
    suits::Array{U,1}
    hlen::T
end

type Deck{T<:Integer}
    cards::Array{T,1}
    design::DeckDesign
end

Deck(n::Integer, des::DeckDesign) = Deck([n], des)

function pokerlayout()
    r = [map(string, 2:10), "J", "Q", "K", "A"]
    r = map(utf8, r)
    s = ["\u2663", "\u2666", "\u2665", "\u2660"]
    DeckDesign(13, 4, r, s, 5)
end

function fresh(des::DeckDesign)
    Deck(collect(1:des.rlen*des.slen), des)
end

Base.isempty(d::Deck) = isempty(d.cards)
Base.empty!(d::Deck) = empty!(d.cards)
Base.length(d::Deck) = length(d.cards)
Base.endof(d::Deck) = endof(d.cards)
Base.shuffle!(d::Deck) = shuffle!(d.cards)
Base.sort!(d::Deck) = sort!(d.cards)
Base.getindex(d::Deck, r) = Deck(getindex(d.cards, r), d.design)
Base.size(d::Deck) = (d.design.rlen, d.design.slen)
function Base.print(d::Deck)
    sz = size(d)
    r = map(x->d.design.ranks[ind2sub(sz, x)[1]], d.cards)
    s = map(x->d.design.suits[ind2sub(sz, x)[2]], d.cards)
    join(r.*s, " ")
end

function deal!{T<:Integer}(d::Deck, hlen::T)
    if hlen < length(d)
        hand = Deck(d.cards[1:hlen], d.design)
        d.cards = d.cards[hlen+1:end]
    else
        hand = d
        empty!(d)
    end
    return hand
end

function deal!(d::Deck)
    deal!(d, d.design.hlen)
end

function pretty(d::Deck)
    s = ""
    llen = d.design.rlen
    dlen = length(d)
    for i in 1:llen:dlen
        j = min(i+llen-1, dlen)
        s *= print(d[i:j])*"\n"
    end
    chop(s)
end

d = fresh(pokerlayout())
println("A new poker deck:")
println(pretty(d))

shuffle!(d)
println()
println("The deck shuffled:")
println(pretty(d))

n = 4
println()
println("Deal ", n, " hands:")
for i in 1:n
    h = deal!(d)
    println(pretty(h))
end

println()
println("And now the deck contains:")
println(pretty(d))

ext: jl

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
#
using DiffMesh
using ModelingToolkit, NeuralPDE
using DiffEqFlux, LinearAlgebra, Plots
using ParameterSchedulers:Scheduler,Sin,SinDecay2,SinExp
using GalacticOptim
#----------------------------------------------------------------------------#
nprt = 2
sdf = (x,y) -> @. (x^2 + y^2) - 1

nprt = 6
pou  = nnPOU2D(nprt,sdf;depth=4,act=relu)

nx = ny = 4
bsp = TensorProdCosines(nprt,nx,ny)
net = MeshNet2D(pou,bsp)

p,re = Flux.destructure(net)
#--------------------------------------#
function utrue(x,y)

    r2 = @. x*x + y*y

    u .= @. 1 - r2
    u .= @. sin(pi*x)
    u .= @. x*y^3 + sin(r2)  + exp(x^2*y)
    return u
end

it = 0 # iter
#----------------------------------------------------------------------------#
function makeplot(phi)

    fpp = phi(xpp,ypp)
    fp  = reshape(fpp,np,np)
    ep  = up .- fp

    p1 = heatmap(Ωp.*fp)
    p2 = heatmap(Ωp.*up)
    p3 = heatmap(Ωp.*ep)

    p1 = plot!(p1,title="Pred")
    p2 = plot!(p2,title="True")
    p3 = plot!(p3,title="Err ")

    plt = plot(p1,p2,p3)

    display(plt)
    return
end
#--------------------------------------#
function callback(p,l)
    phi = re(p)

    leastSq!(phi,xe,ye,ue)

    fpp = phi(xpp,ypp)
    epp = (upp .- fpp) .* Ωpp
    er  = norm(epp,Inf)

    global it += 1
    println("it: $it, loss: $l, pointwise er: $er")
    (it % 10) == 0 && makeplot(phi)
    (l < 1e-6) && return true
    return false
end
#----------------------------------------------------------------------------#
# NeuralPDE setup

@parameters x y
@variables u(..)
Dx  = Differential(x)
Dy  = Differential(y)
Dxx = Dx ∘ Dx
Dyy = Dy ∘ Dy

eq = [Dxx(u(x,y)) + Dyy(u(x,y))  ~ 1f0]
bc = [u(-1f0) ~ 1f0
     ,u( 1f0) ~ 1f0]

domain = [x ∈ ModelingToolkit.Interval(-1f0,1f0),
          y ∈ ModelingToolkit.Interval(-1f0,1f0)]

strategy = GridTraining(0.1)
discretization = PhysicsInformedNN(net,strategy)
pdesys = PDESystem(eq,bc,domain,[x,y],[u],name=:PoissonDisk)
prob = discretize(pdesys,discretization)
#--------------------------------------#
res = GalacticOptim.solve(prob,ADAM(1e-2),cb=callback,maxiters=200)
#----------------------------------------------------------------------------#
return


text:
using SpecialFunctions: lgamma
import SpecialFunctions: lbeta, beta

# Call back functions for `lbeta` and `beta`
# Those in `SpecialFunctions.jl` are implemented for `(::Number, ::Number)` only,
# for which `KnetArray` cannot broadcast through.
# NOTE: `_lbeta` and `_beta` are implemented for explictly test purpose.
# TODO: implement this broadcasting in Knet.jl
_lbeta(α, β) = lgamma.(α) + lgamma.(β) - lgamma.(α + β)
_beta(α, β) = exp.(_lbeta(α, β))
lbeta(α, β) = _lbeta(α, β)
beta(α, β) = _beta(α, β)


text:
module Aligator

export aligator, extract_loop, closed_forms, invariants

using PyCall
using SymPy
using Singular
using AlgebraicDependencies

const AppliedUndef = PyCall.PyNULL()

include("utils.jl")
include("closedform.jl")
include("recurrence.jl")
include("loop.jl")
include("invariants.jl")
include("parse_julia.jl")
include("ideals.jl")
include("singular_imap.jl")

function aligator(str::String)
    _, total = @timed begin

        loop, time = @timed extract_loop(str)
        @debug "Recurrence extraction" time

        cforms, time = @timed closed_forms(loop)
        @debug "Recurrence solving" time

        invs, time = @timed invariants(cforms)
        @debug "Ideal computation" time
    end
    @debug "Total time needed" total

    return InvariantIdeal(invs)
end

function __init__()
    copy!(AppliedUndef, PyCall.pyimport_conda("sympy.core.function", "sympy").AppliedUndef)

    include(joinpath(@__DIR__,"..", "benchmark", "singlepath.jl"))
    include(joinpath(@__DIR__,"..", "benchmark", "multipath.jl"))

    singlepath = [:cohencu, :freire1, :freire2, :(petter(1)), :(petter(2)), :(petter(3)), :(petter(4))]
    multipath = [:divbin, :euclidex, :fermat, :knuth, :lcm, :mannadiv, :wensley]
end

struct InvariantIdeal
    ideal::sideal
end

function Base.show(io::IO, I::InvariantIdeal)
    println(io, "Invariant ideal with $(ngens(I.ideal))-element basis:")
    Base.print_array(io, gens(I.ideal))
end

end # module


text:
_resolve_leaf(tab, s) = SelectionResult(_resolve(tab, s.s), s.r, s.t)

function _resolve(tab, s::OrSelection)
    res = union_results(tab, _resolve_leaf(tab, first(s)), _resolve_leaf(tab, last(s)))
    if bool(s)
        res
    else
        SelectionResult(positive_selection(tab, colnames(res), bool(s)), nothing, nothing)
    end
end

function _resolve(tab, s::AndSelection)
    res = intersect_results(tab, _resolve_leaf(tab, first(s)), _resolve_leaf(tab, last(s)))
    if bool(s)
        res
    else
        SelectionResult(positive_selection(tab, colnames(res), bool(s)), nothing, nothing)
    end
end

function _resolve(tab, s::SubSelection)
    res = sub_results(tab, _resolve_leaf(tab, first(s)), _resolve_leaf(tab, last(s)))
    if bool(s)
        res
    else
        SelectionResult(positive_selection(tab, colnames(res), bool(s)), nothing, nothing)
    end
end

# Union ############################################################################################

function _union_results_inner!(tab, x, y, out)
    length(x) == 0 && return y
    length(y) == 0 && return x
    yset = Set(colnames(y)) # TODO: fix Any
    # yterms = (;(colnames(y) .=> keyfuncs(y) .=> valfuncs(y))...)  # other_cols() cannot be a nt key
    yterms = Dict(colnames(y) .=> keyfuncs(y) .=> valfuncs(y))
    for xterm in x
        xname = colname(xterm)
        if xname in yset
            funs = yterms[xname]
            push!(
                out,
                SelectionTerm(
                    xname,
                    extend(keyfunc(xterm), first(funs)),
                    extend(valfunc(xterm), last(funs))
                )
            )
            delete!(yset, xname)
        else
            push!(out, xterm)
        end
    end
    for yterm in filter(term -> colname(term) in yset, y)
        push!(out, yterm)
    end
    SelectionResult(map(identity, out))
end

function union_results(tab, x::SelectionResult{<:SelectionTerm{ElseSelection}}, y::SelectionResult)
    SelectionResult(
        vcat(
            SelectionResult(
                setdiff(colnames(tab), colnames(y)),
                keyfunc(x[1]),
                valfunc(x[1])
            ),
            y
        )
    )
end

function union_results(tab, x::SelectionResult, y::SelectionResult{<:SelectionTerm{ElseSelection}})
    SelectionResult(
        vcat(
            x,
            SelectionResult(
                setdiff(colnames(tab), colnames(x)),
                keyfunc(y[1]),
                valfunc(y[1])
            )
        )
    )
end

function union_results(
        tab,
        x,
        y
    )
    out = []  # TODO: fix Any
    _union_results_inner!(tab, x, y, out)
end

# Intersection #####################################################################################

function intersect_results(tab, x::SelectionResult, y::SelectionResult)
    length(x) == 0 && return x
    length(y) == 0 && return y

    out = SelectionTerm{Symbol}[]
    yset = Set(colnames(y))
    yterms = (;(colnames(y) .=> keyfuncs(y) .=> valfuncs(y))...)
    for xterm in x
        xname = colname(xterm)
        if xname in yset
            funs = yterms[xname]
            push!(
                out,
                SelectionTerm(
                    xname,
                    extend(keyfunc(xterm), first(funs)),
                    extend(valfunc(xterm), last(funs))
                )
            )
        end
    end
    SelectionResult(out)
end

function intersect_results(tab, x::SelectionResult{<:SelectionTerm{T}}, y::SelectionResult) where {
        T<:Union{ElseSelection,OtherSelection}
    }
    ArgumentError("Cannot intersect with `$(T)`.")
end

function intersect_results(tab, x::SelectionResult, y::SelectionResult{<:SelectionTerm{T}}) where {
        T<:Union{ElseSelection,OtherSelection}
    }
    ArgumentError("Cannot intersect with `$(T)`.")
end

function intersect_results(tab, x::S, y::S) where {
        T<:Union{ElseSelection,OtherSelection},
        S<:SelectionResult{<:SelectionTerm{T}}
    }
    throw(ArgumentError("Cannot resolve two `$(T)`s."))
end

# Setdiff ##########################################################################################

function sub_results(tab, x::SelectionResult, y::SelectionResult)
    yset = Set(colnames(y))
    SelectionResult(filter(term->!(colname(term) in yset), x))
end

function sub_results(tab, x::SelectionResult{<:SelectionTerm{T}}, y::SelectionResult) where {
        T<:Union{ElseSelection,OtherSelection}
    }
    ArgumentError("Cannot subtract `$(T)`.")
end

function sub_results(tab, x::SelectionResult, y::SelectionResult{<:SelectionTerm{T}}) where {
        T<:Union{ElseSelection,OtherSelection}
    }
    ArgumentError("Cannot subtract `$(T)`.")
end

function sub_results(tab, x::S, y::S) where {
        T<:Union{ElseSelection,OtherSelection},
        S<:SelectionResult{<:SelectionTerm{T}}
    }
    throw(ArgumentError("Cannot resolve two `$(T)`s."))
end


text:
# # 2d mixed mesh import
# In this example, we will import a 2d mesh of a Compact Tension (CT) specimen
# from Abaqus. The mesh has both triangular and quadrilateral cells, and we define
# a few sets in Abaqus.
#
#-
# ## Abaqus setup
#
# ![](compact_tension_specimen.svg)
# Figure 1: Geometry and sets, mesh overview and detailed mesh from Abaqus.
#
# We have created sets for the hole (blue line: "Hole"), the symmetry edge
# (red line: "Symmetry"), and the crack zone (green area: "CrackZone").
#
# ## Importing mesh
# The mesh above was created in Abaqus cae, and in this example we
# import the generated input file: [compact_tension.inp](compact_tension.inp)
#
using Ferrite, FerriteMeshParser

grid = get_ferrite_grid("compact_tension.inp")

# We can now inspect this grid, showing that we have different cell types
println(typeof(grid))
println(unique(typeof.(getcells(grid))))    # The different cell types in the grid

# Furthermore, the node and cell sets are imported
println([(key, length(set)) for (key, set) in getnodesets(grid)])
println([(key, length(set)) for (key, set) in getcellsets(grid)])

# As we see, in addition to the sets created in Abaqus, the cellsets also include a set
# for each abaqus element type (useful if you for example defined reduced integration
# in only part of the domain and want to have this in Ferrite). Finally, facesets are
# automatically created by default (can be turned off by `generate_facesets=false`
# argument) based on the nodesets:
println([(key, length(set)) for (key, set) in getfacesets(grid)])
# Clearly, the faceset `"CrackZone"` doesn't make much sense, but unless the mesh is
# very large it doesn't hurt. The facesets can be created manually from each nodeset
# by using the `create_faceset` function:
faceset = create_faceset(grid, getnodeset(grid,"Hole"));
# This can if desired be merged into the grid by
merge!(getfacesets(grid), Dict("HoleManual" => faceset))
println([(key, length(set)) for (key, set) in getfacesets(grid)])

#
#md # ## [Plain Program](@id compact-tension-plain-program)
#md #
#md # Below follows a version of the program without any comments.
#md # The file is also available here: [compact_tension.jl](compact_tension.jl)
#md #
#md # ```julia
#md # @__CODE__
#md # ```


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
# We want two dumb turbines.

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 * 30 * 1 - 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
nUnit = 1  # We go from 0 to 1

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

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

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

# "Actual variables power, fuel, etc."
@variable(m, 0 <= xActualLoad[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])

# You still have a single steam power generation,
# so what is the point of
# going below 30?
@variable(m, 0 <= xPowSteaTur[0:tHorz])
@variable(m, 0 <= xAuxPowSteaT[0:tHorz])

# Unit load

# Disagretated variables
@variable(m, xLoadD[0:tHorz, 0:nUnit, 0:nMod])
@variable(m, 0 <= xPowGtD[0:tHorz, 0:nUnit, 0:nMod])
@variable(m, 0 <= xFuelD[0:tHorz, 0:nUnit, 0:nMod])
@variable(m, 0 <= xEmisD[0:tHorz, 0:nUnit, 0:nMod])
@variable(m, 0 <= xAuxGtD[0:tHorz, 0:nUnit, 0:nMod])
@variable(m, 0 <= xPowHpD[0:tHorz, 0:nUnit, 0:nMod])
@variable(m, 0 <= xPowIpD[0:tHorz, 0:nUnit, 0:nMod])
@variable(m, 0 <= xPccRebD[0:tHorz, 0:nUnit, 0:nMod])
@variable(m, 0 <= xDacSbD[0:tHorz, 0:nUnit, 0:nMod])
@variable(m, 0 <= xAllocD[0:tHorz, 0:nUnit, 0:nMod])
@variable(m, 0 <= xAuxStD[0:tHorz, 0:nUnit, 0:nMod])

# Per unit variables
@variable(m, 0 <= xLoadU[0:tHorz, 0:nUnit] <= 100)
@variable(m, 0 <= xPowGtU[0:tHorz, 0:nUnit])
@variable(m, 0 <= xFuelU[0:tHorz, 0:nUnit])
@variable(m, 0 <= xEmisU[0:tHorz, 0:nUnit])
@variable(m, 0 <= xAuxGtU[0:tHorz, 0:nUnit])
@variable(m, 0 <= xPowHpU[0:tHorz, 0:nUnit])
@variable(m, 0 <= xPowIpU[0:tHorz, 0:nUnit])
@variable(m, 0 <= xPccRebU[0:tHorz, 0:nUnit])
@variable(m, 0 <= xDacSbU[0:tHorz, 0:nUnit])
@variable(m, 0 <= xAllocU[0:tHorz, 0:nUnit])
@variable(m, 0 <= xAuxStU[0:tHorz, 0:nUnit])

# 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])
@variable(m, 0 <= xCo2DacFlueIn[0:tHorz])
@variable(m, 0 <= xPccSteaSlack[0:tHorz])

# Dac-Flue, takes one time slot for adsorption
@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, takes two time slots for adsorption
@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 combination
@constraint(m, cConvxEq[i = 0:tHorz, j = 0:nUnit, mod = 0:nMod],
    sum(lambda[i, j, mod, k] * ep_list[mod + 1][k + 1] for k in extrPoint[mod])
    == xLoadD[i, j, mod]  # There is only a single extreme
    )

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

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

# Overall Disjunction

@constraint(m, gasLeq[i = 0:tHorz, j = 0:nUnit],
    xLoadU[i, j] == sum(xLoadD[i, j, mod] for mod in 0:nMod)
    )

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

# Switch
@constraint(m, switchConEq[i = 1:tHorz, j = 0:nUnit, mod = 0:nMod],
    sum(z[i, j, k, mod] for k in 0:nMod) -
    sum(z[i, j, mod, k] for k in 0:nMod)
    == y[i, j, mod] - y[i - 1, j, mod])

# Forbidden
# from off to full
@constraint(m, forbConOffFullEq[i = 0:tHorz, j = 0:nUnit],
    z[i, j, 0, 2] == 0)
# from full to warm-up
@constraint(m, forbConEq2[i = 0:tHorz, j = 0:nUnit],
    z[i, j, 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]
# good
#@constraints(m, begin
#    rampup1eq[i=0:tHorz-1, j=0:nUnit], xLoadD[i+1, j, 1] - xLoadD[i, j, 1] <= rup[2] + maximum(ep_list[1 + 1]) * (y[i+1, j, 1] - y[i, j, 1])
#    rampup2eq[i=0:tHorz-1, j=0:nUnit], xLoadD[i+1, j, 2] - xLoadD[i, j, 2] <= rup[3] + maximum(ep_list[1 + 2]) * (y[i+1, j, 1] - y[i, j, 1])
#end)

@constraints(m, begin
    rampup1eq[i=0:tHorz-1, j=0:nUnit], xLoadD[i+1, j, 1] - xLoadD[i, j, 1] <= rup[2] + (minimum(ep_list[1 + 1])-1) * (y[i+1, j, 1] - y[i, j, 1])
    rampup2eq[i=0:tHorz-1, j=0:nUnit], xLoadD[i+1, j, 2] - xLoadD[i, j, 2] <= rup[3] + (minimum(ep_list[1 + 2])-1) * (y[i+1, j, 2] - y[i, j, 2])
end)

@constraints(m, begin
    rampdo1eq[i=0:tHorz-1, j=0:nUnit], xLoadD[i, j, 1] - xLoadD[i+1, j, 1] <= rdown[2] + (minimum(ep_list[1 + 1])-1) * (y[i, j, 1] - y[i+1, j, 1])
    rampdo2eq[i=0:tHorz-1, j=0:nUnit], xLoadD[i, j, 2] - xLoadD[i+1, j, 2] <= rdown[3] + (minimum(ep_list[1 + 2])-1) * (y[i, j, 2] - y[i+1, j, 2])
end)

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

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

# Disagregated variables
@constraint(m, powgtdeq[i = 0:tHorz, j = 0:nUnit, mod = 0:nMod],
    xPowGtD[i, j, mod] ==
    (aPowGt[mod] * xLoadD[i, j, mod] + bPowGt[mod] * y[i, j, mod]) * 0.5
    )
@constraint(m, fueldeq[i = 0:tHorz, j = 0:nUnit, mod = 0:nMod],
    xFuelD[i, j, mod] ==
    (aFuel[mod] * xLoadD[i, j, mod] + bFuel[mod] * y[i, j, mod]) * 0.5
    )
@constraint(m, emissdeq[i = 0:tHorz, j = 0:nUnit, mod = 0:nMod],
    xEmisD[i, j, mod] ==
    (aEmis[mod] * xLoadD[i, j, mod] + bEmis[mod] * y[i, j, mod]) * 0.5
    )
@constraint(m, auxgtdeq[i = 0:tHorz, j = 0:nUnit, mod = 0:nMod],
    xAuxGtD[i, j, mod] ==
    (aAuxGt[mod] * xLoadD[i, j, mod] + bAuxGt[mod] * y[i, j, mod]) * 0.5
    )
@constraint(m, powhpdeq[i = 0:tHorz, j = 0:nUnit, mod = 0:nMod],
    xPowHpD[i, j, mod] ==
    (aPowHp[mod] * xLoadD[i, j, mod] + bPowHp[mod] * y[i, j, mod]) * 0.5
    )
@constraint(m, powipdeq[i = 0:tHorz, j = 0:nUnit, mod = 0:nMod],
    xPowIpD[i, j, mod] ==
    (aPowIp[mod] * xLoadD[i, j, mod] + bPowIp[mod] * y[i, j, mod]) * 0.5
    )
@constraint(m, pccrebdeq[i = 0:tHorz, j = 0:nUnit, mod = 0:nMod],
    xPccRebD[i, j, mod] ==
    (aCcReb[mod] * xLoadD[i, j, mod] + bCcReb[mod] * y[i, j, mod]) * 0.5
    )
@constraint(m, dacsbtdeq[i = 0:tHorz, j = 0:nUnit, mod = 0:nMod],
    xDacSbD[i, j, mod] ==
    (aDacSb[mod] * xLoadD[i, j, mod] + bDacSb[mod] * y[i, j, mod]) * 0.5
    )
@constraint(m, allocstdeq[i = 0:tHorz, j = 0:nUnit, mod = 0:nMod],
    xAllocD[i, j, mod] ==
    (aSideS[mod] * xLoadD[i, j, mod] + bSideS[mod] * y[i, j, mod]) * 0.5
    )
@constraint(m, auxstdeq[i = 0:tHorz, j = 0:nUnit, mod = 0:nMod],
    xAuxStD[i, j, mod] ==
    (aAuxSt[mod] * xLoadD[i, j, mod] + bAuxSt[mod] * y[i, j, mod]) * 0.5
    )

# Unit variable constraints
@constraint(m, powgtueq[i = 0:tHorz, j = 0:nUnit],
    xPowGtU[i, j] == sum(xPowGtD[i, j, mod] for mod in 0:nMod)
    )
@constraint(m, fuelueq[i = 0:tHorz, j = 0:nUnit],
    xFuelU[i, j] == sum(xFuelD[i, j, mod] for mod in 0:nMod)
    )
@constraint(m, emissueq[i = 0:tHorz, j = 0:nUnit],
    xEmisU[i, j] == sum(xEmisD[i, j, mod] for mod in 0:nMod)
    )
@constraint(m, auxgtueq[i = 0:tHorz, j = 0:nUnit],
    xAuxGtU[i, j] == sum(xAuxGtD[i, j, mod] for mod in 0:nMod)
    )
@constraint(m, powhpueq[i = 0:tHorz, j = 0:nUnit],
    xPowHpU[i, j] == sum(xPowHpD[i, j, mod] for mod in 0:nMod)
    )
@constraint(m, powipueq[i = 0:tHorz, j = 0:nUnit],
    xPowIpU[i, j] == sum(xPowIpD[i, j, mod] for mod in 0:nMod)
    )
@constraint(m, pccrebueq[i = 0:tHorz, j = 0:nUnit],
    xPccRebU[i, j] == sum(xPccRebD[i, j, mod] for mod in 0:nMod)
    )
@constraint(m, dacsbtueq[i = 0:tHorz, j = 0:nUnit],
    xDacSbU[i, j] == sum(xDacSbD[i, j, mod] for mod in 0:nMod)
    )
@constraint(m, allocstueq[i = 0:tHorz, j = 0:nUnit],
    xAllocU[i, j] == sum(xAllocD[i, j, mod] for mod in 0:nMod)
    )
@constraint(m, auxstueq[i = 0:tHorz, j = 0:nUnit],
    xAuxStU[i, j] == sum(xAuxStD[i, j, mod] for mod in 0:nMod)
    )

# Half the load of each GT
@constraint(m ,eloadeq[i = 0:tHorz],
    xActualLoad[i] == sum(xLoadU[i, j] * 0.5 for j in 0:nUnit)
    )

@constraint(m, powGasTur[i = 0:tHorz],
            xPowGasTur[i] == sum(xPowGtU[i, j] for j in 0:nUnit)
           )
#
@constraint(m, fuelEq[i = 0:tHorz],
            xFuel[i] == sum(xFuelU[i, j] for j in 0:nUnit)
           )
#
@constraint(m, co2FuelEq[i = 0:tHorz],
            xCo2Fuel[i] == sum(xEmisU[i, j] for j in 0:nUnit)
           )

@constraint(m, auxPowGasT[i = 0:tHorz],
            xAuxPowGasT[i] == sum(xAuxGtU[i, j] for j in 0:nUnit)
           )
#
# Steam
#
@constraint(m, powHpEq[i = 0:tHorz],
            xPowHp[i] == sum(xPowHpU[i, j] for j in 0:nUnit)
           )
#
@constraint(m, powIpEq[i = 0:tHorz],
            xPowIp[i] == sum(xPowIpU[i, j] for j in 0:nUnit)
           )

#
@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(xPccRebU[i, j] for j in 0:nUnit)
           )

@constraint(m, dacSteaDutyEq[i = 0:tHorz],
            xDacSteaBaseDuty[i] == sum(xDacSbU[i, j] for j in 0:nUnit)
           )

@constraint(m, sideSteaEloadEq[i = 0:tHorz],
            xAllocSteam[i] == sum(xAllocU[i, j] for j in 0:nUnit)
           )

@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(xAuxStU[i, j] for j in 0:nUnit)
           )

# 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, 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, 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])

oneHourofGas60Load = cNgPerMmbtu * (bFuelEload + 60 * aFuelEload)

shutdowncost = Dict(0 => oneHourofGas60Load * 5, 1 => oneHourofGas60Load)
shutdownfromwarm = cNgPerMmbtu * (bFuelEload + 60 * aFuelEload) * 24
startupcost = Dict(0 => oneHourofGas60Load * 24/2, 1 => oneHourofGas60Load/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[j] * z[i, j, 2, 0]
        + startupcost[j] * z[i, j, 0, 1]
        + shutdownfromwarm * z[i, j, 1, 0]
        for i in 0:tHorz for j in 0:nUnit)
    )

@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("xActualLoad") => Float64[],
                 )

# Pow / hSlice
for i in 0:tHorz
    ygastelecload = value(xActualLoad[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

# Binary variables
# Create the DataFrame for the binary variables.
bvars_names = ["y", "z"]

l0 = []
for v in bvars_names
    for j in 0:nUnit
        for mod in 0:nMod
            if v == "z"
                l = [v * string(j) * string(mod) * string(m2) => Float64[]
                for m2 in 0:nMod]
            else
                l = [v * string(j) * string(mod) => Float64[]]
            end
            global l0 = vcat(l0, l)
        end
    end
end

df_binary = DataFrame(l0)
# Add the values to the DataFrame of binaries.
for i in 0:tHorz
    l0 = []
    for v in bvars_names
        s = Symbol(v)
        if v == "z"
            l = [value(m[s][i, j, m1, m2]) for j in 0:nUnit
             for m1 in 0:nMod  for m2 in 0:nMod]
        else
            l = [value(m[s][i, j, mod]) for j in 0:nUnit for mod in 0:nMod]
        end
        l0 = vcat(l0, l)
    end
    push!(df_binary, l0)
end

# ``Main'' disagregated variables.
dis_names = ["xLoadD", "lambda"]
l0 = []
for v in dis_names
    for j in 0:nUnit # Per unit
        for mod in 0:nMod # Per mod
            if v == "lambda"
                l = [v * string(j) * string(mod) * string(k) => Float64[]
                for k in extrPoint[mod]]
            else
                l = [v * string(j) * string(mod) => Float64[]]
            end
            global l0 = vcat(l0, l)
        end
    end
end

df_mdv = DataFrame(l0)

for i in 0:tHorz
    l0 = []
    for v in dis_names
        s = Symbol(v)
        if v == "lambda"
            l = []
            for j in 0:nUnit
                for mod in 0:nMod
                    l1 = [value(m[s][i, j, mod, k]) for k in extrPoint[mod]]
                    l = vcat(l, l1)
                end
            end
        else
            l = [value(m[s][i, j, mod]) for j in 0:nUnit for mod in 0:nMod]
        end
        l0 = vcat(l0, l)
    end
    #print(l0)
    push!(df_mdv, l0)
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_mdv.csv", df_mdv)

# By Unit Vars
dvars = ["xLoadU", "xPowGtU", "xFuelU", "xEmisU", "xAuxGtU", "xPowHpU",
"xPowIpU", "xPccRebU", "xDacSbU", "xAllocU", "xAuxStU"]

l = []
for v in dvars
    global l = vcat(l, [v * string(j) => Float64[] for j in 0:nUnit])
end

df_unit = DataFrame(l)

for i in 0:tHorz
    l0 = []
    for v in dvars
        s = Symbol(v)
        l = [value(m[s][i, j]) for j in 0:nUnit]
        l0 = vcat(l0, l)
    end
    push!(df_unit, l0)
end

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

# By Disjunction Vars
dvars = ["xLoadD", "xPowGtD", "xFuelD", "xEmisD", "xAuxGtD", "xPowHpD",
"xPowIpD", "xPccRebD", "xDacSbD", "xAllocD", "xAuxStD"]

l = []
for v in dvars
    global l = vcat(l, [v * string(j) * string(mod) => Float64[]
        for j in 0:nUnit for mod in 0:nMod])
end

df_disj = DataFrame(l)

for i in 0:tHorz
    l0 = []
    for v in dvars
        s = Symbol(v)
        l = [value(m[s][i, j, mod]) for j in 0:nUnit for mod in 0:nMod]
        l0 = vcat(l0, l)
    end
    push!(df_disj, l0)
end

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

max_stars_repo_name: dthierry/princetonDacLti

max_issues_repo_name: dthierry/princetonDacLti

max_forks_repo_name: dthierry/princetonDacLti


text:
using Thesis

# search for "optimal foraging" or "optimal feeding"
oft_pubmed = CSV.read("/home/beatriz/Documents/Thesis/Miscellaneous/OFT_PubMed.csv", DataFrame)
oft_WOS = CSV.read("/home/beatriz/Documents/Thesis/Miscellaneous/OFT_WOS.csv", DataFrame)
chsh = CSV.read("/home/beatriz/Documents/Thesis/Miscellaneous/Charnov_Shoener.csv", DataFrame)

@df oft_pubmed scatter(:Year, :Count, xticks = 1910:10:2021,
    label = "OFT appearence in publications", legend = :topleft, markersize = 2.5)
    plot!([1966,1966],[60,1], arrow = true, label = "1966: MacArtur & Pianka; Emlen")
    plot!([1976,1976],[60,5], arrow = true, label = "1976: Charnov Marginal Value Theorem")
savefig(joinpath("/home/beatriz/Documents/Thesis/Figures","OFTappearencePubMed.pdf"))

@df oft_WOS scatter(:Year, :Count, xticks = 1910:10:2021,
    label = "OFT appearence in publications", legend = :topleft, markersize = 5,
    xlabel = "Year", ylabel = "Publication Count")
    plot!([1966,1966],[100,1], arrow = true, label = "1966: MacArtur & Pianka; Emlen", linewidth = 3)
    plot!([1976,1976],[100,10], arrow = true, label = "1976: Charnov Marginal Value Theorem", linewidth = 3)
savefig(joinpath("/home/beatriz/Documents/Thesis/Figures","OFTappearenceWOS.pdf"))

@df chsh plot(:Years, :Count_Ch, xticks = 1925:10:2020,
    label = "Charnov_1976", legend = :topleft, markersize = 2.5)
@df chsh plot!(:Years, :Count_Sh, label = "Shoener_1974")


text:
export MMD, mmd, mmd2_and_variance

struct MMD{K<:AbstractKernel,D<:MetricOrFun} <: PreMetric
    kernel::K
    dist::D
end

function (m::MMD)(x::AbstractArray, y::AbstractArray)
    xx = kernelsum(m.kernel, x, m.dist)
    yy = kernelsum(m.kernel, y, m.dist)
    xy = kernelsum(m.kernel, x, y, m.dist)
    xx + yy - 2xy
end

"""
	mmd(AbstractKernel(γ), x, y)
	mmd(AbstractKernel(γ), x, y, n)

MMD with Gaussian kernel of bandwidth `γ` using at most `n` samples
"""
mmd(k::AbstractKernel, x::AbstractArray, y::AbstractArray, dist=pairwisel2) =
    MMD(k, dist)(x, y)
mmd(k::AbstractKernel, x::AbstractArray, y::AbstractArray, n::Int, dist=pairwisel2) =
    mmd(k, samplecolumns(x,n), samplecolumns(y,n), dist)

"""
    mmd2_and_variance(K_XX, K_XY, K_YY, unit_diagonal=false, biased=false)

calculates the mmd and variance according to https://arxiv.org/pdf/1611.04488.pdf
"""
function mmd2_and_variance(K_XX, K_XY, K_YY, unit_diagonal=false, biased=false)
    @assert size(K_XX, 2) == size(K_XY, 2) == size(K_XX, 2)
    m = size(K_XX, 2)  # Assumes X, Y are same shape

    ### Get the various sums of kernels that we'll use
    # Kts drop the diagonal, but we don't need to compute them explicitly
    diag_X = diag_Y = 1
    sum_diag_X = sum_diag_Y = m
    sum_diag2_X = sum_diag2_Y = m
    if !unit_diagonal
        diag_X = diag(K_XX)
        diag_Y = diag(K_YY)

        sum_diag_X = sum(diag_X)
        sum_diag_Y = sum(diag_Y)

        sum_diag2_X = dot(diag_X,diag_X)
        sum_diag2_Y = dot(diag_Y,diag_Y)
    end

    Kt_XX_sums = sum(K_XX, dims = 1)[:] - diag_X
    Kt_YY_sums = sum(K_YY, dims = 1)[:] - diag_Y
    K_XY_sums_0 = sum(K_XY; dims = 2)
    K_XY_sums_1 = sum(K_XY, dims = 1)[:]

    Kt_XX_sum = sum(Kt_XX_sums)
    Kt_YY_sum = sum(Kt_YY_sums)
    K_XY_sum = sum(K_XY_sums_0)

    Kt_XX_2_sum = sum(K_XX.^2) - sum_diag2_X
    Kt_YY_2_sum = sum(K_YY.^2) - sum_diag2_Y
    K_XY_2_sum  = sum(K_XY.^2)

    mmd2 = if biased
        ((Kt_XX_sum + sum_diag_X) / (m * m)
              + (Kt_YY_sum + sum_diag_Y) / (m * m)
              - 2 * K_XY_sum / (m * m))
    else
        (Kt_XX_sum / (m * (m-1))
              + Kt_YY_sum / (m * (m-1))
              - 2 * K_XY_sum / (m * m))
    end

    var_est = (
          2 / (m^2 * (m-1)^2) * (
              2 * dot(Kt_XX_sums,Kt_XX_sums) - Kt_XX_2_sum
            + 2 * dot(Kt_YY_sums, Kt_YY_sums) - Kt_YY_2_sum)
        - (4*m-6) / (m^3 * (m-1)^3) * (Kt_XX_sum^2 + Kt_YY_sum^2)
        + 4*(m-2) / (m^3 * (m-1)^2) * (
              dot(K_XY_sums_1, K_XY_sums_1)
            + dot(K_XY_sums_0, K_XY_sums_0))
        - 4 * (m-3) / (m^3 * (m-1)^2) * K_XY_2_sum
        - (8*m - 12) / (m^5 * (m-1)) * K_XY_sum^2
        + 8 / (m^3 * (m-1)) * (
              1/m * (Kt_XX_sum + Kt_YY_sum) * K_XY_sum
            - dot(Kt_XX_sums, K_XY_sums_1)
            - dot(Kt_YY_sums, K_XY_sums_0))
    )

    (mmd2, var_est)
end

"""
    crit_mmd2_var(k, X, Y)
    crit_mmd2_var(k, X, Y, n)

Calculates the criterion for selection of the width of the kernel based on
https://arxiv.org/pdf/1611.04488.pdf using at most `n` samples.
"""
crit_mmd2_var(k::GaussianKernel, X, Y, n::Int, dist = pairwisel2) =
    crit_mmd2_var(k, samplecolumns(X, n), samplecolumns(Y, n), dist)

function crit_mmd2_var(k::GaussianKernel, X, Y, dist = pairwisel2)
    KXX = k.(dist(X,X))
    KXY = k.(dist(X,Y))
    KYY = k.(dist(Y,Y))
    m, v = mmd2_and_variance(KXX, KXY, KYY, true)
    m / sqrt(max(v, eps(Float64)))
end

crit_mmd2_var(k::IMQKernel, X, Y, n::Int, dist = pairwisel2) =
    crit_mmd2_var(k, samplecolumns(X, n), samplecolumns(Y, n), dist)

function crit_mmd2_var(k::IMQKernel, X, Y, dist = pairwisel2)
    KXX = k.(dist(X,X))
    KXY = k.(dist(X,Y))
    KYY = k.(dist(Y,Y))
    m, v = mmd2_and_variance(KXX, KXY, KYY, false)
    m / sqrt(max(v, eps(Float64)))
end

"""
    crit_mxy_over_mltpl(k, x, y)

Calculates criterion for selection of the kernel width as MMD(X, Y) / SQRT(MMD(X,X) * MMD(Y, Y))
"""
crit_mxy_over_mltpl(k, x, y, dist = pairwisel2) =
    mmd(k, x, y, dist) / sqrt(abs(mmd(k, split2(x)..., dist) * mmd(k, split2(y)..., dist)))

"""
    crit_mxy_over_sum(k, x, y)

Calculates criterion for selection of the kernel width as MMD(X, Y) / (MMD(X,X) + MMD(Y, Y))
"""
crit_mxy_over_sum(k, x, y, dist = pairwisel2) =
    mmd(k, x, y, dist) / (abs(mmd(k, split2(x)..., dist)) + abs(mmd(k, split2(y)..., dist)))

function mmdfromdist(k::AbstractKernel, d, idx, cidx)
    _z = eltype(d)(0)
	li, lc = length(idx), length(cidx)
	kxx = (mapsum(k, d, idx) - li*k(_z))/(li^2 - li)
	kyy = (mapsum(k, d, cidx) - lc*k(_z))/(lc^2 - lc)
	kyx = mapsum(k, d, idx, cidx)/ (lc*li)
	kxx + kyy -2kyx
end

function mapsum(f, d, idx)
    s = eltype(d)(0)
    @inbounds for j in idx, i in idx
        s+= f(d[i,j])
    end
    return s
end

function mapsum(f, d, rowidx, colidx)
    s = eltype(d)(0)
    @inbounds for i in rowidx, j in colidx
        s+= f(d[i,j])
    end
    return s
end


text:
let
    # Include inside this quote any rules that should have FastMath versions
    fastable_ast = quote
        #  Trig-Basics
        @scalar_rule cos(x) -(sin(x))
        @scalar_rule sin(x) cos(x)
        @scalar_rule tan(x) 1 + Ω ^ 2

        # Trig-Hyperbolic
        @scalar_rule cosh(x) sinh(x)
        @scalar_rule sinh(x) cosh(x)
        @scalar_rule tanh(x) 1 - Ω ^ 2

        # Trig- Inverses
        @scalar_rule acos(x) -(inv(sqrt(1 - x ^ 2)))
        @scalar_rule asin(x) inv(sqrt(1 - x ^ 2))
        @scalar_rule atan(x) inv(1 + x ^ 2)

        # Trig-Multivariate
        @scalar_rule atan(y, x) @setup(u = x ^ 2 + y ^ 2) (x / u, -y / u)
        @scalar_rule sincos(x) @setup((sinx, cosx) = Ω) cosx -sinx

        # exponents
        @scalar_rule cbrt(x) inv(3 * Ω ^ 2)
        @scalar_rule inv(x) -(Ω ^ 2)
        @scalar_rule sqrt(x) inv(2Ω)
        @scalar_rule exp(x::Real) Ω
        @scalar_rule exp10(x) Ω * log(oftype(x, 10))
        @scalar_rule exp2(x) Ω * log(oftype(x, 2))
        @scalar_rule expm1(x) exp(x)
        @scalar_rule log(x) inv(x)
        @scalar_rule log10(x) inv(x) / log(oftype(x, 10))
        @scalar_rule log1p(x) inv(x + 1)
        @scalar_rule log2(x) inv(x) / log(oftype(x, 2))

        # Unary complex functions
        @scalar_rule abs(x::Real) sign(x)
        @scalar_rule abs2(x) 2x
        @scalar_rule angle(x::Real) Zero()
        @scalar_rule conj(x::Real) One()

        # Binary functions
        @scalar_rule hypot(x::Real, y::Real) (x / Ω, y / Ω)
        @scalar_rule x + y (One(), One())
        @scalar_rule x - y (One(), -1)
        @scalar_rule x / y (inv(y), -((x / y) / y))
        #log(complex(x)) is require so it give correct complex answer for x<0
        @scalar_rule(x ^ y,
            (ifelse(iszero(y), zero(Ω), y * x ^ (y - 1)), Ω * log(complex(x))),
        )
        @scalar_rule(
            rem(x, y),
            @setup((u, nan) = promote(x / y, NaN16), isint = isinteger(x / y)),
            (ifelse(isint, nan, one(u)), ifelse(isint, nan, -trunc(u))),
        )
        @scalar_rule max(x, y) @setup(gt = x > y) (gt, !gt)
        @scalar_rule min(x, y) @setup(gt = x > y) (!gt, gt)

        # Unary functions
        @scalar_rule +x One()
        @scalar_rule -x -1

        @scalar_rule sign(x) Zero()

        # product rule requires special care for arguments where `mul` is non-commutative
        function frule((_, Δx, Δy), ::typeof(*), x::Number, y::Number)
            # Optimized version of `Δx .* y .+ x .* Δy`. Also, it is potentially more
            # accurate on machines with FMA instructions, since there are only two
            # rounding operations, one in `muladd/fma` and the other in `*`.
            ∂xy = muladd.(Δx, y, x .* Δy)
            return x * y, ∂xy
        end

        function rrule(::typeof(*), x::Number, y::Number)
            function times_pullback(ΔΩ)
                return (NO_FIELDS,  @thunk(ΔΩ * y), @thunk(x * ΔΩ))
            end
            return x * y, times_pullback
        end
    end

    # Rewrite everything to use fast_math functions, including the type-constraints
    eval(Base.FastMath.make_fastmath(fastable_ast))
    eval(fastable_ast)  # Get original definitions
    # we do this second so it overwrites anything we included by mistake in the fastable
end


text:
#!/usr/bin/env julia

# An example for using Gadfly within GTK

using Gtk
using Cairo
using Compose
using Gadfly
using GtkBuilderAid

mutable struct PlotData
  surface::CairoSurface
  plot
end

function draw_backbuffer(area, udata)
  udata.surface = create_similar_surface(area, Gtk.GEnum(Cairo.CONTENT_COLOR_ALPHA))
  Gadfly.draw(
    Compose.CAIROSURFACE(udata.surface, CairoContext(udata.surface)),
    udata.plot)
end

ft = @GtkFunctionTable begin

  @guarded Cint(1) function plot_area_draw_cb(area, ctx_ptr, qsdata)
    udata = qsdata.userdata
    ctx = CairoContext(ctx_ptr)
    set_source_surface(ctx, udata.surface)
    paint(ctx)
    Cint(0)
  end

  @guarded Cint(0) function plot_area_configure_event_cb(
      area,
      configure_event,
      qsdata)
    # redraw the plot
    draw_backbuffer(area, qsdata.userdata)
    Cint(1)
  end

  @guarded function plot_area_realize_cb(
      area,
      qsdata)
    draw_backbuffer(area, qsdata.userdata)
    nothing
  end

end

start_application(
    ft,
    "io.github.matt5sean3.GtkBuilderAid.third",
    "main_window",
    "resources/main.ui",
    PlotData(
        CairoSurface(C_NULL, -1, -1),
        plot(x=collect(1:10), y=rand(10), Geom.LineGeometry)))


text:
module PyGen

using MacroTools

export @pygen

if VERSION > v"0.6-"

macro pygen(f)

    # generate a new symbol for the channel name and
    # wrapper function
    c = gensym()
    η = gensym()

    # yield(λ) → put!(c, λ)
    f′ = MacroTools.postwalk(f) do x
        @capture(x, yield(λ_)) || return x
        return :(put!($c, $λ))
    end

    # Fetch the function name and args
    @capture(f′, function func_(args__) body_ end)

    # wrap up the η function
    final = quote
        function $func($(args...))
            function $η($c)
                $body
            end
            return Channel(c -> $η(c))
        end
    end

    return esc(final)
end

else  # for versions less than v0.6-

macro pygen(f)
    η = gensym()

    f′ = MacroTools.postwalk(f) do x
        @capture(x, yield(λ_)) || return x
        return :(produce($λ))
    end

    @capture(f′, function func_(args__) body_ end)

    final = quote
        function $func($(args...))
            function $η()
                $body
            end
            return Task($η)
        end
    end

    return esc(final)
end

end

"""
    @pygen

Making julia walk a little *more* like python. Stick
the @pygen macro in front of a function declaration and
you can use `yield(foo)` statements to make a python style
generator!

## Examples

```julia
julia> @pygen function fibonacci()
           n, m = 0, 1
           while true
               yield(m)
               n, m = m, n + m
           end
       end
fibonacci (generic function with 1 method)

julia> for n in fibonacci()
           println(n)
           sleep(1)
       end
1
1
2
3
5
8
13
.
.
.
```
"""
:(@pygen)

end


text:
module Oracles

using Memento
using Distances
using Statistics
using Serialization

import SVDD
import Distributions
import LIBSVM
import MLKernels

import GaussianMixtures: GMM
import MLBase: StratifiedKfold

import ..OneClassActiveLearning:
    ConfusionMatrix,

    cohens_kappa,
    convert_labels_to_learning,
    load_data,
    matthews_corr

include("oracle_base.jl")

include("PoolOracle.jl")
include("NoisyOracle.jl")
include("QuerySynthesisFunctionOracle.jl")
include("QuerySynthesisGMMOracle.jl")
include("QuerySynthesisKNNOracle.jl")
include("QuerySynthesisOCCOracle.jl")
include("QuerySynthesisSVMOracle.jl")
include("QuerySynthesisCVWrapperOracle.jl")

export
    Oracle,
    PoolOracle,
    NoisyOracle,
    QuerySynthesisFunctionOracle,
    QuerySynthesisGMMOracle,
    QuerySynthesisKNNOracle,
    QuerySynthesisOCCOracle,
    QuerySynthesisSVMOracle,
    QuerySynthesisCVWrapperOracle,

    ask_oracle,
    initialize_oracle
end


text:
push!(LOAD_PATH, "../")
using QXTools, Documenter

makedocs(;
    modules=[QXTools],
    authors="QuantEx team",
    repo="https://github.com/JuliaQX/QXTools.jl/blob/{commit}{path}#L{line}",
    sitename="QXTools.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaQX.github.io/QXTools.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Getting Started" => "getting_started.md",
        "User's Guide" => "users_guide.md",
        "Tutorials" => ["basics.md", "circuits.md", "features.md", "distributed.md"],
        "LICENSE" => "license.md"
    ],
)

deploydocs(;
    repo="github.com/JuliaQX/QXTools.jl",
)


text:
module LayeredLayouts
using Dates
using JuMP
using ECOS
using Cbc
using IterTools: IterTools
using LightGraphs
using Random

export LayeredMinDistOne, Zarate
export solve_positions

abstract type AbstractLayout end

include("utils.jl")

include("graph_properties.jl")
include("layering.jl")

include("zarate.jl")

end


text:
using Base.Test, FactCheck
using CloudGraphs_RegisteredTypes

# Creating a connection
print("[TEST] Connecting to the local CloudGraphs instance (Neo4j and Mongo)...");
configuration = CloudGraphs.CloudGraphConfiguration("localhost", 7474, "neo4j", "neo5j", "localhost", 27017, false, "", "");
cloudGraph = connect(configuration);
println("Success!");

facts("Standard SlamInDb type registration") do
    registerGeneralVariableTypes!(cloudGraph)
end


text:
"""
# Quantum algorithms

- Qubits, superposition, and measurement
- The plan
- The quantum Fourier transform
- Periodicity
    - The Fourier transform of a periodic vector
- Quantum circuits
    - Elementary quantum gates
    - Two basic types of quantum circuits
    - Quantum Fourier Transform
    - Classical Functions
- Factoring as periodicity
    - Setting up a periodic superposition
    - Quantum physics meets computation
"""
module Chapter10

end # module