text:
#=
This file actually implements the user interface of the Psychro library.
=#

"""
    volume(DryAir, T, P[, out_unit=u"m^3/kg"])
    volume(Vapor, T[, out_unit=u"m^3/kg"])
    volume(MoistAir, T, HumidityType, hum, P[, out_unit="m^3/kg"])

Calculates the specific volume of a gas. In the case of moist air, parameters specifying the humidity shoud be added and the the property is based on the mass of dry air not mass of the gas (as is the case for `DryAir` or `Vapor`).

If `Unitful` units are used, all parameters of the function (except dimensionless) should have units associated. When no units are provided, all parameters should use SI units (and so does the return type). For further information, checkout the doc for `Psychro` module.

Temperature should be in the range `173.15 K < T < 473.15` and pressure should be below 5 MPa.

# Examples
```julia-repl
julia> volume(DryAir, 293.15, 101325.0)
0.8302353995092402

julia> volume(DryAir, 20.0u"°C", 1.0u"atm")
0.8302353995092402 kg^-1 m^3

julia> volume(Vapor, 293.15)
57.77492045936827

julia> volume(Vapor, 20.0u"°C")
57.77492045936827 kg^-1 m^3

julia> volume(MoistAir, 293.15, RelHum, 0.7, 101325.0)
0.843889817602806

julia> volume(MoistAir, 20.0u"°C", WetBulb, 18.0u"°C", 1.0u"atm")
0.8464079202783964 kg^-1 m^3

julia> volume(MoistAir, 20.0u"°C", WetBulb, 18.0u"°C", 1.0u"atm", u"inch^3/lb")
23428.490579267214 in^3 lb^-1
```
"""
function volume(::Type{DryAir}, Tk, P)
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    @assert 0 <= P < 5e6 "Out of range error: Pressure should be below 5×10⁶ Pa"
    volumeair(Tk, P)
end

"""
    volumem(DryAir, T, P[, out_unit=u"m^3/kg"])
    volumem(Vapor, T[, out_unit=u"m^3/kg"])
    volumem(MoistAir, T, HumidityType, hum, P[, out_unit="m^3/kg"])

Calculates the molar volume of a gas. In the case of moist air, parameters specifying the humidity shoud be added.

If `Unitful` units are used, all parameters of the function (except dimensionless) should have units associated. When no units are provided, all parameters should use SI units (and so does the return type). For further information, checkout the doc for `Psychro` module.

Temperature should be in the range `173.15 K < T < 473.15` and pressure should be below 5 MPa.

# Examples
```julia-repl
julia> volumem(DryAir, 293.15, 101325.0)
0.024046522993685877

julia> volumem(DryAir, 20.0u"°C", 1.0u"atm")
0.024046522993685877 m^3 mol^-1

julia> volumem(Vapor, 293.15)
1.040831369053248

julia> volumem(Vapor, 20.0u"°C")
1.040831369053248 m^3 mol^-1

julia> volumem(MoistAir, 293.15, RelHum, 0.7, 101325.0)
0.02404547233948706

julia> volumem(MoistAir, 20.0u"°C", WetBulb, 18.0u"°C", 1.0u"atm")
0.024045209859305458 m^3 mol^-1

julia> volumem(MoistAir, 20.0u"°C", WetBulb, 18.0u"°C", 1.0u"atm", u"inch^3/mol")
1467.32873315839 in^3 mol^-1
```
"""
function volumem(::Type{DryAir}, Tk, P)
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    @assert 0 <= P < 5e6 "Out of range error: Pressure should be below 5×10⁶ Pa"
    molarvolumeair(Tk, P)
end

"""
    density(DryAir, T, P[, out_unit=u"m^3/kg"])
    density(Vapor, T[, out_unit=u"m^3/kg"])
    density(MoistAir, T, HumidityType, hum, P[, out_unit="m^3/kg"])

Calculates the density volume of a gas. In the case of moist air, parameters specifying the humidity shoud be added.

If `Unitful` units are used, all parameters of the function (except dimensionless) should have units associated. When no units are provided, all parameters should use SI units (and so does the return type). For further information, checkout the doc for `Psychro` module.

Temperature should be in the range `173.15 K < T < 473.15` and pressure should be below 5 MPa.

# Examples
```julia-repl
julia> density(DryAir, 293.15, 101325.0)
1.2044776705391136

julia> density(DryAir, 20.0u"°C", 1.0u"atm")
1.2044776705391136 kg m^-3

julia> density(Vapor, 293.15)
0.017308548277505224

julia> density(Vapor, 20.0u"°C")
0.017308548277505224 kg m^-3

julia> density(MoistAir, 293.15, RelHum, 0.7, 101325.0)
1.1971436091840653

julia> density(MoistAir, 20.0u"°C", WetBulb, 18.0u"°C", 1.0u"atm")
1.195819185774941 kg m^-3

julia> density(MoistAir, 20.0u"°C", WetBulb, 18.0u"°C", 1.0u"atm", u"lb/inch^3")
4.3201708903793606e-5 in^-3 lb

```
"""
function density(::Type{DryAir}, Tk, P)
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    @assert 0 <= P < 5e6 "Out of range error: Pressure should be below 5×10⁶ Pa"
    1.0/volumeair(Tk, P)
end

"""
    enthalpy(DryAir, T, P[, out_unit=u"m^3/kg"])
    enthalpy(Vapor, T[, out_unit=u"m^3/kg"])
    enthalpy(MoistAir, T, HumidityType, hum, P[, out_unit="m^3/kg"])

Calculates the specific enthalpy of a gas. In the case of moist air, parameters specifying the humidity shoud be added and the the property is based on the mass of dry air not mass of the gas (as is the case for `DryAir` or `Vapor`).

If `Unitful` units are used, all parameters of the function (except dimensionless) should have units associated. When no units are provided, all parameters should use SI units (and so does the return type). For further information, checkout the doc for `Psychro` module.

Temperature should be in the range `173.15 K < T < 473.15` and pressure should be below 5 MPa.

# Examples
```julia-repl
julia> enthalpy(DryAir, 293.15, 101325.0)
20121.2722813185

julia> enthalpy(DryAir, 20.0u"°C", 1.0u"atm")
20121.2722813185 kg^-1 J

julia> enthalpy(Vapor, 293.15)
2.5374003352412493e6

julia> enthalpy(Vapor, 20.0u"°C")
2.5374003352412493e6 kg^-1 J

julia> enthalpy(MoistAir, 293.15, RelHum, 0.7, 101325.0)
46142.773129687484

julia> enthalpy(MoistAir, 20.0u"°C", WetBulb, 18.0u"°C", 1.0u"atm")
50944.84575377501 kg^-1 J

julia> enthalpy(MoistAir, 20.0u"°C", WetBulb, 18.0u"°C", 1.0u"atm", u"J/lb")
23108.193324739244 J lb^-1

```
"""
function enthalpy(::Type{DryAir}, Tk, P)
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    @assert 0 <= P < 5e6 "Out of range error: Pressure should be below 5×10⁶ Pa"
    enthalpyair(Tk, P)
end

"""
    enthalpym(DryAir, T, P[, out_unit=u"m^3/kg"])
    enthalpym(Vapor, T[, out_unit=u"m^3/kg"])
    enthalpym(MoistAir, T, HumidityType, hum, P[, out_unit="m^3/kg"])

Calculates the molar enthalpy of a gas. In the case of moist air, parameters specifying the humidity shoud be added.

If `Unitful` units are used, all parameters of the function (except dimensionless) should have units associated. When no units are provided, all parameters should use SI units (and so does the return type). For further information, checkout the doc for `Psychro` module.

Temperature should be in the range `173.15 K < T < 473.15` and pressure should be below 5 MPa.

# Examples
```julia-repl
julia> enthalpym(DryAir, 293.15, 101325.0)
582.7824697199684

julia> enthalpym(DryAir, 20.0u"°C", 1.0u"atm")
582.7824697199684 J mol^-1

julia> enthalpy(Vapor, 293.15)
2.5374003352412493e6

julia> enthalpym(Vapor, 20.0u"°C")
45711.977511464975 J mol^-1

julia> enthalpym(MoistAir, 293.15, RelHum, 0.7, 101325.0)
1314.7744549269437

julia> enthalpym(MoistAir, 20.0u"°C", WetBulb, 18.0u"°C", 1.0u"atm")
1447.2684837312868 J mol^-1

julia> enthalpym(MoistAir, 20.0u"°C", WetBulb, 18.0u"°C", 1.0u"atm", u"kJ/kmol")
1447.2684837312868 kJ kmol^-1

```
"""
function enthalpym(::Type{DryAir}, Tk, P)
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    @assert 0 <= P < 5e6 "Out of range error: Pressure should be below 5×10⁶ Pa"
    molarenthalpyair(Tk, P)
end

"""
    entropy(DryAir, T, P[, out_unit=u"m^3/kg"])
    entropy(Vapor, T[, out_unit=u"m^3/kg"])
    entropy(MoistAir, T, HumidityType, hum, P[, out_unit="m^3/kg"])

Calculates the specific entropy of a gas. In the case of moist air, parameters specifying the humidity shoud be added and the the property is based on the mass of dry air not mass of the gas (as is the case for `DryAir` or `Vapor`).

If `Unitful` units are used, all parameters of the function (except dimensionless) should have units associated. When no units are provided, all parameters should use SI units (and so does the return type). For further information, checkout the doc for `Psychro` module.

Temperature should be in the range `173.15 K < T < 473.15` and pressure should be below 5 MPa.

# Examples
```julia-repl
julia> entropy(DryAir, 293.15, 101325.0)
71.09262577195514

julia> entropy(DryAir, 20.0u"°C", 1.0u"atm")
71.09262577195514 kg^-1 J K^-1

julia> entropy(Vapor, 293.15)
8665.873003613997

julia> entropy(Vapor, 20.0u"°C")
8665.873003613997 kg^-1 J K^-1

julia> entropy(MoistAir, 293.15, RelHum, 0.7, 101325.0)
166.33538729355692

julia> entropy(MoistAir, 20.0u"°C", WetBulb, 18.0u"°C", 1.0u"atm")
182.97140931650105 kg^-1 J K^-1

julia> entropy(MoistAir, 20.0u"°C", WetBulb, 18.0u"°C", 1.0u"atm", u"kJ/lb/°F")
0.04610801955228433 °F^-1 kJ lb^-1

```
"""
function entropy(::Type{DryAir}, Tk, P)
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    @assert 0 <= P < 5e6 "Out of range error: Pressure should be below 5×10⁶ Pa"
    entropyair(Tk, P)
end

"""
    entropym(DryAir, T, P[, out_unit=u"m^3/kg"])
    entropym(Vapor, T[, out_unit=u"m^3/kg"])
    entropym(MoistAir, T, HumidityType, hum, P[, out_unit="m^3/kg"])

Calculates the molar entropy of a gas. In the case of moist air, parameters specifying the humidity shoud be added.

If `Unitful` units are used, all parameters of the function (except dimensionless) should have units associated. When no units are provided, all parameters should use SI units (and so does the return type). For further information, checkout the doc for `Psychro` module.

Temperature should be in the range `173.15 K < T < 473.15` and pressure should be below 5 MPa.

# Examples
```julia-repl
julia> entropym(DryAir, 293.15, 101325.0)
2.0590912665460226

julia> entropym(DryAir, 20.0u"°C", 1.0u"atm")
2.0590912665460226 J K^-1 mol^-1

julia> entropy(Vapor, 293.15)
8665.873003613997

julia> entropym(Vapor, 20.0u"°C")
156.11812860454717 J K^-1 mol^-1

julia> entropym(MoistAir, 293.15, RelHum, 0.7, 101325.0)
4.739496639035868

julia> entropym(MoistAir, 20.0u"°C", WetBulb, 18.0u"°C", 1.0u"atm")
5.197949865380578 J K^-1 mol^-1

julia> entropym(MoistAir, 20.0u"°C", WetBulb, 18.0u"°C", 1.0u"atm", u"kJ/kmol/°F")
2.8877499252114327 °F^-1 kJ kmol^-1

```
"""
function entropym(::Type{DryAir}, Tk, P)
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    @assert 0 <= P < 5e6 "Out of range error: Pressure should be below 5×10⁶ Pa"
    molarentropyair(Tk, P)
end

"""
    compressfactor(DryAir, T, P[, out_unit=u"m^3/kg"])
    compressfactor(Vapor, T[, out_unit=u"m^3/kg"])
    compressfactor(MoistAir, T, HumidityType, hum, P[, out_unit="m^3/kg"])

Calculates the compressibility factor Z of a gas. In the case of moist air, parameters specifying the humidity shoud be added.

If `Unitful` units are used, all parameters of the function (except dimensionless) should have units associated. When no units are provided, all parameters should use SI units (and so does the return type). For further information, checkout the doc for `Psychro` module.

Temperature should be in the range `173.15 K < T < 473.15` and pressure should be below 5 MPa.

# Examples
```julia-repl
julia> compressfactor(DryAir, 293.15, 101325.0)
1.2044776705391136

julia> compressfactor(DryAir, 20.0u"°C", 1.0u"atm")
1.2044776705391136 kg m^-3

julia> compressfactor(Vapor, 293.15)
0.017308548277505224

julia> compressfactor(Vapor, 20.0u"°C")
0.017308548277505224 kg m^-3

julia> compressfactor(MoistAir, 293.15, RelHum, 0.7, 101325.0)
1.1971436091840653

julia> compressfactor(MoistAir, 20.0u"°C", WetBulb, 18.0u"°C", 1.0u"atm")
1.195819185774941 kg m^-3

```
"""
function compressfactor(::Type{DryAir}, Tk, P)
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    @assert 0 <= P < 5e6 "Out of range error: Pressure should be below 5×10⁶ Pa"
    P*molarvolumeair(Tk, P) / (R*Tk)
