text:
# vim: set wrap
#: by David Thierry 2021
#: Set all the eqns to hslice - 1
using JuMP
#using Gurobi
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_pow_c = DataFrame(CSV.File("../resources/test_elec.csv"))
df_ng_c = DataFrame(CSV.File("../resources/natural_gas_price.csv"))

#: Assign parameters
#: Usually two parameters a & b.
bPowGasTeLoad = df_gas[1, 2]
aPowGasTeLoad = df_gas[1, 3]

bFuelEload = df_gas[2, 2]
aFuelEload = df_gas[2, 3]

#: Conversion factor
lbcoToTonneco = 0.4535924 / 1000

bEmissFactEload = df_gas[3, 2] * lbcoToTonneco
aEmissFactEload = df_gas[3, 3] * lbcoToTonneco

bPowHpEload = df_gas[4, 2] / 1000  #: 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

#: Base PCC steam from the full 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]

#: Side steam is the difference btw full steam and power
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

# 7 GJ/tonneCO for regeneration, put a factor here (5)
aSteaUseRateDacAir = 5 * (7 * 1e+06 / 3600) * kwhToMmbtu
aSteaUseRateDacFlue = 5 * (7 * 1e+06 / 3600) * kwhToMmbtu

aPowUseRateDacAir = 500 / 1000
aPowUseRateDacFlue = 250 / 1000

# assume 1 mmol/gSorb
# gCo/gSorb
gCogSorbRatio = 1e-03 * 44.0095

aSorbCo2CapFlue = 1 * gCogSorbRatio
aSorbCo2CapAir = gCogSorbRatio

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

aCapRatePcc = 0.97

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 dac-states have the 0th

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

# Dictionaries with the disjunction parameter, e.g., off, cold start, warm start, dispatchs
slopeFactorV1 = 0.0164678362573099
bFactorV1 = 0.0194931773879139

## Slope
aPowGt = Dict(0 => 1,
              1 => (aPowGasTeLoad * 60 + bPowGasTeLoad) * slopeFactorV1,
              2 => (aPowGasTeLoad * 60 + bPowGasTeLoad) * slopeFactorV1,
              3 => aPowGasTeLoad
              )

aFuel = Dict(0 => 1,
             1 => (aFuelEload * 60 + bFuelEload) * slopeFactorV1,
             2 => (aFuelEload * 60 + bFuelEload) * slopeFactorV1,
             3 => aFuelEload
             )

aEmis = Dict(0 => 1,
             1 => (aEmissFactEload * 60 + bEmissFactEload) * slopeFactorV1,
             2 => (aEmissFactEload * 60 + bEmissFactEload) * slopeFactorV1,
             3 => aEmissFactEload
             )

aAuxGt = Dict(0 => 0,
              1 => 0, #(aAuxRateGas * 60 + bAuxRateGas) * slopeFactorV1,
              2 => 0, #(aAuxRateGas * 60 + bAuxRateGas) * slopeFactorV1,
              3 => aAuxRateGas
              )

aPowHp = Dict(0 => 0,
              1 => 0,
              2 => 0,
              3 => aPowHpEload
              )

aPowIp = Dict(0 => 0,
              1 => 0,
              2 => 0,
              3 => aPowIpEload
              )

aCcReb = Dict(
              0 => 0,
              1 => 0,
              2 => 0,
              3 => aCcRebDutyEload
              )

aDacSb = Dict(0 => 0,
              1 => 0,
              2 => 0,
              3 => aDacSteaBaseEload
              )

aSideS = Dict(0 => 0,
              1 => 0,
              2 => 0,
              3 => aSideSteaEload
              )

aAuxSt = Dict(0 => 0,
              1 => 0,
              2 => 0,
              3 => 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 => (aPowGasTeLoad * 60 + bPowGasTeLoad) * bFactorV1,
              2 => (aPowGasTeLoad * 60 + bPowGasTeLoad) * bFactorV1,
              3 => bPowGasTeLoad
              )
bFuel = Dict(
             0 => 0,
             1 => (aFuelEload * 60 + bFuelEload) * bFactorV1,
             2 => (aFuelEload * 60 + bFuelEload) * bFactorV1,
             3 => bFuelEload
             )
bEmis = Dict(
             0 => 0,
             1 => (aEmissFactEload * 60 + bEmissFactEload) * bFactorV1,
             2 => (aEmissFactEload * 60 + bEmissFactEload) * bFactorV1,
             3 => bEmissFactEload
             )

bAuxGt = Dict(
              0 => 0,
              1 => 0, #(aAuxRateGas * 60 + bAuxRateGas) * bFactorV1,
              2 => 0, #(aAuxRateGas * 60 + bAuxRateGas) * bFactorV1,
              3 => bAuxRateGas
              )
bPowHp = Dict(
              0 => 0,
              1 => 0,
              2 => 0,
              3 => bPowHpEload)
bPowIp = Dict(
              0 => 0,
              1 => 0,
              2 => 0,
              3 => bPowIpEload
    )
bCcReb = Dict(
              0 => 0,
              1 => 0,
              2 => 0,
              3 => bCcRebDutyEload)
bDacSb = Dict(
              0 => 0,
              1 => 0,
              2 => 0,
              3 => bDacSteaBaseEload)
bSideS = Dict(
              0 => 0,
              1 => 0,
              2 => 0,
              3 => bSideSteaEload)

bAuxSt = Dict(
              0 => 0,
              1 => 0,
              2 => 0,
              3 => bAuxRateStea
             )

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

# USD/MWh
# Some relevant price signals
#pow_price =(df_pow_c[!, "Test4"])  # USD/MWh
pow_price =(df_pow_c[!, "MiNg_150_ERCOT"])  # USD/MWh
#pow_price =(df_pow_c[!, "MiNg_150_NYISO"])  # USD/MWh
#pow_price =(df_pow_c[!, "MiNg_150_CAISO"])  # USD/MWh
#pow_price =(df_pow_c[!, "MiNg_150_MISO-W"])  # USD/MWh
#pow_price =(df_pow_c[!, "MMiNg_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
#CAISO	2.26	2.26
#ERCOT	2.64	2.64
#MISO-W	1.69	2.01
#NYISO	1.12	1.14
#PJM-W	1.42	1.43

#: Create model
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

# Capital Cost DAC USD/tCo2/yr
cCostInvDacUsdtCo2yr = 750
cCostFixedDacUsdtCo2yr = 25
cCostVariableDacUsdtCo2yr = 12

#: Number of modes and units
nMod = 3  # We go from 0 to 3, i.e. 4 modes
nUnit = 1  # We go from 0 to 1, i.e. 2 units

#:
coldSidx = 1
warmSidx = 2
dsptchIdx = 3

# 0 -> off
# 1 -> cold-start
# 2 -> warm-start
# 3 -> on

# Extreme points per mode
extrPoint = Dict(0 => 0, 1 => 0, 2 => 0, 3 => 0:1)

# Convex weight
@variable(m,
           0 <= lambda[0:tHorz, 0:nUnit, i=0:nMod, extrPoint[i]; i == dsptchIdx] <= 1)
# mode Off/ColdS/WarmS/etc..
@variable(m, y[-1:tHorz, 0:nUnit, 0:nMod], Bin, start = 0)  # On and off
# Transition j1 to j2
@variable(m, z[0:tHorz, 0:nUnit, j1 = 0:nMod, j2 = 0:nMod; j1 != j2], Bin, start = 0)

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

# "overall variables power, fuel, etc." i.e., U1 + U2"
@variable(m, 0 <= xActualLoad[0:tHorz] <= 100)

@variable(m, 0 <= xPowGasTur[0:tHorz]) # gas turbine power
@variable(m, 0 <= xPowGross[0:tHorz]) # gross power
@variable(m, 0 <= xPowOut[0:tHorz]) # power out

@variable(m, 0 <= xAuxPowGasT[0:tHorz]) # auxiliary power gas turbine

@variable(m, 0 <= xPowHp[0:tHorz]) # high pressure power
@variable(m, 0 <= xPowIp[0:tHorz]) # interm. pressure power
@variable(m, 0 <= xPowLp[0:tHorz]) # low pressure power

@variable(m, 0 <= xFuel[0:tHorz]) # fuel
@variable(m, 0 <= xCo2Fuel[0:tHorz]) # co2 emission
@variable(m, 0 <= xDacSteaDuty[0:tHorz]) # dac steam

@variable(m, 0 <= xCcRebDuty[0:tHorz]) # pcc duty
@variable(m, 0 <= xDacSteaBaseDuty[0:tHorz]) # dac base steam

@variable(m, 0 <= xAllocSteam[0:tHorz]) # heat allocation
@variable(m, 0 <= xSteaPowLp[0:tHorz]) # steam for lp power
@variable(m, 0 <= xSideSteaDac[0:tHorz]) # side steam dac

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

# 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]) # co2 capacity pcc
@variable(m, 0 <= xSteaUsePcc[0:tHorz]) # steam used by pcc
@variable(m, 0 <= xPowUsePcc[0:tHorz]) # power used by pcc
@variable(m, 0 <= xCo2PccOut[0:tHorz]) # co2 out pcc

@variable(m, 0 <= vCo2PccVent[0:tHorz]) # vented co2 by pcc
@variable(m, 0 <= xCo2DacFlueIn[0:tHorz]) # flue co2 to gas
@variable(m, 0 <= xPccSteaSlack[0:tHorz]) # steam slack for pcc

# Dac-Flue, takes one time slot for adsorption
@variable(m, 0 <= xA0Flue[0:tHorz, 0:hSlice-1]) # dac input
@variable(m, 0 <= xA1Flue[0:tHorz, 0:hSlice]) # dac State

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

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

@variable(m, 0 <= xCo2CapDacFlue[0:tHorz, 0:hSlice-1]) # co2 from dac-flue
@variable(m, 0 <= xSteaUseDacFlue[0:tHorz, 0:hSlice-1]) # steam use dac-flue
@variable(m, 0 <= xPowUseDacFlue[0:tHorz, 0:hSlice-1]) # pow use dac-flue
@variable(m, 0 <= xCo2DacVentFlue[0:tHorz]) # co2 vented dac

