text:
# vim: set wrap
#: by David Thierry 2021
using JuMP
using Clp
using DataFrames
using CSV

#: Data frames section
#: Load parameters

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

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

#: Load random floats
df_rf = DataFrame(CSV.File("../resources/1year_2c.csv"))

#: Assign parameters
#: Gas parameters

bPowGasTeLoad = df_gas[1, 2]
aPowGasTeLoad = df_gas[1, 3]

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

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

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

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

#: Steam params
bCcRebDutyEload = df_steam_full_steam[2, 2]
aCcRebDutyEload = df_steam_full_steam[2, 3]
#: Full power gives you the min steam
bDacSteaBaseEload = df_steam_full_power[3, 2]
aDacSteaBaseEload = df_steam_full_power[3, 3]

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

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

aLpSteaToPow = 78.60233832  # MMBtu to kwh

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

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

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

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

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

#: Horizon Lenght
tHorz = 1
rnum = df_rf[1:tHorz, 1]
# USD/MWh
# pow_price =(df_pow_c[!, "MiNg_150_ERCOT"])  # USD/MWh

pow_price =(df_pow_c[!, "MiNg_150_PJM-W"])  # USD/MWh
#: Natural gas price
# 0.056 lb/cuft STP
#std_w_ng1000cuft = 0.056 * 1000
#cNgPerLbUsd = (3.5 / 1000) / 0.056
cNgPerMmbtu = 3.5
m = Model()
#
# aPowUseRateComp = 0.279751187  # MWh/tonneCo2
aPowUseRateComp = 0.076 # MWh/tonneCo2 (Trimeric)
#
cCostInvCombTurb = 1e+02
cCostInvSteaTurb = 1e+02
cCostInvTransInter = 1e+02
cCostInvPcc = 1e+02
cCostInvDac = 1e+03
cCostInvComp = 1e+01

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

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

hSlice = 4  # the number of slices of a given hour
# There are tHorz - 1 slices
# Each slice has hSlice points, but only states have the 0th
@variable(m, 60 <= yGasTelecLoad[0:tHorz, 1:hSlice] <= 100)
@variable(m, 0 <= xPowGasTur[0:tHorz, 1:hSlice])
@variable(m, 0 <= xPowGross[0:tHorz, 1:hSlice])
@variable(m, 0 <= xPowOut[0:tHorz, 1:hSlice])

@variable(m, xAuxPowGasT[0:tHorz, 1:hSlice] >= 0)

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

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

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

@variable(m, 0 <= xSideSteam[0:tHorz, 1:hSlice])
@variable(m, 0 <= xSteaPowLp[0:tHorz, 1:hSlice])
@variable(m, 0 <= xSideSteaDac[0:tHorz, 1:hSlice])

#

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

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

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

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

@variable(m, 0 <= xR0Flue[0:tHorz, 0:hSlice])  # Not-a-State

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

@variable(m, 0 <= vAbsFlue[0:tHorz, 0:hSlice])  # Input
@variable(m, 0 <= vRegFlue[0:tHorz, 0:hSlice])  # Input

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

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

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

@variable(m, 0 <= xR0Air[0:tHorz, 0:hSlice])  # Not-a-State

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

@variable(m, 0 <= vAbsAir[0:tHorz, 0:hSlice])  # Input
@variable(m, 0 <= vRegAir[0:tHorz, 0:hSlice])  # Input

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

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

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

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

# Constraints
# Gas Turbine
@constraint(m, powGasTur[i = 0:tHorz, j = 1:hSlice],
            xPowGasTur[i, j] == aPowGasTeLoad * yGasTelecLoad[i, j]
            + bPowGasTeLoad
           )
#
@constraint(m, fuelEq[i = 0:tHorz, j = 1:hSlice],
            xFuel[i, j] == aFuelEload * yGasTelecLoad[i, j]
            + bFuelEload
           )
#
@constraint(m, co2FuelEq[i = 0:tHorz, j = 1:hSlice],
            xCo2Fuel[i, j] == aEmissFactEload * yGasTelecLoad[i, j]
            + bEmissFactEload
           )

@constraint(m, auxPowGasT[i = 0:tHorz, j = 1:hSlice],
            xAuxPowGasT[i, j] == aAuxRateGas * yGasTelecLoad[i, j]
            + bAuxRateGas
           )
#
# Steam
#
@constraint(m, powHpEq[i = 0:tHorz, j = 1:hSlice],
            xPowHp[i, j] == aPowHpEload * yGasTelecLoad[i, j]
            + bPowHpEload
           )
#
@constraint(m, powIpEq[i = 0:tHorz, j = 1:hSlice],
            xPowIp[i, j] == aPowIpEload * yGasTelecLoad[i, j]
            + bPowIpEload
           )

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

@constraint(m, ccRebDutyEq[i = 0:tHorz, j = 1:hSlice],
            xCcRebDuty[i, j] ==
            aCcRebDutyEload * yGasTelecLoad[i, j]
            + bCcRebDutyEload
           )

@constraint(m, dacSteaDutyEq[i = 0:tHorz, j = 1:hSlice],
            xDacSteaBaseDuty[i, j] == aDacSteaBaseEload * yGasTelecLoad[i, j]
            + bDacSteaBaseEload
           )

@constraint(m, sideSteaEloadEq[i = 0:tHorz, j = 1:hSlice],
            xSideSteam[i, j] == aSideSteaEload * yGasTelecLoad[i, j]
            + bSideSteaEload
           )

@constraint(m, sideSteaRatioEq[i = 0:tHorz, j = 1:hSlice],
            xSideSteam[i, j] == xSideSteaDac[i, j] + xSteaPowLp[i, j]
           )

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

@constraint(m, auxPowSteaTEq[i = 0:tHorz, j = 1:hSlice],
            xAuxPowSteaT[i, j] == aAuxRateStea * yGasTelecLoad[i, j]
            + bAuxRateStea
           )

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

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

# DAC-Flue
@constraint(m, xA0FlueEq[i = 0:tHorz, j=0:hSlice],
            xA0Flue[i, j] == vAbsFlue[i, j]
           )
@constraint(m, xR0FlueEq[i = 0:tHorz, j=0:hSlice],
            xR0Flue[i, j] == vRegFlue[i, j]
           )
#
#: "State equation"
@constraint(m, a1dFlueEq[i = 0:tHorz, j=1:hSlice],
            xA1Flue[i, j] == xA0Flue[i, j-1]
           )
#: "State equation"
@constraint(m, a2dFlueEq[i = 0:tHorz, j=1:hSlice],
            xA2Flue[i, j] == xA1Flue[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] - vAbsFlue[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] - vRegFlue[i, j-1] + xA2Flue[i, j-1]
           )
# Initial conditions
@constraint(m, icXa1FlueEq, xA1Flue[0, 0] == 0.)
@constraint(m, icxA2FlueEq, xA2Flue[0, 0] == 0.)
@constraint(m, icAR1FlueEq, xR1Flue[0, 0] == 0.)
#@constraint(m, capDacFlueEq, xFflue[0] == xSorbFreshFlue)
@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, endxA2FlueEq, xA2Flue[tHorz] == 0.)
#@constraint(m, endAR1FlueEq, xR1Flue[tHorz] == 0.)
#@constraint(m, endDacFlueEq, xFflue[tHorz] == xSorbFreshFlue)
#@constraint(m, endDacFlueEq, xFflue[tHorz] == aSorbAmountFreshFlue)
#@constraint(m, endSsFlueEq, xSflue[tHorz] == 0.)

#
#
@constraint(m, co2CapDacFlueEq[i = 0:tHorz, j = 1:hSlice],
            xCo2CapDacFlue[i, j] == aSorbCo2CapFlue * xR1Flue[i, j]
           )
#
@constraint(m, steamUseDacFlueEq[i = 0:tHorz, j = 1:hSlice],
            xSteaUseDacFlue[i, j] == aSteaUseRateDacFlue * xCo2CapDacFlue[i, j]
           )
#
@constraint(m, powUseDacFlueEq[i = 0:tHorz, j = 1:hSlice],
            xPowUseDacFlue[i, j] == aPowUseRateDacFlue * xCo2CapDacFlue[i, j]
           )
# Equal to the amount vented, at least in flue mode.
@constraint(m, co2DacFlueVentEq[i = 0:tHorz, j = 1:hSlice],
            xCo2DacVentFlue[i, j] == xCo2DacFlueIn[i, j] - xCo2CapDacFlue[i, j]
           )

# DAC-Air
# Bluntly assume we can just take CO2 from air in pure form.
@constraint(m, xA0AirEq[i = 0:tHorz, j = 0:hSlice],
            xA0Air[i, j] == vAbsAir[i, j]
           )
@constraint(m, xR0AirEq[i = 0:tHorz, j = 0:hSlice],
            xR0Air[i, j] == vRegAir[i, j]
           )
# "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] - vAbsAir[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] - vRegAir[i, j-1] + xA2Air[i, j-1]
           )

# Initial conditions - Air
#@constraint(m, capDacAirEq, xFair[0] == xSorbFreshAir)
@constraint(m, capDacAirEq, xFair[0, 0] == aSorbAmountFreshAir)
# aSorbAmountFreshFlue
@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[0] == xSorbFreshAir)
#@constraint(m, endDacAirEq, xFair[0, 0] == aSorbAmountFreshAir)
#@constraint(m, endA1AirEq, xA1Air[0, 0] == 0.)
#@constraint(m, endA2AirEq, xA2Air[0, 0] == 0.)
#@constraint(m, endAR1AirEq, xR1Air[0, 0] == 0.)
#@constraint(m, endSsAirEq, xSair[0, 0] == 0.)

#
# Money, baby.
@constraint(m, co2CapDacAirEq[i = 0:tHorz, j=1:hSlice],
            xCo2CapDacAir[i, j] == aSorbCo2CapAir * xR1Air[i, j]
           )
#
@constraint(m, steamUseDacAirEq[i = 0:tHorz, j=1:hSlice],
            xSteaUseDacAir[i, j] == aSteaUseRateDacAir * xCo2CapDacAir[i, j]
           )
#
@constraint(m, powUseDacAirEq[i = 0:tHorz, j=1:hSlice],
            xPowUseDacAir[i, j] == aPowUseRateDacAir * xCo2CapDacAir[i, j]
           )

@constraint(m, dacSteaSlackEq[i = 0:tHorz, j=1:hSlice],
            xDacSteaSlack[i, j] == xDacSteaDuty[i, j]
            - xCo2DacHr[i]
           )

@constraint(m, dacsteahrEq[i = 0:tHorz],
            xCo2DacHr[i] ==
            sum(xSteaUseDacFlue[i, j] + xSteaUseDacAir[i, j]
                for j in 1:hSlice)
           )

# Co2 Compression
#
@constraint(m, co2CompEq[i = 0:tHorz, j = 1:hSlice],
            xCo2Comp[i, j] == xCo2CapPcc[i, j]
            + sum(xCo2CapDacFlue[i, k] + xCo2CapDacAir[i, k] for k in 1:hSlice)
           )
#
@constraint(m, powUseCompEq[i = 0:tHorz, j = 1:hSlice],
            xPowUseComp[i, j] == aPowUseRateComp * xCo2Comp[i, j]
           )
#
# @constraint(m, powUseCompIn[i = 0:tHorz - 1], xPowUseComp[i] <= vCapComp)

# @constraint(m, co2VentEq[i = 0:tHorz - 1],
# xCo2Vent[i] == vCo2PccVent[i] + xCo2DacVentFlue[i])
@constraint(m, co2VentEq[i = 0:tHorz, j = 1:hSlice],
            xCo2Vent[i, j] == vCo2PccVent[i, j]
            + sum(xCo2DacVentFlue[i, k] - xCo2CapDacAir[i, k] for k in 1:hSlice)
           )

## Overall

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

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

# Piece-wise constant DOF
@constraint(m, pwleq[i = 0:tHorz, j = 2:hSlice],
            yGasTelecLoad[i, 1] == yGasTelecLoad[i, j]
           )
@constraint(m, pwssd[i = 0:tHorz, j = 2:hSlice],
           xSideSteaDac[i, 1] == xSideSteaDac[i, j]
           )
@constraint(m, pwco2vent[i = 0:tHorz, j = 2:hSlice],
            vCo2PccVent[i, 1] == vCo2PccVent[i, j]
           )

# 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, conta2flue[i = 1:tHorz], xA2Flue[i, 0] == xA2Flue[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])