end

function volume(::Type{Vapor}, Tk)
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    volumevapor(Tk)
end

function volumem(::Type{Vapor}, Tk)
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    molarvolvapor(Tk)
end

function density(::Type{Vapor}, Tk)
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    1/volumevapor(Tk)
end

function enthalpy(::Type{Vapor}, Tk)
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    enthalpyvapor(Tk)
end

function enthalpym(::Type{Vapor}, Tk)
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    enthalpyvapor(Tk)*Mv
end

function entropym(::Type{Vapor}, Tk)
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    molarentropyvapor(Tk)
end

function entropy(::Type{Vapor}, Tk)
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    entropyvapor(Tk)
end

function compressfactor(::Type{Vapor}, Tk)
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    Zvapor(Tk)
end

"""
    molarfrac(Tk, HumidityType, x, P)

Calculates the molarfraction of water vapor given a humidity parameter.
The humidity parameter can be one of the following:

 * `MolarFrac` - Identity function basically...
 * `RelHum` - Relative humidity in decimal notation (not percentage!)
 * `DewPoint`- Dew point temperature in K
 * `HumRat` - Humidity ration, kg of vapor / kg of dry air
 * `SpecHum` - Specific Humidity of water vapor = mv / (mv + ma)
 * `WetBulb` - Wet bulb temperature in K- estimated by the adiabatic saturation temperature.

As usual all calculations are performed using plain SI units.
The parameters used by this function are:

 * `Tk` - Temperature in K
 * Humidity type - a type listed above that characterizes the humidity
 * The value of the humidity parameter.
 * P - Pressure in Pa
"""
function molarfrac(Tk, ::Type{MolarFrac}, xv, P)
    xv
end

function molarfrac(Tk, ::Type{HumRat}, w, P)
    w  / (Mv/Ma + w)
end

function molarfrac(Tk, ::Type{SpecHum}, q, P)
    q * Ma / (Mv + r*(Ma - Mv))
end

function molarfrac(Tk, ::Type{RelHum}, rel, P)
    rel * efactor(Tk, P) * Pws(Tk) / P
end

function molarfrac(Tk, ::Type{DewPoint}, D, P)
    efactor(D,P) * Pws(D) / P
end

function molarfrac(Tk, ::Type{WetBulb}, B, P)
    w = calc_W_from_B(Tk, B, P)
    w / (Mv/Ma + w)
end

function volume(::Type{MoistAir}, Tk, ::Type{T}, y, P) where {T<:PsychroProperty}
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15! K"
    @assert 0 <= P < 5e6 "Out of range error: Pressure should be below 5×10⁶ Pa"
    xv = molarfrac(Tk, T, y, P)
    volumemoist(Tk, P, xv)
end

function volumem(::Type{MoistAir}, Tk, ::Type{T}, y, P) where {T<:PsychroProperty}
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15! K"
    @assert 0 <= P < 5e6 "Out of range error: Pressure should be below 5×10⁶ Pa"

    xv = molarfrac(Tk, T, y, P)
    molarvolumemoist(Tk, P, xv)
end

function density(::Type{MoistAir}, Tk, ::Type{T}, y, P) where {T<:PsychroProperty}
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    @assert 0 <= P < 5e6 "Out of range error: Pressure should be below 5×10⁶ Pa"

    xv = molarfrac(Tk, T, y, P)
    M = xv*Mv + (1-xv)*Ma
    vm = molarvolumemoist(Tk, P, xv)
    return M/vm
end

function enthalpy(::Type{MoistAir}, Tk, ::Type{T}, y, P) where {T<:PsychroProperty}
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    @assert 0 <= P < 5e6 "Out of range error: Pressure should be below 5×10⁶ Pa"

    xv = molarfrac(Tk, T, y, P)
    enthalpymoist(Tk, P, xv)
end

function enthalpym(::Type{MoistAir}, Tk, ::Type{T}, y, P) where {T<:PsychroProperty}
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    @assert 0 <= P < 5e6 "Out of range error: Pressure should be below 5×10⁶ Pa"

    xv = molarfrac(Tk, T, y, P)
    molarenthalpymoist(Tk, P, xv)
end

function entropy(::Type{MoistAir}, Tk, ::Type{T}, y, P) where {T<:PsychroProperty}
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    @assert 0 <= P < 5e6 "Out of range error: Pressure should be below 5×10⁶ Pa"

    xv = molarfrac(Tk, T, y, P)
    entropymoist(Tk, P, xv)
end

function entropym(::Type{MoistAir}, Tk, ::Type{T}, y, P) where {T<:PsychroProperty}
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    @assert 0 <= P < 5e6 "Out of range error: Pressure should be below 5×10⁶ Pa"

    xv = molarfrac(Tk, T, y, P)
    molarentropymoist(Tk, P, xv)
end

function compressfactor(::Type{MoistAir}, Tk, ::Type{T}, y, P) where {T<:PsychroProperty}
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    @assert 0 <= P < 5e6 "Out of range error: Pressure should be below 5×10⁶ Pa"

    xv = molarfrac(Tk, T, y, P)
    Zmoist(Tk, P, xv)
end

function calcdewpoint(Tk, P, xv, EPS=1e-9, MAXITER=100)
    # Use Ideal Gas to
    D = Tws(xv*P)
    Dnew = D
    err = 0.0
    niter = 1
    for i = 1:MAXITER
        niter = i
        f = efactor(D, P)
        Dnew = Tws(xv*P/f)
        err = abs(D-Dnew)
        if err < EPS
            return Dnew
        end
        D = Dnew
    end
    throw(ConvergenceError("Dew point calculation failed to converge!", D, niter, err))
    return Dnew
end

"""
# Thermodynamic properties of moist air, dry air and saturated water vapor.

The methods listed below calculate thermodynamic properties of moist air:

    volume(MoistAir, T, HumidityType, humidity, P[, outunit])
    volume(MoistAir, T, HumidityType, humidity, P[, outunit])
    density(MoistAir, T, HumidityType, humidity, P[, outunit])
    enthalpy(MoistAir, T, HumidityType, humidity, P[, outunit])
    enthalpym(MoistAir, T, HumidityType, humidity, P[, outunit])
    entropy(MoistAir, T, HumidityType, humidity, P[, outunit])
    entropym(MoistAir, T, HumidityType, humidity, P[, outunit])
    compressfactor(MoistAir, T, HumidityType, humidity, P[, outunit])
    dewpoint(MoistAir, T, HumidityType, humidity, P[, outunit])
    wetbulb(MoistAir, T, HumidityType, humidity, P[, outunit])
    humrat(MoistAir, T, HumidityType, humidity, P)
    relhum(MoistAir, T, HumidityType, humidity, P)
    humrat(MoistAir, T, HumidityType, humidity, P)
    spechum(MoistAir, T, HumidityType, humidity, P)
    molarfrac(MoistAir, T, HumidityType, humidity, P)

The methods listed above calculate the following thermodynamic properties of moist air:

 * `volume` Specific volume
 * `volumem` Molar volume
 * `density` Density
 * `enthalpy` Specific enthalpy
 * `enthalpym` Molar enthalpy
 * `entropy` Specific entropy
 * `entropym` Molar entropy
 * `compressfactor` Compressibility factor Z
 * `dewpoint` Dew point temperature
 * `wetbulb` Adiabatic saturation temperature
 * `humrat` Humidity ratio
 * `relhum` Relative humidity
 * `spechum` Specific humidity
 * `molarfrac` Molar fraction of water vapor

The humidity is specified using two parameters:

 * How the humidity is specified
 * The actual value of humidity

The following types are used to characterize the humidity.

 * `WetBulb` for wet bulb temperature, actually adiabatic saturation temperature
 * `DewPoint` Dew point temperature
 * `RelHum` Relative humidity
 * `HumRat` Humidity ratio (kg of vapor / kg of dry air)
 * `SpecHum` Specific humidity (kg of vapor / kg of moist air)
 * `MolarFrac` Molar fraction of water vapor.

# Examples
```julia-repl
julia> volume(MoistAir, 293.15, WetBulb, 291.15, 101325.0)
0.8464079202783964

julia> volume(MoistAir, 293.15, DewPoint, 291.15, 101325.0)
0.8475219875187474

julia> volume(MoistAir, 293.15, RelHum, 0.7, 101325.0)
0.843889817602806

julia> volume(MoistAir, 20.0u"°C", DewPoint, 60.0u"°F", 1.0u"atm")
0.8449934929585231 kg^-1 m^3

julia> volumem(MoistAir, 293.15, RelHum, 0.5, 93000.0)
0.026199080086890276

julia> volumem(MoistAir, 20.0u"°C", WetBulb, 17.0u"°C", 93u"kPa", u"inch^3/kmol")
1.598733210336603e6 in^3 kmol^-1

julia> density(MoistAir, 20.0u"°C", WetBulb, 17.0u"°C", 93u"kPa")
1.0976075893895811 kg m^-3

julia> density(MoistAir, 20.0u"°C", WetBulb, 17.0u"°C", 93u"kPa", u"lb/inch^3")
3.965358988338535e-5 in^-3 lb

julia> volumem(MoistAir, 20.0u"°C", WetBulb, 17.0u"°C", 93u"kPa", u"inch^3/kmol")
1.598733210336603e6 in^3 kmol^-1

julia> enthalpy(MoistAir, 20.0u"°C", WetBulb, 17.0u"°C", 93u"kPa")
50667.43014746832 kg^-1 J

julia> enthalpym(MoistAir, 20.0u"°C", WetBulb, 17.0u"°C", 93u"kPa")
1439.6551689935861 J mol^-1

julia> compressfactor(MoistAir, -90.0u"°C", RelHum, 0.01, 4.5u"MPa")
0.8552758629097985

julia> wetbulb(MoistAir, 20.0u"°C", WetBulb, 17.0u"°C", 93u"kPa", u"°C")
17.0 °C

julia> dewpoint(MoistAir, 20.0u"°C", WetBulb, 17.0u"°C", 93u"kPa", u"°C")
15.475836053510477 °C

julia> humrat(MoistAir, 20.0u"°C", WetBulb, 17.0u"°C", 93u"kPa")
0.012032930694441925

julia> relhum(MoistAir, 20.0u"°C", WetBulb, 17.0u"°C", 93u"kPa")
0.7517801524436909

julia> spechum(MoistAir, 20.0u"°C", WetBulb, 17.0u"°C", 93u"kPa")
0.011889860823189923

```

"""
function dewpoint(::Type{MoistAir}, Tk, ::Type{T}, y, P) where {T<:PsychroProperty}
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    @assert 0 <= P < 5e6 "Out of range error: Pressure should be below 5×10⁶ Pa"

    if T==DewPoint
        return y
    end

    xv = molarfrac(Tk, T, y, P)
    D = calcdewpoint(Tk, P, xv)
    return D
end

function relhum(::Type{MoistAir}, Tk, ::Type{T}, y, P) where {T<:PsychroProperty}
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    @assert 0 <= P < 5e6 "Out of range error: Pressure should be below 5×10⁶ Pa"

    if T==RelHum
        return y
    end

    xv = molarfrac(Tk, T, y, P)
    xv * P / (efactor(Tk, P) * Pws(Tk))
end
@doc (@doc dewpoint) relhum

function wetbulb(::Type{MoistAir}, Tk, ::Type{T}, y, P) where {T<:PsychroProperty}
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    @assert 0 <= P < 5e6 "Out of range error: Pressure should be below 5×10⁶ Pa"

    if T==WetBulb
        return y
    end
    xv = molarfrac(Tk, T, y, P)

    B = calcwetbulb(Tk, P, xv)
    return B

end
@doc (@doc dewpoint) wetbulb

humrat(xv) = xv / (1-xv) * (Mv/Ma)
function humrat(::Type{MoistAir}, Tk, ::Type{T}, y, P) where {T<:PsychroProperty}
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    @assert 0 <= P < 5e6 "Out of range error: Pressure should be below 5×10⁶ Pa"

    if T==HumRat
        return y
    end

    xv = molarfrac(Tk, T, y, P)
    return humrat(xv)
end
@doc (@doc dewpoint) humrat

function molarfrac(::Type{MoistAir}, Tk, ::Type{T}, y, P) where {T<:PsychroProperty}
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    @assert 0 <= P < 5e6 "Out of range error: Pressure should be below 5×10⁶ Pa"

    if T==MolarFrac
        return y
    end

    xv = molarfrac(Tk, T, y, P)
    return xv
end
@doc (@doc dewpoint) molarfrac

function spechum(::Type{MoistAir}, Tk, ::Type{T}, y, P) where {T<:PsychroProperty}
    @assert 173.1 < Tk < 473.2 "Out of range error: Temperature should be between 173.15 K and 473.15 K!"
    @assert 0 <= P < 5e6 "Out of range error: Pressure should be below 5×10⁶ Pa"

    if T==SpecHum
        return y
    end
    xv = molarfrac(Tk, T, y, P)
    return xv*Mv / ( (1-xv)*Ma + xv*Mv )

end
@doc (@doc dewpoint) spechum


text:
import .Zygote

struct ZygoteGradientLogDensity{L} <: ADGradientWrapper
    ℓ::L
end

"""
    ADgradient(:Zygote, ℓ)
    ADgradient(Val(:Zygote), ℓ)

Gradient using algorithmic/automatic differentiation via Zygote.
"""
ADgradient(::Val{:Zygote}, ℓ) = ZygoteGradientLogDensity(ℓ)

Base.show(io::IO, ∇ℓ::ZygoteGradientLogDensity) = print(io, "Zygote AD wrapper for ", ∇ℓ.ℓ)