# Dac-Air, takes two time slots for adsorption
@variable(m, 0 <= xA0Air[0:tHorz, 0:hSlice-1]) # dac input
@variable(m, 0 <= xA1Air[0:tHorz, 0:hSlice]) # dac State
@variable(m, 0 <= xA2Air[0:tHorz, 0:hSlice]) # dac State

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

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

@variable(m, 0 <= xCo2CapDacAir[0:tHorz, 0:hSlice-1]) # co2 from air
@variable(m, 0 <= xSteaUseDacAir[0:tHorz, 0:hSlice-1]) # steam use dac-air
@variable(m, 0 <= xPowUseDacAir[0:tHorz, 0:hSlice-1]) # power use dac-air

@variable(m, 0 <= xDacSteaSlack[0:tHorz])
# DAC hourly capacity
#
# CO2 compression
@variable(m, 0 <= xCo2Comp[0:tHorz]) # co2 compression
@variable(m, 0 <= xPowUseComp[0:tHorz]) # power use co2 compression
@variable(m, xCo2Vent[0:tHorz])  # This used to be only positive

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

# Constraints
# Op mode
# Down times
# Up times

#: feasible space (load basis)
extreme_d_down = [0]
extreme_d_cold = [60.]
extreme_d_warm = [60.]
extreme_d_on = [60.0, 100.0] #: dispatch

ep_list = Dict(0 => extreme_d_down,
               1 => extreme_d_cold,
               2 => extreme_d_warm,
               3 => extreme_d_on)
# Convex combination (only for dispatch, the other modes don't actually use this)
@constraint(m, cConvxEq[i = 0:tHorz, j = 0:nUnit],
            sum(lambda[i, j, dsptchIdx, k] * ep_list[dsptchIdx][k + 1] for k in extrPoint[dsptchIdx])
            == xLoadD[i, j, dsptchIdx]  # There is only a single extreme
           )

# Time to transition
# 5sync 5soak at cold
# 3sync 3soak at warm
#: Sync
kColdSync = 5
kWarmSync = 3
#: Soak
kColdSoak = 5
kWarmSoak = 3

#: Trajectory for load
delColdS = (60. - 20.)/kColdSoak
delWarmS = (60. - 20.)/kWarmSoak

#: Cold trajectory
FdefCold = Dict()
j = 0
for i in 0:(kColdSync + kColdSoak)
    if i <= kColdSync - 1
        FdefCold[i] = 0 # flat for sync
    else
        FdefCold[i] = delColdS * j + 20.
        global j += 1
    end
end

#: Warm trajectory
FdefWarm = Dict()
j = 0
for i in 0:(kWarmSync + kWarmSoak)
    if i <= kWarmSync - 1
        FdefWarm[i] = 0 # flat for sync
    else
        FdefWarm[i] = delWarmS * j + 20.
        global j += 1
    end
end

println("FdefCold")
for i in 0:(kColdSync + kColdSoak)
    println(i, "\t", FdefCold[i])
end

println("FdefWarm")
for i in 0:(kWarmSync + kWarmSoak)
    println(i, "\t", FdefWarm[i])
end

Fdef = Dict()
Fdef[0] = Dict(1 => 0)
Fdef[coldSidx] = FdefCold
Fdef[warmSidx] = FdefWarm

#: Overall time of mode
Kmm = Dict()
Kmm[0] = 1
Kmm[coldSidx] = kColdSync + kColdSoak
Kmm[warmSidx] = kWarmSync + kWarmSoak
Kmm[dsptchIdx] = 0

println("Kmm")
for i in 0:nMod
    println(i, "\t", Kmm[i])
end

# Trajectory of load for load cold/warm-start
@constraint(m, cTrajEqWarm[i = 0:tHorz, j = 0:nUnit, mode = [coldSidx, warmSidx]],
            xLoadD[i, j, mode] ==
            sum((z[th, j, 0, mode]) * Fdef[mode][i - th + 1] for th in (i - Kmm[mode] + 1):i if th >= 0)
           )

# off mode
@constraint(m, cOffEq[i = 0:tHorz, j = 0:nUnit], xLoadD[i, j, 0] == 0)

# Lambda constraint (dispatch only)
@constraint(m, xLambdaEq[i = 0:tHorz, j = 0:nUnit, mod = 0:nMod; mod == dsptchIdx],
    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])* y[i, j, mod]
    )

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

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

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

# Forbidden Transitions
# from off to full
@constraint(m, forbConOffFullEq[i = 0:tHorz, j = 0:nUnit],
    z[i, j, 0, dsptchIdx] == 0)

# from dispatch to cold
@constraint(m, forbConEqOnBack1[i = 0:tHorz, j = 0:nUnit],
    z[i, j, dsptchIdx, coldSidx] == 0)
# from distpatch to warm
@constraint(m, forbConEqOnBack2[i = 0:tHorz, j = 0:nUnit],
    z[i, j, dsptchIdx, warmSidx] == 0)

# from cold to warm
@constraint(m, forbConEqCsWs[i = 0:tHorz, j = 0:nUnit],
    z[i, j, coldSidx, warmSidx] == 0)
# from warm to cold
@constraint(m, forbConEqWsCs[i = 0:tHorz, j = 0:nUnit],
    z[i, j, warmSidx, coldSidx] == 0)

# from cold to off
@constraint(m, forbConEqCoff[i = 0:tHorz, j = 0:nUnit],
    z[i, j, coldSidx, 0] == 0)
# from warm to off
@constraint(m, forbConEqWoff[i = 0:tHorz, j = 0:nUnit],
    z[i, j, warmSidx, 0] == 0)

#: Fuel trajectory
## Fuel reference points
f20 = (aFuel[coldSidx] * 20 + bFuel[coldSidx]) * 0.5
f60 = (aFuel[coldSidx] * 60 + bFuel[coldSidx]) * 0.5
f100dsp= (aFuel[dsptchIdx] * 100 + bFuel[dsptchIdx]) * 0.5

println("f20\t", f20)
println("f60\t", f60)
println("f100\t", f100dsp)

## Delta for fuel trajectory
dColdF = (f60 - f20)/kColdSoak
dWarmF = (f60 - f20)/kWarmSoak

FdefCfuel = Dict()
j = 0
for i in 0:(kColdSync + kColdSoak)
    if i <= kColdSync - 1
        FdefCfuel[i] = f20 * 0.1 # flat for sync
    else
        FdefCfuel[i] = dColdF * j + f20
        global j += 1
    end
end

FdefWfuel = Dict()
j = 0
for i in 0:(kWarmSync + kWarmSoak)
    if i <= kWarmSync - 1
        FdefWfuel[i] = f20 * 0.1 # flat for sync
    else
        FdefWfuel[i] = dWarmF * j + f20
        global j += 1
    end
end

println("FdefCold Fuel")
for i in 0:(kColdSync + kColdSoak)
    println(i, "\t", FdefCfuel[i])
end

println("FdefWarm Fuel")
for i in 0:(kWarmSync + kWarmSoak)
    println(i, "\t", FdefWfuel[i])
end

# dictonary with the trajectory
FdefFuel = Dict()
FdefFuel[0] = Dict(1 => 0)
FdefFuel[coldSidx] = FdefCfuel
FdefFuel[warmSidx] = FdefWfuel

# fuel for dispatch
@constraint(m, fuelDspEq[i = 0:tHorz, j = 0:nUnit],
    xFuelD[i, j, dsptchIdx] ==
    (aFuel[dsptchIdx] * xLoadD[i, j, dsptchIdx] + bFuel[dsptchIdx] * y[i, j, dsptchIdx]) * 0.5
    )
# fuel for offline
@constraint(m, fuelOffEq[i = 0:tHorz, j = 0:nUnit],
               xFuelD[i, j, 0] ==
               0
               )
# fuel for cold/warm-start
# FdefFuel already includes the 0.5 factor
@constraint(m, fuelTrajEq[i = 0:tHorz, j = 0:nUnit, mod = [coldSidx, warmSidx]],
            xFuelD[i, j, mod] ==
            sum((z[th, j, 0, mod]) * FdefFuel[mod][i - th + 1] for th in (i - Kmm[mod] + 1):i if th >= 0)
           )
# big-M
bMf = Dict(0=> 0, 1 => f60, 2 => f60, 3 => f100dsp)
@constraint(m, fuelBmEq[i = 0:tHorz, j = 0:nUnit, mod = 0:nMod],
               xFuelD[i, j, mod] <= bMf[mod] * y[i, j, mod]
            )
###
#: Co2 trajectory
#0.061919   0.205418
afuelEm = 0.055907879 #: from the new data
bfuelEm = 0.205417634
# em [=] tonneCO2
# fuel [=] MMBTU / hr

e20 = afuelEm * f20 + bfuelEm
e60 = afuelEm * f60 + bfuelEm
e100dsp = afuelEm * f100dsp + bfuelEm

println("e20\t", e20)
println("e60\t", e60)
println("e100\t", e100dsp)

## Delta co2 emission
dColdE = (e60 - e20)/kColdSoak
dWarmE = (e60 - e20)/kWarmSoak

FdefCem = Dict()
j = 0
for i in 0:(kColdSync + kColdSoak)
    if i <= kColdSync - 1
        FdefCem[i] = e20 * 0.1 # flat during sync
    else
        FdefCem[i] = dColdE * j + e20
        global j += 1
    end
end

FdefWem = Dict()
j = 0
for i in 0:(kWarmSync + kWarmSoak)
    if i <= kWarmSync - 1
        FdefWem[i] = e20 * 0.1 # flat during sync
    else
        FdefWem[i] = dWarmE * j + e20
        global j += 1
    end
end

println("FdefCold Em")
for i in 0:(kColdSync + kColdSoak)
    println(i, "\t", FdefCem[i])
end

println("FdefWarm Em")
for i in 0:(kWarmSync + kWarmSoak)
    println(i, "\t", FdefWem[i])
end

# Dictonary with the trajectories
FdefEm = Dict()
FdefEm[0] = Dict(1 => 0)
FdefEm[coldSidx] = FdefCem
FdefEm[warmSidx] = FdefWem

