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

tHorz = 24
m = Model()

## Parameters
#: Gas Turbine
aPowGasTeLoad = 0.005
aFuelEload = 0.0001
aEmissFactEload = 0.002
aAuxRateGas = 1e-03
#: Steam Turbine Parameters
aPowHpEload = 1e-01
aPowIpEload = 1e-02
aPowLpEload = 1e-03

aDacSteaEload = 1.
aAuxRateStea = 1e-08
aCcRebDutyEload = 1e-04

#
#
aSteaUseRatePcc = 0.001
aPowUseRatePcc = 0.022

aSteaUseRateDacAir = 1e-03
aSteaUseRateDacFlue = 1e-03

aPowUseRateDacAir = 5e-03
aPowUseRateDacFlue = 2e-03

aCo2PccCapRate = .1
aSorbCo2CapFlue = 0.1
pSorbCo2CapAir = 1e-06
aPowUseRateDac = 0.0001

aPowUseRateComp = 0.001
#
cCostInvCombTurb = 1e+02
cCostInvSteaTurb = 1e+02
cCostInvTransInter = 1e+02
cCostInvPcc = 1e+02
cCostInvDac = 1e+03
cCostInvComp = 1e+01

# Cost parameters.
cCostFuel = 1e+02
pEmissionPrice = 1e+04
pCo2TranspPrice = 1e+01
pPowBasePrice = 1e+01
pCo2Credit = 1e+00

#vCapCombTurb = 3.
vCapSteaTurb = 2.
vCapTransInter = 5.
vCapPcc = 20.
vCapComp = 1000.

rnum = rand(Float64, tHorz)
@variable(m, 60 <= yGasTelecLoad[0:tHorz - 1] <= 100)
#@variable(m, 0 <= vCapCombTurb)
#@variable(m, 0 <= xPowGasTur[0:tHorz - 1] <= vCapCombTurb)
@variable(m, 0 <= xPowGasTur[0:tHorz - 1])

#@variable(m, 0 <= xPowSteaTur[0:tHorz - 1] <= vCapSteaTurb)
#@variable(m, 0 <= vCapSteaTurb)

#@variable(m, rnum[i+1] <= xPowOut[i = 0:tHorz - 1] <= vCapTransInter)
#@variable(m, xPowOut[i = 0:tHorz - 1] >= rnum[i + 1])
@variable(m, 0 <= xPowGross[0:tHorz - 1])
@variable(m, 0 <= xPowOut[ 0:tHorz - 1])

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

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

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

@variable(m, 0 <= xCcRebDuty[0:tHorz - 1])
@variable(m, 0 <= xDacSteaDuty[0:tHorz - 1])

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

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

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

# Dac-Flue
@variable(m, 0 <= a0Flue[0:tHorz - 1])  # We could generalize this for any number of hours.
@variable(m, 0 <= a1Flue[0:tHorz]) # State
@variable(m, 0 <= a2Flue[0:tHorz]) # State

@variable(m, 0 <= aR0Flue[0:tHorz - 1])  # Not-a-State
@variable(m, 0 <= aR1Flue[0:tHorz])  # State
@variable(m, 0 <= sFflue[0:tHorz])  # State
@variable(m, 0 <= sSflue[0:tHorz])  # State

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

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

# Dac-Air
@variable(m, 0 <= a0Air[0:tHorz - 1])  # We could generalize this for any number of hours.
@variable(m, 0 <= a1Air[0:tHorz]) # State
@variable(m, 0 <= a2Air[0:tHorz]) # State

@variable(m, 0 <= aR0Air[0:tHorz - 1])  # Not-a-State
@variable(m, 0 <= aR1Air[0:tHorz])  # State
@variable(m, 0 <= sFair[0:tHorz])  # State
@variable(m, 0 <= sSair[0:tHorz])  # State

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

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

@variable(m, 0 <= xDacSteaSlack[0:tHorz - 1])

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

@variable(m, xAuxPow[0:tHorz - 1])

# Constraints
# Gas Turbine
@constraint(m, powGasTur[i = 0:tHorz - 1],
            xPowGasTur[i] == aPowGasTeLoad * yGasTelecLoad[i]
           )
# 4
@constraint(m, fuelEq[i = 0:tHorz - 1],
            xFuel[i] == aFuelEload * yGasTelecLoad[i]
           )
# 5
@constraint(m, co2FuelEq[i = 0:tHorz - 1],
            xCo2Fuel[i] == aEmissFactEload * yGasTelecLoad[i]
           )

@constraint(m, auxPowGasT[i = 0:tHorz - 1],
            xAuxPowGasT[i] == aAuxRateGas * yGasTelecLoad[i])
# 6
# Steam
# 7a
@constraint(m, powHpEq[i = 0:tHorz - 1],
            xPowHp[i] == aPowHpEload * yGasTelecLoad[i])
# 7b
@constraint(m, powIpEq[i = 0:tHorz - 1],
            xPowIp[i] == aPowIpEload * yGasTelecLoad[i])

# 11
@constraint(m, powLpEq[i = 0:tHorz - 1],
            xPowLp[i] == aPowLpEload * yGasTelecLoad[i])
# 12
@constraint(m, powerSteaEq[i = 0:tHorz - 1],
            xPowSteaTur[i] == xPowHp[i] + xPowIp[i] + xPowLp[i])

@constraint(m, ccRebDutyEq[i = 0:tHorz - 1],
            xCcRebDuty[i] == aCcRebDutyEload * yGasTelecLoad[i])

@constraint(m, dacSteaDutyEq[i = 0:tHorz - 1],
            xDacSteaDuty[i] == aDacSteaEload * yGasTelecLoad[i])

@constraint(m, auxPowSteaTEq[i = 0:tHorz - 1],
            xAuxPowSteaT[i] == aAuxRateStea * yGasTelecLoad[i])

# PCC
# 13
#@constraint(m, co2CapPccEq[i = 0:tHorz - 1], xCo2CapPcc[i] == aCo2PccCapRate * xCo2Fuel[i])
@constraint(m, co2CapPccEq[i = 0:tHorz - 1],
            xCo2CapPcc[i] == 0.5 * xCo2Fuel[i])
# 14
@constraint(m, co2PccOutEq[i = 0:tHorz - 1],
            xCo2PccOut[i] == xCo2Fuel[i] - xCo2CapPcc[i])
# 15
@constraint(m, co2DacFlueInEq[i = 0:tHorz - 1],
            xCo2DacFlueIn[i] == xCo2PccOut[i] - vCo2PccVent[i])
# 16
# @constraint(m, co2CapPccIn[i = 0:tHorz - 1], xCo2CapPcc[i] <= vCapPcc)
# 17
@constraint(m, steamUsePccEq[i = 0:tHorz - 1],
            xSteaUsePcc[i] == aSteaUseRatePcc * xCo2CapPcc[i])
# 18
@constraint(m, powerUsePccEq[i = 0:tHorz - 1],
            xPowUsePcc[i] == aPowUseRatePcc * xCo2CapPcc[i])

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

# DAC-Flue
@constraint(m, a0FlueEq[i = 0:tHorz - 1], a0Flue[i] == vAbsFlue[i])
@constraint(m, aR0FlueEq[i = 0:tHorz - 1], aR0Flue[i] == vRegFlue[i])
@constraint(m, a1dFlueEq[i = 0:tHorz - 1], a1Flue[i + 1] == a0Flue[i])
@constraint(m, a2dFlueEq[i = 0:tHorz - 1], a2Flue[i + 1] == a1Flue[i])
@constraint(m, aRdFlueEq[i = 0:tHorz - 1], aR1Flue[i + 1] == aR0Flue[i])
@constraint(m, storeFflueeq[i = 0:tHorz - 1], sFflue[i + 1] == sFflue[i] - vAbsFlue[i] + aR1Flue[i])
@constraint(m, storeSflueeq[i = 0:tHorz - 1], sSflue[i + 1] == sSflue[i] - vRegFlue[i] + a2Flue[i])
@constraint(m, capDacFlueEq, sFflue[0] == 100)
# 22
@constraint(m, icA1FlueEq, a1Flue[0] == 0.)
@constraint(m, icA2FlueEq, a2Flue[0] == 0.)
@constraint(m, icAR1FlueEq, aR1Flue[0] == 0.)
#
@constraint(m, icSsFlueEq, sSflue[0] == 0.)
# 23
@constraint(m, co2StorDacFlueEq[i = 0:tHorz - 1], xCo2StorDacFlue[i] == aSorbCo2CapFlue * sSflue[i])
# 24
@constraint(m, co2CapDacFlueEq[i = 0:tHorz - 1], xCo2CapDacFlue[i] == aSorbCo2CapFlue * aR1Flue[i])
# 25
@constraint(m, steamUseDacFlueEq[i = 0:tHorz - 1], xSteaUseDacFlue[i] == aSteaUseRateDacFlue * xCo2CapDacFlue[i])
# 26
@constraint(m, powUseDacFlueEq[i = 0:tHorz - 1], xPowUseDacFlue[i] == aPowUseRateDacFlue * xCo2CapDacFlue[i])
# Equal to the amount vented, at least in flue mode.
@constraint(m, co2DacFlueVentEq[i = 0:tHorz - 1], xCo2DacVentFlue[i] == xCo2DacFlueIn[i] - xCo2CapDacFlue[i])

# DAC-Air
# Bluntly assume we can just take CO2 from air in pure form.
@constraint(m, a0AirEq[i = 0:tHorz - 1], a0Air[i] == vAbsAir[i])
@constraint(m, aR0AirEq[i = 0:tHorz - 1], aR0Air[i] == vRegAir[i])
@constraint(m, a1dAirEq[i = 0:tHorz - 1], a1Air[i + 1] == a0Air[i])
@constraint(m, a2dAirEq[i = 0:tHorz - 1], a2Air[i + 1] == a1Air[i])
@constraint(m, aRdAirEq[i = 0:tHorz - 1], aR1Air[i + 1] == aR0Air[i])
@constraint(m, storeFairEq[i = 0:tHorz - 1], sFair[i + 1] == sFair[i] - vAbsAir[i] + aR1Air[i])
@constraint(m, storeSaireq[i = 0:tHorz - 1], sSair[i + 1] == sSair[i] - vRegAir[i] + a2Air[i])
@constraint(m, capDacAirEq, sFair[0] == 1000)
#
@constraint(m, icA1AirEq, a1Air[0] == 0.)
@constraint(m, icA2AirEq, a2Air[0] == 0.)
@constraint(m, icAR1AirEq, aR1Air[0] == 0.)
#
@constraint(m, icSsAirEq, sSair[0] == 0.)
#
@constraint(m, co2StorDacAirEq[i = 0:tHorz - 1], xCo2StorDacAir[i] == pSorbCo2CapAir * sSair[i])
# Money, baby.
@constraint(m, co2CapDacAirEq[i = 0:tHorz - 1], xCo2CapDacAir[i] == pSorbCo2CapAir * aR1Air[i])
#
@constraint(m, steamUseDacAirEq[i = 0:tHorz - 1], xSteaUseDacAir[i] == aSteaUseRateDacAir * xCo2CapDacAir[i])
#
@constraint(m, powUseDacAirEq[i = 0:tHorz - 1], xPowUseDacAir[i] == aPowUseRateDacAir * xCo2CapDacAir[i])

@constraint(m, dacSteaSlackEq[i = 0:tHorz - 1],
            xDacSteaSlack[i] == xDacSteaDuty[i] - xSteaUseDacFlue[i] - xSteaUseDacAir[i])

# Co2 Compression
# 27
@constraint(m, co2CompEq[i = 0:tHorz - 1], xCo2Comp[i] == xCo2CapPcc[i])
# 28
@constraint(m, powUseCompEq[i = 0:tHorz - 1], xPowUseComp[i] == aPowUseRateComp * xCo2Comp[i])
# 29
# @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 - 1], xCo2Vent[i] == vCo2PccVent[i] + xCo2DacVentFlue[i] - xCo2CapDacAir[i])

## Overall
#
#
@constraint(m, powGrossEq[i = 0:tHorz - 1],
            xPowGross[i] == xPowGasTur[i] + xPowSteaTur[i]
           )
@constraint(m, auxPowEq[i = 0:tHorz - 1],
            xAuxPow[i] == xAuxPowGasT[i] + xAuxPowSteaT[i])
@constraint(m, powOutEq[i = 0:tHorz - 1],
            xPowOut[i] == xPowGross[i] - xPowUsePcc[i] - xPowUseDacFlue[i] - xPowUseDacAir[i] - xPowUseComp[i] - xAuxPow[i]
           )

@expression(m, eObjfExpr, sum(cCostFuel * xFuel[i] +
                              pEmissionPrice * xCo2Vent[i] +
                              pCo2TranspPrice * xCo2Comp[i] -
                              pPowBasePrice * xPowOut[i] for i in 0:tHorz - 1))
#@expression(m, eObjfExpr, sum(-xPowOut[i] for i in 0:tHorz - 1))

@objective(m, Min, eObjfExpr)

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")
#format::MOI.FileFormats.FileFormat = MOI.FileFormats.FORMAT_AUTOMATIC

# Design decisions.
# vCapCombTurb
# vCapSteaTurb
# vCapTransInter
# vCapPcc
# vCapDac
# vCapComp
#
# Raw materials.
# xFuel
# sF0
# sS0
#
# penalties
# vCo2Vent
#

## vCo2PccVent * 2N

# Actual variables
# xPowGasTur * 2

# Co2 Data Frame
df_co = DataFrame(Symbol("Co2Fuel") => Float64[], # Pairs.
                  Symbol("Co2CapPcc") => Float64[],
                  Symbol("Co2PccOut") => Float64[],
                  Symbol("vCo2PccVent") => Float64[],
                  Symbol("Co2DacFlueIn") => Float64[],
                  Symbol("Co2StoreDacFlue") => Float64[],
                  Symbol("Co2CapDacFlue") => Float64[],
                  Symbol("Co2CapDacAir") => Float64[],
                  Symbol("Co2DacVentFlue") => Float64[],
                  Symbol("Co2Vent") => Float64[],
                 )