@expression(m, eObjfExpr, sum(
                              cNgPerMmbtu * xFuel[i, hSlice]
                              + cEmissionPrice * xCo2Vent[i, hSlice]
                              + cCo2TranspPrice * xCo2Comp[i, hSlice]
                              - pow_price[i+ 1] * xPowOut[i, hSlice]
                              for i in 0:tHorz
                             )
           )

@objective(m, Min, eObjfExpr)
println("The number of variables: ")
println(num_variables(m))
global 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()
println("The number of constraints: ", n)
println()
set_optimizer(m, Clp.Optimizer)
set_optimizer_attribute(m, "LogLevel", 3)
set_optimizer_attribute(m, "PresolveType", 1)

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

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

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

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

# Raw materials.
# xFuel
# sF0
# sS0
#
# penalties
# vCo2Vent

## vCo2PccVent * 2N

# Actual variables
# xPowGasTur * 2

max_stars_repo_name: dthierry/princetonDacLti

max_issues_repo_name: dthierry/princetonDacLti

max_forks_repo_name: dthierry/princetonDacLti


text:
using FATEStools
using Test

@testset "FATEStools.jl" begin

    # Test meta-data add function
    fn_nc = "testfile.nc"
    fn_meta = "testmeta_src.json"
    fn_temp = "testmeta_temp.json"
    ncmeta_add!(fn_nc,fn_meta) # what should this return?

    # Test function to define JSON template with necessary meta-date fields for both global and variable attributes
    ncmeta_template(fn_temp)

    # Test function to check the user-defined JSON file for CF conformance (min requirements?)
    ncmeta_cfcheck(fn_nc) # need to define return states

end

ext: jl


text:
function lyon_passive_ear(signal::AbstractVector; sample_rate = 16000, decimation_factor::Int = 1,
        ear_q = 8, step_factor = ear_q / 32, differ = true, useagc = true, tau_factor = 3)

    ear_filters, _ = design_lyon_filters(sample_rate, ear_q, step_factor)
    n_samples = length(signal)

    n_output_samples = floor(Int, n_samples / decimation_factor)
    n_channels = size(ear_filters, 1) ## ???

    sos_output = zeros(n_channels, decimation_factor)
    sos_state = zeros(n_channels, 2)
    agc_state = zeros(n_channels, 4)
    y = zeros(n_channels, n_output_samples)

    dec_eps = epsilon_from_tau(decimation_factor / sample_rate * tau_factor, sample_rate)

    dec_state = zeros(n_channels, 2)

    _coeffs = [0.0, 0, 1, -2 * (1 - dec_eps), (1 - dec_eps) ^ 2]
    dec_filt = set_gain(_coeffs, 1, 0, sample_rate)

    epses = [epsilon_from_tau(x, sample_rate) for x in [.64, .16, .04, .01]]
    tars = [.0032, .0016, .0008, .0004]

    # move agc_params out of loop
    agc_params = hcat(tars, epses)'

    for i in 1:n_output_samples

        @views window = signal[(i - 1) * decimation_factor + 1 : i * decimation_factor]
        sos_output, sos_state = soscascade(window, ear_filters, sos_state)

        sos_output = sos_output' # BUG why transpose?
        output = clamp!(sos_output, 0, Inf) # Half Wave Rectify
        output[1, 1] = 0 # Test Hack to make inversion easier.
        output[2, 1] = 0

        if useagc
            output, agc_state = agc(output, agc_params, agc_state)
        end
        if differ
            # @views output = vcat(output[1:1, :], output[1:end-1, :] .- output[2:end, :])
            output[2:end, :] .= @view(output[1:end-1, :]) .- @view(output[2:end, :])
            clamp!(output, 0, Inf)
        end
        if decimation_factor > 1

            output, dec_state = sosfilters(
                output,
                reshape(dec_filt, 1, :), #[:, np.newaxis])
                dec_state)

        end

        y[:, i] = output[:, end]
    end

    y[3:end, :]
end


text:
function update_parameters!(net::Net{GPUBackend}, method::Adam,
                            alpha::Float64, epsilon::Float64, beta1::Float64, beta2::Float64,
                            m, v, gradient, param_blob, t, data_type)

  # update biased gradient moment estimates, m and v
  # m_t <- beta1*m_{t-1} + (1-beta1)*g_t

  CuBLAS.scal(net.backend.cublas_ctx, length(m), convert(data_type, beta1), m.ptr, 1)
  CuBLAS.axpy(net.backend.cublas_ctx, length(m), convert(data_type, 1-beta1), gradient.ptr, 1, m.ptr, 1)

  # now we need g2 = g.^2
  # Are we better to preallocate (heavy memory usage) or allocate each time (slower?)
  tmp = make_blob(net.backend, data_type, size(gradient))
  copy!(tmp, gradient)
  CuVec.pow!(net.backend, tmp, 2)
  CuBLAS.scal(net.backend.cublas_ctx, length(v), convert(data_type, beta2), tmp.ptr, 1)
  CuBLAS.axpy(net.backend.cublas_ctx, length(v), convert(data_type, 1-beta2), tmp.ptr, 1, v.ptr, 1)

  # Correct bias and calculate effective stepsize for timestep t (this is just scalar arithmetic)
  alpha_t = alpha * sqrt(1 - beta2^t) / (1 - beta1^t)

  # Update param <- param - alpha_t * m / (sqrt(v) + epsilon)
  copy!(tmp, v) # we'll reuse tmp for this to save reallocation
  CuVec.pow!(net.backend, tmp, convert(data_type, 0.5)) # tmp <- sqrt(v)
  CuVec.add_scal!(net.backend, tmp, convert(data_type, epsilon)) # tmp += epsilon

  CuVec.div2!(net.backend, m, tmp) # tmp = m/tmp

  CuBLAS.axpy(net.backend.cublas_ctx, length(param_blob), convert(data_type, -alpha_t), tmp.ptr, 1, param_blob.ptr, 1)

  destroy(tmp)
end


text:
const MAX_COLOUR = 255
const NUM_COLOURS = MAX_COLOUR + 1

struct Palette
    p::Vector{RGB{N0f8}}
end

const palettes = Dict{Symbol, Palette}()

Base.getindex(pal::Palette, idx::Integer) = pal.p[idx]

"""
    readpalette(stream::IO)
Read palette data from `stream`. The data should be
formatted as up to $NUM_COLOURS lines with three integers
each, separated by white space. The integers must fit into
a `UInt8`; they represent the RGB values of a 24-bit colour.
'#' marks the beginning of a comment which runs to the
end of a line and is ignored.
"""
function readpalette(stream::IO)
    i = 0
    pal = Vector{RGB{N0f8}}(undef, NUM_COLOURS)
    for line in readlines(stream)
        line = replace(strip(line), r"\s*#.*" => "")
        if isempty(line)
            continue
        end
        i += 1
        if i > NUM_COLOURS
            break
        end
        r, g, b = split(line)
        pal[i] = RGB{N0f8}(parse(Float32, r) / MAX_COLOUR,
                           parse(Float32, g) / MAX_COLOUR,
                           parse(Float32, b) / MAX_COLOUR)
    end
    while i < NUM_COLOURS
        i += 1
        pal[i] = RGB{N0f8}(0, 0, 0)
    end
    return Palette(pal)
end

"""
    palette(filename::AbstractString)
Reads and returns a `Palette` from `filename`.
    palette(name::Symbol)
Reads and returns a palette from a user supplied file
or from a standard palette file. Defaults to a palette
of $NUM_COLOURS greys.
"""
palette(filename::AbstractString) = palette(Symbol(filename))

function palette(name::Symbol)
    paths = [
        "$name",
        joinpath(dirname(@__FILE__), "maps", "$name.map")
    ]
    if haskey(palettes, name)
        return palettes[name]
    elseif any(ispath, paths)
        path = paths[findfirst(ispath, paths)]
        if !haskey(palettes, name)
            pal = readpalette(open(path))
            palettes[name] = pal
        end
        return palettes[name]
    else
        if !haskey(palettes, :grey)
            pal = Vector{RGB{N0f8}}(undef, NUM_COLOURS)
            for i in 1:NUM_COLOURS
                g = (i - 1) / MAX_COLOUR
                pal[i] = RGB{N0f8}(g, g, g)
            end
            palettes[:gray] = palettes[:grey] = Palette(pal)
        end
        return palettes[:grey]
    end
end


text:
"""
    StickyMessages(io::IO; ansi_codes=io isa Base.TTY && !Sys.iswindows())

A `StickyMessages` type manages the display of a set of persistent "sticky"
messages in a terminal. That is, messages which are not part of the normal
scrolling output. Each message is identified by a label and may may be added to
the set using `push!(messages, label=>msg)`, and removed using
`pop!(messages, label)`, or `empty!()`.

Only a single StickyMessages object should be associated with a given TTY, as
the object manipulates the terminal scrolling region.
"""
mutable struct StickyMessages
    io::IO
    # Bool for controlling TTY escape codes.
    ansi_codes::Bool
    # Messages is just used as a (short) OrderedDict here
    messages::Vector{Pair{Any,String}}
end

function StickyMessages(io::IO; ansi_codes=io isa Base.TTY &&
        # Give up on Windows for now, as libuv doesn't recognize the scroll region code.
        # Will need to be fixed in libuv and thence julia, but first libuv PR 1884 should merge. See
        # https://github.com/libuv/libuv/pull/1884
        # https://github.com/JuliaLang/libuv/commit/ed3700c849289ed01fe04273a7bf865340b2bd7e
                                    !Sys.iswindows())
    sticky = StickyMessages(io, ansi_codes, Vector{Pair{Any,String}}())
    # Ensure we clean up the terminal
    finalizer(sticky) do sticky
        # See also empty!()
        if !sticky.ansi_codes
            return
        end
        prev_nlines = _countlines(sticky.messages)
        if prev_nlines > 0
            # Clean up sticky lines. Hack: must be async to do the IO outside
            # of finalizer. Proper fix would be an uninstall event triggered by
            # with_logger.
            @async showsticky(sticky.io, prev_nlines, [])
        end
    end
    sticky
end

# Count newlines in a message or sequence of messages
_countlines(msg::String) = sum(c->c=='\n', msg)
_countlines(messages) = length(messages) > 0 ? sum(_countlines, messages) : 0
_countlines(messages::Vector{<:Pair}) = _countlines(m[2] for m in messages)

# Selected TTY cursor and screen control via ANSI codes
# * https://en.wikipedia.org/wiki/ANSI_escape_code#CSI_sequences
# * See man terminfo on linux, eg `tput csr $row1 $row2` and `tput cup $row $col`
# * For windows see https://docs.microsoft.com/en-us/windows/console/console-virtual-terminal-sequences
change_scroll_region!(io, rows::Pair) = write(io, "\e[$(rows[1]);$(rows[2])r")
change_cursor_line!(io, line::Integer) = write(io, "\e[$line;1H")
clear_to_end!(io) = write(io, "\e[J")

function showsticky(io, prev_nlines, messages)
    height,_ = displaysize(io)
    iob = IOBuffer()
    if prev_nlines > 0
        change_cursor_line!(iob, height + 1 - prev_nlines)
        clear_to_end!(iob)
    end
    # Set scroll region to the first N lines.
    #
    # Terminal scrollback buffers seem to be populated with a heuristic which
    # relies on the scrollable region starting at the first row, so to have
    # normal scrollback work we have to position sticky messages at the bottom
    # of the screen.
    linesrequired = _countlines(messages)
    if prev_nlines < linesrequired
        # Scroll screen up to preserve the lines which we will overwrite
        change_cursor_line!(iob, height - prev_nlines)
        write(iob, "\n"^(linesrequired-prev_nlines))
    end
    if prev_nlines != linesrequired
        change_scroll_region!(iob, 1=>height-linesrequired)
    end
    # Write messages. Avoid writing \n of last message to kill extra scrolling
    if !isempty(messages)
        change_cursor_line!(iob, height + 1 - linesrequired)
        for i = 1:length(messages)-1
            write(iob, messages[i][2])
        end
        write(iob, chop(messages[end][2]))
    end
    # TODO: Ideally we'd query the terminal for the line it was on before doing
    # all this and restore it if it's not in the new non-scrollable region.
    change_cursor_line!(iob, height - max(prev_nlines, linesrequired))
    # Write in one block to make the whole operation as atomic as possible.
    write(io, take!(iob))
    nothing
end

function Base.push!(sticky::StickyMessages, message::Pair)
    if !sticky.ansi_codes
        println(sticky.io, rstrip(message[2]))
        return
    end
    label,text = message
    endswith(text, '\n') || (text *= '\n';)
    prev_nlines = _countlines(sticky.messages)
    idx = findfirst(m->m[1] == label, sticky.messages)
    if idx === nothing
        push!(sticky.messages, label=>text)
    else
        sticky.messages[idx] = label=>text
    end
    showsticky(sticky.io, prev_nlines, sticky.messages)