# co2 for dispatch
@constraint(m, emDspEq[i = 0:tHorz, j = 0:nUnit],
    xEmisD[i, j, dsptchIdx] ==
    (aEmis[dsptchIdx] * xLoadD[i, j, dsptchIdx] + bEmis[dsptchIdx] * y[i, j, dsptchIdx]) * 0.5
    )
# co2 for offline
@constraint(m, emOffEq[i = 0:tHorz, j = 0:nUnit],
               xEmisD[i, j, 0] ==
               0
               )
# co2 for coldstart/warmstart
@constraint(m, emTrajEq[i = 0:tHorz, j = 0:nUnit, mod = [coldSidx, warmSidx]],
            xEmisD[i, j, mod] ==
            sum((z[th, j, 0, mod]) * FdefEm[mod][i - th + 1] for th in (i - Kmm[mod] + 1):i if th >= 0)
           )
# big-M
bMe = Dict(0=> 0, 1 => e60, 2 => e60, 3 => e100dsp)
@constraint(m, emBmEq[i = 0:tHorz, j = 0:nUnit, mod = 0:nMod],
               xEmisD[i, j, mod] <= bMe[mod] * y[i, j, mod]
            )
#
##
# The rest of the disagregated variable (0.5 because we want half)
@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, 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)
    )
#: Overall variables (unit1 + unit2 + .., etc...)
@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)
           )
#
@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, 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, auxPowSteaTEq[i = 0:tHorz],
            xAuxPowSteaT[i] == sum(xAuxStU[i, j] for j in 0:nUnit)
           )

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

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

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

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

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

# (arbitrary) costs of shutdown/startup
oneHourofGas60Load = cNgPerMmbtu * (bFuelEload + 60 * aFuelEload)

shutdownCostOn = Dict(0 => oneHourofGas60Load * 2,
                      1 => oneHourofGas60Load)
startupCostCold = Dict(0 => oneHourofGas60Load,
                       1 => oneHourofGas60Load/2)
startupCostWarm = Dict(0 => oneHourofGas60Load,
                       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)
    +
        # Cost per unit
    sum(shutdownCostOn[j] * z[i, j, dsptchIdx, 0]
        + startupCostCold[j] * z[i, j, 0, coldSidx]
        + startupCostWarm[j] * z[i, j, 0, warmSidx]
        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, Gurobi.Optimizer)
set_optimizer(m, SCIP.Optimizer)
#set_optimizer_attribute(m, "branching/preferbinary", true)

# solve
optimize!(m)
println(termination_status(m))
println(solve_time(m))

# Turn on remaining constraints
# State-dependent transitions
# If Warm implies On -> Off in the window of 12 hours
@constraint(m, sdTransOffWarmOn[i = 0:tHorz, j = 0:nUnit],
               z[i, j, 0, warmSidx] <= sum(z[i - k, j, dsptchIdx, 0]
                                            for k in (2-1):12 if (i - k) >= 0)
               )

# Minimum stay
# off to cold
# Now from 0th
# It takes 5 hours to sync
@constraint(m, minStayCold1[i = 0:tHorz, j = 0:nUnit],
    y[i, j, coldSidx] == sum(z[i - k, j, 0, coldSidx] for k in 0:(10 - 1) if (i-k) >= 0)
    )
@constraint(m, minStayCold2[i = 0:tHorz, j = 0:nUnit],
    y[i, j, warmSidx] == sum(z[i - k, j, 0, warmSidx] for k in 0:(6 - 1) if (i-k) >= 0)
    )

# optimize again
optimize!(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

#: Write dataframes

# 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"
                local l = [v * string(j) * string(mod) * string(m2) => Float64[]
                for m2 in 0:nMod if mod != m2]
            else
                local 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
    local l0 = []
    for v in bvars_names
        s = Symbol(v)
        if v == "z"
            local l = [value(m[s][i, j, m1, m2]) for j in 0:nUnit
             for m1 in 0:nMod  for m2 in 0:nMod if m1 != m2]
        else
            local 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"
                if mod == dsptchIdx
                    local l = [v * string(j) * string(mod) * string(k) => Float64[]
                    for k in extrPoint[mod]]
                else
                    local l = []
                end
            else
                local 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
    local l0 = []
    for v in dis_names
        s = Symbol(v)
        if v == "lambda"
            local l = []
            for j in 0:nUnit
                for mod in [dsptchIdx]
                    local 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
    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
    local l0 = []
    for v in dvars
        s = Symbol(v)
        local 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
    local l0 = []
    for v in dvars
        s = Symbol(v)
        local 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)

println(solve_time(m))

max_stars_repo_name: dthierry/princetonDacLti

max_issues_repo_name: dthierry/princetonDacLti

max_forks_repo_name: dthierry/princetonDacLti


text:
include("euler/euler.jl")
using .BigIntegers: BigInteger

compute()::BigInteger = 28_433 * powermod(BigInteger(2), 7_830_457, 10_000_000_000) % 10_000_000_000 + 1

ext: jl

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
export Palette,
    Gray

mutable struct Palette
    # SDL requires Int64 not UInt8
    r::Int64
    g::Int64
    b::Int64
    a::Int64

    function Palette()
        new(0, 0, 0, 255)
    end

    function Palette(r::Integer, g::Integer, b::Integer)
        new(Int64(r), Int64(g), Int64(b), 255)
    end

    function Palette(r::Integer, g::Integer, b::Integer, a::Integer)
        new(Int64(r), Int64(g), Int64(b), Int64(a))
    end

    # 0xrrggbbaa
    function Palette(color::UInt32)
        r = Integer((color & 0xff000000) >> 24)
        g = Integer((color & 0x00ff0000) >> 16)
        b = Integer((color & 0x0000ff00) >> 8)
        a = Integer(color & 0x000000ff)
        new(Int64(r), Int64(g), Int64(b), Int64(a))
    end
end

function set_red!(color::Palette, r::Int64)
    color.r = r
end

function set_green!(color::Palette, g::Int64)
    color.g = g
end

function set_blue!(color::Palette, b::Int64)
    color.b = b
end

function set_alpha!(color::Palette, a::Int64)
    color.a = a
end

function White()
    Palette(255, 255, 255)
end

function Black()
    Palette()
end

function Red()
    Palette(255, 0, 0)
end

function Green()
    Palette(0, 255, 0)
end

function Blue()
    Palette(0, 0, 255)
end

function DarkerGray()
    Palette(64, 64, 64)
end

function DarkGray()
    Palette(80, 80, 80)
end

function LightGray()
    Palette(100, 100, 100)
end

function Gray()
    Palette(0xAAAAAAFF)
end

function Silver()
    Palette(0xDDDDDDFF)
end

# ----------------------------------------------------
# Yellow hues
# ----------------------------------------------------
function Orange()
    Palette(255, 127, 0)
end
function SoftOrange()
    Palette(0xFF851BFF)
end
function Yellow()
    Palette(0xFFDC00FF)
end

function Peach()
    Palette(0xF1C6A7FF)
end

# ----------------------------------------------------
# Green hues
# ----------------------------------------------------
function SoftGreen()
    Palette(0x2ECC40FF)
end
function Olive()
    Palette(0x3D9970FF)
end
function Teal()
    Palette(0x39CCCCFF)
end
function Lime()
    Palette(0x01FF70FF)
end

# ----------------------------------------------------
# Blue hues
# ----------------------------------------------------
function SoftBlue()
    Palette(0x0074D9FF)
end
function Navy()
    Palette(0x001f3fFF)
end
function Aqua()
    Palette(0x7FDBFFFF)
end
function LightPurple()
    Palette(0xaaaaffFF)
end

# ----------------------------------------------------
# Pantones
# https://www.pantone.com/color-finder#/pick?pantoneBook=pantoneSolidCoatedV3M2
# ----------------------------------------------------
function LightNavyBlue()
    Palette(0x85B3D1FF)
end

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:
function ueberpruefeEindeutigkeit2x2(f1,f2,f3,f4,B,m_neu)
mneu = m_neu';
f = [f1; f2; f3; f4];
fhorzvert = ones(size(B,2),1);
mhorzvert = 3*ones(size(B,1),1);
mMatrix = B*f;
e1 = [1;0;0;0];
e2 = [0;1;0;0];
e3 = [0;0;1;0];
e4 = [0;0;0;1];

display0 = 1;
display1 = 1;
display2 = 1;
display3 = 1;
display4 = 1;
display5 = 1;

for Zeile in 1:size(B,1)
    if sum(B[Zeile,:]) > sqrt(2)*2 || sum(B[Zeile,:]) <= 1 #Strahlen kürzer als 1cm oder länger als 3* Wurzel(2)
        if display0 == 1
            println("Schaut euch erneut an durch welche Pixel eure Strahlen wie lange laufen und überprüft,")
            println("ob ihr die Längen den richtigen Matrixeinträgen zugeordent habt.")
            display0 = display0+1;
        end
    else
        counter = 1;
        Nulleintraege = 1;
        for Spalte in 1:size(B,2)
            if B[Zeile,Spalte] == 0 #für nur einen Eintrag, dann 1,3,7,9
                Nulleintraege = Nulleintraege+1;
            end
            if B[Zeile,Spalte] > 0 || B[Zeile,Spalte] < 0 #für max 5 Einträge pro Zeile
                counter = counter+1;
            end
        end
        if counter >= 4 #Mehr als 5 Einträge pro Zeile
            if display1 == 1 && display0 == 1
                println("Schaut euch erneut an durch welche Pixel eure Strahlen wie lange laufen und überprüft,")
                println("ob ihr die Längen den richtigen Matrixeinträgen zugeordent habt.")
                display1 = display1+1;
            end
        else
            if B*fhorzvert == mhorzvert #keine neuen Strahlen verwendet
                if display2 == 1 && display1 == 1 && display0 == 1
                    println("Verwendet neue Strahlen und nicht noch einmal die Horizontalen oder Vertikalen, da so keine zusätlichen")
                    println("Informationen weiter gegeben werden.")
                    display2 = display2+1;
                end
            else
                if Nulleintraege-1 >= size(B,2)-1 #nur ein Eintrag falls 1,3,7,9
                    if sum(B[Zeile,:]) > B[Zeile,:]*e1 || sum(B[Zeile,:]) < B[Zeile,:]*e1 || sum(B[Zeile,:]) > B[Zeile,:]*e2 || sum(B[Zeile,:]) < B[Zeile,:]*e2 || sum(B[Zeile,:]) > B[Zeile,:]*e3 || sum(B[Zeile,:]) < B[Zeile,:]*e3 || sum(B[Zeile,:]) > B[Zeile,:]*e4 || sum(B[Zeile,:]) < B[Zeile,:]*e4 && display4 == 1
                        println("Schaut euch erneut an durch welche Pixel eure Strahlen wie lange laufen und überprüft,")
                        println("ob ihr die Pixel korret nummeriert und die richtigen Werte in die Matrix eingetragen habt.")
                        display4 = display4+1;
                    end
                end
            end
        end
    end