for i in 0:tHorz - 1
    push!(df_co, (
                  value(xCo2Fuel[i]),
               value(xCo2CapPcc[i]),
               value(xCo2PccOut[i]),
               value(vCo2PccVent[i]),
               value(xCo2DacFlueIn[i]),
               value(xCo2StorDacFlue[i]),
               value(xCo2CapDacFlue[i]),
               value(xCo2CapDacAir[i]),
               value(xCo2DacVentFlue[i]),
               value(xCo2Vent[i])))
end

# Power Data Frame.
df_pow = DataFrame(Symbol("PowGasTur") => Float64[], # Pairs.
                  Symbol("PowSteaTurb") => Float64[],
                  Symbol("PowHp") => Float64[],
                  Symbol("PowIp") => Float64[],
                  Symbol("PowLp") => Float64[],
                  Symbol("PowUsePcc") => Float64[],
                  Symbol("PowUseDacFlue") => Float64[],
                  Symbol("PowUseDacAir") => Float64[],
                  Symbol("PowUseComp") => Float64[],
                  Symbol("xAuxPowGasT") => Float64[],
                  Symbol("xAuxPowSteaT") => Float64[],
                  Symbol("PowGross") => Float64[],
                  Symbol("PowOut") => Float64[],
                  Symbol("yGasTelecLoad") => Float64[],
                  Symbol("Demand") => Float64[]
                 )
for i in 0:tHorz-1
    push!(df_pow, (
                   value(xPowGasTur[i]),
                   value(xPowSteaTur[i]),
                   value(xPowHp[i]),
                   value(xPowIp[i]),
                   value(xPowLp[i]),
                   value(xPowUsePcc[i]),
                   value(xPowUseDacFlue[i]),
                   value(xPowUseDacAir[i]),
                   value(xPowUseComp[i]),
                   value(xAuxPowGasT[i]),
                   value(xAuxPowSteaT[i]),
                   value(xPowGross[i]),
                   value(xPowOut[i]),
                   value(yGasTelecLoad[i]),
                   rnum[i + 1]
                  ))
end

df_steam = DataFrame(
                     Symbol("xCcRebDuty") => Float64[],
                     Symbol("xDacSteaDuty") => Float64[],
                     Symbol("SteaUsePccFlue") => Float64[],
                     Symbol("SteaUseDacFlue") => Float64[],
                     Symbol("SteaUseDacAir") => Float64[],
                     Symbol("Fuel") => Float64[]
                    )
for i in 0:tHorz-1
    push!(df_steam, (
                     value(xCcRebDuty[i]),
                     value(xDacSteaDuty[i]),
                     value(xSteaUsePcc[i]),
                     value(xSteaUseDacFlue[i]),
                     value(xSteaUseDacAir[i]),
                     value(xFuel[i])
                    ),
         )
end

# Generate Dataframe of the DAC-Flue variables.
df_dac_flue = DataFrame(:sFflue => Float64[],
                     :sSflue => Float64[],
                     :vAbsFlue => Float64[],
                     :vRegFlue => Float64[],
                     :a0Flue => Float64[],
                     :a1Flue => Float64[],
                     :a2Flue => Float64[],
                     :aR0Flue => Float64[],
                     :aR1Flue => Float64[])
for i in 0:tHorz - 1
    push!(df_dac_flue,(value(sFflue[i]), value(sSflue[i]),
                  value(vAbsFlue[i]), value(vRegFlue[i]),
                  value(a0Flue[i]), value(a1Flue[i]), value(a2Flue[i]),
                  value(aR0Flue[i]), value(aR1Flue[i])))
end

# Generate Dataframe of the DAC-Air variables.
df_dac_air = DataFrame(:sFair => Float64[],
                     :sSair => Float64[],
                     :vAbsAir => Float64[],
                     :vRegAir => Float64[],
                     :a0Air => Float64[],
                     :a1Air => Float64[],
                     :a2Air => Float64[],
                     :aR0Air => Float64[],
                     :aR1Air => Float64[])
for i in 0:tHorz - 1
    push!(df_dac_air,(value(sFair[i]), value(sSair[i]),
                  value(vAbsAir[i]), value(vRegAir[i]),
                  value(a0Air[i]), value(a1Air[i]), value(a2Air[i]),
                  value(aR0Air[i]), value(aR1Air[i])))
end

#print(df_co)
#print(df_pow)
#print(df_steam)
#print(df_dac)

CSV.write("df_co.csv", df_co)
CSV.write("df_pow.csv", df_pow)
CSV.write("df_steam.csv", df_steam)
CSV.write("df_dac_flue.csv", df_dac_flue)
CSV.write("df_dac_air.csv", df_dac_air)

max_stars_repo_name: dthierry/princetonDacLti

max_issues_repo_name: dthierry/princetonDacLti

max_forks_repo_name: dthierry/princetonDacLti


text:
using Base.Test, Color

@test isleaftype(eltype(colormap("Grays")))

@test_throws ArgumentError colormap("Grays", N=10)

ext: jl


text:
# 1
# Comment example
# 2
println("Hello world!")
# 3
function this(is::Int64)
    multiple = "lines"
end
# 0
cellmd0 = """# This should be a multi-line string (Markdown)"""
# 4
println("yay")
# 0


