text:
module Backward

export value_and_gradient, gradient

import Base: ==, !=, <, <=, >, >=, isless, isequal
import Base: +, *, -, /, exp, sqrt, log, log1p
import Base: sin, cos, tan, sec, csc, cot
import Base: atan
import Base: sum
import Base: convert, promote_rule, zero, one

import AutoDiff.StatsFunctions: logsumexp

mutable struct BADNode{T<:Number} <: Number
    tape::Array{BADNode{T},1}
    value::T
    adj::T
    backprop!
end

function zero(::Type{BADNode{T}}) where {T <: Number}
    BADNode(BADNode{T}[], zero(T), zero(T), (n)->())
end

function zero(x::BADNode{T}) where {T <: Number}
    BADNode(BADNode{T}[], zero(T), zero(T), (n)->())
end

function one(::Type{BADNode{T}}) where {T <: Number}
    BADNode(BADNode{T}[], one(T), zero(T), (n)->())
end

function one(x::BADNode{T}) where {T <: Number}
    BADNode(BADNode{T}[], one(T), zero(T), (n)->())
end

function convert(::Type{BADNode{T}}, x::T) where {T <: Number}
    BADNode(BADNode{T}[], x, zero(T), (n)->())
end

function promote_rule(::Type{BADNode{T}}, ::Type{T}) where {T <: Number}
    BADNode{T}
end

for fn in [Symbol("=="), Symbol("!="), Symbol("<"), Symbol("<="), Symbol(">"), Symbol(">="), :isless, :isequal]
    @eval begin
        $fn(a::BADNode{T}, b::BADNode{T}) where {T <: Number} = $fn(a.value, b.value)
    end
end

function +(x::BADNode{T}, y::BADNode{T}) where {T <: Number}
    t = (length(x.tape) == 0 ? y.tape : x.tape)
    function bp!(n::BADNode{T})
        x.adj += n.adj
        y.adj += n.adj
    end
    n = BADNode(t, x.value + y.value, zero(T), bp!)
    push!(t, n)
    n
end

function -(x::BADNode{T}, y::BADNode{T}) where {T <: Number}
    t = (length(x.tape) == 0 ? y.tape : x.tape)
    function bp!(n::BADNode{T})
        x.adj += n.adj
        y.adj -= n.adj
    end
    n = BADNode(t, x.value - y.value, zero(T), bp!)
    push!(t, n)
    n
end

function -(x::BADNode{T}) where {T <: Number}
    t = x.tape
    function bp!(n::BADNode{T})
        x.adj -= n.adj
    end
    n = BADNode(t, -x.value, zero(T), bp!)
    push!(t, n)
    n
end

function *(x::BADNode{T}, y::BADNode{T}) where {T <: Number}
    t = (length(x.tape) == 0 ? y.tape : x.tape)
    function bp!(n::BADNode{T})
        x.adj += y.value*n.adj
        y.adj += x.value*n.adj
    end
    n = BADNode(t, x.value*y.value, zero(T), bp!)
    push!(t, n)
    n
end

function /(x::BADNode{T}, y::BADNode{T}) where {T <: Number}
    t = (length(x.tape) == 0 ? y.tape : x.tape)
    r = x.value / y.value
    function bp!(n::BADNode{T})
        x.adj += n.adj / y.value
        y.adj -= r / y.value * n.adj
    end
    n = BADNode(t, r, zero(T), bp!)
    push!(t, n)
    n
end

function exp(x::BADNode{T}) where {T <: Number}
    t = x.tape
    r = exp(x.value)
    function bp!(n::BADNode{T})
        x.adj += r*n.adj
    end
    n = BADNode(t, r, zero(T), bp!)
    push!(t, n)
    n
end

function sqrt(x::BADNode{T}) where {T <: Number}
    t = x.tape
    r = sqrt(x.value)
    function bp!(n::BADNode{T})
        x.adj += n.adj / (convert(T, 2)*r)
    end
    n = BADNode(t, r, zero(T), bp!)
    push!(t, n)
    n
end

function sin_deriv(x,y)
    cos(x)
end

function cos_deriv(x,y)
    -sin(x)
end

function tan_deriv(x, y)
    z = sec(x)
    z*z
end

function sec_deriv(x, y)
    y*tan(x)
end

function csc_deriv(x, y)
    -y*cot(x)
end

function cot_deriv(x, y)
    z = csc(x)
    -z*z
end

function atan_deriv(x, y)
    one(x) / (one(x) + x*x)
end

for (f, df) in zip([:sin, :cos, :tan, :sec, :csc, :cot, :atan], [:sin_deriv, :cos_deriv, :tan_deriv, :sec_deriv, :csc_deriv, :cot_deriv, :atan_deriv])
    fstr = uppercasefirst(string(f))
    nodename = Symbol("BADNode$(fstr)")
    @eval begin
        function $f(x::BADNode{T}) where {T <: Number}
            t = x.tape
            r = $f(x.value)
            function bp!(n::BADNode{T})
                x.adj += $df(x.value, r)*n.adj
            end
            n = BADNode(t, r, zero(T), bp!)
            push!(t, n)
            n
        end
    end
end

function atan(y::BADNode{T}, x::BADNode{T}) where {T <: Number}
    t = (length(x.tape) == 0 ? y.tape : x.tape)
    function bp!(n::BADNode{T})
        denom = x.value*x.value + y.value*y.value
        y.adj += n.adj * x.value / denom
        x.adj -= n.adj * y.value / denom
    end
    n = BADNode(t, atan(y.value, x.value), zero(T), bp!)
    push!(t, n)
    n
end

function sum(x::Array{BADNode{T}, N}) where {T <: Number, N}
    t = BADNode{T}[]
    for i in eachindex(x)
        if length(x[i].tape) > 0
            t = x[i].tape
            break
        end
    end

    s = zero(T)
    for i in eachindex(x)
        s += x[i].value
    end

    function bp!(n::BADNode{T})
        for i in eachindex(x)
            x[i].adj += n.adj
        end
    end

    n = BADNode(t, s, zero(T), bp!)
    push!(t, n)
    n
end

function log(x::BADNode{T}) where {T <: Number}
    t = x.tape
    function bp!(n::BADNode{T})
        x.adj += n.adj / x.value
    end
    n = BADNode(t, log(x.value), zero(T), bp!)
    push!(t, n)
    n
end

function log1p(x::BADNode{T}) where {T <: Number}
    function bp!(n::BADNode{T})
        x.adj += n.adj / (one(T) + x.value)
    end
    n = BADNode(x.tape, log1p(x.value), zero(T), bp!)
    push!(x.tape, n)
    n
end

function logsumexp(x::BADNode{T}, y::BADNode{T}) where {T <: Number}
    t = (length(x.tape) == 0 ? y.tape : x.tape)
    r = logsumexp(x.value, y.value)
    function bp!(n::BADNode{T})
        x.adj += n.adj * exp(x.value - r)
        y.adj += n.adj * exp(y.value - r)
    end
    n = BADNode(t, r, zero(T), bp!)
    push!(t, n)
    n
end

function logsumexp(x::Array{BADNode{T}, N}) where {T <: Number, N}
    t = BADNode{T}[]
    for i in eachindex(x)
        if length(x[i].tape) > 0
            t = x[i].tape
            break
        end
    end

    xx = Array{T, N}(undef, size(x))
    for i in eachindex(x)
        xx[i] = x[i].value
    end
    r = logsumexp(xx)

    function bp!(n::BADNode{T})
        for i in eachindex(x)
            x[i].adj += n.adj*exp(x[i].value - r)
        end
    end

    n = BADNode(t, r, zero(T), bp!)
    push!(t, n)
    n
end

function value_and_gradient(f)
    function vgf(x::T) where T<:Number
        t = BADNode{T}[]
        y = BADNode(t, x, zero(T), (n)->())
        push!(t, y)
        r = f(y)
        r.adj = one(T)
        while length(t) > 0
            n = pop!(t)
            n.backprop!(n)
        end
        (r.value, y.adj)
    end

    function vgf(xs::Array{T, N}) where {T <: Number, N}
        t = BADNode{T}[]
        ys = Array{BADNode{T}, N}(undef, size(xs))
        for i in eachindex(ys)
            ys[i] = BADNode(t, xs[i], zero(T), (n)->())
            push!(t, ys[i])
        end
        r = f(ys)
        r.adj = one(T)

        while length(t) > 0
            n = pop!(t)
            n.backprop!(n)
        end

        zs = Array{T, N}(undef, size(xs))
        for i in eachindex(xs)
            zs[i] = ys[i].adj
        end
        (r.value, zs)
    end

    function vgf(x::Vararg{T,N}) where {T<:Number, N}
        t = BADNode{T}[]

        y = Array{BADNode{T}, 1}(undef, length(x))
        for i in eachindex(x)
            y[i] = BADNode(t, x[i], zero(T), (n)->())
            push!(t, y[i])
        end

        r = f(y...)
        r.adj = one(T)

        while length(t) > 0
            n = pop!(t)
            n.backprop!(n)
        end

        g = T[n.adj for n in y]
        (r.value, g)
    end

    vgf
end

function gradient(f)
    vgf = value_and_gradient(f)
    function g(x...)
        vgf(x...)[2]
    end

    g
end

end


text:
using MLStyle
using ExproniconLite
using Documenter
using ANSIColoredPrinters
using DocThemeIndigo

indigo = DocThemeIndigo.install(ExproniconLite)

makedocs(;
    modules=[ExproniconLite],
    authors="Roger-luo <rogerluo.rl18@gmail.com> and contributors",
    repo="https://github.com/Roger-luo/ExproniconLite.jl/blob/{commit}{path}#{line}",
    sitename="ExproniconLite.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Roger-luo.github.io/ExproniconLite.jl",
        assets=String[indigo, "assets/default.css"],
    ),
    pages=[
        "Quick Start" => "index.md",
        "Case Study" => "case.md",
        "API Reference" => "reference.md",
    ],
)

deploydocs(;
    repo="github.com/Roger-luo/ExproniconLite.jl",
)

ext: jl


text:
# Numerical Jacobian for a sparsity pattern thta only involves neighbors
using SparseArrays
"""
    sparse(IJ::Vector{<:CartesianIndex}, v, m, n)

Overloads SparseArrays' `sparse` function to work directly with `CartesianIndex` type of index.
"""
function SparseArrays.sparse(IJ::Vector{<:CartesianIndex}, v, m, n)
    IJ′ = reinterpret(Int, reshape(IJ, 1, :))
    return sparse(view(IJ′, 1, :), view(IJ′, 2, :), v, m, n)
end

using DualNumbers
"""
    Compressor

Type, containing all the required objects for fast comoression and decompression.
"""
struct Compressor
    i
    j
    icmp
    jcmp
    ijcmp
    idcmp
    jdcmp
    nwet
    ncol
    ε_colors
    nz
    col_idx
end

"""
    compress(Jac, C::Compressor)

Compresses the Jacobian into a smaller matrix.
"""
function compress(Jac, C::Compressor)
    out = zeros(eltype(Jac), C.nwet, C.ncol)
    out[C.ijcmp] .= Jac.nzval
    return out
end

"""
    decompress(Jac, C::Compressor)

Decompresses the Jacobian back into its full size.
"""
decompress(smallJac, C::Compressor) = sparse(C.i, C.j, smallJac[C.ijcmp], C.nwet, C.nwet)

"""
    numJac(f, x, C)

Evaluates the Jacobian of `f` at `x` numerically using the `DualNumbers` package.

The evaluation is efficient by limiting the number of evaluations of `f`.
"""
function numJac(f, x, C::Compressor)
    xpert = repeat(x, 1, C.ncol) + C.ε_colors
    smallJac = epsilon.(f(xpert))
    J = decompress(smallJac, C)
    dropzeros!(J)
    return J
end

"""
    adjacency_matrix_indices(wet3d, dist)

Builds the adjacency matrix indices for the full 3d grid (wet and dry).

`dist` is the infitiny-norm distance of adjacent boxes.
Color indices thus vary from `1` to `(1+2dist)^2` boxes are adjacent for each box.
(`dist = 2` required for OCIM1.1 and POP2 and means 125 adjacent boxes.)
"""
function adjacency_matrix_indices(wet3d, dist)
    nlat, nlon, ndepth = size(wet3d)
    L = LinearIndices(wet3d)
    R = CartesianIndices(size(wet3d)) # Range of all cartesian indices
    i1 = first(R) # convenient to create R1 for accessing neighbors in the line below
    RAdj = CartesianIndices(UnitRange.(Tuple(-dist * i1), Tuple(dist * i1))) # Range of cartesian indices of neighbors 3 directions
    L1 = LinearIndices(size(RAdj))
    nR, nRAdj = length(R), length(RAdj)
    # Pre-allocate sub indices
    ljnz = zeros(Int64, nR, nRAdj)
    linz = zeros(Int64, nR)
    for i in R # Loop through all boxes
        li = L[i]
        linz[li] = li # Allocate index of each box
        for k in L1 # Loop through all adjacent boxes
            lj = L[valid_index(i + RAdj[k], nlat, nlon, ndepth)]
            ljnz[L[i], L1[k]] = lj # Allocate index of each adjacent box
        end
    end
    return linz, ljnz
end