end

Akorrekt = B;

mkorrekt = m_neu;

# homogene Lösungen berechnen
homogenkorrekt = nullspace(Akorrekt);

# inhomogene Lösung berechnen
inhomogenkorrekt = pinv(Akorrekt)*mkorrekt';

homogen = nullspace(B);
inhomogen = pinv(B)*m_neu';
 u1 = [0; 0; 0; 0];
p = inhomogen;

x = [f1 f2 f3 f4]*B[5,:];

if display4 == 1 && display3 == 1 && display2 == 1 && display1 == 1 && display0 == 1 && mMatrix[1] == mneu[1] && mMatrix[2] == mneu[2] &&  mMatrix[3] == mneu[3] && mMatrix[4] == mneu[4] && almostequal(x[1,1], m_neu[5])
    println("Eure Idee die Anzahl der Strahlen zu vergrößern ist korrekt! Ihr erhaltet Folgendes:");
    q = round.(p;digits=2);
    print("Für f1 lautet die Lösung: ", q[1],"\n")
    print("Für f2 lautet die Lösung: ", q[2],"\n")
    print("Für f3 lautet die Lösung: ", q[3],"\n")
    print("Für f4 lautet die Lösung: ", q[4],"\n\n")
    #show(stdout,"text/plain",round.(p;digits=2))
    #show(stdout,"text/plain",round.(u1;digits=2))
    #println("p = $p")
    #println("u1 = $u1")
    #if B[5,:] == [0, sqrt(2), sqrt(2), 0]
    #    println("****Macht weiter mit Aufgabenblatt 1a.****")
    #else
    #    println("****Macht weiter mit Aufgabenblatt 1b.****")
    #end
    println("Für die nächste Aufgabe benötigt ihr dieses Bild:\n")
    Daten = load("../Computertomographie_Bilder/ObjektvieleStrahlen2.png");

elseif display4 == 1 && display3 == 1 && display2 == 1 && display1 == 1 && display0 == 1 && mMatrix[1] == mneu[1]
    println("Überprüft eure Berechnung von m erneut. Beachtet, dass auch hier die Länge der Pixel berücksichtigt werden müssen.")
    Daten = [];
end
return Daten;
end


text:
using VectorSphericalWaves
using WignerD

wignerd_functions_list = [
    WignerD.wignerdjmn,
    VectorSphericalWaves.wignerdjmn_ELZOUKA,
    VectorSphericalWaves.wignerdjmn_recurrence,
    VectorSphericalWaves.wignerdjmn_recurrence_memoize,
    VectorSphericalWaves.∂wignerdjmn_by_∂θ
]

# checking the limit of s,m,n for wignerd
th = (LinRange(0, pi, Int(10)));

VectorSphericalWaves.∂wignerdjmn_by_∂θ(7,-7,7,th[2])

s_max = 15
for f in wignerd_functions_list
    println()
    println("==============================================================================")
    println("Now testing <<$f>>")
    @time for s = 1:s_max
        for n = 1:s
            for m = -n:n
                try
                    wigner = f.(s, m, n, th);
                catch e
                    println("using $f, here is the arguments that caused error: s=$s, m=$m, n=$n")
                    bt = backtrace()
                    msg = sprint(showerror, e, bt)
                    println(msg)
                    #println(e)
                    println()
                    break
                end

            end
        end
    end
end

"""
sample results
==============================================================================
Now testing <<wignerdjmn_recurrence>>
  0.936349 seconds (6.94 M allocations: 212.850 MiB, 14.74% gc time)
==============================================================================
Now testing <<wignerdjmn_ELZOUKA>>
  1.108147 seconds (10.80 M allocations: 272.367 MiB, 30.52% gc time)
==============================================================================
Now testing <<wignerdjmn_recurrence_memoize>>
  0.050361 seconds (76.13 k allocations: 4.080 MiB)
"""
println()
println()
println("comparing wignerd_functions:")
println(wignerd_functions_list)
# TODO: WignerD.wignerdjmn is different from the functions we have in src/utils.jl near poles (θ=0 and θ=π)
# cross validation of Wignerd
for s = 1:s_max
    for n = 1:s
        for m = -n:n
            try
                error_12 = filter(!isnan, abs.((wignerd_functions_list[1].(s, m, n, th) - wignerd_functions_list[2].(s, m, n, th)) ./ wignerd_functions_list[1].(s, m, n, th)));
                error_13 = filter(!isnan, abs.((wignerd_functions_list[1].(s, m, n, th) - wignerd_functions_list[3].(s, m, n, th)) ./ wignerd_functions_list[1].(s, m, n, th)));
                error_14 = filter(!isnan, abs.((wignerd_functions_list[1].(s, m, n, th) - wignerd_functions_list[4].(s, m, n, th)) ./ wignerd_functions_list[1].(s, m, n, th)));
                println("s=$s, n=$n, m=$m: error_12_max = $(max(error_12...)), error_13_max = $(max(error_13...)), error_14_max = $(max(error_14...))")
            catch
                println("using 00, here is the arguments that caused error: s=$s, n=$n, m=$m")
                println()
                break
            end

        end
    end
end


text:
#
# Circular Restricted Three-body Problem models
#

"""
A `ModelingToolkit.ODESystem` for the Circular Restricted Three-body Problem.

The order of the states follows: `[x, y, z, ẋ, ẏ, ż]`.

The order of the parameters follows: `[μ]`.

# Extended Help
The Circular Restricted Three-body Problem is a simplified dynamical model
describing one small body (spacecraft, etc.) and two celestial
bodies moving in a circle about their common center of mass.
This may seem like an arbitrary simplification, but this assumption
holds reasonably well for the Earth-Moon, Sun-Earth, and many other
systems in our solar system.

### Usage

```julia
model = CR3BP(; stm=true)
```
"""
@memoize function CR3BP(; stm=false, structural_simplify=true, name=:CR3BP)

    @parameters t μ
    @variables x(t) y(t) z(t) ẋ(t) ẏ(t) ż(t)
    δ = Differential(t)
    r = @SVector [x,y,z]
    v = @SVector [ẋ,ẏ,ż]

    eqs = vcat(
        δ.(r) .~ v,
        δ(ẋ)   ~ x + 2ẏ - (μ*(x + μ - 1)*(sqrt(y^2 + z^2 + (x + μ - 1)^2)^-3)) - ((x + μ)*(sqrt(y^2 + z^2 + (x + μ)^2)^-3)*(1 - μ)),
        δ(ẏ)   ~ y - (2ẋ) - (y*(μ*(sqrt(y^2 + z^2 + (x + μ - 1)^2)^-3) + (sqrt(y^2 + z^2 + (x + μ)^2)^-3)*(1 - μ))),
        δ(ż)   ~ z*(-μ*(sqrt(y^2 + z^2 + (x + μ - 1)^2)^-3) - ((sqrt(y^2 + z^2 + (x + μ)^2)^-3)*(1 - μ)))
    )

    if stm
        @variables Φ[1:6,1:6](t)
        Φ = Symbolics.scalarize(Φ)
        A = Symbolics.jacobian(map(el -> el.rhs, eqs), vcat(r,v))

        LHS = map(δ, Φ)
        RHS = map(simplify, A * Φ)

        eqs = vcat(eqs, [LHS[i] ~ RHS[i] for i in 1:length(LHS)])
    end

    if string(name) == "CR3BP" && stm
        modelname = Symbol("CR3BPWithSTM")
    else
        modelname = name
    end
    sys = ODESystem(
        eqs, t, stm  ? vcat(r,v,Φ...) : vcat(r,v), [μ];
        name = modelname
    )

    return structural_simplify ? ModelingToolkit.structural_simplify(sys) : sys
end

"""
Returns an `ODEFunction` for CR3BP dynamics.
Results are cached with `Memoize.jl`.

The order of the states follows: `[x, y, z, ẋ, ẏ, ż]`.

The order of the parameters follows: `[μ]`.

# Extended Help

### Usage

The `stm`, `structural_simplify`, and `name` keyword arguments
are passed to `CR3BP`. All other keyword arguments are passed
directly to `SciMLBase.ODEFunction`.

```julia
f = CR3BPFunction(; stm=false, jac=true)
let u = randn(6), p = randn(1), t = 0
    f(u, p, t)
end
```
"""
@memoize function CR3BPFunction(; stm=false, structural_simplify=true, name=:CR3BP, kwargs...)
    defaults = (; jac=true)
    options  = merge(defaults, kwargs)
    return ODEFunction(
        CR3BP(; stm=stm, structural_simplify=structural_simplify, name=name);
        options...
    )
end


text:
using Test
using Arpack
using Rimu
using LinearAlgebra
using NewRimuOperators

include("rimucompat.jl")
include("terms.jl")

@testset "Potentials" begin
end


text:
fname = :TRMA57
c = Combi(hessian_sparse,PDataMA57,solve_modelTRDiag,preprocessMA57,decreaseFact,Tparam())
include("Template.jl")


text:
using Test
using ClimateMachine.TicToc

function foo()
    @tic foo
    sleep(0.25)
    @toc foo
end

function bar()
    @tic bar
    sleep(1)
    @toc bar
end

