text:
# include("FieldCalc.jl")
# # VecDist(X::Array) = [ X[2,1]-X[1,1], X[2,2]-X[1,2] , X[2,3]-X[1,3]]
# #X is the 2 coordinates to take distance between in format of [x₁,y₁,z₁;x₂,y₂,z₂]

abstract type Toroid end

struct DCore <: Toroid

    FlatHeight::Float64
    MaxHeight::Float64 #Height from CENTER LINE, so the actual core is 2 times as tall overall. This is the half-height effectively
    RadiusAtPeak::Float64 #This is how far from the center axis the peak of the geometry is, not the radius of curvature at the apex
    ID::Float64 #Thner diameter of the core
    OD::Float64 #Thner diameter of the core
    Turns::Int64
    Layers::Int64
    WireLength::Float64
    Resistance::Float64
end

struct Circle <: Toroid

    CoreRadius::Float64
    CenterRadius::Float64
    ID::Float64 #Thner diameter of the core
    OD::Float64 #Thner diameter of the core
    Turns::Int64
    Layers::Int64
    WireLength::Float64
    Resistance::Float64
end
struct GeneralParams <: Toroid
    SingleLayerInductance::Float64
    WireDiameter::Float64
    TargetInductance::Float64

end

mutable struct Geom
    Circ::DCore
    DCore::Circle
    General::GeneralParams
end