function logdensity_and_gradient(∇ℓ::ZygoteGradientLogDensity, x::AbstractVector)
    @unpack ℓ = ∇ℓ
    y, back = Zygote.pullback(_logdensity_closure(ℓ), x)
    y, first(back(Zygote.sensitivity(y)))
end

ext: jl


text:
module METIS

using   DSP, # toolbox for digital signal processing (mostly filtering)
        MAT  #to read matlab .mat files'

# Include types
include("type/kronos.jl")

# Include internal subpackages
include("read/bvr.jl")

end


text:
module Kattis

include("origami.jl")
include("zeroone.jl")
# include("linearprog.jl")
include("maxflow.jl")

end


text:
using FermiCG
using Printf
using Test
using LinearAlgebra
using Profile
using HDF5

using PyCall
pydir = joinpath(dirname(pathof(FermiCG)), "python")
pushfirst!(PyVector(pyimport("sys")."path"), pydir)
ENV["PYTHON"] = Sys.which("python")

#@testset "Clusters" begin

out_radii = []
out_fci = []
out_ref = []
out_nb2 = []
out_cepa = []
out_rhf = []
out_mp2 = []
n_steps = 40
start = 1
stop = 5
stepsize = (stop-start)/n_steps

for step in 1:n_steps
    rad = start + (step-1) * stepsize
    atoms = []
    clusters = []
    na = 0
    nb = 0
    init_fspace = []

    function generate_H_ring(n,radius)
        theta = 2*pi/n

        atoms = []
        for i in 0:n-1
            push!(atoms,Atom(i+1,"H",[radius*cos(theta*i), radius*sin(theta*i), 0]))
        end
        return atoms
    end

    if false
        r = 1
        push!(atoms,Atom(1,"H",[0,0,0*r]))
        push!(atoms,Atom(2,"H",[0,0,1*r]))
        push!(atoms,Atom(3,"H",[0,0,2*r]))
        push!(atoms,Atom(4,"H",[0,0,3*r]))
        push!(atoms,Atom(5,"H",[0,0,4*r]))
        push!(atoms,Atom(6,"H",[0,0,5*r]))
        push!(atoms,Atom(7,"H",[0,0,6*r]))
        push!(atoms,Atom(8,"H",[0,0,7*r]))
        push!(atoms,Atom(9,"H",[0,0,8*r]))
        push!(atoms,Atom(10,"H",[0,0,9*r]))
        push!(atoms,Atom(11,"H",[0,0,10*r]))
        push!(atoms,Atom(12,"H",[0,0,11*r]))

        clusters    = [(1:2),(3:4),(5:6),(7:8)]
        init_fspace = [(1,1),(1,1),(1,1),(1,1)]
        na = 4
        nb = 4
        clusters    = [(1:2),(3:4),(5:6),(7:8),(9:10),(11:12)]
        init_fspace = [(1,1),(1,1),(1,1),(1,1),(1,1),(1,1)]
        na = 6
        nb = 6
    elseif false
        push!(atoms,Atom(1,"H",[-1.30,0,0.00]))
        push!(atoms,Atom(2,"H",[-1.30,0,1.00]))
        push!(atoms,Atom(3,"H",[ 0.00,0,0.00]))
        push!(atoms,Atom(4,"H",[ 0.00,0,1.00]))
        push!(atoms,Atom(5,"H",[ 1.33,0,0.00]))
        push!(atoms,Atom(6,"H",[ 1.30,0,1.00]))

        clusters    = [(1:2),(3:4),(5:6)]
        init_fspace = [(1,1),(1,1),(1,1)]
        na = 3
        nb = 3
    elseif true
        atoms = generate_H_ring(8,rad)
        clusters    = [(1:2),(3:4),(5:6),(7:8)]
        init_fspace = [(1,1),(1,1),(1,1),(1,1)]
        na = 4
        nb = 4

        atoms = generate_H_ring(10,rad)
        clusters    = [(1:2),(3:4),(5:6),(7:8),(9:10)]
        init_fspace = [(1,1),(1,1),(1,1),(1,1),(1,1)]
        na = 5
        nb = 5

        atoms = generate_H_ring(12,rad)
        clusters    = [(1:2),(3:4),(5:6),(7:8),(9:10),(11:12)]
        init_fspace = [(1,1),(1,1),(1,1),(1,1),(1,1),(1,1)]
        clusters    = [(1:4),(5:8),(9:12)]
        init_fspace = [(2,2),(2,2),(2,2)]
        na = 6
        nb = 6
    end

    basis = "6-31g"
    basis = "sto-3g"
    mol     = Molecule(0,1,atoms,basis)

    # get integrals
    mf = FermiCG.pyscf_do_scf(mol)
    push!(out_rhf, mf.e_tot)
    nbas = size(mf.mo_coeff)[1]
    ints = FermiCG.pyscf_build_ints(mol,mf.mo_coeff, zeros(nbas,nbas));
    #e_fci, d1_fci, d2_fci = FermiCG.pyscf_fci(ints, na, nb, conv_tol=1e-10,max_cycle=100, nroots=2)

    #run fci with pyscf
    if false
        pyscf = pyimport("pyscf")
        fci = pyimport("pyscf.fci")
        mp = pyimport("pyscf.mp")
        mp2 = mp.MP2(mf)
        push!(out_mp2, mp2.kernel()[1])
        cisolver = pyscf.fci.direct_spin1.FCI()
        cisolver.max_cycle = 100
        cisolver.conv_tol = 1e-10
        nelec = na + nb
        norb = size(ints.h1)[1]
        e_fci, ci = cisolver.kernel(ints.h1, ints.h2, norb , nelec, ecore=0, nroots = 1, verbose=100)
        e_fci = min(e_fci...)
        @printf(" FCI Energy: %12.8f\n", e_fci)

        push!(out_fci, e_fci + ints.h0)
    end
    push!(out_radii, rad)

    # localize orbitals
    C = mf.mo_coeff
    Cl = FermiCG.localize(mf.mo_coeff,"lowdin",mf)
    FermiCG.pyscf_write_molden(mol,Cl,filename="lowdin.molden")
    S = FermiCG.get_ovlp(mf)
    U =  C' * S * Cl
    println(" Rotate Integrals")
    flush(stdout)
    ints = FermiCG.orbital_rotation(ints,U)
    println(" done.")
    flush(stdout)

    # define clusters

    #fname = "job.scr"
    #fid = h5open(fname, "r")
    #Cl = read(fid["mo_coeffs"])
    #close(fid)
    #ints = FermiCG.pyscf_build_ints(mol,Cl, zeros(nbas,nbas));

    clusters = [Cluster(i,collect(clusters[i])) for i = 1:length(clusters)]
    display(clusters)

    rdm1 = zeros(size(ints.h1))
    #e_cmf, U, Da, Db  = FermiCG.cmf_oo(ints, clusters, init_fspace, rdm1,
    #                                   max_iter_oo=40, verbose=0, gconv=1e-6, method="gd", alpha=1e-1)
    #ints = FermiCG.orbital_rotation(ints,U)

    e_cmf, U, Da, Db  = FermiCG.cmf_oo(ints, clusters, init_fspace, rdm1,
                                       max_iter_oo=40, verbose=0, gconv=1e-6, method="bfgs")
    FermiCG.pyscf_write_molden(mol,Cl*U,filename="cmf.molden")
    ints = FermiCG.orbital_rotation(ints,U)

    #fname = "job.scr"
    #fid = h5open(fname, "w")
    #fid["mo_coeffs"] = Cl*U
    #close(fid)

    #continue
    #cmf_out = FermiCG.cmf_ci(ints, clusters, init_fspace, rdm1, verbose=1)
    #e_ref = cmf_out[1]

    e_ref = e_cmf - ints.h0

    max_roots = 20
    # build Hamiltonian, cluster_basis and cluster ops
    #display(Da)
    #cluster_bases = FermiCG.compute_cluster_eigenbasis(ints, clusters, verbose=2, max_roots=max_roots)
    cluster_bases = FermiCG.compute_cluster_eigenbasis(ints, clusters, verbose=0, max_roots=max_roots,
                                                       init_fspace=init_fspace, rdm1a=Da, rdm1b=Db)
    clustered_ham = FermiCG.extract_ClusteredTerms(ints, clusters)
    cluster_ops = FermiCG.compute_cluster_ops(cluster_bases, ints);

    p_spaces = Vector{FermiCG.ClusterSubspace}()
    q_spaces = Vector{FermiCG.ClusterSubspace}()

    #ci_vector = FermiCG.TuckerState(clusters)
    #FermiCG.add_fockconfig!(ci_vector, [(1,1),(1,1),(1,1)])
    #FermiCG.add_fockconfig!(ci_vector, [(2,1),(0,1),(1,1)])
    #FermiCG.add_fockconfig!(ci_vector, [(0,1),(2,1),(1,1)])
    #FermiCG.add_fockconfig!(ci_vector, [(2,1),(0,1),(1,1)])

    #FermiCG.expand_each_fock_space!(ci_vector, cluster_bases)

    for ci in clusters
        tss = FermiCG.ClusterSubspace(ci)
        tss[(2,2)] = 1:1
        #tss[(2,1)] = 1:1
        #tss[(1,2)] = 1:1
        #tss[(0,1)] = 1:1
        #tss[(1,0)] = 1:1
        push!(p_spaces, tss)
    end

    for tssp in p_spaces
        tss = FermiCG.get_ortho_compliment(tssp, cluster_bases[tssp.cluster.idx])
        push!(q_spaces, tss)
    end

    println(" ================= Cluster P Spaces ===================")
    display.(p_spaces)
    println(" ================= Cluster Q Spaces ===================")
    display.(q_spaces)

    nroots = 1
    ci_vector = FermiCG.TuckerState(clusters, p_spaces, na, nb, nroots=nroots)

    ref_vector = deepcopy(ci_vector)
    if true
        for ci in clusters
            tmp_spaces = copy(p_spaces)
            tmp_spaces[ci.idx] = q_spaces[ci.idx]
            FermiCG.add!(ci_vector, FermiCG.TuckerState(clusters, tmp_spaces, na, nb))
        end
    end
    if true
        for ci in clusters
            for cj in clusters
                ci.idx < cj.idx || continue
                tmp_spaces = copy(p_spaces)
                tmp_spaces[ci.idx] = q_spaces[ci.idx]
                tmp_spaces[cj.idx] = q_spaces[cj.idx]
                FermiCG.add!(ci_vector, FermiCG.TuckerState(clusters, tmp_spaces, na, na))
            end
        end
    end
    if false
        for ci in clusters
            for cj in clusters
                for ck in clusters
                    ci.idx < cj.idx || continue
                    cj.idx < ck.idx || continue
                    tmp_spaces = copy(p_spaces)
                    tmp_spaces[ci.idx] = q_spaces[ci.idx]
                    tmp_spaces[cj.idx] = q_spaces[cj.idx]
                    tmp_spaces[ck.idx] = q_spaces[ck.idx]
                    FermiCG.add!(ci_vector, FermiCG.TuckerState(clusters, tmp_spaces, na, na))
                end
            end
        end
    end
    if false
        for ci in clusters
            for cj in clusters
                for ck in clusters
                    for cl in clusters
                        ci.idx < cj.idx || continue
                        cj.idx < ck.idx || continue
                        ck.idx < cl.idx || continue
                        tmp_spaces = copy(p_spaces)
                        tmp_spaces[ci.idx] = q_spaces[ci.idx]
                        tmp_spaces[cj.idx] = q_spaces[cj.idx]
                        tmp_spaces[ck.idx] = q_spaces[ck.idx]
                        tmp_spaces[cl.idx] = q_spaces[cl.idx]
                        FermiCG.add!(ci_vector, FermiCG.TuckerState(clusters, tmp_spaces, na, nb))
                    end
                end
            end
        end
    end

    #S = FermiCG.dot(q_vector, q_vector)
    ##display(S - I)
    #@test isapprox(S-I, zeros(size(S)), atol=1e-10)

    # initialize with eye
    FermiCG.set_vector!(ref_vector, Matrix(1.0I, length(ref_vector),nroots))
    FermiCG.set_vector!(ci_vector, Matrix(1.0I, length(ci_vector),nroots))

    #FermiCG.randomize!(ci_vector, scale=1e-4)

    #FermiCG.orthogonalize!(ci_vector)

    S = FermiCG.dot(ci_vector, ci_vector)
    @test isapprox(S-I, zeros(size(S)), atol=1e-12)

    if true
        #FermiCG.print_fock_occupations(ci_vector)
        println(" Length of CI Vector: ", length(ci_vector))
        @time e_nb2 = FermiCG.tucker_ci_solve!(ci_vector, cluster_ops, clustered_ham)
        push!(out_nb2, min(e_nb2[1]...) + ints.h0)
        #@time FermiCG.tucker_ci_solve!(ci_vector, cluster_ops, clustered_ham)
        FermiCG.print_fock_occupations(ci_vector)
        #display(ci_vector, thresh=-1)
        #display(FermiCG.get_vector(ci_vector))
    end

    if true
        #FermiCG.print_fock_occupations(ref_vector)
        e_ref = FermiCG.tucker_ci_solve!(ref_vector, cluster_ops, clustered_ham)
        println(" Reference State:" )
        push!(out_ref,  min(e_ref[1]...) + ints.h0)
        FermiCG.print_fock_occupations(ref_vector)
        #FermiCG.print_fock_occupations(ci_vector)

        @time e_cepa = FermiCG.tucker_cepa_solve!(ref_vector, ci_vector, cluster_ops, clustered_ham)
        println(e_cepa[1])
        push!(out_cepa, e_cepa[1] + ints.h0)
        FermiCG.print_fock_occupations(ci_vector)
        display(ci_vector)
    end