"""
    adjacency_matrix(wet3d)

Builds the adjacency matrix for the full 3d grid (wet and dry).

`dist` is the infitiny-norm distance of adjacent boxes.
Color indices thus vary from `1` to `(1+2dist)^2` boxes are adjacent for each box.
(`dist = 2` required for OCIM1.1 and POP2 and means 125 adjacent boxes.)
"""
function adjacency_matrix(wet3d, dist)
    linz, ljnz = adjacency_matrix_indices(wet3d, dist)
    linz = repeat(linz, size(ljnz, 2))
    ljnz = vec(ljnz)
    AM_3d = sparse(linz, ljnz, true)
    return AM_3d
end

"""
    valid_index(j, nlat, nlon, ndepth)

Return a valid cartesian index.

If out of bounds in latitude or depth, return the bound.
If out of bounds in longitude, return the modulo(nlon) value.
(Because longitude is periodic.)
"""
function valid_index(j, nlat, nlon, ndepth)
    y = max(1, min(nlat, j[1]))
    x = mod1(j[2], nlon)
    z = max(1, min(ndepth, j[3]))
    return CartesianIndex(y, x, z)
end

"""
    build_compressor(wet3d, iwet, nwet, dist)

Builds a Compressor object.

A Compressor object contains an adjacency matrix and all the indices required
for compressing and decompressing the numerical Jacobian.
"""
function build_compressor(wet3d, iwet, nwet, dist)
    # Decompression indices
    # We use the indices of the adjacency matrix
    AM_3d = adjacency_matrix(wet3d, dist) # Full wet and dry adjacency adjacency_matrix
    AM = AM_3d[iwet, iwet] # Adjacency matrix (wet points only)
    i, j, _ = findnz(AM) # Sub indices

    # Compression indices
    # We use the color indices
    # The compressed index `icmp` is the same as the decompressed `i`
    # The compressed index `jcmp` is the color index
    icmp = i
    ncol, ε_colors, col_idx = color_groups(wet3d, iwet, nwet, dist)
    jcmp = col_idx[j]

    # matrix for preallocating smallJac (for storage in Compressor)
    #smallJac₀ = sparse(icmp, jcmp, 1.0, nwet, ncol)
    ijcmp = CartesianIndex.(icmp, jcmp)

    # number of non zeros for storage and use in numJac
    nz = length(i)

    # Permuted decompression indices (faster!)
    v = 1:nz # to track the index permutation that happense during compression
    AMv = sparse(icmp, jcmp, v, nwet, ncol)
    _, _, v2 = findnz(AMv) # v2 is the verctor of the index permutation
    idcmp, jdcmp = i[v2], j[v2] # This allows direct use of smallJac.nzval
    # which is much faster than accessing those values via (icmp, jcmp)

    return Compressor(i, j, icmp, jcmp, ijcmp, idcmp, jdcmp, nwet, ncol, ε_colors, nz, col_idx)
end

"""
    color_groups(wet3d, iwet, nwet, dist)

Builds color groups.
Outputs:
ncol, ε_colors, col_idx
`ncol` is the number of colors (should be (1+2dist)^2)
`ε_colors`` is the matrix of 1.0's used to perturb `x` in all directions of each color.
`col_idx` is the vector of the colors of each box.

Inputs:
`dist` is the infitiny-norm distance of adjacent boxes.
(`dist = 2` required for OCIM1.1 and POP2 and means 125 adjacent boxes.)
"""
function color_groups(wet3d, iwet, nwet, dist)
    nlon = size(wet3d, 2)
    if (nlon % (2dist + 1) ≠ 0) & (nlon ≠ 1)
        error("nlon must be a multiple of 2*dist+1!!")
    end
    col_idx_3d = zeros(Int64, size(wet3d))
    color = 0
    n = 1 + 2dist # dist is the maximum distance of adjacent boxes
    for idepth in 1:n, ilon in 1:n, ilat in 1:n
        color += 1
        # @show (ilat, ilon, idepth)
        col_idx_3d[ilat:n:end, ilon:n:end, idepth:n:end] .= color
    end
    nlon == 1 ? col_idx_3d[:] .= collect(1:nwet) : nothing # hack for 4-box-model
    col_idx = col_idx_3d[iwet]

    ncol = n^3 # number of colors
    ε_colors = zeros(nwet, ncol) # for perturbation in numJac
    for color = 1:ncol
        indices_of_color = findall(col_idx .== color)
        ε_colors[indices_of_color, color] .= 1.0
    end
    ε_colors = ε * ε_colors # make it dual-unit valued

    return ncol, ε_colors, col_idx
end

function build_big_compressor(wet3d, iwet, nwet, dist, nTracers)
    # Decompression indices
    # We use the indices of the adjacency matrix
    AM_3d = adjacency_matrix(wet3d, dist) # Full wet and dry adjacency adjacency_matrix
    AM = AM_3d[iwet, iwet] # Adjacency matrix (wet points only)
    i, j, _ = findnz(AM) # Sub indices

    # Compression indices
    # We use the color indices
    # The compressed index `icmp` is the same as the decompressed `i`
    # The compressed index `jcmp` is the color index
    icmp = i
    ncol, ε_colors, col_idx = color_groups(wet3d, iwet, nwet, dist)
    jcmp = col_idx[j]

    # number of non zeros for storage and use in numJac
    nz = length(i)

    # make everything size nTracers x nTracers
    bigAM = kron(sparse(trues(nTracers, nTracers)), AM)
    bigi, bigj, _ = findnz(bigAM) # Sub indices

    big_col_idx = zeros(eltype(col_idx), nTracers * nwet)
    for t in 0:nTracers-1
        big_col_idx[(1:nwet) .+ t * nwet] .= col_idx .+ t * ncol
    end

    bigicmp = bigi
    bigjcmp = big_col_idx[bigj]
    bigijcmp = CartesianIndex.(bigicmp, bigjcmp)

    bigε_colors = kron(Matrix(I, nTracers, nTracers), ε_colors)

    # Permuted decompression indices (faster!)
    bigncol = ncol * nTracers
    bignwet = nwet * nTracers
    bigv = 1:(nz * nTracers^2) # to track the index permutation that happense during compression
    bigAMv = sparse(bigicmp, bigjcmp, bigv, bignwet, bigncol)
    bigicmp2, bigjcmp2, bigv2 = findnz(bigAMv) # v2 is the verctor of the index permutation
    bigidcmp, bigjdcmp = bigi[bigv2], bigj[bigv2] # This allows direct use of smallJac.nzval
    # which is much faster than accessing those values via (icmp, jcmp)

    return Compressor(bigi, bigj, bigicmp, bigjcmp, bigijcmp, bigidcmp, bigjdcmp, bignwet, bigncol, bigε_colors, nz, big_col_idx)
end


text:
# Author: Ritchie Lee, ritchie.lee@sv.cmu.@schedule
# Date: 12/11/2014

module AbstractGenerativeModelImpl

export AbstractGenerativeModel

abstract AbstractGenerativeModel

end #module

module AbstractGenerativeModelInterfaces

export get,
       isterminal

function get() end

function isterminal() end

end #module

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
using PyPlot, DelimitedFiles, Random, Statistics, Revise,
    HiQGA.transD_GP
## model fixed parts, i.e., air
Random.seed!(23)
zfixed   = [-1e5]
ρfixed   = [1e12]
nmax = 100
# Note that the receiver and transmitter need to be in layer 1
zstart = 0.0
extendfrac, dz = 1.06, 2.
zall, znall, zboundaries = transD_GP.setupz(zstart, extendfrac, dz=dz, n=40, showplot=true)
z, ρ, nfixed = transD_GP.makezρ(zboundaries; zfixed=zfixed, ρfixed=ρfixed)
##  geometry and modeling parameters
rRx = 13.
zRx = -42.0
zTx = -40.0
freqlow = 1e-3
include("electronics_halt.jl")
## LM operator
Flm = transD_GP.AEM_VMD_HMD.HFieldDHT(
                      lowpassfcs = lowpassfcs,
                      times  = LM_times,
                      ramp   = LM_ramp,
                      nmax   = nmax,
                      zTx    = zTx,
                      rRx    = rRx,
                      rTx    = rTx,
                      zRx    = zRx,
                      freqlow = freqlow
                      )
## HM operator
Fhm = transD_GP.AEM_VMD_HMD.HFieldDHT(
                      lowpassfcs = lowpassfcs,
                      times  = HM_times,
                      ramp   = HM_ramp,
                      nmax   = nmax,
                      zTx    = zTx,
                      rRx    = rRx,
                      rTx    = rTx,
                      zRx    = zRx,
                      freqlow = freqlow
                      )
## fill in detail in ohm-m
ρ[(z.>=zstart) .& (z.<50)] .= 20.
ρ[(z.>=50) .&(z.<80)] .= 1
ρ[(z.>=80) .&(z.<100)] .= 20
ρ[(z.>=100) .&(z.<200)] .= 50
ρ[(z.>=200) .&(z.<250)] .= 80
ρ[(z.>=250)]            .= 150
## add noise
transD_GP.plotmodelfield_skytem!(Flm, Fhm, z, ρ)
#remember to specify halt noise if it is in pV
dlow, dhigh, σlow, σhigh = transD_GP.addnoise_skytem(Flm, Fhm,
                z, ρ, noisefrac=0.03, halt_LM = LM_noise*1e-12, halt_HM = HM_noise*1e-12,
                dz=dz, extendfrac=extendfrac, nfixed=nfixed)
## create operator
dlow, dhigh, σlow, σhigh = (dlow, dhigh, σlow, σhigh)./transD_GP.SkyTEM1DInversion.μ₀
aem = transD_GP.dBzdt(Flm, Fhm, dlow, dhigh, σlow, σhigh, z=z, ρ=ρ, nfixed=nfixed)

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA

max_forks_repo_name: GeoscienceAustralia/HiQGA


text:
using ValueOrientedRiskManagementInsurance
using DataFrames
using Test

include("ECModel.jl")

println("Start ECModel test")

@test round(total.gross.profit_mean, digits = 2) ≈ 293.08
@test round(total.net.profit_mean, digits = 2) ≈ 249.39
@test round(total.gross.eco_cap, digits = 2) ≈ 1096.71
@test round(total.net.eco_cap, digits = 2) ≈ 823.70
@test round(total.gross.rorac, digits = 3) ≈ 0.267
@test round(total.net.rorac, digits = 3) ≈ 0.303
## risk adjusted pricing:
@test round(ins_input_rp[1, :loss_ratio], digits = 4) ≈ 0.7684
@test round(ins_input_rp[2, :loss_ratio], digits = 4) ≈ 0.7158
@test round(ins_input_rp[3, :loss_ratio], digits = 4) ≈ 0.7500
## capital optimization
@test round(rorac_net_oc[i_oc_opt], digits = 4) ≈ 0.2690
## product mix optimization, chosen quote f
@test round(avg_ceded_ofr, digits = 4) ≈ 0.4928
@test round(rorac_net_or_opt, digits = 4) ≈ 0.3376

println("End ECModel test")


text:
struct Clibrary
	handle::Ptr{Cvoid}

	Clibrary(libName::Union{AbstractString, Nothing} = nothing, flags = Libdl.RTLD_LAZY | Libdl.RTLD_DEEPBIND | Libdl.RTLD_LOCAL) = new(Libdl.dlopen(libName === nothing ? _NullCString() : libName, flags))
end

Base.close(lib::Clibrary) = Libdl.dlclose(lib.handle)

# NOTE:  risky hack to trick Julia dlopen into ccalling the C dlopen with NULL (meaning to dlopen current process rather than a library)
struct _NullCString <: AbstractString end
Base.cconvert(::Type{Cstring}, ::_NullCString) = Cstring(C_NULL)


text:
mutable struct Config
    # phenotype config
    input_nodes::Int
    output_nodes::Int
    hidden_nodes::Int
    fully_connected::Bool
    max_weight::Float64
    min_weight::Float64
    max_weight_tau::Float64
    min_weight_tau::Float64
    feedforward::Bool
    nn_activation::Symbol
    weight_stdev::Float64

    # GA config
    pop_size::Int
    max_fitness_threshold::Float64
    prob_addconn::Float64
    prob_addnode::Float64
    prob_mutatebias::Float64
    bias_mutation_power::Float64
    prob_mutate_weight::Float64 # dynamic mutation rate (future release)
    weight_mutation_power::Float64
    prob_togglelink::Float64
    elitism::Bool

    # genotype compatibility
    compatibility_threshold::Float64
    compatibility_change::Float64
    excess_coeficient::Float64
    disjoint_coeficient::Float64
    weight_coeficient::Float64

    # species
    species_size::Int
    survival_threshold::Float64 # only the best 20% for each species is allowed to mate
    old_threshold::Int
    youth_threshold::Int
    old_penalty::Float64    # always in (0,1)
    youth_boost::Float64    # always in (1,2)
    max_stagnation::Int

    function Config(params::Dict{String,String})

        new(
            # phenotype
            parse(Int,params["input_nodes"]),
            parse(Int,params["output_nodes"]),
            parse(Int,params["hidden_nodes"]),
            Bool(parse(Int,params["fully_connected"])),
            parse(Float64,params["max_weight"]),
            parse(Float64,params["min_weight"]),
            parse(Float64,params["max_weight_tau"]),
            parse(Float64,params["min_weight_tau"]),
            Bool(parse(Int,params["feedforward"])),
            Symbol(params["nn_activation"]),
            parse(Float64,params["weight_stdev"]),

            # GA
            parse(Int,params["pop_size"]),
            parse(Float64,params["max_fitness_threshold"]),
            parse(Float64,params["prob_addconn"]),
            parse(Float64,params["prob_addnode"]),
            parse(Float64,params["prob_mutatebias"]),
            parse(Float64,params["bias_mutation_power"]),
            parse(Float64,params["prob_mutate_weight"]),
            parse(Float64,params["weight_mutation_power"]),
            parse(Float64,params["prob_togglelink"]),
            Bool(parse(Int,params["elitism"])),

            # genotype compatibility
            parse(Float64,params["compatibility_threshold"]),
            parse(Float64,params["compatibility_change"]),
            parse(Float64,params["excess_coeficient"]),
            parse(Float64,params["disjoint_coeficient"]),
            parse(Float64,params["weight_coeficient"]),

            # species
            parse(Int,params["species_size"]),
            parse(Float64,params["survival_threshold"]),
            parse(Int,params["old_threshold"]),
            parse(Int,params["youth_threshold"]),
            parse(Float64,params["old_penalty"]),
            parse(Float64,params["youth_boost"]),
            parse(Int,params["max_stagnation"])
         )
    end