text:
#-------------------------------------------------------------------
#* Model of a tray with reaction
#*-------------------------------------------------------------------
type trayReac
	trayReac()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Type=>"PP"
			)),
			DanaInteger(),
			stream(Dict{Symbol,Any}(
				:Brief=>"Feed stream",
				:PosX=>0,
				:PosY=>0.4932,
				:Symbol=>"_{in}"
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"liquid Sidestream",
				:Hidden=>true,
				:Symbol=>"_{outL}"
			)),
			vapour_stream(Dict{Symbol,Any}(
				:Brief=>"vapour Sidestream",
				:Hidden=>true,
				:Symbol=>"_{outV}"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet liquid stream",
				:PosX=>0.5195,
				:PosY=>0,
				:Symbol=>"_{inL}"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet vapour stream",
				:PosX=>0.4994,
				:PosY=>1,
				:Symbol=>"_{inV}"
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Outlet liquid stream",
				:PosX=>0.8277,
				:PosY=>1,
				:Symbol=>"_{outL}"
			)),
			vapour_stream(Dict{Symbol,Any}(
				:Brief=>"Outlet vapour stream",
				:PosX=>0.8043,
				:PosY=>0,
				:Symbol=>"_{outV}"
			)),
			fill(fraction()),
			fill(mol(Dict{Symbol,Any}(
				:Brief=>"Molar Holdup in the tray"
			)),(NComp)),
			mol(Dict{Symbol,Any}(
				:Brief=>"Molar liquid holdup"
			)),
			mol(Dict{Symbol,Any}(
				:Brief=>"Molar vapour holdup"
			)),
			energy(Dict{Symbol,Any}(
				:Brief=>"Total Energy Holdup on tray"
			)),
			volume_mol(Dict{Symbol,Any}(
				:Brief=>"Liquid Molar Volume"
			)),
			volume_mol(Dict{Symbol,Any}(
				:Brief=>"Vapour Molar volume"
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Height of clear liquid on plate"
			)),
			volume(),
			dens_mass(),
			dens_mass(),
			reaction_mol(Dict{Symbol,Any}(
				:Brief=>"Reaction resulting ethyl acetate",
				:DisplayUnit=>"mol/l/s"
			)),
			fill(conc_mol(Dict{Symbol,Any}(
				:Brief=>"Molar concentration",
				:Lower=>-1
			)),(NComp)),
			[
				:(OutletLiquid.z = vL * C),
				:(r3 = exp(-7150*"K"/OutletLiquid.T)*(4.85e4*C(1)*C(2) - 1.23e4*C(3)*C(4))*"l/mol/s"),
				:(M = ML*OutletLiquid.z + MV*OutletVapour.z),
				:(OutletVapour.T = VapourSideStream.T),
				:(OutletLiquid.T = LiquidSideStream.T),
				:(OutletVapour.P= VapourSideStream.P),
				:(OutletLiquid.P = LiquidSideStream.P),
				:(OutletLiquid.z= LiquidSideStream.z),
				:(OutletVapour.z= VapourSideStream.z),
				:(sum(OutletLiquid.z)= 1.0),
				:(vL = PP.LiquidVolume(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z)),
				:(vV = PP.VapourVolume(OutletVapour.T, OutletVapour.P, OutletVapour.z)),
				:(OutletVapour.T = OutletLiquid.T),
				:(OutletVapour.P = OutletLiquid.P),
				:(Vol = ML*vL),
				:(rhoL = PP.LiquidDensity(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z)),
				:(rhoV = PP.VapourDensity(InletVapour.T, InletVapour.P, InletVapour.z)),
				:(PP.LiquidFugacityCoefficient(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z)*OutletLiquid.z = PP.VapourFugacityCoefficient(OutletVapour.T, OutletVapour.P, yideal)*yideal),
				:(sum(OutletLiquid.z)= sum(OutletVapour.z)),
			],
			[
				"Molar Concentration","Reaction","Molar Holdup","Thermal Equilibrium Vapour Side Stream","Thermal Equilibrium Liquid Side Stream","Mechanical Equilibrium Vapour Side Stream","Mechanical Equilibrium Liquid Side Stream","Composition Liquid Side Stream","Composition Vapour Side Stream","Mol fraction normalisation","Liquid Volume","Vapour Volume","Thermal Equilibrium","Mechanical Equilibrium","","Liquid Density","Vapour Density","Chemical Equilibrium","",
			],
			[:PP,:NComp,],
			[:Inlet,:LiquidSideStream,:VapourSideStream,:InletLiquid,:InletVapour,:OutletLiquid,:OutletVapour,:yideal,:M,:ML,:MV,:E,:vL,:vV,:Level,:Vol,:rhoL,:rhoV,:r3,:C,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	Inlet::stream
	LiquidSideStream::liquid_stream
	VapourSideStream::vapour_stream
	InletLiquid::stream
	InletVapour::stream
	OutletLiquid::liquid_stream
	OutletVapour::vapour_stream
	yideal::Array{fraction}
	M::Array{mol}
	ML::mol
	MV::mol
	E::energy
	vL::volume_mol
	vV::volume_mol
	Level::length
	Vol::volume
	rhoL::dens_mass
	rhoV::dens_mass
	r3::reaction_mol
	C::Array{conc_mol}
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export trayReac
function setEquationFlow(in::trayReac)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
	addEquation(11)
	addEquation(12)
	addEquation(13)
	addEquation(14)
	addEquation(15)
	addEquation(16)
	addEquation(17)
	addEquation(18)
	addEquation(19)
end
function atributes(in::trayReac,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Icon]="icon/Tray"
	fields[:Brief]="Model of a tray with reaction."
	fields[:Info]="== Assumptions ==
* both phases (liquid and vapour) exists all the time;
* thermodymanic equilibrium with Murphree plate efficiency;
* no entrainment of liquid or vapour phase;
* no weeping;
* the dymanics in the downcomer are neglected.

== Specify ==
* the Feed stream;
* the Liquid inlet stream;
* the Vapour inlet stream;
* the Vapour outlet flow (OutletVapour.F);
* the reaction related variables.

== Initial ==
* the plate temperature (OutletLiquid.T)
* the liquid height (Level) OR the liquid flow OutletLiquid.F
* (NoComps - 1) OutletLiquid compositions
"
	drive!(fields,_)
	return fields
end
trayReac(_::Dict{Symbol,Any})=begin
	newModel=trayReac()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using EpiJ
using Test

@testset "EpiJ.jl" begin
    # Write your tests here.
end

max_stars_repo_name: JuliaHealth/EpiJ

max_issues_repo_name: JuliaHealth/EpiJ

max_forks_repo_name: JuliaHealth/EpiJ


text:
module FranklinMinify

include("files.jl")
using PyCall: pyimport
using Conda: pip
using Franklin: path

const rx_tags = Dict(
    "" => r"",
    "<script>" =>r"^<script>",
    "<style>" => r"^<style>",
    "</script>" => r"</script>$",
    "</style>" => r"</style>$",
)

function load_html_minify()
	try
        return pyimport("minify_html")
    catch
        pip("install", "minify-html")
        return pyimport("minify_html")
    end
end

function minify_website(;minify_css=true, minify_js=true, minify_kwargs...)
    site = path(:site)
    @assert isdir(site)
    # remove toplevel paths
    for path in readdir(site)
        f = joinpath(site, path)
        islink(f) && rm(f)
    end
    minify = load_html_minify().minify
    buf = IOBuffer()
    for file in walkfiles(site;
                          exts=Set((".html", ".css", ".js")),
                          subdir=true)
        if endswith(file, ".css")
            opentag = "<style>"
            closetag = "</style>"
        elseif endswith(file, ".js")
            opentag = "<script>"
            closetag = "</script>"
        else
            opentag = ""
            closetag = ""
        end
        println(buf, opentag, read(file, String), closetag)
        min = minify(;code=String(take!(buf)),
                     minify_css,
                     minify_js,
                     minify_kwargs...)
        min = replace(min, rx_tags[opentag] => "")
        min = replace(min, rx_tags[closetag] => "")
        write(file, min)
    end
    close(buf)
end

export minify_website

end


text:
struct RealDerivative{NT,V} <: AbstractDerivative
    fields::NT
    vectorised_data::V
end

RealDerivative(net::NeuralNetwork) = RealDerivative(out_type(net), net)
function RealDerivative(T::Type{<:Number}, net::NeuralNetwork)
    pars = trainable(net)

    vec    = similar(trainable_first(pars), T, _tlen(pars))
    i, fields = weight_tuple(net, vec)
    return RealDerivative(fields, (vec,))
end

# add :fields to the properties
@inline Base.propertynames(s::RealDerivative) = propertynames(getfield(s, :fields))
@inline Base.getindex(s::RealDerivative, val) = getproperty(s, val)
@inline Base.getproperty(s::RealDerivative, val::Symbol) = _getproperty(s, val)
@inline Base.getproperty(s::RealDerivative, val::Int) = _getproperty(s, val)
@inline function _getproperty(s::RealDerivative, val)
    val===:tuple_all_weights && return vec_data(s)
    return getproperty(getfield(s, :fields), val)
end
@inline vec_data(s::RealDerivative) = getfield(s, :vectorised_data)
@inline fields(s::RealDerivative) = getfield(s, :fields)

weights(der::RealDerivative) = der

Base.isapprox(x::RealDerivative, y::RealDerivative; kwargs...) =
    _isapprox(fields(x), fields(y); kwargs...)

Base.show(io::IO, der::RealDerivative) = begin
    pn = propertynames(der)
    str = "{"
    for fn=pn[1:end-1]
        str *= ":$fn, "
    end
    str *= ":" * string(last(pn))*" }"
    print(io,
    "RealDerivative with fields: ", str)
end


text:
using GeneralizedGenerated

binary_operators=[
    :(+),
    :(*),
    :(-),
    :(/),
    :(^),
]
unary_operators=[
    :(sin),
    :(cos),
    :(tan),
]

function get_unary_operator()
    return rand(unary_operators)
end
function get_binary_operator()
    return rand(binary_operators)
end

function create_expression( x::Symbol, p::Symbol, depth=0, pcount=1 )
    expr1=x
    if rand(0:depth)>0
        (expr1,pcount)=create_expression( x, p, depth-1, pcount )
    end
	expr2=:( p[$pcount] )
    if rand(0:depth) == 0
        pcount+=1
    else
        (expr2,pcount)=create_expression( x, p, depth-1, pcount )
    end
	expr=expr1
	if rand(1:2)==1
		op=get_unary_operator()
		expr=Expr( :call, op, expr)
	end
	op=get_binary_operator()
	expr=Expr( :call, op, expr, expr2)
    return expr,pcount
end

function crash()
	step=1
	exception_count=0
	nan_count=0
	θ = π/2
	p = [1,2]
	while true
		(e,pcount) = create_expression( :θ, :p, 1)
		f = mk_function( :( (θ,p) -> $e ) )
		try
			new_θ=f(θ,p)
			println(step," ",exception_count," ",nan_count," : ",e," ",p," ",θ,"->",new_θ)
			if isnan(new_θ)
				nan_count+=1
			end
		catch e
			println(step," ",e)
			exception_count+=1
		end
		step+=1
	end
end

crash()


text:
# # Machine Learning Problem : Housing Dataset
#
# The housing problem functions as a starting point in Machine Learning.
# We'll be demonstrating the use of Julia's [Flux Package](https://fluxml.ai/)
# to do this problem.
#
# The data replicates the housing data example from the Knet.jl readme. Although we
# could have reused more of Flux (see the mnist example), the library's
# abstractions are very lightweight and don't force you into any particular
# strategy.
#
# [This](http://www.mit.edu/~6.s085/notes/lecture3.pdf) might help you know more about the fundamentals of what
# we're about to do. If you don't understand something there which is also not mentioned here in this file,
# you may overlook that (or search it up on google to quench your curiosity :-)

using Flux.Tracker, Statistics, DelimitedFiles
using Flux.Tracker: Params, gradient, update!
using DelimitedFiles, Statistics
using Flux: gpu

# ## Getting the data and other pre-processing.
# We'll start by getting `housing.data` and splitting it into
# training and test sets.
# Training Dataset is the sample of data used to **fit** the model while
# Test Dataset is the sample of data used to provide an unbiased evaluation
# of a final model fit on the training dataset.

# Our aim is to predict the price of the house. In this dataset, the last
# feature is the price and would therefore be our target.

cd(@__DIR__)

isfile("housing.data") ||
  download("https://raw.githubusercontent.com/MikeInnes/notebooks/master/housing.data",
           "housing.data")

rawdata = readdlm("housing.data")'

#-

# Specifying the split ratio and **x** and **y**
split_ratio = 0.1

x = rawdata[1:13,:] |> gpu
y = rawdata[14:14,:] |> gpu

# ### Normalising
# What is the need ?
# Normalization is a technique often applied as part of data preparation for machine learning.
# The goal of normalization is to change the values of numeric columns in the dataset to a common scale,
# without distorting differences in the ranges of values. For machine learning, every dataset does not require normalization.
# It is required only when features have different ranges like in this case.

x = (x .- mean(x, dims = 2)) ./ std(x, dims = 2)

# ### Splitting into test and training sets.

split_index = floor(Int,size(x,2)*split_ratio)
x_train = x[:,1:split_index]
y_train = y[:,1:split_index]
x_test = x[:,split_index+1:size(x,2)]
y_test = y[:,split_index+1:size(x,2)]

# ## The Model
# Here comes everyone's favourite part : implementing a machine learning model.
#
# A ML model is in it's simplest terms a mathematical model which has a number of parameters
# that need to be learned from the data provided. The data has an important task: to fit our model parameters.
# The more data we have, the more we can accurately predict the target.
#
# Hyperparameters aren't learnt during the training process. They can be treated as constants that are fixed for the
# entire process. These parameters express important properties of the model such as its complexity or how fast it should learn.
#
# We'll now define the Weight (W) and the Bias (b) terms. They are our hyperparameter which
# we tune to enhance our predictions during gradient descent.
# To get an intution about how gradientDescent actually works, check out Andrew Ng's awesome explaination
# here. [Video 1: Intution](https://www.youtube.com/watch?v=rIVLE3condE) |
# [Video 2: The Algorithm](https://www.youtube.com/watch?v=yFPLyDwVifc)

W = param(randn(1,13)/10) |> gpu
b = param([0.]) |> gpu

# Here are our prediction and loss functions.
# - The prediction functions returns our prediction of the price of the house as
# suggested by our 2 hyperparameters: W and b.
# - MSE is the average of the squared error that is used as the loss function for least squares regression.
# It is defined as the sum, over all the data points, of the square of the difference between the predicted and actual target
# variables, divided by the number of data points.
#
# Loss functions evaluate how well your algorithm models your dataset.
# If predictions are off, the loss function is high. If they're good, it'll be low.

predict(x) = W*x .+ b
meansquarederror(ŷ, y) = sum((ŷ .- y).^2)/size(y, 2)
loss(x, y) = meansquarederror(predict(x), y)

# ### Gradient Descent
# Optimizing our parameters to get accurate prediction. Learn more from the links mentioned above.

η = 0.1
θ = Params([W, b])

for i = 1:10
  g = gradient(() -> loss(x_train, y_train), θ)
  for x in θ
    update!(x, -g[x]*η)
  end
  @show loss(x_train, y_train)
end

# ## Predictions
# Now we're in a position to know how well our program works on the given data.

err = meansquarederror(predict(x_test),y_test)
println(err)

# The prepared model might not very good for predicting the housing prices and may have high error.
# One can improve the prediction results using many other possible machine learning algorithms and techniques.
# If this was your first ML project in Flux, Congrats!
#
# You should have gotten a gist of basic ML functionality in Flux Package using Julia by now.

# ## References :
# 1. [Introduction to Loss Functions](https://algorithmia.com/blog/introduction-to-loss-functions)
# 2. [Why Data Normalization is necessary for Machine Learning models](https://medium.com/@urvashilluniya/why-data-normalization-is-necessary-for-machine-learning-models-681b65a05029)
# 3. [About Train, Validation and Test Sets in Machine Learning](https://towardsdatascience.com/train-validation-and-test-sets-72cb40cba9e7)
# 4. [How to select the Right Evaluation Metric for Machine Learning Models: Part 1 Regression Metrics](https://towardsdatascience.com/how-to-select-the-right-evaluation-metric-for-machine-learning-models-part-1-regrression-metrics-3606e25beae0)
# 5. [MIT's Notes on Linear Regression](http://www.mit.edu/~6.s085/notes/lecture3.pdf)
# 6. [ML | Hyperparameters: An Understanding](https://www.geeksforgeeks.org/ml-hyperparameter-tuning/)
#

max_stars_repo_name: PseudoCodeNerd/model-zoo

max_issues_repo_name: PseudoCodeNerd/model-zoo

max_forks_repo_name: PseudoCodeNerd/model-zoo


text:
using Layered
using Test
using FileIO
using Colors

generate_images = false
originalimagepath(str) = joinpath("original_images", str)
testimagepath(str) = joinpath("test_images", str)

if !isdir("test_images/")
    mkdir("test_images/")
end

function test_or_generate_png(surface, filename)
    if generate_images
        png(surface, originalimagepath(filename); px_per_pt = 1)
    else
        png(surface, testimagepath(filename); px_per_pt = 1)
        testimg = load(testimagepath(filename))
        originalimg = load(originalimagepath(filename))
        @test testimg == originalimg
    end
end

@testset "geometry" begin
    p1 = Point(0, 0)
    x1 = X(0)
    y1 = Y(0)
    @test p1 == x1
    @test p1 == y1

    p2 = P(3, 5.5)
    x2 = X(3)
    y2 = Y(5.5)
    @test p2 == x2 + y2

    p3 = P(-3, 5.5)
    @test p3 == y2 - x2

    p4 = P(1, 1)
    p5 = P(3, 5)
    v = p5 - p4
    @test v == P(2, 4)
    l = Line(p4, p5)
    @test vector(l) == v
    @test length(l) == magnitude(v) == sqrt(2 ^ 2 + 4 ^ 2)
    @test direction(l) == v ./ length(l)
    @test l * 3 == Line(P(3, 3), P(9, 15))
    @test l / 2 == Line(P(0.5, 0.5), P(1.5, 2.5))
end

@testset "image_lines_1" begin
    c, tl = canvas(100, 100)
    l = layer_in_rect!(tl, c.rect, :w; margin=5)
    line!(l, P(-1,  -1), P(1,  -1)) + Stroke("red") + Linewidth(2)
    line!(l, P(-1,  -0.8), P(1,  -0.8)) + Stroke("blue") + Linewidth(3)
    line!(l, P(-1,  -0.6), P(1,  -0.6)) + Stroke("green") + Linewidth(4)

    filename = "image_lines_1.png"
    test_or_generate_png(c, filename)
end

@testset "image_circle_1" begin
    c, tl = canvas(100, 100)
    circle!(tl, O, 25) + Fill("orange") + Stroke("black") + Linewidth(5)

    filename = "image_circle_1.png"
    test_or_generate_png(c, filename)
end

@testset "petals" begin

    c, tl = canvas(400, 400)

    l = layer_in_rect!(tl, c.rect, :w, margin=30)

    degrees = range(0, 360, length=21)[1:end-1]

    petals = path!.(ang -> begin
        ang
        endpoint = P(ang)
        segments = [Arc(O, endpoint, 0.2), Arc(endpoint, O, 0.2)]
        Path(segments)
    end, l, deg.(degrees)) .+ Stroke.(LCHuv.(70, 50, degrees)) .+ Linewidth(10) .+
        Fill.(LCHuvA.(70, 50, degrees, 0.5)) .+ Operator(:mult)

    circ = circle!(l, O, 0.25) + Visible(false)

    l + Clip(circ, c.rect)

    filename = "image_petals.png"
    test_or_generate_png(c, filename)
end

@testset "gradients" begin

    c, tl = canvas(100, 100)

    l = layer_in_rect!(tl, c.rect, :w) + Stroke(nothing)

    rect!(l, P(-0.5, -0.5), 1, 1, deg(0)) + Fill(Gradient(P(-1, -1), O, "red", "blue"))
    rect!(l, P(0.5, -0.5), 1, 1, deg(0)) + Fill(Gradient(P(1, -1), O, "yellow", "green", "orange"))
    rect!(l, P(-0.5, 0.5), 1, 1, deg(0)) + Fill(RadialGradient(O, 1, RGB(1, 0, 1), RGB(0, 1, 0)))
    rect!(l, P(0.5, 0.5), 1, 1, deg(0)) + Fill(RadialGradient(O, 0.3, 1.41, RGB(1, 0, 0), RGB(0, 1, 1)))

    filename = "image_gradients.png"
    test_or_generate_png(c, filename)

end

@testset "svg to path" begin
    c, tl = canvas(100, 100)

    l = layer_in_rect!(tl, c.rect, :w; margin=10)

    headsvg = open("head.svg") do file
        content = join(readlines(file))
        split(split(content, "<path d=\"")[2], "\"/>")[1]
    end

    headpath = centeredin(Path(headsvg), 2)

    path!(l, headpath) + Fill("black")

    filename = "image_headsvg_path.png"
    test_or_generate_png(c, filename)
end

@testset "gridrects" begin

    c, tl = canvas(100, 100)

    rects!(tl, c.rect) do r
        gridrects(r, 3, 4, 5, 5, 5)
    end + Fill(:frac => f -> LCHuv(60, 60, f * 270))

    filename = "image_gridrects.png"
    test_or_generate_png(c, filename)

end

@testset "linaxis" begin

    c, tl = canvas(100, 100)

    ax = LinAxis((0, 10), (1, -1), 0.1, 0.1)
    xx = 0:0.1:10
    yy = sin.(xx)

    path!(tl, c.rect) do r
        Path(ax.(r, xx, yy))
    end + Stroke("red")

    filename = "image_linaxis.png"
    test_or_generate_png(c, filename)
end

@testset "record" begin
    record("test.mp4", 0:3) do t
        c, l = canvas(50, 50)
        c
    end
end


text:
let
    text = raw"""
        # PlutoStaticHTML.jl

        To quickly build two notebooks (evaluate all the code and convert the output to HTML), use:

        ```
        julia> dir = joinpath("posts", "notebooks");

        julia> files = ["notebook1.jl", "notebook2.jl"];

        julia> parallel_build!(dir, files)
        ```

        and to build all notebooks, use:

        ```
        julia> parallel_build!(dir)
        ```

        See https://rikhuijzer.github.io/PlutoStaticHTML.jl/dev/ for the full documentation.
        """
    @doc text PlutoStaticHTML
end


text:
const λ  = 0.003 #in mm => 150GHz
const θᵦ = (π / 180) * 0.1#   8 primi
const ω₀ = λ / ( π * θᵦ ) # θᵦ is the beam opening angle in radians
const z_atm = 3E4 #depends on the observation site
const Tᵧ = 240 # ground temperature in kelvin
const χ₁_0 = 1 # in m^{-2}
const χ₂_0 = 1 # in m^{-2}
const z_0  = 1E3 # the quote where water vapor vanish
const L₀   = 597 # in m, the turbulance correlation length

max_stars_repo_name: gruffoc/AFO

max_issues_repo_name: gruffoc/AFO

max_forks_repo_name: gruffoc/AFO


text:
	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	#	Using librarys
	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	# GNU Scientific Library
	using GSL ;

	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	#	Special functions related
	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	# Spherical Bessel function first kind j_n(x)
	function Sbessj( i::Int64 , x )
		return sqrt(pi*(1/(2*x)))*besselj(i+1.0/2,x) ;
	end

	# Spherical Bessel function second kind n_n(x) [Neumann]
	function Sneum( i::Int64 ,x )
		sqrt(pi*(1/(2*x)))*bessely(i+1.0/2,x) ;
	end

	# Spherical Bessel function first kind derivative ( Morse & Feshbach definition )
	function Sbessj_p( i::Int64 , x )
 		return (1/(2*i+1))*(i*Sbessj(i-1,x) - (i+1)*Sbessj(i+1,x)) ;
	end

	# Spherical Bessel second kind ( Neumann ) derivative
	function Sneum_p( i::Int64 , x )
		return Sneum(i-1,x) - ((i+1)/x)*Sneum(i,x) ;
	end

	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	#	Liquid Sphere
	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	# Liquid sphere coefficient
	function Coef_LiquidSphere( n::Int64, x1::Float64 )
		return x1*Sbessj_p( n, x1 )/Sbessj( n, x1 ) ;
	end

	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	#	Prolate/Oblate auxiliar
	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	function Index( M::Int64, i::Int64, j::Int64 )
		if j < i # Domain errors
			return	display("ERROR: j < i");
		end
		if i > M + 1 || j > M + 1
			return	display("ERROR: M < i-1, j-1");
		end
		return map( Int64, j + ( M + 1 )*( i - 1 ) - i*( i - 1 )/2 );
	end

	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	#	File I/O
	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	function ReadFileToArrayBF( file::AbstractString, sep::Char, elem::Int )
	#	Read 'file' (as BigFloat) with values 'sep' separated.
	#	A newline in the 'file' means that the retrieved structure will be a matrix
	# 	The 'elem' parameter selects all the entries (elem = 0) or a particular value
	# 	(elem = number)
		if elem == 0 # Return entire vector or matrix
			return map( x -> parse(BigFloat,x), readdlm( file, sep ,ASCIIString) );
		else # Return only the specified element
			return map( x -> parse(BigFloat,x), [readdlm( file, sep ,ASCIIString)[elem]] )[ 1 ] ;
		end
	end

	function ReadFileToArrayBF2( file::AbstractString, sep::Char, row::Any, col::Any )
	#	Read 'file' (as BigFloat) with values 'sep' separated.
	# 	Returned value is a matrix in general.
		return map( x -> parse(BigFloat,x), collect( readdlm( file, sep ,ASCIIString)[row,col] ) ) ;
	end

	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	#	Eigenvalues (prolate)
	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	function pro_lambdamn_approx( c::Any, m::Int64, n::Int64 )
		# Transcription to Julia of the AGD (Adelman, Gumerov & Duraiswami) MATLAB code.
		# View ' pro_calculate_lambdamn_approx.m ' in the AGD source directory
		convert( Float64, c ) ;
		const N = m + n + 200 ; # Size of the matrix
		const A = zeros( Float64, N, N ) ;
		if ( mod( n - m, 2) == 0 )
			r = 0 ; # even
		else
			r = 1 ; # odd
		end
		# Fill in the A matrix
		for i = 1 : N
			if i == 1
				A[ 1, 1 ] = calculate_betar( c, m, r ) ;
				A[ 1, 2 ] = calculate_alphar( c, m, r ) ;
			elseif  i == N
				A[ N, N - 1 ] = calculate_gammar( c, m, r ) ;
				A[ N, N ] = calculate_betar( c, m, r ) ;
			else
				A[ i, i - 1 ] = calculate_gammar( c, m, r ) ;
				A[ i, i ] = calculate_betar( c, m, r ) ;
				A[ i, i + 1 ] = calculate_alphar( c, m, r ) ;
			end
			r = r + 2;
		end
		const d = eigvals!(A);
		if ( mod(n - m, 2) == 0 )
			const lambda_approx = real(d[ Int( (n - m + 2) / 2 ) ]) ;
		else
			const lambda_approx = real(d[ Int( (n - m + 1) / 2 ) ]) ;
		end
		filename = @ sprintf("data/pro_%08d_%03d_%03d_lambda_approx.txt", trunc( Int, c*1000 ), m, n ) ;
		writedlm( filename,lambda_approx );
	end

	function calculate_alphar( c::Float64, m::Int64, r::Int64 )
		return ((( 2 * m + r + 2 )*( 2 * m + r + 1 ))/
		(( 2 * m + 2 * r + 5 )*( 2 * m + 2 * r + 3 )))*( c*c );
	end

	function calculate_betar( c::Float64, m::Int64, r::Int64 )
		return ( m + r )*( m + r + 1 ) + (( 2*( m + r )*( m + r + 1 ) - 2*(m*m) - 1) /
		(( 2*m + 2*r - 1 )*( 2*m + 2 * r + 3 )))*(c*c);
	end

	function calculate_gammar( c::Float64, m::Int64, r::Int64 )
		return (( r * (r - 1) )/(( 2*m + 2*r - 3 )*( 2*m + 2*r - 1 )))*(c*c);
	end

	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	#	Eigenvalues (oblate)
	# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	function obl_lambdamn_approx( c::Any, m::Int64, n::Int64 )
		# Transcription to Julia of the AGD (Adelman, Gumerov & Duraiswami) MATLAB code.
		# View ' obl_calculate_lambdamn_approx.m ' in the AGD source directory
		convert( Float64, c ) ;
		const N = m + n + 200 ; # Size of the matrix
		const A = zeros( Float64, N, N ) ;
		if ( mod( n - m, 2) == 0 )
			r = 0 ; # even
		else
			r = 1 ; # odd
		end
		# Fill in the A matrix
		for i = 1 : N
			if i == 1
				A[ 1, 1 ] = calculate_betar_obl( c, m, r ) ;
				A[ 1, 2 ] = calculate_alphar_obl( c, m, r ) ;
			elseif  i == N
				A[ N, N - 1 ] = calculate_gammar_obl( c, m, r ) ;
				A[ N, N ] = calculate_betar_obl( c, m, r ) ;
			else
				A[ i, i - 1 ] = calculate_gammar_obl( c, m, r ) ;
				A[ i, i ] = calculate_betar_obl( c, m, r ) ;
				A[ i, i + 1 ] = calculate_alphar_obl( c, m, r ) ;
			end
			r = r + 2;
		end
		const d = eigvals!(A);
		if ( mod( n - m, 2 ) == 0 )
			const lambda_approx = real(d[ Int( (n - m + 2) / 2 ) ]) ;
		else
			const lambda_approx = real(d[ Int( (n - m + 1) / 2 ) ]) ;
		end
		filename = @ sprintf("data/obl_%08d_%03d_%03d_lambda_approx.txt", trunc( Int, c*1000 ), m, n ) ;
		writedlm( filename,lambda_approx );
	end

	function calculate_alphar_obl( c::Float64, m::Int64, r::Int64 )
		return ((( 2 * m + r + 2 )*( 2 * m + r + 1 ))/
		(( 2 * m + 2 * r + 5 )*( 2 * m + 2 * r + 3 )))*( -c*c );
	end

	function calculate_betar_obl( c::Float64, m::Int64, r::Int64 )
		return ( m + r )*( m + r + 1 ) + (( 2*( m + r )*( m + r + 1 ) - 2*(m*m) - 1) /
		(( 2*m + 2*r - 1 )*( 2*m + 2 * r + 3 )))*(-c*c);
	end

	function calculate_gammar_obl( c::Float64, m::Int64, r::Int64 )
		return (( r * (r - 1) )/(( 2*m + 2*r - 3 )*( 2*m + 2*r - 1 )))*(-c*c);
	end


text:
## Reflections and projections
#==========================================================================================#

"""
    reflect_xz(vector)
    reflect_xz(line :: Line)
    reflect_xz(horseshoe :: Horseshoe)

Reflect the ``y``-coordinate of a given 3-dimensional vector (or ``Line`` or ``Horseshoe``) about the ``x``-``z`` plane.
"""
reflect_xz(vector) = SVector(vector[1], -vector[2], vector[3])
reflect_xz(horseshoe :: Horseshoe) = (Horseshoe ∘ reflect_xz ∘ bound_leg)(horseshoe)

"""
    project_yz(vector)
    project_yz(line :: Line)

Project a given 3-dimensional vector or `Line` into the ``y``-``z`` plane.
"""
project_yz(vector) = SVector(0, vector[2], vector[3])

"""
    flip_xz(vector)

Reflect the ``x``- and ``z``- coordinates of a given 3-dimensional vector about the ``y``-``z`` and ``x``-``y`` planes respectively for the representation in body axes.
"""
flip_xz(vector) = SVector(-vector[1], vector[2], -vector[3])

## Axis transformations
#==========================================================================================#

abstract type AbstractAxisSystem end

struct Geometry <: AbstractAxisSystem end
struct Body      <: AbstractAxisSystem end
struct Stability <: AbstractAxisSystem end
struct Wind      <: AbstractAxisSystem end

Base.show(io :: IO, :: Geometry)  = print(io, "Geometry")
Base.show(io :: IO, :: Body)      = print(io, "Body")
Base.show(io :: IO, :: Stability) = print(io, "Stability")
Base.show(io :: IO, :: Wind)      = print(io, "Wind")

geometry_to_body_axes(coords) = flip_xz(coords)

"""
    geometry_to_stability_axes(coords, α)

Convert coordinates from geometry to stability axes with angle ``α``.
"""
geometry_to_stability_axes(coords, α :: T) where T <: Real = RotY{T}(α) * coords

"""
    geometry_to_stability_axes(coords, α)

Convert coordinates from stability to geometry axes with angle ``α``.
"""
stability_to_geometry_axes(coords, α :: T) where T <: Real = geometry_to_stability_axes(coords, -α)

"""
    geometry_to_wind_axes(coords, α, β)
    geometry_to_wind_axes(line :: Line, α, β)

Convert coordinates from geometry axes to wind axes for given angles of attack ``α`` and sideslip ``\\beta.``
"""
geometry_to_wind_axes(coords, α, β) = let T = promote_type(eltype(α), eltype(β)); RotZY{T}(β, α) * coords end

function geometry_to_wind_axes(horseshoe :: Horseshoe, α, β)
    T = promote_type(eltype(α), eltype(β))
    transform(horseshoe, LinearMap(RotZY{T}(β, α)))
end

geometry_to_wind_axes(coords, fs :: Freestream) = geometry_to_wind_axes(coords, fs.alpha, fs.beta)
geometry_to_wind_axes(horseshoe :: Horseshoe, fs :: Freestream) = geometry_to_wind_axes(horseshoe, fs.alpha, fs.beta)

"""
    wind_to_geometry_axes(coords, α, β)

Convert coordinates from wind axes to geometry axes for given angles of attack ``α`` and sideslip \\beta.``
"""
wind_to_geometry_axes(coords, α, β) = let T = promote_type(eltype(α), eltype(β)); RotYZ{T}(-α, -β) * coords end

function wind_to_geometry_axes(horseshoe :: Horseshoe, α, β)
    T = promote_type(eltype(α), eltype(β))
    transform(horseshoe, LinearMap(RotZY{T}(-α, -β)))
end

# function rotate_zy(θ₁, θ₂)
#     sinθ₁, cosθ₁ = sincos(θ₁)
#     sinθ₂, cosθ₂ = sincos(θ₂)
#     z            = zero(sinθ₁)

#     @SMatrix [ cosθ₁* cosθ₂  -sinθ₁  cosθ₁ * sinθ₂
#                sinθ₁* cosθ₂   cosθ₁  sinθ₁ * sinθ₂
#                  -sinθ₂         z        cosθ₂     ]
# end


text:
export constraint2coarse

"""
Adapts constraint definitions for the original grid, to coarser grids.
Some of these 'rules' are empirical or user defined. Some constraints are not adapted as
they are often considered invariant under subsampling the grid for imaging applications.
"""
function constraint2coarse(constraint,comp_grid,coarsening_factor)

  nr_constraints = length(constraint)

  #bound constraints: same as on fine grid

  #rank: same as on fine grid (limit by minimum dimension)
  for i=1:nr_constraints
    if constraint[i].set_type == "rank"
      constraint[i].max = min(constraint[i].max,minimum(comp_grid.n))
    end
  end

  #cardinality in a transform-domain: same as on fine grid (limit by number of elements)
  for i=1:nr_constraints
    if constraint[i].set_type == "cardinality"
      constraint[i].max = min(constraint[i].max,prod(comp_grid.n))
    end
  end

  #define below what to do with norms (l1, l2 and nuclear)
  # will probably depend on constraints and interpolation type that is used
  # below we use simple heuristics, better methods may improve performance of the multilevel scheme

if length(comp_grid.n)==3 && comp_grid.n[3]>1 #use 3D

  # #for point-wise bound constraints in a transform-domain on a 3D grid if the transform-domain is a derivative operator:
  # for i=1:3
  #   if haskey(constraint,string("use_TD_bounds_",i)) && (constraint[string("use_TD_bounds_",i)]==true)
  #     if (constraint[string("TDB_operator_",i)] in ["D_x","D_y","D_z","DFT"])
  #       constraint[string("TD_LB_",i)]=constraint[string("TD_LB_",i)].*coarsening_factor
  #       constraint[string("TD_UB_",i)]=constraint[string("TD_UB_",i)].*coarsening_factor
  #     end
  #   end
  # end

  #for l1 norm in a transform-domain: on a 3D grid: ||.||_1(coarse)=||.||_1(fine)/(coarsening factor^3)
  for i=1:nr_constraints
    if constraint[i].set_type == "l1"
      constraint[i].max = constraint[i].max/(coarsening_factor^3)
    end
  end

  #for l2 norm in a transform-domain: on a 3D grid: ||.||_2(coarse)=||.||_2(fine)/sqrt(coarsening factor^3)
  for i=1:nr_constraints
    if constraint[i].set_type == "l2"
      constraint[i].max = constraint[i].max/(sqrt(coarsening_factor^3))
    end
  end

else #use 2D

  # #for point-wise bound constraints in a transform-domain on a 2D grid if the transform-domain is a derivative operator:
  # for i=1:3
  #   if haskey(constraint,string("use_TD_bounds_",i)) && (constraint[string("use_TD_bounds_",i)]==true)
  #     if (constraint[string("TDB_operator_",i)] in ["D_x","D_y","D_z","DFT"])
  #       constraint[string("TD_LB_",i)]=constraint[string("TD_LB_",i)].*coarsening_factor
  #       constraint[string("TD_UB_",i)]=constraint[string("TD_UB_",i)].*coarsening_factor
  #     end
  #   end
  # end

    #for l1 norm in a transform-domain: on a 2D grid: ||.||_1(coarse)=||.||_1(fine)/(coarsening factor^2)
    for i=1:nr_constraints
      if constraint[i].set_type == "l1"
        constraint[i].max = constraint[i].max/(coarsening_factor^2)
      end
    end

    #for l2 norm in a transform-domain: on a 2D grid: ||.||_2(coarse)=||.||_2(fine)/(coarsening factor)
    for i=1:nr_constraints
      if constraint[i].set_type == "l2"
      constraint[i].max = constraint[i].max/(coarsening_factor)
      end
    end

    #nuclear norm:
    for i=1:nr_constraints
      if constraint[i].set_type == "nuclear"
        constraint[i].max = constraint[i].max/2.7
      end
    end

end #END if 2D or 3D

  # #subspace constraint
  # if constraint[i].set_type == "subspace"
  #   error("still need to define how to map subspace to a coarser grid in: function constraint2coarse")
  #   #maybe just use coarsening of matrix to start simple
  # end

return constraint
end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
# Input
print("Enter a binary number (up to 8 binary digits): ")
binary = readline()
# Validation
is_valid = match(r"[01]{1,8}", binary).match == binary
if !is_valid
    error("Invalid binary format given. Please check the constraints!")
end
# Parsing
decimal = parse(Int, binary, base = 2)
# Output
println("The decimal number is $decimal.")

max_stars_repo_name: feluelle/get-to-know-julia

max_issues_repo_name: feluelle/get-to-know-julia

max_forks_repo_name: feluelle/get-to-know-julia


text:
function rm_count(deck)
    r = 0
    m = 0
    for elem in deck
        if elem.prints[1].rarity == 3
            r += elem.amount
        elseif elem.prints[1].rarity == 4
            m += elem.amount
        end
    end
    (r = r, m = m)
end

max_stars_repo_name: asieev/GP-Lille

max_issues_repo_name: asieev/GP-Lille

max_forks_repo_name: asieev/GP-Lille


text:
using Test
include("pythagorean-triplet.jl")

@testset "triplets whose sum is 12" begin
    @test pythagorean_triplets(12) == [(3, 4, 5)]
end

@testset "triplets whose sum is 108" begin
    @test pythagorean_triplets(108) == [(27, 36, 45)]
end

@testset "triplets whose sum is 1000" begin
    @test pythagorean_triplets(1000) == [(200, 375, 425)]
end

@testset "triplets whose sum is 1001" begin
    @test pythagorean_triplets(1001) == []
end

@testset "returns all matching triplets" begin
    @test pythagorean_triplets(90) == [(9, 40, 41), (15, 36, 39)]
end

@testset "several matching triplets" begin
    @test pythagorean_triplets(840) == [
        (40, 399, 401),
        (56, 390, 394),
        (105, 360, 375),
        (120, 350, 370),
        (140, 336, 364),
        (168, 315, 357),
        (210, 280, 350),
        (240, 252, 348),
    ]
end

@testset "triplets for large number" begin
    @test pythagorean_triplets(30000) == [
        (1200, 14375, 14425),
        (1875, 14000, 14125),
        (5000, 12000, 13000),
        (6000, 11250, 12750),
        (7500, 10000, 12500),
    ]
end


text:
# This file contains some implementations of prox-operators

# L1 - PROX - OPERATOR
# Param : vector to
function L1prox(x, lambda)
    prox = zeros(size(x,1))
    for index in 1:size(x,1)
        if x[index] < -lambda
            prox[index] = x[index] + lamda
        elseif x[index] > lambda
            prox[index] = x[index] - lambda
        else
            prox[index] = 0
        end
    end
    return prox
end

# L2 - PROX - OPERATOR
# Param :
#   - x = vector to apply proxity
#   - λ = regularization parameter
function L2prox(x,lambda)
    # create empty prox vector
    prox = zeros(size(x,1))
    l2norm = norm(x,2)
    if l2norm > lambda
        return (1-lambda/l2norm)*x
    else
        return zeros(size(x,1))
    end
end

max_stars_repo_name: zagazao/num-opt

max_issues_repo_name: zagazao/num-opt

max_forks_repo_name: zagazao/num-opt


text:
using Random
using MultiAgentPathFinding
using MultiAgentAllocationTransit

## Hard-code a bunch of depots and many more sites
const DEPOT1 = LatLonCoords((lat = 37.762892, lon = -122.472193))
const DEPOT2 = LatLonCoords((lat = 37.751751, lon = -122.410654))
const DEPOT3 = LatLonCoords((lat = 37.718779, lon = -122.462401))

const SITE1 = LatLonCoords((lat = 37.789290, lon = -122.426797)) # Top right
const SITE2 = LatLonCoords((lat = 37.739611, lon = -122.492203)) # Way left
const SITE3 = LatLonCoords((lat = 37.780593, lon = -122.434555)) # Under SITE 1
const SITE4 = LatLonCoords((lat = 37.739011, lon = -122.430059)) # Near the bottom
const SITE5 = LatLonCoords((lat = 37.749018, lon = -122.462966)) # 2nd most left

depots = [DEPOT1, DEPOT2, DEPOT3]
sites = [SITE1, SITE2, SITE3, SITE4, SITE5]
depot_sites = vcat(depots, sites)

n_depots = 3
n_agents = 5
n_sites = 5

cost_fn(i, j) = MultiAgentAllocationTransit.distance_lat_lon_euclidean(depot_sites[i], depot_sites[j])
agent_tours = task_allocation(n_depots, n_sites, n_agents, depot_sites, cost_fn)
agent_tasks = get_agent_task_set(agent_tours, n_depots, n_sites)


text:
#######################################################################
# Analysing Julia files and packages for lower bounds
###############################
#
# Functions for gathering lower-bound statistics
# for individual files, packages, and lists of packages
#######################################################################

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Analysis of files
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# Checks if `filePath` is a Julia file
isJuliaFile(filePath :: String) :: Bool = endswith(filePath, ".jl")

# Assuming that `pkgPath` is a prefix of `filePath`,
# carves out the path to the file within ``pkgPath``
filePathWithinPkg(filePath::String, pkgPath::String) :: String =
    filePath[length(pkgPath)+2:end]

# Reads file `filePath` located in `pkgPath`
# and updates `pkgStat` accordingly with the file's lower-bounds stats
# - keeps track of failed and interesting files
processFile(filePath::String, pkgPath::String, pkgStat::PackageStat) = begin
    fileInfo = nothing
    try
        # file might not be accessible
        fileInfo = lbFileInfo(filePath; isPath=true)
    catch e
        isa(e, Base.IOError) ? @debugonly(@warn e) : @error e
    end
    if fileInfo === nothing || fileInfo.err !== nothing
        pkgStat.failedFilesNum += 1
    elseif nonVacuous(fileInfo.lbStat)
        pkgStat.nonVacFilesNum += 1
        # cut pkgPath from file name for readability
        pkgStat.filesInfo[filePathWithinPkg(filePath, pkgPath)] = fileInfo
    end
end

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Analysis of a package
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# Walks `files` in `root`
# and adds their lower-bounds statistics to `pkgStat`
processFilesInDir(root, files, pkgsPath::String, pkgStat::PackageStat) = begin
    # we are only interested in Julia files
    files = filter(isJuliaFile, files)
    pkgStat.totalFilesNum += length(files)
    map(file -> processFile(joinpath(root, file), pkgsPath, pkgStat),
        files)
end

# Returns cumulative lower-bound statistics
# for all files in `filesInfo`
cumulativeLBStat(filesInfo :: FilesLBInfo) :: LBStat = begin
    totalLBFreq = LBValsFreq()
    for fInfo in values(filesInfo)
        fInfo.lbStat === nothing ||
            unionMergeWith!(totalLBFreq, fInfo.lbStat.lbsFreq)
    end
    LBStat(totalLBFreq)
end

# Walks `src` directory of package `pkgName` located at `pkgPath`
# and collects its lower-bounds statistics
processPkg(pkgPath :: String, pkgName :: String) :: PackageStat = begin
    srcPath = joinpath(pkgPath, "src")
    # we assume that correct Julia packages have `src` folder,
    # otherwise we are not interested
    pkgStat = PackageStat(pkgName, isdir(srcPath))
    pkgStat.hasSrc || return pkgStat
    # recursively walk all files in [src]
    for (root, _, files) in walkdir(srcPath)
        processFilesInDir(root, files, pkgPath, pkgStat)
    end
    # package summary statistics
    if pkgStat.nonVacFilesNum > 0
        pkgStat.lbStat = cumulativeLBStat(pkgStat.filesInfo)
    end
    pkgStat
end

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Analysis of packages
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pkgHasSrc(pkgStat :: PackageStat) :: Bool = pkgStat.hasSrc

# Returns cumulative lower-bound statistics
# for all packages in `pkgsStat`
cumulativeLBStat(pkgsStat :: Vector{PackageStat}) :: LBStat = begin
    totalLBFreq = LBValsFreq()
    for pkgStat in pkgsStat
        unionMergeWith!(totalLBFreq, pkgStat.lbStat.lbsFreq)
    end
    LBStat(totalLBFreq)
end

# String → (Vector{PackageStat}, Vector{PackageStat}, LBStat)
# Processes every folder in `pkgsPath` as a package folder
#   and computes lower-bounds statistics for it.
# Returns failed packages, successfully processed packages with their stats,
#   and cumulative statstics of good packages.
processPkgsDir(pkgsPath :: String) = begin
    pathsAndNames = subdirPathsWithNames(pkgsPath)
    processPkgPathAndName(pn) = begin
        @statusb "Processing package $(pn[2])"
        processPkg(pn[1], pn[2])
    end

    pkgsNum = length(pathsAndNames)
    pkgsStats = PackageStat[]
    if VERBOSE
        pkgsStats = Vector{PackageStat}(undef, pkgsNum)
        for pkgi in 1:pkgsNum
            pkgsStats[pkgi] = processPkgPathAndName(pathsAndNames[pkgi])
            pkgi % PKGS_NUM_STEP == 0 && @info "$pkgi PKGS PROCESSED"
        end
    else
        pkgsStats = map(processPkgPathAndName, pathsAndNames)
    end

    # sort packages information from most interesting to less interesting
    goodPkgs = filter(pkgHasSrc, pkgsStats)
    badPkgs  = filter(!pkgHasSrc, pkgsStats)
    (badPkgs, goodPkgs, cumulativeLBStat(goodPkgs))
end

max_stars_repo_name: prl-julia/julia-sub

max_issues_repo_name: prl-julia/julia-sub

max_forks_repo_name: prl-julia/julia-sub


text:
#=

  Arch Friends puzzle (Dells Logic Puzzles) in Julia ConstraintSolver.jl

  From http://brownbuffalo.sourceforge.net/ArchFriendsClues.html
  """
  Title: Arch Friends
  Author: Mark T. Zegarelli
  Publication: Dell Logic Puzzles
  Issue: April, 1998
  Page: 7
  Stars: 1

  Harriet, upon returning from the mall, is happily describing her four
  shoe purchases to her friend Aurora. Aurora just loves the four
  different kinds of shoes that Harriet bought
     (ecru espadrilles,fuchsia flats, purple pumps, and suede sandals),
  but Harriet can't recall at which different store
     (Foot Farm, Heels in a Handcart, The Shoe Palace, or Tootsies)
  she got each pair. Can you help these two figure out the order in
  which Harriet bought each pair of shoes, and where she bought each?

  1. Harriet bought fuchsia flats at Heels in a Handcart.
  2. The store she visited just after buying her purple pumps was not
     Tootsies.
  3. The Foot Farm was Harriet's second stop.
  4. Two stops after leaving The Shoe Place, Harriet bought her suede
     sandals.

  Determine: Order - Shoes - Store
  """

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function arch_friends(print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))
    n = 4

    @variable(model, 1 <= shoes[1:n] <= n, Int)
    EcruEspadrilles, FuchsiaFlats, PurplePumps, SuedeSandals = shoes

    @variable(model, 1 <= stores[1:n] <= n, Int)
    FootFarm, HeelsInAHandcart, TheShoePalace, Tootsies = stores

    @constraint(model, shoes in CS.AllDifferentSet())
    @constraint(model, stores in CS.AllDifferentSet())

    # 1. Harriet bought fuchsia flats at Heels in a Handcart.
    @constraint(model, FuchsiaFlats == HeelsInAHandcart)

    # 2. The store she visited just after buying her purple pumps was not
    #    Tootsies.
    @constraint(model, PurplePumps + 1 != Tootsies)

    # 3. The Foot Farm was Harriet's second stop.
    @constraint(model, FootFarm == 2)

    # 4. Two stops after leaving The Shoe Place, Harriet bought her suede
    # sandals.
    @constraint(model, TheShoePalace + 2 == SuedeSandals)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                shoes_val = convert.(Integer,JuMP.value.(shoes; result=sol))
                stores_val = convert.(Integer,JuMP.value.(stores; result=sol))
                println("shoes:$shoes_val")
                println("stores:$stores_val")
                println()

            end
        end
    else
        println("status:$status")
    end

    return status