end

#end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
function chebyQuad(f,x)
    n = length(x)
    f[:] = zeros(n);
    for i = 2:2:n
        f[i] = n/(i*i-1)
    end
    for l = 1:n
        factt = 4.0*x[l] - 2.0
        ti2 = 1.0
        ti1 = 0.5*factt
        f[1] = ti1 + f[1]
        for i = 2:n
            ti = factt*ti1-ti2
            f[i] = ti + f[i]
            ti2 = ti1
            ti1 = ti
        end
    end
    return nothing
end

function chebyQuadJac(J,x)
    n = length(x)
    for j = 1:n
        factt = 4.0*x[j] - 2.0
        ti2 = 1.0
        ti1 = 0.5*factt
        tabli2 = 0.0
        tabli1 = 2.0
        J[1,j] = tabli1
        for i = 2:n
            ti = factt*ti1-ti2
            tabli = 4.0*ti1+factt*tabli1-tabli2
            J[i,j] = tabli
            ti2 = ti1
            ti1 = ti
            tabli2 = tabli1
            tabli1 = tabli
        end
    end
    return nothing
end


text:
using Plots,FFTW,JLD2

function main()
    sampling_rate = 0.0
    window_size = 64*1024;
    w_m_vals::Array{Float64} = []

    jldopen(ARGS[1], "r") do f
        sampling_rate = f["samplingrate"]
        w_m_vals = f["w_m"]
    end

    f_t = open("time_fft.dat", "w")
    f_i = open("I_fft.dat", "w")
    f_q = open("Q_fft.dat", "w")
    f_am_p = open("am_fft.dat", "w")

    duration = 2500000
    dt = 1/sampling_rate
    for t in 1:duration
        I::Float64 = 0
        Q::Float64 = 0
        spectrum = fft(w_m_vals[t:window_size+t-1])|>fftshift
        I = real(spectrum[32*1024+626])/window_size
        Q = -imag(spectrum[32*1024+626])/window_size
        A = 2.0*sqrt(I^2+Q^2)
        write(f_t, t*dt)
        write(f_i, I)
        write(f_q, Q)
        write(f_am_p, A)
    end
    close(f_am_p)
    close(f_q)
    close(f_i)
    close(f_t)

    timeline::Array{Float64}=[]
    power::Array{Float64}=[]
    open("time_fft.dat","r") do f_t
        open("am_fft.dat","r") do f_am
            for step in 1:1000*1000
                t = read(f_t, Float64)
                p = read(f_am, Float64)
                push!(timeline, t)
                push!(power, p)
            end
        end
    end
    plot_graph=plot(timeline, power, fmt=:png)
    savefig(plot_graph, "demodulate_fft_am.png")

    jldopen("am_demodulate_fft.jld2","w") do file
        file["samplingrate"] = sampling_rate

        vals = []
        open("time.dat","r") do f
            while !eof(f_t)
                push!(vals, read(f, Float64))
            end
        end
        file["times"] = vals

        vals = []
        open("I.dat","r") do f
            while !eof(f)
                push!(vals, read(f, Float64))
            end
        end
        file["I"] = vals

        vals = []
        open("Q.dat","r") do f
            while !eof(f)
                push!(vals, read(f, Float64))
            end
        end
        file["Q"] = vals

        vals = []
        open("am.dat","r") do f
            while !eof(f)
                push!(vals, read(f, Float64))
            end
        end
        file["demodulated_signal"] = vals

    end

    rm("time_fft.dat")
    rm("I_fft.dat")
    rm("Q_fft.dat")
    rm("am_fft.dat")

end

main()

max_stars_repo_name: seisyuu-hantatsushi/julia

max_issues_repo_name: seisyuu-hantatsushi/julia

max_forks_repo_name: seisyuu-hantatsushi/julia


text:
using SnakeLadder
using Base.Test

include("squaresubstitutiontest.jl")

max_stars_repo_name: qlete/markov-decision

max_issues_repo_name: qlete/markov-decision

max_forks_repo_name: qlete/markov-decision


text:
import Nemo
import Nemo: acb, arb, AcbField,  onei
export NemoCharFunc
#support for Nemo

struct NemoCharFunc{MT} <: CharFunc{MT,acb} #model type, return type (e.g. Complex or acb)
    model::MT
    field::AcbField
end

@inline CharFuncPricing.model(cf::NemoCharFunc) = cf.model
@inline CharFuncPricing.oneim(cf::NemoCharFunc) = onei(cf.field)
@inline CharFuncPricing.const_pi(cf::NemoCharFunc) = real(Nemo.const_pi(cf.field))
@inline Base.expm1(x::acb) = exp(x) - 1 #for Nemo - now in latest Nemo 0.19.2
@inline Base.sign(x::arb) = x < 0 ? -1 : 1


text:
"""
documentation tools.
"""
module Doc

using Comonicon
using Pkg
using LiveServer: servedocs

"""
build documentation.

# Args

- `path`: path of the project.
- `args...`: command line arguments for `docs/make.jl`
"""
@cast function build(path::String=pwd(), args...)
    project_dir = dirname(Base.current_project(path))
    docs_dir = joinpath(project_dir, "docs")

    # with project
    old_project = Pkg.project().path
    Pkg.activate(docs_dir)
    Pkg.develop(path=project_dir)
    Pkg.instantiate()
    # with ARGS
    old_ARGS = copy(ARGS)
    copy!(ARGS, collect(args))

    Main.include(joinpath(docs_dir, "make.jl"))
    copy!(ARGS, old_ARGS)
    Pkg.activate(old_project)
    return
end

"""
serve documentation.

# Options

- `--foldername <name>`: specify the name of the content folder if different than "docs".
- `--literate <path>`: is the path to the folder containing the literate scripts, if
    left empty, it will be assumed that they are in docs/src.

# Flags

- `--verbose`: show verbose log.
- `--doc-env`: is a boolean switch to make the server start by activating the
    doc environment or not (i.e. the Project.toml in docs/).
- `-l, --launch-browser`: open documentation in default browser.
"""
@cast function serve(;verbose::Bool=false, literate="", foldername="docs", launch_browser::Bool=false)
    servedocs(;verbose=verbose, literate=literate, doc_env=true, foldername=abspath(foldername), launch_browser=launch_browser)
end

end

@cast Doc


text:
function statisticalModel(
    nDays,
    GPP_daily_mu,
    GPP_daily_sigma,
    ER_daily_mu,
    ER_daily_sigma,
    k600_meanlog,
    k600_sdlog,
    K02_conv,
    err_obs_iid_sigma_scale,
    err_proc_iid_sigma_scale)

    return eval(quote
        @model function f(::Type{T} = Float64; DOdata, ODE, times) where {T}

            # Priors on daily GPP and ER, k600, observation error, and process error
            GPP_daily ~ MvNormal(FillArrays.Fill($GPP_daily_mu, $nDays), $GPP_daily_sigma)
            ER_daily ~ MvNormal(FillArrays.Fill($ER_daily_mu, $nDays), $ER_daily_sigma)
            k600 ~ LogNormal($k600_meanlog, $k600_sdlog)
            kGas = $K02_conv*k600
            err_obs_iid_sigma ~ truncated(Cauchy(0, $err_obs_iid_sigma_scale), 0, Inf)
            err_proc_iid_sigma ~ truncated(Cauchy(0, $err_proc_iid_sigma_scale), 0, Inf)

            # Specify model initial condition and parameters
            p = Dict(
                "PAR" => ODE.p["PAR"],
                "OSat" => ODE.p["OSat"],
                "day" => ODE.p["day"],
                "GPP_daily"=> GPP_daily,
                "ER_daily" => ER_daily,
                "kGas" => kGas,
                "Q" => ODE.p["Q"],
                "V" => ODE.p["V"])

            DO_true = Vector{T}(undef, (length(times)))
            DO_true[1] ~ Normal(DOdata[1], err_obs_iid_sigma)

            for i = 2:length(times)
                problem = remake(ODE, p = p, u0 = [DO_true[i-1]], tspan = (times[i-1], times[i]))
                predicted = solve(problem, alg = Euler(), dt = 0.1)
                DO_true[i] ~ Normal(predicted(times[i])[1], err_proc_iid_sigma)
                DOdata[i] ~ Normal(DO_true[i], err_obs_iid_sigma)
            end
        end
    end)
end

max_stars_repo_name: elizabethjmohr/mesocosm-metabolism

max_issues_repo_name: elizabethjmohr/mesocosm-metabolism

max_forks_repo_name: elizabethjmohr/mesocosm-metabolism


text:
using TrajectoryOptimization
using RobotDynamics
using StaticArrays
using LinearAlgebra

const TO = TrajectoryOptimization
const RD = RobotDynamics


text:
using LinearAlgebra
import Base: *
import LinearAlgebra: mul!

struct BroydenOperator{R,TH}
    H::TH
    theta_bar::R
end

function BroydenOperator(::Type{T}, n::Integer, theta_bar=T(0.2)) where T
    H = Matrix{T}(I, n, n)
    BroydenOperator(H, theta_bar)
end

function BroydenOperator(x::AbstractVector{T}, theta_bar=T(0.2)) where T
    BroydenOperator(T, length(x), theta_bar)
end

_sign(x::R) where {R} = x == 0 ? R(1) : sign(x)

function update!(L::BroydenOperator{R,TH}, s, y) where {R,TH}
    Hy = L.H * y
    sH = s' * L.H
    delta = dot(Hy, s) / norm(s)^2
    theta = if abs(delta) >= L.theta_bar
        R(1)
    else
        (1 - _sign(delta) * L.theta_bar) / (1 - delta)
    end
    L.H .+= (s - Hy) / dot(s, (1 / theta - 1) * s + Hy) * sH
end

function reset!(L::BroydenOperator{R,TH}) where {R,TH}
    L.H .= 0
    L.H[diagind(L.H)] .= 1
end

function (*)(L::BroydenOperator, v)
    w = similar(v)
    return mul!(w, L, v)
end

function mul!(d::T, L::BroydenOperator, v::T) where {T}
    mul!(d, L.H, v)
    return d
end

Base.@kwdef struct Broyden{R}
    theta_bar::R = 0.2
end

acceleration_style(::Type{<:Broyden}) = QuasiNewtonStyle()

function initialize(broyden::Broyden, x::AbstractVector{R}) where R
    return BroydenOperator(x, broyden.theta_bar)
end


text:
import HTTP: request
import JSON: parse

"""
    symbols(conn::Connection; format::String="JSON")

This call returns an array of symbols the Investors Exchange supports for
trading. This list is updated daily as of 7:45 a.m. ET. Symbols may be added or
removed by the Investors Exchange after the list was produced
"""
function symbols(conn::Connection; format::String = "JSON")
    url = conn.url * conn.version * "/ref-data/symbols/"
    r = request("GET", url;
        query = Dict("token" => conn.token, "format" => format)
    )
    return parse(String(r.body))
end

max_stars_repo_name: GodolphinCapitalManagement/Iex

max_issues_repo_name: GodolphinCapitalManagement/Iex

max_forks_repo_name: GodolphinCapitalManagement/Iex


text:
""" Extends the matrices u,v,η,sst with a halo of ghost points for boundary conditions."""
function add_halo(  u::Array{T,2},
                    v::Array{T,2},
                    η::Array{T,2},
                    sst::Array{T,2},
                    S::ModelSetup) where {T<:AbstractFloat}

    @unpack nx,ny,nux,nuy,nvx,nvy = S.grid
    @unpack halo,haloη,halosstx,halossty = S.grid

    # Add zeros to satisfy kinematic boundary conditions
    u = cat(zeros(T,halo,nuy),u,zeros(T,halo,nuy),dims=1)
    u = cat(zeros(T,nux+2*halo,halo),u,zeros(T,nux+2*halo,halo),dims=2)

    v = cat(zeros(T,halo,nvy),v,zeros(T,halo,nvy),dims=1)
    v = cat(zeros(T,nvx+2*halo,halo),v,zeros(T,nvx+2*halo,halo),dims=2)

    η = cat(zeros(T,haloη,ny),η,zeros(T,haloη,ny),dims=1)
    η = cat(zeros(T,nx+2*haloη,haloη),η,zeros(T,nx+2*haloη,haloη),dims=2)

    sst = cat(zeros(T,halosstx,ny),sst,zeros(T,halosstx,ny),dims=1)
    sst = cat(zeros(T,nx+2*halosstx,halossty),sst,zeros(T,nx+2*halosstx,halossty),dims=2)

    # SCALING
    @unpack scale,scale_sst = S.constants
    u *= scale
    v *= scale
    sst *= scale_sst

    ghost_points!(u,v,η,S)
    ghost_points_sst!(sst,S)
    return u,v,η,sst
end

"""Cut off the halo from the prognostic variables."""
function remove_halo(   u::Array{T,2},
                        v::Array{T,2},
                        η::Array{T,2},
                        sst::Array{T,2},
                        S::ModelSetup) where {T<:AbstractFloat}

    @unpack halo,haloη,halosstx,halossty = S.grid
    @unpack scale_inv,scale_sst = S.constants

    # undo scaling as well
    @views ucut = scale_inv*u[halo+1:end-halo,halo+1:end-halo]
    @views vcut = scale_inv*v[halo+1:end-halo,halo+1:end-halo]
    @views ηcut = η[haloη+1:end-haloη,haloη+1:end-haloη]
    @views sstcut = sst[halosstx+1:end-halosstx,halossty+1:end-halossty]/scale_sst

    return ucut,vcut,ηcut,sstcut
end