end

function Base.pop!(sticky::StickyMessages, label)
    sticky.ansi_codes || return
    idx = findfirst(m->m[1] == label, sticky.messages)
    if idx !== nothing
        prev_nlines = _countlines(sticky.messages)
        deleteat!(sticky.messages, idx)
        showsticky(sticky.io, prev_nlines, sticky.messages)
    end
    nothing
end

function Base.empty!(sticky::StickyMessages)
    sticky.ansi_codes || return
    prev_nlines = _countlines(sticky.messages)
    empty!(sticky.messages)
    showsticky(sticky.io, prev_nlines, sticky.messages) # Resets scroll region
    nothing
end


text:
using Statistics: mean
using ColorSchemes: ColorScheme, colorschemes
using Luxor: @draw, background, circle, fontsize, line, Point, sethue, text
using CairoMakie: Axis, Colorbar, DataAspect, Figure,
    hidedecorations!, hidespines!, lines!, scatterlines!, text!,
    xlims!, ylims!

# TODO
# - text font for bases
# - legend for base_colors gradient
# - choose png or svg output
# - be able to ignore base_colors if not desired,
#   with base_colorscheme == flag_id it is white by default
# - make `constants` block configurable

"""
    plot_structure(structure; [sequence, layout_type, base_colors, base_colorscheme])

Plot a secondary structure to a PNG image.
"""
function plot_structure(structure::AbstractString;
                        sequence::AbstractString=" "^length(structure),
                        layout_type::Symbol=:simple,
                        base_colors::Vector=zeros(length(structure)),
                        base_colorscheme::ColorScheme=colorschemes[:flag_id])
    # constants
    base_radius = 10
    font_size = 20
    background_color = "white"
    base_circle_color = "black"
    base_text_color = "black"
    backbone_color = "black"
    basepair_color = "blue"
    scale_plot_coords = if layout_type == :circular
        200.0
    elseif layout_type == :turtle || layout_type == :puzzler
        1.3
    else
        2.5
    end
    x_pad = 5
    y_pad = 5

    function line_connect_bases(xa, ya, xb, yb, color)
        sethue(color)
	# TODO: line goes over the edge of a base, this is fixed in
	# draw_base by first overwriting with background color
	line(Point(xa, ya), Point(xb, yb), :stroke)
    end
    draw_backbone(xa, ya, xb, yb) =
	line_connect_bases(xa, ya, xb, yb, backbone_color)
    draw_basepair(xa, ya, xb, yb) =
	line_connect_bases(xa, ya, xb, yb, basepair_color)
    function draw_base(x, y, txt::AbstractString, fillcolor)
        p = Point(x, y)
        # overwrite with background
	sethue(background_color)
	circle(Point(x, y), base_radius, :fill)
        # base_color
        sethue(fillcolor)
	circle(p, base_radius, :fill)
        # base circle
	sethue(base_circle_color)
	circle(p, base_radius, :stroke)
        # text
        sethue(base_text_color)
	text(txt, p, halign=:center, valign=:middle)
    end

    n = length(structure)
    # TODO: we need the structure as a String for plot_coords,
    # and as a Pairtable to find basepairs
    pt = Pairtable(structure)
    length(sequence) == n ||
        throw(ArgumentError("structure and sequence must have same length"))
    length(base_colors) == n ||
        throw(ArgumentError("base_colors vector must have same length as structure"))

    xs, ys = plot_coords(structure; plot_type=layout_type)
    # center coords, and calculate width and height of image
    # TODO: not sure if optimal this way, seems ok-ish
    xs .*= scale_plot_coords
    ys .*= scale_plot_coords
    x_origin = mean(xs)
    y_origin = mean(ys)
    xs .-= x_origin
    ys .-= y_origin
    x_min, x_max = extrema(xs)
    y_min, y_max = extrema(ys)
    x_width = 2 * (ceil(max(abs(x_min), abs(x_max))) + base_radius + x_pad)
    y_width = 2 * (ceil(max(abs(y_min), abs(y_max))) + base_radius + y_pad)

    # Note:
    # in Luxor, (0,0) is the center,
    #           (-width/2, -height/2) is top left
    #       and (width/2, height/2) is lower right
    @draw begin
	background(background_color)
	fontsize(font_size)
	# backbone
	for i = 1:n
	    if i != n
                draw_backbone(xs[i], ys[i], xs[i+1], ys[i+1])
	    end
	end
	# basepairs
	for (i,j) in basepairs(pt)
            draw_basepair(xs[i], ys[i], xs[j], ys[j])
	end
        # bases (these must come last, as we overwrite too long
        # backbone and basepair lines here)
	for i = 1:n
            draw_base(xs[i], ys[i], string(sequence[i]),
                      base_colorscheme[base_colors[i]])
	end
    end x_width y_width
end

plot_structure(pt::Pairtable;
               sequence::AbstractString=" "^length(pt),
               layout_type::Symbol=:naview,
               base_colors::Vector=zeros(length(pt)),
               base_colorscheme::ColorScheme=colorschemes[:flag_id]) =
                   plot_structure(String(pt);
                                  sequence, base_colors, base_colorscheme)

"""
    plot_structure_makie(structure; [sequence, filepath, layout_type, colorscheme])

Plot a secondary structure to a PNG image or PDF file depending on `filepath` ending.
"""
function plot_structure_makie(
    structure::AbstractString;
    sequence::AbstractString=" "^length(structure),
    filepath::String = "",
    layout_type::Symbol=:simple,
    colorscheme::Symbol=:lightrainbow)

    fc = FoldCompound(sequence)
    partfn(fc)
    pt = Pairtable(structure)
    x_coords, y_coords = plot_coords(structure; plot_type=layout_type)

    markersize = 100 / sqrt(length(structure))
    positions = [(a, b) for (a, b) in zip(x_coords, y_coords)]
    pairs = basepairs(pt)
    f = Figure()
    ax = Axis(f[1, 1])
    xlims!(
        round(Int, minimum(x_coords)) - 20,
        round(Int, maximum(x_coords)) + 20)
    ylims!(
        round(Int, minimum(y_coords)) - 20,
        round(Int, maximum(y_coords)) + 20)
    hidedecorations!(ax)
    hidespines!(ax)

    for (i, j) in pairs
        lines!(
            [x_coords[i], x_coords[j]],
            [y_coords[i], y_coords[j]],
            color = :black,
            linestyle = :dot,
            linewidth = 1,
        )
    end
    scatterlines!(ax, x_coords, y_coords,
        markersize = markersize,
        markercolor = prob_of_basepairs(fc, pt),
        markercolormap = colorscheme,
        markercolorrange = (0, 1), # probabilities [0, 1]
        linewidth = 3,
    )
    text!(string.(collect(sequence)),
        position = positions,
        align = (:center, :center),
        textsize = markersize / 2,
    )
    Colorbar(
        f[2, 1],
        vertical = false,
        colormap = colorscheme,
        width = 500,
        height = 12,
    )
    ax.aspect = DataAspect()
    if ! isempty(filepath)
        save(filepath, f)
    end
    return f
end


text:
using DiffEqEnvironments
using Test
using SafeTestsets

@time @safetestset "DiffEqEnvironments.jl" begin
    @time @safetestset "Scalar ODE" begin include("test_scalar_autonomous_ode.jl") end
    @time @safetestset "LQR control" begin include("test_lqr_control.jl") end
    @time @safetestset "LTIQuadraticEnv" begin include("test_lti.jl") end
    @time @safetestset "IC sampler" begin include("test_ic_sampler.jl") end
end


text:
import Base.Random: rand, srand
import RNG: AbstractRNG, gen_seed, split_uint

"""
```julia
AbstractXorshift128 <: AbstractRNG{UInt64}
```

The base abstract type for `Xorshift128`, `Xorshift128Star` and `Xorshift128Plus`.
"""
abstract AbstractXorshift128 <: AbstractRNG{UInt64}

for (star, plus) in (
        (false, false),
        (false, true),
        (true, false),
    )
    rng_name = Symbol(string("Xorshift128", star ? "Star" : plus ? "Plus" :""))
    @eval begin
        type $rng_name <: AbstractXorshift128
            x::UInt64
            y::UInt64
            function $rng_name(seed::NTuple{2, UInt64}=gen_seed(UInt64, 2))
                r = new(0, 0)
                srand(r, seed)
                r
            end
        end

        $rng_name(seed::Integer) = $rng_name(split_uint(seed % UInt128))

        @inline function xorshift_next(r::$rng_name)
            t = r.x $ r.x << 23
            r.x = r.y
            r.y = t $ (t >> 3) $ r.y $ (r.y >> 24)
            $(star ? :(r.y * 2685821657736338717) :
              plus ? :(r.y + r.x) : :(r.y))
        end
    end
end

srand(r::AbstractXorshift128, seed::Integer) = srand(r, split_uint(seed % UInt128))

function srand(r::AbstractXorshift128, seed::NTuple{2, UInt64}=gen_seed(UInt64, 2))
    r.x = seed[1]
    r.y = seed[2]
    xorshift_next(r)
    xorshift_next(r)
    r
end

@inline rand(r::AbstractXorshift128, ::Type{UInt64}) = xorshift_next(r)


text:
module CARMA

include("Celerite.jl")
include("Kepler.jl")
include("CARMAKepler.jl")
include("Kalman.jl")
include("Green.jl")
include("LowPass.jl")

using .Kepler
using .CARMAKepler
using .Celerite
using .Green
using .Kalman
using .LowPass

export Green, Kalman, Celerite, CARMAKepler, Kepler, LowPass

end


text:
using ONNX, Flux, ProtoBuf
include("ops_tests.jl")

#test reshape one dim:
main_test("$ONNX_TEST_PATH/test_reshape_one_dim",
          read_output("$ONNX_TEST_PATH/test_reshape_extended_dims"),
          read_input("$ONNX_TEST_PATH/test_reshape_extended_dims")[1],
          read_input("$ONNX_TEST_PATH/test_reshape_extended_dims")[2])

#test reshape extended dim:
main_test("$ONNX_TEST_PATH/test_reshape_extended_dims",
          read_output("$ONNX_TEST_PATH/test_reshape_extended_dims"),
          read_input("$ONNX_TEST_PATH/test_reshape_extended_dims")[1],
          read_input("$ONNX_TEST_PATH/test_reshape_extended_dims")[2])

#test reshape reordered dim:
main_test("$ONNX_TEST_PATH/test_reshape_reordered_dims",
          read_output("$ONNX_TEST_PATH/test_reshape_reordered_dims"),
          read_input("$ONNX_TEST_PATH/test_reshape_reordered_dims")[1],
          read_input("$ONNX_TEST_PATH/test_reshape_reordered_dims")[2])

#test reshape reduced dim
main_test("$ONNX_TEST_PATH/test_reshape_reduced_dims",
          read_output("$ONNX_TEST_PATH/test_reshape_reduced_dims"),
          read_input("$ONNX_TEST_PATH/test_reshape_reduced_dims")[1],
          read_input("$ONNX_TEST_PATH/test_reshape_reduced_dims")[2])

## Transpose test:

main_test("$ONNX_TEST_PATH/test_transpose_all_permutations_0",
          read_output("$ONNX_TEST_PATH/test_transpose_all_permutations_0"),
          read_input("$ONNX_TEST_PATH/test_transpose_all_permutations_0")[1])

main_test("$ONNX_TEST_PATH/test_transpose_all_permutations_1",
          read_output("$ONNX_TEST_PATH/test_transpose_all_permutations_1"),
          read_input("$ONNX_TEST_PATH/test_transpose_all_permutations_1")[1])

main_test("$ONNX_TEST_PATH/test_transpose_all_permutations_2",
          read_output("$ONNX_TEST_PATH/test_transpose_all_permutations_2"),
          read_input("$ONNX_TEST_PATH/test_transpose_all_permutations_2")[1])

main_test("$ONNX_TEST_PATH/test_transpose_all_permutations_3",
          read_output("$ONNX_TEST_PATH/test_transpose_all_permutations_3"),
          read_input("$ONNX_TEST_PATH/test_transpose_all_permutations_3")[1])

main_test("$ONNX_TEST_PATH/test_transpose_all_permutations_4",
          read_output("$ONNX_TEST_PATH/test_transpose_all_permutations_4"),
          read_input("$ONNX_TEST_PATH/test_transpose_all_permutations_4")[1])