function ToroidOptimizer(
    Dia,
    LTarget;
    NumLayers = 2,
    CoreMu = 1,
    Alpha = 2,
    CuFillFactor = 1,
)
    ## Summary: To use this function enter the wire diameter and the target
    #inductance (in meters and Henries). The calculations assume "NumLayers" fully
    #wound layers as it is a good balance between heat dissipation and
    #geometric efficiency (2 is a good option, unless heat is an issue). The output is a struct containing the geometric
    #information for a toroid with a D shaped core, and a circular core.

    ## Constants
    μ_o = 4 * pi * 10^(-7)
    μ = μ_o * CoreMu
    L0 = μ * Dia / (2 * pi)
    L = LTarget / NumLayers^2 #Use the assumption that there will be 'NumLayers' fully wound layers...
    #If there are two fully wound layers, it can be calculated at just one
    #layer with 1/4 the inductance. Then double it when winding
    LperL0 = L / L0 #An intermediate dimensionless parameter that comes from the papers which these formulas are derived.
    SingleLayerInductance = L
    WireDiameter = Dia
    TargetInductance = LTarget
    General =
        GeneralParams(SingleLayerInductance, WireDiameter, TargetInductance)

    ## D shaped core
    #Alpha = 2; #User can change this to an integer between 2 and 5 (inclusive) varies the aspect ratio of the D
    #A high alpha value such as 5 will be a taller/bigger core with fewer
    #turns. It is theoretically more efficient but will probably have more
    #leakage flux.

    # AlphaMat Key. The columns represent [alpha e h p s t]

    AlphaMat = [
        2 0.26 0.645 3.6 0.72 1.77
        3 0.85 1.5 8.0 2.74 4.42
        4 1.6 2.4 12.8 5.76 8.09
        5 2.45 3.4 17.9 9.61 12.6
    ] #Table 1 From  "D-Shaped toroidal cage inductors" P.N. Murgatroyd & D. Belahrache,1989

    T = AlphaMat[Alpha-1, 6] #Picking a value in the table for clarity
    S = AlphaMat[Alpha-1, 5] #Also Selecting value for clarity
    P = AlphaMat[Alpha-1, 4] # As above

    KFunk_D(k) =
        (2 * pi)^0.5 * (S / (P^(3 / 2))) * k .^ (3 / 2) + 0.25 * k - LperL0 #From "The Optimal Form for Coreless Inductor", P. Murgatroyd IEEE TMI, 25 No. 3 1989
    K = mybisect(KFunk_D, 0, 1e5, 50) #solving for ideal K with bisection method
    WireLength_D = K * Dia #Equation 6 in "Economic designs for..."
    Turns_D = (2 * pi * K / P)^(1 / 2) # Also from P. Murgatroyd, section 4, Eq. 20 in "Economic designs for single-layer toroidal inductors"
    # TurnsD = sqrt(2*pi*L/(T*μ*B)); #Turns needed in a D shaped core toroid
    # ID_DToroid = 2*B;
    # B = WireLength_D/(Turns_D*P); #B=inner radius, so B = total length/(number of turns*Perimeter of each turn)
    B = Dia / 2 + Turns_D * Dia / (2 * pi) #Making it so that the wires all tough on the inner edge of the core
    ## For the following reference Fig 4 in "D-Shaped toroidal cage inductors" P.N. Murgatroyd & D. Belahrache,1989
    FlatHeight = B * AlphaMat[Alpha-1, 2]
    MaxHeight = B * AlphaMat[Alpha-1, 3] #Height from CENTER LINE, so the actual core is 2 times as tall overall. This is the half-height effectively
    RadiusAtPeak = Alpha^0.5 * B #This is how far from the center axis the peak of the geometry is, not the radius of curvature at the apex
    ID = B * 2#Thner diameter of the core
    OD = 2 * Alpha * B
    Turns = round(Turns_D * NumLayers)
    Layers = NumLayers
    WireLength = WireLength_D * NumLayers
    D_WireR = Length2Resist(WireLength, Dia; FillFac = CuFillFactor)
    DCoreGeom = DCore(
        FlatHeight,
        MaxHeight,
        RadiusAtPeak,
        ID,
        OD,
        round(Turns),
        Layers,
        WireLength,
        D_WireR,
    )

    ## Two-layer circular core

    KFunk(k) = 0.2722 * k .^ (3 / 2) + 0.25 * k - LperL0 #From "The Optimal Form for Coreless Inductor", P. Murgatroyd IEEE TMI, 25 No. 3 1989
    K = mybisect(KFunk, 0, 1e5, 50) #solving for ideal K with bisection method
    WireLength = K * Dia
    Turns = 0.8165 * K^(1 / 2) # Also from P. Murgatroyd, section 2 in "Economic designs for single-layer toroidal inductors"
    R = WireLength / (2 * pi * Turns)# Equation 2 in "Economic designs for single-layer toroidal inductors",P. Murgatroyd
    T = Dia / (sin(pi / Turns) * 2) + R# Equation 3 in "Economic designs for single-layer toroidal inductors"

    Turns = Turns * NumLayers
    WireLength = WireLength * NumLayers
    Layers = NumLayers
    ID = 2 * (T - R)
    OD = 2 * (T + R)
    CoreRadius = R
    CenterRadius = T
    Circ_WireR = Length2Resist(WireLength, Dia; FillFac = CuFillFactor)
    CircleGeom = Circle(
        CoreRadius,
        CenterRadius,
        ID,
        OD,
        round(Turns),
        Layers,
        WireLength,
        Circ_WireR,
    )

    ## Back-checking the work with two approximations.

    Approx_EnclosedArea_DCore =
        DCoreGeom.MaxHeight * 2 * (DCoreGeom.OD - DCoreGeom.ID) / 2 * 0.75#This ROUGH approximation is assuming the core fills 75# of the circumscribed rectangle around the core (roughly a semi-cirle)
    SanityCheckInductance =
        μ * DCoreGeom.Turns^2 * Approx_EnclosedArea_DCore /
        (pi * (DCoreGeom.OD + DCoreGeom.ID) / 2)

    MeanPathLength = (DCoreGeom.OD + DCoreGeom.ID) * pi

    CoreFlux(
        MeanPathLength,
        Approx_EnclosedArea_DCore,
        DCoreGeom.Turns;
        CoreMu = CoreMu,
        Current = 1e-3,
        AirGapSize = 0,
    )

    Laa = μ * CircleGeom.Turns^2 * R^2 / (2 * T)#Another formula for inductance of circular
    # toroid from hyperphysics
    println("Target inductance is: $(round(LTarget*1e6)) μH")
    println(
        "Sanity check inductance is: $(round(SanityCheckInductance*1e6)) μH",
    )
    println("
    FlatHeight     :$(round(DCoreGeom.FlatHeight*1000;sigdigits=2)) mm
    MaxHeight      :$(round(DCoreGeom.MaxHeight*1000;sigdigits=2)) mm
    RadiusAtPeak   :$(round(DCoreGeom.RadiusAtPeak*1000;sigdigits=2)) mm
    ID             :$(round(DCoreGeom.ID*1000;sigdigits=2)) mm
    OD             :$(round(DCoreGeom.OD*1000;sigdigits=2)) mm
    Turns          :$(DCoreGeom.Turns) Turns
    Layers         :$(DCoreGeom.Layers) Turns
    WireLength     :$(round(DCoreGeom.WireLength;sigdigits=2)) meters
    WireResist     :$(round(DCoreGeom.Resistance;sigdigits=2)) Ω")
    CoilGeom = Geom(DCoreGeom, CircleGeom, General)
    return CoilGeom
end

function CoreFlux(Le, Ae, Turns; CoreMu = 1, Current = 1, AirGapSize = 0)
    ## This assumes a simple geometry where the gap area equals Ferrite area
    μ_o = 4 * pi * 1e-7
    μ = μ_o * CoreMu
    Φ =
        μ_o * Turns * Current /
        ((Le - AirGapSize) / (CoreMu * Ae) + AirGapSize / Ae)
    L = Turns * Φ
    Flux = Φ / Ae
    println("The total B field within the core given $(Current*1000)mA
     is $(round(Flux*1000;sigdigits=2))mT")
    println("~~~~~~~~~~~~")
    return Flux

end

function Length2Resist(L, Diam; ρ = 1.68e-8, FillFac = 1)
    if Diam > 0.1
        println(
            "Diameter must be in METERS, it was automatically multiplied by 1e-3 in case you entered it in mm",
        )
        Diam = Diam * 1e-3
    end
    A = pi / 4 .* Diam .^ 2
    return L * ρ * FillFac / A
end

"""
This function takes in an inner and outer radius and makes the ideal "D-shaped" toroid core

It returns a list of coordinates (r,z) where r is the radial location
on the toroid, and z is the height from the center axis. In a cross-section
view this is essentially x and y.

r₁ - inner radius
r₂ - outer radius
dr (kwarg) is the dr between each sampled location

"""
function DCoreGeom(r₁, r₂; dr = 0.0001,NPts = nothing,PlotOn=false,UpsamplePoints = 1e4)
    dZdr(r) = log.(sqrt(r₁ * r₂) ./ r) ./ sqrt.(log.(r ./ r₁) .* log.(r₂ ./ r))
    if ~(NPts===nothing)
        dr = (r₂ - r₁) / UpsamplePoints
    end
    rᵢ = r₁ + dr / 10#starting at a small fraction above zero, to avoid an inf.
    r = [rᵢ]
    zᵢ = 0.0
    z = [zᵢ] # initializing vector
    II = 1 #iterator
    NumIters = floor((r₂ - r₁) / dr)

    while (NumIters > II)
        dz = dZdr(rᵢ)
        zᵢ +=dz*dr
        rᵢ += dr

        push!(z, zᵢ)
        push!(r, rᵢ)
        II += 1
    end
    z .-= z[end] #the boundary condition at the end must be fixed to equal 0

    z[1] = 0 #because of the inf initial dz/dr, the initial boundary must be manually set to zero
    ZFlatHeight = z[2]
    StartPts = Int(ceil(ZFlatHeight/dr))
    for i in 1:StartPts
        insert!(z,1+i,z[1+i]*i/StartPts)
        insert!(r,1+i,r[1])
    end

    z = vcat(z[:],-1 .*reverse(z[2:end-2]))
    r = vcat(r[:],reverse(r[2:end-2]))
    DownsampleIndex = 1:Int(round(UpsamplePoints/NPts)):length(z)
    z = z[DownsampleIndex]
    r = r[DownsampleIndex]
    if PlotOn
        plot(r, z)
    end
    return r, z
end

function DCore_PointPath(r₁, r₂; NPts = 100)
    r, z = DCoreGeom(r₁, r₂; NPts = NPts*10,PlotOn=false)
    CumulativeDist = 0
    for i in 1:(length(r)-1)
        CumulativeDist += sqrt((r[i+1]-r[i])^2+(z[i+1]-z[i])^2)
    end
    TargetDist = CumulativeDist/(NPts+1)
    rᵣ = [r[1]]
    zᵣ = [z[1]]
    CumulativeDist2 = 0
    for i in 1:(length(r)-1)
        CumulativeDist2 += sqrt((r[i+1]-r[i])^2+(z[i+1]-z[i])^2)
        if CumulativeDist2>TargetDist
            push!(rᵣ,r[i])
            push!(zᵣ,z[i])
            CumulativeDist2 = 0
        end
    end

    PointPathZeros = zeros(size(rᵣ))
    PointPath = hcat(rᵣ,zᵣ,PointPathZeros)
    return PointPath
end

function FieldMapDToroid(r₁, r₂; NPts = 100, TestLayers = 20)

    PointPath = DCore_PointPath(r₁, r₂;NPts = NPts)
    FieldMapPointPath(PointPath, TestLayers; WeightRadius = true, InvWeights = true)
end

"""
RectCoords=[ 1   2.5 0
         5   2.5 0
         5  -2.5 0
         1  -2.5 0]

         Dr₁ = 1
         Dr₂ = 5

         Cr₁ = 2
         Cr₂ = 2
         CircCenter = [3,0,0]
         CompareD_Circ_Rect_Toroid(Dr₁, Dr₂,RectCoords,Cr₁, Cr₂,CircCenter,20)

"""
function CompareD_Circ_Rect_Toroid(Dr₁, Dr₂,RectCoords,Cr₁, Cr₂,CircCenter,TestLayers = 20)
    figure(75)
    FieldMapDToroid(Dr₁,Dr₂)

    figure(76)
    R_PointPath = MakeRectPointPath(RectCoords;NElement = 50, PlotOn=false)
    FieldMapPointPath(R_PointPath, 20; WeightRadius = true,InvWeights = true)

     figure(77)
    C_PointPath = MakeEllip(Cr₁,Cr₂;Center = CircCenter,NPts=100)
    C_PointPath = reverse(C_PointPath,dims=1)
FieldMapPointPath(C_PointPath, 20; WeightRadius = true,InvWeights = true)

end

"""
This function takes in the parameters of a circular Rogowski coil
    N: Num turns,
    ID: ID of circle in meters
    OD: OD of circle in meters
    and kwargs:
    Current: Current in the wire being sensed
    μ: permeability of the material
    ω: frequency (rad/sec) of current
"""
function Rogowski_Calc(N,ID,OD;Current=1,μ=4*π*1e-7,ω = 25e3*2*π)

    # IStraightWire(r) = μ*Current/(2*π*r)
    # CircWidth(x,ID,OD) = √(((OD-ID)/2)^2-(x-(ID+OD)/2)^2)

    # dr =  (OD-ID)/1e4
    # MutualInductance = N*sum([CircWidth(xx,ID,OD).*(IStraightWire(xx)).*dr for xx in (ID+dr):dr:OD][:])./Current #Derivation by integrating flux over area
    M₂₁ =μ*N/2*(ID+OD-2*√(ID*OD)) #From multiple sources.

    V = ω*M₂₁*Current
    return (M₂₁,V)
end

"""
This function takes in the parameters of a rectagular Rogowski coil
    N: Num turns,
    h: height of rectangle in meters
    ID: ID of rectangle in meters
    OD: OD of rectangle in meters
    and kwargs:
    Current: Current in the wire being sensed
    μ: permeability of the material
    ω: frequency (rad/sec) of current
"""
function Rogowski_Calc(N,h,ID,OD;Current=1,μ=4*π*1e-7,ω = 25e3*2*π)

    M₂₁ =μ*N/(2*π)*h*(log(OD/ID)) #From multiple sources.

    V = ω*M₂₁*Current
    return (M₂₁,V)
end

function DCore_DetermineIdealInduct(ID,N,α)
    #Equation 7
     # AlphaMat Key. The columns represent [alpha e h p s t]

     AlphaMat = [
        2 0.26  0.645 3.6   0.72  1.77
        3 0.85  1.5   8.0   2.74  4.42
        4 1.6   2.4   12.8  5.76  8.09
        5 2.45  3.4   17.9  9.61  12.6
        6 3.4   4.5   23.3  14.2  17.8
        7 4.4   5.6   28.9  19.4  23.7
    ] #Table 1 From  "D-Shaped toroidal cage inductors" P.N. Murgatroyd & D. Belahrache,1989

    T = AlphaMat[α-1, 6] #Picking a value in the table for clarity
    S = AlphaMat[α-1, 5] #Also Selecting value for clarity
    #S₁ = AlphaMat[α, 5] #Also Selecting value for clarity
    P = AlphaMat[α-1, 4] # As above
    E = AlphaMat[α-1, 2] # As above
    B = ID/2 #For consistency with the paper
    #L = μ₀*N^2*B / (2*π)*((2*E+1)/4 + 2/3*S + S₁)
    L = μ₀*N^2*B / (2*π) * T
    return L
end

function CircCore_DetermineIdealInduct(IRad,N,α)
    IRad = IRad*100
    ORad = IRad*α
    R = (ORad+IRad)/2
    a = (ORad-IRad)/2
    return 1e-6*(0.01257*N^2*(R-√*(R^2-a^2)))
end


text:
lines = map(collect, readlines("data_x99.txt"))
data = permutedims(hcat(lines...))

movements = [(1,1),(3,1),(5,1),(7,1),(1,2)]

function countTrees(data, movement)
    c_move = movement[1]
    r_move = movement[2]
    trees = 0

    rows = size(data)[1]
    cols = size(data)[2]
    col = 1
    for row in 1:r_move:rows
        #println("row=", row, " col=", col)
        if(data[row,col] == '#')
            #println(data[row,col])
            trees += 1
        end
        col += c_move
        if(col > cols)
            break
        end
    end
    println("Slope with movement ", movement, " has ", trees, " trees")
    return trees
end

function getResult(data, movements)
    result = 1
    for movement in movements
        result *= countTrees(data,movement)
    end
    return result
end

println("result = ", getResult(data, movements))

# 1,1 = 68
# 3,1 = 203
# 5,1 = 78
# 7,1 = 77
# 1,2 = 40

3316272960

ext: jl


text:
@symbol_func function phi_RU(cur_reactor::AbstractReactor)
  cur_phi = L_P(cur_reactor)

  cur_phi *= f_IN(cur_reactor) * cur_reactor.I_P

  cur_phi *= 1e6

  cur_phi
end


text:
######### CmdStan batch program example  ###########

using CmdStan, StatsPlots

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  eightschools ="
  data {
    int<lower=0> J; // number of schools
    real y[J]; // estimated treatment effects
    real<lower=0> sigma[J]; // s.e. of effect estimates
  }
  parameters {
    real mu;
    real<lower=0> tau;
    real eta[J];
  }
  transformed parameters {
    real theta[J];
    for (j in 1:J)
      theta[j] <- mu + tau * eta[j];
  }
  model {
    eta ~ normal(0, 1);
    y ~ normal(theta, sigma);
  }
  "

  schools8data = Dict("J" => 8,
      "y" => [28,  8, -3,  7, -1,  1, 18, 12],
      "sigma" => [15, 10, 16, 11,  9, 11, 10, 18],
      "tau" => 25
    )

  global stanmodel, rc, chns, cnames
  stanmodel = Stanmodel(name="schools8", model=eightschools,
    output_format=:mcmcchains);
  rc, chn, cnames = stan(stanmodel, schools8data, ProjDir, CmdStanDir=CMDSTAN_HOME)

  if rc == 0

    chns = set_section(chn, Dict(
      :parameters => ["mu", "tau"],
      :thetas => ["theta.$i" for i in 1:8],
      :etas => ["eta.$i" for i in 1:8],
      :internals => ["lp__", "accept_stat__", "stepsize__", "treedepth__", "n_leapfrog__",
        "divergent__", "energy__"]
      )
    )

    if isdefined(Main, :StatsPlots)
      p1 = plot(chns)
      savefig(p1, "traceplot.pdf")
      #p2 = plot(chns, [:thetas])
      #savefig(p2, "thetas.pdf")
    end

    show(chns)
    println("\n")
    summarize(chns, sections=[:thetas])

  end
end # cd


text:
using InterpretableValidation
using Test

@test not_inv(:anybool, Random.GLOBAL_RNG) == :anybool
@test not_inv(true, Random.GLOBAL_RNG) == false
@test not_inv(false, Random.GLOBAL_RNG) == true

@test and_inv(true, Random.GLOBAL_RNG) == (true, true)
@test and_inv(false, Random.GLOBAL_RNG) in [(false, :anybool), (:anybool, false)]
@test and_inv(:anybool, Random.GLOBAL_RNG) == (:anybool, :anybool)

@test or_inv(true, Random.GLOBAL_RNG) in [(true, :anybool), (:anybool, true)]
@test or_inv(false, Random.GLOBAL_RNG) == (false, false)
@test or_inv(:anybool, Random.GLOBAL_RNG) == (:anybool, :anybool)

@test all_inv(true, 3, Random.GLOBAL_RNG) == ([true, true, true],)
@test all_inv(false, 2, Random.GLOBAL_RNG) in [([false, :anybool],), ([:anybool, false],)]
@test all_inv(:anybool, 3, Random.GLOBAL_RNG) == ([:anybool, :anybool, :anybool],)

@test any_inv(true, 2, Random.GLOBAL_RNG) in [([true, :anybool],), ([:anybool, true],)]
@test any_inv(false, 3, Random.GLOBAL_RNG) == ([false, false, false],)
@test any_inv(:anybool, 3, Random.GLOBAL_RNG) == ([:anybool, :anybool, :anybool],)

@test bitwise_and_inv([true, true], Random.GLOBAL_RNG) == ([true, true], [true, true])
@test bitwise_and_inv([:anybool, false], Random.GLOBAL_RNG) in [([:anybool, :anybool], [:anybool, false]),
                                    ([:anybool, false], [:anybool, :anybool])]

@test bitwise_or_inv([false, false], Random.GLOBAL_RNG) == ([false, false], [false, false])
@test bitwise_or_inv([:anybool, true], Random.GLOBAL_RNG) in [([:anybool, :anybool], [:anybool, true]),
                                        ([:anybool, true], [:anybool, :anybool])]

@test all_before([true, false, true], 2) == false
@test all_before([true, false, true], 1) == true
@test all_before_inv(true, 2, 3, Random.GLOBAL_RNG) == ([true, true, :anybool],)
@test all_before_inv(false, 2, 3, Random.GLOBAL_RNG)[1] in ([:anybool, false, :anybool], [false, :anybool, :anybool] )

@test all_after([true, false, true], 2) == false
@test all_after([true, false, true], 3) == true
@test all_after_inv(true, 2, 3, Random.GLOBAL_RNG) == ([:anybool, true, true],)
@test all_after_inv(false, 2, 3, Random.GLOBAL_RNG)[1] in ([:anybool, false, :anybool], [:anybool, :anybool, false] )

@test !all_between([false, true, true, false], 2, 4)
@test !all_between([false, true, true, false], 1, 3)
@test all_between([false, true, true, false], 2, 3)
@test all_between_inv(true, 2, 3, 4, Random.GLOBAL_RNG) == ([:anybool, true, true, :anybool],)
@test all_between_inv(false, 2, 3, 4, Random.GLOBAL_RNG)[1] in ([:anybool, false, :anybool, :anybool],[:anybool, :anybool, false, :anybool])

@test !any_between([true, false, false, false], 2, 4)
@test !any_between([false, true, true, false], 4, 4)
@test any_between([false, true, false, false], 2, 3)
@test any_between_inv(true, 2, 3, 4, Random.GLOBAL_RNG)[1] in ([:anybool, :anybool, true, :anybool], [:anybool, true, :anybool, :anybool])
@test any_between_inv(false, 2, 3, 4, Random.GLOBAL_RNG) == ([:anybool, false, false, :anybool],)

constraints = sample_constraints(Meta.parse("all(x .>= 1.) && any(x .>= 10)"), 10, Random.GLOBAL_RNG)
@test length(constraints) == 2
@test constraints[1][1] == Meta.parse("x .>= 1.")
@test constraints[1][2] == fill(true, 10)
@test constraints[2][1] == Meta.parse("x .>= 10")
@test sum(constraints[2][2] .== :anybool) == 9
@test sum(constraints[2][2] .== true) == 1

# Test nested temporal operators
expr = Meta.parse("any_between( all(x .>= 1.), 1,2)")
constraints = sample_constraints(expr, 10, Random.GLOBAL_RNG)

expr = Meta.parse("all( any_between(x .>= 1., 1, 3))")
x = collect(1:30.)
constraints = sample_constraints(expr, 30, Random.GLOBAL_RNG)
mvts = MvTimeseriesDistribution(:x => IID(x, Uniform(0,1)))
constrain_timeseries!(mvts,constraints)


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -37.14009556780016 <= q <= 48.365016418887166)

                     add_NL_constraint(m, :(swish(0.24902435571631587 + -0.02312729681609138*swish(-0.46870412880299783 + 0.1700038462852933*swish(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.44738483682757924*swish(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.5112184633122538*swish(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.17714072857386487*swish(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9529511507705823*swish(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.5418072450548097*swish(0.5019192867350393 + 0.2997113620022893*swish(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.5721484689942011*swish(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.4073624372378828*swish(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + -0.07630883358086304*swish(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.3103910161241421*swish(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.9789829058919253*swish(-0.652369829242545 + -0.998242239199818*swish(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.041041478696029365*swish(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.4348093907280841*swish(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.7773011179238121*swish(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7537599741964742*swish(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.13446361488468206*swish(-0.40141742757451926 + 0.7289092837988314*swish(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + 0.16077592157977172*swish(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.6378484913808045*swish(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9326240829024273*swish(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9638370498324353*swish(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.3069132295129271*swish(0.17680227033398044 + 0.8293506390949057*swish(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.09356004658531791*swish(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.16358082594109513*swish(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9029174038261001*swish(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7318887196209252*swish(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3])))) + swish(-0.017552134571214406 + -0.41004797399293524*swish(-0.46870412880299783 + 0.1700038462852933*swish(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.44738483682757924*swish(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.5112184633122538*swish(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.17714072857386487*swish(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9529511507705823*swish(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.47019448283071075*swish(0.5019192867350393 + 0.2997113620022893*swish(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.5721484689942011*swish(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.4073624372378828*swish(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + -0.07630883358086304*swish(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.3103910161241421*swish(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.35795867979675133*swish(-0.652369829242545 + -0.998242239199818*swish(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.041041478696029365*swish(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.4348093907280841*swish(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.7773011179238121*swish(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7537599741964742*swish(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.44715886126593407*swish(-0.40141742757451926 + 0.7289092837988314*swish(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + 0.16077592157977172*swish(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.6378484913808045*swish(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9326240829024273*swish(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9638370498324353*swish(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.22638221307132733*swish(0.17680227033398044 + 0.8293506390949057*swish(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.09356004658531791*swish(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.16358082594109513*swish(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9029174038261001*swish(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7318887196209252*swish(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3])))) + swish(-0.4029592795849051 + -0.6197606146665708*swish(-0.46870412880299783 + 0.1700038462852933*swish(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.44738483682757924*swish(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.5112184633122538*swish(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.17714072857386487*swish(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9529511507705823*swish(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.8464388391219866*swish(0.5019192867350393 + 0.2997113620022893*swish(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.5721484689942011*swish(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.4073624372378828*swish(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + -0.07630883358086304*swish(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.3103910161241421*swish(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.6262663939568953*swish(-0.652369829242545 + -0.998242239199818*swish(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.041041478696029365*swish(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.4348093907280841*swish(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.7773011179238121*swish(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7537599741964742*swish(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.06090941414894013*swish(-0.40141742757451926 + 0.7289092837988314*swish(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + 0.16077592157977172*swish(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.6378484913808045*swish(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9326240829024273*swish(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9638370498324353*swish(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.6699649674430836*swish(0.17680227033398044 + 0.8293506390949057*swish(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.09356004658531791*swish(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.16358082594109513*swish(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9029174038261001*swish(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7318887196209252*swish(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3])))) + swish(0.6579858782021573 + -0.21911632564812766*swish(-0.46870412880299783 + 0.1700038462852933*swish(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.44738483682757924*swish(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.5112184633122538*swish(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.17714072857386487*swish(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9529511507705823*swish(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.7292607357990395*swish(0.5019192867350393 + 0.2997113620022893*swish(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.5721484689942011*swish(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.4073624372378828*swish(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + -0.07630883358086304*swish(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.3103910161241421*swish(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.23882038501163416*swish(-0.652369829242545 + -0.998242239199818*swish(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.041041478696029365*swish(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.4348093907280841*swish(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.7773011179238121*swish(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7537599741964742*swish(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.6855544316424598*swish(-0.40141742757451926 + 0.7289092837988314*swish(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + 0.16077592157977172*swish(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.6378484913808045*swish(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9326240829024273*swish(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9638370498324353*swish(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.2868272308560216*swish(0.17680227033398044 + 0.8293506390949057*swish(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.09356004658531791*swish(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.16358082594109513*swish(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9029174038261001*swish(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7318887196209252*swish(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3])))) + swish(-0.10147744083820065 + -0.7175110849623545*swish(-0.46870412880299783 + 0.1700038462852933*swish(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.44738483682757924*swish(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.5112184633122538*swish(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.17714072857386487*swish(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9529511507705823*swish(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.9583278917345623*swish(0.5019192867350393 + 0.2997113620022893*swish(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.5721484689942011*swish(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.4073624372378828*swish(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + -0.07630883358086304*swish(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.3103910161241421*swish(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.7313133539296679*swish(-0.652369829242545 + -0.998242239199818*swish(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.041041478696029365*swish(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.4348093907280841*swish(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.7773011179238121*swish(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7537599741964742*swish(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.008349755673621395*swish(-0.40141742757451926 + 0.7289092837988314*swish(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + 0.16077592157977172*swish(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.6378484913808045*swish(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9326240829024273*swish(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9638370498324353*swish(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.13663535205769417*swish(0.17680227033398044 + 0.8293506390949057*swish(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.09356004658531791*swish(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.16358082594109513*swish(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9029174038261001*swish(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7318887196209252*swish(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# Copyright 2018 Erik Edin
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

abstract type ExecutionEnvironment end

struct NoExecutionEnvironment <: ExecutionEnvironment end
beforescenario(::NoExecutionEnvironment, ::StepDefinitionContext, ::Gherkin.Scenario) = nothing
afterscenario(::NoExecutionEnvironment, ::StepDefinitionContext, ::Gherkin.Scenario) = nothing
beforefeature(::NoExecutionEnvironment, ::Gherkin.Feature) = nothing
afterfeature(::NoExecutionEnvironment, ::Gherkin.Feature) = nothing
beforeall(::NoExecutionEnvironment) = nothing
afterall(::NoExecutionEnvironment) = nothing

module GlobalExecEnv
    envs = Dict{Symbol, Function}()

    function clear()
        global envs
        envs = Dict{Symbol, Function}()
    end
end

macro beforescenario(ex::Expr)
    envdefinition = :( $ex )
    quote
        GlobalExecEnv.envs[:beforescenario] = $(esc(envdefinition))
    end
end

macro afterscenario(ex::Expr)
    envdefinition = :( $ex )
    quote
        GlobalExecEnv.envs[:afterscenario] = $(esc(envdefinition))
    end
end

macro beforefeature(ex::Expr)
    envdefinition = :( $ex )
    quote
        GlobalExecEnv.envs[:beforefeature] = $(esc(envdefinition))
    end
end

macro afterfeature(ex::Expr)
    envdefinition = :( $ex )
    quote
        GlobalExecEnv.envs[:afterfeature] = $(esc(envdefinition))
    end
end

macro beforeall(ex::Expr)
    envdefinition = :( $ex )
    quote
        GlobalExecEnv.envs[:beforeall] = $(esc(envdefinition))
    end
end

macro afterall(ex::Expr)
    envdefinition = :( $ex )
    quote
        GlobalExecEnv.envs[:afterall] = $(esc(envdefinition))
    end
end

struct FromSourceExecutionEnvironment <: ExecutionEnvironment
    envdefinitions::Dict{Symbol, Function}

    function FromSourceExecutionEnvironment(source::String)
        include_string(Main, source)

        this = new(GlobalExecEnv.envs)
        GlobalExecEnv.clear()
        this
    end
end

function invokeenvironmentmethod(
        executionenv::FromSourceExecutionEnvironment,
        methodsym::Symbol,
        args...)

    if haskey(executionenv.envdefinitions, methodsym)
        method = executionenv.envdefinitions[methodsym]
        Base.invokelatest(method, args...)
    end
end

beforescenario(executionenv::FromSourceExecutionEnvironment,
               context::StepDefinitionContext,
               scenario::Gherkin.Scenario) = invokeenvironmentmethod(executionenv, :beforescenario, context, scenario)

afterscenario(executionenv::FromSourceExecutionEnvironment,
              context::StepDefinitionContext,
              scenario::Gherkin.Scenario) = invokeenvironmentmethod(executionenv, :afterscenario, context, scenario)

beforefeature(executionenv::FromSourceExecutionEnvironment,
              feature::Gherkin.Feature) = invokeenvironmentmethod(executionenv, :beforefeature, feature)

afterfeature(executionenv::FromSourceExecutionEnvironment,
             feature::Gherkin.Feature) = invokeenvironmentmethod(executionenv, :afterfeature, feature)

beforeall(executionenv::FromSourceExecutionEnvironment) = invokeenvironmentmethod(executionenv, :beforeall)

afterall(executionenv::FromSourceExecutionEnvironment) = invokeenvironmentmethod(executionenv, :afterall)


text:
function version_this_pkg()
    v = PkgVersion.Version(LVServer)
    return (; major=v.major, minor=v.minor, patch=v.patch)
end

function trunc_v(v)
    return VersionNumber(v.major, v.minor, v.patch)
end

"""
    check_pkg_versions(;LV_v, LVServer_current, LVServer_min)

Compare the (LV) client and (Julia) server package software versions. \
Error on incompatible versions, warn if not the most recent version known to the opposite side.
# Arguments
- `LV_v`: version of `LabVIEW0` on the client as reported by `LV2Julia_core.lvlib`
- `LVServer_current`: most current `LVServer` version known to `LabVIEW0`/`LV2Julia_core.lvlib`
- `LVServer_min`: the least `LVServer` version compatible with the `LabVIEW0`/`LV2Julia_core.lvlib` client
# Examples
```julia-repl

julia> lvnt = (; LV_v=(major=0, minor=2, patch=1), LVServer_current=(major=1, minor=5, patch=6), LVServer_min=(major=0, minor=2, patch=0));
julia> check_pkg_versions(; lvnt...)
```
"""
function check_pkg_versions(;LV_v, LVServer_current, LVServer_min)

    LV_v = VersionNumber(NamedTuple(LV_v)...)
    LVServer_min = VersionNumber(NamedTuple(LVServer_min)...)
    LVServer_remote_current = VersionNumber(NamedTuple(LVServer_current)...)

    LVServer_localv = PkgVersion.Version(LVServer)

    # @show LV_v LVServer_localv LVServer_min LVServer_remote_current

    LVServer_incompat =  LVServer_localv < LVServer_min
    LV_incompat = LV_v < LabVIEW0_LOCAL_MIN

    LVServer_OK = trunc_v(LVServer_localv) >= LVServer_remote_current
    LV_OK = LV_v >= trunc_v(LabVIEW0_LOCAL_CURRENT)

    if LVServer_incompat
        err_msg = "The LVServer package version on this Julia server is $LVServer_localv. \
        At least $LVServer_min is expected by LabVIEW0 client."
    end

    if LV_incompat
        err_msg = "The version of LabVIEW0 client (or of the currently used LV2Julia_core.lvlib) is $LV_v. \
        At least $LabVIEW0_LOCAL_MIN is expected by LVServer package version on this Julia server"
    end

    if (LVServer_incompat | LV_incompat)
        @warn err_msg # for local display on the server
        throw(ErrorException(err_msg)) # returned to the LV client
    end

    if !LVServer_OK
        warn_msg = "The LVServer package version on this Julia server is $LVServer_localv. \
        According to LabVIEW0 client, $LVServer_remote_current is available, please consider LVServer update."
    end

    if !LV_OK
        warn_msg = "The version of LabVIEW0 client (or of the currently used LV2Julia_core.lvlib) is $LV_v. \
        According to LVServer package on this Julia server, $LabVIEW0_LOCAL_CURRENT  is available, \
        please consider LabVIEW0 update"
    end

    all_current = (LVServer_OK & LV_OK)
    if !all_current
        @warn warn_msg # for local display on the server
    end
    v = LVServer_localv
    return (;all_current, LVServer_v=(; major=v.major, minor=v.minor, patch=v.patch))
end

versionfunctions = (; version_this_pkg, check_pkg_versions,)

# const LabVIEW0_LOCAL_CURRENT = v"0.1.1"
# const LabVIEW0_LOCAL_MIN = v"0.1.0"
# version = "0.2.2"
#
# lvnt_e1 = (; LV_v=(major=0, minor=1, patch=1), LVServer_current=(major=4, minor=5, patch=6), LVServer_min=(major=1, minor=3, patch=0));
# LVServer below min
#
# lvnt_e2 = (; LV_v=(major=0, minor=0, patch=3), LVServer_current=(major=0, minor=2, patch=2), LVServer_min=(major=0, minor=2, patch=0));
# LabVIEW0 client below min
#
# lvnt_w1 = (; LV_v=(major=0, minor=1, patch=0), LVServer_current=(major=0, minor=2, patch=1), LVServer_min=(major=0, minor=2, patch=0));
# LabVIEW0 client not the most current
#
# lvnt_w2 = (; LV_v=(major=0, minor=1, patch=1), LVServer_current=(major=0, minor=2, patch=5), LVServer_min=(major=0, minor=2, patch=0));
# LVServer not the most current
#
# lvnt_ok = (; LV_v=(major=0, minor=1, patch=1), LVServer_current=(major=0, minor=2, patch=2), LVServer_min=(major=0, minor=2, patch=0));
# lvjs = JSON3.write(lvnt);
# lvnt1 = lvnt1=JSON3.read(lvjs);


text:
# PATH uses the Float64 value 1e20 to represent +infinity.
const INFINITY = 1e20

###
### License.h
###

function c_api_License_SetString(license::String)
    return ccall((:License_SetString, PATH_SOLVER), Cint, (Ptr{Cchar},), license)
end

###
### Output_Interface.h
###

const c_api_Output_Log     = Cint(1 << 0)
const c_api_Output_Status  = Cint(1 << 1)
const c_api_Output_Listing = Cint(1 << 2)

mutable struct OutputInterface
    output_data::Ptr{Cvoid}
    print::Ptr{Cvoid}
    flush::Ptr{Cvoid}
end

function _c_flush(data::Ptr{Cvoid}, mode::Cint)
    io = unsafe_pointer_to_objref(data)::IO
    flush(io)
    return
end

function _c_print(data::Ptr{Cvoid}, mode::Cint, msg::Ptr{Cchar})
    if (
        mode & c_api_Output_Log == c_api_Output_Log ||
        # TODO(odow): decide whether to print the Output_Status. It has a lot of
        # information...
        # mode & c_api_Output_Status == c_api_Output_Status ||
        mode & c_api_Output_Listing == c_api_Output_Listing
    )
        io = unsafe_pointer_to_objref(data)::IO
        print(io, unsafe_string(msg))
    end
    return
end

function OutputInterface(io)
    _C_PRINT = @cfunction(_c_print, Cvoid, (Ptr{Cvoid}, Cint, Ptr{Cchar}))
    _C_FLUSH = @cfunction(_c_flush, Cvoid, (Ptr{Cvoid}, Cint))
    return OutputInterface(pointer_from_objref(io), _C_PRINT, _C_FLUSH)
end

function c_api_Output_SetInterface(o::OutputInterface)
    return ccall(
        (:Output_SetInterface, PATH_SOLVER), Cvoid, (Ref{OutputInterface},), o
    )
end

###
### Options.h
###

mutable struct Options
    ptr::Ptr{Cvoid}
    function Options(ptr::Ptr{Cvoid})
        o = new(ptr)
        finalizer(c_api_Options_Destroy, o)
        return o
    end
end

Base.cconvert(::Type{Ptr{Cvoid}}, x::Options) = x
Base.unsafe_convert(::Type{Ptr{Cvoid}}, x::Options) = x.ptr

function c_api_Options_Create()
    ptr = ccall((:Options_Create, PATH_SOLVER), Ptr{Cvoid}, ())
    return Options(ptr)
end

function c_api_Options_Destroy(o::Options)
    return ccall((:Options_Destroy, PATH_SOLVER), Cvoid, (Ptr{Cvoid},), o)
end

function c_api_Options_Default(o::Options)
    return ccall((:Options_Default, PATH_SOLVER), Cvoid, (Ptr{Cvoid},), o)
end

function c_api_Options_Display(o::Options)
    return ccall((:Options_Display, PATH_SOLVER), Cvoid, (Ptr{Cvoid},), o)
end

function c_api_Options_Read(o::Options, filename::String)
    return ccall(
        (:Options_Read, PATH_SOLVER), Cvoid, (Ptr{Cvoid}, Ptr{Cchar}), o, filename
    )
end

function c_api_Path_AddOptions(o::Options)
    return ccall((:Path_AddOptions, PATH_SOLVER), Cvoid, (Ptr{Cvoid},), o)
end

###
### MCP_Interface.h
###

mutable struct InterfaceData
    n::Cint
    nnz::Cint
    F::Function
    J::Function
    lb::Vector{Cdouble}
    ub::Vector{Cdouble}
    z::Vector{Cdouble}
    variable_names::Vector{String}
    constraint_names::Vector{String}
end

function _c_problem_size(
    id_ptr::Ptr{Cvoid}, n_ptr::Ptr{Cint}, nnz_ptr::Ptr{Cint}
)
    id_data = unsafe_pointer_to_objref(id_ptr)::InterfaceData
    n = unsafe_wrap(Array{Cint}, n_ptr, 1)
    n[1] = id_data.n
    nnz = unsafe_wrap(Array{Cint}, nnz_ptr, 1)
    nnz[1] = id_data.nnz
    return
end

function _c_bounds(
    id_ptr::Ptr{Cvoid},
    n::Cint,
    z_ptr::Ptr{Cdouble},
    lb_ptr::Ptr{Cdouble},
    ub_ptr::Ptr{Cdouble},
)
    id_data = unsafe_pointer_to_objref(id_ptr)::InterfaceData
    z = unsafe_wrap(Array{Cdouble}, z_ptr, n)
    lb = unsafe_wrap(Array{Cdouble}, lb_ptr, n)
    ub = unsafe_wrap(Array{Cdouble}, ub_ptr, n)
    for i = 1:n
        z[i] = id_data.z[i]
        lb[i] = id_data.lb[i]
        ub[i] = id_data.ub[i]
    end
    return
end

function _c_function_evaluation(
    id_ptr::Ptr{Cvoid}, n::Cint, x_ptr::Ptr{Cdouble}, f_ptr::Ptr{Cdouble}
)
    id_data = unsafe_pointer_to_objref(id_ptr)::InterfaceData
    x = unsafe_wrap(Array{Cdouble}, x_ptr, n)
    f = unsafe_wrap(Array{Cdouble}, f_ptr, n)
    err = id_data.F(n, x, f)
    return err
end

function _c_jacobian_evaluation(
    id_ptr::Ptr{Cvoid},
    n::Cint,
    x_ptr::Ptr{Cdouble},
    wantf::Cint,
    f_ptr::Ptr{Cdouble},
    nnz_ptr::Ptr{Cint},
    col_ptr::Ptr{Cint},
    len_ptr::Ptr{Cint},
    row_ptr::Ptr{Cint},
    data_ptr::Ptr{Cdouble}
)
    id_data = unsafe_pointer_to_objref(id_ptr)::InterfaceData
    x = unsafe_wrap(Array{Cdouble}, x_ptr, n)
    err = Cint(0)
    if wantf > 0
        f = unsafe_wrap(Array{Cdouble}, f_ptr, n)
        err += id_data.F(n, x, f)
    end
    nnz = unsafe_wrap(Array{Cint}, nnz_ptr, 1)
    col = unsafe_wrap(Array{Cint}, col_ptr, n)
    len = unsafe_wrap(Array{Cint}, len_ptr, n)
    row = unsafe_wrap(Array{Cint}, row_ptr, nnz[1])
    data = unsafe_wrap(Array{Cdouble}, data_ptr, nnz[1])
    err += id_data.J(n, nnz[1], x, col, len, row, data)
    nnz[1] = sum(len)
    return err
end

function _c_variable_name(
    id_ptr::Ptr{Cvoid},
    i::Cint,
    buf_ptr::Ptr{UInt8},
    buf_size::Cint
)
    id_data = unsafe_pointer_to_objref(id_ptr)::InterfaceData
    name = collect(codeunits(id_data.variable_names[i]))
    n = length(name)
    resize!(name, buf_size)
    name[min(buf_size, n + 1):end] .= UInt8('\0')
    GC.@preserve name begin
        unsafe_copyto!(buf_ptr, pointer(name), buf_size)
    end
    return
end

function _c_constraint_name(
    id_ptr::Ptr{Cvoid},
    i::Cint,
    buf_ptr::Ptr{UInt8},
    buf_size::Cint
)
    id_data = unsafe_pointer_to_objref(id_ptr)::InterfaceData
    name = collect(codeunits(id_data.constraint_names[i]))
    n = length(name)
    resize!(name, buf_size)
    name[min(buf_size, n + 1):end] .= UInt8('\0')
    GC.@preserve name begin
        unsafe_copyto!(buf_ptr, pointer(name), buf_size)
    end
    return
end

"""
    MCP_Interface

A storage struct that is used to pass problem-specific functions to PATH.
"""
mutable struct MCP_Interface
    interface_data::Ptr{Cvoid}
    problem_size::Ptr{Cvoid}
    bounds::Ptr{Cvoid}
    function_evaluation::Ptr{Cvoid}
    jacobian_evaluation::Ptr{Cvoid}
    # TODO(odow): the .h files I have don't include the hessian evaluation in
    #             MCP_Interface, but Standalone_Path.c includes it. Ask M.
    #             Ferris to look at the source.
    # Answer: there is an #ifdef to turn it on or off. In the GAMS builds it
    #         appears to be on, but we should be careful when updating PATH
    #         versions.
    hessian_evaluation::Ptr{Cvoid}
    start::Ptr{Cvoid}
    finish::Ptr{Cvoid}
    variable_name::Ptr{Cvoid}
    constraint_name::Ptr{Cvoid}
    basis::Ptr{Cvoid}

    function MCP_Interface(interface_data::InterfaceData)
        _C_PROBLEM_SIZE = @cfunction(
            _c_problem_size,
            Cvoid,
            (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cint})
        )
        _C_BOUNDS = @cfunction(
            _c_bounds,
            Cvoid,
            (Ptr{Cvoid}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble})
        )
        _C_FUNCTION_EVALUATION = @cfunction(
            _c_function_evaluation,
            Cint,
            (Ptr{Cvoid}, Cint, Ptr{Cdouble}, Ptr{Cdouble})
        )
        _C_JACOBIAN_EVALUATION = @cfunction(
            _c_jacobian_evaluation,
            Cint,
            (
                Ptr{Cvoid}, Cint, Ptr{Cdouble}, Cint, Ptr{Cdouble},
                Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}
            )
        )
        if isempty(interface_data.variable_names)
            _C_VARIABLE_NAME = C_NULL
        else
            _C_VARIABLE_NAME = @cfunction(
                _c_variable_name,
                Cvoid,
                (Ptr{Cvoid}, Cint, Ptr{Cuchar}, Cint)
            )
        end
        if isempty(interface_data.constraint_names)
            _C_CONSTRAINT_NAME = C_NULL
        else
            _C_CONSTRAINT_NAME = @cfunction(
                _c_constraint_name,
                Cvoid,
                (Ptr{Cvoid}, Cint, Ptr{Cuchar}, Cint)
            )
        end

        return new(
            pointer_from_objref(interface_data),
            _C_PROBLEM_SIZE,
            _C_BOUNDS,
            _C_FUNCTION_EVALUATION,
            _C_JACOBIAN_EVALUATION,
            C_NULL,
            C_NULL,  # See TODO note in definition of fields above.
            C_NULL,
            _C_VARIABLE_NAME,
            _C_CONSTRAINT_NAME,
            C_NULL
        )
    end
end

mutable struct MCP
    n::Int
    ptr::Ptr{Cvoid}
    id_data::Union{Nothing, InterfaceData}
    function MCP(n::Int, ptr::Ptr{Cvoid})
        m = new(n, ptr, nothing)
        finalizer(c_api_MCP_Destroy, m)
        return m
    end
end

Base.cconvert(::Type{Ptr{Cvoid}}, x::MCP) = x
Base.unsafe_convert(::Type{Ptr{Cvoid}}, x::MCP) = x.ptr

function c_api_MCP_Create(n::Int, nnz::Int)
    ptr = ccall((:MCP_Create, PATH_SOLVER), Ptr{Cvoid}, (Cint, Cint), n, nnz)
    return MCP(n, ptr)
end

function c_api_MCP_Destroy(m::MCP)
    if m.ptr === C_NULL
        return
    end
    ccall((:MCP_Destroy, PATH_SOLVER), Cvoid, (Ptr{Cvoid},), m)
    return
end

function c_api_MCP_SetInterface(m::MCP, interface::MCP_Interface)
    ccall(
        (:MCP_SetInterface, PATH_SOLVER),
        Cvoid,
        (Ptr{Cvoid}, Ref{MCP_Interface}),
        m, interface
    )
    return
end

function c_api_MCP_GetX(m::MCP)
    ptr = ccall((:MCP_GetX, PATH_SOLVER), Ptr{Cdouble}, (Ptr{Cvoid},), m)
    return copy(unsafe_wrap(Array{Cdouble}, ptr, m.n))
end

###
### Types.h
###

@enum(
    MCP_Termination,
    MCP_Solved = 1,
    MCP_NoProgress,
    MCP_MajorIterationLimit,
    MCP_MinorIterationLimit,
    MCP_TimeLimit,
    MCP_UserInterrupt,
    MCP_BoundError,
    MCP_DomainError,
    MCP_Infeasible,
    MCP_Error,
    MCP_LicenseError,
    MCP_OK
)

mutable struct Information
    # Double residual;	     /* Value of residual at final point             */
    residual::Cdouble
    # Double distance;	     /* Distance between initial and final point     */
    distance::Cdouble
    # Double steplength;	     /* Steplength taken                             */
    steplength::Cdouble
    # Double total_time;	     /* Amount of time spent in the code             */
    total_time::Cdouble

    # Double maximum_distance;   /* Maximum distance from init point allowed     */
    maximum_distance::Cdouble

    # Int major_iterations;	     /* Major iterations taken                       */
    major_iterations::Cint
    # Int minor_iterations;	     /* Minor iterations taken                       */
    mainor_iterations::Cint
    # Int crash_iterations;	     /* Crash iterations taken                       */
    crash_iterations::Cint
    # Int function_evaluations;  /* Function evaluations performed               */
    function_evaluations::Cint
    # Int jacobian_evaluations;  /* Jacobian evaluations performed               */
    jacobian_evaluations::Cint
    # Int gradient_steps;	     /* Gradient steps taken                         */
    gradient_steps::Cint
    # Int restarts;		     /* Restarts used                                */
    restarts::Cint

    # Int generate_output;       /* Mask where output can be displayed.          */
    generate_output::Cint
    # Int generated_output;      /* Mask where output displayed.                 */
    generated_output::Cint

    # Boolean forward;	     /* Move forward?                                */
    forward::Bool
    # Boolean backtrace;	     /* Back track?                                  */
    backtrace::Bool
    # Boolean gradient;	     /* Take gradient step?                          */
    gradient::Bool

    # Boolean use_start;         /* Use the starting point provided?             */
    use_start::Bool
    # Boolean use_basics;        /* Use the basis provided?                      */
    use_basics::Bool

    # Boolean used_start;        /* Was the starting point given used?           */
    used_start::Bool
    # Boolean used_basics;       /* Was the initial basis given used?            */
    used_basics::Bool

    function Information(;
        use_start::Bool = true,
    )
        return new(
            0.0, 0.0, 0.0, 0.0, 0.0,
            0, 0, 0, 0, 0, 0, 0,
            0, 0,
            false, false, false,
            use_start, false,
            false, false
        )
    end
end

###
### Path.h
###

"""
    c_api_Path_CheckLicense(n::Int, nnz::Int)

Check that the current license (stored in the environment variable
`PATH_LICENSE_STRING` if present) is valid for problems with `n` variables and
`nnz` non-zeros in the Jacobian.

Returns a nonzero value on successful completion, and a zero value on failure.
"""
function c_api_Path_CheckLicense(n::Int, nnz::Int)
    return ccall((:Path_CheckLicense, PATH_SOLVER), Cint, (Cint, Cint), n, nnz)
end

"""
    c_api_Path_Version()

Return a string of the PATH version.
"""
function c_api_Path_Version()
    ptr = ccall((:Path_Version, PATH_SOLVER), Ptr{Cchar}, ())
    return unsafe_string(ptr)
end

"""
    c_api_Path_Solve(m::MCP, info::Information)

Returns a MCP_Termination status.
"""
function c_api_Path_Solve(m::MCP, info::Information)
    return ccall(
        (:Path_Solve, PATH_SOLVER), Cint, (Ptr{Cvoid}, Ref{Information}), m, info
    )
end

###
### Standalone interface
###

"""
    solve_mcp(
        F::Function,
        J::Function
        lb::Vector{Cdouble},
        ub::Vector{Cdouble},
        z::Vector{Cdouble};
        nnz::Int = length(lb)^2,
        variable_name::Vector{String}=String[],
        constraint_name::Vector{String}=String[],
        silent::Bool = false,
        kwargs...
    )

Mathematically, the mixed complementarity problem is to find an x such that
for each i, at least one of the following hold:

   1.  F_i(x) = 0, lb_i <= (x)_i <= ub_i
   2.  F_i(x) > 0, (x)_i = lb_i
   3.  F_i(x) < 0, (x)_i = ub_i

where F is a given function from R^n to R^n, and lb and ub are prescribed
lower and upper bounds.

`z` is an initial starting point for the search.

`F` is a function `F(n::Cint, x::Vector{Cdouble}, f::Vector{Cdouble})` that
should calculate the function F(x) and store the result in `f`.
"""
function solve_mcp(
    F::Function,
    J::Function,
    lb::Vector{Cdouble},
    ub::Vector{Cdouble},
    z::Vector{Cdouble};
    nnz::Int = length(lb)^2,
    variable_names::Vector{String} = String[],
    constraint_names::Vector{String} = String[],
    silent::Bool = false,
    kwargs...
)
    @assert length(z) == length(lb) == length(ub)
    if !silent && isimmutable(stdout)
        @warn(
            "Disabling output because PATH.jl does not support writing to " *
            "this `stdout`. Are you using Atom-Hydrogen or Jupyter? See: " *
            "https://github.com/odow/PATH.jl/issues/10."
        )
        silent = true
    end
    out_io = silent ? IOBuffer() : stdout
    GC.@preserve out_io begin
    c_api_Output_SetInterface(OutputInterface(out_io))

    n = length(z)
    if n == 0
        return MCP_Solved, nothing, nothing
    end

    # Convert `Int` to `Float64` for check to avoid overflow.
    dnnz = min(1.0 * nnz, 1.0 * n^2)
    if dnnz > typemax(Cint)
        return MCP_Error, nothing, nothing
    end
    nnz = Int(dnnz + 1)

    o = c_api_Options_Create()
    c_api_Path_AddOptions(o)
    c_api_Options_Default(o)
    m = c_api_MCP_Create(n, nnz)
    m.id_data = InterfaceData(
        Cint(n), Cint(nnz),
        F, J, lb, ub, z,
        variable_names, constraint_names
    )
    m_interface = MCP_Interface(m.id_data)
    c_api_MCP_SetInterface(m, m_interface)
    if length(kwargs) > 0
        mktemp() do path, io
            println(io, "* Automatically generated by PATH.jl. Do not edit.")
            for (key, val) in kwargs
                println(io, key, " ", val)
            end
            close(io)
            c_api_Options_Read(o, path)
        end
    end
    c_api_Options_Display(o)
    info = Information(use_start = true)
    status = c_api_Path_Solve(m, info)
    end  # GC.@preserve
    X = c_api_MCP_GetX(m)
    # TODO(odow): I don't know why, but manually calling MCP_Destroy was
    # necessary to avoid a segfault on Julia 1.0 when using LUSOL. I guess it's
    # something to do with the timing of when things need to get freed on the
    # PATH side? i.e., MCP_Destroy before other things?
    c_api_MCP_Destroy(m)
    m.ptr = C_NULL
    return MCP_Termination(status), X, info
end

function _linear_function(M::AbstractMatrix, q::Vector)
    if size(M, 1) != size(M, 2)
        error("M not square! size = $(size(M))")
    elseif size(M, 1) != length(q)
        error("q is wrong shape. Expected $(size(M, 1)), got $(length(q)).")
    end
    return (n::Cint, x::Vector{Cdouble}, f::Vector{Cdouble}) -> begin
        f .= M * x .+ q
        return Cint(0)
    end
end

function _linear_jacobian(M::SparseArrays.SparseMatrixCSC{Cdouble, Cint})
    if size(M, 1) != size(M, 2)
        error("M not square! size = $(size(M))")
    end
    return (
        n::Cint,
        nnz::Cint,
        x::Vector{Cdouble},
        col::Vector{Cint},
        len::Vector{Cint},
        row::Vector{Cint},
        data::Vector{Cdouble}
    ) -> begin
        @assert n == length(x) == length(col) == length(len) == size(M, 1)
        @assert nnz == length(row) == length(data)
        @assert nnz >= SparseArrays.nnz(M)
        for i = 1:n
            col[i] = M.colptr[i]
            len[i] = M.colptr[i + 1] - M.colptr[i]
        end
        for (i, v) in enumerate(SparseArrays.rowvals(M))
            row[i] = v
        end
        for (i, v) in enumerate(SparseArrays.nonzeros(M))
            data[i] = v
        end
        return Cint(0)
    end
end

"""
    solve_mcp(;
        M::SparseArrays.SparseMatrixCSC{Cdouble, Cint},
        q::Vector{Cdouble},
        lb::Vector{Cdouble},
        ub::Vector{Cdouble},
        z::Vector{Cdouble};
        kwargs...
    )

Mathematically, the mixed complementarity problem is to find an x such that
for each i, at least one of the following hold:

   1.  F_i(x) = 0, lb_i <= (x)_i <= ub_i
   2.  F_i(x) > 0, (x)_i = lb_i
   3.  F_i(x) < 0, (x)_i = ub_i

where F is a function `F(x) = M * x + q` from R^n to R^n, and lb and ub are
prescribed lower and upper bounds.

`z` is an initial starting point for the search.
"""
function solve_mcp(
    M::SparseArrays.SparseMatrixCSC{Cdouble, Cint},
    q::Vector{Cdouble},
    lb::Vector{Cdouble},
    ub::Vector{Cdouble},
    z::Vector{Cdouble};
    kwargs...
)
    return solve_mcp(
        _linear_function(M, q),
        _linear_jacobian(M),
        lb,
        ub,
        z;
        nnz = SparseArrays.nnz(M),
        kwargs...
    )
end


text:
module Copeland

function level_of_dominance(v1::Array{Int,1}, v2::Array{Int,1})::Int64
    lod = zero(Int64)
    n = length(v1)
    for i = 1:n
        if v1[i] < v2[i]
            lod += 1
        end
    end
    return lod
end

function dominance_scores(ordering_mat::Array{Int, 2})::Array{Int, 2}
    n, _ = size(ordering_mat)
    lev_dom = zeros(Int64, n, n)
    for i in 1:n
        for j in 1:n
            lev_dom[i, j] = level_of_dominance(ordering_mat[i,:], ordering_mat[j,:])
        end
    end
    return lev_dom
end

function winloss_scores(dommat::Array{Int, 2})::Array{Int, 2}
    n, _ = size(dommat)
    winlossmat = zeros(Int64, n, n)
    for i in 1:n
        for j in 1:n
            winlossmat[i, j] = Int(sign(dommat[i, j] - dommat[j, i]))
        end
    end
    return winlossmat
end

function copeland(ordering_mat::Array{Int, 2})::Array{Int, 1}
    winlosses = ordering_mat |> dominance_scores |> winloss_scores
    n, _ = size(winlosses)
    scores = map(i -> Int(sum(winlosses[i, :])), 1:n)
    return scores
end

end # end of module

max_stars_repo_name: jbytecode/jMCDM

max_issues_repo_name: jbytecode/jMCDM

max_forks_repo_name: jbytecode/jMCDM


text:
"""
    set_value(p::ParameterRef, val::Real)

Sets the parameter `p` to the new value `val`.
"""
function JuMP.set_value(p::ParameterRef, val::Real)
    params = _getparamdata(p)::_ParameterData
    params.sync = false
    params.future_values[index(p)] = val
    return val
end

"""
    value(p::ParameterRef)

Return the current value of the parameter `p`.
"""
function JuMP.value(p::ParameterRef)
    data = _getparamdata(p)::_ParameterData
    return data.future_values[index(p)]
end

# interface continues

JuMP.owner_model(p::ParameterRef) = p.model

struct ParameterNotOwned <: Exception
    parameter::ParameterRef
end

function JuMP.check_belongs_to_model(p::ParameterRef, model::AbstractModel)
    if owner_model(p) !== model
        throw(ParameterNotOwned(p))
    end
end

Base.iszero(::ParameterRef) = false
Base.copy(p::ParameterRef) = ParameterRef(p.ind, p.model)

"""
    is_valid(model::Model, parameter::ParameterRef)

Return `true` if `parameter` refers to a valid parameter in `model`.
"""
function JuMP.is_valid(model::Model, parameter::ParameterRef)
    return model === owner_model(parameter)
end

# The default hash is slow. It's important for the performance of AffExpr to
# define our own.
# https://github.com/JuliaOpt/MathOptInterface.jl/issues/234#issuecomment-366868878
function Base.hash(p::ParameterRef, h::UInt)
    return hash(objectid(owner_model(p)), hash(p.ind, h))
end
function Base.isequal(p1::ParameterRef, p2::ParameterRef)
    return owner_model(p1) === owner_model(p2) && p1.ind == p2.ind
end

JuMP.name(p::ParameterRef) = get(_getparamdata(p).names, p, "")
JuMP.set_name(p::ParameterRef, s::String) = _getparamdata(p).names[p] = s

function parameter_by_name(model::Model, name::String)
    # can be improved with a lazy rev_names map
    dict = _getparamdata(model).names
    for (par, n) in dict
        if n == name
            return par
        end
    end
    return nothing
end


text:
module MIME

using Hiccup
using ..Chess

export html

const BOARD_SIZE = 280
const DARK_SQUARE_COLOR = "#a4704b"
const LIGHT_SQUARE_COLOR = "#ffd7a6"
const HIGHLIGHT_COLOR = "#47d18b"

function squarecolor(file, rank)
    if iseven(file + rank)
        LIGHT_SQUARE_COLOR
    else
        DARK_SQUARE_COLOR
    end
end

function imgurl(p::Piece)::String
    c = pcolor(p)
    t = ptype(p)
    "https://raw.githubusercontent.com/romstad/Chess.jl/master/img/" *
    tochar(c) *
    lowercase(tochar(t)) *
    ".svg"
end

function squarehighlight(s)
    f = file(s).val - 1
    r = rank(s).val - 1
    Node(
        :circle,
        Dict(
            :cx => f + 0.5,
            :cy => r + 0.5,
            :r => 0.4,
            :fill => HIGHLIGHT_COLOR,
            :opacity => 0.5,
        ),
    )
end

function squarehighlights(ss::SquareSet)
    Node(:g, map(squarehighlight, Chess.squares(ss)))
end

function square(file::Int, rank::Int, piece)
    Node(
        :g,
        Node(
            :rect,
            Dict(
                :fill => squarecolor(file, rank),
                :x => file,
                :y => rank,
                :width => 1,
                :height => 1,
            ),
        ),
        if piece == EMPTY
            Node(:g)
        else
            Node(
                :image,
                Dict(
                    Symbol("xlink:href") => imgurl(piece),
                    :x => file,
                    :y => rank,
                    :width => 1,
                    :height => 1,
                ),
            )
        end,
    )
end

function square(s::Square, piece)
    f = file(s).val - 1
    r = rank(s).val - 1
    square(f, r, piece)
end

function squares(board)
    Node(:g, map(s -> square(Square(s), pieceon(board, Square(s))), 1:64))
end

function svg(; board = emptyboard(), highlight = SS_EMPTY)
    Node(
        :svg,
        Dict(
            :style => "float: left; margin-right: 20px",
            :viewBox => "0 0 8 8",
            :width => BOARD_SIZE,
            :height => BOARD_SIZE,
        ),
        [squares(board), squarehighlights(highlight)],
    )
end

function lichesslink(board)
    Node(:a, Dict(:href => lichessurl(board), :target => "_blank"), "Open in lichess")
end

function description(board)
    Node(
        :div,
        [
            Node(:p, sidetomove(board) == WHITE ? "White to move" : "Black to move"),
            board.castlerights == 0 ? "" :
            Node(:p, "Castle rights: " * Chess.castlestring(board)),
            epsquare(board) == SQ_NONE ? "" :
            Node(:p, "En passant square: " * tostring(epsquare(board))),
            Node(:p, lichesslink(board)),
        ],
    )
end

function html(board::Board; highlight = SS_EMPTY)
    Node(
        :div,
        Dict(:class => "chessboard"),
        [svg(board = board, highlight = highlight), description(board)],
    )
end

function html(ss::SquareSet)
    Node(:div, Dict(:class => "chessboard"), svg(highlight = ss))
end

end


text:
using Iris
using NLsolve
using Test
using LinearAlgebra

@testset "Linear Scattering" begin
    @testset "1D" begin
        @testset "PML 2-sided" begin
            ω=40
            sim = Iris.Examples.slab_pml(ω=ω, dx=.00035)
            ls = @test_nowarn HelmholtzLS(sim, ω, 1)
            @test_nowarn @inferred HelmholtzLS(sim, ω, 1)
            f(x) = x.operator
            @test_nowarn @inferred f(ls)
            f(x) = x.equivalent_source
            @test_nowarn @inferred f(ls)
            f(x) = x.j
            @test_nowarn @inferred f(ls)
            f(x) = x.solved
            @test_nowarn @inferred f(ls)
            f(x) = x.converged
            @test_nowarn @inferred f(ls)
            f(x) = x.solution
            @test_nowarn @inferred f(ls)
            f(x) = x.sim
            @test_nowarn @inferred f(ls)
            f(x) = x.ω
            @test_nowarn @inferred f(ls)
            f(x) = x.a
            @test_nowarn @inferred f(ls)
            f(x) = x.H
            @test_nowarn @inferred f(ls)
            f(x) = x.src
            @test_nowarn @inferred f(ls)
            f(x) = x.sol
            @test_nowarn @inferred f(ls)
            es = ls.equivalent_source
            f(es) = es.incoming_mask
            @test_nowarn @inferred f(es)
            f(es) = es.outgoing_mask
            @test_nowarn @inferred f(es)
            f(es) = es.field
            @test_nowarn @inferred f(es)
            f(es) = es.simulation
            @test_nowarn @inferred f(es)
            f(es) = es.ω
            @test_nowarn @inferred f(es)
            f(es) = es.a
            @test_nowarn @inferred f(es)
            f(es) = es.channelflux
            @test_nowarn @inferred f(es)
            sol = ls.solution
            f(sol) = sol.total
            @test_nowarn @inferred f(sol)
            f(sol) = sol.incident
            @test_nowarn @inferred f(sol)
            f(sol) = sol.scattered
            @test_nowarn @inferred f(sol)
            f(sol) = sol.ω
            @test_nowarn @inferred f(sol)
            f(sol) = sol.a
            @test_nowarn @inferred f(sol)
            @test_nowarn @inferred HelmholtzLS(sim, ω, 1)
            @test_nowarn scattering!(ls)
            @test ls.converged[]
            @test ls.helmholtz.A*ls.solution.total≈ls.j
            @test isapprox(ls.solution.scattered[288],0.07314178272991415 + 0.06502554346709624im; atol=1e-3)
        end
        @testset "Matching 2-sided" begin
            ω=40
            sim = Iris.Examples.slab_out(ω=ω, dx=.00035)
            ls = @test_nowarn HelmholtzLS(sim, ω, 1)
            @test_nowarn scattering!(ls)
            @test ls.converged[]
            @test ls.helmholtz.A*ls.solution.total≈ls.j
            @test isapprox(ls.solution.scattered[288],0.07314178272991415 + 0.06502554346709624im; atol=1e-3)
        end
        @testset "PML 1-sided" begin
            ω0 = 40
            lat = Lattices.Cartesian(.00035)
            bnd = Boundary(Interval(-.7,.7),PML{1}(),DirichletBC)
            cell = LatticeDomain(bnd,lat)
            ndcavity = NondispersiveDomain(Interval(-.5,.5),3)
            sim = Simulation(ω0, ndcavity, cell)
            smooth!(sim)
            ls = @test_nowarn HelmholtzLS(sim, ω0, 1)
            @test_nowarn scattering!(ls)
            @test ls.converged[]
            @test ls.helmholtz.A*ls.solution.total≈ls.j
            @test isapprox(ls.solution.scattered[288],-0.1290389888137668 + 0.27289178643785555im; atol=1e-3)
        end
        @testset "Matched 1-sided" begin
            ω0 = 40
            lat = Lattices.Cartesian(.00035)
            bnd = Boundary(Interval(-.7,.7),PML{1}(),MatchedBC{1}(out=1),DirichletBC{2}())
            cell = LatticeDomain(bnd,lat)
            ndcavity = NondispersiveDomain(Interval(-.5,.5),3)
            sim = Simulation(ω0, ndcavity, cell)
            smooth!(sim)
            ls = @test_nowarn HelmholtzLS(sim, ω0, 1)
            @test_nowarn scattering!(ls)
            @test ls.converged[]
            @test ls.helmholtz.A*ls.solution.total≈ls.j
            @test isapprox(ls.solution.scattered[288],-0.1290389888137668 + 0.27289178643785555im; atol=1e-3)
        end
    end
end

@testset "Nonlinear Scattering" begin
    @testset "1D" begin
        @testset "PML 2-sided" begin
            ω=40
            sim = Iris.Examples.slab_pml(ω=ω, dx=.00035, D=.1, γ=2)
            nls = @test_nowarn HelmholtzNLS(sim, ω, 1)
            @test_nowarn @inferred HelmholtzNLS(sim,ω,1)
            f(x) = x.linearscatter
            @test_nowarn @inferred f(nls)
            f(x) = x.ψ
            @test_nowarn @inferred f(nls)
            f(x) = x.residual
            @test_nowarn @inferred f(nls)
            f(x) = x.x
            @test_nowarn @inferred f(nls)
            f(x) = x.fixedpoint
            @test_nowarn @inferred f(nls)
            f(x) = x.lupack
            @test_nowarn @inferred f(nls)
            f(x) = x.a
            @test_nowarn @inferred f(nls)
            f(x) = x.equivalent_source
            @test_nowarn @inferred f(nls)
            f(x) = x.helmholtz
            @test_nowarn @inferred f(nls)
            f(x) = x.operator
            @test_nowarn @inferred f(nls)
            f(x) = x.j
            @test_nowarn @inferred f(nls)
            f(x) = x.simulation
            @test_nowarn @inferred f(nls)
            f(x) = x.solution
            @test_nowarn @inferred f(nls)
            f(x) = x.ω
            @test_nowarn @inferred f(nls)
            f(x) = x.results
            @test_nowarn scattering!(nls)
            @test_nowarn @inferred f(nls)
            nls = @test_nowarn HelmholtzNLS(sim, ω, 1)
            @test_nowarn fixedpoint(nls)
            @test nls.converged[]
            @test nls.helmholtz.A*nls.solution.total≈nls.j
            nls = @test_nowarn HelmholtzNLS(sim, ω, 1)
            @test_nowarn nlsolve(nls)
            @test nls.converged[]
            @test nls.helmholtz.A*nls.solution.total≈nls.j
            nls = @test_nowarn HelmholtzNLS(sim, ω, 1)
            @test_nowarn scattering!(nls)
            @test nls.converged[]
            @test nls.helmholtz.A*nls.solution.total≈nls.j
            @test isapprox(nls.solution.scattered[288],0.1881169235362705 + 0.15566101877014im; atol=1e-3)
            nls = @test_nowarn HelmholtzNLS(sim, ω, 10)
            @test_nowarn scattering!(nls)
            @test nls.converged[]
            @test nls.helmholtz.A*nls.solution.total≈nls.j
            @test isapprox(nls.solution.scattered[288], 1.1192308728340543 + 0.7280855538638448im; atol=1e-3)
        end
        @testset "Matching 2-sided" begin
            ω0=40
            sim = Iris.Examples.slab_out(ω=ω0, dx=.00035, D=.1, γ=2)
            nls = @test_nowarn HelmholtzNLS(sim, ω0, 1)
            @test_nowarn fixedpoint(nls)
            @test nls.converged[]
            @test nls.helmholtz.A*nls.solution.total≈nls.j
            nls = @test_nowarn HelmholtzNLS(sim, ω0, 1)
            @test_nowarn nlsolve(nls)
            @test nls.converged[]
            @test nls.helmholtz.A*nls.solution.total≈nls.j
            nls = @test_nowarn HelmholtzNLS(sim, ω0, 1)
            @test_nowarn scattering!(nls)
            @test nls.converged[]
            @test nls.helmholtz.A*nls.solution.total≈nls.j
            @test isapprox(nls.solution.scattered[288],0.18811652728844638 + 0.15566252089682026im; atol=1e-3)
            nls = @test_nowarn HelmholtzNLS(sim, ω0, 10)
            @test_nowarn scattering!(nls)
            @test nls.converged[]
            @test nls.helmholtz.A*nls.solution.total≈nls.j
            @test isapprox(nls.solution.scattered[288],1.1192537618187042 + 0.7280915535384698im; atol=1e-3)
        end
        @testset "PML 1-sided" begin
            ω0 = 40
            lat = Lattices.Cartesian(.00035)
            bnd = Boundary(Interval(-1,1),PML{1}(),MatchedBC{1}(out=1),DirichletBC{2}())
            cell = LatticeDomain(bnd,lat)
            ndcavity = NondispersiveDomain(Interval(-.5,.5),3)
            dcavity = DispersiveDomain(Interval(-.3,.3),TwoLevelSystem(ω0,.1,2))
            sim = Simulation(ω0, ndcavity, cell, dcavity)
            smooth!(sim)
            nls = @test_nowarn HelmholtzNLS(sim, ω0, 1)
            @test_nowarn fixedpoint(nls)
            @test nls.converged[]
            @test nls.helmholtz.A*nls.solution.total≈nls.j
            nls = @test_nowarn HelmholtzNLS(sim, ω0, 1)
            @test_nowarn nlsolve(nls)
            @test nls.converged[]
            @test nls.helmholtz.A*nls.solution.total≈nls.j
            nls = @test_nowarn HelmholtzNLS(sim, ω0, 1)
            @test_nowarn scattering!(nls)
            @test nls.converged[]
            @test nls.helmholtz.A*nls.solution.total≈nls.j
            @test isapprox(nls.solution.scattered[288],-0.07570911698903976 + 0.3176555836644436im; atol=1e-3)
            nls = @test_nowarn HelmholtzNLS(sim, ω0, 10)
            @test_nowarn scattering!(nls)
            @test nls.converged[]
            @test nls.helmholtz.A*nls.solution.total≈nls.j
            @test isapprox(nls.solution.scattered[288],0.5081023860501906 + 2.6581334068616136im; atol=1e-3)
        end
        @testset "Matched 1-sided" begin
            ω0 = 40
            lat = Lattices.Cartesian(.00035)
            bnd = Boundary(Interval(-1,1),PML{1}(),MatchedBC{1}(out=1),DirichletBC{2}())
            cell = LatticeDomain(bnd,lat)
            ndcavity = NondispersiveDomain(Interval(-.5,.5),3)
            dcavity = DispersiveDomain(Interval(-.3,.3),TwoLevelSystem(ω0,.1,2))
            sim = Simulation(ω0, ndcavity, cell, dcavity)
            smooth!(sim)
            nls = @test_nowarn HelmholtzNLS(sim, ω0, 1)
            @test_nowarn fixedpoint(nls)
            @test nls.converged[]
            @test nls.helmholtz.A*nls.solution.total≈nls.j
            nls = @test_nowarn HelmholtzNLS(sim, ω0, 1)
            @test_nowarn nlsolve(nls)
            @test nls.converged[]
            @test nls.helmholtz.A*nls.solution.total≈nls.j
            nls = @test_nowarn HelmholtzNLS(sim, ω0, 1)
            @test_nowarn scattering!(nls)
            @test nls.converged[]
            @test nls.helmholtz.A*nls.solution.total≈nls.j
            @test isapprox(nls.solution.scattered[288],-0.07570911698903976 + 0.3176555836644436im; atol=1e-3)
            nls = @test_nowarn HelmholtzNLS(sim, ω0, 10)
            @test_nowarn scattering!(nls)
            @test nls.converged[]
            @test nls.helmholtz.A*nls.solution.total≈nls.j
            @test isapprox(nls.solution.scattered[288],0.5081023860501906 + 2.6581334068616136im; atol=1e-3)
        end
    end

    @testset "SPA" begin
        ω=100.4
        sim = Iris.Examples.slab_pml(ω=ω, dx=.00035, D=-1, γ=2)
        nls = @test_nowarn HelmholtzNLS(sim, ω, 50, 10)
        scattering!(nls)
        @test_nowarn @inferred spa(nls; nev=2)
    end
end


text:
function uniform_probabilities(T, npoints::Integer)
    pmin = 1 / (2 * T(npoints))
    return pmin:(1 / T(npoints)):(1 - pmin)
end

as_array(x::AbstractArray) = x
as_array(x) = [x]

missing_to_nan(x::AbstractArray{>:Missing}) = replace(x, missing => NaN)
missing_to_nan(::Missing) = NaN
missing_to_nan(x) = x

# dimension corresponding to parameters
function param_dim(x)
    @assert ndims(x) > 1
    return 1
end

# view of first draw from first chain (i.e. vector of parameters)
function first_draw(x::AbstractArray)
    dims = Base.setindex(ntuple(one, ndims(x)), :, param_dim(x))
    return view(x, dims...)
end

# view of all draws
param_draws(x::AbstractArray, i::Int) = selectdim(x, param_dim(x), i)

# dimensions corresponding to draws and chains
function sample_dims(x::AbstractArray)
    d = param_dim(x)
    return filter(!=(d), ntuple(identity, ndims(x)))
end
sample_dims(::AbstractVector) = Colon()


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""Represents a Glusterfs mount that lasts the lifetime of a pod. Glusterfs volumes do not support ownership management or SELinux relabeling.

    IoK8sApiCoreV1GlusterfsVolumeSource(;
        endpoints=nothing,
        path=nothing,
        readOnly=nothing,
    )

    - endpoints::String : EndpointsName is the endpoint name that details Glusterfs topology. More info: https://examples.k8s.io/volumes/glusterfs/README.md#create-a-pod
    - path::String : Path is the Glusterfs volume path. More info: https://examples.k8s.io/volumes/glusterfs/README.md#create-a-pod
    - readOnly::Bool : ReadOnly here will force the Glusterfs volume to be mounted with read-only permissions. Defaults to false. More info: https://examples.k8s.io/volumes/glusterfs/README.md#create-a-pod
"""
mutable struct IoK8sApiCoreV1GlusterfsVolumeSource <: SwaggerModel
    endpoints::Any # spec type: Union{ Nothing, String } # spec name: endpoints
    path::Any # spec type: Union{ Nothing, String } # spec name: path
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly

    function IoK8sApiCoreV1GlusterfsVolumeSource(;endpoints=nothing, path=nothing, readOnly=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1GlusterfsVolumeSource, Symbol("endpoints"), endpoints)
        setfield!(o, Symbol("endpoints"), endpoints)
        validate_property(IoK8sApiCoreV1GlusterfsVolumeSource, Symbol("path"), path)
        setfield!(o, Symbol("path"), path)
        validate_property(IoK8sApiCoreV1GlusterfsVolumeSource, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        o
    end
end # type IoK8sApiCoreV1GlusterfsVolumeSource

const _property_map_IoK8sApiCoreV1GlusterfsVolumeSource = Dict{Symbol,Symbol}(Symbol("endpoints")=>Symbol("endpoints"), Symbol("path")=>Symbol("path"), Symbol("readOnly")=>Symbol("readOnly"))
const _property_types_IoK8sApiCoreV1GlusterfsVolumeSource = Dict{Symbol,String}(Symbol("endpoints")=>"String", Symbol("path")=>"String", Symbol("readOnly")=>"Bool")
Base.propertynames(::Type{ IoK8sApiCoreV1GlusterfsVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1GlusterfsVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1GlusterfsVolumeSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1GlusterfsVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1GlusterfsVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1GlusterfsVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1GlusterfsVolumeSource)
    (getproperty(o, Symbol("endpoints")) === nothing) && (return false)
    (getproperty(o, Symbol("path")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1GlusterfsVolumeSource }, name::Symbol, val)
end


text:
using LinearAlgebra
using Plots

#%
#% The 'illustrative example' from Section 1.4, Chapter 1.
#%
#% written by John Bardsley 2016.
#%
#% Define the matrix A using the SVD
v1 = [1/sqrt(2) 1/sqrt(2)]';
v2 = [-1/sqrt(2) 1/sqrt(2)]';
s1 = 1;
s2 = 1e-2;
A  = s1*v1*v1'+s2*v2*v2';

#% Define true x and noise free b
x   = [1 1]';
b_e = A*x;

#% Create realization from the data model
nsamp  = 10000;
σ  = 0.1;
b_samp = repeat(b_e,1,nsamp)+σ*randn(2,nsamp);

#figure(1),
  scatter(b_samp[1,:],b_samp[2,:]) #,'k*')
#% Create the corresponding least squares solutions
x_LS = A\b_samp;
#figure(2)
  scatter(x_LS[1,:],x_LS[2,:]) #,'k*')

#% Check ill-posedness inequality (1.28)
ill_posed_lhs = 0.5*(1/s1+1/s2)
ill_posed_rhs = norm(x)/norm(b_e)