""" Copy ghost points for u from inside to the halo in the nonperiodic case. """
function ghost_points_u_nonperiodic!(u::AbstractMatrix{T},one_minus_α::T) where T
    n,m = size(u)

    # tangential boundary condition
    @inbounds for i in 1:n
        u[i,1] = one_minus_α*u[i,4]
        u[i,2] = one_minus_α*u[i,3]
        u[i,end-1] = one_minus_α*u[i,end-2]
        u[i,end] = one_minus_α*u[i,end-3]
    end
end

""" Copy ghost points for u from inside to the halo in the periodic case. """
function ghost_points_u_periodic!(::Type{Tcomm},u::AbstractMatrix{T},one_minus_α::T) where {Tcomm,T}
    n,m = size(u)

    # periodic bc - mimic MPI communication with (reduced precision) type Tcomm
    @inbounds for j in 1:m
        u[1,j] = T(Tcomm(u[end-3,j]))
        u[2,j] = T(Tcomm(u[end-2,j]))
        u[end-1,j] = T(Tcomm(u[3,j]))
        u[end,j] = T(Tcomm(u[4,j]))
    end

    # tangential bc
    # ghost_points_u_nonperiodic!(u,one_minus_α)
end

""" Copy ghost points for v from inside to the halo in the nonperiodic case. """
function ghost_points_v_nonperiodic!(v::AbstractMatrix{T},one_minus_α::T) where T
    n,m = size(v)

    # tangential boundary condition
    @inbounds for j in 1:m
        v[1,j] = one_minus_α*v[4,j]
        v[2,j] = one_minus_α*v[3,j]
        v[end-1,j] = one_minus_α*v[end-2,j]
        v[end,j] = one_minus_α*v[end-3,j]
    end
end

""" Copy ghost points for v from inside to the halo in the periodic case. """
function ghost_points_v_periodic!(::Type{Tcomm},v::Array{T,2}) where {Tcomm,T}
    n,m = size(v)

    # mimic MPI communication with (reduced precision) type Tcomm
    @inbounds for j in 1:m
        v[1,j] = T(Tcomm(v[end-3,j]))
        v[2,j] = T(Tcomm(v[end-2,j]))
        v[end-1,j] = T(Tcomm(v[3,j]))
        v[end,j] = T(Tcomm(v[4,j]))
    end
end

""" Copy ghost points for η from inside to the halo in the nonperiodic case. """
function ghost_points_η_nonperiodic!(η::AbstractMatrix)
    n,m = size(η)

    # assume no gradients of η across solid boundaries
    # the 4 corner points are copied twice, but it's faster!
    @inbounds for j in 1:m
        η[1,j] = η[2,j]
        η[end,j] = η[end-1,j]
    end

    @inbounds for i in 1:n
        η[i,1] = η[i,2]
        η[i,end] = η[i,end-1]
    end
end

""" Copy ghost points for η from inside to the halo in the periodic case. """
function ghost_points_η_periodic!(::Type{Tcomm},η::Array{T,2}) where {Tcomm,T}
    n,m = size(η)

    # corner points are copied twice, but it's faster!
    # mimic MPI communication with (reduced precision) type Tcomm
    @inbounds for j in 1:m
        η[1,j] = T(Tcomm(η[end-1,j]))
        η[end,j] = T(Tcomm(η[2,j]))
    end

    @inbounds for i in 1:n
        η[i,1] = η[i,2]
        η[i,end] = η[i,end-1]
    end
end

""" Copy ghost points for η from inside to the halo in the nonperiodic case. """
function ghost_points_sst_nonperiodic!(sst::AbstractMatrix,
                                        halosstx::Int,
                                        halossty::Int)
    n,m = size(sst)

    # assume no gradients of η across solid boundaries
    # the 4 corner points are copied twice, but it's faster!
    @inbounds for j in 1:m
        for i ∈ 1:halosstx
            sst[i,j] = sst[halosstx+1,j]
            sst[end-i+1,j] = sst[end-halosstx,j]
        end
    end

    @inbounds for j ∈ 1:halossty
        for i in 1:n
            sst[i,j] = sst[i,halossty+1]
            sst[i,end-j+1] = sst[i,end-halossty]
        end
    end
end

""" Copy ghost points for η from inside to the halo in the periodic case. """
function ghost_points_sst_periodic!(::Type{Tcomm},
                                    sst::Array{T,2},
                                    halosstx::Int,
                                    halossty::Int) where {Tcomm,T}
    n,m = size(sst)

    # corner points are copied twice, but it's faster!
    @inbounds for j in 1:m
        for i ∈ 1:halosstx
            sst[i,j] = T(Tcomm(sst[end-2*halosstx+i,j]))
            sst[end-halosstx+i,j] = T(Tcomm(sst[halosstx+i,j]))
        end
    end

    @inbounds for j ∈ 1:halossty
        for i in 1:n
            sst[i,j] = sst[i,halossty+1]
            sst[i,end-j+1] = sst[i,end-halossty]
        end
    end
end

"""Decide on boundary condition P.bc which ghost point function to execute."""
function ghost_points!( u::AbstractMatrix,
                        v::AbstractMatrix,
                        η::AbstractMatrix,
                        S::ModelSetup)

    @unpack bc,Tcomm = S.parameters
    @unpack one_minus_α = S.constants

    if bc == "periodic"
        ghost_points_u_periodic!(Tcomm,u,one_minus_α)
        ghost_points_v_periodic!(Tcomm,v)
        ghost_points_η_periodic!(Tcomm,η)
    else
        ghost_points_u_nonperiodic!(u,one_minus_α)
        ghost_points_v_nonperiodic!(v,one_minus_α)
        ghost_points_η_nonperiodic!(η)
    end
end

"""Decide on boundary condition P.bc which ghost point function to execute."""
function ghost_points_uv!( u::AbstractMatrix,
                        v::AbstractMatrix,
                        S::ModelSetup)

    @unpack bc,Tcomm = S.parameters
    @unpack one_minus_α = S.constants

    if bc == "periodic"
        ghost_points_u_periodic!(Tcomm,u,one_minus_α)
        ghost_points_v_periodic!(Tcomm,v)
    else
        ghost_points_u_nonperiodic!(u,one_minus_α)
        ghost_points_v_nonperiodic!(v,one_minus_α)
    end
end

"""Decide on boundary condition P.bc which ghost point function to execute."""
function ghost_points_η!(   η::AbstractMatrix,
                            S::ModelSetup)

    @unpack bc, Tcomm = S.parameters

    if bc == "periodic"
        ghost_points_η_periodic!(Tcomm,η)
    else
        ghost_points_η_nonperiodic!(η)
    end
end

"""Decide on boundary condition P.bc which ghost point function to execute."""
function ghost_points_sst!(sst::AbstractMatrix,S::ModelSetup)

    @unpack bc,Tcomm = S.parameters
    @unpack halosstx,halossty = S.grid

    if bc == "periodic"
        ghost_points_sst_periodic!(Tcomm,sst,halosstx,halossty)
    else
        ghost_points_sst_nonperiodic!(sst,halosstx,halossty)
    end
end


text:
# This is free and unencumbered software released into the public domain.

module DRYlib

################################################################################
# Type Definitions

# Boolean (true or false)
const Bool = Core.Bool

# Character (Unicode code point)
const Char = Core.Char

# Complex number (arbitrary size)
const Complex = Base.Complex{Base.MPFR.BigFloat}

# Floating-point number (native size)
const Float = Core.Float64

# Floating-point number (32-bit single-precision)
const Float32 = Core.Float32

# Floating-point number (64-bit double-precision)
const Float64 = Core.Float64

# Integer number (native size)
const Int = if (Core.Int === Core.Int64) Core.Int64 else Core.Int32 end

# Integer number (8-bit)
const Int8 = Core.Int8

# Integer number (16-bit)
const Int16 = Core.Int16

# Integer number (32-bit)
const Int32 = Core.Int32

# Integer number (64-bit)
const Int64 = Core.Int64

# Integer number (128-bit)
const Int128 = Core.Int128

# Integer number (arbitrary size)
const Integer = Base.GMP.BigInt

# Natural number (arbitrary size)
const Natural = Integer

# Rational number (arbitrary size)
const Rational = Base.Rational{Base.GMP.BigInt}

# Real number (arbitrary size)
const Real = Base.MPFR.BigFloat

# Machine word (native size)
const Word = if (Core.UInt === Core.UInt64) Core.UInt64 else Core.UInt32 end

# Machine word (8-bit)
const Word8 = Core.UInt8

# Machine word (16-bit)
const Word16 = Core.UInt16

# Machine word (32-bit)
const Word32 = Core.UInt32

# Machine word (64-bit)
const Word64 = Core.UInt64

# Machine word (128-bit)
const Word128 = Core.UInt128

################################################################################

end # module DRYlib

#println("OK") # julia drylib.jl


text:
"""
    struct WendlandC2{T} <: SPHKernel
        n_neighbours::Int64
        norm_1D::T
        norm_2D::T
        norm_3D::T
    end
"""
struct WendlandC2{T} <: SPHKernel
    n_neighbours::Int64
    norm_1D::T
    norm_2D::T
    norm_3D::T
end

"""
    WendlandC2(T::DataType=Float64, n_neighbours::Integer=100)

Set up a `WendlandC2` with a given DataType `T`.
"""
WendlandC2(T::DataType=Float64, n_neighbours::Integer=100) = WendlandC2{T}(n_neighbours, T(5.0/4.0), T(7.0/π), T(21.0/(2π)))

"""
    kernel_value_1D(kernel::WendlandC2{T}, u::Real, h_inv::Real) where T

Evaluate WendlandC2 spline at position ``u = \\frac{x}{h}``.
"""
@inline function kernel_value_1D(kernel::WendlandC2{T}, u::Real, h_inv::Real) where T

    if u < 1
        n = kernel.norm_1D * h_inv
        t1 = 1 - u
        t3 = t1*t1*t1
        return ( t3 * ( 1 + 3u )) * n |> T
    else
        return 0.0 |> T
    end

end

"""
    kernel_deriv_1D(kernel::WendlandC2{T}, u::Real, h_inv::Real) where T

Evaluate the derivative of the WendlandC2 spline at position ``u = \\frac{x}{h}``.
"""
@inline function kernel_deriv_1D(kernel::WendlandC2{T}, u::Real, h_inv::Real) where T

    if u < 1
        n = kernel.norm_1D * h_inv^2
        t1 = 1 - u
        return ( -12u * t1  ) * n |> T
    else
        return 0.0 |> T
    end

end

"""
    bias_correction_1D(kernel::WendlandC2{T}, density::Real, m::Real, h_inv::Real) where T

Corrects the density estimate for the kernel bias. See Dehnen&Aly 2012, eq. 18+19.
"""
@inline function bias_correction_1D(kernel::WendlandC2{T}, density::Real, m::Real, h_inv::Real) where T
    n = kernel.norm_1D * h_inv^3
    wc_correction = T(0.0294) * ( kernel.n_neighbours * T(0.01) )^T(-0.977) * m * n
    return density - wc_correction |> T
end

"""
    kernel_value_2D(kernel::WendlandC2{T}, u::Real, h_inv::Real) where T

Evaluate WendlandC2 spline at position ``u = \\frac{x}{h}``.
"""
@inline function kernel_value_2D(kernel::WendlandC2{T}, u::Real, h_inv::Real) where T

    if u < 1
        n = kernel.norm_2D * h_inv^2
        t1 = 1 - u
        t4 = t1 * t1 * t1 * t1
        return ( t4 * ( 1 + 4u ) ) * n |> T
    else
        return 0.0 |> T
    end

end

"""
    kernel_deriv_2D(kernel::WendlandC2{T}, u::Real, h_inv::Real) where T

Evaluate the derivative of the WendlandC2 spline at position ``u = \\frac{x}{h}``.
"""
@inline function kernel_deriv_2D(kernel::WendlandC2{T}, u::Real, h_inv::Real) where T

    if u < 1
        n = kernel.norm_2D * h_inv^3
        t1 = 1 - u
        t3 = t1 * t1 * t1
        return ( -20u * t3 ) * n |> T
    else
        return 0.0 |> T
    end

end

"""
    bias_correction_2D(kernel::WendlandC2{T}, density::Real, m::Real, h_inv::Real) where T

Corrects the density estimate for the kernel bias. See Dehnen&Aly 2012, eq. 18+19.
"""
@inline function bias_correction_2D(kernel::WendlandC2{T}, density::Real, m::Real, h_inv::Real) where T
    n = kernel.norm_2D * h_inv^3
    wc_correction = T(0.0294) * ( kernel.n_neighbours * T(0.01) )^T(-0.977) * m * n
    return density - wc_correction |> T
end

"""
    kernel_value_3D(kernel::WendlandC2{T}, u::Real, h_inv::Real) where T

Evaluate WendlandC2 spline at position ``u = \\frac{x}{h}``.
"""
@inline function kernel_value_3D(kernel::WendlandC2{T}, u::Real, h_inv::Real) where T

    if u < 1
        n  = kernel.norm_3D * h_inv^3
        t1 = 1 - u
        t4 = t1 * t1 * t1 * t1
        return ( t4 * ( 1 + 4u ) ) * n |> T
    else
        return 0.0 |> T
    end

end

"""
    kernel_deriv_3D(kernel::WendlandC2{T}, u::Real, h_inv::Real) where T

Evaluate the derivative of the WendlandC2 spline at position ``u = \\frac{x}{h}``.
"""
@inline function kernel_deriv_3D(kernel::WendlandC2{T}, u::Real, h_inv::Real) where T

    if u < 1
        n = kernel.norm_3D * h_inv^2 * h_inv^2
        t1 = 1 - u
        t3 = t1 * t1 * t1
        return ( -20u * t3 ) * n |> T
    else
        return 0.0 |> T
    end

end