end

function loadConfig(file::String)
    str = open(file) do f
        read(f, String)
    end

    str = replace(str, r"\r(\n)?" => '\n')
    ls = Base.split(str, "\n")

    ls = filter(l->length(l)>0 && l[1] != '#', ls)
    lsMap = map(x->Base.split(x,'='),ls)
    params = Dict{String,String}()
    for i = 1:length(lsMap)
#         println((i,lsMap[i][1],lsMap[i][2]))
        params[rstrip(lsMap[i][1])] = lstrip(lsMap[i][2])
    end

    return params
end


text:
# ---------------------------------------------------------------------------
# Construction gloutonne d'une solution admissible de SPP - porter de DM1

function GreedyConstruction(C, A)
    x = zeros(Int64, length(C))
    S = 1:length(C)
    M = 1:size(A, 1)
    E = Array{Array{Int64,1},1}(undef, length(C))
    Einit = Array{Array{Int64,1},1}(undef, length(C))
    F = Array{Array{Int64,1},1}(undef, length(C))
    Cout = Array{Float64,1}(undef, length(C))
    firstloop = true
    while (length(S) != 0)
        fill!(E, Int64[])
        fill!(F, Int64[])
        fill!(Cout, -Inf)
        if length(S) == 1
            imax = S[1]
        else
            for j in M
                current = intersect(findall(isequal(1), A[j, :]), S)
                for i in current
                    E[i] = union(E[i], setdiff(current, i))
                    F[i] = union(F[i], j)
                end
            end
            for i in S
                if length(E[i]) != 0
                    Cout[i] = C[i] - sum(C[u] for u in E[i])
                end
            end
            if maximum(Cout) == -Inf        #pas de conflit
                imax = S[1]
            else
                imax = findfirst(isequal(maximum(Cout)), Cout)
            end
        end
        if firstloop
            Einit = copy(E)
            firstloop = false
        end
        x[imax] = 1
        S = setdiff(S, union(E[imax], imax))
        M = setdiff(M, F[imax])
    end
    z = sum(C[i]x[i] for i in 1:length(C))
    return x, z, Einit
end

# ---------------------------------------------------------------------------
# le metaheuristique Tabou

function Tabou(C,A,iterMax,iterPenaliser,length_tm,xInit,zInit,Einit)
    iter  = 1
    nIter = 1
    x = copy(xInit)
    z = zInit
    xBest = copy(x)
    zBest = z
    tb_list = Array{Tuple{Int64,Int64},1}(undef, length_tm)
    fill!(tb_list,(0,0))
    memory_longtime = zeros(Int64,2,length(C))          #1ere ligne est 0, 2eme ligne est 1
    while iter < iterMax
        eqZero = findall(isequal(0),x)
        eqOne  = findall(isequal(1),x)
        for j in eqOne
            setdiff!(eqZero,Einit[j])
        end
        changeable = union(eqZero,eqOne)    #table d'indice de x[j] qui peut être fait add ou drop (neighbor search)
        neighborBest = (0,-Inf)
        foundAdd = false                    #var boolean indiqué si on a trouvé un voisin valid avec mouvement add

        #local search
        for j in changeable
            if nIter < iterPenaliser        #pas de penalisation
                if x[j] == 0                #mouvement add
                    zPrime = z + C[j]
                    if (zPrime > zBest) || (!((j,1) in tb_list) && zPrime > neighborBest[2])           #override tabou or not tabou
                        neighborBest = (j,zPrime)
                        foundAdd = true
                    end
                elseif !foundAdd            #pas besoin de faire le mouvement drop si on a deja trouvé un voisin valid avec mouvement add (en cas de sans penaliser)
                    zPrime = z - C[j]
                    if zPrime > neighborBest[2] && !((j,0) in tb_list)                                #not tabou
                        neighborBest = (j,zPrime)
                    elseif length(changeable) == 1    #cas special, une seule variable x[j] = 1 qui a conflit avec tous les autres mais ne peut pas faire mouvement drop car tabou
                        neighborBest = (j,zPrime)   #override tabou
                    end
                end
            else                            #avec la penalisation pour chercher sur un autre voisinage
                if x[j] == 0
                    zPrime = z + C[j] - memory_longtime[1,j]
                    if (zPrime + memory_longtime[1,j] > zBest)                                         #override tabou and memory longtime
                        neighborBest = (j,zPrime + memory_longtime[1,j])
                    elseif (zPrime > neighborBest[2] && !((j,1) in tb_list))                           #not tabou
                        neighborBest = (j,zPrime)
                    end
                else
                    zPrime = z - C[j] - memory_longtime[2,j]
                    if zPrime > neighborBest[2] && !((j,0) in tb_list)                                #not tabou
                        neighborBest = (j,zPrime)
                    elseif length(changeable) == 1    #cas special, une seule variable x[j] = 1 qui a conflit avec tous les autres mais ne peut pas faire mouvement drop car tabou
                        neighborBest = (j,zPrime)   #override tabou
                    end
                end
            end
        end

        #modifier x, z, taboulist, memoire longtemps, xBest, zBest, nIter, iter
        if neighborBest[1] == 0             #le cas de blocage car |TM| trop grande
            println("Blocked ! Every mouvement is TABOU ! Reduce |TM| by 1 ...")
            xBest2, zBest2, length_tm2 = Tabou(C, A, iterMax, iterPenaliser, length_tm - 1, xInit, zInit, Einit)
            if zBest2 > zBest       #test pour choisir la meilleure solution entre avant et apres modifier |TM|
                return xBest2, zBest2, length_tm2
            else
                return xBest, zBest, length_tm
            end
        else
            before = x[neighborBest[1]]
            after  = (before == 0 ? 1 : 0)
            if (nIter >= iterPenaliser)
                if (neighborBest[2] < zBest)        #le statut tabou n'est pas dominé
                    neighborBest = (neighborBest[1], neighborBest[2] + memory_longtime[before+1,neighborBest[1]])        #retouner à la valeur sans penaliser
                end
                nIter = 0       #reinitialiser nIter apres changer de voisinage
            end
            x[neighborBest[1]] = after
            z = neighborBest[2]
            tb_list = push!(tb_list[2:length_tm],(neighborBest[1],before))
            memory_longtime[before+1,neighborBest[1]] += 1
            if z > zBest
                zBest = z
                xBest = copy(x)
                nIter = 1       #reinitialiser nIter si trouver la nouvelle meilleure solution
                println("New zBest = ",z," found at iter = ",iter," !")
            else
                nIter += 1
            end
            iter += 1
        end
    end
    println("|TM| final = ", length_tm)
    return xBest, zBest, length_tm
end

max_stars_repo_name: ducanh-le/metaSPPstu

max_issues_repo_name: ducanh-le/metaSPPstu

max_forks_repo_name: ducanh-le/metaSPPstu


text:
"""
    servedocs_callback!(docwatcher, filepath, path2makejl, literate)

Custom callback used in [`servedocs`](@ref) triggered when the file corresponding to `filepath`
is changed. If that file is `docs/make.jl`, the callback will check whether any new files have
been added in the `docs/src` folder and add them to the watched files, it will also remove any
file that may have been deleted or renamed.
Otherwise, if the modified file is in `docs/src` or is `docs/make.jl`, a pass of Documenter is
triggered to regenerate the documents, subsequently the LiveServer will render the produced pages
in `docs/build`.
"""
function servedocs_callback!(dw::SimpleWatcher, fp::AbstractString, makejl::AbstractString,
                             literate::String="")
    # if the file that was changed is the `make.jl` file, assume that maybe new files are # referenced and so refresh the vector of watched files as a result.
    if fp == makejl
        # it's easier to start from scratch (takes negligible time)
        empty!(dw.watchedfiles)
        scan_docs!(dw, literate)
    end
    fext = splitext(fp)[2]
    P1 = fext ∈ (".md", ".jl")
    # the second condition is for CSS files, we want to track it but not the output
    # if we track the output then there's an infinite loop being triggered (see docstring)
    if P1 || (fext == ".css" && !occursin(joinpath("docs", "build", "assets"), fp))
        Main.include(makejl)
        file_changed_callback(fp)
    end
    return nothing
end

"""
    scan_docs!(dw::SimpleWatcher, literate="")

Scans the `docs/` folder in order to recover the path to all files that have to be watched and add
those files to `dw.watchedfiles`. The function returns the path to `docs/make.jl`. A list of
folders and file paths can also be given for files that should be watched in addition to the
content of `docs/src`.
"""
function scan_docs!(dw::SimpleWatcher, literate::String="")
    src = joinpath("docs", "src")
    if !(isdir("docs") && isdir(src))
        @error "I didn't find a docs/ or docs/src/ folder."
    end
    makejl = joinpath("docs", "make.jl")
    push!(dw.watchedfiles, WatchedFile(makejl))
    if isdir("docs")
        # add all files in `docs/src` to watched files
        for (root, _, files) ∈ walkdir(joinpath("docs", "src")), file ∈ files
            push!(dw.watchedfiles, WatchedFile(joinpath(root, file)))
        end
    end
    if !isempty(literate)
        isdir(literate) || @error "I didn't find the provided literate folder $literate."
        for (root, _, files) ∈ walkdir(literate), file ∈ files
            push!(dw.watchedfiles, WatchedFile(joinpath(root, file)))
        end
    end

    # When using literate.jl, we should only watch the source file otherwise we would double
    # trigger: first when the script.jl is modified then again when the script.md is created
    # which would cause an infinite loop if both `script.jl` and `script.md` are watched.
    # So here we remove from the watchlist all files.md that have a files.jl with the same path.
    remove = Int[]
    if isempty(literate)
        # assumption is that the scripts are in `docs/src/...` and that the generated markdown
        # goes in exactly the same spot so for instance:
        # docs
        # └── src
        #     ├── index.jl
        #     └── index.md
        for wf ∈ dw.watchedfiles
            spath = splitext(wf.path)
            spath[2] == ".jl" || continue
            k = findfirst(e -> splitext(e.path) == (spath[1], ".md"), dw.watchedfiles)
            k === nothing || push!(remove, k)
        end
    else
        # assumption is that the scripts are in `literate/` and that the generated markdown goes
        # in `docs/src` with the same relative paths so for instance:
        # docs
        # ├── lit
        # │   └── index.jl
        # └── src
        #     └── index.md
        for (root, _, files) ∈ walkdir(literate), file ∈ files
            spath = splitext(joinpath(root, file))
            spath[2] == ".jl" || continue
            path = replace(spath[1], Regex("^$literate") => joinpath("docs", "src"))
            k = findfirst(e -> splitext(e.path) == (path, ".md"), dw.watchedfiles)
            k === nothing || push!(remove, k)
        end
    end
    deleteat!(dw.watchedfiles, remove)
    return makejl
end

"""
    servedocs(; verbose=false, literate="")

Can be used when developing a package to run the `docs/make.jl` file from Documenter.jl and
then serve the `docs/build` folder with LiveServer.jl. This function assumes you are in the
directory `[MyPackage].jl` with a subfolder `docs`.

* `verbose` is a boolean switch to make the server print information about file changes and
connections.
* `literate` is the path to the folder containing the literate scripts, if left empty, it will be
assumed that they are in `docs/src`.
"""
function servedocs(; verbose::Bool=false, literate::String="")
    # Custom file watcher: it's the standard `SimpleWatcher` but with a custom callback.
    docwatcher = SimpleWatcher()
    set_callback!(docwatcher, fp->servedocs_callback!(docwatcher, fp, makejl, literate))

    # Retrieve files to watch
    makejl = scan_docs!(docwatcher, literate)

    # trigger a first pass of Documenter (& possibly Literate)
    Main.include(makejl)

    # note the `docs/build` exists here given that if we're here it means the documenter
    # pass did not error and therefore that a docs/build has been generated.
    serve(docwatcher, dir=joinpath("docs", "build"), verbose=verbose)
    return nothing
end

#
# Miscellaneous utils
#

"""
    setverbose(b)

Set the verbosity of LiveServer to either `true` (showing messages upon events) or `false` (default).
"""
setverbose(b::Bool) = (VERBOSE.x = b)

"""
    example()

Simple function to copy an example website folder to the current working directory that can be
watched by the LiveServer to get an idea of how things work.

### Example

```julia
LiveServer.example()
cd("example")
serve()
```
"""
function example(; basedir="")
    isempty(basedir) && (basedir = pwd())
    cp(joinpath(dirname(dirname(pathof(LiveServer))), "example"), joinpath(basedir, "example"))