if TicToc.tictoc_enabled
    @testset "TicToc" begin
        @test tictoc() >= 2
        foo_i = findfirst(s -> s == :tictoc__foo, TicToc.timing_info_names)
        bar_i = findfirst(s -> s == :tictoc__bar, TicToc.timing_info_names)
        @test foo_i != nothing
        @test bar_i != nothing
        foo()
        foo()
        @test TicToc.timing_infos[foo_i].ncalls == 2
        @test TicToc.timing_infos[bar_i].ncalls == 0
        @test TicToc.timing_infos[foo_i].time >= 5e8
        bar()
        @test TicToc.timing_infos[bar_i].ncalls == 1
        @test TicToc.timing_infos[bar_i].time >= 1e9
        buf = IOBuffer()
        old_stdout = stdout
        try
            rd, = redirect_stdout()
            TicToc.print_timing_info()
            Libc.flush_cstdio()
            flush(stdout)
            write(buf, readavailable(rd))
        finally
            redirect_stdout(old_stdout)
        end
        str = String(take!(buf))
        @test findnext("tictoc__foo", str, 1) != nothing
        @test findnext("tictoc__bar", str, 1) != nothing
    end
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
include("problem.jl")
include("options.jl")

include("debug_options.jl")
include("record_options.jl")

include("cost_plan.jl")
include("gradient_plan.jl")
include("conjugate_gradient_plan.jl")
include("quasi_newton_plan.jl")
include("hessian_plan.jl")
include("proximal_plan.jl")
include("subgradient_plan.jl")
include("particle_swarm_plan.jl")
include("primal_dual_plan.jl")
include("stochastic_gradient_plan.jl")

include("stopping_criterion.jl")
include("stepsize.jl")


text:
#---------------------------------------------------------
# # [Sound overview](@id 01-overview)
#---------------------------------------------------------

#=
This page illustrates the Julia package
[`Sound`](https://github.com/JeffFessler/Sound.jl).

This page was generated from a single Julia file:
[01-overview.jl](@__REPO_ROOT_URL__/01-overview.jl).
=#

#md # In any such Julia documentation,
#md # you can access the source code
#md # using the "Edit on GitHub" link in the top right.

#md # The corresponding notebook can be viewed in
#md # [nbviewer](http://nbviewer.jupyter.org/) here:
#md # [`01-overview.ipynb`](@__NBVIEWER_ROOT_URL__/01-overview.ipynb),
#md # and opened in [binder](https://mybinder.org/) here:
#md # [`01-overview.ipynb`](@__BINDER_ROOT_URL__/01-overview.ipynb).

# ### Setup

# Packages needed here.

using Sound
using SampledSignals: SampleBuf
using InteractiveUtils: versioninfo

# ### Overview

# The primary purpose of this package, at least initially,
# was to provide a simple way to hear audio signals in Julia.

S = 8192 # sampling rate in Hz
x = 0.7*cos.(2π*(1:S÷2)*440/S)
y = 0.8*sin.(2π*(1:S÷2)*660/S)
sound(x, S) # monophonic

# stereo:
sound([x y], S)

# scale to unit amplitude:
soundsc([x y], S)

# Using `SampleBuf` may provide some convenience.

S = 8192 # sampling rate in Hz
x = 0.7*cos.(2π*(1:S÷2)*440/S)
y = 0.8*sin.(2π*(1:S÷2)*660/S)
sb = SampleBuf([x y], S) # stereo data
sound(sb)

# ### Reproducibility

# This page was generated with the following version of Julia:

io = IOBuffer(); versioninfo(io); split(String(take!(io)), '\n')

# And with the following package versions

import Pkg; Pkg.status()


text:
module TestConstraintQuadToSOC

import LinearAlgebra
import SparseArrays
using Test

using MathOptInterface
const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

include("../utilities.jl")

function test_error_for_nonconvex_quadratic_constraints()
    mock = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    bridged_mock = MOI.Bridges.Constraint.QuadtoSOC{Float64}(mock)
    x = MOI.add_variable(bridged_mock)
    @test_throws(
        MOI.UnsupportedConstraint,
        MOI.add_constraint(
            bridged_mock,
            MOI.ScalarQuadraticFunction(
                [MOI.ScalarQuadraticTerm(1.0, x, x)],
                MOI.ScalarAffineTerm{Float64}[],
                0.0,
            ),
            MOI.GreaterThan(0.0),
        )
    )
    @test_throws(
        MOI.UnsupportedConstraint,
        MOI.add_constraint(
            bridged_mock,
            MOI.ScalarQuadraticFunction(
                [MOI.ScalarQuadraticTerm(-1.0, x, x)],
                MOI.ScalarAffineTerm{Float64}[],
                0.0,
            ),
            MOI.LessThan(0.0),
        )
    )
    return
end

function test_quadratic_constraints_with_2_variables()
    mock = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.QuadtoSOC{Float64}(mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [0.5, 0.5]),
        ),
    )
    MOI.Test.test_constraint_qcp_duplicate_diagonal(bridged_mock, config)
    MOI.empty!(bridged_mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [0.5, (√13 - 1) / 4]),
        ),
    )
    MOI.Test.test_constraint_qcp_duplicate_off_diagonal(bridged_mock, config)
    ci = first(
        MOI.get(
            mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.RotatedSecondOrderCone,
            }(),
        ),
    )
    x, y = MOI.get(mock, MOI.ListOfVariableIndices())
    # The matrix is
    # 2 1
    # 1 2
    # for which the cholesky factorization is U' * U with U =
    # √2 √2/2
    #  . √3/√2
    expected = MOI.VectorAffineFunction{Float64}(
        MOI.VectorAffineTerm.(
            [3, 3, 4],
            MOI.ScalarAffineTerm.([√2, √2 / 2, √3 / √2], [x, y, y]),
        ),
        [1.0, 1.0, 0.0, 0.0],
    )
    @test MOI.get(mock, MOI.ConstraintFunction(), ci) ≈ expected
end

function test_qcp()
    mock = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.QuadtoSOC{Float64}(mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [1 / 2, 7 / 4],
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) =>
                [zeros(2)],
            (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone) => [[0.25, 1.0, -1 / √2]],
        ),
    )
    MOI.Test.test_quadratic_constraint_integration(bridged_mock, config)
    MOI.empty!(bridged_mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [√2],
            (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone) => [[1 / √2, 1 / (2 * √2), -1 / √2]],
        ),
    )
    MOI.Test.test_quadratic_constraint_basic(bridged_mock, config)
    MOI.empty!(bridged_mock)
    MOI.Test.test_quadratic_constraint_minimize(bridged_mock, config)
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.ScalarQuadraticFunction{Float64},
                MOI.LessThan{Float64},
            }(),
        ),
    )
    _test_delete_bridge(
        bridged_mock,
        ci,
        1,
        ((MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone, 0),),
    )
    MOI.empty!(bridged_mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [1.0, 1.0],
            (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone) => [[1.0, 1 / 3, -1 / √2, -1 / √6]],
        ),
    )
    MOI.Test.test_quadratic_constraint_LessThan(bridged_mock, config)
    MOI.empty!(bridged_mock)
    MOI.Test.test_quadratic_constraint_GreaterThan(bridged_mock, config)
    return
end

function test_fill_reducing_permutation()
    model = MOI.Utilities.Model{Float64}()
    bridge = MOI.Bridges.Constraint.QuadtoSOC{Float64}(model)
    x = MOI.add_variables(bridge, 3)
    Q = Float64[2 1 1; 1 2 0; 1 0 2]
    f = 0.5 * x' * Q * x
    MOI.add_constraint(bridge, f, MOI.LessThan(2.0))
    indices = MOI.get(
        model,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.RotatedSecondOrderCone,
        }(),
    )
    F = LinearAlgebra.cholesky(Q)
    # Test that the sparse cholesky pivot is permuted
    SF = LinearAlgebra.cholesky(SparseArrays.sparse(Q))
    @test SF.p == [3, 2, 1]
    U = Matrix(F.U)
    @test MOI.get(model, MOI.ConstraintFunction(), indices[1]) ≈
          MOI.VectorAffineFunction(
        [
            MOI.VectorAffineTerm(3, MOI.ScalarAffineTerm(U[1, 1], x[1])),
            MOI.VectorAffineTerm(3, MOI.ScalarAffineTerm(U[1, 2], x[2])),
            MOI.VectorAffineTerm(3, MOI.ScalarAffineTerm(U[1, 3], x[3])),
            MOI.VectorAffineTerm(4, MOI.ScalarAffineTerm(U[2, 2], x[2])),
            MOI.VectorAffineTerm(4, MOI.ScalarAffineTerm(U[2, 3], x[3])),
            MOI.VectorAffineTerm(5, MOI.ScalarAffineTerm(U[3, 3], x[3])),
        ],
        [1.0, 2.0, 0.0, 0.0, 0.0],
    )
    return
end

end  # module

TestConstraintQuadToSOC.runtests()


text:
function Vector3DFormat()
    return Vector3DFormat(())
end

function Vector3DFormat(arg0::JString, arg1::JString, arg2::JString)
    return Vector3DFormat((JString, JString, JString), arg0, arg1, arg2)
end

function Vector3DFormat(arg0::JString, arg1::JString, arg2::JString, arg3::NumberFormat)
    return Vector3DFormat((JString, JString, JString, NumberFormat), arg0, arg1, arg2, arg3)
end

function Vector3DFormat(arg0::NumberFormat)
    return Vector3DFormat((NumberFormat,), arg0)
end

function format(obj::Vector3DFormat, arg0::JVector, arg1::StringBuffer, arg2::FieldPosition)
    return jcall(obj, "format", StringBuffer, (JVector, StringBuffer, FieldPosition), arg0, arg1, arg2)
end

function format(obj::VectorFormat, arg0::JVector)
    return jcall(obj, "format", JString, (JVector,), arg0)
end

function get_available_locales(::Type{VectorFormat})
    return jcall(VectorFormat, "getAvailableLocales", Vector{Locale}, ())
end

function get_format(obj::VectorFormat)
    return jcall(obj, "getFormat", NumberFormat, ())
end

function get_prefix(obj::VectorFormat)
    return jcall(obj, "getPrefix", JString, ())
end

function get_separator(obj::VectorFormat)
    return jcall(obj, "getSeparator", JString, ())
end

function get_suffix(obj::VectorFormat)
    return jcall(obj, "getSuffix", JString, ())
end