main_test("$ONNX_TEST_PATH/test_transpose_all_permutations_5",
          read_output("$ONNX_TEST_PATH/test_transpose_all_permutations_5"),
          read_input("$ONNX_TEST_PATH/test_transpose_all_permutations_5")[1])

## Test concat

#Test concat 1d axis 0
#ip = read_input("$ONNX_TEST_PATH/test_concat_1d_axis_0")
#main_test("$ONNX_TEST_PATH/test_concat_1d_axis_0",
#   read_output("$ONNX_TEST_PATH/test_concat_1d_axis_0"),
#        read_input("$ONNX_TEST_PATH/test_concat_1d_axis_0")[1],
#            read_input("$ONNX_TEST_PATH/test_concat_1d_axis_0")[2])

#Test concat 2d axis 0
#ip = read_input("$ONNX_TEST_PATH/test_concat_2d_axis_0")
#main_test("$ONNX_TEST_PATH/test_concat_2d_axis_0",
#    read_output("$ONNX_TEST_PATH/test_concat_2d_axis_0"),
#        read_input("$ONNX_TEST_PATH/test_concat_2d_axis_0")[1],
#            read_input("$ONNX_TEST_PATH/test_concat_2d_axis_0")[2])

#Test concat 2d axis 1
#ip = read_input("$ONNX_TEST_PATH/test_concat_2d_axis_1")
#main_test("$ONNX_TEST_PATH/test_concat_2d_axis_1",
#    read_output("$ONNX_TEST_PATH/test_concat_2d_axis_1"),
#        read_input("$ONNX_TEST_PATH/test_concat_2d_axis_1")[1],
#            read_input("$ONNX_TEST_PATH/test_concat_2d_axis_1")[2])

#Test concat 3d axis 0
#ip = read_input("$ONNX_TEST_PATH/test_concat_3d_axis_0")
#main_test("$ONNX_TEST_PATH/test_concat_3d_axis_0",
#    read_output("$ONNX_TEST_PATH/test_concat_3d_axis_0"),
#        read_input("$ONNX_TEST_PATH/test_concat_3d_axis_0")[1],
#            read_input("$ONNX_TEST_PATH/test_concat_3d_axis_0")[2])

#Test concat 3d axis 1
#ip = read_input("$ONNX_TEST_PATH/test_concat_3d_axis_1")
#main_test("$ONNX_TEST_PATH/test_concat_3d_axis_1",
#    read_output("$ONNX_TEST_PATH/test_concat_3d_axis_1"),
#        read_input("$ONNX_TEST_PATH/test_concat_3d_axis_1")[1],
#            read_input("$ONNX_TEST_PATH/test_concat_3d_axis_1")[2])

#Test concat 3d axis 2
#ip = read_input("$ONNX_TEST_PATH/test_concat_3d_axis_2")
#main_test("$ONNX_TEST_PATH/test_concat_3d_axis_2",
#    read_output("$ONNX_TEST_PATH/test_concat_3d_axis_2"),
#        read_input("$ONNX_TEST_PATH/test_concat_3d_axis_2")[1],
#            read_input("$ONNX_TEST_PATH/test_concat_3d_axis_2")[2])


text:
#
# Auxiliary Functions
#
# Not exactly related, but not exactly not

typealias Assignment Dict{Symbol, Any}

####################################################################################################
#                                   Dimensions

# make sure all the values are valid (in d)
_check_values_valid(v, d::Dimension) = (v in d) || not_in_dimension_error(v, d)

@inline function _check_values_valid(vs::AbstractArray, d::Dimension)
    for v in vs
        (v in d) || not_in_dimension_error(v, d)
    end
end

# make sure the Assignment only specifies one value per Dimension
_check_assignment_singleton(i::Int, ::Symbol) = nothing
_check_assignment_singleton(::Array, name::Symbol) = nonsingleton_assignment_error(name)

####################################################################################################
#                                   Factors

_dims_sizes{D<:Dimension}(dims::Vector{D}) = ntuple(k -> length(dims[k]), length(dims))

# dims are unique
_check_dims_unique{D<:Dimension}(dims::Vector{D}) =
    (length(dims) == 0) || _check_dims_unique(map(name, dims))

_check_dims_unique(dims::Vector{Symbol}) = allunique(dims) || non_unique_dims_error()

# dimensions need at least 1 value
_check_dims_singleton(d::Dimension) = isempty(d) && empty_support_error(name(d))

@inline function _check_dims_singleton{D<:Dimension}(dims::Vector{D})
    for d in dims
        isempty(d) && empty_support_error(name(d))
    end
end

# make sure all dims are valid (in φ)
_check_dims_valid(d::Symbol, φ::Factor) = (d in φ) || not_in_factor_error(d, φ)

@inline function _check_dims_valid(dims::Vector{Symbol}, φ::Factor)
    for d in dims
        (d in φ) || not_in_factor_error(d, φ)
    end
end

"""
    duplicate(A, dims::Integer...)

Repeates an array only through higer dimensions `dims`.

Custom version of `Base.repeat`, but only outer repetition, and only duplicates
the array for the number of times specified in `dims` for dimensions greater
than `ndims(A)`. If `dims` is empty, returns a copy of `A`.

```jldoctest
julia> duplicate(collect(1:3), 2)
3×2 Array{Int64,2}:
 1  1
 2  2
 3  3

julia> duplicate([1 3; 2 4], 3)
2×2×3 Array{Int64,3}:
[:, :, 1] =
 1  3
 2  4

[:, :, 2] =
 1  3
 2  4

[:, :, 3] =
 1  3
 2  4
```
"""
@inline function duplicate(A::Array, dims::Vararg{Integer})
    size_in = size(A)

    length_in = length(A)
    size_out = (size_in..., dims...)::Dims

    B = similar(A, size_out)

    # zero in dims means no looping
    @simd for index in 1:prod(dims)
        unsafe_copy!(B, (index - 1) * length_in + 1, A, 1, length_in)
    end

    return B
end


text:
struct Dataset{T, S}
    # if T is a Matrix, each of its columns should correspond to each of input data
    input::T
    target::S
end

function update_fitness!(fitness, pop)
    for i in 1:length(pop)
        fitness[i] = entry_pool[pop[i]].train_fitness
    end
end

function update_test_fitness!(fitness, pop)
    for i in 1:length(pop)
        fitness[i] = entry_pool[pop[i]].test_fitness
    end
end

max_stars_repo_name: kogad/JuliaGSGP

max_issues_repo_name: kogad/JuliaGSGP

max_forks_repo_name: kogad/JuliaGSGP


text:
@testset "single.jl" begin
    @testset "probabilities" begin
        probs = rand(10)
        consistent = Consistent(probs)
        sampler = Random.Sampler(Random.GLOBAL_RNG, consistent, Val(1))
        @test sampler isa Random.SamplerTrivial{<:Consistent{Float64,Bool}}

        # single sample
        for spl in (consistent, sampler)
            sample = @inferred(rand(spl))
            @test sample isa Tuple{Float64,Bool}
            @test first(sample) in probs
        end

        # multiple samples
        samples_oop = @inferred(rand(sampler, 100_000))
        @test samples_oop isa Vector{Tuple{Float64,Bool}}
        @test length(samples_oop) == 100_000
        samples_iip = similar(samples_oop)
        rand!(samples_iip, sampler)

        for samples in (samples_oop, samples_iip)
            @test all(first(s) in probs for s in samples)

            for prob in probs
                @test count(s -> first(s) == prob, samples) ≈ 10_000 atol = 300
                @test mean(t for (p, t) in samples if p == prob) ≈ prob atol = 2e-2
            end
        end
    end

    @testset "categorical" begin
        predictions = map(1:10) do _
            prediction = rand(5)
            lmul!(inv(sum(prediction)), prediction)
            return prediction
        end
        consistent = Consistent(predictions)
        sampler = Random.Sampler(Random.GLOBAL_RNG, consistent, Val(1))
        @test sampler isa Random.SamplerTrivial{<:Consistent{Vector{Float64},Int}}

        # single sample
        for spl in (consistent, sampler)
            sample = @inferred(rand(spl))
            @test sample isa Tuple{Vector{Float64},Int}
            @test first(sample) in predictions
            @test 1 ≤ last(sample) ≤ 5
        end

        # multiple samples
        samples_oop = @inferred(rand(sampler, 100_000))
        @test samples_oop isa Vector{Tuple{Vector{Float64},Int}}
        @test length(samples_oop) == 100_000
        samples_iip = similar(samples_oop)
        rand!(samples_iip, sampler)

        for samples in (samples_oop, samples_iip)
            @test all(first(s) in predictions for s in samples)
            @test all(1 ≤ last(s) ≤ 5 for s in samples)

            for prediction in predictions
                @test count(s -> first(s) == prediction, samples) ≈ 10_000 atol = 300
                @test proportions([t for (p, t) in samples if p == prediction], 1:5) ≈ prediction atol =
                    2e-2
            end
        end
    end

    @testset "Distributions" begin
        predictions = map(Normal, randn(10), rand(10))
        consistent = Consistent(predictions)
        sampler = Random.Sampler(Random.GLOBAL_RNG, consistent, Val(1))
        @test sampler isa Random.SamplerTrivial{<:Consistent{Normal{Float64},Float64}}

        # single sample
        for spl in (consistent, sampler)
            sample = @inferred(rand(spl))
            @test sample isa Tuple{Normal{Float64},Float64}
            @test first(sample) in predictions
        end

        # multiple samples
        samples_oop = @inferred(rand(sampler, 100_000))
        @test samples_oop isa Vector{Tuple{Normal{Float64},Float64}}
        @test length(samples_oop) == 100_000
        samples_iip = similar(samples_oop)
        rand!(samples_iip, sampler)

        for samples in (samples_oop, samples_iip)
            @test all(first(s) in predictions for s in samples)

            for prediction in predictions
                @test count(s -> first(s) == prediction, samples) ≈ 10_000 atol = 300
                @test mean(t for (p, t) in samples if p == prediction) ≈ mean(prediction) atol =
                    2e-2
                @test std(t for (p, t) in samples if p == prediction) ≈ std(prediction) atol =
                    2e-2
            end
        end
    end
end


text:
using Images
using TestImages
using Documenter
using ChaoticEncryption