end

@time arch_friends()

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
using RelocatableFolders, Test

module M

using RelocatableFolders

const DIR = folder"path"

end

@testset "RelocatableFolders" begin
    tests = function ()
        @test isdir(M.DIR)
        @test sort(readdir(M.DIR)) == ["file.jl", "subfolder", "text.md"]
        @test readdir(joinpath(M.DIR, "subfolder")) == ["other.jl"]

        @test read(joinpath(M.DIR, "file.jl"), String) == "# file.jl"
        @test read(joinpath(M.DIR, "text.md"), String) == "text.md"
        @test read(joinpath(M.DIR, "subfolder", "other.jl"), String) == "# other.jl"

        @test length(M.DIR.files) == 3
        @test M.DIR.mod == M
        @test M.DIR.path == joinpath(@__DIR__, "path")
    end
    from, to = joinpath.(Ref(@__DIR__), ("path", "moved"))
    try
        tests()
        @test String(M.DIR) == joinpath(@__DIR__, "path")

        # Remove the referenced folder `DIR`.
        @test isdir(from)
        @test !isdir(to)
        mv(from, to)
        @test isdir(to)
        @test !isdir(from)

        let path = String(M.DIR)
            rm(path; recursive = true)
            @test !isdir(path)
        end
        tests()
        @test String(M.DIR) != joinpath(@__DIR__, "path")

        # Return the referenced folder `DIR`.
        @test !isdir(from)
        @test isdir(to)
        mv(to, from)
        @test isdir(from)
        @test !isdir(to)

        tests()
        @test String(M.DIR) == joinpath(@__DIR__, "path")
    finally
        isdir(from) || mv(to, from)
    end