end


text:
using Gadfly

set_default_plot_size(8inch, 3inch)

dodge1 = plot(x=["a","a"], y=[3,2], color=["red","blue"],
      Geom.bar(position=:dodge));
dodge2 = plot(x=["a","a","b","b"], y=[3,2,1,2], color=["red","blue","red","blue"],
      Geom.bar(position=:dodge));
stack1 = plot(x=["a","a"], y=[3,2], color=["red","blue"],
      Geom.bar(position=:stack));
stack2 = plot(x=["a","a","b","b"], y=[3,2,1,2], color=["red","blue","red","blue"],
      Geom.bar(position=:stack));

gridstack([dodge1 dodge2; stack1 stack2])


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Data: 06/04/2014

module Util

include("format_string.jl")

export outputGFormatString

end


text:
module PlusMinus

export plusminus, ±, minusplus, ∓

plusminus(x::Number,y::Number) = (x+y,x-y)
±(x,y) = plusminus(x,y)

minusplus(x::Number,y::Number) = (x-y,x+y)
∓(x,y) = minusplus(x,y)

end # module


text:
module ExcludeFile

include_dependency(joinpath(dirname(@__DIR__), "deps", "dependency.txt"))
include("f.jl")

end # module


text:
@testset "Bitstring" begin

    @testset for f in (-64.0,-1.0,-1e-2,1e-1,1.0,64.0)
        for T in (Posit8,Posit8_2)
            T(f) == T(parse(UInt8,bitstring(T(f)),base=2))
        end

        for T in (Posit16,Posit16_2)
            T(f) == T(parse(UInt16,bitstring(T(f)),base=2))
        end

        # don't test for Posit24_2 as there is no 24bit UInt format.
        # rely on bitstring(p," ") testing

        Posit32(f) == Posit32(parse(UInt32,bitstring(Posit32(f)),base=2))
    end

    # Posit8
    @test "0 10 00000" == bitstring(Posit8(1)," ")
    @test "0 111110 1" == bitstring(Posit8(24)," ")
    @test "0 1111110" == bitstring(Posit8(32)," ")
    @test "0 1111111" == bitstring(Posit8(64)," ")

    # Posit16
    @test "0 111111111111110" == bitstring(Posit16(1e8)," ")
    @test "0 11111111111110 1" == bitstring(Posit16(5e7)," ")
    @test "0 1111111111110 1 0" == bitstring(Posit16(1e7)," ")
    @test "0 10 0 000000000000" == bitstring(Posit16(1)," ")

    # Posit32
    @test "0 1111111111111111111111111111110" == bitstring(Posit32(1e35)," ")
    @test "0 111111111111111111111111111110 0" == bitstring(Posit32(1e34)," ")
    @test "0 11111111111111111111111111110 10" == bitstring(Posit32(1e33)," ")
    @test "0 1111111111111111111111111110 10 0" == bitstring(Posit32(1e32)," ")
    @test "0 10 00 000000000000000000000000000" == bitstring(Posit32(1)," ")

    # Posit8_2
    @test "0 10 00 000" == bitstring(Posit8_2(1)," ")
    @test "0 11110 01" == bitstring(Posit8_2(1e4)," ")
    @test "0 111110 0" == bitstring(Posit8_2(1e5)," ")
    @test "0 1111110" == bitstring(Posit8_2(1e6)," ")

    # Posit16_2
    @test "0 10 00 00000000000" == bitstring(Posit16_2(1)," ")
    @test "0 111111111111110" == bitstring(Posit16_2(1e16)," ")
    @test "0 11111111111110 1" == bitstring(Posit16_2(1e15)," ")
    @test "0 1111111111110 10" == bitstring(Posit16_2(1e14)," ")

    # Posit24_2
    @test "0 11111111111111111111110" == bitstring(Posit24_2(1e25)," ")
    @test "0 1111111111111111111110 0" == bitstring(Posit24_2(1e24)," ")
    @test "0 111111111111111111110 00" == bitstring(Posit24_2(1e23)," ")
    @test "0 10 00 0000000000000000000" == bitstring(Posit24_2(1)," ")

end


text:
# test adding new infinite domain
@testset "Infinite Domains" begin
    # load in the extension
    include("./extensions/infinite_domain.jl")

    # initialize model
    m = InfiniteModel()

    # test definition
    @test MyNewDomain(0, 1) isa MyNewDomain
    domain = MyNewDomain(0, 1)
    @test @infinite_parameter(m, domain = domain) isa GeneralVariableRef
    @test infinite_domain(first(all_parameters(m))) == domain
    @test @infinite_parameter(m, par in domain, supports = [0, 1]) isa GeneralVariableRef
    infinite_domain(m[:par]) == domain
    supports(m[:par]) == [0., 1.]
    @test @infinite_parameter(m, par2 in domain, num_supports = 3) isa GeneralVariableRef
    supports(m[:par2]) == [0., 0.5, 1.]
    @test @infinite_parameter(m, [1:2], domain = domain, num_supports = 3) isa Vector

    # set support methods
    par = m[:par]
    @test_throws ErrorException add_supports(par, 2)

    # test bound methods
    @test has_lower_bound(m[:par2])
    @test lower_bound(m[:par2]) == 0
    @test set_lower_bound(m[:par2], -1) isa Nothing
    @test lower_bound(m[:par2]) == -1
    @test has_upper_bound(m[:par2])
    @test upper_bound(m[:par2]) == 1
    @test set_upper_bound(m[:par2], 0) isa Nothing
    @test upper_bound(m[:par2]) == 0

    # add variables
    @test @variable(m, x >= 0, Infinite(par)) isa GeneralVariableRef
    x = m[:x]

    # test measures
    @test @integral(x^2 + par, par, num_nodes = 2, eval_method = GaussLegendre()) isa GeneralVariableRef
    @test expect(x, par) isa GeneralVariableRef

    # test constraints
    @test @constraint(m, x + par <= 0) isa InfOptConstraintRef

    # transcribe the model
    @test build_optimizer_model!(m) isa Nothing
    @test num_variables(optimizer_model(m)) == 4
end

# test using the new derivative evaluation method
@testset "Derivative Method" begin
    # load in the extension
    include("./extensions/derivative_method.jl")

    # initialize model and objects
    m = InfiniteModel()
    @infinite_parameter(m, t in [0, 10], num_supports = 3)
    @variable(m, x, Infinite(t))
    d = deriv(x, t)
    method = MyDerivMethod(0.5)

    # test definition and addition to a parameter
    @test method isa MyDerivMethod
    @test set_derivative_method(t, method) isa Nothing
    @test derivative_method(t) isa MyDerivMethod

    # test preliminaries to evaluation
    @test generative_support_info(method) == UniformGenerativeInfo([0.5], InternalLabel)
    @test support_label(method) == InternalLabel

    # test evaluation
    @test InfiniteOpt.evaluate_derivative(d, method, m) isa Vector
    @test has_generative_supports(t)
    @test num_supports(t, label = All) == 5
    @test evaluate(d) isa Nothing
    @test length(derivative_constraints(d)) == 4
    @test num_constraints(m) == 4
end

# Test extensions of measure data
@testset "Measure Data" begin
    # load in the extension
    include("./extensions/measure_data.jl")

    # setup the model
    m = InfiniteModel()
    @infinite_parameter(m, t in [0, 10])
    @variable(m, x >= 0, Infinite(t))
    @variable(m, z)
    data1 = DiscreteMeasureData(t, [2.5, 2.5], [2., 4.], lower_bound = 0., upper_bound = 5.)
    data2 = DiscreteMeasureData(t, [1., 1.], [3., 6.], lower_bound = 1., upper_bound = 10.)

    # test definition
    @test NewMeasureData("test", data1) isa NewMeasureData
    new_data1 = NewMeasureData("test", data1)
    new_data2 = NewMeasureData("test", data2)

    # test data queries
    @test parameter_refs(new_data1) == t
    @test parameter_refs(new_data2) == t
    @test supports(new_data1) == supports(data1)
    @test supports(new_data2) == supports(data2)
    @test coefficients(new_data1) == coefficients(data1)
    @test coefficients(new_data2) == coefficients(data2)
    @test weight_function(new_data1) == weight_function(data1)
    @test weight_function(new_data2) == weight_function(data2)

    # test measure definition
    @test measure(x + z, new_data1) isa GeneralVariableRef
    @test new_measure(x^2, t, 0, 4, num_supports = 4) isa GeneralVariableRef
    @test_throws ErrorException new_measure(x^2 + z, t, 6, 10)

    # test expansion
    pvrefs = [GeneralVariableRef(m, i, PointVariableIndex) for i in 1:2]
    @test expand(measure(x + z, new_data1)) == 2.5 * (pvrefs[1] + pvrefs[2]) + 5z

    # test transcription
    @test @constraint(m, z == measure(x, new_data1)) isa InfOptConstraintRef
    @test build_optimizer_model!(m) isa Nothing
    @test num_variables(optimizer_model(m)) == 6

    # test deletion
    @test_throws ErrorException delete(m, t)
    indices = collect(keys(m.measures))
    for index in indices
        @test delete(m, MeasureRef(m, index)) isa Nothing
    end
    @test delete(m, x) isa Nothing
    @test delete(m, t) isa Nothing
end

# Test extensions of measure evaluation methods
@testset "Measure Evaluation" begin
    # load in the extension
    include("./extensions/measure_eval.jl")

    # set up the model
    m = InfiniteModel()
    @infinite_parameter(m, t in [0, 5])
    @infinite_parameter(m, x[1:2] in [0, 1], independent = true)
    @infinite_parameter(m, p[1:2] in [0, 1])
    @infinite_parameter(m, xi ~ Normal(0., 1.))
    @variable(m, y >= 0, Infinite(t))
    @variable(m, f, Infinite(x))
    mref = integral(y^2 + t, t, 0, 4, num_supports = 5,
                    eval_method = NewUniEvalMethod())
    @test supports(measure_data(mref)) == Array([0., 1., 2., 3., 4.])
    warn = "The method is implemented for independent multivariate parameters."
    @test_logs (:warn, warn) integral(f, p, num_supports = 3,
                                      eval_method = NewMultiEvalMethod())
    mref2 = integral(f, x, num_supports = 3, eval_method = NewMultiEvalMethod())
    @test supports(measure_data(mref2)) == Float64[0 0.5 1; 0 0.5 1]
end

# Test extensions of generative support info
@testset "Generative Support Info" begin
    # load in the extension
    include("./extensions/generative_info.jl")

    # set up the model
    m = InfiniteModel()
    @infinite_parameter(m, t in [0, 5], num_supports = 4)
    pref = dispatch_variable_ref(t)

    # test creation and basic functions
    @test MyGenInfo(2) isa AbstractGenerativeInfo
    @test support_label(MyGenInfo(2)) == MyGenLabel
    @test length(make_generative_supports(MyGenInfo(2), t, supports(t))) == 6

    # test incorporation into infinite parameter
    @test InfiniteOpt._set_generative_support_info(pref, MyGenInfo(2)) isa Nothing
    @test add_generative_supports(t) isa Nothing
    @test num_supports(t, label = All) == 10
end