"""
    bias_correction_3D(kernel::WendlandC2, density::Real, m::Real, h_inv::Real)

Corrects the density estimate for the kernel bias. See Dehnen&Aly 2012, eq. 18+19.
"""
@inline function bias_correction_3D(kernel::WendlandC2{T}, density::Real, m::Real, h_inv::Real) where T

    n = kernel.norm_3D * h_inv^3
    wc_correction = T(0.0294) * ( kernel.n_neighbours * T(0.01) )^T(-0.977) * m * n
    return density - wc_correction |> T
end


text:
module contactForceLaw_ballWithBox_45Deg

using Modia3D
using Modia3D.StaticArrays
import Modia3D.ModiaMath

vmat1 = Modia3D.Material(color="LightBlue" , transparency=0.5)    # material of SolidFileMesh
vmat2 = deepcopy(vmat1)                                           # material of convex decomposition of SolidFileMesh
vmat2.transparency = 0.7

cmat = "Steel"

@assembly ThreeDFiles begin
  world = Modia3D.Object3D(visualizeFrame=true)

  sphereMoving     = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidSphere(0.5) , "Steel", vmat1; contactMaterial = cmat); r=[0.0, 0.0, 0.0],  fixed=false, v_start=[2.0,0.0,-3.0]) # R=ModiaMath.rot2(-pi/2),

  box     = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidBox(20.0,4.0,1.0) , "Steel", vmat1; contactMaterial = cmat); r=[0.0, 0.0, -3.0], fixed=true ) #, R=ModiaMath.rot2(-pi/3))
end

gravField = Modia3D.UniformGravityField(g=9.81, n=[0,0,-1])
threeD = ThreeDFiles(sceneOptions=Modia3D.SceneOptions(gravityField=gravField,visualizeFrames=true, defaultFrameLength=0.7,nz_max = 100, enableContactDetection=true, visualizeContactPoints=true, visualizeSupportPoints=true))

model = Modia3D.SimulationModel( threeD )
ModiaMath.print_ModelVariables(model)
result = ModiaMath.simulate!(model; stopTime=4.0, tolerance=1e-8,interval=0.001, log=false)
ModiaMath.plot(result, ["sphereMoving.r","sphereMoving.v"])

println("... success of contactForceLaw_ballWithBox_45Deg.jl!")

end


text:
"""
    Reluplex(optimizer, eager::Bool)

Reluplex uses binary tree search to find an activation pattern that maps a feasible input to an infeasible output.

# Problem requirement
1. Network: any depth, ReLU activation
2. Input: hyperrectangle
3. Output: PolytopeComplement

# Return
`CounterExampleResult`

# Method
Binary search of activations (0/1) and pruning by optimization.

# Property
Sound and complete.

# Reference
[G. Katz, C. Barrett, D. L. Dill, K. Julian, and M. J. Kochenderfer,
"Reluplex: An Efficient SMT Solver for Verifying Deep Neural Networks," in
*International Conference on Computer Aided Verification*, 2017.](https://arxiv.org/abs/1702.01135)
"""
@with_kw struct Reluplex <: Solver
    optimizer = GLPK.Optimizer
end

function solve(solver::Reluplex, problem::Problem)
    model = Model(solver)
    encode(solver, model, problem)
    return reluplex_step(solver, model)
end

function find_relu_to_fix(ẑ, z)
    for i in 1:length(z), j in 1:length(z[i])
        ẑᵢⱼ = value(ẑ[i][j])
        zᵢⱼ = value(z[i][j])

        if type_one_broken(ẑᵢⱼ, zᵢⱼ) ||
           type_two_broken(ẑᵢⱼ, zᵢⱼ)
            return (i, j)
        end
    end
    return (0, 0)
end

 # TODO consider renaming to `inactive_broken` and `active_broken`
type_one_broken(ẑᵢⱼ, zᵢⱼ) = (zᵢⱼ > 0.0 + TOL[])  && (!(-TOL[] < ẑᵢⱼ - zᵢⱼ < TOL[])) # (zᵢⱼ > 0) && (ẑᵢⱼ != zᵢⱼ)
type_two_broken(ẑᵢⱼ, zᵢⱼ) = (-TOL[] < zᵢⱼ < TOL[]) && (ẑᵢⱼ > 0.0 + TOL[]) # (zᵢⱼ == 0) && (ẑᵢⱼ > 0)

# Corresponds to a ReLU that shouldn't be active but is
function type_one_repair!(model, ẑᵢⱼ, zᵢⱼ)
    con_one = @constraint(model, ẑᵢⱼ == zᵢⱼ)
    con_two = @constraint(model, ẑᵢⱼ >= 0.0)
    return con_one, con_two
end
# Corresponds to a ReLU that should be active but isn't
function type_two_repair!(model, ẑᵢⱼ, zᵢⱼ)
    con_one = @constraint(model, ẑᵢⱼ <= 0.0)
    con_two = @constraint(model, zᵢⱼ == 0.0)
    return con_one, con_two
end

function activation_constraint!(model, ẑᵢ, zᵢ, act::ReLU)
    # ReLU ensures that the variable after activation is always
    # greater than before activation and also ≥0
    @constraint(model, zᵢ .>= ẑᵢ)
    @constraint(model, zᵢ .>= 0.0)
end

function activation_constraint!(model, ẑᵢ, zᵢ, act::Id)
    @constraint(model, zᵢ .== ẑᵢ)
end

function encode(solver::Reluplex, model::Model,  problem::Problem)
    layers = problem.network.layers
    ẑ = init_vars(model, layers, :ẑ, with_input=true) # before activation
    z = init_vars(model, layers, :z, with_input=true) # after activation

    # Each layer has an input set constraint associated with it based on the bounds.
    # Additionally, consective variables zᵢ, ẑᵢ₊₁ are related by a constraint given
    # by the affine map encoded in the layer Lᵢ.
    # Finally, the before-activation-variables and after-activation-variables are
    # related by the activation function. Since the input layer has no activation,
    # its variables are related implicitly by identity.
    activation_constraint!(model, ẑ[1], z[1], Id())
    bounds = get_bounds(problem)
    for (i, L) in enumerate(layers)
        @constraint(model, affine_map(L, z[i]) .== ẑ[i+1])
        add_set_constraint!(model, bounds[i], z[i])
        activation_constraint!(model, ẑ[i+1], z[i+1], L.activation)
    end
    # Add the bounds on your output layer
    add_set_constraint!(model, last(bounds), last(z))
    # Add the complementary set defind as part of the problem
    add_complementary_set_constraint!(model, problem.output, last(z))
    feasibility_problem!(model)

    return nothing
end

function reluplex_step(solver::Reluplex, model::Model)

    optimize!(model)

    ẑ, z = model[:ẑ], model[:z]
    # If the problem is optimally solved, this could potentially be a counterexample.
    # Branch by repair type (inactive or active) and redetermine if this is a valid
    # counterexample. If the problem is infeasible or unbounded. The property holds.
    if termination_status(model) == OPTIMAL
        i, j = find_relu_to_fix(ẑ, z)

        # In case no broken relus could be found, return the input as a counterexample
        i == 0 && return CounterExampleResult(:violated, value(first(ẑ)))

        for repair! in (type_one_repair!, type_two_repair!)
            # Add the constraints associated with the ReLU being fixed
            new_constraints = repair!(model, ẑ[i][j], z[i][j])

            # Recurse with the ReLU i, j fixed to active or inactive
            result = reluplex_step(solver, model)

            # Return (all the way to top level) if violated.
            result.status == :violated && return result

            # Reset the model when we're done with this ReLU
            delete.(model, new_constraints)
        end
    end
    return CounterExampleResult(:holds)
end


text: myfun(x) = x^2 + x

max_stars_repo_name: Sarah-W/athing

max_issues_repo_name: Sarah-W/athing

max_forks_repo_name: Sarah-W/athing


text:
# Alphabet
# ========
#
# Alphabet of biological symbols.
#
# Subtypes of `Alphabet` represent a domain of biological characters. For
# example, `DNAAlphabet{2}` has a domain of unambiguous nucleotides (i.e. A, C,
# G, and T). These types are used for parameterizing biological sequences and so
# on. A pair of encoder and decoder is associated with an alphabet, which maps
# values between binary and Julia-level representation.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

"""
Alphabet of biological characters.
"""
abstract Alphabet

"""
DNA nucleotide alphabet.
"""
immutable DNAAlphabet{n} <: Alphabet end

"""
RNA nucleotide alphabet.
"""
immutable RNAAlphabet{n} <: Alphabet end

"""
Amino acid alphabet.
"""
immutable AminoAcidAlphabet <: Alphabet end

"""
General character alphabet.
"""
immutable CharAlphabet <: Alphabet end

"""
Void alphabet (internal use only).
"""
immutable VoidAlphabet <: Alphabet end

"""
The number of bits to represent the alphabet.
"""
function bitsof end

for n in (2, 4)
    @eval begin
        bitsof(::Type{DNAAlphabet{$n}}) = $n
        bitsof(::Type{RNAAlphabet{$n}}) = $n
    end
end
bitsof(::Type{AminoAcidAlphabet}) = 8
bitsof(::Type{CharAlphabet}) = 32
bitsof(::Type{VoidAlphabet}) = 0

Base.eltype(::Type{DNAAlphabet}) = DNANucleotide
Base.eltype(::Type{RNAAlphabet}) = RNANucleotide
Base.eltype{n}(::Type{DNAAlphabet{n}}) = DNANucleotide
Base.eltype{n}(::Type{RNAAlphabet{n}}) = RNANucleotide
Base.eltype(::Type{AminoAcidAlphabet}) = AminoAcid
Base.eltype(::Type{CharAlphabet}) = Char
Base.eltype(::Type{VoidAlphabet}) = Void

alphabet(::Type{DNAAlphabet{2}}) = DNA_A:DNA_T
alphabet(::Type{RNAAlphabet{2}}) = RNA_A:RNA_U
alphabet(::Type{DNAAlphabet{4}}) = alphabet(DNANucleotide)
alphabet(::Type{RNAAlphabet{4}}) = alphabet(RNANucleotide)
alphabet(::Type{AminoAcidAlphabet}) = alphabet(AminoAcid)
# TODO: this alphabet includes invalid Unicode scalar values
alphabet(::Type{CharAlphabet}) = typemin(Char):typemax(Char)
alphabet(::Type{VoidAlphabet}) = nothing

# Encoders & Decoders
# -------------------

"""
Encode biological characters to binary representation.
"""
function encode end

immutable EncodeError{A<:Alphabet,T} <: Exception
    val::T
end

EncodeError{A,T}(::Type{A}, val::T) = EncodeError{A,T}(val)

function Base.showerror{A}(io::IO, err::EncodeError{A})
    print(io, "cannot encode ", err.val, " in ", A)
end

"""
Decode biological characters from binary representation.
"""
function decode end

immutable DecodeError{A<:Alphabet,T} <: Exception
    val::T
end

DecodeError{A,T}(::Type{A}, val::T) = DecodeError{A,T}(val)

function Base.showerror{A}(io::IO, err::DecodeError{A})
    print(io, "cannot decode ", err.val, " in ", A)
end

for (A, T, U, ub) in [
        (DNAAlphabet{2},    DNANucleotide, UInt8,  DNA_T     ),
        (DNAAlphabet{4},    DNANucleotide, UInt8,  DNA_Gap   ),
        (RNAAlphabet{2},    RNANucleotide, UInt8,  RNA_U     ),
        (RNAAlphabet{4},    RNANucleotide, UInt8,  RNA_Gap   ),
        (AminoAcidAlphabet, AminoAcid,     UInt8,  AA_Gap    ),
        (CharAlphabet,      Char,          UInt32, '\U10ffff')]
    @assert sizeof(T) == sizeof(U)
    @assert isa(ub, T)
    @eval begin
        @inline function encode(::Type{$A}, x::$T)
            if x > $ub
                throw(EncodeError($A, x))
            end
            return reinterpret($U, x)
        end
        @inline function decode(::Type{$A}, x::$U)
            if x > $(reinterpret(U, ub))
                throw(DecodeError($A, x))
            end
            return reinterpret($T, x)
        end
        decode(::Type{$A}, x::Unsigned) = decode($A, $U(x))
    end
end


text:
# derived_domain.jl
# Generic code for a domain that is implemented in terms of an underlying domain

"""
DerivedDomain is an abstract supertype for domains that are implemented in terms
of another domain using composition. DerivedDomain transfers the interface of
Domain to this superdomain.

A concrete subtype that inherits from DerivedDomain and stores a `superdomain`
is functionally equivalent to that superdomain. Any properties of the superdomain
can be modified by overriding a suitable function. For example, `in` of the
concrete domain may be implemented in terms of the `in` of the superdomain.
"""
abstract type DerivedDomain{T} <: Domain{T}
end

# We assume that the underlying domain is stored in a field called superdomain
superdomain(d::DerivedDomain) = d.superdomain

"Return the eltype of the superdomain."
supereltype(d::DerivedDomain) = eltype(superdomain(d))

indomain(x, d::DerivedDomain) = in(x, superdomain(d))


text:
const SNP_INFO_KEYS = [:chromosome, :snpid, :genetic_distance, :position, :allele1, :allele2]
const PERSON_INFO_KEYS = [:fid, :iid, :father, :mother, :sex, :phenotype]

"""
    SnpData
    SnpData(plknm; famnm::AbstractString=plinknm*".fam", bimnm::AbstractString=plinknm*".bim" )

Type to store SNP and person information along with the SnpArray.
"""
struct SnpData
    people::Int
    snps::Int
    snparray::SnpArray
    snp_info::DataFrame
    person_info::DataFrame
    srcbed::AbstractString
    srcbim::AbstractString
    srcfam::AbstractString
end

