text:
# Based on crying baby problem (Kochenderfer 2015)
# POMDP, Robust POMDP, and Robust ρPOMDP formulations

struct BabyPOMDP <: POMDP{Symbol, Symbol, Symbol}
    r_feed::Float64
    r_hungry::Float64
    discount::Float64
end

struct BabyRPOMDP <: RPOMDP{Symbol, Symbol, Symbol}
    r_feed::Float64
    r_hungry::Float64
    discount::Float64
    usize::Float64
end

struct BabyIPOMDP <: IPOMDP{Symbol, Symbol, Symbol}
    r_feed::Float64
    r_hungry::Float64
    discount::Float64
    inforeward::Vector{Vector{Float64}}
end

struct BabyRIPOMDP <: RIPOMDP{Symbol, Symbol, Symbol}
    r_feed::Float64
    r_hungry::Float64
    discount::Float64
    usize::Float64
    inforeward::Vector{Vector{Float64}}
end

const states_baby = [:hungry, :full]
const actions_baby = [:feed, :nothing]
const observations_baby = [:quiet, :crying, :yelling]
const inforeward_baby = [[1.0, -1.0], [-1.0, 1.0], [0.2, 0.2]]

BabyPOMDP(r_feed, r_hungry) = BabyPOMDP(r_feed, r_hungry, 0.9)
BabyPOMDP() = BabyPOMDP(-5.0, -10.0, 0.9)

BabyRPOMDP(r_feed, r_hungry) = BabyRPOMDP(r_feed, r_hungry, 0.9, 0.025)
BabyRPOMDP(err) = BabyRPOMDP(-5.0, -10.0, 0.9, err)
BabyRPOMDP() = BabyRPOMDP(-5.0, -10.0, 0.9, 0.025)

BabyIPOMDP(r_feed, r_hungry) = BabyIPOMDP(r_feed, r_hungry, 0.9, inforeward_baby)
BabyIPOMDP(alphas) = BabyIPOMDP(-5.0, -10.0, 0.9, alphas)
BabyIPOMDP() = BabyIPOMDP(-5.0, -10.0, 0.9, inforeward_baby)

BabyRIPOMDP(r_feed, r_hungry) = BabyRIPOMDP(r_feed, r_hungry, 0.9, 0.025, inforeward_baby)
BabyRIPOMDP(err::Float64, alphas::Vector{Vector{Float64}}) = BabyRIPOMDP(-5.0, -10.0, 0.9, err, alphas)
BabyRIPOMDP() = BabyRIPOMDP(-5.0, -10.0, 0.9, 0.025, inforeward_baby)

states(::Union{BabyPOMDP, BabyRPOMDP, BabyIPOMDP, BabyRIPOMDP}) = states_baby
actions(::Union{BabyPOMDP, BabyRPOMDP, BabyIPOMDP, BabyRIPOMDP}) = actions_baby
observations(::Union{BabyPOMDP, BabyRPOMDP, BabyIPOMDP, BabyRIPOMDP}) = observations_baby

n_states(::Union{BabyPOMDP, BabyRPOMDP, BabyIPOMDP, BabyRIPOMDP}) = 2
n_actions(::Union{BabyPOMDP, BabyRPOMDP, BabyIPOMDP, BabyRIPOMDP}) = 2
n_observations(::Union{BabyPOMDP, BabyRPOMDP, BabyIPOMDP, BabyRIPOMDP}) = 3

state_index(::Union{BabyPOMDP, BabyRPOMDP, BabyIPOMDP, BabyRIPOMDP}, s::Symbol) = s == :hungry ? 1 : 2
action_index(::Union{BabyPOMDP, BabyRPOMDP, BabyIPOMDP, BabyRIPOMDP}, a::Symbol) = a == :feed ? 1 : 2
observation_index(::Union{BabyPOMDP, BabyRPOMDP, BabyIPOMDP, BabyRIPOMDP}, z::Symbol) = z == :quiet ? 1 : z == :crying ? 2 : 3
obs_index(prob::Union{BabyPOMDP, BabyRPOMDP, BabyIPOMDP, BabyRIPOMDP}, z::Symbol) = observation_index(prob, z)

initial_state_distribution(::Union{BabyPOMDP, BabyRPOMDP, BabyIPOMDP, BabyRIPOMDP}) = SparseCat([:hungry, :full], [0.0, 1.0])
initial_belief(::Union{BabyPOMDP, BabyRPOMDP, BabyIPOMDP, BabyRIPOMDP}) = [0.0, 1.0]

# Transition functions
# Nominal transition function distributions
const pϵ_baby = 1e-6
const tdist_baby = [SparseCat([:hungry, :full], [0.0, 1.0]) SparseCat([:hungry, :full], [1.0, 0.0]);
                    SparseCat([:hungry, :full], [0.0, 1.0]) SparseCat([:hungry, :full], [0.1, 0.9])]
# Nominal transition function array tarray_baby[s,a,s'] = Pr(s' | s ,a)
const tarray_baby = cat(3, [0.0 1.0; 0.0 0.1], [1.0 0.0; 1.0 0.9]) # P(hungry | s, a), P(full | s, a)

# Nominal transitions
transition(prob::Union{BabyIPOMDP, BabyPOMDP}, s::Symbol, a::Symbol) = tdist_baby[state_index(prob, s), action_index(prob, a)]
transition(prob::Union{BabyIPOMDP, BabyPOMDP}, s::Symbol, a::Symbol, sp::Symbol) = tarray_baby[state_index(prob, s), action_index(prob, a), state_index(prob, sp)]

# Robust transitions
# P(z,t|s,a) = P(z|s,a,t)P(t|s,a)
# PztLower = PzLower * PtLower
# PtLower = Plower(t|s,a) = t_lower[s,a,sp]
function transition(prob::Union{BabyRPOMDP, BabyRIPOMDP}, s::Symbol, a::Symbol)
    si, ai = state_index(prob, s), action_index(prob, a)
    plower = max.(tarray_baby[si,ai,:] - prob.usize, 0.0 + pϵ_baby)
    pupper = min.(tarray_baby[si,ai,:] + prob.usize, 1.0 - pϵ_baby)
    plower, pupper
end
function transition(prob::Union{BabyRPOMDP, BabyRIPOMDP}, s::Symbol, a::Symbol, sp::Symbol)
    si, ai, spi = state_index(prob, s), action_index(prob, a), state_index(prob, sp)
    plower = max.(tarray_baby[si, ai, spi] - prob.usize, 0.0 + pϵ_baby)
    pupper = min.(tarray_baby[si, ai, spi] + prob.usize, 1.0 - pϵ_baby)
    plower, pupper
end

# Nominal observation function distributions
const odist_baby = [SparseCat([:quiet, :crying, :yelling], [0.3, 0.6, 0.1])
                    SparseCat([:quiet, :crying, :yelling], [0.8, 0.1, 0.1])]

# Nominal observation function array oarray_baby[a,sp,z] = Pr(z|a,sp)
const oarray_baby = cat(3, [0.3 0.8; 0.3 0.8], [0.6 0.1; 0.6 0.1], [0.1 0.1; 0.1 0.1]) # Pr(quiet|a,sp), # Pr(crying|a,sp), # Pr(yelling|a,sp)

# Nominal observations
observation(prob::Union{BabyIPOMDP, BabyPOMDP}, a::Symbol, sp::Symbol) = odist_baby[state_index(prob, sp)]
observation(prob::Union{BabyIPOMDP, BabyPOMDP}, a::Symbol, sp::Symbol, z::Symbol) = oarray_baby[action_index(prob, a), state_index(prob, sp), observation_index(prob, z)]

# Robust observations
# P(z,t|s,a) = P(z|s,a,t)P(t|s,a)
# PztLower = PzLower * PtLower
# PzLower = Plower(z|s,a,t) = o_lower[a,sp,z]
function observation(prob::Union{BabyRPOMDP, BabyRIPOMDP}, a::Symbol, sp::Symbol)
    ai, spi = action_index(prob, a), state_index(prob, sp)
    plower = max.(oarray_baby[ai,spi,:] - prob.usize, 0.0 + pϵ_baby)
    pupper = min.(oarray_baby[ai,spi,:] + prob.usize, 1.0 - pϵ_baby)
    plower, pupper
end
function observation(prob::Union{BabyRPOMDP, BabyRIPOMDP}, a::Symbol, sp::Symbol, z::Symbol)
    ai, spi, zi = action_index(prob, a), state_index(prob, sp), observation_index(prob, z)
    plower = max.(oarray_baby[ai, spi, zi] - prob.usize, 0.0 + pϵ_baby)
    pupper = min.(oarray_baby[ai, spi, zi] + prob.usize, 1.0 - pϵ_baby)
    plower, pupper
end

function dynamics(prob::Union{BabyPOMDP,BabyIPOMDP})
    ns = n_states(prob)
    nz = n_observations(prob)
    na = n_actions(prob)
    tarr = states(prob)
    sarr = states(prob)
    aarr = actions(prob)
    zarr = observations(prob)
    p = Array{Float64}(ns,nz,ns,na)
    plower = Array{Float64}(ns,nz,ns,na)
    pupper = Array{Float64}(ns,nz,ns,na)
    for tind in 1:ns, zind in 1:nz, sind in 1:ns, aind in 1:na
        pt = transition(prob, sarr[sind], aarr[aind], tarr[tind])
        pz = observation(prob, aarr[aind], tarr[tind], zarr[zind])
        p[tind, zind, sind, aind] = pz * pt
    end
    p
end

function dynamics(prob::Union{BabyRPOMDP,BabyRIPOMDP})
    ns = n_states(prob)
    nz = n_observations(prob)
    na = n_actions(prob)
    tarr = states(prob)
    sarr = states(prob)
    aarr = actions(prob)
    zarr = observations(prob)
    p = Array{Float64}(ns,nz,ns,na)
    plower = Array{Float64}(ns,nz,ns,na)
    pupper = Array{Float64}(ns,nz,ns,na)
    for tind in 1:ns, zind in 1:nz, sind in 1:ns, aind in 1:na
        ptl, ptu = transition(prob, sarr[sind], aarr[aind], tarr[tind])
        pzl, pzu = observation(prob, aarr[aind], tarr[tind], zarr[zind])
        plower[tind, zind, sind, aind] = pzl * ptl
        pupper[tind, zind, sind, aind] = pzu * ptu
    end
    plower, pupper
end

function reward(prob::Union{BabyPOMDP, BabyRPOMDP}, s::Symbol, a::Symbol)
    r = 0.0
    if s == :hungry
        r += prob.r_hungry
    end
    if a == :feed
        r += prob.r_feed
    end
    return r
end
reward(prob::Union{BabyPOMDP,BabyRPOMDP}, b::Vector{Float64}, s::Symbol, a::Symbol, sp::Symbol) = reward(prob,s,a)

function reward(prob::Union{BabyIPOMDP, BabyRIPOMDP}, b::Vector{Float64}, a::Symbol)
    rmax = -Inf
    for α in prob.inforeward
        rmax = max(rmax, dot(α, b))
    end
    rmax
end
reward(prob::Union{BabyIPOMDP,BabyRIPOMDP}, b::Vector{Float64}, s::Symbol, a::Symbol, sp::Symbol) = reward(prob,b,a)

function rewardalpha(prob::Union{BabyIPOMDP, BabyRIPOMDP}, b::Vector{Float64}, a::Symbol)
    rmax = -Inf
    ralpha = nothing
    for α in prob.inforeward
        rnext = dot(α, b)
        if rnext > rmax
            rmax = rnext
            ralpha = copy(α)
        end
    end
    ralpha
end

discount(p::Union{BabyPOMDP,BabyRPOMDP,BabyIPOMDP,BabyRIPOMDP}) = p.discount

# Simulation Functions
function initial_state(prob::Union{BabyPOMDP,BabyIPOMDP,BabyRPOMDP,BabyRIPOMDP}, rng::AbstractRNG)
    d = initial_state_distribution(prob)
    return rand(rng, d)
end

function generate_sor(prob::Union{BabyPOMDP,BabyIPOMDP}, b, s, a, rng::AbstractRNG)
    sp = rand(rng, transition(prob, s, a))
    o = rand(rng, observation(prob, a, sp))
    r = reward(prob, b, s, a, sp)
    sp, o, r
end

function generate_sor(prob::Union{BabyRPOMDP,BabyRIPOMDP}, b, s, a, rng::AbstractRNG)
    tdist = SparseCat(states(prob), psample(transition(prob, s, a)...))
    sp = rand(rng, tdist)
    odist = SparseCat(observations(prob), psample(observation(prob, a, sp)...))
    o = rand(rng, odist)
    r = reward(prob, b, s, a, sp)
    sp, o, r
end


text:
# Simple Adaptive Metropolis Proposal
# The proposal at each step is equal to a scalar multiple
# of the empirical posterior covariance plus a fixed, small covariance
# matrix epsilon which is also used for initial exploration.
#
# Reference:
#   H. Haario, E. Saksman, and J. Tamminen, "An adaptive Metropolis algorithm",
#   Bernoulli 7(2): 223-242 (2001)
mutable struct AMProposal{FT <: Real, CT <: AbstractMatrix{FT},
                          MNT <: AbstractMvNormal} <: Proposal{MNT}
    epsilon::CT
    scalefactor::FT
    proposal::MNT
    μ::Vector{FT}
    M::Matrix{FT}
    δ::Vector{FT}
    N::Int
end

function AMProposal(
    epsilon::AbstractMatrix{FT}, scalefactor=FT(2.38^2 / size(epsilon, 1))
) where {FT <: Real}
    sym = PDMat(Symmetric(epsilon))
    proposal = MvNormal(zeros(size(sym, 1)), sym)
    AMProposal(sym, scalefactor, proposal, zeros(size(sym,1)), zeros(size(sym)...),
               zeros(size(sym,1)), 0)
end