# Test otpimizer model extensions
@testset "Optimizer Model" begin
    # load in the extension
    include("./extensions/optimizer_model.jl")

    # setup the infinite model
    m = InfiniteModel()
    @infinite_parameter(m, par in [0, 1])
    @variable(m, x, Infinite(par))
    @variable(m, x0, Point(x, 0))
    @variable(m, y)
    meas = integral(x, par)
    @constraint(m, c1, x + y - 2 <= 0)
    @constraint(m, c2, meas == 0)
    @constraint(m, c3, x0 + y == 5)
    @objective(m, Min, y)

    # test optimizer model constructor
    @test NewReformModel() isa Model
    @test NewReformModel().ext[:ReformData] isa NewReformData

    # test data extraction
    @test_throws ErrorException reform_data(optimizer_model(m))
    @test reform_data(NewReformModel()) isa NewReformData

    # test replacing the optimizer model
    @test set_optimizer_model(m, NewReformModel()) isa Nothing
    @test haskey(optimizer_model(m).ext, :ReformData)

    # test making InfiniteModel with the new optimizer model
    @test InfiniteModel(OptimizerModel = NewReformModel) isa InfiniteModel
    @test optimizer_model_key(InfiniteModel(OptimizerModel = NewReformModel)) == :ReformData

    # test retrival errors
    @test_throws ErrorException optimizer_model_variable(x)
    @test_throws ErrorException optimizer_model_variable(x0)
    @test_throws ErrorException optimizer_model_variable(y)
    @test_throws ErrorException optimizer_model_constraint(c1)
    @test_throws ErrorException optimizer_model_constraint(c2)
    @test_throws ErrorException optimizer_model_constraint(c3)
    @test_throws ErrorException supports(x)
    @test_throws ErrorException supports(c1)
    @test_throws ErrorException supports(c2)
    @test_throws ErrorException supports(c3)

    # test build_optimizer_model!
    @test build_optimizer_model!(m, my_kwarg = true) isa Nothing
    @test haskey(optimizer_model(m).ext, :ReformData)
    @test num_variables(optimizer_model(m)) == 13

    # test retrivals
    @test optimizer_model_variable(x, my_kwarg = true) isa Vector{VariableRef}
    @test optimizer_model_variable(x0, my_kwarg = true) isa VariableRef
    @test optimizer_model_variable(y, my_kwarg = true) isa VariableRef
    @test optimizer_model_variable(meas, my_kwarg = true) isa Vector{VariableRef}
    @test optimizer_model_constraint(c1, my_kwarg = true) isa Vector{<:ConstraintRef}
    @test optimizer_model_constraint(c2, my_kwarg = true) isa Vector{<:ConstraintRef}
    @test optimizer_model_constraint(c3, my_kwarg = true) isa Vector{<:ConstraintRef}
    @test optimizer_model_expression(x^2) == zero(AffExpr)
    @test supports(x, my_kwarg = true) == [(0.,), (1.,)]
    @test supports(y) == ()
    @test supports(meas) == [(-1.,), (-2.,)]
    @test supports(c1, my_kwarg = true) == [(2.,), (3.,)]
    @test supports(x + y) == [(-42.,), (1.,)]

    # test optimization with rebuild
    mockoptimizer = () -> MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()),
                                             eval_objective_value=false)
    @test set_optimizer_model_ready(m, false) isa Nothing
    @test set_optimizer(m, mockoptimizer) isa Nothing
    @test set_silent(m)
    @test set_time_limit_sec(m, 42) == 42
    @test optimize!(m) isa Nothing
    @test get_optimizer_attribute(m, MOI.Silent())

    # prepare for result queries
    mockoptimizer = JuMP.backend(m).optimizer.model
    MOI.set(mockoptimizer, MOI.TerminationStatus(), MOI.OPTIMAL)
    MOI.set(mockoptimizer, MOI.RawStatusString(), "solver specific string")
    MOI.set(mockoptimizer, MOI.ResultCount(), 2)
    MOI.set(mockoptimizer, MOI.PrimalStatus(1), MOI.FEASIBLE_POINT)
    MOI.set(mockoptimizer, MOI.DualStatus(1), MOI.FEASIBLE_POINT)
    MOI.set(mockoptimizer, MOI.SolveTime(), 0.42)
    MOI.set(mockoptimizer, MOI.ObjectiveValue(), -1.0)
    MOI.set(mockoptimizer, MOI.ObjectiveBound(), 2.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(),
            JuMP.optimizer_index(optimizer_model_variable(x)[1]), -1.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(),
            JuMP.optimizer_index(optimizer_model_variable(x)[2]), 0.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(),
            JuMP.optimizer_index(optimizer_model_variable(y)), 1.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(),
            JuMP.optimizer_index(optimizer_model_variable(x0)), 42.)
    MOI.set(mockoptimizer, MOI.VariablePrimal(),
            JuMP.optimizer_index(optimizer_model_variable(meas)[1]), 2.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(),
            JuMP.optimizer_index(optimizer_model_variable(meas)[2]), -2.0)
    MOI.set(mockoptimizer, MOI.ConstraintDual(),
            JuMP.optimizer_index(optimizer_model_constraint(c1)[1]), -1.0)
    MOI.set(mockoptimizer, MOI.ConstraintDual(),
            JuMP.optimizer_index(optimizer_model_constraint(c1)[2]), -1.0)
    MOI.set(mockoptimizer, MOI.ConstraintDual(),
            JuMP.optimizer_index(optimizer_model_constraint(c2)[1]), 0.0)
    MOI.set(mockoptimizer, MOI.ConstraintDual(),
            JuMP.optimizer_index(optimizer_model_constraint(c2)[2]), -1.0)
    MOI.set(mockoptimizer, MOI.ConstraintDual(),
            JuMP.optimizer_index(optimizer_model_constraint(c3)[1]), 0.0)
    MOI.set(mockoptimizer, MOI.ConstraintDual(),
            JuMP.optimizer_index(optimizer_model_constraint(c3)[2]), -1.0)

    # test result queries
    @test termination_status(m) == MOI.OPTIMAL
    @test raw_status(m) == "solver specific string"
    @test result_count(m) == 2
    @test primal_status(m) == MOI.FEASIBLE_POINT
    @test dual_status(m) == MOI.FEASIBLE_POINT
    @test solve_time(m) == 0.42
    @test has_values(m)
    @test has_duals(m)
    @test objective_bound(m) == 2
    @test objective_value(m) == -1
    @test value(x) == [-1, 0]
    @test value(x0) == 42.
    @test value(y) == 1
    @test value(meas) == [2., -2.]
    @test value(x + y) == 0.
    @test dual(c1) == [-1, -1]
    @test dual(c2) == [0., -1.]
    @test dual(c3) == [0., -1.]
    @test optimizer_index(x) == optimizer_index.(optimizer_model_variable(x))
    @test optimizer_index(x0) == optimizer_index(optimizer_model_variable(x0))
    @test optimizer_index(y) == optimizer_index(optimizer_model_variable(y))
    @test optimizer_index(c1) == optimizer_index.(optimizer_model_constraint(c1))
    @test optimizer_index(c2) == optimizer_index.(optimizer_model_constraint(c2))
    @test optimizer_index(c3) == optimizer_index.(optimizer_model_constraint(c3))
    @test shadow_price(c1) == [1, 1]
    @test shadow_price(c2) == [-0., 1.]
    @test shadow_price(c3) == [-0., 1.]
end


text:
"""
Retrieve the composition of a user-defined portfolio as a data frame.
The name of the portfolio is something like U12345678-1 (see PRTU <GO>).
Four kinds of requests exist: members (default), positions, weights, data.
Use the date argument to retrieve historical compositions (YYYYMMDD format).

#Examples
```julia
bport("U12345678-1")
bport("U12345678-1",retrieve="data")
bport("U12345678-1",date="20170930")
```
"""
function bport(name::String;retrieve::String="members",date::String="",overrides::Dict=Dict())
    global session
    global services
    field = @Match.match retrieve begin
        "members" => "PORTFOLIO_MEMBERS"
        "positions" => "PORTFOLIO_MPOSITION"
        "weights" => "PORTFOLIO_MWEIGHT"
        "data" => "PORTFOLIO_DATA"
        _ => error("the retrieve parameter has to be one of [members, positions, weights, data]")
    end
    service = services["refdata"]
    request = @with_pointer request blpapi_Service_createRequest(service,request,"PortfolioDataRequest")
    elements = blpapi_Request_elements(request)
    try
        securities_ = @with_pointer securities_ blpapi_Element_getElement(elements,securities_,"securities",C_NULL)
        blpapi_Element_setValueString(securities_,name*" Client",-1)
        fields_ = @with_pointer fields_ blpapi_Element_getElement(elements,fields_,"fields",C_NULL)
        blpapi_Element_setValueString(fields_,field,-1)
        if length(date) > 0 || length(overrides) > 0
            overrides_ = @with_pointer overrides_ blpapi_Element_getElement(elements,overrides_,"overrides",C_NULL)
            new_ = @with_pointer new_ blpapi_Element_appendElement(overrides_,new_)
            if length(date) > 0
                blpapi_Element_setElementString(new_, "fieldId", C_NULL, "REFERENCE_DATE")
                blpapi_Element_setElementString(new_, "value", C_NULL, date)
            end
            for (k,v) in overrides
                new_ = @with_pointer new_ blpapi_Element_appendElement(overrides_,new_)
                blpapi_Element_setElementString(new_, "fieldId", C_NULL, string(k))
                blpapi_Element_setElementString(new_, "value", C_NULL, string(v))
            end
        end
        correlation_id = blpapi_CorrelationId(0,0)
        blpapi_Session_sendRequest(session,request,Ref(correlation_id),C_NULL,C_NULL,C_NULL,0)
    finally
        blpapi_Request_destroy(request)
    end
    process_tree(retrieve_response())[1][2]
end


text:
module Relations

using Compat
using Reexport
@reexport using NullableArrays
@reexport using StructuredQueries
const SQ = StructuredQueries
@reexport using CSV

export  Attributed,
        Relation

include("attributed.jl")
include("relation.jl")
include("utils.jl")
include("show.jl")

include("collect/attributive/select.jl")
include("collect/attributive/filter.jl")
include("collect/attributive/groupby.jl")
include("collect/attributive/summarize.jl")

include("collect/attributive/grouped/summarize.jl")

end # module


text:
#Contains functions to update the networks

#Simply iterate trough the dataset and update:
function do_td3_update(Qf_1::Any, Qf_tar_1::Any, Qf2::Any, Qf_tar_2::Any, Pl::Any, Pl_tar::Any,
    buffer::Traj_Data,params::Params, opt_para_pl, opt_para_qf1, opt_para_qf2)

    n_traj = buffer.n_traj
    shuffle_arr = shuffle(1:n_traj)
    τ = params.τ

    #Parameter selecting when to update the target and policy networks:
    d = params.d_update

    for i_true = 2:n_traj
        i_traj = shuffle_arr[i_true]
        traj = buffer.traj[i_traj]

        #Only every d-steps: Delayed Policy Update
        if mod(i_true, d) == 0
            pl_update = true
        else
            pl_update = false
        end

        #Get the gradients:
        ĝ_qf1, ĝ_qf2, ĝ_pl = grad_all_networks(pl_update, traj, Qf_1, Qf_tar_1, Qf_2, Qf_tar_2, Pl, Pl_tar, params)

        #Always update the normal Q-functions
        Knet.update!(Qf_1, ĝ_qf1, opt_para_qf1)
        Knet.update!(Qf_2, ĝ_qf2, opt_para_qf2)

        #Do the delayed Policy update:
        if pl_update == true
            #update the policy function
            Knet.update!(Pl, ĝ_pl, opt_para_pl)

            #Update all Target functions
            Qf_tar_1 = τ .* Qf_1 .+ (1-τ) .* Qf_tar_1
            Qf_tar_2 = τ .* Qf_2 .+ (1-τ) .* Qf_tar_2
            Pl_tar   = τ .* Pl   .+ (1-τ) .* Pl_tar
        end
    end
end

#This function will update the networks in parallel:
function do_td3_update_parallel!(Qf_1::Any, Qf_tar_1::Any, Qf_2::Any, Qf_tar_2::Any, Pl::Any, Pl_tar::Any,
    buffer::Traj_Data,params::Params, opt_para_pl, opt_para_qf1, opt_para_qf2)

    """
    Setup the problem
    """
    #Take all partial traj and shuffle their order!
    n_traj = buffer.n_traj
    shuffle_arr = shuffle(1:n_traj)
    τ = params.τ

    active_workers = workers()
    n_workers = length(active_workers)

    ĝ_qf1 = copy(Qf_1)
    ĝ_qf2 = copy(Qf_2)
    ĝ_pl  = copy(Pl)

    # Create Channels storing results and starting values
    input_channel  = RemoteChannel(()->Channel{Tuple}(length(workers())+1)) # will contain networks and traj
    output_channel = RemoteChannel(()->Channel{Tuple}(length(workers())+1)) # will contain networks

    # Create problem on every worker
    begin
        for p in workers()
            remote_do(distributed_td3_update, p, params, input_channel, output_channel)
        end
    end

    """
    Do the work in parallel
    """
    partial_step = 0.1
    n_parallel   = Int(ceil((n_traj*partial_step) / n_workers))
    #prog       = Progress(n_parallel,1)

    #Parameter selecting when to update the target and policy networks:
    d = params.d_update

    for i_parallel = 1 : n_parallel

        #Only every d-steps: Delayed Policy Update
        if mod(i_parallel, d) == 0
            pl_update = true
        else
            pl_update = false
        end

        #1. Send the work tasks:
        for i_w = 1 : n_workers
            i_true = ((i_parallel-1)*n_workers + i_w) #get true order
            if i_true <= n_traj
                i_traj  = shuffle_arr[i_true] #shuffle order
                traj    = buffer.traj[i_traj]
                @async put!(input_channel, (pl_update, traj, Qf_1, Qf_tar_1, Qf_2, Qf_tar_2, Pl, Pl_tar))
            end
        end

        # 2. Take the results
        ĝ_qf1 *= 0
        ĝ_qf2 *= 0
        ĝ_pl  *= 0

        for i_w = 1: n_workers
            i_true = ((i_parallel-1)*n_workers + i_w)
            if i_true <= n_traj
                #Take from Channel
                gradients = take!(output_channel)
                #Reorder
                ĝ_qf1 += gradients[1]
                ĝ_qf2 += gradients[2]

                #Only when policy-gradient was calculated
                if pl_update == true
                    ĝ_pl  += gradients[3]
                end

            end
        end

        #Always update the normal Q-functions
        Knet.update!(Qf_1, ĝ_qf1, opt_para_qf1)
        Knet.update!(Qf_2, ĝ_qf2, opt_para_qf2)

        #Only frequently update target networks and Policy
        if pl_update == true
            #update the policy network
            Knet.update!(Pl, ĝ_pl, opt_para_pl)

            #Update all target networks
            Qf_tar_1 = τ .* Qf_1 .+ (1-τ) .* Qf_tar_1
            Qf_tar_2 = τ .* Qf_2 .+ (1-τ) .* Qf_tar_2
            Pl_tar   = τ .* Pl   .+ (1-τ) .* Pl_tar
        end

        #Update Progress
        #ProgressMeter.update!(prog,i_parallel) #update progress bar
    end

    """
    Close the channels
    """
    # Tell workers that its done
    NaN_mat = zeros(params.batch_size)
    fill!(NaN_mat,NaN)
    ending_traj = init_traj(params,params.batch_size)
    ending_traj.r_t = NaN_mat

    for j = 1 : n_workers
        @async put!(input_channel, (true, ending_traj, Qf_1, Qf_tar_1, Qf_2, Qf_tar_2, Pl, Pl_tar))
    end

    #Close Channels
    sleep(0.1)
    close(input_channel)
    close(output_channel)

    return nothing