end


text:
l_bottom_vertex = 2
l_top_vertex = 5
r_top_vertex = 7
r_bottom_vertex = 11

mf = TrapezoidalMF(l_bottom_vertex, l_top_vertex, r_top_vertex, r_bottom_vertex)

# Evaluation tests
@assert Fuzzy.eval(mf, l_top_vertex) == Fuzzy.eval(mf, r_top_vertex) == 1
@assert Fuzzy.eval(mf, l_bottom_vertex) == Fuzzy.eval(mf, r_bottom_vertex) == 0
@assert Fuzzy.eval(mf, (l_bottom_vertex + l_top_vertex) / 2) == Fuzzy.eval(mf, (r_top_vertex + r_bottom_vertex) / 2) == 0.5

# Mean finding tests
@assert Fuzzy.mean_at(mf, 1) == (l_top_vertex + r_top_vertex) / 2
@assert Fuzzy.mean_at(mf, 0) == (l_bottom_vertex + r_bottom_vertex) / 2


text:
# this code currently only uses one custom type

mutable struct pfracinfo{E}
  xalphas           ::Vector{E}                 # [x_i - alpha_i for i = 1..r]
  betas             ::Vector{Vector{E}}         # matrix of beta evaluations
  inv_prod_betas1   ::Vector{E}                 # info for univariate solver
  prod_betas_coeffs ::Vector{Vector{Vector{E}}} # matrix of taylor coeffs
  deltas            ::Vector{Vector{E}}         # answers
  delta_coeffs      ::Vector{Vector{Vector{E}}} # taylor coeffs of answer