logratio_proposal_density(::AMProposal, params_prev, params) = 0

# When the proposal is initialised the empirical posterior covariance is zero
function trackstep!(proposal::AMProposal, params)
    proposal.μ .= params
    proposal.M .= 0
    proposal.proposal = MvNormal(zeros(size(proposal.μ)), proposal.epsilon)
    proposal.N = 1
end

# Recompute the empirical posterior covariance matrix
function trackstep!(proposal::AMProposal, params, ::Union{Val{true}, Val{false}})
    proposal.N += 1
    proposal.δ .= params .- proposal.μ
    proposal.μ .+= proposal.δ ./ proposal.N

    mul!(proposal.M, params .- proposal.μ, proposal.δ', true, true)

    prop_cov = proposal.M .* proposal.scalefactor ./ proposal.N .+ proposal.epsilon
    proposal.proposal = MvNormal(mean(proposal.proposal), Symmetric(prop_cov))
end

function propose(rng::Random.AbstractRNG, p::AMProposal, ::DensityModel)
    return rand(rng, p.proposal)
end

function propose(rng::Random.AbstractRNG, p::AMProposal, ::DensityModel, t)
    return t + rand(rng, p.proposal)
end

ext: jl


text:
module Optinum

using LinearAlgebra
using Test
using Markdown
using TestOptinum

include("Algorithme_De_Newton.jl")
export  Algorithme_De_Newton

include("Pas_De_Cauchy.jl")
export Pas_De_Cauchy

include("Gradient_Conjugue_Tronque.jl")
export Gradient_Conjugue_Tronque

include("Regions_De_Confiance.jl")
export Regions_De_Confiance

include("Lagrangien_Augmente.jl")
export Lagrangien_Augmente

end # module

max_stars_repo_name: faicaltoubali/ENSEEIHT

max_issues_repo_name: faicaltoubali/ENSEEIHT

max_forks_repo_name: faicaltoubali/ENSEEIHT


text:
function generate_get_expr(domain::Domain, state::State, term::Const,
                           varmap=Dict{Var,Any}(), state_var=:state)
    return :($state_var.$(term.name))
end

function generate_get_expr(domain::Domain, state::State, term::Var,
                           varmap=Dict{Var,Any}(), state_var=:state)
    return :($(varmap[term]))
end

function generate_get_expr(domain::Domain, state::State, term::Compound,
                           varmap=Dict{Var,Any}(), state_var=:state)
    indices = generate_fluent_ids(domain, state, term,
                                  get_fluent(domain, term.name), varmap, state_var)
    return :($state_var.$(term.name)[$(indices...)])
end

function generate_set_expr(domain::Domain, state::State, term::Const,
                           val, varmap=Dict{Var,Any}(), state_var=:state)
    if domain isa AbstractedDomain && domain.interpreter.autowiden
        prev_val = generate_get_expr(domain, state, term, varmap, :prev_state)
        return :(setfield!($state_var, $(QuoteNode(term.name)), widen($prev_val, $val)))
    else
        return :(setfield!($state_var, $(QuoteNode(term.name)), $val))
    end
end

function generate_set_expr(domain::Domain, state::State, term::Compound,
                           val, varmap=Dict{Var,Any}(), state_var=:state)
    indices = generate_fluent_ids(domain, state, term,
                                  get_fluent(domain, term.name), varmap, state_var)
    if domain isa AbstractedDomain && domain.interpreter.autowiden
        prev_val = generate_get_expr(domain, state, term, varmap, :prev_state)
        return :($state_var.$(term.name)[$(indices...)] = widen($prev_val, $val))
    else
        return :($state_var.$(term.name)[$(indices...)] = $val)
    end
end

function generate_fluent_ids(domain::Domain, state::State, term::Term,
                             sig::Signature, varmap=Dict{Var,Any}(),
                             state_var=:state)
    if get_requirements(domain)[:typing]
        ids = map(zip(term.args, sig.argtypes)) do (arg, type)
            if arg isa Var
                :(objectindex($state_var, $(QuoteNode(type)), $(varmap[arg])))
            else
                objects = sort(get_objects(domain, state, type), by=x->x.name)
                findfirst(isequal(arg), objects)
            end
        end
    else
        objects = sort(collect(get_objects(state)), by=x->x.name)
        ids = map(term.args) do arg
            if arg isa Var
                :(objectindex($state_var, $(varmap[arg])))
            else
                findfirst(isequal(arg), objects)
            end
        end
    end
    return ids
end

function generate_fluent_dims(domain::Domain, state::State, pred::Signature)
    if get_requirements(domain)[:typing]
        dims = [length(get_objects(domain, state, ty)) for ty in pred.argtypes]
    else
        dims = fill(length(get_objects(state)), length(pred.args))
    end
    return dims
end


text:
## Metsa-Oy Forest Products Supply Chain
## Edward J. Xu <edxu96@outlook.com>
## Jan 24th, 2020

function print_result_3(
    result_h_t, result_y_i, result_y_j, result_y_paper, result_z_ik,
    result_z_jk, result_z_paper_k, result_h_nmt, result_y_nmi, result_y_nmj,
    result_y_paper_nm, result_z_nmik, result_z_nmjk, result_z_paper_nmk,
    alpha_t, beta_t, a_i, b_i, e_i, c_i, r_saw, r_plywood, a_j,
    c_j, r_j, a_paper_t, b_paper_j, c_paper, r_paper, gamma_ik, delta_ik,
    gamma_jk, delta_jk, gamma_paper_k, delta_paper_k, p_fuel, sigma, omega_i,
    omega_j, omega_paper, nu_saw, nu_plywood, nu_j, nu_paper, o_saw, o_plywood,
    o_j, o_paper, i_t1, i_t2, j_t2, pi_n, rho_nm, result_x_saw_nm,
    result_x_plywood_nm, result_x_nmj, result_x_paper_nm)

  h_t = result_h_t
  y_i = result_y_i
  y_j = result_y_j
  y_paper = result_y_paper
  z_ik = result_z_ik
  z_jk = result_z_jk
  z_paper_k = result_z_paper_k
  h_nmt = result_h_nmt
  y_nmi = result_y_nmi
  y_nmj = result_y_nmj
  y_paper_nm = result_y_paper_nm
  z_nmik = result_z_nmik
  z_nmjk = result_z_nmjk
  z_paper_nmk = result_z_paper_nmk
  x_saw_nm = result_x_saw_nm
  x_plywood_nm = result_x_plywood_nm
  x_nmj = result_x_nmj
  x_paper_nm = result_x_paper_nm

  cal_dis_profit_3_first_year(h_t, y_i, y_j, y_paper, z_ik,
    z_jk, z_paper_k, h_nmt, y_nmi, y_nmj,
    y_paper_nm, z_nmik, z_nmjk, z_paper_nmk,
    alpha_t, beta_t, a_i, b_i, e_i, c_i, r_saw, r_plywood, a_j,
    c_j, r_j, a_paper_t, b_paper_j, c_paper, r_paper, gamma_ik, delta_ik,
    gamma_jk, delta_jk, gamma_paper_k, delta_paper_k, p_fuel, sigma, omega_i,
    omega_j, omega_paper, nu_saw, nu_plywood, nu_j, nu_paper, o_saw, o_plywood,
    o_j, o_paper, i_t1, i_t2, j_t2, pi_n, rho_nm, x_saw_nm,
    x_plywood_nm, x_nmj, x_paper_nm)

  cal_dis_profit_3_second_year(h_t, y_i, y_j, y_paper, z_ik,
    z_jk, z_paper_k, h_nmt, y_nmi, y_nmj,
    y_paper_nm, z_nmik, z_nmjk, z_paper_nmk,
    alpha_t, beta_t, a_i, b_i, e_i, c_i, r_saw, r_plywood, a_j,
    c_j, r_j, a_paper_t, b_paper_j, c_paper, r_paper, gamma_ik, delta_ik,
    gamma_jk, delta_jk, gamma_paper_k, delta_paper_k, p_fuel, sigma, omega_i,
    omega_j, omega_paper, nu_saw, nu_plywood, nu_j, nu_paper, o_saw, o_plywood,
    o_j, o_paper, i_t1, i_t2, j_t2, pi_n, rho_nm, x_saw_nm,
    x_plywood_nm, x_nmj, x_paper_nm)

  cal_dis_profit_3_third_year(h_t, y_i, y_j, y_paper, z_ik,
    z_jk, z_paper_k, h_nmt, y_nmi, y_nmj,
    y_paper_nm, z_nmik, z_nmjk, z_paper_nmk,
    alpha_t, beta_t, a_i, b_i, e_i, c_i, r_saw, r_plywood, a_j,
    c_j, r_j, a_paper_t, b_paper_j, c_paper, r_paper, gamma_ik, delta_ik,
    gamma_jk, delta_jk, gamma_paper_k, delta_paper_k, p_fuel, sigma, omega_i,
    omega_j, omega_paper, nu_saw, nu_plywood, nu_j, nu_paper, o_saw, o_plywood,
    o_j, o_paper, i_t1, i_t2, j_t2, pi_n, rho_nm, x_saw_nm,
    x_plywood_nm, x_nmj, x_paper_nm)
end

function cal_dis_profit_3_first_year(h_t, y_i, y_j, y_paper, z_ik,
    z_jk, z_paper_k, h_nmt, y_nmi, y_nmj,
    y_paper_nm, z_nmik, z_nmjk, z_paper_nmk,
    alpha_t, beta_t, a_i, b_i, e_i, c_i, r_saw, r_plywood, a_j,
    c_j, r_j, a_paper_t, b_paper_j, c_paper, r_paper, gamma_ik, delta_ik,
    gamma_jk, delta_jk, gamma_paper_k, delta_paper_k, p_fuel, sigma, omega_i,
    omega_j, omega_paper, nu_saw, nu_plywood, nu_j, nu_paper, o_saw, o_plywood,
    o_j, o_paper, i_t1, i_t2, j_t2, pi_n, rho_nm, x_saw_nm,
    x_plywood_nm, x_nmj, x_paper_nm)

  dis_profit =
    ## Objective fucntions regarding here-and-now decision variables.
    -sum(alpha_t[t] * h_t[t] + beta_t[t] * h_t[t]^2 for t = 1:6) +
    -sum(c_i[i] * y_i[i] for i=1:5) +
    -sum(c_j[j] * y_j[j] for j=1:2) - c_paper * y_paper +
    sum(alpha_t[t] * (
      h_t[t] - sum(a_i[i] * y_i[i] for i in i_t1[t])
      ) for t=1:3) +
    alpha_t[4] * (h_t[4] + sum(b_i[i] * y_i[i] for i in
      i_t2[1]) - sum(a_j[j] * y_j[j] for j in j_t2[1]) -
      a_paper_t[1] * y_paper) +
    alpha_t[6] * (h_t[6] + sum(b_i[i] * y_i[i] for i in
      i_t2[3]) - sum(a_j[j] * y_j[j] for j in j_t2[3]) -
      a_paper_t[3] * y_paper) +
    ## when t = 5
    alpha_t[5] * (h_t[5] + sum(b_i[i] * y_i[i] for i in
      i_t2[2]) - a_paper_t[2] * y_paper) +
    sum(e_i[i] * y_i[i] for i=1:5) * p_fuel +
    sum(gamma_ik[i, k] * z_ik[i, k] - delta_ik[i, k] * z_ik[i, k]^2
      for i=1:5, k=1:4) +
    sum(gamma_jk[j, k] * z_jk[j, k] - delta_jk[j, k] * z_jk[j, k]^2
      for j=1:2, k=1:4) +
    sum(gamma_paper_k[k] * z_paper_k[k] -
      delta_paper_k[k] * z_paper_k[k]^2
      for k=1:4) +

    ## here-and-now capacity in the second year.
    ## those in the first scenario are used because they are the same.
    - (o_saw * (x_saw_nm[1, 2] - r_saw) +
      o_plywood * (x_plywood_nm[1, 2] - r_plywood) +
      sum(o_j[j] * (x_nmj[1, 2, j] - r_j[j]) for j=1:2) +
      o_paper * (x_paper_nm[1, 2] - r_paper))

  println("dis_profit in the first year = $(dis_profit) ;")
end

function cal_dis_profit_3_second_year(h_t, y_i, y_j, y_paper, z_ik,
    z_jk, z_paper_k, h_nmt, y_nmi, y_nmj,
    y_paper_nm, z_nmik, z_nmjk, z_paper_nmk,
    alpha_t, beta_t, a_i, b_i, e_i, c_i, r_saw, r_plywood, a_j,
    c_j, r_j, a_paper_t, b_paper_j, c_paper, r_paper, gamma_ik, delta_ik,
    gamma_jk, delta_jk, gamma_paper_k, delta_paper_k, p_fuel, sigma, omega_i,
    omega_j, omega_paper, nu_saw, nu_plywood, nu_j, nu_paper, o_saw, o_plywood,
    o_j, o_paper, i_t1, i_t2, j_t2, pi_n, rho_nm, x_saw_nm,
    x_plywood_nm, x_nmj, x_paper_nm)
  m = 2

  dis_profit = sum(pi_n[n] * (sigma^(m-1) * (
    ## 1, cost of timber procurement
    -sum(alpha_t[t] * h_nmt[n, m, t] + beta_t[t] * h_nmt[n, m, t]^2
      for t = 1:6) +
    ## 2, cost of wood production
    -sum(c_i[i] * y_nmi[n, m, i] for i=1:5) +
    ## 3, cost of pulp and paper production
    -sum(c_j[j] * y_nmj[n, m, j] for j=1:2) +
    -c_paper * y_paper_nm[n, m] +
    ## 4, profit of left timbers selling
    sum(alpha_t[t] * (h_nmt[n, m, t] - sum(a_i[i] * y_nmi[n, m, i]
      for i in i_t1[t])) for t=1:3) +
    sum(alpha_t[t] * (h_nmt[n, m, t] + sum(b_i[i] * y_nmi[n, m, i]
      for i in i_t2[t-3]) - sum(a_j[j] * y_nmj[n, m, j] for j in
      j_t2[t-3]) - a_paper_t[t-3] * y_paper_nm[n, m]) for t in [4, 6]) +
    ## when t = 5
    alpha_t[5] * (h_nmt[n, m, 5] + sum(b_i[i] * y_nmi[n, m, i]
      for i in i_t2[2]) - a_paper_t[2] * y_paper_nm[n, m]) +
    ## 5, profit of fuel wood selling
    p_fuel * sum(e_i[i] * y_nmi[n, m, i] for i=1:5) +
    ## 6, profit of wood selling
    sum(omega_i[i]^(m-1) * rho_nm[n, m] * gamma_ik[i, k] *
      z_nmik[n, m, i, k] - omega_i[i]^(m-1) * rho_nm[n, m] *
      delta_ik[i, k] * z_nmik[n, m, i, k]^2 for i=1:5, k=1:4) +
    ## 7, profit of pulp selling
    sum(omega_j[j]^(m-1) * rho_nm[n, m] * gamma_jk[j, k] *
      z_nmjk[n, m, j, k] - omega_j[j]^(m-1) * rho_nm[n, m] *
      delta_jk[j, k] * z_nmjk[n, m, j, k]^2 for j=1:2, k=1:4) +
    ## 8, profit of paper selling
    sum(omega_paper^(m-1) * rho_nm[n, m] * gamma_paper_k[k] *
      z_paper_nmk[n, m, k] - omega_paper^(m-1) * rho_nm[n, m] *
      delta_paper_k[k] * z_paper_nmk[n, m, k]^2 for k=1:4)
    ) +

    ## Cost of capacity expansion in the third year
    - sigma * (
      o_saw * (x_saw_nm[1, 3] - x_saw_nm[1, 2]) +
      o_plywood * (x_plywood_nm[1, 3] - x_plywood_nm[1, 2]) +
      sum(o_j[j] * (x_nmj[1, 3, j] - x_nmj[1, 2, j]) for j=1:2) +
      o_paper * (x_paper_nm[1, 3] - x_paper_nm[1, 2]))
    ) for n=1:4)

  println("dis_profit in the second year = $(dis_profit) ;")
end

function cal_dis_profit_3_third_year(h_t, y_i, y_j, y_paper, z_ik,
    z_jk, z_paper_k, h_nmt, y_nmi, y_nmj,
    y_paper_nm, z_nmik, z_nmjk, z_paper_nmk,
    alpha_t, beta_t, a_i, b_i, e_i, c_i, r_saw, r_plywood, a_j,
    c_j, r_j, a_paper_t, b_paper_j, c_paper, r_paper, gamma_ik, delta_ik,
    gamma_jk, delta_jk, gamma_paper_k, delta_paper_k, p_fuel, sigma, omega_i,
    omega_j, omega_paper, nu_saw, nu_plywood, nu_j, nu_paper, o_saw, o_plywood,
    o_j, o_paper, i_t1, i_t2, j_t2, pi_n, rho_nm, x_saw_nm,
    x_plywood_nm, x_nmj, x_paper_nm)

  m = 3

  ## Objective fucntions regarding wait-and-see decision variables.
  dis_profit = sum(
    pi_n[n] * (
      sigma^(m-1) * (
        ## 1, cost of timber procurement
        -sum(alpha_t[t] * h_nmt[n, m, t] + beta_t[t] * h_nmt[n, m, t]^2
          for t = 1:6) +
        ## 2, cost of wood production
        -sum(c_i[i] * y_nmi[n, m, i] for i=1:5) +
        ## 3, cost of pulp and paper production
        -sum(c_j[j] * y_nmj[n, m, j] for j=1:2) +
        -c_paper * y_paper_nm[n, m] +
        ## 4, profit of left timbers selling
        sum(alpha_t[t] * (h_nmt[n, m, t] - sum(a_i[i] * y_nmi[n, m, i]
          for i in i_t1[t])) for t=1:3) +
        sum(alpha_t[t] * (h_nmt[n, m, t] + sum(b_i[i] * y_nmi[n, m, i]
          for i in i_t2[t-3]) - sum(a_j[j] * y_nmj[n, m, j] for j in
          j_t2[t-3]) - a_paper_t[t-3] * y_paper_nm[n, m]) for t in [4, 6]) +
        ## when t = 5
        alpha_t[5] * (h_nmt[n, m, 5] + sum(b_i[i] * y_nmi[n, m, i]
          for i in i_t2[2]) - a_paper_t[2] * y_paper_nm[n, m]) +
        ## 5, profit of fuel wood selling
        p_fuel * sum(e_i[i] * y_nmi[n, m, i] for i=1:5) +
        ## 6, profit of wood selling
        sum(omega_i[i]^(m-1) * rho_nm[n, m] * gamma_ik[i, k] *
          z_nmik[n, m, i, k] - omega_i[i]^(m-1) * rho_nm[n, m] *
          delta_ik[i, k] * z_nmik[n, m, i, k]^2 for i=1:5, k=1:4) +
        ## 7, profit of pulp selling
        sum(omega_j[j]^(m-1) * rho_nm[n, m] * gamma_jk[j, k] *
          z_nmjk[n, m, j, k] - omega_j[j]^(m-1) * rho_nm[n, m] *
          delta_jk[j, k] * z_nmjk[n, m, j, k]^2 for j=1:2, k=1:4) +
        ## 8, profit of paper selling
        sum(omega_paper^(m-1) * rho_nm[n, m] * gamma_paper_k[k] *
          z_paper_nmk[n, m, k] - omega_paper^(m-1) * rho_nm[n, m] *
          delta_paper_k[k] * z_paper_nmk[n, m, k]^2 for k=1:4)
        )
      ) for n=1:4
    )

  println("dis_profit in the third year = $(dis_profit) ;")
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/optimize/optimize_nonconvex.jl
# Contains the optimize! routine and subroutines needed in the branch and
# bound routine called by EAGO.
#############################################################################

function set_evaluator_flags!(d, is_post, is_intersect, is_first_eval, interval_intersect)

    d.is_post = is_post
    d.is_intersect = is_intersect
    d.is_first_eval = is_first_eval
    d.interval_intersect = interval_intersect

    return nothing
end

function reset_relaxation!(m::Optimizer)

    m._working_problem._relaxed_evaluator.is_first_eval = true
    fill!(m._working_problem._relaxed_evaluator.subexpressions_eval, false)

    m._new_eval_objective = true
    m._new_eval_constraint = true

    delete_nl_constraints!(m)
    delete_objective_cuts!(m)

    return nothing
end

"""
$(TYPEDSIGNATURES)

Creates an initial node with initial box constraints and adds it to the stack.
"""
function create_initial_node!(m::Optimizer)

    branch_variable_count = m._branch_variable_count

    variable_info = m._working_problem._variable_info
    lower_bound = zeros(Float64, branch_variable_count)
    upper_bound = zeros(Float64, branch_variable_count)
    branch_count = 1

    for i = 1:m._working_problem._variable_count
        vi = variable_info[i]
        if vi.branch_on === BRANCH
            lower_bound[branch_count] = vi.lower_bound
            upper_bound[branch_count] = vi.upper_bound
            branch_count += 1
        end
    end

    n = NodeBB(lower_bound, upper_bound, -Inf, Inf, 1, 1)
    push!(m._stack, n)
    m._node_count = 1
    m._maximum_node_id += 1

    return nothing
end

"""
$(TYPEDSIGNATURES)

Loads variables, linear constraints, and empty storage for first nlp and
quadratic cut.
"""
function load_relaxed_problem!(m::Optimizer)

    relaxed_optimizer = m.relaxed_optimizer

    # add variables and indices and constraints
    wp = m._working_problem
    branch_variable_count = 0

    variable_count = wp._variable_count
    for i = 1:variable_count

        relaxed_variable_indx = MOI.add_variable(relaxed_optimizer)
        relaxed_variable = SV(relaxed_variable_indx)
        push!(m._relaxed_variable_index, relaxed_variable_indx)

        vinfo =  wp._variable_info[i]

        is_branch_variable =  m._branch_variables[i]
        vinfo.branch_on = is_branch_variable ? BRANCH : NO_BRANCH
        is_branch_variable && (branch_variable_count += 1)

        if vinfo.is_integer

        elseif vinfo.is_fixed
            ci_sv_et = MOI.add_constraint(relaxed_optimizer, relaxed_variable, ET(vinfo.lower_bound))
            if is_branch_variable
                push!(m._relaxed_variable_eq, (ci_sv_et, branch_variable_count))
                wp._var_eq_count += 1
            end

        else
            if vinfo.has_lower_bound
                ci_sv_gt = MOI.add_constraint(relaxed_optimizer, relaxed_variable, GT(vinfo.lower_bound))
                if is_branch_variable
                    push!(m._relaxed_variable_gt, (ci_sv_gt, branch_variable_count))
                    wp._var_geq_count += 1
                end
            end

            if vinfo.has_upper_bound
                ci_sv_lt = MOI.add_constraint(relaxed_optimizer, relaxed_variable, LT(vinfo.upper_bound))
                if is_branch_variable
                    push!(m._relaxed_variable_lt, (ci_sv_lt, branch_variable_count))
                    wp._var_leq_count += 1
                end
            end
        end
    end

    # set node index to single variable constraint index maps
    m._node_to_sv_leq_ci = fill(CI{SV,LT}(-1), branch_variable_count)
    m._node_to_sv_geq_ci = fill(CI{SV,GT}(-1), branch_variable_count)
    for i = 1:wp._var_leq_count
        ci_sv_lt, branch_index = m._relaxed_variable_lt[i]
        m._node_to_sv_leq_ci[branch_index] = ci_sv_lt
    end
    for i = 1:wp._var_geq_count
        ci_sv_gt, branch_index = m._relaxed_variable_gt[i]
        m._node_to_sv_geq_ci[branch_index] = ci_sv_gt
    end

    # set number of variables to branch on
    m._branch_variable_count = branch_variable_count

    # add linear constraints
    add_linear_constraints!(m, relaxed_optimizer)

    # sets relaxed problem objective sense to Min as all problems
    # are internally converted in Min problems in EAGO
    MOI.set(relaxed_optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    return nothing
end

function presolve_global!(t::ExtensionType, m::Optimizer)

    load_relaxed_problem!(m)
    create_initial_node!(m)

    branch_variable_count = m._branch_variable_count

    m._current_xref             = fill(0.0, branch_variable_count)
    m._candidate_xref           = fill(0.0, branch_variable_count)
    m._current_objective_xref   = fill(0.0, branch_variable_count)
    m._prior_objective_xref     = fill(0.0, branch_variable_count)
    m._lower_lvd                = fill(0.0, branch_variable_count)
    m._lower_uvd                = fill(0.0, branch_variable_count)

    # populate in full space until local MOI nlp solves support constraint deletion
    # uses input model for local nlp solves... may adjust this if a convincing reason
    # to use a reformulated upper problem presents itself
    m._lower_solution      = zeros(Float64, m._working_problem._variable_count)
    m._cut_solution        = zeros(Float64, m._working_problem._variable_count)
    m._continuous_solution = zeros(Float64, m._working_problem._variable_count)
    m._upper_solution      = zeros(Float64, m._working_problem._variable_count)
    m._upper_variables     = fill(VI(-1), m._working_problem._variable_count)

    # add storage for fbbt
    m._lower_fbbt_buffer   = zeros(Float64, m._working_problem._variable_count)
    m._upper_fbbt_buffer   = zeros(Float64, m._working_problem._variable_count)

    # add storage for obbt ( perform obbt on all relaxed variables, potentially)
    m._obbt_working_lower_index = fill(false, branch_variable_count)
    m._obbt_working_upper_index = fill(false, branch_variable_count)
    m._old_low_index            = fill(false, branch_variable_count)
    m._old_upp_index            = fill(false, branch_variable_count)
    m._new_low_index            = fill(false, branch_variable_count)
    m._new_upp_index            = fill(false, branch_variable_count)
    m._lower_indx_diff          = fill(false, branch_variable_count)
    m._upper_indx_diff          = fill(false, branch_variable_count)
    m._obbt_variable_count      = branch_variable_count

    # add storage for objective cut if quadratic or nonlinear
    wp = m._working_problem
    obj_type = wp._objective_type
    if obj_type === SCALAR_QUADRATIC
        wp._objective_saf.terms = copy(wp._objective_sqf.saf.terms)
    elseif obj_type === NONLINEAR
        wp._objective_saf.terms = copy(wp._objective_nl.saf.terms)
    end

    # set subgradient refinement flag
    wp._relaxed_evaluator.is_post = m._parameters.subgrad_tighten
    wp._relaxed_evaluator.subgrad_tighten = m._parameters.subgrad_tighten
    wp._relaxed_evaluator.reverse_subgrad_tighten =  m._parameters.reverse_subgrad_tighten

    m._presolve_time = time() - m._parse_time

    return nothing
end

"""
$(SIGNATURES)

Selects node with the lowest lower bound in stack.
"""
function node_selection!(t::ExtensionType, m::Optimizer)

    m._node_count -= 1
    m._current_node = popmin!(m._stack)

    return nothing

end

"""
$(SIGNATURES)

Creates two nodes from `current_node` using information available the `x`
and stores them to the stack. By default, relative width bisection is perfomed
at a point `branch_pnt` which is a convex combination
(parameter: `branch_cvx_factor`) of the solution to the relaxation and
the midpoint of the node. If this solution lies within `branch_offset/width` of
a bound then the branch point is moved to a distance of `branch_offset/width`
from the bound.
"""
function branch_node!(t::ExtensionType, m::Optimizer)

    n = m._current_node

    lvbs = n.lower_variable_bounds
    uvbs = n.upper_variable_bounds

    max_pos = 0
    max_val = -Inf
    temp_max = 0.0

    flag = true
    for i = 1:m._branch_variable_count
        si = m._branch_to_sol_map[i]
        vi = m._working_problem._variable_info[si]
        if vi.branch_on === BRANCH
            temp_max =  uvbs[i] - lvbs[i]
            temp_max /= vi.upper_bound - vi.lower_bound
            if temp_max > max_val
                max_pos = i
                max_val = temp_max
            end
        end
    end

    lvb  = lvbs[max_pos]
    uvb  = uvbs[max_pos]
    si   = m._branch_to_sol_map[max_pos]
    lsol = m._lower_solution[si]

    cvx_f = m._parameters.branch_cvx_factor
    cvx_g = m._parameters.branch_offset

    branch_pnt = cvx_f*lsol + (1.0 - cvx_f)*(lvb + uvb)/2.0
    if branch_pnt < lvb*(1.0 - cvx_g) + cvx_g*uvb
        branch_pnt = (1.0 - cvx_g)*lvb + cvx_g*uvb
    elseif branch_pnt > cvx_g*lvb + (1.0 - cvx_g)*uvb
        branch_pnt = cvx_g*lvb + (1.0 - cvx_g)*uvb
    end

    # rounds into branch points, which in turn prevents the
    # solution at the branch point from being discarded
    N1::Interval{Float64} = Interval{Float64}(lvb, branch_pnt)
    N2::Interval{Float64} = Interval{Float64}(branch_pnt, uvb)
    lvb_1 = copy(lvbs)
    uvb_1 = copy(uvbs)
    lvb_2 = copy(lvbs)
    uvb_2 = copy(uvbs)
    lvb_1[max_pos] = N1.lo
    uvb_1[max_pos] = N1.hi
    lvb_2[max_pos] = N2.lo
    uvb_2[max_pos] = N2.hi

    lower_bound = max(n.lower_bound, m._lower_objective_value)
    upper_bound = min(n.upper_bound, m._upper_objective_value)
    new_depth = n.depth + 1

    m._maximum_node_id += 1
    X1 = NodeBB(lvb_1, uvb_1, lower_bound, upper_bound, new_depth, m._maximum_node_id)
    m._maximum_node_id += 1
    X2 = NodeBB(lvb_2, uvb_2, lower_bound, upper_bound, new_depth, m._maximum_node_id)

    push!(m._stack, X1)
    push!(m._stack, X2)

    m._node_repetitions = 1
    m._node_count += 2

    return nothing
end

"""
$(SIGNATURES)

Stores the current node to the stack after updating lower/upper bounds.
"""
function single_storage!(t::ExtensionType, m::Optimizer)
    y = m._current_node
    m._node_repetitions += 1
    m._node_count += 1
    lower_bound = max(y.lower_bound, m._lower_objective_value)
    upper_bound = min(y.upper_bound, m._upper_objective_value)
    push!(m._stack, NodeBB(y.lower_variable_bounds, y.upper_variable_bounds,
                           lower_bound, upper_bound, y.depth, y.id))

    return nothing
end

"""
$(SIGNATURES)

Selects and deletes nodes from stack with lower bounds greater than global
upper bound.
"""
function fathom!(t::ExtensionType, m::Optimizer)

    upper = m._global_upper_bound
    continue_flag = !isempty(m._stack)

    while continue_flag
        max_node = maximum(m._stack)
        max_check = (max_node.lower_bound > upper)

        if max_check
            popmax!(m._stack)
            m._node_count -= 1
            if isempty(m._stack)
                continue_flag = false
            end

        else
            if !max_check
                continue_flag = false
            elseif isempty(m._stack)
                continue_flag = false
            end

        end
    end

    return nothing
end

"""
$(SIGNATURES)

Checks to see if current node should be reprocessed.
"""
function repeat_check(t::ExtensionType, m::Optimizer)
    return false
end

relative_gap(L::Float64, U::Float64) = ((L > -Inf) && (U < Inf)) ?  abs(U - L)/(max(abs(L), abs(U))) : Inf
relative_tolerance(L::Float64, U::Float64, tol::Float64) = relative_gap(L, U)  > tol || ~(L > -Inf)

"""
$(SIGNATURES)

Checks for termination of algorithm due to satisfying absolute or relative
tolerance, infeasibility, or a specified limit, returns a boolean valued true
if algorithm should continue.
"""
function termination_check(t::ExtensionType, m::Optimizer)

    node_in_stack = length(m._stack)
    L = m._global_lower_bound
    U = m._global_upper_bound

    if node_in_stack === 0

        if m._first_solution_node > 0
            m._termination_status_code = MOI.OPTIMAL
            m._result_status_code = MOI.FEASIBLE_POINT
            (m._parameters.verbosity >= 3) && println("Empty Stack: Exhaustive Search Finished")

        else
            m._termination_status_code = MOI.INFEASIBLE
            m._result_status_code = MOI.INFEASIBILITY_CERTIFICATE
            (m._parameters.verbosity >= 3) && println("Empty Stack: Infeasible")
        end

    elseif node_in_stack >= m._parameters.node_limit

        m._termination_status_code = MOI.NODE_LIMIT
        m._result_status_code = MOI.UNKNOWN_RESULT_STATUS
        (m._parameters.verbosity >= 3) && println("Node Limit Exceeded")

    elseif m._iteration_count >= m._parameters.iteration_limit

        m._termination_status_code = MOI.ITERATION_LIMIT
        m._result_status_code = MOI.UNKNOWN_RESULT_STATUS
        (m._parameters.verbosity >= 3) && println("Maximum Iteration Exceeded")

    elseif ~relative_tolerance(L, U, m._parameters.relative_tolerance)

        m._termination_status_code = MOI.OPTIMAL
        m._result_status_code = MOI.FEASIBLE_POINT
        (m._parameters.verbosity >= 3) && println("Relative Tolerance Achieved")

    elseif (U - L) < m._parameters.absolute_tolerance

        m._termination_status_code = MOI.OPTIMAL
        m._result_status_code = MOI.FEASIBLE_POINT
        (m._parameters.verbosity >= 3) && println("Absolute Tolerance Achieved")

    elseif m._run_time > m._parameters.time_limit

        m._termination_status_code = MOI.TIME_LIMIT
        m._result_status_code = MOI.UNKNOWN_RESULT_STATUS
        (m._parameters.verbosity >= 3) && println("Time Limit Exceeded")

    else

        return false

    end

    return true
end

"""
$(SIGNATURES)

Checks for convergence of algorithm with respect to absolute and/or relative
tolerances.
"""
function convergence_check(t::ExtensionType, m::Optimizer)

  L = m._lower_objective_value
  U = m._global_upper_bound
  t = (U - L) <= m._parameters.absolute_tolerance
  if (U < Inf) && (L > Inf)
      t |= (abs(U - L)/(max(abs(L), abs(U))) <= m._parameters.relative_tolerance)
  end

  if t && m._min_converged_value < Inf
      m._min_converged_value = min(m._min_converged_value, L)
  else
      m._min_converged_value = L
  end

  return t
end

"""
$(SIGNATURES)

Takes an `MOI.TerminationStatusCode` and a `MOI.ResultStatusCode` and returns
the tuple `(valid_result::Bool, feasible::Bool)`. The value `valid_result` is
`true` if the pair of codes prove that either the subproblem solution was solved
to global optimality or the subproblem solution is infeasible. The value of
`feasible` is true if the problem is feasible and false if the problem is infeasible.
"""
function is_globally_optimal(t::MOI.TerminationStatusCode, r::MOI.ResultStatusCode)

    feasible = false
    valid_result = false

    if (t === MOI.INFEASIBLE && r === MOI.INFEASIBILITY_CERTIFICATE)
        valid_result = true

    elseif (t === MOI.INFEASIBLE && r === MOI.NO_SOLUTION)
        valid_result = true

    elseif (t === MOI.INFEASIBLE && r === MOI.UNKNOWN_RESULT_STATUS)
        valid_result = true

    elseif (t === MOI.OPTIMAL && r === MOI.FEASIBLE_POINT)
        valid_result = true
        feasible = true

    elseif (t === MOI.INFEASIBLE_OR_UNBOUNDED && r === MOI.NO_SOLUTION)
        valid_result = true
        feasible = false

    end

    return valid_result, feasible
end

"""
$(SIGNATURES)

Takes an `MOI.TerminationStatusCode` and a `MOI.ResultStatusCode` and returns `true`
if this corresponds to a solution that is proven to be feasible.
Returns `false` otherwise.
"""
function is_feasible_solution(t::MOI.TerminationStatusCode, r::MOI.ResultStatusCode)

    termination_flag = false
    result_flag = false

    (t === MOI.OPTIMAL) && (termination_flag = true)
    (t === MOI.LOCALLY_SOLVED) && (termination_flag = true)

    # This is default solver specific... the acceptable constraint tolerances
    # are set to the same values as the basic tolerance. As a result, an
    # acceptably solved solution is feasible but non necessarily optimal
    # so it should be treated as a feasible point
    if (t === MOI.ALMOST_LOCALLY_SOLVED) && (r === MOI.NEARLY_FEASIBLE_POINT)
        termination_flag = true
        result_flag = true
    end

    (r === MOI.FEASIBLE_POINT) && (result_flag = true)

    return (termination_flag && result_flag)
end

"""
$(SIGNATURES)

Retrieves the lower and upper duals for variable bounds from the
`relaxed_optimizer` and sets the appropriate values in the
`_lower_lvd` and `_lower_uvd` storage fields.
"""
function set_dual!(m::Optimizer)

    relaxed_optimizer = m.relaxed_optimizer
    relaxed_variable_lt = m._relaxed_variable_lt
    relaxed_variable_gt = m._relaxed_variable_gt

    for i = 1:m._working_problem._var_leq_count
        ci_lt, i_lt = @inbounds relaxed_variable_lt[i]
        @inbounds m._lower_uvd[i_lt] = MOI.get(relaxed_optimizer, MOI.ConstraintDual(), ci_lt)
    end
    for i = 1:m._working_problem._var_geq_count
        ci_gt, i_gt = @inbounds relaxed_variable_gt[i]
        @inbounds m._lower_lvd[i_gt] = MOI.get(relaxed_optimizer, MOI.ConstraintDual(), ci_gt)
    end

    return nothing
end

"""
$(SIGNATURES)

Runs interval, linear, quadratic contractor methods followed by obbt and a
constraint programming walk up to tolerances specified in
`EAGO.Optimizer` object.
"""
function preprocess!(t::ExtensionType, m::Optimizer)

    reset_relaxation!(m)

    wp = m._working_problem
    params = m._parameters

    # Sets initial feasibility
    feasible_flag = true
    m._obbt_performed_flag = false

    # compute initial volume
    m._initial_volume = prod(upper_variable_bounds(m._current_node) -
                             lower_variable_bounds(m._current_node))

    if params.fbbt_lp_depth >= m._iteration_count
        load_fbbt_buffer!(m)
        for i = 1:m._parameters.fbbt_lp_repetitions
            if feasible_flag
                for j = 1:wp._saf_leq_count
                    !feasible_flag && break
                    saf_leq =  wp._saf_leq[j]
                    feasible_flag &= fbbt!(m, saf_leq)
                end
                !feasible_flag && break

                for j = 1:wp._saf_eq_count
                    !feasible_flag && break
                    saf_eq = wp._saf_eq[j]
                    feasible_flag &= fbbt!(m, saf_eq)
                end
                !feasible_flag && break
            end
        end
        unpack_fbbt_buffer!(m)
    end

    # done after cp to prevent using cp specific flags in cut generation
    set_first_relax_point!(m)

    cp_walk_count = 0
    perform_cp_walk_flag = feasible_flag
    perform_cp_walk_flag &= (params.cp_depth >= m._iteration_count)
    perform_cp_walk_flag &= (cp_walk_count < m._parameters.cp_repetitions)
    while perform_cp_walk_flag
        feasible_flag &= set_constraint_propagation_fbbt!(m)
        !feasible_flag && break
        cp_walk_count += 1
        perform_cp_walk_flag = (cp_walk_count < m._parameters.cp_repetitions)
    end

    obbt_count = 0
    perform_obbt_flag = feasible_flag
    perform_obbt_flag &= (params.obbt_depth >= m._iteration_count)
    perform_obbt_flag &= (obbt_count < m._parameters.obbt_repetitions)

    while perform_obbt_flag
        feasible_flag &= obbt!(m)
        m._obbt_performed_flag = true
        !feasible_flag && break
        obbt_count += 1
        perform_obbt_flag     = (obbt_count < m._parameters.obbt_repetitions)
    end

    m._final_volume = prod(upper_variable_bounds(m._current_node) -
                           lower_variable_bounds(m._current_node))

    m._preprocess_feasibility = feasible_flag

    return nothing
end

"""
$(SIGNATURES)

Updates the relaxed constraint by setting the constraint set of `v == x*`` ,
`xL_i <= x_i`, and `x_i <= xU_i` for each such constraint added to the relaxed
optimizer.
"""
function update_relaxed_problem_box!(m::Optimizer)

    opt = m.relaxed_optimizer
    wp = m._working_problem

    n = m._current_node
    lower_bound = n.lower_variable_bounds
    upper_bound = n.upper_variable_bounds

    relaxed_variable_eq = m._relaxed_variable_eq
    for i = 1:wp._var_eq_count
        constr_indx, node_indx =  relaxed_variable_eq[i]
        MOI.set(opt, MOI.ConstraintSet(), constr_indx, ET( lower_bound[node_indx]))
    end

    relaxed_variable_lt = m._relaxed_variable_lt
    for i = 1:wp._var_leq_count
        constr_indx, node_indx =  relaxed_variable_lt[i]
        MOI.set(opt, MOI.ConstraintSet(), constr_indx, LT( upper_bound[node_indx]))
    end

    relaxed_variable_gt = m._relaxed_variable_gt
    for i = 1:wp._var_geq_count
        constr_indx, node_indx =  relaxed_variable_gt[i]
        MOI.set(opt, MOI.ConstraintSet(), constr_indx, GT( lower_bound[node_indx]))
    end

    return nothing
end

function interval_objective_bound(m::Optimizer, n::NodeBB)

    interval_objective_bound = bound_objective(m)

    if interval_objective_bound > m._lower_objective_value
        m._lower_objective_value = interval_objective_bound
        fill!(m._lower_lvd, 0.0)
        fill!(m._lower_uvd, 0.0)
        m._cut_add_flag = false
        return true

    end

    return false
end

"""
$(SIGNATURES)

A fallback lower bounding problem that consists of an natural interval extension
calculation. This is called when the optimizer used to compute the lower bound
does not return a termination and primal status code indicating that it
successfully solved the relaxation to a globally optimal point.
"""
function fallback_interval_lower_bound!(m::Optimizer, n::NodeBB)

    feasible_flag = true

    if !cp_condition(m)
        for i = 1:m._working_problem._saf_leq_count
            saf_leq =  m._working_problem._saf_leq[i]
            feasible_flag &= (lower_interval_bound(m, saf_leq, n) <= 0.0)
            !feasible_flag && break
        end

        if feasible_flag
            for i = 1:m._working_problem._saf_eq_count
                saf_eq =  m._working_problem._saf_eq[i]
                lower_value, upper_value = interval_bound(m, saf_eq, n)
                feasible_flag &= (lower_value <= 0.0 <= upper_value)
                !feasible_flag && break
            end
        end

        if feasible_flag
            for i = 1:m._working_problem._sqf_leq_count
                sqf_leq =  m._working_problem._sqf_leq[i]
                feasible_flag &= (lower_interval_bound(m, sqf_leq, n) <= 0.0)
                !feasible_flag && break
            end
        end

        if feasible_flag
            for i = 1:m._working_problem._sqf_eq_count
                sqf_eq =  m._working_problem._sqf_eq[i]
                lower_value, upper_value = interval_bound(m, sqf_eq, n)
                feasible_flag &= (lower_value <= 0.0 <= upper_value)
                !feasible_flag && break
            end
        end

        if feasible_flag
            for i = 1:m._working_problem._nonlinear_count
                nl_constr =  m._working_problem._nonlinear_constr[i]
                lower_value, upper_value = interval_bound(m, nl_constr, n)
                feasible_flag &= upper_value < nl_constr.lower_bound
                feasible_flag &= lower_value > nl_constr.upper_bound
                !feasible_flag && break
            end
        end
    end

    if feasible_flag
        interval_objective_used = interval_objective_bound(m, n)
        @__dot__ m._current_xref = 0.5*(n.upper_variable_bounds + n.lower_variable_bounds)
        unsafe_check_fill!(isnan, m._current_xref, 0.0, length(m._current_xref))
    else
        m._lower_objective_value = -Inf
    end
    m._lower_feasibility = feasible_flag

    return nothing
end

"""
$(SIGNATURES)

Constructs and solves the relaxation using the default EAGO relaxation scheme
and optimizer on node `y`.
"""
function lower_problem!(t::ExtensionType, m::Optimizer)

    n = m._current_node

    m._working_problem._relaxed_evaluator.is_post = m._parameters.subgrad_tighten
    if !m._obbt_performed_flag
        if m._nonlinear_evaluator_created
            set_node!(m._working_problem._relaxed_evaluator, n)
            set_node_flag!(m)
            set_reference_point!(m)
            fill!(m._working_problem._relaxed_evaluator.subexpressions_eval, false)
        end
        update_relaxed_problem_box!(m)
    end
    m._working_problem._objective_nl.has_value = false
    m._working_problem._relaxed_evaluator.interval_intersect = false

    if !m._obbt_performed_flag
        relax_constraints!(m, 1)
    end
    relax_objective!(m, 1)

    # Optimizes the object
    relaxed_optimizer = m.relaxed_optimizer

    MOI.set(relaxed_optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    MOI.optimize!(relaxed_optimizer)

    m._lower_termination_status = MOI.get(relaxed_optimizer, MOI.TerminationStatus())
    m._lower_result_status = MOI.get(relaxed_optimizer, MOI.PrimalStatus())
    valid_flag, feasible_flag = is_globally_optimal(m._lower_termination_status, m._lower_result_status)

    if valid_flag && feasible_flag
        set_dual!(m)
        m._cut_add_flag = true
        m._lower_feasibility = true
        m._lower_objective_value = MOI.get(relaxed_optimizer, MOI.ObjectiveValue())
        for i = 1:m._working_problem._variable_count
             m._lower_solution[i] = MOI.get(relaxed_optimizer, MOI.VariablePrimal(), m._relaxed_variable_index[i])
        end

    elseif valid_flag
        m._cut_add_flag = false
        m._lower_feasibility  = false
        m._lower_objective_value = -Inf

    else
        fallback_interval_lower_bound!(m, n)
    end

    return nothing
end

"""
$(SIGNATURES)

Updates the internal storage in the optimizer after a valid feasible cut is added.
"""
function cut_update!(m::Optimizer)

    m._cut_feasibility = true

    relaxed_optimizer = m.relaxed_optimizer
    obj_val = MOI.get(relaxed_optimizer, MOI.ObjectiveValue())
    prior_obj_val = (m._cut_iterations == 2) ? m._lower_objective_value : m._cut_objective_value

    m._cut_add_flag = true
    m._lower_termination_status = m._cut_termination_status
    m._lower_result_status = m._cut_result_status
    m._cut_solution[:] = MOI.get(relaxed_optimizer, MOI.VariablePrimal(), m._relaxed_variable_index)

    if prior_obj_val < obj_val
        m._cut_objective_value = obj_val
        m._lower_objective_value = obj_val
        set_dual!(m)
        copyto!(m._lower_solution, m._cut_solution)

    else
        m._cut_objective_value = prior_obj_val
        m._lower_objective_value = prior_obj_val
        m._cut_add_flag = false
    end

    return nothing
end

"""
$(SIGNATURES)

Checks if a cut should be added and computes a new reference point to add the
cut at. If no cut should be added the constraints not modified in place are
deleted from the relaxed optimizer and the solution is compared with the
interval lower bound. The best lower bound is then used.
"""
function cut_condition(t::ExtensionType, m::Optimizer)

    # always add cut if below the minimum iteration limit, otherwise add cut
    # the number of cuts is less than the maximum and the distance between
    # prior solutions exceeded a tolerance.
    continue_cut_flag = m._cut_add_flag
    continue_cut_flag &= (m._cut_iterations < m._parameters.cut_max_iterations)

    # compute distance between prior solutions and compare to tolerances
    n = m._current_node
    ns_indx = m._branch_to_sol_map

    cvx_factor =  m._parameters.cut_cvx
    xsol = (m._cut_iterations > 1) ? m._cut_solution[ns_indx] : m._lower_solution[ns_indx]
    xnew = (1.0 - cvx_factor)*mid(n) + cvx_factor*xsol

    continue_cut_flag &= (norm((xsol - xnew)/diam(n), 1) > m._parameters.cut_tolerance)
    continue_cut_flag |= (m._cut_iterations < m._parameters.cut_min_iterations)

    # update reference point for new cut
    if continue_cut_flag
        copyto!(m._current_xref, xnew)
        if m._nonlinear_evaluator_created
            set_reference_point!(m)
            fill!(m._working_problem._relaxed_evaluator.subexpressions_eval, false)
        end
    end

    # check to see if interval bound is preferable and replaces the objective
    # value with the interval value if so. Any available dual values are then
    # set to zero since the interval bounds are by definition constant
    if m._lower_feasibility && !continue_cut_flag
        objective_lo = -Inf
        obj_type = m._working_problem._objective_type
        if obj_type === SINGLE_VARIABLE
            var_index = m._working_problem._objective_sv.variable.value
            if m._branch_variables[var_index]
                obj_indx = m._sol_to_branch_map[var_index]
                lower_variable_bnd = n.lower_variable_bounds[obj_indx]
                if !isinf(lower_variable_bnd)
                    objective_lo = lower_variable_bnd
                end
            end

        elseif obj_type === SCALAR_AFFINE
            objective_lo = lower_interval_bound(m, m._working_problem._objective_saf_parsed, n)

        elseif obj_type === SCALAR_QUADRATIC
            objective_lo = lower_interval_bound(m, m._working_problem._objective_sqf, n)

        elseif obj_type === NONLINEAR
            objective_lo = lower_interval_bound(m, m._working_problem._objective_nl, n)

        end

        if objective_lo > m._lower_objective_value
            m._lower_objective_value = objective_lo
            fill!(m._lower_lvd, 0.0)
            fill!(m._lower_uvd, 0.0)
        end
    end

    m._cut_iterations += 1

    return continue_cut_flag
end

"""
$(SIGNATURES)

Adds a cut for each constraint and the objective function to the subproblem.
"""
function add_cut!(t::ExtensionType, m::Optimizer)

    fill!(m._working_problem._relaxed_evaluator.subexpressions_eval, false)
    m._working_problem._relaxed_evaluator.is_first_eval = true
    m._working_problem._relaxed_evaluator.is_intersect = false
    m._new_eval_objective = true
    m._new_eval_constraint = true

    relax_constraints!(m, m._cut_iterations)
    relax_objective!(m, m._cut_iterations)

    # Optimizes the object
    relaxed_optimizer = m.relaxed_optimizer
    MOI.optimize!(relaxed_optimizer)

    m._cut_termination_status = MOI.get(relaxed_optimizer, MOI.TerminationStatus())
    m._cut_result_status = MOI.get(relaxed_optimizer, MOI.PrimalStatus())
    valid_flag, feasible_flag = is_globally_optimal(m._cut_termination_status, m._cut_result_status)

    if valid_flag && feasible_flag
        cut_update!(m)

    elseif valid_flag
        m._cut_add_flag = false
        m._lower_feasibility  = false
        m._lower_objective_value = -Inf

    else
        m._cut_add_flag = false
    end

    return nothing
end

"""
$(SIGNATURES)

Default check to see if the upper bounding problem should be run. By default,
The upper bounding problem is run on every node up to depth `upper_bounding_depth`
and is triggered with a probability of `0.5^(depth - upper_bounding_depth)`
afterwards.
"""
function default_nlp_heurestic(m::Optimizer)
    bool = false
    ubd_limit = m._parameters.upper_bounding_depth
    depth = m._current_node.depth
    bool |= (depth <= ubd_limit)
    bool |= (rand() < 0.5^(depth - m._parameters.upper_bounding_depth))
    return bool
end

"""
$(SIGNATURES)

Default upper bounding problem which simply calls `solve_local_nlp!` to solve
the nlp locally.
"""
function upper_problem!(t::ExtensionType, m::Optimizer)

    if !default_nlp_heurestic(m)
        m._upper_feasibility = false
        m._upper_objective_value = Inf

    else
        solve_local_nlp!(m)

    end

    return nothing
end

"""
$(SIGNATURES)

Default postprocess perfoms duality-based bound tightening on the `y`.
"""
function postprocess!(t::ExtensionType, m::Optimizer)

    if m._parameters.dbbt_depth > m._iteration_count
        variable_dbbt!(m._current_node, m._lower_lvd, m._lower_uvd,
                       m._lower_objective_value, m._global_upper_bound,
                       m._branch_variable_count)
    end

    return nothing
end

"""
$(SIGNATURES)

Provides a hook for extensions to EAGO as opposed to standard global, local,
or linear solvers.
"""
optimize_hook!(t::ExtensionType, m::Optimizer) = nothing

function store_candidate_solution!(m::Optimizer)

    if m._upper_feasibility && (m._upper_objective_value < m._global_upper_bound)

        m._feasible_solution_found = true
        m._first_solution_node = m._maximum_node_id
        m._solution_value = m._upper_objective_value
        m._global_upper_bound = m._upper_objective_value
        @__dot__ m._continuous_solution = m._upper_solution

    end
    return nothing
end

function set_global_lower_bound!(m::Optimizer)

    if !isempty(m._stack)

        min_node = minimum(m._stack)
        lower_bound = min_node.lower_bound
        if m._global_lower_bound < lower_bound
            m._global_lower_bound = lower_bound
        end

    end

    return nothing
end

# wraps subroutine call to isolate ExtensionType
parse_global!(m::Optimizer) = parse_global!(m.ext_type, m)
presolve_global!(m::Optimizer) = presolve_global!(m.ext_type, m)
termination_check(m::Optimizer) = termination_check(m.ext_type, m)
cut_condition(m::Optimizer) = cut_condition(m.ext_type, m)
convergence_check(m::Optimizer) = convergence_check(m.ext_type, m)
repeat_check(m::Optimizer) = repeat_check(m.ext_type, m)
node_selection!(m::Optimizer) = node_selection!(m.ext_type, m)
preprocess!(m::Optimizer) = preprocess!(m.ext_type, m)
lower_problem!(m::Optimizer) = lower_problem!(m.ext_type, m)
add_cut!(m::Optimizer) = add_cut!(m.ext_type, m)
upper_problem!(m::Optimizer) = upper_problem!(m.ext_type, m)
postprocess!(m::Optimizer) = postprocess!(m.ext_type, m)
single_storage!(m::Optimizer) = single_storage!(m.ext_type, m)
branch_node!(m::Optimizer) = branch_node!(m.ext_type, m)
fathom!(m::Optimizer) = fathom!(m.ext_type, m)
revert_adjusted_upper_bound!(m::Optimizer) = revert_adjusted_upper_bound!(m.ext_type, m)

"""
$(TYPEDSIGNATURES)

Solves the branch and bound problem with the input EAGO optimizer object.
"""
function global_solve!(m::Optimizer)

    m._iteration_count = 1
    m._node_count = 1

    parse_global!(m)
    presolve_global!(m)

    logging_on = m._parameters.log_on
    verbosity = m._parameters.verbosity

    # terminates when max nodes or iteration is reach, or when node stack is empty
    while !termination_check(m)

        # Selects node, deletes it from stack, prints based on verbosity
        node_selection!(m)
        (verbosity >= 3) && print_node!(m)

        # Performs prepocessing and times
        logging_on && (start_time = time())
        preprocess!(m)
        if logging_on
            m._last_preprocess_time = time() - start_time
        end

        if m._preprocess_feasibility

            # solves & times lower bounding problem
            logging_on && (start_time = time())
            m._cut_iterations = 1
            lower_problem!(m)
            while cut_condition(m)
                add_cut!(m)
            end
            if logging_on
                m._last_lower_problem_time = time() - start_time
            end
            print_results!(m, true)
            print_results_post_cut!(m)

            # checks for infeasibility stores solution
            if m._lower_feasibility
                if !convergence_check(m)

                    logging_on && (start_time = time())
                    upper_problem!(m)
                    if logging_on
                        m._last_upper_problem_time = time() - start_time
                    end
                    print_results!(m, false)
                    store_candidate_solution!(m)
                    if m._input_problem._optimization_sense === MOI.FEASIBILITY_SENSE
                        if !m.feasible_local_continue || m.local_solve_only
                            break
                        end
                    end

                    # Performs and times post processing
                    logging_on && (start_time = time())
                    postprocess!(m)
                    if logging_on
                        m._last_postprocessing_time = time() - start_time
                    end

                    # Checks to see if the node
                    if m._postprocess_feasibility
                        if repeat_check(m)
                            single_storage!(m)
                        else
                            branch_node!(m)
                        end
                    end
                end
            end
            fathom!(m)
        else
            m._lower_objective_value = -Inf
            m._lower_feasibility = false
            m._upper_feasibility = false
        end
        set_global_lower_bound!(m)
        m._run_time = time() - m._start_time
        m._time_left = m._parameters.time_limit - m._run_time
        log_iteration!(m)
        print_iteration!(m)
        m._iteration_count += 1
    end

    revert_adjusted_upper_bound!(m)
    m._objective_value = m._global_upper_bound

    # Prints the solution
    print_solution!(m)

    return nothing
end

optimize!(::Val{MINCVX}, m::Optimizer) = global_solve!(m)


text:
using sem, Test, SafeTestsets

@time @safetestset "Examples" begin include("examples/example_tests.jl") end


text:
module CuBatchedRoutines

using GPUArrays, CuArrays
using Reexport
@reexport using BatchedRoutines
import BatchedRoutines: batched_gemm!, batched_gemm, batched_tr!, batched_tr

include("blas.jl")
include("linalg.jl")

end # module


text:
function build_application(m, options)
    isnothing(options.application) && return
    build_dir = application_build_dir(m, options)
    ispath(build_dir) || mkpath(build_dir)
    @info "application options: " options.application build_dir

    exec_file = String[pkgdir(m, x) for x in options.application.precompile.execution_file]
    stmt_file = String[pkgdir(m, x) for x in options.application.precompile.statements_file]

    create_app(
        pkgdir(m),
        build_dir;
        executables = [options.name => "julia_main"],
        precompile_execution_file = exec_file,
        precompile_statements_file = stmt_file,
        incremental = options.application.incremental,
        filter_stdlibs = options.application.filter_stdlibs,
        force = true,
        include_lazy_artifacts = options.application.include_lazy_artifacts,
        cpu_target = options.application.cpu_target,
        c_driver_program = options.application.c_driver_program,
    )

    # build completions
    build_application_completion(m, options)

    # bundle assets
    bundle_assets(m, options)
    return
end

function build_application_completion(m::Module, options::Configs.Comonicon)
    isempty(options.application.shell_completions) && return
    build_dir = application_build_dir(m, options)
    completion_dir = joinpath(build_dir, "completions")
    if !ispath(completion_dir)
        @info "creating path: $completion_dir"
        mkpath(completion_dir)
    end

    for sh in options.application.shell_completions
        script = completion_script(m, options, sh)
        write(joinpath(completion_dir, "$sh.completion"), script)
    end
    return
end

function bundle_assets(m::Module, options::Configs.Comonicon)
    isempty(options.application.assets) && return
    build_dir = application_build_dir(m, options)
    share_dir = joinpath(build_dir, "share")
    ispath(share_dir) || mkpath(share_dir)
    for asset in options.application.assets
        if isnothing(asset.package)
            dst = joinpath(build_dir, "share", asset.path)
        else
            dst = joinpath(build_dir, "share", asset.package, asset.path)
        end
        ispath(dirname(dst)) || mkpath(dirname(dst))
        cp(Configs.get_path(m, asset), dst; force = true, follow_symlinks = true)
    end
    return
end

function build_application_tarball(m, options)
    tarball = joinpath(pwd(), tarball_name(m, options.name, "application"))
    cd(options.application.path) do
        run(`tar -czvf $tarball $(options.name)`)
    end
end

function application_build_dir(m, options)
    return abspath(options.application.path, options.name)
end


text:
using Base.Test
using Lens

function f()
  x = 10;
  lens(:mylens,x=x)
  y = 20
  lens(:his,y)
  z = 30
  lens(:yourlens,z=z,y=y)
  40
end

value, results = capture(f, [(:mylens,[:x])])
@test value == 40
@test get(results)[1] == 10
value, results = capture(f, :his)
@test get(results)[1] == 20

value, results = capture(f, [(:yourlens, [:z,:y])])
@test get(results;capturename=:z)[1] == 30
@test get(results;capturename=:y)[1] == 20

## Parallel Tests
addprocs(1)
@everywhere using Lens

function parbench()
  a = zeros(100)
  tot = 0
  lens(:sum,tot)
  for proc in procs()
    Bref = @spawn (s = sum(rand(1000,1000)^2); lens(:sum,s); s)
    tot += fetch(Bref)
  end
  tot
end

## Test Macro
@capture((x = 3 ;x + lens(:tlens, x^x)),:tlens)


text:
module Stats
    import Base.quantile

    export

    # scalar_stats
    gmean,
    hmean,
    variation,
    sem,
    mad,
    skewness,
    kurtosis,
    minmax,
    midrange,
    range,
    tiedrank,
    percentile,
    quantile,
    quartile,
    quintile,
    decile,
    iqr,
    table,
    mode,
    modes,
    describe,

    # weighted stats
    wmean,

    # corr
    autocor,
    cor_spearman,
    cor_kendall,

    # others
    ## Types
    StatisticalModel,
    RegressionModel,

    ## Fcuntions
    coef,
    coeftable,
    confint,
    ecdf,
    findat,
    inverse_rle,
    loglikelihood,
    nobs,
    predict,
    residuals,
    model_response,
    rle,
    stderr,
    vcov

    include("scalar_stats.jl")
    include("weighted_stats.jl")
 #   include("corr.jl")
    include("others.jl")

end # module


text:
export GP

"""
    GP{Tm<:MeanFunction, Tk<:Kernel}

A Gaussian Process (GP) with known mean `m` and kernel `k`, coordinated by `gpc`.
"""
struct GP{Tm<:MeanFunction, Tk<:Kernel} <: AbstractGP
    m::Tm
    k::Tk
    n::Int
    gpc::GPC
    function GP{Tm, Tk}(m::Tm, k::Tk, gpc::GPC) where {Tm, Tk}
        gp = new{Tm, Tk}(m, k, next_index(gpc), gpc)
        gpc.n += 1
        return gp
    end
end
GP(m::Tm, k::Tk, gpc::GPC) where {Tm<:MeanFunction, Tk<:Kernel} = GP{Tm, Tk}(m, k, gpc)

GP(f, k::Kernel, gpc::GPC) = GP(CustomMean(f), k, gpc)
GP(m::Real, k::Kernel, gpc::GPC) = GP(ConstMean(m), k, gpc)
GP(k::Kernel, gpc::GPC) = GP(ZeroMean(), k, gpc)
GP(k::Kernel, m, gpc::GPC) = GP(m, k, gpc)

mean_vector(f::GP, x::AV) = ew(f.m, x)

cov(f::GP, x::AV) = pw(f.k, x)
cov_diag(f::GP, x::AV) = ew(f.k, x)

cov(f::GP, x::AV, x′::AV) = pw(f.k, x, x′)
cov_diag(f::GP, x::AV, x′::AV) = ew(f.k, x, x′)

function cov(f::GP, f′::GP, x::AV, x′::AV)
    return f === f′ ? cov(f, x, x′) : zeros(length(x), length(x′))
end
function cov_diag(f::GP, f′::GP, x::AV, x′::AV)
    return f === f′ ? cov_diag(f, x, x′) : zeros(length(x))
end


text:
using LinearAlgebra, Statistics, Roots, Optim, Plots, Random; Random.seed!(1234)
n = 1000 # number of samples to draw

# the true regression function
r(x) = 2 + 1/50*x*(30-x)
# the true density function
σy = 3/2
f(x,y) = 3/4000 * 1/√(2π*σy^2) * x*(20-x)*exp(-1/(2σy^2)*(y-r(x))^2)

# x values and y values for a grid
xs = 0:1/2^3:20
ys = 0:1/2^3:12

# F(t) is the CDF of X
F(t) = -t^3/4000 + 3t^2/400
# Inverse CDF sampling
function sampleX()
    U = rand()
    find_zero(t->F(t)-U,(0,20),Bisection())
end

function sampleXY(r,σ)
    X = sampleX()
    Y = r(X) + σ*randn()
    (X,Y)
end

samples = [sampleXY(r,σy) for i=1:n]

D(u) = abs(u) < 1 ? 70/81*(1-abs(u)^3)^3 : 0 # tri-cube function
D(λ,u) = 1/λ*D(u/λ) # scaled tri-cube
K(λ,x,y) = D(λ,x) * D(λ,y) # kernel
kde(λ,x,y,samples) = sum(K(λ,x-Xi,y-Yi) for (Xi,Yi) in samples)/length(samples)

# `optimize` takes functions which accept vector arguments, so we treat
# λ as a one-entry vector
L(λ) = sum((f(x,y) - kde(λ,x,y,samples))^2 for x=xs,y=ys)*step(xs)*step(ys)

# minimize L using the BFGS method
λ_best = optimize(λ->L(first(λ)),[1.0],BFGS())

using AsyPlots
#plot = Plots.plot
#heatmap = Plots.heatmap
    # background colors
softseagreen = NamedColor("SoftSeaGreen",(0.96, 0.995, 0.98))
softblue = NamedColor("SoftBlue",(0.92,0.95,0.99))
softyellow = NamedColor("SoftYellow",(0.98, 0.98, 0.9))

densities = [f(x,y) for x=xs,y=ys]
xmax = xs[end]
ymax = ys[end]

# plot accoutrements
extras = Plot([Path([0 10; -0.5 10]),Label("10",(-1.5,10))]) +
    Plot([Path([20 0; 20 -0.66]),Label("20",(20,-1.5))]) +
    Plot(Path([0 0; 0 ymax],arrow=Arrow())) +
    Plot(Path([0 0; 1.05*xmax 0],arrow=Arrow())) +
    Plot(Label("hours studied",(10,-1))) +
    Plot(Label("score",(-1,5);rotation=90);
         width=256,height=128,ignoreaspect=true,bgfill=false)
function sampleXY()
    X = sampleX()
    Y = r(X) + σy*randn()
    if !(0 < Y < 10)
        sampleXY()
    else
        (X,Y)
    end
end
Random.seed!(1234)
n = 1000
samples = [sampleXY() for i=1:n]
samplesfig = Plot([Point(point;color="MidnightBlue",linewidth=2)
                   for point in samples]) + extras
save("figures/exam-samples.pdf",samplesfig)

densityfig = AsyPlots.heatmap(densities,(0,0),(xmax,ymax);
                  colors=[NamedColor("white"),NamedColor("MidnightBlue")],
              smooth=true,alpha=fill(4/5,size(densities)...)) + extras

save("exam-density.pdf",densityfig)

densitywithline = AsyPlots.heatmap(densities,(0,0),(xmax,ymax);
                                        colors=[NamedColor("white"),NamedColor("MidnightBlue")],
                                        smooth=true,alpha=fill(4/5,size(densities)...)) +
                   Plot(Path([(x,r(x)) for x=0:1/2^4:20],color="DarkRed",linewidth=2)) + extras

save("exam-density-line.pdf",densitywithline)

cd(@__DIR__)

pwd()

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


text:
@testset "Sitemap gen" begin
    f = joinpath(p, "basic", "__site", "sitemap.xml")
    @test isfile(f)
    fc = prod(readlines(f, keep=true))

    @test occursin(raw"""
        <?xml version="1.0" encoding="utf-8" standalone="yes" ?>
        <urlset xmlns="http://www.sitemaps.org/schemas/sitemap/0.9">""", fc)
    # check pages
    for pg in ("", "menu1", "menu2", "menu3")
        pgg = joinpath(pg, "index.html")
        @test occursin("""
            <loc>https://tlienart.github.io/FranklinTemplates.jl/$pgg</loc>""", fc)
    end
end

@testset "Robots.txt gen" begin
    f = joinpath(p, "basic", "__site", "robots.txt")
    @test isfile(f)
    fc = prod(readlines(f, keep=true))

    @test occursin(raw"""
        Sitemap: https://tlienart.github.io/FranklinTemplates.jl/sitemap.xml""", fc)
end


text:
using Pkg
Pkg.activate(@__DIR__)
using Conda
using PyCall

# check what is contained in the julia env
Pkg.status()

# check what is contained in the conda env
Conda.list()

# check which python pycall is using
@show pyimport("sys").executable

max_stars_repo_name: jmaces/fw-rde

max_issues_repo_name: ZIB-IOL/fw-rde

max_forks_repo_name: ZIB-IOL/fw-rde


text:
@testset "Boundary utilities" begin
    rd = RefElemData(Tri(), N=3)
    md = MeshData(uniform_mesh(Tri(), 1)..., rd)
    on_bottom_boundary(xy, tol = 1e-13) = abs(xy[2]+1) < tol
    on_top_boundary(xy, tol = 1e-13) = abs(xy[2]-1) < tol

    boundary_dict = tag_boundary_faces(md, Dict(:bottom => on_bottom_boundary, :top => on_top_boundary))
    @test boundary_dict == Dict(:bottom=>[1],:top=>[4])

    boundary_dict = tag_boundary_faces(md, nothing)
    @test boundary_dict == Dict(:entire_boundary => [1,2,4,5])

    # test named tuple version
    boundary_dict = tag_boundary_faces(md,(; :bottom => on_bottom_boundary, :top => on_top_boundary))
    @test boundary_dict == (; :bottom=>[1],:top=>[4])
end


text:
using Makie

function draw_truss!(scene, X::Matrix{Float64}, T::Matrix{Int}, S::Matrix{Int}; area=undef, stress=undef, line_width::Float64=1.0, draw_supp::Bool=true, supp_scale::Float64=0.1, xaxis_label::String="x", plot_limits=undef)
    if area != undef
        @assert(size(T, 1) == length(area))
        a = reshape([area area]', 2*length(area))
        # a ./= maximum(a)
        a .*= line_width
    else
        a = ones(2*size(T,1))
        a .*= line_width
    end
    if stress != undef
        # stress color
        q_color = Array{RGBAf0, 1}(undef, length(stress))
        for i=1:length(stress)
            if stress[i] < 0
                q_color[i] = RGBAf0(0,0,1,0.8)
            else
                q_color[i] = RGBAf0(1,0,0,0.8)
            end
        end
        color = reshape([q_color q_color]', 2*length(stress))
    else
        color = :black
    end

    # seg_ids = reshape([T[:,1] T[:,2]], length(T[:,1])+length(T[:,2]))
    seg_ids = reshape(T', prod(size(T)))
    if plot_limits == undef
        max_xlim = maximum(X[:,1]) - minimum(X[:,1])
        max_ylim = maximum(X[:,2]) - minimum(X[:,2])
        max_lim = max(max_xlim, max_ylim)
        plot_limits = FRect(minimum(X[:,1]), minimum(X[:,2]),
                       max_lim, max_lim)
    end

    linesegments!(scene, X[seg_ids, 1], X[seg_ids, 2], linewidth = a,
                  color = color,
                  limits = plot_limits,
                  axis = (names = (axisnames = (xaxis_label, "y"),),
                          grid = (linewidth = (1, 1),),
                          )
                  )
    if draw_supp
        fix_ids = findall(x->x==1, S)
        for s in fix_ids
            if 1 == s[2]
                # x dir
                arrows!(scene, [X[s[1], 1]], [X[s[1], 2]], [supp_scale], [0], linecolor = :green, arrowcolor = :green, limits = plot_limits,
                # linewidth = supp_scale*10,
                arrowsize = 0.1,
                )
            end
            if 2 == s[2]
                # y dir
                arrows!(scene, [X[s[1], 1]], [X[s[1], 2]], [0], [supp_scale], linecolor = :green, arrowcolor = :green, limits = plot_limits,
                # linewidth = supp_scale*10,
                arrowsize = 0.1,
                )
            end
            # if 1 == S[i,4]
            #     scatter!(scene, [X[S[i,1], 1]], [X[S[i,1], 2]], color = :green, markersize=supp_scale, marker=:x, limits = plot_limits)
            # else
            #     scatter!(scene, [X[S[i,1], 1]], [X[S[i,1], 2]], color = :green, markersize=supp_scale, marker=:circle, limits = plot_limits)
            # end
        end
    end
    # axis = scene[Axis]
    # axis[:grid][:linewidth] = (1, 1)
end

function draw_load!(scene, X::Matrix{Float64}, F::Matrix{Float64}; load_scale::Float64=0.1, xaxis_label::String="x", plot_limits=undef)
    if plot_limits == undef
        max_xlim = maximum(X[:,1]) - minimum(X[:,1])
        max_ylim = maximum(X[:,2]) - minimum(X[:,2])
        max_lim = max(max_xlim, max_ylim)
        plot_limits = FRect(minimum(X[:,1]), minimum(X[:,2]),
                       max_lim, max_lim)
    end

    for i=1:size(F,1)
        if any(i -> i!=0, F[i,:])
            arrows!(scene, [X[i, 1]],
                           [X[i, 2]],
                           [F[i,1]*load_scale],
                           [F[i,2]*load_scale],
                           linecolor = :orange, arrowcolor = :orange,
                           # linewidth = load_scale*10,
                           arrowsize = 0.1,
                           limits = plot_limits,
                           axis = (names = (axisnames = (xaxis_label, "y"),),)
                           )
         end
    end
    # axis = scene[Axis]
    # axis[:grid][:linewidth] = (1, 1)
end


text:
function satisfy(interpreter::ConcreteInterpreter,
                 domain::Domain, state::GenericState,
                 terms::AbstractVector{<:Term})
     # Quick check that avoids SLD resolution unless necessary
     sat = all(check(interpreter, domain, state, t) for t in terms)
     if sat !== missing return sat end
     # Call SLD resolution only if there are free variables or axioms
     return !isempty(satisfiers(interpreter, domain, state, terms))
end

function satisfy(interpreter::ConcreteInterpreter,
                 domain::Domain, state::GenericState, term::Term)
    # Quick check that avoids SLD resolution unless necessary
    sat = check(interpreter, domain, state, term)
    if sat !== missing return sat end
    # Call SLD resolution only if there are free variables or axioms
    return !isempty(satisfiers(interpreter, domain, state, term))
end

function satisfiers(interpreter::ConcreteInterpreter,
                    domain::Domain, state::GenericState,
                    terms::AbstractVector{<:Term})
    # Initialize Julog knowledge base
    clauses = Clause[get_clauses(domain); collect(state.types); collect(state.facts)]
    # Pass in fluents and function definitions as a dictionary of functions
    funcs = merge(GLOBAL_FUNCTIONS, state.values, get_funcdefs(domain))
    return resolve(collect(terms), clauses; funcs=funcs, mode=:all)[2]
end

function satisfiers(interpreter::ConcreteInterpreter,
                    domain::Domain, state::GenericState, term::Term)
    return satisfiers(interpreter, domain, state, [term])
end

function check(interpreter::ConcreteInterpreter,
               domain::Domain, state::GenericState, term::Compound)
    sat = if term.name == :and
        all(check(interpreter, domain, state, a) for a in term.args)
    elseif term.name == :or
        any(check(interpreter, domain, state, a) for a in term.args)
    elseif term.name == :imply
        !check(interpreter, domain, state, term.args[1]) |
        check(interpreter, domain, state, term.args[2])
    elseif term.name == :not
        !check(interpreter, domain, state, term.args[1])
    elseif is_quantifier(term)
        missing
    elseif !is_ground(term)
        missing
    elseif is_derived(term, domain)
        missing
    elseif is_type(term, domain)
        if has_subtypes(term, domain)
            missing
        elseif term in state.types
            true
        elseif !isempty(get_constants(domain))
            domain.constypes[term.args[1]] == term.name
        else
            false
        end
    elseif is_global_pred(term)
        evaluate(interpreter, domain, state, term)::Bool
    elseif is_func(term, domain)
        evaluate(interpreter, domain, state, term)::Bool
    else
        term = partialeval(domain, state, term)
        term in state.facts
    end
    return sat
end

function check(interpreter::ConcreteInterpreter,
               domain::Domain, state::GenericState, term::Const)
    if term in state.facts || term in state.types || term.name == true
        return true
    elseif is_func(term, domain) || is_derived(term, domain)
        return missing
    else
        return false
    end
end

function check(interpreter::ConcreteInterpreter,
               domain::Domain, state::GenericState, term::Var)
    return missing
end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

n = 200 #3000 #1000 #20000
t = 0.35 #0.15 #0.4 #0.15

function randlines(n=1000,t=0.2)
	V = zeros(Float64,2,2*n)
	EV = [zeros(Int64,2) for k=1:n]
	for k=1:n
		v1 = rand(Float64,2)
		v2 = rand(Float64,2)
		vm = (v1+v2)/2
		transl = rand(Float64,2)
		V[:,k] = (v1-vm)*t + transl
		V[:,n+k] = (v2-vm)*t + transl
		EV[k] = [k,n+k]
	end
	V = GL.normalize2(V)
	model = (V,EV)
	return model
end

model = randlines(n,t)
Sigma = Lar.spaceindex(model)
GL.VIEW([ GL.GLLines(model...)] )

W,EW = Lar.fraglines(1.5,1.5,1.5)(model)
GL.VIEW([ GL.GLLines(W,EW,GL.COLORS[1])] );

W,EW = Lar.fragmentlines(model)
V,EVs = Lar.biconnectedComponent((W,EW::Lar.Cells)) # 2-connected components (H & T)

graphs = [ GL.GLLines(V,EVs[k],GL.COLORS[(k-1)%12+1]) for k=1:length(EVs) ]
GL.VIEW(graphs);


text:
using Luxor, Test, Colors, Random

Random.seed!(42)

function noisetest(fname)
    Drawing(800, 800, fname)
    background("chartreuse4")
    origin()

    initnoise(rand(1:12))

    @test 0.0 < noise(.5) <= 1.0
    @test 0.0 < noise(0.1, 0.5) <= 1.0
    @test 0.0 < noise(0.1, 0.1, -0.1) <= 1.0
    @test 0.0 < noise(0.5, 2.0, -2.0, 0.1) <= 1.0

    freq = 0.02
    tiles = Tiler(800, 800, 150, 150)
    for k in 1:2:10
        for (pos, n) in tiles
            f, d = .01, k
            ns = noise(pos.x * freq, pos.y * freq, detail=3, persistence=.3)
            ns1 = noise(ns, detail = 3, persistence=2)
            ns2 = noise(pos.x * freq, pos.y * freq, ns, detail = 2)
            ns3 = noise(pos.x * freq, pos.y * freq, ns2)
            setopacity(ns3)
            sethue(LCHab(80ns, 100 * ns1, 360 * ns2))
            box(pos, tiles.tilewidth, tiles.tileheight, :fill)
        end
    end

    # test a custom rng
    rng = MersenneTwister(1234)
    initnoise(rng)

    @test 0.0 < noise(.5) <= 1.0
    @test 0.0 < noise(0.1, 0.5) <= 1.0
    @test 0.0 < noise(0.1, 0.1, -0.1) <= 1.0
    @test 0.0 < noise(0.5, 2.0, -2.0, 0.1) <= 1.0

    @test finish() == true
end

fname = "noise-test.png"
noisetest(fname)
println("...finished noise test: output in $(fname)")


text:
using VolterraGP
using Test
using BenchmarkTools
using LinearAlgebra
using Flux

P = 1

X = [collect(-1:0.1:1), collect(-1:0.1:1)]
Y = [sin.(X[1]), cos.(X[2])]

Xd = [collect(-1:0.1:1), collect(0:0.2:1.5)]
Yd = [sin.(Xd[1]), cos.(Xd[2])]

t = [collect(-1.5:0.3:1.5), collect(-1.5:0.3:1.5)]

data = Data(X, Y)
datad = Data(Xd, Yd)

dpars1 = DiffableParameters([0.1, 0.1], ones((2, sum(1:1), P)), [0.1])
gp1 = GaussianProcess(threeEQs, 2, 1, P, data, dpars1)

dpars2 = DiffableParameters([0.1, 0.1], ones((2, sum(1:2), P)), [0.1])
gp2 = GaussianProcess(threeEQs, 2, 2, P, data, dpars2)

dparsrand = DiffableParameters([0.3, 0.4], 0.1 .+ rand(Float64, (2, sum(1:2), P)), [0.1])
gprand = GaussianProcess(threeEQs, 2, 2, P, data, dparsrand)

gpd = GaussianProcess(threeEQs, 2, 2, P, datad, dpars2)

dparssca = DiffableParameters([0.1, 0.1], 0.5ones((2, sum(1:2), 2)), [0.1])
gpscaled = GaussianProcess(scaledEQs, 2, 2, P, data, dparssca)

test_psd = K -> minimum(eigvals(K)) > -sqrt(eps())
test_herm = K -> maximum(K' - K) < sqrt(eps())

# gp3

@testset "VolterraGP.jl" begin
    @test negloglikelihood(gp1) ≈ 475.6091937696483
    @test negloglikelihood(gp2) ≈ 476.28491875918456
    @test negloglikelihood(gpd) ≈ 48.92347237044826
    @test negloglikelihood(gpscaled) ≈ 477.9933646808806
end

@testset "volterra.jl" begin
    # when C=1 running thru volterra should be same as base
    @test VolterraGP.kernel(1., 2., 1, 2, gp1) ≈ threeEQs(1., 2., gp1.dpars.G[1, 1, :], gp1.dpars.G[2, 1, :], gp1.dpars.u)
    @test VolterraGP.kernel(1., 2., 1, 1, gp1) ≈ threeEQs(1., 2., gp1.dpars.G[1, 1, :], gp1.dpars.G[1, 1, :], gp1.dpars.u)

    # test for symmerty in kernel arguments
    @test VolterraGP.kernel(1., 2., 1, 2, gprand) ≈ VolterraGP.kernel(1., 2., 2, 1, gprand)
    @test VolterraGP.kernel(1., 2., 1, 1, gprand) ≈ VolterraGP.kernel(2., 1., 1, 1, gprand)

    # examples from paper
    @test VolterraGP.kan_rv_prod(ones(Float64, (4, 4))) ≈ 3.
    @test VolterraGP.kan_rv_prod(ones(Float64, (2, 2))) ≈ 1.

    phi_E = VolterraGP.get_phi_E(0.2, 2, 1, gprand)
    @test test_psd(phi_E)
    @test test_herm(phi_E)

    phi_K = VolterraGP.get_phi_cov(0.1, 0.2, 2, 2, 2, 2, gprand)
    @test test_psd(phi_K)
    @test test_herm(phi_K)

end

@testset "gp.jl" begin
    for gp in (gp1, gp2, gpd, gpscaled)

        K = VolterraGP.fill_K(t, t, gp)
        μp, Kp = posterior(t, gp)

        @test test_psd(K)
        @test test_herm(K)
        @test test_psd(Kp)
        @test test_herm(Kp)
    end

end

@testset "gradients" begin
    # test gradient works
    g1 = gradient(negloglikelihood, gp1)
    @test g1 == g1
    g2 = gradient(negloglikelihood, gp2)
    @test g2 == g2
    gr = gradient(negloglikelihood, gprand)
    @test gr == gr
    gd = gradient(negloglikelihood, gpd)
    @test gd == gd
    gr = gradient(negloglikelihood, gpscaled)
    @test gr == gr

    # test forward vs reverse gradients
    g_fwd = VolterraGP.fwd_grad_like(gpscaled)
    g_bwd = Flux.gradient(Flux.params(gpscaled.dpars.σ, gpscaled.dpars.G, gpscaled.dpars.u)) do
        negloglikelihood(gpscaled)
    end

    @test all(g_bwd[gpscaled.dpars.σ] .≈ g_fwd[1])
    @test all(g_bwd[gpscaled.dpars.G] .≈ g_fwd[2])
    @test all(g_bwd[gpscaled.dpars.u] .≈ g_fwd[3])

end


text:
"""
SNOPT wrapper
"""

"""
    minimize_solver_specific(solver::SNOPT, cache, x0, lx, ux, lg, ug, rows, cols,
        outputfile::Bool, lencw::Int=500)

Minimize function interfacing to snopta() routine

# Arguments
    - `cache::AbstractSolver`: cache created using _create_cache()
"""
function minimize_snopt(options::Dict, cache, x0, lx, ux, lg, ug, rows, cols,
            outputfile::Bool, lencw::Int=500, iSumm::Int=6)

    # define objective function for SNOPT
    function fsnopt!(g, df, dg, x, deriv)
        f = evaluate!(g, df, dg, x, cache)
        fail = false  # check?
        # potentially provide df if user-defined
        return f, fail
    end

    # call snopta
    start = ColdStart(x0, length(lg)+1)
    xstar, fstar, info, out = snopta(fsnopt!, start, lx, ux, lg, ug, rows, cols;
                                    options=options, lencw=lencw)
    return xstar, fstar, info
end

max_stars_repo_name: Yuricst/joptimise

max_issues_repo_name: Yuricst/joptimise

max_forks_repo_name: Yuricst/joptimise


text:
"
This script contains the functions to evaluate the expectation values ⟨Â⟩ for the three oscillators

The parameters to be given are: ω, λ, g, γ, α, and a list of sampled times t

"

export exv_a, exv_a₁, exv_a₂

## Definitions of the parameters and function parameters

ωₘ(ω, λ) = ω - λ
ωₚ(ω, λ) = ω + λ

Ω(ω, λ, g) = (1/2) * (ω + ωₚ(ω, λ) - √(8 * g^2 + (ωₚ(ω, λ) - ω)^2))
Ω₂(ω, λ, g) = (1/2) * (ω + ωₚ(ω, λ) + √(8 * g^2 + (ωₚ(ω, λ) - ω)^2))

c2ϕ(ω, λ, g) = (ωₚ(ω, λ) - ω)/√(8 * g^2 + (ωₚ(ω, λ) - ω)^2)
s2ϕ(ω, λ, g) = (2 * √2 * g)/√(8 * g^2 + (ωₚ(ω, λ) - ω)^2)

cc2ϕ(ω, λ, g) = (1/2) * (1 + c2ϕ(ω, λ, g))
ss2ϕ(ω, λ, g) = (1/2) * (1 - c2ϕ(ω, λ, g))

##
## Definition of the exponentials

expΩΩ₂(ω, λ, g, γ, t) = exp(big(γ * t * exp((1im/γ) * (Ω(ω, λ, g) - Ω₂(ω, λ, g)))))

expωₘΩ₂(ω, λ, g, γ, t) = exp(big(γ * t * exp((1im/γ) * (ωₘ(ω, λ) - Ω₂(ω, λ, g)))))

expωₘΩ(ω,λ, g, γ, t) = exp(big(γ * t * exp((1im/γ) * (ωₘ(ω, λ) - Ω(ω, λ, g)))))

##
## Definition of the functions who gives the expectation values

"`exv_a(ω, λ, g, γ, α, t)` gives the expectation value over time of the oscillator ⟨â⟩(t). The function will be waiting for the parameters {ω, λ, g, γ, α} and a list of sampled times."
exv_a(ω, λ, g, γ, α, t) = abs2(α) * cc2ϕ(ω, λ, g) * ss2ϕ(ω, λ, g) * (1 - (1/2) * exp(-γ * t) * (expΩΩ₂(ω, λ, g, γ, t) + conj(expΩΩ₂(ω, λ, g, γ, t))))

"`exv_a₁(ω, λ, g, γ, α, t)` gives the expectation value over time of the oscillator ⟨â₁⟩(t). The function will be waiting for the parameters {ω, λ, g, γ, α} and a list of sampled times."
exv_a₁(ω, λ, g, γ, α, t) = (1/4) * abs2(α) * (1 + cc2ϕ(ω, λ, g)^2 + ss2ϕ(ω, λ, g)^2 + exp(-γ * t) * ((expΩΩ₂(ω, λ, g, γ, t) + conj(expΩΩ₂(ω, λ, g, γ, t))) * cc2ϕ(ω, λ, g) * ss2ϕ(ω, λ, g) + (expωₘΩ₂(ω, λ, g, γ, t) + conj(expωₘΩ₂(ω, λ, g, γ, t))) * cc2ϕ(ω, λ, g) + (expωₘΩ(ω, λ, g, γ, t) + conj(expωₘΩ(ω, λ, g, γ, t))) * ss2ϕ(ω, λ, g)))

"`exv_a₂(ω, λ, g, γ, α, t)` gives the expectation value over time of the oscillator ⟨â₂⟩(t). The function will be waiting for the parameters {ω, λ, g, γ, α} and a list of sampled times."
exv_a₂(ω, λ, g, γ, α, t) = (1/4) * abs2(α) * (1 + cc2ϕ(ω, λ, g)^2 + ss2ϕ(ω, λ, g)^2 + exp(-γ * t) * ((expΩΩ₂(ω, λ, g, γ, t) + conj(expΩΩ₂(ω, λ, g, γ, t))) * cc2ϕ(ω, λ, g) * ss2ϕ(ω, λ, g) - (expωₘΩ₂(ω, λ, g, γ, t) + conj(expωₘΩ₂(ω, λ, g, γ, t))) * cc2ϕ(ω, λ, g) - (expωₘΩ(ω, λ, g, γ, t) + conj(expωₘΩ(ω, λ, g, γ, t))) * ss2ϕ(ω, λ, g)))

##

"WARNING: There's overflows issues for when the exponentials goes to very large numbers, their evaluations returns numbers below the machine precision."

max_stars_repo_name: rurz/IntrinsicDecoherence

max_issues_repo_name: rurz/IntrinsicDecoherence

max_forks_repo_name: rurz/IntrinsicDecoherence


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

struct AppsPermissionsResourcesApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_appsPermissionsResourcesList(_api::AppsPermissionsResourcesApi; cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/apps.permissions.resources.list", ["slackAuth"])
    Swagger.set_param(_ctx.query, "cursor", cursor)  # type String
    Swagger.set_param(_ctx.query, "limit", limit)  # type Int32
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Returns list of resource grants this app has on a team.
Param: cursor::String
Param: limit::Int32
Return: Dict{String, Any}
"""
function appsPermissionsResourcesList(_api::AppsPermissionsResourcesApi; cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_appsPermissionsResourcesList(_api; cursor=cursor, limit=limit, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function appsPermissionsResourcesList(_api::AppsPermissionsResourcesApi, response_stream::Channel; cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_appsPermissionsResourcesList(_api; cursor=cursor, limit=limit, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export appsPermissionsResourcesList


text:
type MHRemoteSegments{T<:AbstractFloat} <: AbstractRemoteSegments
    @compat remote::Vector{Future}
    numproposalspersegment::Int
    numsegments::Int

    acceptances::Vector
    prop2collected::Dict{Int,Tuple{Int,Int}}
    collectedsamples::Vector{AbstractSample}
    @compat MHRemoteSegments(r::Vector{Future},npps::Int,ns::Int,::Type{T}) = new(r,npps,ns,zeros(T,npps*ns),Dict{Int,Tuple{Int,Int}}(),Vector{AbstractSample}(ns))
end

@inline _numjobsegments(policy_::MHRuntimePolicy,nproposals::Int) = min(nproposals,_numjobsegments(traittype(policy_.jobsegments)))
@inline _processnumbers(policy_::MHRuntimePolicy,njobsegments::Int) = collect(take(cycle(_processnumbers(traittype(policy_.jobsegments))),njobsegments))

function _remotesegments(policy_::MHRuntimePolicy,model_::AbstractModel,sampler_::AbstractSampler,nproposals::Int)
    njobsegments = _numjobsegments(policy_,nproposals)
    nproposalspersegment = _numproposalspersegment(nproposals,njobsegments)
    @compat r = Future[remotecall(segment,i,policy_,model_,sampler_,nproposalspersegment) for i in _processnumbers(policy_,njobsegments)]
    MHRemoteSegments{policy_.calculationtype}(r,nproposalspersegment,njobsegments,policy_.calculationtype)
end

_prop2seg(segments_::MHRemoteSegments,j::Int) = ind2sub((segments_.numproposalspersegment,segments_.numsegments),j)

function _insegmentindex(segments_::MHRemoteSegments,sampleindex_::AbstractVector)
    boolindex_ = [fill(false,segments_.numproposalspersegment) for k=1:segments_.numsegments]
    ntp = numtotalproposals(segments_)
    for i=1:length(sampleindex_)
        @inbounds i1 = sampleindex_[i]-1
        i1==ntp?continue:nothing #if this the final
        si = div(i1,segments_.numproposalspersegment) + 1
        pi = mod(i1,segments_.numproposalspersegment) + 1
        @inbounds boolindex_[si][pi] = true
    end
    map(find,boolindex_)
end

function _insegmentindex2collected!(segments_::MHRemoteSegments,insegmentindex::Array{Array{Int,1},1})
    empty!(segments_.prop2collected)
    d = (segments_.numproposalspersegment,segments_.numsegments)
    for i=1:segments_.numsegments
        for j=1:length(insegmentindex[i])
            segments_.prop2collected[sub2ind(d,insegmentindex[i][j],i)] = (j,i)
        end
    end
    segments_.prop2collected
end

#start the jobs on each process
function iterate!(segments_::MHRemoteSegments,indicatorstate::AbstractSamplerState)
    @compat a = Array{Future}(segments_.numsegments)
    @sync begin
        for j=1:segments_.numsegments
            a[j] = remotecall(iterate!,segments_.remote[j].where,segments_.remote[j],indicatorstate)
        end
    end #@sync means wait for all processes to finish
    iCounter = 1
    for j=1:segments_.numsegments
        copy!(segments_.acceptances,iCounter,fetch(a[j]),1,segments_.numproposalspersegment)
        iCounter += segments_.numproposalspersegment
    end
    segments_.acceptances
end

#call prepare in order to copy over the new indicator state
function prepare!(segments_::MHRemoteSegments,indicatorstate::AbstractSamplerState,j::Int)
    p,s = _prop2seg(segments_,j)
    r = segments_.remote[s]
    remotecall_fetch(prepare!,r.where,r,indicatorstate,p)
end

function retrievesamples!(segments_::MHRemoteSegments,sampleindex_::AbstractVector)
    insegmentindex = _insegmentindex(segments_,sampleindex_)
    segments_.prop2collected = _insegmentindex2collected!(segments_,insegmentindex)
    @sync begin
        map!((r,i)->(~isempty(i)?remotecall_fetch(getsamples,r.where,r,i):samples(:base,0,0,Float64,Float64)),segments_.collectedsamples,segments_.remote,insegmentindex)
    end
    segments_.collectedsamples
end

function getsamples(segments_::MHRemoteSegments,j::Int)
    if !haskey(segments_.prop2collected,j)
        p,s = _prop2seg(segments_,j)
        r = segments_.remote[s]
        return remotecall_fetch(getsamples,r.where,r,p)
    else
        p,s = segments_.prop2collected[j]
        return copy(segments_.collectedsamples[s],p)
    end
end

function getsamplerstatevars(segments_::MHRemoteSegments)
    map((r)->remotecall_fetch(getsamplerstatevars,r.where,r),segments_.remote)
end

function store!(segments_::MHRemoteSegments,chain_::AbstractChain,j::Int)
    p,s = segments_.prop2collected[j]
    store!(chain_,segments_.collectedsamples[s],p)
end

function tune!(segments_::MHRemoteSegments,tvals...)
    @sync begin
        for i=1:segments_.numsegments
            r = segments_.remote[i]
            remotecall_wait(tune!,r.where,r,tvals...)
        end
    end
end

function show(io::IO,r::MHRemoteSegments)
    println(io,"RemoteSegments with $(r.numsegments) segment",r.numsegments==1?"":"s"," and $(r.numproposalspersegment) proposal",r.numproposalspersegment==1?"":"s"," per segment.")
    println(io,"Additional fields: :remote, :acceptances, :collectedsamples, :prop2collected")
    nothing
end