end

#Function running on the worker
@everywhere function distributed_td3_update(params::Params, input_channel::RemoteChannel, output_channel::RemoteChannel)
   # Start loop where work is done
   while true
       # break loop if NaN state received
       try
           #Take the new Input
           inputs = take!(input_channel)

           #Reorder data:
           pl_update  = inputs[1]
           traj       = inputs[2]
           Qf_1       = inputs[3]
           Qf_tar_1   = inputs[4]
           Qf_2       = inputs[5]
           Qf_tar_2   = inputs[6]
           Pl         = inputs[7]
           Pl_tar     = inputs[8]

           #Check if u reached final
           if (sum(isnan, traj.r_t) > 0)
               GC.gc()
               break
           end

           #DO THE WORK: (in ./td3/gradients)
           #-----------------------------------------------------------------------
           ĝ_qf1, ĝ_qf2, ĝ_pl = grad_all_networks(pl_update, traj, Qf_1, Qf_tar_1, Qf_2, Qf_tar_2, Pl, Pl_tar, params)
           #-----------------------------------------------------------------------

           if pl_update == true
               put!(output_channel, (ĝ_qf1, ĝ_qf2, ĝ_pl))
           else
               put!(output_channel, (ĝ_qf1, ĝ_qf2))
           end

       catch
           println("Some problem occured, worker $(myid()) exiting trajectory optimization function.")
           break
       end
   end

   return nothing
end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
#
# linesearches.jl --
#
# Line search methods for OptimPackNextGen.
#
# -----------------------------------------------------------------------------
#
# This file is part of OptimPackNextGen.jl which is licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2015-2019, Éric Thiébaut.
# <https://github.com/emmt/OptimPackNextGen.jl>.
#

module LineSearches

import ...OptimPackNextGen
import OptimPackNextGen.Float

export
    start!,
    iterate!,
    getreason,
    getstatus,
    getstep,
    gettask,
    usederivatives,
    LineSearch,
    ArmijoLineSearch,
    MoreToraldoLineSearch,
    MoreThuenteLineSearch

const REASONS = Dict{Symbol,String}(
    :NOT_STARTED => "Algorithm not yet started",
    :F_LE_FMIN => "F(X) ≤ FMIN",
    :COMPUTE_FG => "Caller must compute f(x) and g(x)",
    :NEW_X => "A new iterate is available for examination",
    :FINAL_X => "A solution has been found within tolerances",
    :FATOL_TEST_SATISFIED => "FATOL test satisfied",
    :FRTOL_TEST_SATISFIED => "FRTOL test satisfied",
    :GATOL_TEST_SATISFIED => "GATOL test satisfied",
    :GRTOL_TEST_SATISFIED => "GRTOL test satisfied",
    :FIRST_WOLFE_SATISFIED =>"First Wolfe condition satisfied",
    :STRONG_WOLFE_SATISFIED => "Strong Wolfe conditions both satisfied",
    :STP_EQ_STPMIN => "Step at lower bound",
    :STP_EQ_STPMAX => "Step at upper bound",
    :XTOL_TEST_SATISFIED => "XTOL test satisfied",
    :ROUNDING_ERRORS => "Rounding errors prevent progress",
    :NOT_DESCENT => "Search direction is not a descent direction",
    :INIT_STP_LE_ZERO => "Initial step ≤ 0"
)

"""
## Line search methods

Line search methods are instances of types derived from the abstract type
`LineSearch{T}` which is parameterized by the floating point type `T` for the
computations.  Assuming `SomeLineSearch` is a concrete line search type, a
typical line search is performed as follows:

    # Create an instance of the line search method:
    ls = SomeLineSearch(T)

    # Start the line search and loop until a step satisfying
    # some conditions is found:
    x0 = ...            # initial variables
    f0 = func(x)        # function value at x0
    g0 = grad(x)        # gradient at x0
    d = ...             # search direction
    dtg0 = vdot(d, g0)  # directional derivative at x0
    stp = ...           # initial step
    stpmin = ...        # lower bound for the step (usually zero)
    stpmax = ...        # upper bound for the step (usually a large number)
    task = start!(ls, f0, dtg0, stp; stpmin = stpmin, stpmax = stpmax)
    while task == :SEARCH
        stp = getstep(ls) # get step length to try
        x = x0 + stp*d    # compute trial point
        f = func(x)       # function value at x
        g = grad(x)       # gradient at x
        dtg = vdot(d, g)  # directional derivative at x
        task = iterate!(ls, stp, f, dtg)
    end
    task = gettask(ls)
    if task != :CONVERGENCE
        if task == :WARNING
            @warn getreason(ls)
        else
            error(getreason(ls))
        end
    end

Note that the same line search instance may be re-used for subsequent line
searches (with the same settings).

"""
abstract type LineSearch{T<:AbstractFloat} end

"""
    getstep(ls)

yields the length of the next step to try in the line search implemented by
`ls`.
"""
getstep(ls::LineSearch) = ls.stp

"""
    gettask(ls)

yields the current pending task in the line search implemented by `ls`.
"""
gettask(ls::LineSearch) = ls.task

"""
    getstatus(ls)

yields the current status of the line search implemented by `ls`.
"""
getstatus(ls::LineSearch) = ls.status

# Helper function.
function report!(ls::LineSearch, task::Symbol, status::Symbol)
    ls.status = status
    ls.task = task
    return task
end

"""
    getreason(ls)

yields a textual explanation for the current state of the line search
implemented by `ls`.
"""
getreason(ls::LineSearch) = REASONS[ls.status]

"""
    usederivative(ls)

indicates whether the line search instance `ls` requires the derivative of the
objective function when calling `iterate!`.  Alternatively `ls` can also be the
line search type.  Note that the derivative is always needed by the `start!`
method.
"""
usederivatives(::T) where {T<:LineSearch} = usederivatives(T)

"""
    start!(ls, f0, g0, stp; stpmin = 0, stpmax = Inf) -> task

starts a new line search with the method implemented by the line search
instance `ls`, the returned value is the next task to perform (see `iterate!`
for the interpretation of this value).  Arguments `f0` and `g0` are the value
and the directional derivative of the objective function at the start of the
search (that is for a step length equal to zero).  Argument `stp` is the length
of the next step to try and must be nonnegative.  Keywords `stpmin` and
`stpmax` can be used to specify bounds for the length of the step.
"""
function start!(ls::LineSearch{T},
                f0::Real, g0::Real, stp::Real;
                stpmin::Real = zero(T),
                stpmax::Real = typemax(T)) where {T<:AbstractFloat}
    return start!(ls, convert(T, f0), convert(T, g0), convert(T, stp),
                  convert(T, stpmin), convert(T, stpmax))
end

"""
    iterate!(ls, stp, f, g) -> task

performs one iteration of the line search implemented by `ls` with `f` and `g`
the value and directional derivative of the objective function for the step
`stp` (which must be equal to `getstep(ls)`).

The returned `task` is one of the following symbolic values:

* `:SEARCH` to indicate that the caller should take the new trial step, given
  by `getstep(ls)`, and call `iterate!` with the updated values of the
  objective function and of its directional derivative.

* `:CONVERGENCE` to indicate that the line search criterion holds.

* `:WARNING` to indicate that no further progresses are possible.

""" iterate!

#-------------------------------------------------------------------------------

mutable struct ArmijoLineSearch{T<:AbstractFloat} <: LineSearch{T}
    task::Symbol
    status::Symbol
    ftol::T
    finit::T
    ginit::T
    stp::T
    stpmin::T
end

"""

    ArmijoLineSearch(T; ftol = 1e-4) -> ls

yields a line search instance for floating-point type `T` which implements
Armijo's method.  Keyword `ftol` can be used to specify a nonnegative tolerance
for the sufficient decrease condition.

## Description

The objective of Armijo's method is to find a step `stp` that satisfies the
sufficient decrease condition (1st Wolfe condition):

    f(stp) ≤ f(0) + ftol*stp*f'(0),

where `stp` is smaller or equal the initial step (backtracking).  The method
consists in reducing the step (by a factor 2) until the criterion holds.

## References

* L. Armijo, "*Minimization of functions having Lipschitz continuous first
  partial derivatives*" in Pacific Journal of Mathematics, vol. 16, pp. 1–3
  (1966).

"""
function ArmijoLineSearch(::Type{T} = Float;
                          ftol::Real = 1e-4) where {T<:AbstractFloat}
    @assert 0 < ftol ≤ 1/2
    ArmijoLineSearch{T}(:START, :NOT_STARTED, ftol, 0, 0, 0, 0)
end

# Armijo's line search does not use the directional derivative to refine the
# step.
usederivatives(::Type{<:ArmijoLineSearch}) = false

# start! method for Armijo's line search has the same implementation as the
# Moré & Toraldo line search.

function iterate!(ls::ArmijoLineSearch, stp::Real, f::Real, g::Real)
    @assert ls.task == :SEARCH
    @assert stp == ls.stp
    if f ≤ ls.finit + ls.ftol*stp*ls.ginit
        return report!(ls, :CONVERGENCE, :FIRST_WOLFE_CONDITION_SATISFIED)
    elseif stp > ls.stpmin
        ls.stp = max(stp/2, ls.stpmin)
        return report!(ls, :SEARCH, :COMPUTE_FG)
    else
        ls.stp = ls.stpmin
        return report!(ls, :WARNING, :STP_EQ_STPMIN)
    end
end

#-------------------------------------------------------------------------------

mutable struct MoreToraldoLineSearch{T<:AbstractFloat} <: LineSearch{T}
    task::Symbol
    status::Symbol
    gamma1::T
    gamma2::T
    ftol::T
    finit::T
    ginit::T
    stp::T
    stpmin::T
    strict::Bool
end

"""

    MoreToraldoLineSearch(T; ...) -> ls

yields a backtracking line search proposed by Moré & Toraldo (1991) and which
is based on a safeguarded quadratic interpolation.  Argument `T` is the
floating-point type for the computations.

Keyword `ftol` can be used to specify a nonnegative tolerance for the
sufficient decrease condition.

Keyword `gamma` can be used to specify a tuple of two values to safeguard the
quadratic interpolation of the step and such that `0 < gamma[1] < gamma[2] <
1`).  In Moré & Toraldo (1991, GPCG algorithm) `gamma = (0.01,0.5)` while in
Birgin et al. (2000, SPG2 algorithm) `gamma = (0.1,0.9)`.  The default settings
are `gamma = (0.1, 0.5)`.

If keyword `strict` is set to `false`, a variant of the method is used which
takes a bissection step whenever the curvature of the model is not striclty
positive.  By default, `strict` is `true`.

## Description

The objective of Moré & Toraldo method is to find a step `stp` that satisfies
the sufficient decrease condition (1st Wolfe condition):

    f(stp) ≤ f(0) + ftol*stp*f'(0),

where `stp` is smaller or equal the initial step (backtracking).  The principle
of the method is to reduce the step by a safeguarded quadratic interpolation
(or by a bissection if the minimum of the quadratic interpolation falls outside
the current search interval) until the criterion is met.

The new step is computed as follows:

    newstp = clamp(fact, gamma[1], gamma[2])*stp

where `fact` is the factor by which the current step should be multiplied to
reach the minimum of the quadratic interpolation of the objective function.

## References

* J.J. Moré & G. Toraldo, "*On the Solution of Large Quadratic Programming
  Problems with Bound Constraints*", SIAM J. Optim., vol. 1, pp. 93–113 (1991).

* E.G. Birgin, J.M. Martínez & M. Raydan, "*Nonmonotone Spectral Projected
  Gradient Methods on Convex Sets*", SIAM J. Optim., vol. 10, pp. 1196–1211
  (2000).

"""
function MoreToraldoLineSearch(::Type{T} = Float;
                               ftol::Real = 1e-4,
                               strict::Bool = false,
                               gamma::NTuple{2,Real} = (0.1,0.5)
                               ) where {T<:AbstractFloat}
    @assert 0 < ftol ≤ 1/2
    @assert 0 < gamma[1] < gamma[2] < 1
    MoreToraldoLineSearch{T}(:START, :NOT_STARTED, gamma[1], gamma[2], ftol,
                             0, 0, 0, 0, strict)
end

# Backtracking line search does not use the directional derivative to refine
# the step but it does use the initial derivative.
usederivatives(::Type{<:MoreToraldoLineSearch}) = false

function start!(ls::Union{MoreToraldoLineSearch,ArmijoLineSearch},
                f0::Real, g0::Real, stp::Real, stpmin::Real, stpmax::Real)
    # Check the input arguments for errors.
    if stpmax < stpmin
        throw(ArgumentError("STPMAX < STPMIN"))
    end
    if stpmin < 0
        throw(ArgumentError("STPMIN < 0"))
    end
    if stp > stpmax
         throw(ArgumentError("STP > STPMAX"))
    end
    if stp < stpmin
         throw(ArgumentError("STP < STPMIN"))
    end
    if g0 ≥ 0
        throw(ArgumentError("Not a descent direction"))
        #return report!(ls, :ERROR, :NOT_DESCENT)
    end

    # Store objective function and directional derivative at start of line
    # search.
    ls.finit = f0
    ls.ginit = g0
    ls.stp = stp
    ls.stpmin = stpmin
    return report!(ls, :SEARCH, :COMPUTE_FG)