end

########### various helpers for expansions in terms of (x - alpha) ############

# evaluation at x = alpha should use divrem. Don't use evaluate for this!!!

# set t to the coefficients of q when expanded in powers of (x - alpha)
function taylor_get_coeffs!(t::Vector{E}, q::E, xalpha::E) where E
  R = parent(xalpha)
  empty!(t)
  while !iszero(q)
    q, r = divrem(q, xalpha)
    push!(t, r)
  end
end

function taylor_get_coeffs(q::E, xalpha::E) where E
  t = E[]
  taylor_get_coeffs(t, q, xalpha)
  return t
end

# opposite of taylor_get_coeffs
function taylor_from_coeffs(B::Vector{E}, xalpha::E) where E
  R = parent(xalpha)
  a = R(0)
  for i in length(B):-1:1
    a = a*xalpha + B[i]
  end
  return a
end

function taylor_set_coeff!(t::Vector{E}, i::Int, c::E) where E
  R = parent(c)
  while 1 + i > length(t)
    push!(t, R(0))
  end
  t[1 + i] = c
  while length(t) > 0 && iszero(t[end])
    popback!(t)
  end
end

#################### basic manipulation of Fac ################################

# a *= b^e
function mulpow!(a::Fac{T}, b::T, e::Int) where T
  @assert e >= 0
  if e == 0
    return
  end
  if isconstant(b)
    a.unit *= b^e
  elseif haskey(a.fac, b)
    a.fac[b] += e
  else
    a.fac[b] = e
  end
end

function mulpow!(a::Fac{T}, b::Fac{T}, e::Int) where T
  @assert !(a === b)
  @assert e >= 0
  if e == 0
    return
  end
  a.unit *= b.unit^e
  if isempty(a.fac) && e == 1
    a.fac = b.fac
  else
    for i in b.fac
      mulpow!(a, i[1], i[2]*e)
    end
  end
end

################ basic manipulation of polys #################################

function to_univar(a, v, Kx)
  return evaluate(a, [i == v ? gen(Kx) : Kx(0) for i in 1:nvars(parent(a))])
end

function from_univar(a, v, Kxyz)
  return evaluate(a, gen(Kxyz, v))
end

function gcdcofactors(a, b)
  g = gcd(a, b)
  if iszero(g)
    @assert iszero(a)
    @assert iszero(b)
    return (g, a, b)
  elseif isone(g)
    return (g, a, b)
  else
    return (g, divexact(a, g), divexact(b, g))
  end
end

# remove the content with respect to variable v
function primitive_part(a::E, v::Int) where E
  R = parent(a)
  d = degree(a, v)
  g = R(0)
  for i in 0:d
    ci = coeff(a, Int[v], Int[i])
    if !iszero(ci)
      g = gcd(g, ci)::E
    end
  end
  if iszero(g)
    return a
  elseif isone(g)
    return a
  else
    return divexact(a, g)
  end
end

function get_lc(a::E, v::Int) where E
  d = degree(a, v)
  @assert d >= 0
  return coeff(a, Int[v], Int[d])
end

function set_lc(a::E, v::Int, b::E) where E
  R = parent(a)
  d = degree(a, v)
  @assert d >= 0
  diff = (b - coeff(a, Int[v], Int[d]))
  if iszero(diff)
    return a
  else
    return a + diff*gen(R, v)^d
  end
end

function make_monic(a::E) where E
  if length(a) < 1
    return a
  end
  lc = coeff(a, 1)
  if isone(lc)
    return a
  else
    return inv(lc) * a
  end
end

function eval_one(a::E, v, alpha) where E
  R = parent(a)
  return divrem(a, gen(R, v) - alpha)[2]
end

################# generic partial fractions ###################################

#=
  R is a multivariate ring. We will be dealing with polynomials in l + 1
  variables, i.e. only the variables X, x_1, ..., x_l appear,
  where X = gen(R, mainvar) and x_i = gen(R, minorvar[i]).

  For fixed B[1], ..., B[r], precompute info for solving

    t/(B[1]*...*B[r]) = delta[1]/B[1] + ... + delta[r]/B[r] mod (x_i - alpha_i)^(deg_i + 1)

  for delta[1], ..., delta[r] given t in X, x_1, ..., x_l.

  alpha is an array of evaluation points
    x_1 = alpha[1]
    ...
    x_l = alpha[l]

  If, after evaluation at all x_j = alpha[j], the B[i] did not drop degree
  in X and are pairwise relatively prime in K[X], then the function returns
  true and a suitable pfracinfo. Otherwise, it returns false and junk

  Upon success, we are ready to solve for the delta[i] using pfrac.
  The only divisions that pfrac does in K are by the leading coefficients
  of the inv_prod_betas1 member of pfracinfo
=#

function pfracinit(
  B::Vector{E},
  mainvar,
  minorvars::Vector{Int},
  alphas::Vector
) where E

  r = length(B)       # number of factors
  l = length(minorvars)   # number of evaluated variables
  @assert r > 1
  @assert length(alphas) >= l

  R = parent(B[1])
  K = base_ring(R)

  # betas is a 1+l by r matrix with the original B[i] in the last row
  # and successive evaluations in the previous rows
  betas = [E[R(0) for i in 1:r] for j in 0:l]

  # prod_betas is like betas, but has B[1]*...*B[r]/B[i] in the i^th column
  # we don't need prod_betas, but rather
  #   prod_betas_coeffs = taylor of prod_betas in x_(1+j) - alpha[j]
  prod_betas_coeffs = [[E[] for i in 1:r] for j in 0:l]

  # inv_prod_betas1[i] will hold invmod(B[1]*...*B[r]/B[i], B[i]) in R[X]
  # after all x_j = alpha[j] have been evaluated away
  inv_prod_betas1 = E[R(0) for i in 1:r]

  # working space that need not be reallocated on each recursive call
  deltas = [E[R(0) for i in 1:r] for j in 0:l]
  delta_coeffs = [[E[] for i in 1:r] for j in 0:l]

  xalphas = E[gen(R, minorvars[j]) - R(alphas[j]) for j in 1:l]

  I = pfracinfo{E}(xalphas, betas, inv_prod_betas1, prod_betas_coeffs,
                                                          deltas, delta_coeffs)
  for i in 1:r
    I.betas[1+l][i] = deepcopy(B[i])
  end
  for j in l:-1:1
    for i in 1:r
      _, I.betas[j][i] = divrem(I.betas[j + 1][i], xalphas[j])
    end
  end

  for i in 1:r
    if degree(I.betas[1][i], mainvar) != degree(I.betas[l + 1][i], mainvar)
      # degree drop in univariates in X
      return false, I
    end
  end

  # univariate ring for gcdx
  S, x = PolynomialRing(K, "x")
  sub = [k == mainvar ? x : S(1) for k in 1:nvars(R)]

  for j in 0:l
    for i in 1:r
      p = R(1)
      for k in 1:r
        if k != i
          p *= I.betas[1 + j][k]
        end
      end
      if j > 0
        taylor_get_coeffs!(I.prod_betas_coeffs[1 + j][i], p, I.xalphas[j])
      else
        s = evaluate(betas[1][i], sub)
        t = evaluate(p, sub)
        g, s1, t1 = gcdx(s, t)
        if !isconstant(g)
          # univariates are not pairwise prime
          return false, I
        end
        I.inv_prod_betas1[i] = evaluate(divexact(t1, g), gen(R, mainvar))
        # leave prod_betas_coeffs[1][i] as zero, it doesn't matter
      end
    end
  end

  return true, I
end

#=
  Try to solve

    t/(B[1]*...*B[r]) = delta[1]/B[1] + ... + delta[r]/B[r] mod I

  where I = <(x_i - alpha_i)^(degs[i] + 1)>, and
      the x_i and alpha_i are those parameters of pfracinit, and
      the degs[i] are parameters of this function.

  With check = false, a solution to the above problem is generated without fail.
  With check = true, the solutions delta[j] are limited by

    deg_{x_i} delta[j] <= degs[i] - deg_{x_i}(B[1]*...*B[r]/B[j])

  The second return of pfrac (the array delta) is owned by I.
  So, at least its length should not be mutated.
=#
function pfrac(I::pfracinfo{E}, t::E, degs::Vector{Int}, check::Bool) where E
  return pfrac_recursive(I, t, degs, length(I.xalphas), check)
end

function pfrac_recursive(
  I::pfracinfo{E},
  t::E,
  degs::Vector{Int},
  lev::Int,
  check::Bool
) where E

  @assert 0 <= lev <= length(I.xalphas)
  @assert lev <= length(degs)

  r = length(I.inv_prod_betas1)

  deltas = I.deltas[1 + lev]

  if lev == 0
    for i in 1:r
      Q, deltas[i] = divrem(t*I.inv_prod_betas1[i], I.betas[1][i])
      # TODO check the denominator of deltas[i] and possibly return false
    end
    return true, deltas
  end

  newdeltas = I.deltas[lev]
  delta_coeffs = I.delta_coeffs[1 + lev]
  pbetas = I.prod_betas_coeffs[1 + lev]

  for i in 1:r
    empty!(delta_coeffs[i])
  end

  for l in 0:degs[lev]
    t, newt = divrem(t, I.xalphas[lev])
    for i in 1:r
      for s in 0:l-1
        if s < length(delta_coeffs[i]) && l - s < length(pbetas[i])
          newt -= delta_coeffs[i][1 + s] * pbetas[i][1 + l - s]
        end
      end
    end

    if iszero(newt)
      continue
    end

    ok, newdeltas = pfrac_recursive(I, newt, degs, lev - 1, check)
    if !ok
      return false, deltas
    end

    for i in 1:r
      if iszero(newdeltas[i])
        continue
      end
      if check && l + length(pbetas[i]) - 1 > degs[lev]
        return false, deltas
      end
      taylor_set_coeff!(delta_coeffs[i], l, newdeltas[i])
    end
  end

  for i in 1:r
    deltas[i] = taylor_from_coeffs(delta_coeffs[i], I.xalphas[lev])
  end

  return true, deltas