function get_vector3d_format(::Type{Vector3DFormat})
    return jcall(Vector3DFormat, "getVector3DFormat", Vector3DFormat, ())
end

function get_vector3d_format(::Type{Vector3DFormat}, arg0::Locale)
    return jcall(Vector3DFormat, "getVector3DFormat", Vector3DFormat, (Locale,), arg0)
end

function parse(obj::Vector3DFormat, arg0::JString)
    return jcall(obj, "parse", Vector3D, (JString,), arg0)
end

function parse(obj::Vector3DFormat, arg0::JString, arg1::ParsePosition)
    return jcall(obj, "parse", Vector3D, (JString, ParsePosition), arg0, arg1)
end


text:
#---------------------------------------------------------------------#
#This function warps a 2D grid
#Written by F.X. Giraldo on 5/2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function rotate_grid!(coord,Npoin,DFloat)

    α=π/4
    for i=1:Npoin
        x, y = coord[1,i], coord[2,i]
        coord[1,i]=cos(α)*x - sin(α)*y
        coord[2,i]=sin(α)*x + cos(α)*y
    end

end


text:
export get_invite,
    delete_invite

"""
    get_invite(c::Client, invite::AbstractString; kwargs...} -> Invite

Get an [`Invite`](@ref) to a [`Guild`](@ref).
More details [here](https://discordapp.com/developers/docs/resources/invite#get-invite).
"""
get_invite(c::Client, invite::AbstractString; kwargs...) = Response{Invite}(c, :GET, "/invites/$invite"; kwargs...)

"""
    delete_invite(c::Client, invite::AbstractString) -> Invite

Delete an [`Invite`](@ref) to a [`Guild`](@ref).
"""
delete_invite(c::Client, invite::AbstractString) = Response{Invite}(c, :DELETE, "/invites/$invite")


text:
# Obsoleted samplers

# These samplers are not actually used by anyone, and they have not been properly tested.
# The codes are still kept here in case we want to pick them up in future.
#

##### Draw Table #####

# Store an alias table
immutable DiscreteDistributionTable <: Sampler{Univariate,Discrete}
    table::Vector{Vector{Int64}}
    bounds::Vector{Int64}
end

# TODO: Test if bit operations can speed up Base64 mod's and fld's
function DiscreteDistributionTable{T <: Real}(probs::Vector{T})
    # Cache the cardinality of the outcome set
    n = length(probs)

    # Convert all Float64's into integers
    vals = Array(Int64, n)
    for i in 1:n
        vals[i] = round(Int, probs[i] * 64^9)
    end

    # Allocate digit table and digit sums as table bounds
    table = Array(Vector{Int64}, 9)
    bounds = zeros(Int64, 9)

    # Special case for deterministic distributions
    for i in 1:n
        if vals[i] == 64^9
            table[1] = Array(Int64, 64)
            for j in 1:64
                table[1][j] = i
            end
            bounds[1] = 64^9
            for j in 2:9
                table[j] = Array(Int64, 0)
                bounds[j] = 64^9
            end
            return DiscreteDistributionTable(table, bounds)
        end
    end

    # Fill tables
    multiplier = 1
    for index in 9:-1:1
        counts = Array(Int64, 0)
        for i in 1:n
            digit = mod(vals[i], 64)
            # vals[i] = fld(vals[i], 64)
            vals[i] >>= 6
            bounds[index] += digit
            for itr in 1:digit
                push!(counts, i)
            end
        end
        bounds[index] *= multiplier
        table[index] = counts
        # multiplier *= 64
        multiplier <<= 6
    end

    # Make bounds cumulative
    bounds = cumsum(bounds)

    return DiscreteDistributionTable(table, bounds)
end

function rand(table::DiscreteDistributionTable)
    # 64^9 - 1 == 0x003fffffffffffff
    i = rand(1:(64^9 - 1))
    # if i == 64^9
    #   return table.table[9][rand(1:length(table.table[9]))]
    # end
    bound = 1
    while i > table.bounds[bound] && bound < 9
        bound += 1
    end
    if bound > 1
        index = fld(i - table.bounds[bound - 1] - 1, 64^(9 - bound)) + 1
    else
        index = fld(i - 1, 64^(9 - bound)) + 1
    end
    return table.table[bound][index]
end

Base.show(io::IO, table::DiscreteDistributionTable) = @printf io "DiscreteDistributionTable"

##### Huffman Table ######

abstract HuffmanNode{T} <: Sampler{Univariate,Discrete}

immutable HuffmanLeaf{T} <: HuffmanNode{T}
    value::T
    weight::UInt64
end

immutable HuffmanBranch{T} <: HuffmanNode{T}
    left::HuffmanNode{T}
    right::HuffmanNode{T}
    weight::UInt64
end
HuffmanBranch{T}(ha::HuffmanNode{T},hb::HuffmanNode{T}) = HuffmanBranch(ha, hb, ha.weight + hb.weight)

Base.isless{T}(ha::HuffmanNode{T}, hb::HuffmanNode{T}) = isless(ha.weight,hb.weight)
Base.show{T}(io::IO, t::HuffmanNode{T}) = show(io,typeof(t))

function Base.getindex{T}(h::HuffmanBranch{T},u::UInt64)
    while isa(h,HuffmanBranch{T})
        if u < h.left.weight
            h = h.left
        else
            u -= h.left.weight
            h = h.right
        end
    end
    h.value
end

# build the huffman tree
# could be slightly more efficient using a Deque.
function huffman{T}(values::AbstractVector{T},weights::AbstractVector{UInt64})
    leafs = [HuffmanLeaf{T}(values[i],weights[i]) for i = 1:length(weights)]
    sort!(leafs; rev=true)

    branches = Array(HuffmanBranch{T},0)

    while !isempty(leafs) || length(branches) > 1
        left = isempty(branches) || (!isempty(leafs) && first(leafs) < first(branches)) ? pop!(leafs) : pop!(branches)
        right = isempty(branches) || (!isempty(leafs) && first(leafs) < first(branches)) ? pop!(leafs) : pop!(branches)
        unshift!(branches,HuffmanBranch(left,right))
    end

    pop!(branches)
end

function rand{T}(h::HuffmanNode{T})
    w = h.weight
    # generate uniform UInt64 objects on the range 0:(w-1)
    # unfortunately we can't use Range objects, as they don't have sufficient length
    u = rand(UInt64)
    if (w & (w-1)) == 0
        # power of 2
        u = u & (w-1)
    else
        m = typemax(UInt64)
        lim = m - (rem(m,w)+1)
        while u > lim
            u = rand(UInt64)
        end
        u = rem(u,w)
    end
    h[u]
end


text:
import CLIMAParameters

@testset "Override defaults" begin
    struct EarthParameterSet <: AbstractEarthParameterSet end

    ps = EarthParameterSet()
    CLIMAParameters.Planet.grav(::EarthParameterSet) = 2.0
    @test grav(ps) ≈ 2.0

end

max_stars_repo_name: climate-machine/CLIMAParameters

max_issues_repo_name: climate-machine/CLIMAParameters

max_forks_repo_name: climate-machine/CLIMAParameters


text:
#=
trie:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-27
=#

module MyTrie

    # trie node class
    mutable struct TrieNode
        children
        is_end
    end

    const alphabet_num = 26

    function get_node()
        return TrieNode([TrieNode([], false) for i in 1:alphabet_num], false)
    end

    function char_to_idx(ch)
        if Int(ch) - Int('a') == 0
            return 1
        else
            return Int(ch) - Int('a')
        end
    end

    function insert(key)
        # if not present, inserts key into trie
        # if the key is prefix of trie node,
        # just marks leaf node
        p_crawl = root
        len_key = length(key)
        for level in 1:len_key
            idx = char_to_idx(key[level])
            if isempty(p_crawl.children[idx].children) == true
                p_crawl.children[idx] = get_node()
            end
            p_crawl = p_crawl.children[idx]
        end
        p_crawl.is_end = true
    end

    function search(key)
        p_crawl = root
        len_key = length(key)
        for level in 1:len_key
            idx = char_to_idx(key[level])
            if isempty(p_crawl.children[idx].children) == true
                return false
            end
            p_crawl = p_crawl.children[idx]
        end
        return isempty(p_crawl.children) == false && p_crawl.is_end == true
    end

    function show_result(flag)
        if flag == true
            return "Present"
        else
            return "Not Present"
        end
    end

    function main()
        # input keys
        # use only "a" through "z" and lower case
        keys = ["the", "a", "there", "anaswe", "any", "by", "their"]

        println("Input Keys: $(keys)")

        # init root
        global root = get_node()

        # construct trie
        for key in keys
            insert(key)
        end

        # search different keys
        println("Search Key: the = $(show_result(search("the")))")
        println("Search Key: these = $(show_result(search("these")))")
        println("Search Key: their = $(show_result(search("their")))")
        println("Search Key: thaw = $(show_result(search("thaw")))")
        println("Search Key: answer = $(show_result(search("answer")))")
        println("Search Key: by = $(show_result(search("by")))")
    end

end

if abspath(PROGRAM_FILE) == @__FILE__
    using .MyTrie
    MyTrie.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
################################################################################
#                                                                              #
#                        DIRCOL VARIABLES STRUCTURE                            #
#                                                                              #
################################################################################

struct DircolVars{T}
    Z::Vector{T}
    X::SubArray{T}
    U::SubArray{T}
end

function DircolVars(Z::Vector,n::Int,m::Int,N::Int)
    z = reshape(Z,n+m,N)
    X = view(z,1:n,:)
    U = view(z,n+1:n+m,:)
    DircolVars(Z,X,U)
end

function DircolVars(n::Int,m::Int,N::Int)
    Z = zeros((n+m)N)
    z = reshape(Z,n+m,N)
    X = view(z,1:n,:)
    U = view(z,n+1:n+m,:)
    DircolVars(Z,X,U)
end

function DircolVars(X::Matrix,U::Matrix)
    Z = packZ(X,U)
    n,m,N = get_sizes(X,U)
    DircolVars(Z,n,m,N)
end

function DircolVars(res::SolverIterResults)
    DircolVars(to_array(res.X), [to_array(res.U) res.U[end]])
end