end

function iterate!(ls::MoreToraldoLineSearch, stp::Real, f::Real, g::Real)
    @assert ls.task == :SEARCH
    @assert stp == ls.stp
    if f ≤ ls.finit + ls.ftol*stp*ls.ginit
        return report!(ls, :CONVERGENCE, :FIRST_WOLFE_CONDITION_SATISFIED)
    elseif stp > ls.stpmin
        q = -ls.ginit*stp
        r = f - (ls.finit - q)
        if (ls.strict ? r > 0 && q > 0 : r != 0)
            # Safeguarded quadratic interpolation.
            stp *= clamp(q/(r + r), ls.gamma1, ls.gamma2)
        else
            # Bissection.
            stp /= 2
        end
        ls.stp = max(stp, ls.stpmin)
        return report!(ls, :SEARCH, :COMPUTE_FG)
    else
        ls.stp = ls.stpmin
        return report!(ls, :WARNING, :STP_EQ_STPMIN)
    end
end

#-------------------------------------------------------------------------------

"""

`MoreThuenteLineSearch{T}` defines the work-space for Moré & Thuente line
search.

Line search methods which need to compute a cubic step via `cstep!` need to
store some parameters in an instance of `LineSearchInterval` which has the
following members:

* `lower` and `upper` contain respectively a lower and an upper bound for the
  step;

* `stx`, `fx`, `dx` contain the values of the step, function, and derivative at
  the best step so far;

* `sty`, `fy`, `dy` contain the value of the step, function, and derivative at
  the other endpoint`sty`;

* `brackt` indicates whether a minimum is bracketed in the interval
  `(stx,sty)`.

At the start of the line search, `brackt = false`, `stx = sty = 0` while `fx =
fy = f0` and `dx = dy = g0` the value of the function and its derivative for
the step `stp = 0`.  Then `cstep!` can be called (typically in the `iterate!`
method) to compute a new step (based on cubic or quadratic interpolation) and
to maintain the interval of search.

"""
mutable struct MoreThuenteLineSearch{T<:AbstractFloat} <: LineSearch{T}
    task::Symbol
    status::Symbol
    ftol::T
    gtol::T
    xtol::T
    stpmin::T
    stpmax::T
    finit::T
    ginit::T
    stp::T
    stx::T
    fx::T
    gx::T
    sty::T
    fy::T
    gy::T
    lower::T     # current lower bound for the step
    upper::T     # current upper bound for the step
    width::T     # current width of the interval
    oldwidth::T  # previous width of the interval
    stage::Int
    brackt::Bool # minimum has been bracketed?
end

function MoreThuenteLineSearch(::Type{T} = Float;
                               ftol::Real = 0.001,
                               gtol::Real = 0.9,
                               xtol::Real = 0.1) where {T<:AbstractFloat}
    @assert 0 < ftol < 1
    @assert 0 < gtol < 1
    @assert 0 < xtol < 1
    MoreThuenteLineSearch{T}(:START, :NOT_STARTED,
                             ftol, gtol, xtol, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                             0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                             0, false)
end

# Moré & Thuente line search does use the directional derivative to refine the
# step.
usederivatives(::Type{<:MoreThuenteLineSearch}) = true

# A few constants.
const XTRAPL = 1.1
const XTRAPU = 4.0
const STPMAX = 1e20

function start!(ls::MoreThuenteLineSearch{T},
                f0::Real, g0::Real, stp::Real;
                ftol::Real = ls.ftol,
                gtol::Real = ls.gtol,
                xtol::Real = ls.xtol,
                stpmin::Real = zero(T),
                stpmax::Real = stp*STPMAX) where {T<:AbstractFloat}
    start!(ls, f0, g0, stp, ftol, gtol, xtol, stpmin, stpmax)
end

function start!(ls::MoreThuenteLineSearch,
                f0::Real, g0::Real, stp::Real,
                ftol::Real, gtol::Real, xtol::Real,
                stpmin::Real, stpmax::Real)

    # Check the input arguments for errors.
    if stpmax < stpmin
        throw(ArgumentError("STPMAX < STPMIN"))
    end
    if stpmin < 0
        throw(ArgumentError("STPMIN < 0"))
    end
    if xtol < 0
         throw(ArgumentError("XTOL < 0"))
    end
    if ftol ≤ 0
         throw(ArgumentError("FTOL ≤ 0"))
    end
    if gtol ≤ 0
         throw(ArgumentError("GTOL ≤ 0"))
    end
    if g0 ≥ 0
         throw(ArgumentError("Not a descent direction"))
    end
    if stp > stpmax
         throw(ArgumentError("STP > STPMAX"))
    end
    if stp < stpmin
         throw(ArgumentError("STP < STPMIN"))
    end

    # Initialize local variables.
    # The variables STX, FX, GX contain the values of the step,
    # function, and derivative at the best step.
    # The variables STY, FY, GY contain the value of the step,
    # function, and derivative at STY.
    # The variables STP, F, G contain the values of the step,
    # function, and derivative at STP.
    ls.ftol      = ftol
    ls.gtol      = gtol
    ls.xtol      = xtol
    ls.stpmin    = stpmin
    ls.stpmax    = stpmax
    ls.finit     = f0
    ls.ginit     = g0
    ls.stp       = stp
    ls.stx       = 0
    ls.fx        = f0
    ls.gx        = g0
    ls.sty       = 0
    ls.fy        = f0
    ls.gy        = g0
    ls.lower     = 0
    ls.upper     = stp + stp*XTRAPU
    ls.width     = stpmax - stpmin
    ls.oldwidth  = 2*(stpmax - stpmin)
    ls.stage     = 1
    ls.brackt    = false
    report!(ls, :SEARCH, :COMPUTE_FG)
end

function iterate!(ls::MoreThuenteLineSearch{T},
                  stp::Real, f::Real, g::Real) where {T<:AbstractFloat}
    return iterate!(ls, convert(T, stp), convert(T, f), convert(T, g))
end

function iterate!(ls::MoreThuenteLineSearch{T},
                  stp::T, f::T, g::T) where {T<:AbstractFloat}

    @assert ls.task == :SEARCH
    @assert stp == ls.stp

    P66 = T(0.66)
    HALF = T(0.5)
    ZERO = zero(T)

    # If psi(stp) ≤ 0 and f'(stp) ≥ 0 for some step, then the algorithm enters
    # the second stage.
    gtest = ls.ftol*ls.ginit
    ftest = ls.finit + stp*gtest
    if ls.stage == 1 && f ≤ ftest && g ≥ ZERO
        ls.stage = 2
    end

    # Test for termination: convergence or warnings.
    if f ≤ ftest && abs(g) ≤ -ls.gtol*ls.ginit
        return report!(ls, :CONVERGENCE, :STRONG_WOLFE_SATISFIED)
    end
    if stp == ls.stpmin && (f > ftest || g ≥ gtest)
        return report!(ls, :WARNING, :STP_EQ_STPMIN)
    end
    if stp == ls.stpmax && f ≤ ftest && g ≤ gtest
        return report!(ls, :WARNING, :STP_EQ_STPMAX)
    end
    if ls.brackt
        if ls.upper - ls.lower ≤ ls.xtol*ls.upper
            return report!(ls, :WARNING, :XTOL_TEST_SATISFIED)
        end
        if stp ≤ ls.lower || stp ≥ ls.upper
            return report!(ls, :WARNING, :ROUNDING_ERRORS)
        end
    end

    # A modified function is used to predict the step during the first stage if
    # a lower function value has been obtained but the decrease is not
    # sufficient.

    if ls.stage == 1 && f ≤ ls.fx && f > ftest

        # Call CSTEP to update STX, STY, and to compute the new step for the
        # modified function and its derivatives.
        info, ls.brackt,
        ls.stx, fxm, gxm,
        ls.sty, fym, gym,
        stp = cstep(ls.brackt, ls.lower, ls.upper,
                    ls.stx, ls.fx - ls.stx*gtest, ls.gx - gtest,
                    ls.sty, ls.fy - ls.sty*gtest, ls.gy - gtest,
                    stp, f - stp*gtest, g - gtest)

        # Reset the function and derivative values for F.
        ls.fx = fxm + ls.stx*gtest
        ls.fy = fym + ls.sty*gtest
        ls.gx = gxm + gtest
        ls.gy = gym + gtest

    else

        # Call CSTEP to update STX, STY, and to compute the new step.
        info, ls.brackt,
        ls.stx, ls.fx, ls.gx,
        ls.sty, ls.fy, ls.gy,
        stp = cstep(ls.brackt, ls.lower, ls.upper,
                    ls.stx, ls.fx, ls.gx,
                    ls.sty, ls.fy, ls.gy,
                    stp, f, g)

    end

    # Decide if a bisection step is needed.
    if ls.brackt
        wcur = abs(ls.sty - ls.stx)
        if wcur ≥ P66*ls.oldwidth
            stp = ls.stx + HALF*(ls.sty - ls.stx)
        end
        ls.oldwidth = ls.width
        ls.width = wcur
    end

    # Set the minimum and maximum steps allowed for STP.
    if ls.brackt
        if ls.stx ≤ ls.sty
            ls.lower = ls.stx
            ls.upper = ls.sty
        else
            ls.lower = ls.sty
            ls.upper = ls.stx
        end
    else
        ls.lower = stp + XTRAPL*(stp - ls.stx)
        ls.upper = stp + XTRAPU*(stp - ls.stx)
    end

    # Force the step to be within the bounds STPMAX and STPMIN.
    stp = clamp(stp, ls.stpmin, ls.stpmax)

    # If further progress is not possible, let STP be the best point obtained
    # during the search.
    if (ls.brackt && (stp ≤ ls.lower || stp ≥ ls.upper ||
                      ls.upper - ls.lower ≤ ls.xtol*ls.upper))
        stp = ls.stx
    end

    # Save next step to try.
    ls.stp = stp

    # Obtain another function and derivative.
    return report!(ls, :SEARCH, :COMPUTE_FG)
end

"""
# Compute a safeguarded cubic step.

    cstep(brackt, stpmin, stpmax,
          stx, fx, dx,
          sty, fy, dy,
          stp, fp, dp) ->  info, brackt, stx,fx,dx, sty,fy,dy, stp

The function `cstep` computes a safeguarded step for a search procedure and
updates an interval that contains a step that satisfies a sufficient decrease
and a curvature condition.  The algorithm is described in:

* J.J. Moré and D.J. Thuente, "*Line search algorithms with guaranteed
  sufficient decrease*" in ACM Transactions on Mathematical Software, vol. 20,
  pp. 286–307 (1994).

The parameter `stx` contains the step with the least function value.  If
`brackt` is set to true (*i.e.,* non-zero) then a minimizer has been bracketed
in an interval with endpoints `stx` and `sty`.  The parameter `stp` contains
the current step.  The subroutine assumes that if `brackt` is true then:

     min(stx,sty) < stp < max(stx,sty),

and that the derivative at `stx` is negative in the direction of the step.

On output, the updated parameters are returned.

Parameter `brackt` specifies if a minimizer has been bracketed.  Initially
`brackt` must be set to false.  The returned value of `brackt` indicates if a
minimizer has been bracketed.

Parameters `stpmin` and `stpmax` specify the lower and the upper bounds for the
step.

Parameters `stx`, `fx` and `dx` specify the step, the function and the
derivative at the best step obtained so far.  The derivative must be negative
in the direction of the step, that is, `dx` and `stp - stx` must have opposite
signs.  On return, these parameters are updated appropriately.

Parameters `sty`, `fy` and `dy` specify the step, the function and the
derivative at the other endpoint of the interval of uncertainty.  On return,
these parameters are updated appropriately.

Parameters `stp`, `fp` and `dp` specify the step, the function and the
derivative at the current step.  If `brackt` is true, then `stp` must be
between `stx` and `sty`.  On return, these parameters are updated
appropriately.  The returned value of `stp` is the next trial step.

The returned value `info` indicates which case occured for computing the new
step.

"""
function cstep(brackt::Bool, stpmin::Real, stpmax::Real,
               stx::Real, fx::Real, dx::Real,
               sty::Real, fy::Real, dy::Real,
               stp::Real, fp::Real, dp::Real)
    cstep(brackt,
          promote(stpmin, stpmax, stx, fx, dx, sty, fy, dy, stp, fp, dp)...)
end