end

################# generic hensel lifting ######################################

#=
  Try to lift a factorization of A in 1 variable to 1 + n variables without
  any leading coefficient information.
=#
function hlift_without_lcc(
  A::E,
  Auf::Vector{E},       # univariate factorization in mainvar
  mainvar::Int,
  minorvars::Vector{Int},   # length n
  alphas::Vector
) where E

  R = parent(A)
  n = length(minorvars)
  r = length(Auf)
  @assert n > 0
  @assert r > 0

  if r < 2
    return true, [primitive_part(A, mainvar)]
  end

  lcA = get_lc(A, mainvar)
  A *= lcA^(r-1)

  degs = [degree(A, i) for i in minorvars]

  evals = E[R(0) for i in 1:n+1]
  lcevals = E[R(0) for i in 1:n+1]

  evals[n + 1] = A
  lcevals[n + 1] = lcA
  for i in n:-1:1
    evals[i] = eval_one(evals[i + 1], minorvars[i], alphas[i])
    lcevals[i] = eval_one(lcevals[i + 1], minorvars[i], alphas[i])
  end

  fac = [divexact(lcevals[1], get_lc(Auf[i], mainvar))*Auf[i] for i in 1:r]

  for k in 1:n
    ok, fac = hliftstep(fac, mainvar, minorvars[1:k], degs, alphas, evals[k+1], true)
    if !ok
      return false, fac
    end
  end

  return true, [make_monic(primitive_part(fac[i], mainvar)) for i in 1:r]
end

#=
  Try to lift a factorization of A in 1 variable to 1 + n variables with the
  assumption that divs[i] is a divisor of the leading coefficient of
  the i^th liftedfactor.
=#
function hlift_with_lcc(
  A::E,
  Auf::Vector{E},       # univariate factorization in mainvar
  divs::Vector{E},      # lead coeff divisors
  mainvar::Int,
  minorvars::Vector{Int},   # length n
  alphas::Vector,
) where E

  R = parent(A)
  K = base_ring(R)
  n = length(minorvars)
  r = length(Auf)
  @assert n > 0
  @assert r > 1

  fac = copy(Auf)

  m = get_lc(A, mainvar)
  for i in divs
    ok, m = divides(m, i)
    if !ok
      # lcc was definitely wrong
      return false, fac
    end
  end

  lcs = zeros(R, n + 1, r)
  for j in 1:r
    lcs[n + 1, j] = divs[j]*m
    for i in n:-1:1
      lcs[i, j] = eval_one(lcs[i + 1, j], minorvars[i], alphas[i])
    end
  end

  evals = zeros(R, n + 1)
  evals[n + 1] = A*m^(r - 1)
  for i in n:-1:1
    evals[i] = eval_one(evals[i + 1], minorvars[i], alphas[i])
  end

  tdegs = degrees(evals[n + 1])
  liftdegs = [tdegs[minorvars[i]] for i in 1:n]

  for j in 1:r
    @assert isconstant(lcs[1, j])
    fac[j] *= divexact(lcs[1, j], get_lc(fac[j], mainvar))
  end

  for i in 2:n+1
    tfac = [set_lc(fac[j], mainvar, lcs[i, j]) for j in 1:r]
    ok, fac = hliftstep(tfac, mainvar, minorvars[1:i-1], liftdegs, alphas, evals[i], true)
    if !ok
      return false, fac
    end
  end

  if !isconstant(m)
    fac = [primitive_part(i, mainvar) for i in fac]
  end

  return true, fac
end

function hliftstep(
  fac::Vector{E},
  mainvar::Int,
  minorvars::Vector{Int},
  liftdegs::Vector{Int},
  alphas::Vector,
  a::E,
  check::Bool
) where E

  r = length(fac)
  @assert r >= 2

  if r == 2
    return hliftstep_quartic2(fac, mainvar, minorvars, liftdegs, alphas, a, check)
  elseif r < 30
    return hliftstep_quartic(fac, mainvar, minorvars, liftdegs, alphas, a, check)
  else
    return hliftstep_quintic(fac, mainvar, minorvars, liftdegs, alphas, a, check)
  end
end

function hliftstep_quartic2(
  fac::Vector{E},
  mainvar::Int,
  minorvars::Vector{Int},
  liftdegs::Vector{Int},
  alphas::Vector,
  a::E,
  check::Bool
) where E

  R = parent(a)
  m = length(minorvars)

  @assert m > 0
  @assert length(liftdegs) >= m
  @assert length(alphas) >= m
  @assert length(fac) == 2

  xalpha = gen(R, minorvars[m]) - alphas[m]

  newfac = E[R(0), R(0)]

  Blen = Int[0, 0]
  B = [E[], E[]]
  for i in 1:2
    taylor_get_coeffs!(B[i], fac[i], xalpha)
    Blen[i] = length(B[i])
    # push extra stuff to avoid annoying length checks
    while length(B[i]) < liftdegs[m] + 1
      push!(B[i], R(0))
    end
  end

  ok, I = pfracinit([B[1][1], B[2][1]], mainvar, minorvars[1:m-1], alphas)
  @assert ok

  a, t = divrem(a, xalpha)
  @assert t == B[1][1] * B[2][1]

  M = zeros(R, liftdegs[m] + 1)

  for j in 1:liftdegs[m]
    a, t = divrem(a, xalpha)
    for i in 0:j
      p1 = B[1][1 + i] * B[2][1 + j - i]
      if !iszero(p1)
        t -= p1
      end
    end

    if iszero(t)
      continue
    end

    ok, deltas = pfrac(I, t, liftdegs, check)
    if !ok
      return false, newfac
    end

    tdeg = 0
    for i in 1:2
      B[i][1 + j] += deltas[i]
      if !iszero(B[i][j + 1])
        Blen[i] = max(Blen[i], j + 1)
      end
      tdeg += Blen[i] - 1
    end

    if check && tdeg > liftdegs[m]
##      println("high degree: ", tdeg, " > ", liftdegs[m])
      return false, newfac
    end

    M[1 + j] = B[1][1 + j]*B[2][1 + j]
  end

  for i in 1:2
    while length(B[i]) > Blen[i]
      @assert iszero(B[i][end])
      pop!(B[i])
    end
    newfac[i] = taylor_from_coeffs(B[i], xalpha)
  end

  return true, newfac
end

function hliftstep_quartic(
  fac::Vector{E},
  mainvar::Int,
  minorvars::Vector{Int},
  liftdegs::Vector{Int},
  alphas::Vector,
  a::E,
  check::Bool
) where E

  R = parent(a)
  r = length(fac)
  m = length(minorvars)

  @assert m > 0
  @assert length(liftdegs) >= m
  @assert length(alphas) >= m
  @assert r > 2

  xalpha = gen(R, minorvars[m]) - alphas[m]

  newfac = E[R(0) for i in 1:r]

  Blen = zeros(Int, r)
  B = [E[] for i in 1:r]
  U = [E[R(0) for j in 1:liftdegs[m] + 1] for i in 1:r]
  for i in 1:r
    taylor_get_coeffs!(B[i], fac[i], xalpha)
    Blen[i] = length(B[i])
    # push extra stuff to avoid anoying length checks
    while length(B[i]) < liftdegs[m] + 1
      push!(B[i], R(0))
    end
  end

  ok, I = pfracinit([B[i][1] for i in 1:r], mainvar, minorvars[1:m-1], alphas)
  @assert ok

  # lets pretend julia is 0-based for this awful mess

  k = r - 2
  U[1+k][1+0] = B[1+k][1+0] * B[1+k+1][1+0]
  k -= 1
  while k >= 1
    U[1+k][1+0] = B[1+k][1+0] * U[1+k+1][1+0]
    k -= 1
  end

  a, t = divrem(a, xalpha)
  @assert t == prod(B[i][1] for i in 1:r)

  for j in 1:liftdegs[m]
    k = r - 2

    U[1+k][1+j] = R(0)
    for i in 0:j
      U[1+k][1+j] += B[1+k][1+i] * B[1+k+1][1+j-i]
    end
    k -= 1
    while k >= 1
      U[1+k][1+j] = R(0)
      for i in 0:j
        U[1+k][1+j] += B[1+k][1+i] * U[1+k+1][1+j-i]
      end
      k -= 1
    end

    a, t = divrem(a, xalpha)
    for i in 0:j
      t -= B[1+0][1+i]*U[1+1][1+j-i]
    end

    if iszero(t)
      continue
    end

    ok, deltas = pfrac(I, t, liftdegs, check)
    if !ok
      return false, newfac
    end

    tdeg = 0
    for i in 1:r
      B[i][1+j] += deltas[i]
      if !iszero(B[i][1+j])
        Blen[i] = max(Blen[i], j + 1)
      end
      tdeg += Blen[i] - 1
    end

    if check && tdeg > liftdegs[m]
      return false, fac
    end

    k = r - 2
    t = B[1+k][1+0]*deltas[1+k+1] + deltas[1+k]*B[1+k+1][1+0]
    U[1+k][1+j] += t
    k -= 1
    while k >= 1
      t = B[1+k][1+0]*t + deltas[1+k]*U[1+k+1][1+0]
      U[1+k][1+j] += t
      k -= 1
    end
  end

  for i in 1:r
    while length(B[i]) > Blen[i]
      @assert iszero(B[i][end])
      pop!(B[i])
    end
    newfac[i] = taylor_from_coeffs(B[i], xalpha)
  end

  return true, newfac
end

function hliftstep_quintic(
  fac::Vector{E},
  mainvar::Int,
  minorvars::Vector{Int},
  liftdegs::Vector{Int},
  alphas::Vector,
  a::E,
  check::Bool
) where E

  R = parent(a)
  r = length(fac)
  m = length(minorvars)

  @assert m > 0
  @assert length(liftdegs) >= m
  @assert length(alphas) >= m
  @assert r > 2

  xalpha = gen(R, minorvars[m]) - alphas[m]

  betas = [eval_one(f, minorvars[m], alphas[m]) for f in fac]
  ok, I = pfracinit(betas, mainvar, minorvars[1:m-1], alphas)
  @assert ok

  newfac = [f for f in fac]
  err = a - prod(f for f in newfac)

  for j in 1:liftdegs[m]
    if iszero(err)
      return true, newfac
    end
    _, t = divrem(divexact(err, xalpha^j), xalpha)
    ok, deltas = pfrac(I, t, liftdegs, check)
    if !ok
      return false, newfac
    end
    for i in 1:r
      newfac[i] += deltas[i]*xalpha^j
    end
    err = a - prod(f for f in newfac)
  end

  return (!check || iszero(err)), newfac
end

###################### generic factoring ######################################

#=
  Factor a into irreducibles assuming a depends only on variable var
  Return (true, monic factors of a) if a is squarefree
         (false, junk)              if a is not squarefree
=#
function mfactor_irred_univar(a::E, var::Int) where E
  R = parent(a)
  K = base_ring(R)
  Kx, _ = PolynomialRing(K, "x")
  F = Hecke.factor(to_univar(a, var, Kx))
  res = E[]
  ok = true
  for f in F.fac
    ok = ok && (f[2] == 1)
    push!(res, make_monic(from_univar(f[1], var, R)))
  end
  return ok, res
end

#=
  Factor a into irreducibles assuming
    a depends only on variables xvar and yvar, and
    a is squarefree and primitive wrt xvar
  Return monic squarefree factors of a
=#
function mfactor_irred_bivar_char_zero(a::E, xvar::Int, yvar::Int) where E

  xdeg = degree(a, xvar)

  alpha = 0
  @goto got_alpha

@label next_alpha

  if alpha > 0
    alpha = -alpha
  else
    alpha = 1 - alpha
  end

@label got_alpha

  u = eval_one(a, yvar, alpha)
  if degree(u, xvar) != xdeg
    @goto next_alpha
  end

  sqrfree, ufac = mfactor_irred_univar(u, xvar)
  if !sqrfree
    @goto next_alpha
  end

  ok, cont, res = hlift_bivar_combine(a, xvar, yvar, alpha, ufac)
  if !ok
    @goto next_alpha
  end

  @assert isconstant(cont)

  return map(make_monic, res)
end

#=
  return (true, cont(a, x), array of factors) or
         (false, junk, junk) if the ufacs don't lift straight to bivariate ones

  TODO: allow this function to do a bit of zassenhaus and add a parameter to
        limit the size of the subsets.