DocMeta.setdocmeta!(ChaoticEncryption, :DocTestSetup, :(using ChaoticEncryption); recursive = true)
makedocs(
    sitename = "ChaoticEncryption",
    format = Documenter.HTML(sidebar_sitename=false),
    modules = [ChaoticEncryption, Images, TestImages],
    doctest = true,
    strict = :doctest,
    pages = [
        "Home" => "index.md",
        "Tutorials" => [
            "Pseudo-Random Number Generators" => "tutorials/prng.md",
            "Substitution algorithm" => "tutorials/substitution.md"
        ],
        "API Documentation" => [
            "Pseudo-Random Number Generators" => "apidocs/prngs.md",
            "Encryption/decryption algorithms" => "apidocs/algorithms.md"
        ],
        "Developer Documentation" => [
            "Encryption/decryption algorithms" => "devdocs/algorithms.md",
        ]
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/Saransh-cpp/ChaoticEncryption.jl.git",
    push_preview = true,
    devurl = "dev",
    versions = ["stable" => "v^", "v#.#.#", "dev" => "dev"],
)


text:
# TODO decide on full interface that should be supported by subtypes
# An AbstractTimeGF G should support at least:
# * G[k,l,t1,t2] where k,l are orbital indices and t1,t2 are TimeGridPoints
# * TimeDomain(G)
# * similar(G)
# * zero(G)
# * norbitals(G)
# * G.grid
"""
    AbstractTimeGF{T, scalar}

Abstract supertype for Keldysh Green's function objects.\\
An `AbstractTimeGF` object `G` supports indexing with `TimeGridPoint`s and interpolation with `BranchPoint`s.
"""
abstract type AbstractTimeGF{T, scalar} end

@enum GFSignEnum fermionic=-1 bosonic=1

Base.eltype(::Type{<:AbstractTimeGF{T, scalar}}) where {T, scalar} = T
Base.eltype(X::AbstractTimeGF) = eltype(typeof(X))

is_scalar(::Type{<:AbstractTimeGF{T, scalar}}) where {T, scalar} = scalar
is_scalar(X::AbstractTimeGF) = is_scalar(typeof(X))

Base.getindex(X::AbstractTimeGF{T,true}, t1::TimeGridPoint, t2::TimeGridPoint, greater=true) where {T} = X[1,1,t1,t2,greater]
Base.getindex(X::AbstractTimeGF{T,false}, t1::TimeGridPoint, t2::TimeGridPoint, greater=true) where {T} = X[:,:,t1,t2,greater]
Base.setindex!(X::AbstractTimeGF{T,true}, v, t1::TimeGridPoint, t2::TimeGridPoint) where {T} = X[1,1,t1,t2] = v
Base.setindex!(X::AbstractTimeGF{T,false}, v, t1::TimeGridPoint, t2::TimeGridPoint) where {T} = X[:,:,t1,t2] = v

function Base.getindex(G::AbstractTimeGF, b1::BranchEnum, b2::BranchEnum)
  grid = G.grid
  @assert b1 ∈ grid.contour && b2 ∈ grid.contour
  x = b1 == backward_branch ? reverse(grid[b1]) : grid[b1]
  y = b2 == backward_branch ? reverse(grid[b2]) : grid[b2]

  norb = norbitals(G)
  n = length(x)
  m = length(y)

  out = zeros(eltype(G), norb, norb, n, m)

  for (i,t1) in enumerate(x)
    for (j,t2) in enumerate(y)
      out[:,:,i,j] .= G[t1, t2]
    end
  end

  return norb == 1 ? out[1,1,:,:] : out
end

function Base.getindex(G::AbstractTimeGF, component::Symbol)
  grid = G.grid
  norb = norbitals(G)

  if component == :greater
    return G[backward_branch, forward_branch]
  elseif component == :lesser
    return G[forward_branch, backward_branch]
  elseif component == :matsubara
    X = real(-im .* G[imaginary_branch, imaginary_branch])
    return norb == 1 ? X[:,1] : X[:,:,:,1]
  elseif component == :retarded
    # TODO include heaviside function?
    ret = G[:greater] - G[:lesser]
    return ret
  elseif component == :advanced
    # TODO include heaviside function?
    adv = G[:lesser] - G[:greater]
    return adv
  elseif component == :leftmixing
    G[backward_branch, imaginary_branch]
  else
    throw(ArgumentError("component $component not recognized"))
  end
end

function interpolate(G::AbstractTimeGF{T, true}, t1::BranchPoint, t2::BranchPoint) where T
  grid = G.grid

  t1_l = find_lower(grid, t1)
  t1_r = grid[t1_l.cidx + 1]

  t2_l = find_lower(grid, t2)
  t2_r = grid[t2_l.cidx + 1]

  has_cut = t1.domain == t2.domain
  greater = heaviside(t1, t2)

  do_triangular_interp = has_cut && (t1_l.cidx == t2_l.cidx && t1_r.cidx == t2_r.cidx)

  if t1 == t2
    return _line_interp(t1.val,
                        t1_l.bpoint.val,
                        t1_r.bpoint.val,
                        G[t1_l, t1_l, greater],
                        G[t1_r, t1_r, greater])
  end

  if !do_triangular_interp
    return _square_interp(t1.val,
                          t1_l.bpoint.val,
                          t1_r.bpoint.val,
                          t2.val,
                          t2_l.bpoint.val,
                          t2_r.bpoint.val,
                          G[t1_l, t2_l, greater],
                          G[t1_l, t2_r, greater],
                          G[t1_r, t2_l, greater],
                          G[t1_r, t2_r, greater])
  else
    if greater # t1 >= t2
      return _tri_interp(t1.val,
                         t1_l.bpoint.val,
                         t1_r.bpoint.val,
                         t2.val,
                         t2_l.bpoint.val,
                         t2_r.bpoint.val,
                         G[t1_l, t2_l, greater],
                         G[t1_r, t2_l, greater],
                         G[t1_r, t2_r, greater])
    else
      return _tri_interp(t2.val,
                         t2_l.bpoint.val,
                         t2_r.bpoint.val,
                         t1.val,
                         t1_l.bpoint.val,
                         t1_r.bpoint.val,
                         G[t1_l, t2_l, greater],
                         G[t1_l, t2_r, greater],
                         G[t1_r, t2_r, greater])
    end
  end
end

function interpolate!(x, G::AbstractTimeGF{T, false}, t1::BranchPoint, t2::BranchPoint) where T
  grid = G.grid

  t1_l = find_lower(grid, t1)
  t1_r = grid[t1_l.cidx + 1]

  t2_l = find_lower(grid, t2)
  t2_r = grid[t2_l.cidx + 1]

  has_cut = t1.domain == t2.domain
  greater = heaviside(t1, t2)

  do_triangular_interp = has_cut && (t1_l.cidx == t2_l.cidx && t1_r.cidx == t2_r.cidx)

  if t1 == t2
    return _line_interp!(x,
                        t1.val,
                        t1_l.bpoint.val,
                        t1_r.bpoint.val,
                        G[t1_l, t1_l, greater],
                        G[t1_r, t1_r, greater])
  end

  if !do_triangular_interp
    return _square_interp!(x,
                          t1.val,
                          t1_l.bpoint.val,
                          t1_r.bpoint.val,
                          t2.val,
                          t2_l.bpoint.val,
                          t2_r.bpoint.val,
                          G[t1_l, t2_l, greater],
                          G[t1_l, t2_r, greater],
                          G[t1_r, t2_l, greater],
                          G[t1_r, t2_r, greater])
  else
    if greater # t1 >= t2
      return _tri_interp!(x,
                         t1.val,
                         t1_l.bpoint.val,
                         t1_r.bpoint.val,
                         t2.val,
                         t2_l.bpoint.val,
                         t2_r.bpoint.val,
                         G[t1_l, t2_l, greater],
                         G[t1_r, t2_l, greater],
                         G[t1_r, t2_r, greater])
    else
      return _tri_interp!(x,
                         t2.val,
                         t2_l.bpoint.val,
                         t2_r.bpoint.val,
                         t1.val,
                         t1_l.bpoint.val,
                         t1_r.bpoint.val,
                         G[t1_l, t2_l, greater],
                         G[t1_l, t2_r, greater],
                         G[t1_r, t2_r, greater])
    end
  end
end

function (G::AbstractTimeGF{T, true})(t1::BranchPoint, t2::BranchPoint) where T
  return interpolate(G, t1, t2)
end

function (G::AbstractTimeGF{T, false})(t1::BranchPoint, t2::BranchPoint) where T
  norb = norbitals(G)
  x = zeros(T, norb, norb)
  return interpolate!(x, G, t1, t2)
end

include("gf/generic.jl")
include("gf/full.jl")
include("gf/time_invariant_full.jl")
include("gf/keldysh.jl")
include("gf/time_invariant_keldysh.jl")
include("gf/imag.jl")


text:
using JuMP, NPZ, Tulip
using Printf
using Statistics
using LinearAlgebra

function Solve_func(year, tol)

	println("\nSolving weights for $year ...\n\n")

	# we only solve the weights for years where the targets have changed. If the
	# targets have not changed, we don't write the _input.npz file
	if isfile(string(year, "_input.npz"))
		array = npzread(string(year, "_input.npz"))
	else
		println("Skipping solver for $year \n")
		return nothing
	end

	A1 = array["A1"]
	A2 = array["A2"]
	b = array["b"]

	N = size(A1)[2]

        # scaling: determine a scaling vector with one value per constraint
	#  - the goal is to keep coefficients reasonably near 1.0
	#  - multiply each row of A1 and A2 by its specific scaling constant
	#  - multiply each element of the b target vector by its scaling constant
	#  - current approach: choose scale factors so that the sum of absolute values in each row of
	#    A1 and of A2 will equal the total number of records / 1000; maybe we can improve on this

	scale = (N / 1000.) ./ sum(abs.(A1), dims=2)

        A1s = scale .* A1
	A2s = scale .* A2
	bs = scale .* b

	model = Model(Tulip.Optimizer)
	set_optimizer_attribute(model, "OutputLevel", 1)  # 0=disable output (default), 1=show iterations
	set_optimizer_attribute(model, "IPM_IterationsLimit", 100)  # default 100 seems to be enough

	# r and s must each fall between 0 and the tolerance
	@variable(model, 0 <= r[1:N] <= tol)
	@variable(model, 0 <= s[1:N] <= tol)

	@objective(model, Min, sum(r[i] + s[i] for i in 1:N))

	# Ax = b  - use the scaled matrices and vector
	@constraint(model, [i in 1:length(bs)], sum(A1s[i,j] * r[j] + A2s[i,j] * s[j]
		                          for j in 1:N) == bs[i])

	optimize!(model)

	println("Termination status: ", termination_status(model))
	@printf "Objective = %.4f\n" objective_value(model)

	r_vec = value.(r)
	s_vec = value.(s)

	npzwrite(string(year, "_output.npz"), Dict("r" => r_vec, "s" => s_vec))

	println("\n")

	# quick checks on results

	# Did we satisfy constraints?
	rs = r_vec - s_vec
	b_calc = sum(rs' .* A1, dims=2)
	check = vec(b_calc) ./ b

	q = (0, .1, .25, .5, .75, .9, 1)
	println("Quantiles used below: ", q)

	println("\nQuantiles of ratio of calculated targets to intended targets: ")
	println(quantile!(check, q))

	# Are the ratios of new weights to old weights in bounds (within tolerances)?
	x = 1.0 .+ r_vec - s_vec  # note the .+
	println("\nQuantiles of ratio of new weight to initial weight: ")
	println(quantile!(x, q))

end

year_list = [x for x in 2012:2031]
tol_list = [0.40, 0.38, 0.35, 0.33, 0.30,
 			0.45, 0.45, 0.45, 0.45, 0.45,
			0.45, 0.45, 0.45, 0.45,	0.45,
			0.45, 0.45, 0.45, 0.45, 0.45]

# Run solver function for all years and tolerances (in order)
for i in zip(year_list, tol_list)
	Solve_func(i[1], i[2])
end

max_stars_repo_name: OpenSourcePolicyCenter/taxdata

max_issues_repo_name: open-source-economics/taxdata

max_forks_repo_name: open-source-economics/taxdata


text:
facts("Utilities") do

    T = rand_kruskal3(2, (10, 20, 30), true)
    @fact size(T) --> (10, 20, 30)

    context("_row_unfold()") do
        res = TensorDecompositions._row_unfold(T, 1)
        @fact size(res) --> (10, 600)

        res = TensorDecompositions._row_unfold(T, 2)
        @fact size(res) --> (20, 300)

        res = TensorDecompositions._row_unfold(T, 3)
        @fact size(res) --> (30, 200)
    end

    context("_col_unfold()") do
        res = TensorDecompositions._col_unfold(T, 1)
        @fact size(res) --> (600, 10)

        res = TensorDecompositions._col_unfold(T, 2)
        @fact size(res) --> (300, 20)

        res = TensorDecompositions._col_unfold(T, 3)
        @fact size(res) --> (200, 30)
    end

    context("tensorcontractmatrices()") do
        factors = TensorDecompositions._random_factors(size(T), (5, 2, 6))
        res = tensorcontractmatrices(T, factors)
        @fact size(res) --> (5, 2, 6)
    end

end


text:
struct SecondaryLimbDark{LD1 <: AbstractLimbDark,LD2 <: AbstractLimbDark,T} <: AbstractLimbDark
    primary_driver::LD1
    secondary_driver::LD2
    brightness_ratio::T
end

@doc raw"""
    SecondaryLimbDark(primary::AbstractLimbDark,
                      secondary::AbstractLimbDark;
                      brightness_ratio=1)
    SecondaryLimbDark(u_p::AbstractVector, u_s=u_p; kwargs...)

Compose two limb darkening laws together to add a secondary eclipse. If vectors of coefficients are provided, laws will automatically be constructed using [`PolynomialLimbDark`](@ref). The surface brightness ratio is given in terms of the host; e.g., if the companion is half as bright as the host, the ratio would be 0.5.

!!! note "Interface"
    [`SecondaryLimbDark`](@ref) only works with an orbit, since the companion's reference frame needs to be calculated. This means you can't call it using an impact parameter like `ld(b, r)` directly.

# Mathematical form
```math
f(t, r) = \frac{2f_p(t, r) + \eta r^2 f_s(t', r')}{1 + f_p(t, r) + \eta r^2 f_s(t', r')}
```
where $f_p$ is to the primary flux, $f_s$ is to the secondary flux, and $\eta$ is the surface brightness ratio. $t'$ and $r'$ correspond to the time and radius ratio from the companion's reference frame.

# Examples
```jldoctest second
# equal size and limb darkening
r = 1.0
u = [0.4, 0.26]
# companion is 1/10 as bright
brightness_ratio = 0.1
ld = SecondaryLimbDark(u; brightness_ratio)
orbit = SimpleOrbit(period=2, duration=0.5)
fp = ld(orbit, 0, r) # primary egress
fs = ld(orbit, 1, r) # secondary egress

fp ≈ brightness_ratio * fs

# output
true
```
"""
function SecondaryLimbDark(driver1::AbstractLimbDark, driver2::AbstractLimbDark; brightness_ratio=1.0)
    return SecondaryLimbDark(driver1, driver2, brightness_ratio)
end

function SecondaryLimbDark(u1::AbstractVector, u2=u1; kwargs...)
    driver1 = PolynomialLimbDark(u1)
    driver2 = PolynomialLimbDark(u2)
    return SecondaryLimbDark(driver1, driver2; kwargs...)
end

function compute(ld::SecondaryLimbDark, orbit::AbstractOrbit, t, ror; kwargs...)
    f1 = ld.primary_driver(orbit, t, ror; kwargs...)
    orbit2 = Orbits.flip(orbit, ror)
    f2 = ld.secondary_driver(orbit2, t, inv(ror); kwargs...)
    flux_ratio = ld.brightness_ratio * ror^2
    return (f1 + flux_ratio * f2) / (1 + flux_ratio)
end

function Base.show(io::IO, ld::SecondaryLimbDark{L1,L2}) where {L1,L2}
    r = ld.brightness_ratio
    np = L1.name.name
    ns = L2.name.name
    print(io, "SecondaryLimbDark($np, $ns, $r)")
end

function Base.show(io::IO, ::MIME"text/plain", ld::SecondaryLimbDark)
    p = ld.primary_driver
    s = ld.secondary_driver
    r = ld.brightness_ratio
    print(io, "SecondaryLimbDark\n primary: $p\n secondary: $s\n ratio: $r")
end


text:
using BinaryBuilder

name = "Mingw32"
version = v"5.0.4"

# sources to build, such as mingw32, our patches, etc....
sources = [
    "https://sourceforge.net/projects/mingw-w64/files/mingw-w64/mingw-w64-release/mingw-w64-v$(version).tar.bz2" =>
	"5527e1f6496841e2bb72f97a184fc79affdcd37972eaa9ebf7a5fd05c31ff803",
    "patches",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/mingw-*/

# We need newer configure scripts
update_configure_scripts

# This is where we'll install everything
sysroot=${prefix}/${target}/sys-root

# Patch mingw to build 32-bit cross compiler with GCC 7.1+ (This no longer needed with mingw 5.0.3+,
# so we let it fail, but leave it in incase we ever need to build an older mingw)
patch -p1 < $WORKSPACE/srcdir/patches/mingw_gcc710_i686.patch || true

# Install mingw headers
cd $WORKSPACE/srcdir/mingw-*/mingw-w64-headers
./configure --prefix=/ \
    --enable-sdk=all \
    --enable-secure-api \
    --host=${target}
make install DESTDIR=${sysroot}

# Build CRT
mkdir -p $WORKSPACE/srcdir/mingw_crt_build
cd $WORKSPACE/srcdir/mingw_crt_build
MINGW_CONF_ARGS=""

# If we're building a 32-bit build of mingw, add `--disable-lib64`
if [[ "${target}" == i686-* ]]; then
	MINGW_CONF_ARGS="${MINGW_CONF_ARGS} --disable-lib64"
else
    MINGW_CONF_ARGS="${MINGW_CONF_ARGS} --disable-lib32"
fi

$WORKSPACE/srcdir/mingw-*/mingw-w64-crt/configure \
    --prefix=/ \
    --host=${target} \
    ${MINGW_CONF_ARGS}
make -j${nproc}
make install DESTDIR=${sysroot}

# Build winpthreads
mkdir -p $WORKSPACE/srcdir/mingw_winpthreads_build
cd $WORKSPACE/srcdir/mingw_winpthreads_build
$WORKSPACE/srcdir/mingw-*/mingw-w64-libraries/winpthreads/configure \
    --prefix=/ \
    --host=${target} \
    --enable-static \
    --enable-shared

make -j${nproc}
make install DESTDIR=${sysroot}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Windows(:i686),
    Windows(:x86_64),
]

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

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

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

max_stars_repo_name: staticfloat/MingwBuilder

max_issues_repo_name: staticfloat/MingwBuilder

max_forks_repo_name: staticfloat/MingwBuilder


max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using Test
using Pkg.Artifacts
using NCDatasets

using RRTMGP
using RRTMGP.Device: array_type, array_device
using RRTMGP.Vmrs
using RRTMGP.LookUpTables
using RRTMGP.AtmosphericStates
using RRTMGP.Optics
using RRTMGP.Sources
using RRTMGP.BCs
using RRTMGP.Fluxes
using RRTMGP.AngularDiscretizations
using RRTMGP.RTE
using RRTMGP.RTESolver

import CLIMAParameters as CP
include(joinpath(pkgdir(RRTMGP), "parameters", "create_parameters.jl"))
# overriding some parameters to match with RRTMGP FORTRAN code
overrides = (; grav = 9.80665, molmass_dryair = 0.028964, molmass_water = 0.018016)
param_set = create_insolation_parameters(FT, overrides)

include("reference_files.jl")
include("read_rfmip_clear_sky.jl")
#---------------------------------------------------------------
function lw_rfmip(
    ::Type{OPC},
    ::Type{SRC},
    ::Type{VMR},
    ::Type{FT},
    ::Type{I},
    ::Type{DA},
) where {FT <: AbstractFloat, I <: Int, OPC, SRC, VMR, DA}
    opc = Symbol(OPC)
    lw_file = get_ref_filename(:lookup_tables, :clearsky, λ = :lw) # lw lookup tables
    lw_input_file = get_ref_filename(:atmos_state, :clearsky)      # clear-sky atmos state
    # reference data files for comparison
    flux_up_file = get_ref_filename(:comparison, :clearsky, λ = :lw, flux_up_dn = :flux_up, opc = opc)
    flux_dn_file = get_ref_filename(:comparison, :clearsky, λ = :lw, flux_up_dn = :flux_dn, opc = opc)

    FTA1D = DA{FT, 1}
    FTA2D = DA{FT, 2}
    max_threads = Int(256)
    exp_no = 1
    n_gauss_angles = 1

    # reading longwave lookup data
    ds_lw = Dataset(lw_file, "r")
    lookup_lw, idx_gases = LookUpLW(ds_lw, I, FT, DA)
    close(ds_lw)
    # reading rfmip data to atmospheric state
    ds_lw_in = Dataset(lw_input_file, "r")

    (as, sfc_emis, _, _, _) = setup_rfmip_as(ds_lw_in, idx_gases, exp_no, lookup_lw, FT, DA, VMR, max_threads)
    close(ds_lw_in)

    ncol, nlay, ngpt = as.ncol, as.nlay, lookup_lw.n_gpt
    nlev = nlay + 1
    op = OPC(FT, ncol, nlay, DA)                                   # allocating optical properties object
    src_lw = SRC(FT, DA, nlay, ncol)                               # allocating longwave source function object
    bcs_lw = LwBCs{FT, typeof(sfc_emis), Nothing}(sfc_emis, nothing)            # setting up boundary conditions
    fluxb_lw = FluxLW(ncol, nlay, FT, DA)                          # flux storage for bandwise calculations
    flux_lw = FluxLW(ncol, nlay, FT, DA)                           # longwave fluxes

    # initializing RTE solver
    slv = Solver(as, op, src_lw, nothing, bcs_lw, nothing, fluxb_lw, nothing, flux_lw, nothing)

    solve_lw!(slv, max_threads, lookup_lw)
    println("Timing ==================================================")
    for i in 1:5
        @time solve_lw!(slv, max_threads, lookup_lw)
        println("*************************************************")
    end
    # comparing with data from RRTMGP FORTRAN code
    flip_ind = nlev:-1:1

    ds_flux_up = Dataset(flux_up_file, "r")
    comp_flux_up = Array{FT}(ds_flux_up["rlu"][:][flip_ind, :, exp_no])
    close(ds_flux_up)

    ds_flux_dn = Dataset(flux_dn_file, "r")
    comp_flux_dn = Array{FT}(ds_flux_dn["rld"][:][flip_ind, :, exp_no])
    close(ds_flux_dn)

    max_err_flux_up = FT(maximum(abs.(Array(slv.flux_lw.flux_up) .- comp_flux_up)))
    max_err_flux_dn = FT(maximum(abs.(Array(slv.flux_lw.flux_dn) .- comp_flux_dn)))

    println("=======================================")
    println("Clear-sky longwave test - $opc")
    println("max_err_flux_up = $max_err_flux_up")
    println("max_err_flux_dn = $max_err_flux_dn")
    println("=======================================")
    toler = FT(1e-4)

    @test max_err_flux_up ≤ toler
    @test max_err_flux_dn ≤ toler
    return nothing
end

lw_rfmip(OneScalar, SourceLWNoScat, VmrGM, Float64, Int, array_type())
lw_rfmip(TwoStream, SourceLW2Str, VmrGM, Float64, Int, array_type())

max_stars_repo_name: climate-machine/JRRTMGP

max_issues_repo_name: climate-machine/JRRTMGP

max_forks_repo_name: climate-machine/JRRTMGP


text:
using Pkg
using PyCall
using Conda

try
    pyimport("andes")
catch
    @warn "PyCall is not configured to an existing Python env."
    @warn "Andes.jl will use Conda for PyCall and install andes."
    ENV["PYTHON"] = Conda.PYTHONDIR
    Pkg.build("PyCall")

    Conda.add("andes", channel="conda-forge")
end


text:
using ITensors
using ITensorInfiniteMPS

# H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ
function ising_mpo(s; J, h)
  N = length(s)
  a = AutoMPO()
  for n in 1:(N - 1)
    a .+= -J, "X", n, "X", n + 1
  end
  for n in 1:N
    a .+= -h, "Z", n
  end
  return H = MPO(a, s)
end

function ising_infinitempo_tensor(s⃗, l, r; J, h)
  s = only(s⃗)
  dₗ = 3 # The link dimension of the TFI
  Hmat = fill(ITensor(s', dag(s)), dₗ, dₗ)
  Hmat[1, 1] = op("Id", s)
  Hmat[2, 2] = op("Id", s)
  Hmat[3, 3] = op("Id", s)
  Hmat[2, 1] = -J * op("X", s)
  Hmat[3, 2] = -J * op("X", s)
  Hmat[3, 1] = -h * op("Z", s)
  H = emptyITensor()
  for i in 1:dₗ, j in 1:dₗ
    H += Hmat[i, j] * setelt(l => i) * setelt(r => j)
  end
  return H
end

# H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ
function ising_infinitempo(s; J, h)
  N = length(s)
  # Use a finite MPO to create the infinite MPO
  H = InfiniteMPO(N)
  # For the finite unit cell
  s¹ = addtags.(s, "c=1")
  l¹ = [Index(3, "Link,l=$n,c=1") for n in 1:N]
  for n in 1:N
    # TODO: use a CelledVector here to make the code logic simpler
    l¹ₙ = n == 1 ? replacetags(dag(l¹[N]), "c=1" => "c=0") : l¹[n - 1]
    r¹ₙ = l¹[n]
    H[n] = ising_infinitempo_tensor(s¹[n], l¹ₙ, r¹ₙ; J=J, h=h)
  end
  return H
end

# H = -J X₁X₂ - h Z₁
function ising_localop(s1::Index, s2::Index; J, h)
  return -J * op("X", s1) * op("X", s2) - h * op("Z", s1) * op("Id", s2)
end
# Version accepting IndexSet
ising_localop(s1, s2; kwargs...) = ising_localop(Index(s1), Index(s2); kwargs...)

function ising_infinitesumlocalops(s; J, h)
  N = length(s)
  H = InfiniteSumLocalOps(N)
  s∞ = CelledVector(s)
  return InfiniteSumLocalOps([ising_localop(s∞[n], s∞[n + 1]; J=J, h=h) for n in 1:N])
end

let
  N = 100
  s = siteinds("S=1/2", N)

  J = 1.0
  h = 1.5
  H = ising_mpo(s; J=J, h=h)
  ψ0 = randomMPS(s)

  sweeps = Sweeps(10)
  maxdim!(sweeps, 10)
  cutoff!(sweeps, 1E-10)
  energy, ψ = dmrg(H, ψ0, sweeps)
  @show energy / N

  n1 = N ÷ 2
  n2 = n1 + 1

  ψ = orthogonalize(ψ, n1)

  println("\nFinite MPS energy on bond (n1, n2) = $((n1, n2))")

  ϕ12 = ψ[n1] * ψ[n2]
  ham12 = ising_localop(s[n1], s[n2]; J=J, h=h)
  @show (ϕ12 * ham12 * prime(dag(ϕ12), commoninds(ϕ12, ham12)))[]

  # Get an infinite MPS that approximates
  # the finite MPS in the bulk.
  # `nsites` sets the number of sites in the unit
  # cell of the infinite MPS.
  nsites = 1
  ψ∞ = infinitemps_approx(ψ; nsites=nsites, nsweeps=10)

  println("\nInfinite MPS energy on a few different bonds")
  println("Infinite MPS has unit cell of size $nsites")

  # Compute a few energies
  for n1 in 1:3
    n2 = n1 + 1
    @show n1, n2
    ϕ12 = ψ∞.AL[n1] * ψ∞.AL[n2] * ψ∞.C[n2]
    s1 = siteind(ψ∞, n1)
    s2 = siteind(ψ∞, n2)
    ham12 = ising_localop(s1, s2; J=J, h=h)
    ϕ12ᴴ = prime(dag(ϕ12), commoninds(ϕ12, ham12))
    @show (ϕ12 * ham12 * ϕ12ᴴ)[] / norm(ϕ12)
  end
end


text:
using Conductor, IfElse, OrdinaryDiffEq, Unitful, ModelingToolkit
import Unitful: mV, mS, cm, pA, nA, ms, nS
import Conductor: Na, K

Vₘ = MembranePotential()

nav_kinetics = [
    Gate(AlphaBetaRates,
         αₘ = IfElse.ifelse(Vₘ == -40.0, 1.0, (0.1*(Vₘ + 40.0))/(1.0 - exp(-(Vₘ + 40.0)/10.0))),
         βₘ = 4.0*exp(-(Vₘ + 65.0)/18.0),
         p = 3)
    Gate(AlphaBetaRates,
         αₕ = 0.07*exp(-(Vₘ+65.0)/20.0),
         βₕ = 1.0/(1.0 + exp(-(Vₘ + 35.0)/10.0)))]

kdr_kinetics = [
    Gate(AlphaBetaRates,
         αₙ = IfElse.ifelse(Vₘ == -55.0, 0.1, (0.01*(Vₘ + 55.0))/(1.0 - exp(-(Vₘ + 55.0)/10.0))),
         βₙ = 0.125 * exp(-(Vₘ + 65.0)/80.0),
         p = 4)]

@named NaV = IonChannel(Sodium, nav_kinetics, 120mS/cm^2)
@named Kdr = IonChannel(Potassium, kdr_kinetics, 36mS/cm^2)
@named leak = PassiveChannel(Leak, 0.3mS/cm^2)

gradients = Equilibria([Na   =>  50.0mV, K    => -77.0mV, Leak => -54.4mV])
area = 0.629e-3cm^2
@named neuron1 = Soma([NaV,Kdr,leak], gradients, holding = 5000pA, area = ustrip(Float64, cm^2, area));
@named neuron2 = Soma([NaV,Kdr,leak], gradients, area = ustrip(Float64, cm^2, area));

# Synaptic model
syn∞ = 1/(1 + exp((-35 - Vₘ)/5))
τsyn = (1 - syn∞)/(1/40)
syn_kinetics = Gate(SteadyStateTau, :s, syn∞, τsyn, 1)
EGlut = Equilibrium{Conductor.Mixed}(0mV, :Glut)

@named Glut = Conductor.SynapticChannel(Leak, [syn_kinetics], EGlut, 30nS);
topology = [neuron1 => (neuron2, Glut)];
network = Conductor.Network([neuron1, neuron2], topology)

t = 250.
simul = Simulation(network, time = t*ms)

# Hand written problem
# Steady-state functions
minf(V0) = inv((4.0 + 0.1V0)*inv(1.0 - exp(-4.0 - 0.1V0)) + 4.0*exp(-3.6111111111111107 - 0.05555555555555555V0)) *
           (4.0 + 0.1V0) * inv(1.0 - exp(-4.0 - 0.1V0))
ninf(V0) = (0.55 + 0.01V0) * inv(1.0 - exp(-5.5 - 0.1V0)) * inv(0.125exp(-0.8125 - 0.0125V0) +
                                                                (0.55 + 0.01V0) * inv(1.0 - exp(-5.5 - 0.1V0)))
hinf(V0) = 0.07*exp(-3.25 - 0.05V0) * inv(0.07*exp(-3.25 - 0.05V0) + inv(1.0 + exp(-3.5 - 0.1V0)) )

# Initial conditions
V0 = -65.0
u0 = zeros(Float64, 12)
p = zeros(Float64, 18)

# Serialized MTK function initial conditions
# Neuron 1
u0[1]  = 0        # neuron1₊Isyn
u0[2]  = 0.005    # neuron1₊Iapp
u0[3]  = V0       # neuron1₊Vₘ
u0[4]  = minf(V0) # neuron1₊NaV₊m
u0[5]  = hinf(V0) # neuron1₊NaV₊h
u0[6]  = ninf(V0) # neuron1₊Kdr₊n

# Neuron 2
u0[7]  = 0.0      # neuron2₊Iapp
u0[8]  = V0       # neuron2₊Vₘ
u0[9]  = minf(V0) # neuron2₊NaV₊m
u0[10] = hinf(V0) # neuron2₊NaV₊h
u0[11] = ninf(V0) # neuron2₊Kdr₊n

# Synapse
u0[12] = 0.0      # Glut1₊s

# Parameters
# Neuron 1
p[2]   =  0.001     # neuron1₊cₘ
p[3]   =  0.000629  # neuron1₊aₘ
p[4]   =  50.0      # neuron1₊ENa
p[5]   =  -77.0     # neuron1₊EK
p[6]   =  -54.4     # neuron1₊El
p[7]   =  120.0     # neuron1₊NaV₊gbar
p[8]   =  36.0      # neuron1₊Kdr₊gbar
p[9]   =  0.3       # neuron1₊leak₊g

# Neuron 2
p[10]  =  0.001     # neuron2₊cₘ
p[11]  =  0.000629  # neuron2₊aₘ
p[12]  =  50.0      # neuron2₊ENa
p[13]  =  -77.0     # neuron2₊EK
p[14]  =  -54.4     # neuron2₊El
p[15]  =  120.0     # neuron2₊NaV₊gbar
p[16]  =  36.0      # neuron2₊Kdr₊gbar
p[17]  =  0.3       # neuron2₊leak₊g

# Synapse
p[1]   =  0.0       # EGlut
p[18]  =  3.0e-5    # Glut1₊gbar

function simple_synapse!(du, u, p, t)
    let neuron1₊Isyn        = @inbounds(u[1]),
        neuron1₊Iapp        = @inbounds(u[2]),
        neuron1₊Vₘ          = @inbounds(u[3]),
        neuron1₊NaV₊m       = @inbounds(u[4]),
        neuron1₊NaV₊h       = @inbounds(u[5]),
        neuron1₊Kdr₊n       = @inbounds(u[6]),
        neuron2₊Iapp        = @inbounds(u[7]),
        neuron2₊Vₘ          = @inbounds(u[8]),
        neuron2₊NaV₊m       = @inbounds(u[9]),
        neuron2₊NaV₊h       = @inbounds(u[10]),
        neuron2₊Kdr₊n       = @inbounds(u[11]),
        Glut1₊s             = @inbounds(u[12]),

        EGlut               = @inbounds(p[1]),
        neuron1₊cₘ          = @inbounds(p[2]),
        neuron1₊aₘ          = @inbounds(p[3]),
        neuron1₊ENa         = @inbounds(p[4]),
        neuron1₊EK          = @inbounds(p[5]),
        neuron1₊El          = @inbounds(p[6]),
        neuron1₊NaV₊gbar    = @inbounds(p[7]),
        neuron1₊Kdr₊gbar    = @inbounds(p[8]),
        neuron1₊leak₊g      = @inbounds(p[9]),
        neuron2₊cₘ          = @inbounds(p[10]),
        neuron2₊aₘ          = @inbounds(p[11]),
        neuron2₊ENa         = @inbounds(p[12]),
        neuron2₊EK          = @inbounds(p[13]),
        neuron2₊El          = @inbounds(p[14]),
        neuron2₊NaV₊gbar    = @inbounds(p[15]),
        neuron2₊Kdr₊gbar    = @inbounds(p[16]),
        neuron2₊leak₊g      = @inbounds(p[17]),
        Glut1₊gbar          = @inbounds(p[18])

        let
            du[1] = 0
            du[2] = 0
            du[3] = (-1 * neuron1₊aₘ * neuron1₊leak₊g * (neuron1₊Vₘ + -1neuron1₊El) +
                     -1 * neuron1₊Kdr₊gbar * neuron1₊aₘ * neuron1₊Kdr₊n^4.0 * (neuron1₊Vₘ + -1neuron1₊EK) +
                     -1 * neuron1₊NaV₊gbar * neuron1₊aₘ * neuron1₊NaV₊h * neuron1₊NaV₊m^3.0 * (neuron1₊Vₘ + -1neuron1₊ENa) +
                     neuron1₊Iapp + -1neuron1₊Isyn) * inv(neuron1₊aₘ) * inv(neuron1₊cₘ)

            du[4] = if neuron1₊Vₘ == -40.0
                1.0
            else
                (4.0 + 0.1neuron1₊Vₘ) * inv(1.0 + -1 * exp(-4.0 + -0.1neuron1₊Vₘ))
            end * (1 + -1neuron1₊NaV₊m) + -4.0 * neuron1₊NaV₊m *
                    exp(-3.6111111111111107 + -0.05555555555555555neuron1₊Vₘ)

            du[5] = -1 * neuron1₊NaV₊h * inv(1.0 + exp(-3.5 + -0.1neuron1₊Vₘ)) +
                    0.07 * exp(-3.25 + -0.05neuron1₊Vₘ) * (1 + -1neuron1₊NaV₊h)

            du[6] = (1 + -1neuron1₊Kdr₊n) * if neuron1₊Vₘ == -55.0
                0.1
            else
                (0.55 + 0.01neuron1₊Vₘ) * inv(1.0 + -1 * exp(-5.5 + -0.1neuron1₊Vₘ))
            end + -0.125 * neuron1₊Kdr₊n * exp(-0.8125 + -0.0125neuron1₊Vₘ)

            du[7] = 0
            du[8] = (-1 * Glut1₊gbar * Glut1₊s * (neuron2₊Vₘ + -1 * EGlut) +
                     -1 * neuron2₊aₘ * neuron2₊leak₊g * (neuron2₊Vₘ + -1neuron2₊El) +
                     -1 * neuron2₊Kdr₊gbar * neuron2₊aₘ * neuron2₊Kdr₊n^4.0 * (neuron2₊Vₘ + -1neuron2₊EK) +
                     -1 * neuron2₊NaV₊gbar * neuron2₊aₘ * neuron2₊NaV₊h * neuron2₊NaV₊m^3.0 * (neuron2₊Vₘ + -1neuron2₊ENa) +
                     neuron2₊Iapp) * inv(neuron2₊aₘ) * inv(neuron2₊cₘ)

            du[9] = if neuron2₊Vₘ == -40.0
                1.0
            else
                (4.0 + 0.1neuron2₊Vₘ) * inv(1.0 + -1 * exp(-4.0 + -0.1neuron2₊Vₘ))
            end * (1 + -1neuron2₊NaV₊m) + -4.0 * neuron2₊NaV₊m *
                    exp(-3.6111111111111107 + -0.05555555555555555neuron2₊Vₘ)

            du[10] = -1 * neuron2₊NaV₊h * inv(1.0 + exp(-3.5 + -0.1neuron2₊Vₘ)) +
                     0.07 * exp(-3.25 + -0.05neuron2₊Vₘ) * (1 + -1neuron2₊NaV₊h)

            du[11] = (1 + -1neuron2₊Kdr₊n) * if neuron2₊Vₘ == -55.0
                0.1
            else
                (0.55 + 0.01neuron2₊Vₘ) * inv(1.0 + -1 * exp(-5.5 + -0.1neuron2₊Vₘ))
            end + -0.125 * neuron2₊Kdr₊n * exp(-0.8125 + -0.0125neuron2₊Vₘ)

            du[12] = (-1Glut1₊s + inv(1 + exp(-7 // 1 + -1 // 5 * neuron1₊Vₘ))) *
                     inv(40.0 + -40.0 * inv(1 + exp(-7 // 1 + -1 // 5 * neuron1₊Vₘ)))
            nothing
        end
    end
end


text:
#####
## Static definitions for use with TintiNet.jl
#

## 1. Sequence (FASTA) Dictionary

const fastaCategories = [ "-",                  # Padding/unknown
                          "A", "C", "D", "E",
                          "F", "G", "H", "I",
                          "K", "L", "M", "N",
                          "P", "Q", "R", "S",
                          "T", "V", "W", "Y",   # 20 natural aminoacids
                          "X" ]                 # any aminoacid

# Create a hashing table for FastaCategories from Transformers.Basic.Vocabulary
const fastaVocab = Vocabulary(fastaCategories, "-")

## 2. Structure (SS3) Dictionary

const ss3Categories = [ "-",                    # Padding/ Unknown
                        "H", "E", "C" ]         # 3-state SS

# Create a hashing table for ss3Categories from Transformers.Basic.Vocabulary
const ss3Vocab = Vocabulary(ss3Categories, "-")

## 3. Structure (SS8) Dictionary

const ss8Categories = [ "-",                    # Padding/ Unknown
                        "H", "G", "I", "P",     # HELIX types
                        "E", "B",               # STRAND Types
                        "C", "T", "S" ]         # COIL types ## 8-state SS ( ! recent DSSP revision includes a 9th state)

# Create a hashing table for ss8Categories from Transformers.Basic.Vocabulary
const ss8Vocab = Vocabulary(ss8Categories, "-")

## 4. G-R-G ASA table
# RSA is the residue's ASA in the protein divided by the residue's ASA calculated when the residue is in an unfolded state,
# which is modeled as a Gly-R-Gly tripeptide.
const rsa_denominator_dictionary = Dict(
    "A" => 113,
    "R" => 241,
    "N" => 158,
    "D" => 151,
    "C" => 140,
    "Q" => 189,
    "E" => 183,
    "G" => 85,
    "H" => 194,
    "I" => 182,
    "L" => 180,
    "K" => 211,
    "M" => 204,
    "F" => 218,
    "P" => 143,
    "S" => 122,
    "T" => 146,
    "W" => 259,
    "Y" => 229,
    "V" => 160
)


text:
using Tables
using MLJBase
using StatsBase: sample, mean, quantile
using Statistics
using CategoricalArrays
using Distributions
using EvoTrees
using EvoTrees: logit, sigmoid
# import EvoTrees: EvoTreeRegressor, EvoTreeClassifier, EvoTreeCount, EvoTreeGaussian

##################################################
### Regression - small data
##################################################
features = rand(10_000) .* 5 .- 2
X = reshape(features, (size(features)[1], 1))
Y = sin.(features) .* 0.5 .+ 0.5
Y = logit(Y) + randn(size(Y))
Y = sigmoid(Y)
y = Y
X = Tables.table(X)

# @load EvoTreeRegressor
# linear regression
tree_model = EvoTreeRegressor(max_depth=5, η=0.05, nrounds=10)
# logistic regression
tree_model = EvoTreeRegressor(loss=:logistic, max_depth=5, η=0.05, nrounds=10)
# quantile regression
tree_model = EvoTreeRegressor(loss=:quantile, α=0.75, max_depth=5, η=0.05, nrounds=10)

tree = machine(tree_model, X, y)
train, test = partition(eachindex(y), 0.7, shuffle=true); # 70:30 split
fit!(tree, rows=train, verbosity=1)

tree.model.nrounds += 10
fit!(tree, rows=train, verbosity=1)

# predict on train data
pred_train = predict(tree, selectrows(X, train))
println(mean(abs.(pred_train - selectrows(Y, train))))

# predict on test data
pred_test = predict(tree, selectrows(X, test))
println(mean(abs.(pred_test - selectrows(Y, test))))

##################################################
### classif
##################################################
# features = rand(10_000) .* 5 .- 2
# X = reshape(features, (size(features)[1], 1))
# Y = sin.(features) .* 0.5 .+ 0.5
# Y = logit(Y) + randn(size(Y))
# Y = sigmoid(Y)
# y = Int.(round.(Y)) .+ 1
# y = _levels[y]
# # y = string.(y)
# y = CategoricalArray(y, ordered=false)
# X = Tables.table(X)
# X_matrix = MLJBase.matrix(X)

X, y = @load_crabs
X_matrix = matrix(X)

# define hyperparameters
tree_model = EvoTreeClassifier(max_depth=4, η=0.05, λ=0.0, γ=0.0, nbins=32, nrounds=20, metric=:none)
model = fit_evotree(tree_model, X_matrix, y);
pred = predict(model, X_matrix)
pred_cat = pred .> 0.5
sum((y .== "B") .== pred_cat[:,1]) / length(y)

# @load EvoTreeRegressor
tree = machine(tree_model, X, y)
train, test = partition(eachindex(y), 0.7, shuffle=true); # 70:30 split
fit!(tree, rows=train, verbosity=1)
report(tree)[:feature_importances]

tree.model.nrounds += 10
fit!(tree, rows=train, verbosity=1)
report(tree)[:feature_importances]

pred_train = predict(tree, selectrows(X, train))
pred_train_mode = predict_mode(tree, selectrows(X, train))
println(cross_entropy(pred_train, selectrows(y, train)) |> mean)
println(sum(pred_train_mode .== y[train]))

pred_test = predict(tree, selectrows(X, test))
pred_test_mode = predict_mode(tree, selectrows(X, test))
println(cross_entropy(pred_test, selectrows(y, test)) |> mean)
println(sum(pred_test_mode .== y[test]))
pred_test_mode = predict_mode(tree, selectrows(X, test))

# using LossFunctions, Plots
# evo_model = EvoTreeClassifier(max_depth=6, η=0.05, λ=1.0, γ=0.0, nrounds=10, nbins=64)
# evo = machine(evo_model, X, y)
# r = range(evo_model, :nrounds, lower=1, upper=500)
# @time curve = learning_curve!(evo, range=r, resolution=10, measure=HingeLoss())
# plot(curve.parameter_values, curve.measurements)

##################################################
### regression - Larger data
##################################################
features = rand(1_000_000, 100)
# features = rand(100, 10)
X = features
Y = rand(size(X, 1))
𝑖 = collect(1:size(X, 1))

# train-eval split
𝑖_sample = sample(𝑖, size(𝑖, 1), replace=false)
train_size = 0.8
𝑖_train = 𝑖_sample[1:floor(Int, train_size * size(𝑖, 1))]
𝑖_eval = 𝑖_sample[floor(Int, train_size * size(𝑖, 1)) + 1:end]

X_train, X_eval = X[𝑖_train, :], X[𝑖_eval, :]
Y_train, Y_eval = Y[𝑖_train], Y[𝑖_eval]

# @load EvoTreeRegressor
tree_model = EvoTreeRegressor(
    loss=:linear, metric=:mae,
    nrounds=10,
    λ=0.0, γ=0.0, η=0.1,
    max_depth=6, min_weight=1.0,
    rowsample=0.5, colsample=0.5, nbins=32)

X = Tables.table(X);
X = Tables.rowtable(X);
# X = Tables.columntable(X)
# X_matrix = MLJBase.matrix(X)

# typeof(X)
@time tree = machine(tree_model, X, Y);
train, test = partition(eachindex(Y), 0.8, shuffle=true); # 70:30 split
@time fit!(tree, rows=train, verbosity=1, force=true)

using LossFunctions
using MLJ
r = range(tree_model, :nrounds, lower=1, upper=100)
m = rms
@time curve = learning_curve!(evo, range=r, resolution=100, measure=m)

tree.model.nrounds += 1
@time update(tree.model, 0, tree.fitresult, tree.cache, X, Y);

tree.model.nrounds += 1
@time fit!(tree, rows=train, verbosity=1)
# @time MLJBase.fit!(tree, rows=train, verbosity=1)

# yhat = MLJBase.predict(tree.model, tree.fitresult, MLJ.selectrows(X,test))
pred_train = predict(tree, selectrows(X, train))
mean(abs.(pred_train - selectrows(Y, train)))

##################################################
### count - Larger data
##################################################
features = rand(100_000, 100)
# features = rand(100, 10)
X = features
Y = rand(UInt8, size(X, 1))
𝑖 = collect(1:size(X, 1))

# train-eval split
𝑖_sample = sample(𝑖, size(𝑖, 1), replace=false)
train_size = 0.8
𝑖_train = 𝑖_sample[1:floor(Int, train_size * size(𝑖, 1))]
𝑖_eval = 𝑖_sample[floor(Int, train_size * size(𝑖, 1)) + 1:end]

X_train, X_eval = X[𝑖_train, :], X[𝑖_eval, :]
Y_train, Y_eval = Y[𝑖_train], Y[𝑖_eval]

# @load EvoTreeRegressor
tree_model = EvoTreeCount(
    loss=:poisson, metric=:poisson,
    nrounds=10,
    λ=0.0, γ=0.0, η=0.1,
    max_depth=6, min_weight=1.0,
    rowsample=0.5, colsample=0.5, nbins=32)

X = Tables.table(X)
X = Tables.rowtable(X)
X = Tables.columntable(X)
X_matrix = MLJBase.matrix(X)

# typeof(X)
@time tree = machine(tree_model, X, Y)
train, test = partition(eachindex(Y), 0.8, shuffle=true); # 70:30 split
@time fit!(tree, rows=train, verbosity=1, force=true)

tree.model.nrounds += 10
@time MLJBase.update(tree.model, 0, tree.fitresult, tree.cache, X, Y)

tree.model.nrounds += 10
@time fit!(tree, rows=train, verbosity=1)
# @time MLJBase.fit!(tree, rows=train, verbosity=1)

# yhat = MLJBase.predict(tree.model, tree.fitresult, MLJ.selectrows(X,test))
pred = predict(tree, selectrows(X, train))
pred_mean = predict_mean(tree, selectrows(X, train))
pred_mode = predict_mode(tree, selectrows(X, train))

##################################################
### Gaussian - Larger data
##################################################
features = rand(100_000, 100)
# features = rand(100, 10)
X = features
Y = rand(size(X, 1))
𝑖 = collect(1:size(X, 1))

# train-eval split
𝑖_sample = sample(𝑖, size(𝑖, 1), replace=false)
train_size = 0.8
𝑖_train = 𝑖_sample[1:floor(Int, train_size * size(𝑖, 1))]
𝑖_eval = 𝑖_sample[floor(Int, train_size * size(𝑖, 1)) + 1:end]

X_train, X_eval = X[𝑖_train, :], X[𝑖_eval, :]
Y_train, Y_eval = Y[𝑖_train], Y[𝑖_eval]

# @load EvoTreeRegressor
tree_model = EvoTreeGaussian(
    loss=:gaussian, metric=:gaussian,
    nrounds=10,
    λ=0.0, γ=0.0, η=0.1,
    max_depth=6, min_weight=1.0,
    rowsample=0.5, colsample=0.5, nbins=32)

X = Tables.table(X)
X_matrix = MLJBase.matrix(X)

# typeof(X)
@time tree = machine(tree_model, X, Y)
train, test = partition(eachindex(Y), 0.8, shuffle=true); # 70:30 split
@time fit!(tree, rows=train, verbosity=1, force=true)

tree.model.nrounds += 10
@time MLJBase.update(tree.model, 0, tree.fitresult, tree.cache, X, Y)

tree.model.nrounds += 10
@time fit!(tree, rows=train, verbosity=1)
# @time MLJBase.fit!(tree, rows=train, verbosity=1)

# yhat = MLJBase.predict(tree.model, tree.fitresult, MLJ.selectrows(X,test))
pred = predict(tree, selectrows(X, train))
pred_mean = predict_mean(tree, selectrows(X, train))
pred_mode = predict_mode(tree, selectrows(X, train))
mean(abs.(pred_mean - selectrows(Y, train)))

q_20 = quantile.(pred, 0.20)
q_20 = quantile.(pred, 0.80)

#########################################
# MLJ2 test
#########################################
using EvoTrees
using MLJModelInterface
using MLJBase
using StatsBase: sample, mean, quantile
using Tables

X = rand(1_000_000, 100);
Y = rand(size(X, 1))

# @load EvoTreeRegressor
tree_model = EvoTreeRegressor(
    loss=:linear, metric=:mae,
    nrounds=10,
    λ=0.0, γ=0.0, η=0.1,
    max_depth=6, min_weight=1.0,
    rowsample=0.5, colsample=0.5, nbins=32)

X = Tables.table(X);
# X = Tables.rowtable(X);
# X = Tables.columntable(X);
# X_matrix = MLJBase.matrix(X);

# typeof(X)
@time tree = machine(tree_model, X, Y);
train, test = partition(eachindex(Y), 0.8, shuffle=true); # 70:30 split
@time fit!(tree, rows=train, verbosity=1, force=false)

tree.model.nrounds += 1
@time fit!(tree, rows=train, verbosity=1)