function SnpData(
    plink_file::AbstractString, args...;
    famnm::AbstractString=plink_file*".fam",
    bimnm::AbstractString=plink_file*".bim",
    kwargs...)

    # load snp info
    snp_info = DataFrame(CSV.File(bimnm,  delim='\t', header=SNP_INFO_KEYS,
        types=[String, String, Float64, Int, String, String]); copycols = false)

    # load person info
    dlmobj = readdlm(famnm, AbstractString)
    person_info = DataFrame(table(dlmobj, header = Symbol.(:x, axes(dlmobj, 2))))
    # person_info = convert(DataFrame, readdlm(famnm, AbstractString))
    rename!(person_info, collect(f => t for (f, t) = zip(names(person_info),
                    PERSON_INFO_KEYS)))

    # load snp array
    snparray = SnpArray(string(plink_file, ".bed"), args...; famnm=famnm, kwargs...)
    people, snps = size(snparray)

    SnpData(people, snps, snparray, snp_info, person_info, plink_file*".bed", bimnm, famnm)
end

Base.size(x::SnpData) = Base.size(x.snparray)

"""
    show(io::IO, x::SnpData)

Simple string representation of SnpData
"""
function Base.show(io::IO, x::SnpData)
    print(io, "SnpData(people: $(x.people), snps: $(x.snps),\n" *
        join([
                "snp_info: \n$((join(split(string(first(x.snp_info, 6)), "\n")[2:end], "\n")))",
                x.snps > 6 ? "…," : ",",
                "person_info: \n$((join(split(string(first(x.person_info, 6)), "\n")[2:end], "\n")))",
                x.people > 6 ? "…," : ",",
                "srcbed: $(x.srcbed)\nsrcbim: $(x.srcbim)\nsrcfam: $(x.srcfam)"],
                "\n") *
        "\n)")
end

"""
    filter(s::SnpData, rowinds::AbstractVector{<:Integer}, colinds::AbstractVector{<:Integer}; des::AbstractString)
    filter(s; des::AbstractString, f_person::Function, f_snp::Function)

Filter `s` according to `f_person` and `f_snp`. The resultiing plink files are saved at `des.[bed|bim|fam]`.
"""
@inline _trueftn(x) = true
function filter(s::SnpData, rowinds::AbstractVector{<:Integer}, colinds::AbstractVector{<:Integer}; des::AbstractString=split(s.srcbed, ".bed")[1] * ".filtered")
    if eltype(rowinds) == Bool
        rmask = rowinds
    else
        rmask = falses(size(s.snparray, 1))
        rmask[rowinds] .= true
    end
    if eltype(colinds) == Bool
        cmask = colinds
    else
        cmask = falses(size(s.snparray, 2))
        cmask[colinds] .= true
    end

    snps = count(cmask)
    people = count(rmask)
    snparray = SnpArrays.filter(s.srcbed, s.srcbim, s.srcfam, rmask, cmask; des=des)
    snp_info = s.snp_info[cmask, :]
    person_info = s.person_info[rmask, :]
    SnpData(people, snps, snparray, snp_info, person_info,
            des * ".bed", des * ".bim", des * ".fam")
end

function filter(s::SnpData; des::AbstractString = split(s.srcbed, ".bed")[1] * ".filtered", f_person::Function = _trueftn, f_snp::Function = _trueftn)
    f_person == _trueftn && f_snp == _trueftn && @warn "No nontrivial function provided. Just copying."
    colinds = collect(f_snp(r)::Bool for r in eachrow(s.snp_info))
    rowinds = collect(f_person(r)::Bool for r in eachrow(s.person_info))
    SnpArrays.filter(s::SnpData, rowinds::Vector{Bool}, colinds::Vector{Bool}; des = des::AbstractString)
end

filter(s::AbstractString; des::AbstractString = split(s.srcbed, ".bed")[1] * ".filtered",
    f_person::Function = _trueftn, f_snp::Function = _trueftn) = SnpArrays.filter(SnpData(s); des = des, f_person = f_person, f_snp = f_snp)

"""
    split_plink(s, key; prefix)

Split data `s` according to chromosome, sex, or phenotype. Returns a dictionary of splitted data.
"""
function split_plink(s::SnpData, key::Symbol = :chromosome; prefix = split(s.srcbed, ".bed")[1] * string(key))
    key in [:chromosome, :sex, :phenotype] || throw(ArgumentError("key should be one of :chromosome, :sex, or :phenotype"))
    r = Dict{AbstractString, SnpData}()
    if key == :chromosome
        for chr in unique(s.snp_info[!, key])
            r[chr] = SnpArrays.filter(s; des = prefix * chr, f_snp = x -> (x[key] == chr))
        end
        r
    else
        for val in unique(s.person_info[!, key])
            r[val] = SnpArrays.filter(s; des = prefix * string(val), f_person = x -> (x[key] == val))
        end
        r
    end
end
split_plink(src::AbstractString, key::Symbol = :chromosome; prefix = src * string(key)) = split_plink(SnpData(src), key; prefix = prefix)

"""
    vcat(A...;des="tmp_vcat_" * string(vcat_counter))

Concatenate SnpData along dimension 1.
"""
vcat_counter = 1
function Base.vcat(A::SnpData...; des="tmp_vcat_" * string(vcat_counter))
    global vcat_counter
    if des == "tmp_vcat_" * string(vcat_counter)
        vcat_counter = vcat_counter + 1
    end

    # vcat person_info
    person_info = vcat([x.person_info for x in A]...)

    # get snp_info
    @assert allequal([x.snp_info for x in A]) "snp_info are different"
    snp_info = A[1].snp_info

    # vcat snparray
    snparray = vcat([x.snparray for x in A]...; des=des)

    people, snps = size(person_info,1), size(snp_info, 1)

    bimfile = des * ".bim"
    famfile = des * ".fam"

    writedlm(bimfile, hcat([snp_info[!, k] for k in SNP_INFO_KEYS]...))
    writedlm(famfile, hcat([person_info[!, k]
                                for k in PERSON_INFO_KEYS]...))

    SnpData(people, snps, snparray, snp_info, person_info, des * ".bed", bimfile, famfile)
end

"""
    hcat(A...;des="tmp_hcat_" * string(hcat_counter))

Concatenate SnpData along dimension 2.
"""
hcat_counter = 1
function Base.hcat(A::SnpData...; des="tmp_hcat_" * string(hcat_counter))
    global hcat_counter
    if des == "tmp_hcat_" * string(hcat_counter)
        hcat_counter = hcat_counter + 1
    end

    # get person_info
    @assert allequal([x.person_info for x in A]) "person_info are different"
    person_info = A[1].person_info

    # vcat snp_info
    snp_info = vcat([x.snp_info for x in A]...)

    # hcat snparray
    snparray = hcat([x.snparray for x in A]...; des=des)

    people, snps = size(person_info, 1), size(snp_info, 1)

    bimfile = des * ".bim"
    famfile = des * ".fam"

    writedlm(bimfile, hcat([snp_info[!, k] for k in SNP_INFO_KEYS]...))
    writedlm(famfile, hcat([person_info[!, k]
                                for k in PERSON_INFO_KEYS]...))

    SnpData(people, snps, snparray, snp_info, person_info, des * ".bed", bimfile, famfile)
end

"""
   hvcat(rows::Tuple{Vararg{Int}}, values...; des="tmp_hvcat_" * string(hvcat_counter))

Horizontal and vertical concatenation in one call.
"""
hvcat_counter = 1
function Base.hvcat(rows::Tuple{Vararg{Int}}, A::SnpData...; des="tmp_hvcat" * string(hvcat_counter))
    global hvcat_counter
    if des == "tmp_hvcat_" * string(hvcat_counter)
        hvcat_counter = hvcat_counter + 1
    end

    num_block_rows = length(rows)

    # collect person_info
    a = 1

    person_info = A[1].person_info

    for i = 2:num_block_rows
        person_info = vcat(person_info, A[a].person_info)
        a += rows[i]
    end
    # collect snp_info
    snp_info = vcat([x.snp_info for x in A[1:rows[1]]]...)

    # hvcat snparray
    snparray = hvcat(rows, [x.snparray for x in A]...; des=des)

    people, snps = size(person_info, 1), size(snp_info, 1)

    bimfile = des * ".bim"
    famfile = des * ".fam"

    writedlm(bimfile, hcat([snp_info[!, k] for k in SNP_INFO_KEYS]...))
    writedlm(famfile, hcat([person_info[!, k]
                                for k in PERSON_INFO_KEYS]...))

    SnpData(people, snps, snparray, snp_info, person_info, des * ".bed", bimfile, famfile)
end

"""
    merge_plink(d)
    merge_plink(des, d)

Merge the SnpData of the splitted plink files.
Returns: merged SnpData. If `des` is given, it is written on that destination.
"""
@inline function isless_chromosome(x::AbstractString, y::AbstractString)
    x_int = try parse(Int32, x)
    catch
        typemax(Int32)
    end
    y_int = try parse(Int32, y)
    catch
        typemax(Int32)
    end
    x_int == y_int ? x < y : x_int < y_int
end

function merge_plink(d::Dict{AbstractString, SnpData})
    ks = sort(collect(keys(d)), lt=isless_chromosome)

    # get person_info
    person_info = d[ks[1]].person_info

    # vcat snp_info
    @time snp_info = Base.vcat([d[k].snp_info for k in ks]...)
    # hcat snparray
    @time data = Base.hcat([d[k].snparray.data for k in ks]...)
    rowcounts = d[ks[1]].snparray.rowcounts
    @time columncounts = Base.hcat([d[k].snparray.columncounts for k in ks]...)
    snparray = SnpArray(data, rowcounts, columncounts, size(data)[1])

    people, snps = size(person_info,1), size(snp_info, 1)
    SnpData(people, snps, snparray, snp_info, person_info, "", "", "")
end

merge_plink(des::AbstractString, d::Dict{AbstractString, SnpData}) = write_plink(des, merge_plink(d))

"""
    merge_plink(prefix; des::AbstractString = prefix * ".merged")

merge the plink files beginning with `prefix`.
"""
function merge_plink(prefix::AbstractString; des::AbstractString = prefix * ".merged")
    l = glob(prefix * "*.bed")
    matching_srcs = map(x -> splitext(x)[1], l)
    d = Dict{AbstractString, SnpData}()
    @time for fn in matching_srcs
        chrsnpdata = SnpData(fn)
        chr = chrsnpdata.snp_info[!, :chromosome][1]
        @assert all(chr .== chrsnpdata.snp_info[!, :chromosome]) "Not all chrs are the same in $fn.bim."
        d[String(chr)] = chrsnpdata
    end
    merge_plink(des, d)
end

"""
    write_plink(filename, snpdata)

Write SnpData to Plink bed, bim, fam files.
Returns: `snpdata`
"""
function write_plink(filename::AbstractString, snpdata::SnpData)
    bimfile = filename * ".bim"
    bedfile = filename * ".bed"
    famfile = filename * ".fam"
    # isfile(bimfile) && error("($bimfile) already exists.")
    # isfile(bedfile) && error("($bedfile) already exists.")
    # isfile(famfile) && error("($famfile) already exists.")

    # write bim file
    snp_info = snpdata.snp_info
    writedlm(bimfile, hcat([snp_info[!, k] for k in SNP_INFO_KEYS]...))

    # write fam file
    person_info = snpdata.person_info
    writedlm(famfile, hcat([person_info[!, k]
                                for k in PERSON_INFO_KEYS]...))

    # write bed file
    open(bedfile, "w") do io
        write(io, 0x1b6c)
        write(io, 0x01)
        write(io, snpdata.snparray.data)
    end
    SnpData(snpdata.people, snpdata.snps, snpdata.snparray, snpdata.snp_info, snpdata.person_info, bimfile, bedfile, famfile)
end


text:
export LgbmModel, LgbmRegressor, LgbmClassifier

@with_kw mutable struct LgbmModel <: BaseEstimator
    pyo::PyObject = PyNULL()
    so::Vector{UInt8} = UInt8[]
    max_bin::Int = 16
    objective::String = "binary:logistic"
    warm_start::Bool = false
    num_leaves::Int = 31
    min_child_samples::Int = 1000
    n_estimators::Int = 50
    reg_lambda::Float32 = 0.1
    subsample::Float32 = 1
    colsample_bytree::Float32 = 1
end

is_classifier(m::LgbmModel) = occursin(r"binary|multi", m.objective)

function gridparams(m::LgbmModel)
    grid = [
        "num_leaves" => [15, 31],
        "min_child_samples" => [1000, 100, 5000],
        "n_estimators" => [20, 50],
        "reg_lambda" => [0.1, 1, 10],
        "subsample" => [1, 0.8],
        "colsample_bytree" => [1, 0.8],
    ]
    params = gridparams(grid)
end