=#
function hlift_bivar_combine(
  a::E,
  xvar::Int,
  yvar::Int,
  alpha,
  ufacs::Vector{E}  # factorization of a(x, y = alpha)
) where E

  R = parent(a)
  K = base_ring(R)
  r = length(ufacs)
  tfac = E[R(0) for i in 1:r]

  if r < 2
    tfac[1] = primitive_part(a, xvar)
    return true, divexact(a, tfac[1]), tfac
  end

  xdeg = degree(a, xvar)
  ydeg = degree(a, yvar)

  Ky, y = PolynomialRing(K, "x")

  yalpha = gen(R, yvar) - R(alpha)
  yalphapow = yalpha^(ydeg + 1)

  lc = coeff(a, Int[xvar], Int[xdeg])
  lc = evaluate(lc, [k == yvar ? y : Ky(0) for k in 1:nvars(R)])
  g, s, lcinv = gcdx((y - alpha)^(ydeg + 1), lc)
  @assert isone(g)
  lcinv = evaluate(lcinv, gen(R, yvar))

  monica = divrem(lcinv*a, yalphapow)[2]
  @assert isone(coeff(monica, Int[xvar], Int[xdeg]))

  ok, bfacs = hliftstep(ufacs, xvar, [yvar], [ydeg], [alpha], monica, false)
  if !ok
    return false, a, tfac
  end
  for i in 1:r
    @assert degree(a, xvar) >= 0
    tfac[i] = coeff(a, Int[xvar], Int[degree(a, xvar)]) * bfacs[i]
    tfac[i] = primitive_part(divrem(tfac[i], yalphapow)[2], xvar)
    a, r = divrem(a, tfac[i])
    if !iszero(r)
      return false, a, tfac
    end
  end

  return true, a, tfac
end

#=
  a and b are both factorizations. Make the bases coprime without changing
  the values of factorizations. TODO this is probably done somewhere else.
=#
function make_bases_coprime!(a::Array{Pair{E, Int}}, b::Array{Pair{E, Int}}) where E
  lena = length(a)
  lenb = length(b)
  for i in 1:lena
    for j in 1:lenb
      ai = a[i].first
      bj = b[j].first
      (g, ai, bi) = gcdcofactors(ai, bj)
      if !isconstant(g)
        a[i] = ai => a[i].second
        b[i] = bi => b[i].second
        push!(a, g => a[i].second)
        push!(b, g => b[i].second)
      end
    end
  end
  filter!(t->!isconstant(t.first), a)
  filter!(t->!isconstant(t.first), b)
end

# Return A/b^bexp
function divexact_pow(A::Fac{E}, b::E, bexp::Int) where E

  R = parent(A.unit)
  a = collect(A.fac)
  abases = E[t.first for t in a]
  aexps = Int[t.second for t in a]

  i = 1 # index strickly before which everthing is coprime to b

  while i <= length(abases) && !isconstant(b)
    abase_new, abases[i], b = gcdcofactors(abases[i], b)
    if isconstant(abase_new)
      i += 1
      continue
    end
    aexp_new = aexps[i] - bexp
    if aexp_new < 0
      error("non-exact division in divexact_pow")
    elseif aexp_new > 0
      push!(abases, abase_new)
      push!(aexps, aexp_new)
    end
    if isconstant(abases[i])
      deleteat!(abases, i)
      deleteat!(aexps, i)
    else
      i += 1
    end
  end

  if !isconstant(b)
    error("non-exact division in divexact_pow")
  end

  #return Fac{E}(R(1), Dict(abases .=> aexps))
  f = Fac{E}()
  f.unit = R(1)
  f.fac = Dict(abases .=> aexps)
  return f
end

function lcc_kaltofen_step!(
  divs::Vector{E},  # modifed
  Af::Fac{E},       # unmodified, possibly new one is returned
  Au::Vector{E},    # univariates in gen(v) from the lc's of bvar factors
  v::Int,           # the main variable for this step
  minorvars::Vector{Int},
  alphas::Vector
) where E

  R = parent(Af.unit)
  r = length(Au)
  @assert r == length(divs)
  Kx, _ = PolynomialRing(base_ring(R), string(gen(R,v)))

  Auf = [collect(Hecke.factor_squarefree(to_univar(Au[i], v, Kx)).fac) for i in 1:r]

  Afdegv = 0
  Afp = R(1)
  for i in Af.fac
    thisdeg = degree(i.first, v)
    Afdegv += thisdeg
    if thisdeg != 0
      Afp *= i.first
    end
  end

  t = Afp
  for i in 1:length(minorvars)
    t = eval_one(t, minorvars[i], alphas[i])
  end
  t = to_univar(make_monic(t), v, Kx)

  if degree(t) != Afdegv || Afdegv < 1
    return false, Af
  end

  for i in 1:r-1
    for j in i+1:r
      make_bases_coprime!(Auf[i], Auf[j])
    end
  end

  f = Set{elem_type(Kx)}()
  for i in Auf
    for j in i
      push!(f, make_monic(j.first))
    end
  end
  f = collect(f)

  if t != prod(i for i in f)
    return false, Af
  end

  liftarg = [from_univar(i, v, R) for i in f]
  ok, lifts = hlift_without_lcc(Afp, liftarg, v, minorvars, alphas)
  if !ok
    return false, Af
  end

  for i in 1:r
    for j in Auf[i]
      for k in 1:length(f)
        if f[k] == j.first
          Af = divexact_pow(Af, lifts[k], j.second)
          divs[i] *= lifts[k]^j.second
        end
      end
    end
  end
  return true, Af
end

#=
  Try to determine divisors of the leading coefficients of the factors of A.
  This is accomplished by looking at the bivariate factoration of A when
  all but one of the minor variables are evaluated away. The resulting
  univariate leading coefficients are lifted against the supplied
  leading_coefficient(A) factorization. return is Tuple{::Boole, ::Vector{E}}
  If the bool is true, then the method can be considered to have fully found
  the leading coefficients, otherwise not. In any case, the second return
  is a tentative guess at the leading coefficients.
=#
function lcc_kaltofen(
  lcAf::Fac{E},       # square free factorization of lc(A), not modified
  A::E,
  mainvar::Int,
  maindeg::Int,
  minorvars::Vector{Int},
  alphas::Vector,
  ufacs::Vector{E}    # univariate factors of (A evaluated at minor vars)
) where E

  R = parent(A)
  r = length(ufacs)

  @assert r > 1

  divs = E[R(1) for i in 1:r]
  ulcs = E[R() for i in 1:r]

  for vi in 1:length(minorvars)
    if isempty(lcAf.fac)
      break
    end

    other_minorvars = deleteat!(copy(minorvars), vi)
    other_alphas = deleteat!(copy(alphas), vi)

    beval = A
    for i in 1:length(other_minorvars)
      beval = eval_one(beval, other_minorvars[i], other_alphas[i])
    end
    @assert degree(beval, mainvar) == maindeg

    ok, cont, bfac = hlift_bivar_combine(beval, mainvar, minorvars[vi],
                                                             alphas[vi], ufacs)
    if !ok
      return false, divs
    end
    if !isconstant(cont)
      continue
    end

    divides_ok = true
    for k in 1:r
      ueval = divs[k]
      for i in 1:length(other_minorvars)
        ueval = eval_one(ueval, other_minorvars[i], other_alphas[i])
      end
      ok, ulcs[k] = divides(get_lc(bfac[k], mainvar), ueval)
      divides_ok = divides_ok && ok
    end
    if !divides_ok
      continue
    end
    ok, lcAf = lcc_kaltofen_step!(divs, lcAf, ulcs, minorvars[vi],
                                                 other_minorvars, other_alphas)
  end

  return isempty(lcAf.fac), divs

end

# factor a truely multivariate A in at least three variables
function mfactor_irred_mvar_char_zero(
  A::E,             # squarefree, primitve wrt mainvar, monic
  mainvar::Int,
  minorvars::Vector{Int}
) where E

  n = length(minorvars)
  R = parent(A)
  K = base_ring(R)
  @assert length(A) > 0

  evals = zeros(R, n + 1)
  alphas = zeros(K, n)
  alpha_modulus = 0
  lcc_fails_remaining = 3

  maindeg = degree(A, mainvar)

@label next_alpha

  alpha_modulus += 1
  for i in 1:n
    alphas[i] = K(rand(Int) % alpha_modulus)
  end

  evals[n + 1] = A
  for i in n:-1:1
    evals[i] = eval_one(evals[i + 1], minorvars[i], alphas[i])
    if degree(evals[i], mainvar) != maindeg
      @goto next_alpha
    end
  end

  # make sure univar is squarefree. TODO also zassenhaus pruning here
  ok, fac = mfactor_irred_univar(evals[1], mainvar)
  if !ok
    @goto next_alpha
  end

  r = length(fac)
  if r < 2
    return [A]
  end

  lcAf = mfactor_squarefree_char_zero(get_lc(A, mainvar))

  ok, divs = lcc_kaltofen(lcAf, A, mainvar, maindeg, minorvars, alphas, fac)
  if !ok
    lcc_fails_remaining -= 1
    if lcc_fails_remaining >= 0
      @goto next_alpha
    end
  end

  ok, fac = hlift_with_lcc(A, fac, divs, mainvar, minorvars, alphas)
  if !ok
    @goto next_alpha
  end

  return map(make_monic, fac)
end

########## factorization in three steps #######################################

# take out the content and lowest power wrt variable v
function mfactor_primitive(f::E, v::Int) where E
  R = parent(f)
  d = degree(f, v)
  g = R(0)
  smallest = d
  for i in 0:d
    ci = coeff(f, Int[v], Int[i])
    if !iszero(ci)
      g = gcd(g, ci)::E
      smallest = min(smallest, i)
    end
  end
  @assert !iszero(g)
  return (g, divexact(f, g*gen(R, v)^smallest)::E, smallest)
end

# assume a is primitive wrt each variable appearing in it
# return squarefree factors
function mfactor_sqrfree_char_zero(a::E) where E
  R = parent(a)
  @assert length(a) >= 1
  res = Fac{E}()
  res.unit = R(1)
  for v in 1:nvars(R)
    Sp = derivative(a, v)
    if !iszero(Sp)
      (Sm, Ss, Y) = gcdcofactors(a, Sp)
      k = 1
      while begin Z = Y - derivative(Ss, v); !iszero(Z) end
        (S, Ss, Y) = gcdcofactors(Ss, Z)
        mulpow!(res, S, k)
        k += 1
      end
      mulpow!(res, Ss, k)
      return res
    end
  end
  @assert isconstant(a)
  mulpow!(res, a, 1)
  return res
end

# assume a is primitive wrt each variable appearing in it and squarefree
# return irreducible factors
function mfactor_irred_char_zero(a::E) where E
  R = parent(a)
  K = base_ring(R)
  @assert length(a) > 0

  res = Fac{E}()
  res.unit = R(1)

  lc = coeff(a, 1)
  if !isone(lc)
    res.unit = lc
    a *= inv(lc)
  end

  degs = degrees(a)
  vars = Int[]    # variables that actually appear
  for v in 1:nvars(R)
    if degs[v] > 0
      push!(vars, v)
    end
  end

  if isempty(vars)
    @assert isone(a)
    return res
  end

  sort!(vars, by = (v -> degs[v]), alg=InsertionSort)

  if degs[1] == 1
    # linear is irreducible by assumption
    mulpow!(res, a, 1)
    return res
  end

  local f::Vector{E}

  if length(vars) == 1
    sqrfree, f = mfactor_irred_univar(a, vars[1])
    @assert sqrfree
  elseif length(vars) == 2
    f = mfactor_irred_bivar_char_zero(a, vars[1], vars[2])
  else
    f = mfactor_irred_mvar_char_zero(a, vars[1], vars[2:end])
  end

  for i in f
    mulpow!(res, i, 1)
  end
  return res
end

function mfactor_squarefree_char_zero(a::E) where E
  R = parent(a)
  K = base_ring(R)

  res = Fac{E}()
  tres = Fac{E}()

  if iszero(a)
    res.unit = R(0)   # not a unit :-(
    return res
  end

  # start with a monic version of a
  lc = coeff(a, 1)
  res.unit = R(lc)
  res.fac = Dict{E, Int}(1//lc * a => 1)

  # pure variable powers in the final factorization
  var_powers = zeros(Int, nvars(R))

  # ensure factors are primitive wrt any variable
  for v in 1:nvars(R)
    tres.unit = res.unit
    empty!(tres.fac)
    for p in res.fac
      (content, primitive, power) = mfactor_primitive(p[1], v)
      var_powers[v] += power
      mulpow!(tres, content, 1)
      mulpow!(tres, primitive, 1)
    end
    res, tres = tres, res
  end

  # ensure factors are squarefree
  tres.unit = res.unit
  empty!(tres.fac)
  for i in res.fac
    mulpow!(tres, mfactor_sqrfree_char_zero(i[1]), i[2])
  end
  res, tres = tres, res

  # put pure variable powers back in
  for v in 1:nvars(R)
    mulpow!(res, gen(R, v), var_powers[v])
  end

  return res
end

# factor a multivariate over an exact field of characteristic 0
function mfactor_char_zero(a::E) where E
  tres = mfactor_squarefree_char_zero(a)
  # ensure factors are irreducible
  res = Fac{E}()
  res.unit = tres.unit
  empty!(res.fac)
  for i in tres.fac
    mulpow!(res, mfactor_irred_char_zero(i[1]), i[2])
  end
  return res
end

function factor(a::MPolyElem)
  R = parent(a)
  K = base_ring(R)
  if elem_type(K) <: AbstractAlgebra.FieldElem && iszero(characteristic(K))
    return mfactor_char_zero(a)
  else
    error("factor not implemented for the ring $R")
  end
end

function factor_squarefree(a::MPolyElem)
  R = parent(a)
  K = base_ring(R)
  if elem_type(K) <: AbstractAlgebra.FieldElem && iszero(characteristic(K))
    return mfactor_squarefree_char_zero(a)
  else
    error("factor_squarefree not implemented for the ring $R")
  end
end


text:
using DashBootstrapComponents, DashHtmlComponents

pagination = html_div([
    dbc_pagination(max_value = 5, first_last = true),
    dbc_pagination(max_value = 5, previous_next = true),
    dbc_pagination(max_value = 5, first_last = true, previous_next = true),
])

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
# module Ecto

# in schema.jl
export Schema, put_meta, rem, |>

# in changeset.jl
export Changeset

# in query.jl
export Query, in, from, and, set

# in repo.jl
export Repo