DiffFloat = Union{Float64,ForwardDiff.Dual}
struct DircolResults <: SolverIterResults
    vars::DircolVars
    # vars_::DircolVars
    Z::Vector{Float64}
    X::SubArray{Float64}
    U::SubArray{Float64}
    fVal::Matrix{Float64}
    X_::Union{Matrix{Float64},SubArray{Float64}}
    U_::Union{Matrix{Float64},SubArray{Float64}}
    fVal_::Matrix{Float64}
    A::Array{Float64,3}
    B::Array{Float64,3}
    weights::Vector{Float64}
    # c_colloc::Matrix{Float64}
    # ceq::Matrix{Float64}
    # c::Matrix{Float64}
    # DircolResults(Z,X,U,fVal,X_,U_,fVal_,A,B,weights) =
    #     new(Z,X,U,fVal,X_,U_,fVal_,A,B,weights)
end
#
# abstract type DircolMethod end
#
# type HermiteSimpson <: DircolMethod end

function DircolResults(solver::Solver,method::Symbol)
    DircolResults(get_sizes(solver)...,method)
end

function DircolResults(n::Int,m::Int,N::Int,method::Symbol)
    N,N_ = get_N(N,method)
    Z = zeros(N*(n+m))
    vars = DircolVars(Z,n,m,N)
    X,U = vars.X,vars.U
    fVal = zero(X)
    if method == :hermite_simpson
        X_ = zeros(n,N_)
        U_ = zeros(m,N_)
        fVal_ = zeros(n,N_)
        A = zeros(n,n,N_)
        B = zeros(n,m,N_)
    elseif method == :midpoint
        X_ = zeros(n,N) # midpoints plus terminal
        U_ = U
        fVal_ = zero(X_)
        N_ = size(X_,2)
        A = zeros(n,n,N_)
        B = zeros(n,m,N_)
    else
        X_ = X
        U_ = U
        fVal_ = fVal
        N_ = size(X_,2)
        A = zeros(n,n+m,N_) # These methods conveniently use the gradient of Z
        B = zeros(0,0,N_)
    end

    weights = get_weights(method,N_)
    DircolResults(vars,Z,X,U,fVal,X_,U_,fVal_,A,B,weights)
end


text:
function sample_time(
    rng::AbstractRNG, proj::Project, overall_start=:start, overall_finish=:finish
)
    realizations = Dict{Symbol, Float64}()
    realization(task_name::Symbol) = get!(realizations, task_name) do
        task_dist = proj.task_durations[task_name]
        if task_dist isa Real  # it is some constant
            return task_dist
        else  # it is a collection or a Distribution
            rand(rng, task_dist)
        end
    end

    function time_taken(start, finish)
        start == finish && return 0.0
        dependencies = (pre for (pre, this) in proj.links if this==finish)
        time_to_begin = maximum(time_taken.(start, dependencies))
        time_to_finish = realization(finish)
        return time_to_begin + time_to_finish
    end

    return time_taken(overall_start, overall_finish)
end

Base.eltype(::Type{<:Project}) = Float64
Base.eltype(::Type{<:Random.SamplerTrivial{<:Project}}) = Float64
function Random.Sampler(rng::AbstractRNG, proj::Project, ::Union{Val{1}, Val{Inf}})
    return Random.SamplerTrivial(proj)
end

Random.rand(rng::AbstractRNG, d::Random.SamplerTrivial{<:Project}) = sample_time(rng, d[])

# For Plots.jl users: should really be called with StatsPlot's `density` function
@recipe function f(proj::Project, nsamples=100_000, rng=Random.default_rng())
    title --> "Completion Time"
    rand(rng, proj, nsamples)
end

"""
    density(proj::Project, nsamples=100_000, rng=Random.default_rng())

Displayed a probability density function for how long it will take to complete the project.
`nsamples` controls how many samples are used to estimate the distribution, increasing this makes it take longer, but be more accurate.
`rng` controls the random number generator to use, you shouldn't normally need to touch this.

[Gadfly.jl](https://github.com/GiovineItalia/Gadfly.jl) is used for the displaying.
A [RecipesBase.jl](https://github.com/JuliaPlots/RecipesBase.jl) recipe is also provided and so this can be used with [StatsPlots.jl](https://github.com/JuliaPlots/StatsPlots.jl) similarly: as `StatsPlots.density(proj)`.
"""
function density(proj::Project, nsamples=100_000, rng=Random.default_rng(); kwargs...)
    # Our own Gadfly based density function
    if !isempty(kwargs)
        @warn "ProjectMangement.density does not support these keyword arguments. You may have intended to use StatsPlot.density." kwargs
    end
    samples = rand(rng, proj, nsamples)
    Gadfly.plot(
        x=samples,
        Gadfly.Geom.density,
        Gadfly.Guide.title("Completion Time PDF"),
        Gadfly.Guide.xlabel("Completed By"),
        Gadfly.Guide.ylabel("Probability Density"),
    )
end


text:
using Combinatorics
using LinearAlgebra
using ArgCheck

function is_a_square_three_tensor(W::Array)

	length(size(W)) == 3 && all(size(W) .== size(W)[1])

end

function naive_tensor(W::Array)

    if is_a_square_three_tensor(W)
        n = size(W)[1]

        sum_tensor_diag(sigma,rho, W) = prod(W[i,sigma[i],rho[i]] for i = 1:n)
        return sum(sum_tensor_diag(sigma,rho, W) for sigma in permutations(collect(1:n)) for rho in permutations(collect(1:n)))
    else
        error("tensor permanent implemented only for square 3-indices tensors")
    end

end

function ryser_tensor(W::Array)

	"""implements the tensor permanent according to Ryser's factorization
	without using Gray factoring

	corresponds to Eq. 21 of https://arxiv.org/abs/1410.7687"""

	@argcheck is_a_square_three_tensor(W) "tensor permanent implemented only for square 3-indices tensors"

	n = size(U)[1]
    nstring = collect(1:n)
    sub_nstring = collect(powerset(nstring))
    sub_nstring = sub_nstring[2:length(sub_nstring)]
    res = 0

    function delta_set(S1, S2)
        if S1 == S2
            return 1
        else
            return 0
        end
    end

    for r = 1:length(sub_nstring)
        for s = r:length(sub_nstring)
            R = sub_nstring[r]
            S = sub_nstring[s]

            t = prod(
                sum(
					W[rr,ss,j] for rr in R for ss in S
                ) for j = 1:n
            )
            res +=
                (2 - delta_set(S, R)) * (-1)^(length(S) + length(R)) * real(t)
        end
    end
    return res

end

function multi_dim_ryser(U, gram_matrix)

    # https://arxiv.org/pdf/1410.7687.pdf

    n = size(U)[1]
    nstring = collect(1:n)
    sub_nstring = collect(powerset(nstring))
    sub_nstring = sub_nstring[2:length(sub_nstring)]
    res = 0

    function delta_set(S1, S2)
        if S1 == S2
            return 1
        else
            return 0
        end
    end

    for r = 1:length(sub_nstring)
        for s = r:length(sub_nstring)
            R = sub_nstring[r]
            S = sub_nstring[s]

            t = prod(
                sum(
                    U[ss, j] * conj(U[rr, j]) * gram_matrix[rr, ss] for rr in R for ss in S
                ) for j = 1:n
            )
            res +=
                (2 - delta_set(S, R)) * (-1)^(length(S) + length(R)) * real(t)
        end
    end
    return res
end


text:
#########################################################################
#
#   Expression to graph conversion
#
#########################################################################

##########  Parameterized type to ease AST exploration  ############
type ExH{H}
    head::Symbol
    args::Vector
    typ::Any
end
toExH(ex::Expr) = ExH{ex.head}(ex.head, ex.args, ex.typ)
toExpr(ex::ExH) = Expr(ex.head, ex.args...)

typealias ExEqual    ExH{:(=)}
typealias ExDColon   ExH{:(::)}
typealias ExColon    ExH{:(:)}
typealias ExPEqual   ExH{:(+=)}
typealias ExMEqual   ExH{:(-=)}
typealias ExTEqual   ExH{:(*=)}
typealias ExTrans    ExH{Symbol("'")}
typealias ExCall     ExH{:call}
typealias ExBlock    ExH{:block}
typealias ExLine     ExH{:line}
typealias ExVcat     ExH{:vcat}
typealias ExVect     ExH{:vect}
typealias ExCell1d   ExH{:cell1d}
typealias ExCell     ExH{:cell1d}
typealias ExFor      ExH{:for}
typealias ExRef      ExH{:ref}
typealias ExIf       ExH{:if}
typealias ExComp     ExH{:comparison}
typealias ExDot      ExH{:.}
typealias ExTuple    ExH{:tuple}
typealias ExReturn   ExH{:return}
typealias ExBody     ExH{:body}
typealias ExQuote    ExH{:QuoteNode}