function fit!(m::LgbmModel, x, y, w = nothing; columns = nothing)
    columns = something(columns, string.(1:size(x, 1)))
    @unpack max_bin, objective, warm_start, num_leaves, min_child_samples = m
    @unpack n_estimators, reg_lambda, subsample, colsample_bytree = m
    dump_lgbm_data(x, y, w, max_bin = max_bin, columns = columns)
    if occursin(r"binary|multi", objective)
        num_class = length(unique(y))
    end
    if occursin("binary", objective) && num_class > 2
        objective = "multi:softmax"
    end
    min_child_samples = min(length(y) ÷ 10, min_child_samples)
    # get node lists
    nodes = unique(pmap(n -> gethostname(), workers()))
    nnode, hosts = length(nodes), join(nodes, ',')
    # make cmd args
    args = String[]
    push!(args, "task=train")
    objective = denormalize_objective(objective, "lightgbm")
    push!(args, "objective=$objective")
    occursin("multiclass", objective) &&
    push!(args, "num_class=$num_class")
    push!(args, "tree_learner=data")
    push!(args, "n_estimators=$n_estimators")
    push!(args, "n_jobs=20")
    push!(args, "train_data=train.lgbm")
    push!(args, "num_machines=$nnode")
    push!(args, "output_model=lgbm")
    warm_start && isfile("lgbm") &&
    push!(args, "input_model=lgbm")
    push!(args, "num_leaves=$num_leaves")
    push!(args, "min_child_samples=$min_child_samples")
    push!(args, "reg_lambda=$reg_lambda")
    push!(args, "subsample=$subsample")
    push!(args, "colsample_bytree=$colsample_bytree")
    push!(args, "max_bin=$(max_bin - 1)")
    if isnothing(Sys.which("mpiexec"))
        run(`$LIGHTGBM $args`)
    else
        run(`mpiexec --host $hosts $LIGHTGBM $args`)
    end
    @from lightgbm imports Booster
    @from treelite imports Model
    pyo = Booster(model_file = "lgbm")
    em = eachmatch(r"Tree=(\d+)", read("lgbm", String))
    ntree = collect(em)[end].captures[1]
    if parse(Int, ntree) > 0
        model = Model.load("lgbm", model_format = "lightgbm")
        rm("lgbm.so", force = true)
        model.export_lib(toolchain = "gcc", libpath = "lgbm.so")
        so = read("lgbm.so")
    else
        so = UInt8[]
    end
    @pack! m = pyo, so
    BSON.bson("model.bson", model = m)
    visualize(m, columns)
    return m
end

function predict_lgbm(m::LgbmModel, x; treelite)
    @unpack pyo, so = m
    if treelite && !isempty(so)
        write("lgbm.so", so)
        @imports treelite.runtime as rt
        batch = rt.Batch.from_npy2d(pymat(x))
        pyo = rt.Predictor("lgbm.so")
        pyo.predict(batch)
    else
        pyo.predict(pymat(x))
    end
end

function predict_proba(m::LgbmModel, x; treelite = false)
    ŷ = predict_lgbm(m, x, treelite = treelite)
    ndims(ŷ) == 1 ? hcat(1 .- ŷ,  ŷ) : ŷ
end

function predict(m::LgbmModel, x; treelite = false)
    @unpack objective = m
    ŷ = predict_lgbm(m, x, treelite = treelite)
    !is_classifier(m) && return ŷ
    if ndims(ŷ) == 1
        signone.(ŷ .- 0.5)
    else
        [I[2] - 1 for I in argmax(ŷ, dims = 2)]
    end
end

function visualize(m::LgbmModel, columns)
    @unpack pyo = m
    gbm = GbmModel(pyo = pyo)
    visualize(gbm, columns)
end

reset!(m::LgbmModel) = foreach(c -> rm(c, force = true), ["train.lgbm", "lgbm", "lgbm.so", "ref.bson"])

istree(m::LgbmModel) = true

modelhash(m::LgbmModel) = hash(m.pyo)

function dump_lgbm_data(x, y, w = nothing; max_bin = 16, columns)
    w = isnothing(w) ? fill!(similar(y), 1) : w
    if isfile("ref.bson") && isfile("train.lgbm")
        BSON.@load "ref.bson" xseg yoff woff
    else
        xseg, yoff, woff = zeros(Float32), -1, -1
    end
    if xseg == x[1:min(end, 100000)] && length(unique(y)) <= 2
        open("train.lgbm", "r+") do fid
            @assert yoff > 0
            seek(fid, yoff - 1)
            write(fid, y)
            if woff > 0
                seek(fid, woff - 1)
                write(fid, w)
            end
        end
    else
        xseg = x[1:min(end, 100000)]
        x, y, w = pymat(x), vec(y), vec(w)
        @from lightgbm imports Dataset
        @from unidecode imports unidecode
        dset = Dataset(data = x, label = y, weight = w, params = Dict("max_bin" => max_bin - 1))
        columns′ = [filter(!isspace, unidecode(c)) for c in columns]
        dset.set_feature_name(columns′)
        isfile("train.lgbm") && rm("train.lgbm")
        @time dset.save_binary("train.lgbm")
        @gc dset
        data = Mmap.mmap("train.lgbm")
        yoff = align(data, reinterpret(UInt8, y))
        woff = align(data, reinterpret(UInt8, w))
        BSON.@save "ref.bson" xseg yoff woff
    end
end

LgbmRegressor(;ka...) = LgbmModel(;objective = "reg:mse", ka...)

LgbmClassifier(;ka...) = LgbmModel(;objective = "binary:logistic", ka...)


text:
@ms include("instance.jl")
@ms include("solvers/greedy.jl")

function main_greedy(args)
    println("="^70)
    println("Début de l'action greedy")

    filename = args[:infile]
    instance = Instance(filename)
    solution = Solution(instance)

    ms_start = ms()
    # initial_sort!(solution)
    greedy!(solution)
    ms_stop = ms()
    nb_sec = round(ms_stop - ms_start, digits=3)
    print_sol(solution)
    println("Fin de l'algorithme glouton")
end

main_greedy(Args.args)


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************

"""
    AbsorptionChiller

struct with inner constructor:
```julia
function AbsorptionChiller(;
    min_ton::Real = 0.0,
    max_ton::Real = 0.0,
    chiller_cop::Real,
    chiller_elec_cop::Real = 14.1,
    installed_cost_per_ton::Real,
    om_cost_per_ton::Real,
    macrs_option_years::Real = 0,
    macrs_bonus_pct::Real = 0
)
```
"""
struct AbsorptionChiller <: AbstractThermalTech
    min_ton::Real
    max_ton::Real
    chiller_cop::Real
    chiller_elec_cop::Real
    installed_cost_us_dollars_per_ton::Real
    om_cost_us_dollars_per_ton::Real
    macrs_option_years::Real
    macrs_bonus_pct::Real
    min_kw::Real
    max_kw::Real
    installed_cost_per_kw::Real
    om_cost_per_kw::Real
    function AbsorptionChiller(;
        min_ton::Real = 0.0,
        max_ton::Real = 0.0,
        chiller_cop::Real,
        chiller_elec_cop::Real = 14.1,
        installed_cost_per_ton::Real,
        om_cost_per_ton::Real,
        macrs_option_years::Real = 0,
        macrs_bonus_pct::Real = 0,
        )

        min_kw = min_ton * KWH_THERMAL_PER_TONHOUR
        max_kw = max_ton * KWH_THERMAL_PER_TONHOUR
        installed_cost_per_kw = installed_cost_per_ton / KWH_THERMAL_PER_TONHOUR
        om_cost_per_kw = om_cost_per_ton / KWH_THERMAL_PER_TONHOUR

        new(
            min_ton,
            max_ton,
            chiller_cop,
            chiller_elec_cop,
            installed_cost_per_ton,
            om_cost_per_ton,
            macrs_option_years,
            macrs_bonus_pct,
            min_kw,
            max_kw,
            installed_cost_per_kw,
            om_cost_per_kw
        )
    end
end

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
# Implementation of the AbstractMOP interface
# Formerly, this was the main problem type used internally.
# But it is not strongly typed and mutable.
# For performance, we convert it to a `StaticMOP` when
# the user is done.
# (see file `AbstractMOPInterface.jl`)

@with_kw mutable struct MixedMOP <: AbstractMOP{true}
    vars :: Vector{Int} = [];

    vector_of_objectives :: Vector{ <:AbstractObjective } = AbstractObjective[]
    objf_output_mapping :: Dict{Any,Vector{Int}} = Dict{Any,Vector{Int}}()

    lb :: Dict{ Int, Real } = Dict{Int,Real}();
    ub :: Dict{ Int, Real } = Dict{Int,Real}();

    # change whenever vars or constraints are added
    var_state :: UUIDs.UUID = UUIDs.uuid4();      # TODO maybe use automatic hashing? pro: no need to implement everything, cons: slower?
    objf_state :: UUIDs.UUID = UUIDs.uuid4();
end

# legacy constructors
function MixedMOP( lb :: Vec, ub :: Vec )
    @assert length(lb) == length(ub);
    @assert all( lb .<= ub );
    @assert length(lb) > 0;     # for unconstrained problems, we need the number of variables

    LB = Dict{Int,Real}();
    for (var_int, b) ∈ enumerate(lb)
        if !isinf( b )
            LB[var_int] = b
        end
    end

    UB = Dict{Int,Real}();
    for (var_int, b) ∈ enumerate(ub)
        if !isinf( b )
            UB[var_int] = b
        end
    end

    MixedMOP(;
        vars = collect( 1 : length(lb) ),
        lb = LB,
        ub = UB,
    )
end

function MixedMOP(n_vars :: Int)
    return MixedMOP(;
        vars = collect(1:n_vars)
    );
end

list_of_objectives( mop :: MixedMOP ) = mop.vector_of_objectives;

function output_indices( objf :: AbstractObjective, mop :: MixedMOP )
    return mop.objf_output_mapping[ objf ]
end

function _del!( mop :: MixedMOP, objf :: AbstractObjective )
    position = _objf_index( objf, mop );
    deleteat!( mop.vector_of_objectives, position );
    nothing
end

function _add!( mop :: MixedMOP, objf :: AbstractObjective, out_indices :: Union{Nothing, Vector{Int}} = nothing )
    if isnothing( out_indices )
        out_indices = let n_out = num_objectives( mop );
            collect( n_out + 1 : n_out + num_outputs(objf) )
        end
    end
    if length(out_indices) != num_outputs(objf)
        error("Number of objective outputs does not match length of output indices!");
    end
    push!(list_of_objectives(mop), objf);
    mop.objf_output_mapping[objf] = out_indices;
    mop.objf_state = UUIDs.uuid4();
    nothing
end

function MOI.add_variable( mop :: MixedMOP )
    new_var = length(mop.vars) + 1;
    push!( mop.vars, new_var );
    mop.state_hash = UUIDs.uuid4();
    return MOI.VariableIndex( new_var )
end

MOI.get( mop :: MixedMOP, :: MOI.NumberOfVariables ) = length(mop.vars) :: Int;
MOI.get( mop :: MixedMOP, :: MOI.ListOfVariableIndices) = MOI.VariableIndex.(mop.vars);

function MOI.get( mop :: MixedMOP, :: MOI.ListOfConstraints ) :: Vector{Tuple}
    constraints = [];
    for (i, var_int) ∈ enumerate(mop.vars)
        if haskey( mop.lb, var_int ) && haskey( mop.ub, var_int ) &&
            !isinf(mop.lb[var_int]) && !isinf(mop.ub[var_int])
            # TODO: what to do with half-open intervals?
            push!( constraints, ( MOI.SingleVariable, MOI.Interval ))
        end
    end
    return constraints;
end

function MOI.add_variables( mop :: MixedMOP, N :: Int )
    new_vars = [length(mop.vars) + i for i = 1 : N];
    push!( mop.vars, new_vars... );
    mop.var_state = UUIDs.uuid4();
    return MOI.VariableIndex.( new_vars )
end

function MOI.supports_constraint( ::MixedMOP, ::Type{MOI.SingleVariable}, ::Type{MOI.Interval} )
    true;
end

function MOI.add_constraint( mop :: MixedMOP, var ::F, bounds ::MOI.Interval ) where {F<:MOI.SingleVariable}
    var_int =  var.variable.value;
    mop.lb[ var_int ] = bounds.lower;
    mop.ub[ var_int ] = bounds.upper;
    mop.var_state = UUIDs.uuid4();
    return MOI.ContstraintIndex{F, MOI.Interval}
end

# TODO support half-open intervals

@memoize ThreadSafeDict function _full_lb( mop :: MixedMOP, hash :: UUIDs.UUID )
    [ haskey(mop.lb, var_int) ? mop.lb[var_int] : -Inf for var_int ∈ mop.vars ]
end

function full_lower_bounds( mop :: MixedMOP )
    _full_lb( mop, mop.var_state )
end

@memoize ThreadSafeDict function _full_ub( mop :: MixedMOP, hash :: UUIDs.UUID )
    [ haskey(mop.ub, var_int) ? mop.ub[var_int] : Inf for var_int ∈ mop.vars ]
end

function full_upper_bounds( mop :: MixedMOP )
    _full_ub( mop, mop.var_state )
end

# overwrite these derived methods for efficacy
@memoize ThreadSafeDict function _full_lb_internal(mop :: MixedMOP, hash :: UUIDs.UUID )
    [ isinf(l) ? l : 0.0 for l ∈ full_lower_bounds(mop) ];
end

@memoize ThreadSafeDict function _full_ub_internal(mop :: MixedMOP, hash :: UUIDs.UUID )
    [ isinf(u) ? u : 1.0 for u ∈ full_upper_bounds(mop) ];
end

full_lower_bounds_internal( mop :: MixedMOP ) = _full_lb_internal( mop , mop.var_state );
full_upper_bounds_internal( mop :: MixedMOP ) = _full_ub_internal( mop , mop.var_state );

# Evaluating and sorting objectives

# overwrite to exploit memoization
@memoize ThreadSafeDict function _output_indices( mop :: MixedMOP, hash :: UUIDs.UUID )
    all_outputs = Int[];
    for objf ∈ list_of_objectives( mop )
        push!( all_outputs, output_indices( objf, mop )...);
    end
    return all_outputs;
end

function output_indices( mop :: MixedMOP )
    return _output_indices(mop, mop.objf_state )
end

@memoize ThreadSafeDict function _reverse_internal_sorting_indices( mop :: MixedMOP, hash :: UUIDs.UUID )
    internal_indices = output_indices(mop);
    return sortperm( internal_indices );
end

function reverse_internal_sorting_indices(mop :: MixedMOP)
    return _reverse_internal_sorting_indices(mop, mop.objf_state);
end