function cstep(brackt::Bool, stpmin::T, stpmax::T,
               stx::T, fx::T, dx::T,
               sty::T, fy::T, dy::T,
               stp::T, fp::T, dp::T) where {T<:AbstractFloat}

    ZERO = zero(T)
    ONE = one(T)
    TWO = convert(T, 2)
    THREE = convert(T, 3)
    P66 = convert(T, 0.66)
    local theta::T, gamma::T, p::T, q::T, r::T, s::T, t::T

    info = 0

    # Check the input parameters for errors.
    if brackt && (stp ≤ min(stx, sty) || stp ≥ max(stx, sty))
        throw(ArgumentError("STP outside bracket (STX,STY)"))
    end
    if dx*(stp - stx) ≥ ZERO
        throw(ArgumentError("descent condition violated"))
    end
    if stpmax < stpmin
        throw(ArgumentError("STPMAX < STPMIN"))
    end

    # Determine if the derivatives have opposite sign.
    opposite = ((dx < ZERO && dp > ZERO) || (dx > ZERO && dp < ZERO))

    if fp > fx
        # First case.  A higher function value.  The minimum is bracketed.  If
        # the cubic step is closer to STX than the quadratic step, the cubic
        # step is taken, otherwise the average of the cubic and quadratic steps
        # is taken.
        info = 1
        theta = THREE*(fx - fp)/(stp - stx) + dx + dp
        s = max(abs(theta), abs(dx), abs(dp))
        gamma = s*sqrt((theta/s)^2 - (dx/s)*(dp/s))
        if stp < stx
            gamma = -gamma
        end
        p = (gamma - dx) + theta
        q = ((gamma - dx) + gamma) + dp
        r = p/q
        stpc = stx + r*(stp - stx)
        stpq = stx + ((dx/((fx - fp)/(stp - stx) + dx))/TWO)*(stp - stx)
        if abs(stpc - stx) < abs(stpq - stx)
            stpf = stpc
        else
            stpf = stpc + (stpq - stpc)/TWO
        end
        brackt = true
    elseif opposite
        # Second case.  A lower function value and derivatives of opposite
        # sign.  The minimum is bracketed.  If the cubic step is farther from
        # STP than the secant (quadratic) step, the cubic step is taken,
        # otherwise the secant step is taken.
        info = 2
        theta = THREE*(fx - fp)/(stp - stx) + dx + dp
        s = max(abs(theta), abs(dx), abs(dp))
        gamma = s*sqrt((theta/s)^2 - (dx/s)*(dp/s))
        if stp > stx
            gamma = -gamma
        end
        p = (gamma - dp) + theta
        q = ((gamma - dp) + gamma) + dx
        r = p/q
        stpc = stp + r*(stx - stp)
        stpq = stp + (dp/(dp - dx))*(stx - stp)
        if abs(stpc - stp) > abs(stpq - stp)
            stpf = stpc
        else
            stpf = stpq
        end
        brackt = true
    elseif abs(dp) < abs(dx)
        # Third case.  A lower function value, derivatives of the same sign,
        # and the magnitude of the derivative decreases.  The cubic step is
        # computed only if the cubic tends to infinity in the direction of the
        # step or if the minimum of the cubic is beyond STP.  Otherwise the
        # cubic step is defined to be the secant step.
        info = 3
        theta = THREE*(fx - fp)/(stp - stx) + dx + dp
        s = max(abs(theta), abs(dx), abs(dp))
        # The case GAMMA = 0 only arises if the cubic does not tend to infinity
        # in the direction of the step.
        t = (theta/s)^2 - (dx/s)*(dp/s)
        gamma = (t > ZERO ? s*sqrt(t) : ZERO)
        if stp > stx
            gamma = -gamma
        end
        p = (gamma - dp) + theta
        #q = ((gamma - dp) + gamma) + dx
        q = (gamma + (dx - dp)) + gamma
        r = p/q
        if r < ZERO && gamma != ZERO
            stpc = stp + r*(stx - stp)
        elseif stp > stx
            stpc = stpmax
        else
            stpc = stpmin
        end
        stpq = stp + (dp/(dp - dx))*(stx - stp)
        if brackt
            # A minimizer has been bracketed.  If the cubic step is closer to
            # STP than the secant step, the cubic step is taken, otherwise the
            # secant step is taken.
            stpf = abs(stpc - stp) < abs(stpq - stp) ? stpc : stpq
            t = stp + P66*(sty - stp)
            if stp > stx ? stpf > t : stpf < t
                stpf = t
            end
        else
            # A minimizer has not been bracketed. If the cubic step is farther
            # from stp than the secant step, the cubic step is taken, otherwise
            # the secant step is taken.
            stpf = abs(stpc - stp) > abs(stpq - stp) ? stpc : stpq
            stpf = max(stpmin, min(stpf, stpmax))
        end
    else
        # Fourth case.  A lower function value, derivatives of the same sign,
        # and the magnitude of the derivative does not decrease.  If the
        # minimum is not bracketed, the step is either STPMIN or STPMAX,
        # otherwise the cubic step is taken.
        info = 4
        if brackt
            theta = THREE*(fp - fy)/(sty - stp) + dy + dp
            s = max(abs(theta), abs(dy), abs(dp))
            t = theta/s
            gamma = s*sqrt((theta/s)^2 - (dy/s)*(dp/s))
            if stp > sty
                gamma = -gamma
            end
            p = (gamma - dp) + theta
            q = ((gamma - dp) + gamma) + dy
            r = p/q
            stpc = stp + r*(sty - stp)
            stpf = stpc
        elseif stp > stx
            stpf = stpmax
        else
            stpf = stpmin
        end
    end

    # Update the interval which contains a minimizer and guess for next step.
    if fp > fx
        sty = stp
        fy = fp
        dy = dp
    else
        if opposite
            sty = stx
            fy = fx
            dy = dx
        end
        stx = stp
        fx = fp
        dx = dp
    end
    stp = stpf
    return (info, brackt, stx,fx,dx, sty,fy,dy, stp)
end

end # module


text:
export Non_negative
export Map_0_1
export Piecewise_positive
export Pow4_positive
export Abs_positive

 # All these functions return a mapping function and
 # the inverse of it
 # they are used to map the real numbers to non-negative real numbers

"""
    Non_negative()

Returns a function and an inverse function inverse function
to map numbers to non-negative numbers.
We use a parabola.

# Examples
```julia-repl
julia> p, p_inv = Non_negative()
(DeconvOptim.var"#5#7"(), DeconvOptim.var"#6#8"())

julia> x = [-1, 2, -3]
3-element Array{Int64,1}:
 -1
  2
 -3

julia> p(x)
3-element Array{Int64,1}:
 1
 4
 9

julia> p_inv(p(x))
3-element Array{Float64,1}:
 1.0
 2.0
 3.0
```
"""
function Non_negative()
    #return x -> map(abs2, x), parab_inv
    return parab, parab_inv
end

parab(x) = abs2.(x)
parab_inv(x) = sqrt.(x)
 # define custom adjoint for parab because of
 # slow broadcasting
function ChainRulesCore.rrule(::typeof(parab), x)
    Y = parab(x)
    function aux_pullback(barx)
        return zero(eltype(Y)), (2 .* barx) .* x
    end
    return Y, aux_pullback
end

"""
    Map_0_1()

Returns a function and an inverse function
to map numbers to an interval between 0 and 1.
via an exponential function.

"""
function Map_0_1()
    return f01, f01_inv
end

f01(x) = 1 .- exp.(.- x.^2)
f01_inv(y) = sqrt.(.- log.(1 .- y))

"""
    Piecewise_positive()

Returns a function and an inverse function
to map numbers to larger than 0 via
two function stitched together.
"""
function Piecewise_positive()
    return f_pw_pos, f_pw_pos_inv
end

f_pw_pos(x) = ifelse.(x .> 0, one(eltype(x)) .+ x,one(eltype(x))./(one(eltype(x)).-x))
f_pw_pos_grad(x) = ifelse.(x .> 0, one(eltype(x)) , one(eltype(x))./abs2.(one(eltype(x)).-x))
f_pw_pos_inv(y) = ifelse.(y .> 1, y .- one(eltype(y)), one(eltype(y)) .- one(eltype(y))./y)

function ChainRulesCore.rrule(::typeof(f_pw_pos), x)
    Y = f_pw_pos(x)
    function aux_pullback(barx)
        return zero(eltype(Y)), barx .* f_pw_pos_grad(x)
    end
    return Y, aux_pullback
end

"""
    Pow4_positive()

Returns a function and an inverse function
to map numbers to larger than 0 with `abs2.(abs2.(x))`

"""
function Pow4_positive()
    return f_pow4, f_pow4_inv
end

f_pow4(x) = abs2.(abs2.(x))
f_pow4_inv(y) = sqrt.(sqrt.(y))

"""
    Abs_positive()

Returns a function and an inverse function
to map numbers to larger than 0.
"""
function Abs_positive()
    return f_abs, f_abs_inv
end

f_abs(x) = abs.(x)
f_abs_inv(y) = abs.(y)


text:
# Autogenerated wrapper script for libsvm_jll for x86_64-w64-mingw32
export libsvm, svm_predict, svm_scale, svm_train

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("libsvm")
JLLWrappers.@declare_library_product(libsvm, "libsvm.dll")
JLLWrappers.@declare_executable_product(svm_predict)
JLLWrappers.@declare_executable_product(svm_scale)
JLLWrappers.@declare_executable_product(svm_train)
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libsvm,
        "bin\\libsvm.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        svm_predict,
        "bin\\svm-predict.exe",
    )

    JLLWrappers.@init_executable_product(
        svm_scale,
        "bin\\svm-scale.exe",
    )

    JLLWrappers.@init_executable_product(
        svm_train,
        "bin\\svm-train.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    StochasticModel

A mathematical model of a stochastic optimization problem.
"""
struct StochasticModel{N, P <: NTuple{N, StageParameters}}
    parameters::P
    generator::Function

    function StochasticModel(generator::Function, parameters::Vararg{StageParameters,N}) where N
        return new{N,typeof(parameters)}(parameters, generator)
    end
end
nstages(::StochasticModel{N}) where N = N
"""
    @stochastic_model(def)

Define a stochastic model capable of instantiating stochastic programs, using the syntax
```julia
sm = @stochastic_model begin
    ...
    @stage x begin
      ...
    end
    ...
end
```
where the inner blocks are [`@stage`](@ref) blocks. At least two stages must be specified in consecutive order. A stochastic model object can later be used to [`instantiate`](@ref) stochastic programs using a given set of scenarios or to create [`SAA`](@ref) models using samplers.

## Examples

The following defines a stochastic model consisitng of the first stage model given by:
```math
  minimize 100x₁ + 150x₂
    s.t  x₁ + x₂ ≤ 120
         x₁ ≥ 40
         x₂ ≥ 20
```
and the second-stage model given by:
```math
  minimize q₁(ξ)y₁ + q₂(ξ)y₂
    s.t  6y₁ + 10y₂ ≤ 60x₁
         8y₁ + 5y₂ ≤ 60x₂
         0 ≤ y₁ ≤ d₁(ξ)
         0 ≤ y₂ ≤ d₂(ξ)
```
where ``q₁(ξ), q₂(ξ), d₁(ξ), d₂(ξ)`` depend on the scenario ``ξ``.

```julia
sm = @stochastic_model begin
    @stage 1 begin
        @variable(model, x₁ >= 40)
        @variable(model, x₂ >= 20)
        @objective(model, Min, 100*x₁ + 150*x₂)
        @constraint(model, x₁ + x₂ <= 120)
    end
    @stage 2 begin
        @decision x₁ x₂
        @uncertain q₁ q₂ d₁ d₂
        @variable(model, 0 <= y₁ <= d₁)
        @variable(model, 0 <= y₂ <= d₂)
        @objective(model, Min, q₁*y₁ + q₂*y₂)
        @constraint(model, 6*y₁ + 10*y₂ <= 60*x₁)
        @constraint(model, 8*y₁ + 5*y₂ <= 80*x₂)
    end
end
```

See also: [`@stage`](@ref), [`@parameters`](@ref), [`@decision`](@ref), [`@uncertain`](@ref)
"""
macro stochastic_model(def)
    stage = 0
    scenariodef = Expr(:block)
    paramdefs = Vector{Expr}()
    def = prewalk(prettify(def)) do x
        x = if @capture(x, @stage n_ arg_)
            if @capture(arg, sp_ = def_)
                x
            else
                stage == 0 && (n == 1 || error("A first stage must be defined."))
                stage == n - 1 || error("Define the stages in coherent order.")
                stage += 1
                push!(paramdefs, :(StageParameters()))
                return :(@stage $n sp = $arg)
            end
        else
            x
        end
        paramdef = if @capture(x, @parameters arg_)
            names = Vector{Symbol}()
            default = Vector{Expr}()
            for paramdef in prettify(arg).args
                if @capture(paramdef, key_Symbol = val_)
                    push!(names, key)
                    push!(default, paramdef)
                elseif @capture(paramdef, key_Symbol)
                    push!(names, key)
                else
                    error("Incorrect @parameters syntax. Specify parameter names, possibly with a default value.")
                end
            end
            :(StageParameters($names; $(default...)))
        elseif @capture(x, @parameters args__)
            args = convert(Vector{Symbol}, args)
            :(StageParameters($args))
        else
            nothing
        end
        if paramdef != nothing
            paramdefs[stage] = paramdef
        end
        scenariodef = if @capture(x, @uncertain var_Symbol::t_Symbol = def_)
            esc(@q begin
                @scenario $t = $def
            end)
        else
            scenariodef
        end
        return x
    end
    stage >= 2 || error("Define at least two stages.")
    code = @q begin
        $scenariodef
        StochasticModel($(esc.(paramdefs)...)) do $(esc(:sp))
            $(esc(def))
        end
    end
    return prettify(code)
end

# Printing #
# ========================== #
function Base.show(io::IO, stochasticmodel::StochasticModel)
    println(io, "Multi-stage Stochastic Model")
end
function Base.show(io::IO, stochasticmodel::StochasticModel{2})
    modelstr = "minimize f₀(x) + 𝔼[f(x,ξ)]
  x∈𝒳

where

f(x,ξ) = min  f(y; x, ξ)
              y ∈ 𝒴 (x, ξ)"
    print(io, "Two-Stage Stochastic Model\n\n")
    println(io, modelstr)
end
# ========================== #