#  s     : expression to convert
#  svars : vars set since the toplevel graph (helps separate globals / locals)
function tograph(s, evalmod=Main, svars=Any[])

    explore(ex::Any)       = error("[tograph] unmanaged type $ex ($(typeof(ex)))")
    explore(ex::Expr)      = explore(toExH(ex))
    explore(ex::ExH)       = error("[tograph] unmanaged expr type $(ex.head) in ($ex)")

    explore(ex::ExLine)         = nothing     # remove line info
    explore(ex::LineNumberNode) = nothing     # remove line info
    explore(ex::LabelNode)      = nothing     # remove label info
    explore(ex::GotoNode)       = error("[tograph] GotoNode found") # this should be caught before in frdiff.jl
    explore(ex::QuoteNode)      = addnode!(g, NConst(ex.value))  # consider as constant

    explore(ex::ExReturn)  = explore(ex.args[1]) # focus on returned statement

    explore(ex::ExVcat)    = explore(Expr(:call, :vcat, ex.args...) )  # translate to vcat() call, and explore
    explore(ex::ExVect)    = explore(Expr(:call, :vcat, ex.args...) )  # translate to vcat() call, and explore
    explore(ex::ExCell1d)  = explore(Expr(:call, :(Base.cell_1d), ex.args...) )  # translate to cell_1d() call, and explore
    explore(ex::ExTrans)   = explore(Expr(:call, :transpose, ex.args[1]) )  # translate to transpose() and explore
    explore(ex::ExColon)   = explore(Expr(:call, :colon, ex.args...) )  # translate to colon() and explore
    explore(ex::ExTuple)   = explore(Expr(:call, :tuple, ex.args...) )  # translate to tuple() and explore

    explore(ex::ExPEqual)  = (args = ex.args ; explore( Expr(:(=), args[1], Expr(:call, :+, args[1], args[2])) ) )
    explore(ex::ExMEqual)  = (args = ex.args ; explore( Expr(:(=), args[1], Expr(:call, :-, args[1], args[2])) ) )
    explore(ex::ExTEqual)  = (args = ex.args ; explore( Expr(:(=), args[1], Expr(:call, :*, args[1], args[2])) ) )

    explore(ex::Real)      = addnode!(g, NConst(ex))

    explore(ex::ExBlock)   = map( explore, ex.args )[end]
    explore(ex::ExBody)    = map( explore, ex.args )[end]

    explore(ex::ExComp)    = addnode!(g, NComp(ex.args[2], [explore(ex.args[1]), explore(ex.args[3])]))

    # explore(ex::ExDot)     = addnode!(g, NDot(ex.args[2],     [ explore(ex.args[1]) ]))
    explore(ex::ExDot)     = explore(Expr(:call, :getfield, ex.args...))
    explore(ex::ExRef)     = explore(Expr(:call, :getindex, ex.args...))

    function explore(ex::Symbol)
        hassym(g.seti, ex)       && return getnode(g.seti, ex)
        hassym(g.exti, ex)       && return getnode(g.exti, ex)

        nn = addnode!(g, NExt(ex))    # create external node for this var
        g.exti[nn] = ex
        return nn
    end

    function explore(ex::ExCall)
        sf  = ex.args[1]

        # catch comparisons (they are :call since julia 0.5)
        if sf in [:(>), :(<), :(<=), :(>=)]
            return addnode!(g, NComp(ex.args[1], [explore(ex.args[2]), explore(ex.args[3])]))
        end

        # catch getindex, etc. qualified by a module
        sf2 = if isa(sf, Expr) && sf.head == :. && isa(sf.args[2], QuoteNode)
                sf.args[2].value
              else
                nothing
              end

        if sf == :getindex || sf2 == :getindex
            nv = explore(ex.args[2])
            ps = indexspec(nv, ex.args[3:end])
            return addnode!(g, NRef(:getidx, vcat([nv], ps)))

        elseif sf == :setindex! || sf2 == :setindex!
            isa(ex.args[2], Symbol) ||
                error("[tograph] setindex! only allowed on variables, $(ex.args[2]) found")

            nv  = explore(ex.args[2]) # node whose subpart is assigned
            ps  = indexspec(nv, ex.args[4:end])
            rhn = addnode!(g, NSRef(:setidx,
                                    [ nv,                               # var modified in pos #1
                                      explore(ex.args[3]),              # value affected in pos #2
                                      ps...] ))                         # dims

            rhn.precedence = filter(n -> nv in n.parents && n != rhn, g.nodes)
            g.seti[rhn] = ex.args[2]

            return nothing

        elseif sf == :getfield || sf2 == :getfield
            return addnode!(g, NDot(ex.args[3], [ explore(ex.args[2]) ]))

        elseif sf == :setfield! || sf2 == :setfield!
            isa(ex.args[2], Symbol) ||
                error("[tograph] setfield! only allowed on variables, $(ex.args[2]) found")

            nv  = explore(ex.args[2]) # node whose subpart is assigned
            rhn = addnode!(g, NSDot(ex.args[3],
                                    [ nv,                               # var modified in pos #1
                                      explore(ex.args[4])]))            # value affected in pos #2

            rhn.precedence = filter(n -> nv in n.parents && n != rhn, g.nodes)
            g.seti[rhn] = ex.args[2]

            return nothing

        else
            return addnode!(g, NCall(  :call,
                                        map(explore, ex.args[1:end]) ))
        end
    end

    function explore(ex::ExEqual)
        lhs = ex.args[1]

        if isSymbol(lhs)  # x = ....
            lhss = lhs

            # set before ? call explore
            if lhss in union(svars, collect(syms(g.seti)))
                vn = explore(lhss)
                rhn  = addnode!(g, NSRef(:setidx,
                                         [ vn,                     # var modified in pos #1
                                           explore(ex.args[2]) ])) # value affected in pos #2
                rhn.precedence = filter(n -> vn in n.parents && n != rhn, g.nodes)

            else # never set before ? assume it is created here
                rhn = explore(ex.args[2])

                # we test if RHS has already a Symbol
                # if it does, to avoid loosing it, we create an NIn node
                if hasnode(g.seti, rhn)
                    rhn = addnode!(g, NIn(lhss, [rhn]))
                end
            end
            g.seti[rhn] = lhss

        elseif isRef(lhs)   # x[i] = ....
            lhss = lhs.args[1]
            explore( Expr(:call, :setindex!, lhss, ex.args[2], lhs.args[2:end]...) )

        elseif isDot(lhs)   # x.field = ....
            lhss = lhs.args[1]
            explore( Expr(:call, :setfield!, lhss, lhs.args[2], ex.args[2]) )

        else
            error("[tograph] $(toExpr(ex)) not allowed on LHS of assigment")
        end

        return nothing
        # return rhn
    end

    function explore(ex::ExFor)
        is = ex.args[1].args[1]
        isa(is, Symbol) ||
            error("[tograph] for loop using several indexes : $is ")

        # explore the index range
        nir = explore(ex.args[1].args[2])

        # explore the for block as a separate graph
        nsvars = union(svars, collect(syms(g.seti)))
        g2 = tograph(ex.args[2], evalmod, nsvars)

        # create "for" node
        nf = addnode!(g, NFor( Any[ is, g2 ] ))
        nf.parents = [nir]  # first parent is indexing range fo the loop

        # create onodes (node in parent graph) for each exti
        for (k, sym) in g2.exti.kv
            sym==is  && continue # loop index should be excluded
            pn = explore(sym)  # look in setmap, externals or create it
            g2.exto[pn] = sym
            push!(nf.parents, pn) # mark as parent of for loop
        end

        # create onodes and 'Nin' nodes for each seti
        #  will be restricted to variables that are defined in parent
        #   (others are assumed to be local to the loop)
        for (k, sym) in g2.seti.kv
            if sym in nsvars && sym != is # only for variables set in parent scope
                pn = explore(sym)                   # create node if needed
                rn = addnode!(g, NIn(sym, [nf]))    # exit node for this var in this graph
                g.seti[rn] = sym                    # signal we're setting the var
                g2.seto[rn] = sym

                append!(nf.precedence, filter(n -> pn in n.parents && n != nf, g.nodes))

                # create corresponding exti if it's not already done
                if !hassym(g2.exto, sym)
                    g2.exto[pn] = sym
                    push!(nf.parents, pn) # mark as parent of for loop
                end
            end
        end
    end

    #### translates ':' and 'end' special symbols in getindex / setindex!
    function indexspec(nv, as)
        p  = ExNode[]
        for (i,na) in enumerate(as)
            if length(as)==1 # single dimension
                ns = addgraph!(:( length(x) ), g, Dict(:x => nv) )
            else # several dimensions
                ns = addgraph!(:( size(x, $i) ), g, Dict(:x => nv) )
            end

            na==:(:) && (na = Expr(:(:), 1, :end) )  # replace (:) with (1:end)

            # explore the dimension expression
            nsvars = union(svars, collect(syms(g.seti)))
            ng = tograph(na, evalmod, nsvars)

            # find mappings for onodes, including :end
            vmap = Dict()
            for (k, sym) in ng.exti.kv
                vmap[sym] = sym == :end ? ns : explore(sym)
            end

            nd = addgraph!(ng, g, vmap)
            push!(p, nd)
        end
        p
    end

    #  top level graph
    g = ExGraph()

    exitnode = explore(s)
    # exitnode = nothing if only variable assigments in expression
    #          = ExNode of last calc otherwise

    # id is 'nothing' for unnassigned last statement
    exitnode!=nothing && ( g.seti[exitnode] = nothing )

    # Resolve external symbols that are Functions, DataTypes or Modules
    # and turn them into constants
    for en in filter(n -> isa(n, NExt) & !in(n.main, svars) , keys(g.exti))
        if isdefined(evalmod, en.main)  # is it defined
            tv = evalmod.eval(en.main)
            isa(tv, TypeConstructor) && error("[tograph] TypeConstructors not supported: $ex $(tv), use DataTypes")
            if isa(tv, DataType) || isa(tv, Module) || isa(tv, Function)
                delete!(g.exti, en)
                nc = addnode!(g, NConst( tv ))
                fusenodes(g, nc, en)
            end
        end
    end

    g
end


text:
using MPSGE
using MPSGE.JuMP.Containers

m = Model()

goods = [:x, :y]
factors = [:l, :k]

factor = DenseAxisArray(Float64[50 50; 20 30], goods, factors)
supply = DenseAxisArray(Float64[100, 50], goods)

@parameter(m, endow, 1.0)

Y = add!(m, Sector(:Y, indices=(goods,)))
U = add!(m, Sector(:U))

PC = add!(m, Commodity(:PC, indices=(goods,)))
PU = add!(m, Commodity(:PU))
PF = add!(m, Commodity(:PF, indices=(factors,)))

RA = add!(m, Consumer(:RA, benchmark=150.))

for i in goods
    @production(m, Y[i], 1, [Output(PC[i], supply[i])], [Input(PF[:l], factor[i,:l]), Input(PF[:k], factor[i,:k])])
end

# @production(m, [i in goods], Y[i], 1, PC[i], supply[i],  [Input(PF[f], factor[i,f]) for f in factors])

@production(m, U, 1, [Output(PU, 150)], [Input(PC[:x], 100), Input(PC[:y], 50)])

@demand(m, RA, [Demand(PU, 150)], [Endowment(PF[:l], :(70 * $endow)), Endowment(PF[:k], 80.)])

solve!(m, cumulative_iteration_limit=0)

set_fixed!(PC[:x], true)
set_value(endow, 1.1)

solve!(m)

algebraic_version(m)