text:
## Abstract Metropolized Sampliing
## ===============================

"""Abstract Metropolized Sampliing samples events in preimage uniformly in
convergence of the Markov Chain.  This algorithm is useful for high dimensions"""
immutable AMS <: MCMCAlgorithm end

"Proposes a box by querying satisfiability solver multiple times"
function proposebox_tl{D <: Domain}(
    X::RandVar, box::D;
    split::Function = weighted_partial_split,
    maxdepth::Int = 1000,
    precision::Float64 = DEFAULT_PREC,
    args...)

  niters = 0
  depth = 0
  logq = 0.0 # == log(1.0)
  prevolfrac = 1.0

  A::D = box
  lens(:refine,time_ns())
  @show image::AbstractBool = call(X, A)
  while (niters <= 1000) && (depth <= maxdepth)
    @show depth
    @show measure(A)
    A, image
    if issmall(A, precision)
      lens(:depth, depth)
      lens(:refine,time_ns())
      println("Got Sample")
      return A, logq, 1.0  # Assume boxes are full
    # else if  isequal(image,t)
    #   lens(:proposing, depth=depth, niters=niters)
    #   return A, logq, 1.0  # Assume boxes are full
    elseif isequal(image, tf)
      @compat children::Vector{Tuple{Domain,Float64}} = split(A, depth)
      statuses = AbstractBool[]

      # Due to bug in DReal we're getting both a query and its negation unsat
      for child in children
        try
          child_status = call(X, child[1])
          push!(statuses,child_status)
        catch e
          rethrow(e)
          println("caught exception $e")
          println(A)
          println(X)
          [println(child) for child in children]
          push!(statuses,f)
        end
      end
      weights = pnormalize([isequal(statuses[i],f) ? 0.0 : children[i][2] for i = 1:length(children)])

      # Sometimes all the children are false, even thoughrand pnarent is true, bug?
      if all([isequal(status,f) for status in statuses])
        @show "Found Bug $depth"
        @show statuses
        @show children
        @show A
        @show typeof(X)
        lens(:depth, depth)
        lens(:refine,time_ns())
        return A, logq, 1.0
      end

      # Choose a random child
      # @show statuses
      # @show weights
      # @show A
      # @show children
      rand_index = rand(Categorical(weights))
      A = children[rand_index][1]
      status = statuses[rand_index]

      # Shouldn't go to empty children
      assert(!isequal(status,f))
      logq += log(weights[rand_index])

      depth += 1; niters += 1
      lens(:children_split, weights, children, depth, niters)
    elseif isequal(image, f) # Condition is unsatisfiable
      error("Cannot condition on unsatisfiable events")
    end
  end
  error("Unexpected Branch")
end

"Propose `n` boxes"
function propose_boxes{D <: Domain}(
    Y::SymbolicRandVar{Bool},
    init_box::D,
    nsamples::Integer;
    RandVarType::Type = default_randvar(),
    precision::Float64 = default_precision(),
    args...)

  Y_conv = convert(RandVarType{Bool}, Y)
  set_precision!(Y_conv, precision)

  # DReal.opensmt_print_expr(Y_conv.ex.e)
  @compat samples = Array(Tuple{D,Float64,Float64}, nsamples)
  for i = 1:nsamples
    before = time_ns()
    samples[i] = preimage_proposal(Y_conv, init_box; args...)
    after = time_ns()
    lens(:sat_check, after-before)
  end
  cleanup(Y_conv)
  samples
end

## Parallel
## ========
"Generate a stack of `n` boxes in parallel"
function propose_boxes_parallel{D<:Domain}(
    Y::SymbolicRandVar,
    init_box::D,
    nsamples::Integer;
    ncores = nprocs() - 1,
    args...)

  println("Using $ncores cores to generate $nsamples")
  samplespercore = div(nsamples, ncores) + 1
  lst = [i for i = 1:ncores]
  g = _ -> propose_boxes(Y, init_box, samplespercore; args...)
  a = pmaplm(g, lst; ncores = min(nprocs(),ncores))
  samples = Tuple{D,Float64,Float64}[]
  for sampleset in a
    append!(samples, sampleset)
  end
  return samples
end

"From set of boxes run a markov chain in batch mode"
function run_chain{D <: Domain}(boxes::Vector{@compat Tuple{D,Float64,Float64}})
  box, logq, prevolfrac = boxes[1]
  logp = logmeasure(box) + log(prevolfrac)
  chain = Array(D, length(boxes))
  chain[1] = box

  naccepted = 0
  for i = 2:length(boxes)
    nextbox, nextlogq, prevolfrac = boxes[i]
    nextlogp = logmeasure(nextbox) + log(prevolfrac)

    loga = nextlogp + logq - logp - nextlogq
    a = exp(loga)

    # MH accept/reject step
    if a >= 1 || rand() < a
      naccepted += 1
      box = nextbox
      logp = nextlogp
      logq = nextlogq
    end
    chain[i] = box
  end
  lens(:mh_ratio, naccepted/length(boxes))

  chain
end

"Uniform sample of subset of preimage of `Y` under `f` unioned with `X`."
function pre_mc(
    Y::SymbolicRandVar{Bool},
    nsamples::Integer,
    ::Type{AMS};
    RandVarType::Type = default_randvar(),
    parallel::Bool = true,
    args...)

  init_box = unit_box(LazyBox{Float64}, dims(Y))
  if parallel && nprocs() > 1
    boxes = propose_boxes_parallel(Y, init_box, nsamples; RandVarType = RandVarType, args...)
  else
    boxes = propose_boxes(Y, init_box, nsamples; RandVarType = RandVarType, args...)
  end

  run_chain(boxes)
end


text:
module MPCCmod

using NLPModels

import Relaxation
import ParamSetmod
import AlgoSetmod

"""
Definit le type MPCC :
min_x f(x)
l <= x <= u
lb <= c(x) <= ub
0 <= G(x) _|_ H(x) >= 0

liste des constructeurs :
MPCC(f::Function,x0::Vector,G::Function,H::Function,nb_comp::Int64,lvar::Vector,uvar::Vector,c::Function,y0::Vector,lcon::Vector,ucon::Vector)
MPCC(f::Function,x0::Vector,G::Function,H::Function,nb_comp::Int64,lvar::Vector,uvar::Vector,c::Function,y0::Vector,lcon::Vector,ucon::Vector,prec::Float64)
MPCC(mp::NLPModels.AbstractNLPModel,G::Function,H::Function,nb_comp::Int64)

liste des accesseurs :
addInitialPoint(mod::MPCC,x0::Vector)

liste des fonctions :
viol_contrainte_norm(mod::MPCCmod.MPCC,x::Vector,yg::Vector,yh::Vector)
viol_contrainte_norm(mod::MPCCmod.MPCC,x::Vector)
viol_contrainte(mod::MPCCmod.MPCC,x::Vector,yg::Vector,yh::Vector)
viol_contrainte(mod::MPCCmod.MPCC,x::Vector)
viol_comp(mod::MPCCmod.MPCC,x::Vector)
viol_cons(mod::MPCCmod.MPCC,x::Vector)

MPCCtoRelaxNLP(mod::MPCC, r::Float64, s::Float64, t::Float64, relax::AbstractString)
"""

# TO DO List
#Major :
# - MPCCtoRelaxNLP : bug à corriger
# - appel de la hessienne d'un SimpleNLPModel ? NLPModels.hess(mod.H,x)

type MPCC
 mp::NLPModels.AbstractNLPModel
 G::NLPModels.AbstractNLPModel
 H::NLPModels.AbstractNLPModel
 xj::Vector #itéré courant
 nb_comp::Int64

 algoset::AlgoSetmod.AlgoSet
 paramset::ParamSetmod.ParamSet
end

#Constructeurs supplémentaires :
function MPCC(f::Function,x0::Vector,G::NLPModels.AbstractNLPModel,H::NLPModels.AbstractNLPModel,nb_comp::Int64,lvar::Vector,uvar::Vector,c::Function,lcon::Vector,ucon::Vector)

 mp=ADNLPModel(f, x0, lvar=lvar, uvar=uvar, c=c, lcon=lcon, ucon=ucon)
 nbc=length(mp.meta.lvar)+length(mp.meta.uvar)+length(mp.meta.lcon)+length(mp.meta.ucon)+2*nb_comp

 return MPCC(mp,G,H,x0,nb_comp,AlgoSetmod.AlgoSet(),ParamSetmod.ParamSet(nbc))
end

function MPCC(mp::NLPModels.AbstractNLPModel,G::NLPModels.AbstractNLPModel,H::NLPModels.AbstractNLPModel,nb_comp)
 nbc=length(mp.meta.lvar)+length(mp.meta.uvar)+length(mp.meta.lcon)+length(mp.meta.ucon)+2*nb_comp

 return MPCC(mp,G,H,mp.meta.x0,nb_comp,AlgoSetmod.AlgoSet(),ParamSetmod.ParamSet(nbc))
end

function MPCC(mp::NLPModels.AbstractNLPModel)

 G=SimpleNLPModel(x->0, [0.0])
 H=SimpleNLPModel(x->0, [0.0])

 nb_comp=0
 nbc=length(mp.meta.lvar)+length(mp.meta.uvar)+length(mp.meta.lcon)+length(mp.meta.ucon)+2*nb_comp

 return MPCC(mp,G,H,mp.meta.x0,nb_comp,
		AlgoSetmod.AlgoSet(),ParamSetmod.ParamSet(nbc))
end

function MPCC(mp::NLPModels.AbstractNLPModel,algo::AlgoSetmod.AlgoSet)

 G=SimpleNLPModel(x->0, [0.0])
 H=SimpleNLPModel(x->0, [0.0])

 nb_comp=0
 nbc=length(mp.meta.lvar)+length(mp.meta.uvar)+length(mp.meta.lcon)+length(mp.meta.ucon)+2*nb_comp

 return MPCC(mp,G,H,mp.meta.x0,nb_comp,
		AlgoSetmod.AlgoSet(),ParamSetmod.ParamSet(nbc))
end

function MPCC(mp::NLPModels.AbstractNLPModel,G::NLPModels.AbstractNLPModel,H::NLPModels.AbstractNLPModel;nb_comp::Float64=NaN)

 nb_comp=isnan(nb_comp)?length(NLPModels.cons(G,mp.meta.x0)):nb_comp
 nbc=length(mp.meta.lvar)+length(mp.meta.uvar)+length(mp.meta.lcon)+length(mp.meta.ucon)+2*nb_comp

 return MPCC(mp,G,H,mp.meta.x0,nb_comp,
		AlgoSetmod.AlgoSet(),ParamSetmod.ParamSet(nbc))
end

"""
Accesseur : modifie le point initial
"""

function addInitialPoint(mod::MPCC,x0::Vector)
 mod.xj=x0
 return mod
end

"""
Donne la norme 2 de la violation des contraintes avec slack

note : devrait appeler viol_contrainte
"""
function viol_contrainte_norm(mod::MPCCmod.MPCC,x::Vector,yg::Vector,yh::Vector;tnorm::Real=2)
 return norm(viol_contrainte(mod,x,yg,yh),tnorm)
end

function viol_contrainte_norm(mod::MPCCmod.MPCC,x::Vector;tnorm::Real=2) #x de taille n+2nb_comp

 n=length(mod.mp.meta.x0)
 if length(x)==n
  resul=max(viol_comp(mod,x),viol_cons(mod,x))
 else
  resul=viol_contrainte_norm(mod,x[1:n],x[n+1:n+mod.nb_comp],x[n+mod.nb_comp+1:n+2*mod.nb_comp],tnorm=tnorm)
 end
 return resul
end

"""
Donne le vecteur de violation des contraintes dans l'ordre : G(x)-yg ; H(x)-yh ; lvar<=x ; x<=uvar ; lvar<=c(x) ; c(x)<=uvar
"""
function viol_contrainte(mod::MPCCmod.MPCC,x::Vector,yg::Vector,yh::Vector)

 c=NLPModels.cons(mod.mp,x)
 if mod.nb_comp>0
  G=NLPModels.cons(mod.G,x)
  H=NLPModels.cons(mod.H,x)
  return [G-yg;H-yh;max.(mod.mp.meta.lvar-x,0);max.(x-mod.mp.meta.uvar,0);max.(mod.mp.meta.lcon-c,0);max.(c-mod.mp.meta.ucon,0)]
 else
  return [yg;yh;max.(mod.mp.meta.lvar-x,0);max.(x-mod.mp.meta.uvar,0);max.(mod.mp.meta.lcon-c,0);max.(c-mod.mp.meta.ucon,0)]
 end
end

function viol_contrainte(mod::MPCCmod.MPCC,x::Vector) #x de taille n+2nb_comp
 n=length(mod.mp.meta.x0)

 return viol_contrainte(mod,x[1:n],x[n+1:n+mod.nb_comp],x[n+mod.nb_comp+1:n+2*mod.nb_comp])
end

"""
Donne la norme de la violation de la complémentarité min(G,H)
"""
function viol_comp(mod::MPCCmod.MPCC,x::Vector;tnorm::Real=2)

 n=length(mod.mp.meta.x0)
 x=length(x)==n?x:x[1:n]
 G=NLPModels.cons(mod.G,x)
 H=NLPModels.cons(mod.H,x)

 return mod.nb_comp>0?norm(G.*H./(G+H+1),tnorm):0
end

"""
Donne la norme de la violation des contraintes \"classiques\"
"""
function viol_cons(mod::MPCCmod.MPCC,x::Vector;tnorm::Real=2)

 n=length(mod.mp.meta.x0)
 x=length(x)==n?x:x[1:n]
 feas=norm([max.(mod.mp.meta.lvar-x,0);max.(x-mod.mp.meta.uvar,0)],tnorm)

 if mod.mp.meta.ncon !=0
  c=NLPModels.cons(mod.mp,x)
  #feas=max(max(maximum(mod.mp.meta.lcon-c),maximum(c-mod.mp.meta.ucon)),feas)
  feas=max(norm([max.(mod.mp.meta.lcon-c,0);max.(c-mod.mp.meta.ucon,0)],tnorm),feas)
 end

 return feas
end

"""
MPCCtoRelaxNLP(mod::MPCC, t::Float64)
mod : MPCC
return : le MPCC en version NL pour un t donné
"""
function MPCCtoRelaxNLP(mod::MPCC, r::Float64, s::Float64, t::Float64, relax::AbstractString)
 G(x)=NLPModels.cons(mod.G,x)
 H(x)=NLPModels.cons(mod.H,x)

 #concatène les contraintes de complémentarité + positivité :
 nl_constraint(x)=emptyfunc
 if relax=="SS"
  nl_constraint(x)=[mod.mp.c(x);G(x).*H(x)-t;G(x);H(x)]
 elseif relax=="KDB" #(G(x)-s)(H(x)-s)<=0, G(x)>=-r, H(x)>=-r
  nl_constraint(x)=[mod.mp.c(x);(G(x)-s).*(H(x)-s);G(x)+r;H(x)+r]
 elseif relax=="KS" #si G(x)-s+H(x)-s>=0 ? (G(x)-s)(H(x)-s)<=0 : -1/2*((G(x)-s)^2+(H(x)-s)^2), G(x)>=0, H(x)>=0
  KS(x)= G(x)-s+H(x)-s>=0 ? (G(x)-s).*(H(x)-s) : -0.5*((G(x)-s).^2+(H(x)-s).^2)
  nl_constraint(x)=[mod.mp.c(x);KS(x);G(x);H(x)]
 elseif relax=="Butterfly"
# On devrait appeler Relaxation et pas Thetamod
#  FG(x)=mod.G(x)-s-t*Thetamod.theta(mod.H(x)-s,r) Bug à corriger
#  FH(x)=mod.H(x)-s-t*Thetamod.theta(mod.G(x)-s,r) Bug à corriger
  FG(x)=G(x)-s-t*(H(x)-s)
  FH(x)=H(x)-s-t*(G(x)-s)
  B(x)= FG(x)+FH(x)>=0 ? FG(x).*FH(x) : -0.5*(FG(x).^2+FH(x).^2)
  nl_constraint(x)=[mod.mp.c(x);B(x);G(x)+r;H(x)+r]
 else
  println("No matching relaxation name. Default : No relaxation. Try : SS, KDB, KS or Butterfly")
  nl_constraint(x)=[mod.mp.c(x);G(x).*H(x);G(x);H(x)]
 end

 lcon=[mod.mp.meta.lcon;-Inf*ones(mod.nb_comp);zeros(mod.nb_comp*2)]
 ucon=[mod.mp.meta.ucon;zeros(mod.nb_comp);Inf*ones(mod.nb_comp*2)]
 y0=[mod.mp.meta.y0;zeros(3*mod.nb_comp)]

 #appel au constructeur NLP que l'on souhaite utiliser.
 nlp = ADNLPModel(mod.mp.f, mod.mp.meta.x0, lvar=mod.mp.meta.lvar, uvar=mod.mp.meta.uvar, y0=y0, c=nl_constraint, lcon=lcon, ucon=ucon)

 return nlp
end

#end du module
end

ext: jl


text:
module FunctionZeros
import SpecialFunctions
import Roots

export besselj_zero, bessely_zero

# Asymptotic formula for zeros of Bessel J function of order nu
"""
    besselj_zero_asymptotic(nu, n)

Asymptotic formula for the `n`th zero of the the Bessel J function of order `nu`.
"""
besselj_zero_asymptotic(nu, n) = pi * (n - 1 + nu / 2 + 3//4)

"""
    bessely_zero_asymptotic(nu, n)

Asymptotic formula for the `n`th zero of the the Bessel Y function of order `nu`.
"""
bessely_zero_asymptotic(nu, n) = pi * (n + nu / 2 - 3//4)

# Use the asymptotic values as starting values.
# These find the correct zeros even for n = 1,2,...
# Order 0 is 6 times slower and 50-100 times less accurate
# than higher orders, with other parameters constant.
"""
    besselj_zero(nu, n)

`n`th zero of the Bessel J function of order `nu`,
for `n` = `1,2,...`.
"""
besselj_zero(nu, n) = Roots.fzero((x) -> SpecialFunctions.besselj(nu, x),
                                  besselj_zero_asymptotic(nu, n); order=1)

"""
    bessely_zero(nu, n)

`n`th zero of the Bessel Y function of order `nu`,
for `n` = `1,2,...`.
"""
bessely_zero(nu, n) = Roots.fzero((x) -> SpecialFunctions.bessely(nu, x),
                                 bessely_zero_asymptotic(nu, n); order=1)

end # module FunctionZeros


text:
rosenbrock2d(x) =  (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
RANDOM_SEED = 809230
rng = MersenneTwister(RANDOM_SEED)

@testset "SimulatedAnnealing" begin

    # * Single-run tests for user-defined functions
    @test begin
        dimension = 2
        ground_truth = 0.0
        val = SimulatedAnnealing(
            rosenbrock2d, -5.0, 5.0, dimension, rng;
            low_temp=5000
        )
        @show val
        isapprox(val.min, ground_truth, atol=1e-2)
    end

    # * Single-run tests for benchmark functions
    @test begin
        dimension = 3
        ground_truth = 0.0
        val = SimulatedAnnealing(
            Rosenbrock(), -5.0, 5.0, dimension, rng;
            low_temp=20000
        )
        @show val
        isapprox(val.min, ground_truth, atol=1e-2)
    end

    # * Tests for random seeds
    @test begin
        dimension = 3
        ground_truth = 0.0
        val = SimulatedAnnealing(
            Rosenbrock(), -5.0, 5.0, dimension, rng;
            low_temp=20000
        )
        @show val
        isapprox(val.min, ground_truth, atol=1e-2)
    end
end

@testset "GeneralizedSimulatedAnnealing" begin
    # * Single-run tests for user-defined functions
    @test begin
        dimension = 2
        ground_truth = 0.0
        val = GeneralSimulatedAnnealing(
            rosenbrock2d, -5.0, 5.0, dimension, rng;
            low_temp=35000
        )
        @show val
        isapprox(val.min, ground_truth, atol=1e-2)
    end

    # * Single-run tests for different Tsallis' parameters
    @test begin
        dimension = 2
        ground_truth = 0.0
        val = GeneralSimulatedAnnealing(
            rosenbrock2d,
            -5.0,
            5.0,
            dimension,
            rng;
            low_temp=30000,
            qv=2.5,
            qa=1.0000001
        )
        @show val
        isapprox(val.min, ground_truth, atol=1e-2)
    end

    # * Single-run tests for benchmark functions
    @test begin
        dimension = 3
        ground_truth = 0.0
        val = GeneralSimulatedAnnealing(
            Ackley(),
            -3.0,
            3.0,
            dimension,
            rng;
            low_temp=10000
        )
        @show val
        isapprox(val.min, ground_truth, atol=1e-2)
    end

    # * Test for random seeds
    @test begin
        dimension = 3
        ground_truth = 0.0
        val = GeneralSimulatedAnnealing(
            Ackley(),
            -3.0,
            3.0,
            dimension,
            rng;
            low_temp=10000
        )
        @show val
        isapprox(val.min, ground_truth, atol=1e-2)
    end
end


text:
using Test

using ARCHModels
using GLM
using DataFrames
using StableRNGs

T = 10^4;

@testset "lgamma" begin
    @test ARCHModels.lgamma(1.0f0) == 0.0f0
end

@testset "TGARCH" begin
    @test ARCHModels.nparams(TGARCH{1, 2, 3}) == 7
    @test ARCHModels.presample(TGARCH{1, 2, 3}) == 3
    spec = TGARCH{1,1,1}([1., .05, .9, .01]);
    str = sprint(show, spec)
    if VERSION < v"1.5.5"
        @test startswith(str, "TGARCH{1,1,1} specification.\n\n─────────────────────────────────\n               ω    γ₁   β₁    α₁\n─────────────────────────────────\nParameters:  1.0  0.05  0.9  0.01\n─────────────────────────────────\n")
    else
        @test startswith(str, "TGARCH{1, 1, 1} specification.\n\n─────────────────────────────────\n               ω    γ₁   β₁    α₁\n─────────────────────────────────\nParameters:  1.0  0.05  0.9  0.01\n─────────────────────────────────\n")
    end
    am = simulate(spec, T, rng=StableRNG(1));
    am = selectmodel(TGARCH, am.data; meanspec=NoIntercept(), show_trace=true, maxlags=2)
    @test all(isapprox.(coef(am), [1.3954654215590847,
                                   0.06693040956623193,
                                   0.8680818765441008,
                                   0.006665140784151278], rtol=1e-4))
   #everything below is just pure GARCH, in fact
    spec = GARCH{1, 1}([1., .9, .05])
    am0 = simulate(spec, T; rng=StableRNG(1));
    am00 = deepcopy(am0)
    am00.data .= 0.
    simulate!(am00, rng=StableRNG(1))
    @test all(am00.data .== am0.data)
    am00 = simulate(am0; rng=StableRNG(1))
    @test all(am00.data .== am0.data)
    am000 = simulate(am0, nobs(am0); rng=StableRNG(1))
    @test all(am000.data .== am0.data)
    am = selectmodel(GARCH, am0.data; meanspec=NoIntercept(), show_trace=true)
    @test isfitted(am) == true
    @test all(isapprox.(coef(am), [1.116707484875346,
                                   0.8920705288828562,
                                   0.05103227915762242], rtol=1e-4))
    @test all(isapprox.(stderror(am), [ 0.22260057264313066,
                                        0.016030182299773734,
                                        0.006460941055580745], rtol=1e-3))
    @test sum(volatilities(am0)) ≈ 44285.00568611553
    @test sum(abs, residuals(am0)) ≈ 7964.585890843087
    @test sum(abs, residuals(am0, standardized=false)) ≈ 35281.71207401529
    am2 = UnivariateARCHModel(spec, am0.data)
    @test isfitted(am2) == false
    io = IOBuffer()
    str = sprint(io -> show(io, am2))
    if VERSION < v"1.5.5"
        @test startswith(str, "\nTGARCH{0,1,1}")
    else
        @test startswith(str, "\nGARCH{1, 1}")
    end
    fit!(am2)
    @test isfitted(am2) == true
    io = IOBuffer()
    str = sprint(io -> show(io, am2))
    if VERSION < v"1.5.5"
        @test startswith(str, "\nTGARCH{0,1,1}")
    else
        @test startswith(str, "\nGARCH{1, 1}")
    end
    am3 = fit(am2)
    @test isfitted(am3) == true
    @test all(am2.spec.coefs .== am.spec.coefs)
    @test all(am3.spec.coefs .== am2.spec.coefs)
end
@testset "ARCH" begin
    spec = ARCH{2}([1., .3, .4]);
    am = simulate(spec, T; rng=StableRNG(1));
    @test selectmodel(ARCH, am.data).spec.coefs == fit(ARCH{2}, am.data).spec.coefs
    spec = ARCH{0}([1.]);
    am = simulate(spec, T, rng=StableRNG(1));
    fit!(am)
    @test all(isapprox.(coef(am),  0.991377950108106, rtol=1e-4))

end

@testset "EGARCH" begin
    @test ARCHModels.nparams(EGARCH{1, 2, 3}) == 7
    @test ARCHModels.presample(EGARCH{1, 2, 3}) == 3
    am = simulate(EGARCH{1, 1, 1}([.1, 0., .9, .1]), T; meanspec=Intercept(3), rng=StableRNG(1))
    am7 = selectmodel(EGARCH, am.data; maxlags=2, show_trace=true)
    @test all(isapprox(coef(am7), [ 0.1240152087585493,
                                   -0.010544394266072957,
                                    0.874501604519596,
                                    0.10762246065941368,
                                    3.0008464829419053], rtol=1e-4))

    @test coefnames(EGARCH{2, 2, 2}) == ["ω", "γ₁", "γ₂", "β₁", "β₂", "α₁", "α₂"]
    @test_throws Base.ErrorException predict.(am7, :variance, 1:3)
end
@testset "StatisticalModel" begin
    #not implemented: adjr2, deviance, mss, nulldeviance, r2, rss, weights
    spec = GARCH{1, 1}([1., .9, .05])
    am = simulate(spec, T; rng=StableRNG(1))
    fit!(am)
    @test loglikelihood(am) ==  ARCHModels.loglik!(Float64[],
                                                                Float64[],
                                                                Float64[],
                                                                Float64[],
                                                                typeof(spec),
                                                                StdNormal{Float64},
                                                                NoIntercept{Float64}(),
                                                                am.data,
                                                                spec.coefs
                                                                )
    @test nobs(am) == T
    @test dof(am) == 3
    @test coefnames(GARCH{1, 1}) == ["ω", "β₁", "α₁"]
    @test aic(am) ≈ 57949.19500673284 rtol=1e-4
    @test bic(am) ≈ 57970.82602784877 rtol=1e-4
    @test aicc(am) ≈ 57949.19740769323 rtol=1e-4

    @test all(coef(am) .== am.spec.coefs)
    @test all(isapprox(confint(am), [ 0.680418   1.553;
                                      0.860652   0.923489;
                                      0.0383691  0.0636955],
                       rtol=1e-4)
                       )
    @test all(isapprox(informationmatrix(am; expected=false)/T, [ 0.125032   2.33319   2.07012;
                                                                2.33319   44.6399   40.8553;
                                                                2.07012   40.8553   41.2192],
                       rtol=1e-4)
                       )
    @test_throws ErrorException informationmatrix(am)
    @test all(isapprox(score(am), [0. 0. 0.], atol=1e-3))
    @test islinear(am::UnivariateARCHModel) == false
    @test predict(am) ≈ 4.296827552671104
    @test predict(am, :variance) ≈ 18.46272701739355
    @test predict(am, :return) == 0.0
    @test predict(am, :VaR) ≈ 9.995915642276554
    for what in [:return, :variance]
        @test predict.(am, what, 1:3) == [predict(am, what, h) for h in 1:3]
    end
    @test_throws Base.ErrorException predict.(am, :VaR, 1:3)
    @test_throws Base.ErrorException predict.(am, :volatility, 1:3)
end

@testset "MeanSpecs" begin
    spec = GARCH{1, 1}([1., .9, .05])
    am = simulate(spec, T; meanspec=Intercept(0.), rng=StableRNG(1))
    fit!(am)
    @test all(isapprox(coef(am), [1.1176635890968043,
                                  0.8919906787166815,
                                  0.05106346071866704,
                                  0.00952591461710004], rtol=1e-4))
    @test ARCHModels.coefnames(Intercept(0.)) == ["μ"]
    @test ARCHModels.nparams(Intercept) == 1
    @test ARCHModels.presample(Intercept(0.)) == 0
    @test ARCHModels.constraints(Intercept{Float64}, Float64) == (-Float64[Inf], Float64[Inf])
    @test typeof(NoIntercept()) == NoIntercept{Float64}
    @test ARCHModels.coefnames(NoIntercept()) == []
    @test ARCHModels.constraints(NoIntercept{Float64}, Float64) == (Float64[], Float64[])
    @test ARCHModels.nparams(NoIntercept) == 0
    @test ARCHModels.presample(NoIntercept()) == 0
    @test ARCHModels.uncond(NoIntercept()) == 0
    @test mean(zeros(5), zeros(5), zeros(5), zeros(5), NoIntercept(), zeros(5), 4) == 0.
    ms = ARMA{2, 2}([1., .5, .2, -.1, .3])
    @test ARCHModels.nparams(typeof(ms)) == length(ms.coefs)
    @test ARCHModels.presample(ms) == 2
    @test ARCHModels.coefnames(ms) == ["c", "φ₁", "φ₂", "θ₁", "θ₂"]
    spec = GARCH{1, 1}([1., .9, .05])
    am = simulate(spec, T; meanspec=ms, rng=StableRNG(1))
    fit!(am)
    @test all(isapprox(coef(am), [ 1.1375727511714622,
                                   0.8903853180079492,
                                   0.05158067874765809,
                                   1.0091192373639755,
                                   0.482666588367849,
                                   0.21802258440272837,
                                  -0.08390300941364812,
                                   0.28868236034111855], rtol=1e-4))
    @test predict(am, :return) ≈ 2.335436537249963 rtol = 1e-6
    am = selectmodel(ARCH, BG96;  meanspec=AR, maxlags=2);
    @test all(isapprox(coef(am), [0.1191634087516343,
                                  0.31568628680702837,
                                  0.18331803992648235,
                                 -0.006857008709781168,
                                  0.035836278501164005], rtol=1e-4))
    @test typeof(Regression([1 2; 3 4])) == Regression{2, Float64}
    @test typeof(Regression([1. 2.; 3. 4.])) == Regression{2, Float64}
    @test typeof(Regression{Float32}([1 2; 3 4])) == Regression{2, Float32}
    @test typeof(Regression([1 2; 3 4])) == Regression{2, Float64}
    @test typeof(Regression([1, 2], [1 2; 3 4.0f0])) ==  Regression{2, Float32}
    @test typeof(Regression([1, 2.], [1 2; 3 4.0f0])) ==  Regression{2, Float64}
    @test typeof(Regression([1], [1, 2, 3, 4.0f0])) ==  Regression{1, Float32}
    @test typeof(Regression([1, 2, 3, 4.0f0])) ==  Regression{1, Float32}
    @test ARCHModels.nparams(Regression{2, Float64}) == 2

    rng = StableRNG(1)
    beta = [1, 2]
    reg = Regression(beta, rand(rng, 2000, 2))
    u = randn(rng, 2000)*.1
    y = reg.X*reg.coefs+u
    @test ARCHModels.coefnames(reg) == ["β₀", "β₁"]
    @test ARCHModels.presample(reg) == 0
    @test ARCHModels.constraints(typeof(reg), Float64) == ([-Inf, -Inf], [Inf, Inf])
    @test all(isapprox(ARCHModels.startingvals(reg, y),
        [0.992361089980835, 2.003646964507331], rtol=1e-4))
    @test ARCHModels.uncond(reg) === 0.
    am = simulate(GARCH{1, 1}([1., .9, .05]), 2000; meanspec=reg, warmup=0, rng=StableRNG(1))
    fit!(am)
    @test_throws Base.ErrorException predict(am, :return)

    @test all(isapprox(coef(am), [1.098632569628791,
                                  0.8866288812154145,
                                  0.05770241980639491,
                                  0.7697476790102007,
                                  2.403750061921962], rtol=1e-4))
    am = simulate(GARCH{1, 1}([1., .9, .05]), 1999; meanspec=reg, warmup=0, rng=StableRNG(1))
    @test predict(am, :return) ≈ 2.3760239544958175
    data = DataFrame(X=ones(1974), Y=BG96)
    model = lm(@formula(Y ~ -1 + X), data)
    am = fit(GARCH{1, 1}, model)
    @test all(isapprox(coef(am), coef(fit(GARCH{1, 1}, BG96, meanspec=Intercept)), rtol=1e-4))
    @test coefnames(am)[end] == "X"
    @test all(isapprox(coef(am), coef(fit(GARCH{1, 1}, model.model)), rtol=1e-4))
    @test sum(coef(fit(ARMA{1, 1}, BG96))) ≈ 0.21595383060382695
    @test sum(coef(selectmodel(ARMA, BG96; minlags=2, maxlags=3))) ≈ 0.25404112526031397
end

@testset "VaR" begin
    am = fit(GARCH{1, 1}, BG96)
    @test sum(VaRs(am)) ≈ 2077.0976454790807
end
@testset "Errors" begin
    #with unconditional as presample:
    #@test_warn "Fisher" stderror(UnivariateARCHModel(GARCH{3, 0}([1., .1, .2, .3]), [.1, .2, .3, .4, .5, .6, .7]))
    #@test_warn "non-positive" stderror(UnivariateARCHModel(GARCH{3, 0}([1., .1, .2, .3]), -5*[.1, .2, .3, .4, .5, .6, .7]))

    # the following are temporarily disabled while we use FiniteDiff for Hessians:
    #@test_logs (:warn, "Fisher information is singular; vcov matrix is inaccurate.") stderror(UnivariateARCHModel(GARCH{1, 0}( [1.0, .1]), [0., 1.]))
    #@test_logs (:warn, "non-positive variance encountered; vcov matrix is inaccurate.") stderror(UnivariateARCHModel(GARCH{1, 0}( [1.0, .1]), [1., 1.]))
    e = @test_throws ARCHModels.NumParamError ARCHModels.loglik!(Float64[], Float64[], Float64[], Float64[], GARCH{1, 1}, StdNormal{Float64},
                                                     NoIntercept{Float64}(), zeros(T),
                                                     [0., 0., 0., 0.]
                                                     )
    str = sprint(showerror, e.value)
    @test startswith(str, "incorrect number of parameters")
    @test_throws ARCHModels.NumParamError GARCH{1, 1}([.1])
    e = @test_throws ErrorException predict(UnivariateARCHModel(GARCH{0, 0}([1.]), zeros(10)), :blah)
    str = sprint(showerror, e.value)
    @test startswith(str, "Prediction target blah unknown")
    @test_throws ARCHModels.NumParamError ARMA{1, 1}([1.])
    @test_throws ARCHModels.NumParamError Intercept([1., 2.])
    @test_throws ARCHModels.NumParamError NoIntercept([1.])
    @test_throws ARCHModels.NumParamError StdNormal([1.])
    @test_throws ARCHModels.NumParamError StdT([1., 2.])
    @test_throws ARCHModels.NumParamError StdSkewT([2.])
    @test_throws ARCHModels.NumParamError StdGED([1., 2.])
    @test_throws ARCHModels.NumParamError Regression([1], [1 2; 3 4])
    at = zeros(10)
    data = rand(StableRNG(1), 10)
    reg = Regression(data[1:5])
    @test_throws ErrorException mean(at, at, at, data, reg, [0.], 6)
end

@testset "Distributions" begin
    a=rand(StableRNG(1), StdT(3))

    b=rand(StableRNG(1), StdT(3), 1)[1]
    @test a==b

    @test rand(StableRNG(1), StdNormal()) ≈ -0.5325200748641231
    @testset "Gaussian" begin
        data = rand(StableRNG(1), T)
        @test typeof(StdNormal())==typeof(StdNormal(Float64[]))
        @test fit(StdNormal, data).coefs == Float64[]
        @test coefnames(StdNormal) == String[]
        @test ARCHModels.distname(StdNormal) == "Gaussian"
        @test quantile(StdNormal(), .05) ≈ -1.6448536269514724
        @test ARCHModels.constraints(StdNormal{Float64}, Float64) == (Float64[], Float64[])
    end
    @testset "Student" begin
        data = rand(StableRNG(1), StdT(4), T)
        spec = GARCH{1, 1}([1., .9, .05])
        @test fit(StdT, data).coefs[1] ≈ 4.027241574961463 rtol=1e-4
        @test coefnames(StdT) == ["ν"]
        @test ARCHModels.distname(StdT) == "Student's t"
        @test quantile(StdT(3), .05) ≈ -1.3587150125838563
        datat = simulate(spec, T; dist=StdT(4), rng=StableRNG(1)).data
        datam = simulate(spec, T; dist=StdT(4), meanspec=Intercept(3), rng=StableRNG(1)).data
        am4 = selectmodel(GARCH, datat; dist=StdT, meanspec=NoIntercept{Float64}(), show_trace=true)
        am5 = selectmodel(GARCH, datam; dist=StdT, show_trace=true)
        @test coefnames(am5) == ["ω", "β₁", "α₁", "ν", "μ"]
        @test all(coeftable(am4).cols[2] .== stderror(am4))
        @test all(isapprox(coef(am4), [0.7400801158793416,
                                       0.9185237789418755,
                                       0.04062109019124998,
                                       4.267591319191636], rtol=1e-4))
        @test all(isapprox(coef(am5), [0.7402313332187656,
                                       0.9184701329706532,
                                       0.040654605524829274,
                                       4.269395399912836,
                                       3.035298011833236], rtol=1e-4))
    end
    @testset "HansenSkewedT" begin
       data = rand(StableRNG(1), StdSkewT(4,-0.3), T)
       spec = GARCH{1, 1}([1., .9, .05])
       c = fit(StdSkewT, data).coefs
       @test c[1] ≈ 3.990671630456716 rtol=1e-4
       @test c[2] ≈ -0.3136773995478942 rtol=1e-4
       @test typeof(StdSkewT(3,0)) == typeof(StdSkewT(3.,0)) == typeof(StdSkewT([3,0.0]))
       @test coefnames(StdSkewT) == ["ν", "λ"]
       @test ARCHModels.nparams(StdSkewT) == 2
       @test ARCHModels.distname(StdSkewT) == "Hansen's Skewed t"
       @test ARCHModels.constraints(StdNormal{Float64}, Float64) == (Float64[], Float64[])
       @test quantile(StdSkewT(3,0), 0.5) == 0
       @test quantile(StdSkewT(3,0), .05) ≈ -1.3587150125838563
       @test ARCHModels.constraints(StdSkewT{Float64}, Float64) == (Float64[20/10, -one(Float64)], Float64[Inf,one(Float64)])
       dataskt = simulate(spec, T; dist=StdSkewT(4,-0.3), rng=StableRNG(1)).data
       datam = simulate(spec, T; dist=StdSkewT(4,-0.3), meanspec=Intercept(3), rng=StableRNG(1)).data
       am4 = selectmodel(GARCH, dataskt; dist=StdSkewT, meanspec=NoIntercept{Float64}(), show_trace=true)
       am5 = selectmodel(GARCH, datam; dist=StdSkewT, show_trace=true)
       @test coefnames(am5) == ["ω", "β₁", "α₁", "ν", "λ", "μ"]
       @test all(coeftable(am4).cols[2] .== stderror(am4))
       @test all(isapprox(coef(am4), [ 1.0123398035363282,
                                       0.9010308454299863,
                                       0.042335307040165894,
                                       4.24455990918083,
                                      -0.3115002211205442], rtol=1e-4))
       @test all(isapprox(coef(am5), [ 1.0151845148616474,
                                       0.9009908899358181,
                                       0.04243949895951436,
                                       4.241005415020919,
                                      -0.3124667515252298,
                                       2.9931917146031144], rtol=1e-4))
    end
    @testset "GED" begin
        @test typeof(StdGED(3)) == typeof(StdGED(3.)) == typeof(StdGED([3.]))
        data = rand(StableRNG(1), StdGED(1), T)
        @test fit(StdGED, data).coefs[1] ≈ 1.0062771307312837 rtol=1e-4
        @test coefnames(StdGED) == ["p"]
        @test ARCHModels.nparams(StdGED) == 1
        @test ARCHModels.distname(StdGED) == "GED"
        @test quantile(StdGED(1), .05) ≈ -1.6281735335151468
    end
    @testset "Standardized" begin
        using Distributions
        @test eltype(StdNormal{Float64}()) == Float64
        MyStdT=Standardized{TDist}
        @test typeof(MyStdT([1.])) == typeof(MyStdT(1.))
        @test ARCHModels.logconst(MyStdT, [0]) == 0.
        @test coefnames(MyStdT{Float64}) == ["ν"]
        @test ARCHModels.distname(MyStdT{Float64}) == "TDist"
        @test all(isapprox.(ARCHModels.startingvals(MyStdT, [0.]), eps()))
        @test quantile(MyStdT(3.), .1) ≈ quantile(StdT(3.), .1)
        ARCHModels.startingvals(::Type{<:MyStdT}, data::Vector{T}) where T = T[3.]
        am = simulate(GARCH{1, 1}([1, 0.9, .05]), 1000, dist=MyStdT(3.); rng=StableRNG(1))
        @test  loglikelihood(fit(am)) ≈ -2543.7820751346744
    end
end
@testset "tests" begin
    am = fit(GARCH{1, 1}, BG96)
    LM = ARCHLMTest(am)
    @test pvalue(LM) ≈ 0.1139758664282619
    str = sprint(show, LM)
    @test startswith(str, "ARCH LM test for conditional heteroskedasticity")
    @test ARCHModels.testname(LM) == "ARCH LM test for conditional heteroskedasticity"

    vars = VaRs(am, 0.01)
    DQ = DQTest(BG96, VaRs(am), 0.01)
    @test pvalue(DQ) ≈ 2.3891461144184955e-11
    str = sprint(show, DQ)
    @test startswith(str, "Engle and Manganelli's (2004) DQ test (out of sample)")
    @test ARCHModels.testname(DQ) == "Engle and Manganelli's (2004) DQ test (out of sample)"
end
@testset "multivariate" begin
    am1 = fit(DCC, DOW29[:, 1:2])
    am2 = fit(DCC, DOW29[:, 1:2]; method=:twostep)
    am3 = MultivariateARCHModel(DCC{1, 1}([1. 0.; 0. 1.], [0., 0.], [GARCH{1, 1}([1., 0., 0.]), GARCH{1, 1}([1., 0., 0.])]), DOW29[:, 1:2]) # not fitted
    am4 = fit(DCC, DOW29[1:20, 1:29]) # shrinkage n<p
    @test all(fit(am1).spec.coefs .== am1.spec.coefs)
    @test all(isapprox(am1.spec.coefs, [0.8912884521017908, 0.05515419379547665], rtol=1e-3))
    @test all(isapprox(am2.spec.coefs,    [0.8912161306136979, 0.055139392936998946], rtol=1e-3))
    @test all(isapprox(am4.spec.coefs, [0.8935938309400944, 6.938893903907228e-18], atol=1e-3))
    @test all(isapprox(stderror(am1)[1:2], [0.0434344187103969, 0.020778846682313102], rtol=1e-3))
    @test all(isapprox(stderror(am2)[1:2], [0.030405542205923865, 0.014782869078355866], rtol=1e-4))
    @test all(isapprox(predict(am1; what=:correlation)[:], [1.0, 0.4365129466277069, 0.4365129466277069, 1.0], rtol=1e-4))
    @test all(isapprox(predict(am1; what=:covariance)[:], [6.916591739333349, 1.329392154000225, 1.329392154000225,  1.340972349032465], rtol=1e-4))
    @test_throws ErrorException predict(am1; what=:bla)
    @test residuals(am1)[1, 1] ≈ 0.5107042609407892
    @test_throws ErrorException fit(DCC, DOW29; method=:bla)
    @test_throws ARCHModels.NumParamError DCC{1, 1}([1. 0.; 0. 1.], [1., 0., 0.], [GARCH{1, 1}([1., 0., 0.]), GARCH{1, 1}([1., 0., 0.])])
    @test_throws AssertionError DCC{1, 1}([1. 0.; 0. 1.], [0., 0.], [GARCH{1, 1}([1., 0., 0.]), GARCH{1, 1}([1., 0., 0.])]; method=:bla)
    @test coefnames(am1) == ["β₁", "α₁", "ω₁", "β₁₁", "α₁₁", "μ₁", "ω₂", "β₁₂", "α₁₂", "μ₂"]
    @test ARCHModels.nparams(DCC{1, 1}) == 2
    ARCHModels.nparams(DCC{1, 1, GARCH{1, 1}, Float64, 2}) == 8
    @test ARCHModels.presample(DCC{1, 2, GARCH{3, 4}}) == 4
    @test ARCHModels.presample(DCC{1, 2, GARCH{3, 4, Float64}, Float64, 2}) == 4
    io = IOBuffer()
    str = sprint(io -> show(io, am1))
    @test startswith(str, "\n2-dim")
    str = sprint(io -> show(io, am3))
    @test startswith(str, "\n2-dim")
    str = sprint(io -> show(io, am3.spec))
    @test startswith(str, "DCC{1, 1")
    str = sprint(io -> show(IOContext(io, :se=>true), am1))
    @test occursin("Std.Error", str)
    @test_throws ErrorException fit(DCC, DOW29[1:11, :]) # shrinkage requires n>=12
    @test loglikelihood(am1) ≈ -9810.905799585276

    @test ARCHModels.nparams(MultivariateStdNormal) == 0
    @test typeof(MultivariateStdNormal{Float64, 3}()) == typeof(MultivariateStdNormal{Float64, 3}(Float64[]))
    @test typeof(MultivariateStdNormal(Float64, 3)) == typeof(MultivariateStdNormal{Float64, 3}(Float64[]))
    @test typeof(MultivariateStdNormal(Float64[], 3)) == typeof(MultivariateStdNormal{Float64, 3}(Float64[]))
    @test typeof(MultivariateStdNormal{Float64}(3)) == typeof(MultivariateStdNormal{Float64, 3}(Float64[]))
    @test typeof(MultivariateStdNormal(3)) == typeof(MultivariateStdNormal{Float64, 3}(Float64[]))
    @test all(isapprox(rand(StableRNG(1), MultivariateStdNormal(2)), [-0.5325200748641231,  0.098465514284785], rtol=1e-6))
    @test coefnames(MultivariateStdNormal) == String[]
    @test ARCHModels.distname(MultivariateStdNormal) == "Multivariate Normal"

    am = am1
    am.spec.coefs .= [.7, .2]
    ams  = simulate(am; rng=StableRNG(1))
    @test isfitted(ams) == false
    fit!(ams)
    @test isfitted(ams) == true
    @test all(isapprox(ams.spec.coefs, [0.6611103068430052, 0.23089471530783906], rtol=1e-4))
    simulate!(ams; rng=StableRNG(2))
    @test ams.fitted == false
    fit!(ams)
    @test all(isapprox(ams.spec.coefs, [0.6660369039914371, 0.2329752007155509], rtol=1e-4))
    amc = fit(DCC{1, 2, GARCH{3, 2}}, DOW29[:, 1:4]; meanspec=AR{3})
    ams = simulate(amc, T; rng=StableRNG(1))
    fit!(ams)
    @test all(isapprox(ams.meanspec[1].coefs, [-0.1040426570178552, 0.03639191550146291, 0.033657970110476075, -0.020300480179225668], rtol=1e-4))
    ame = fit(DCC{1, 2, EGARCH{1, 1, 1}}, DOW29[:, 1:4])
    ams = simulate(ame, T; rng=StableRNG(1))
    fit!(ams)
    @test all(isapprox(ams.spec.univariatespecs[1].coefs, [0.05335407349997172, -0.08008165178490954,  0.9627467601623543,  0.22652855417695117], rtol=1e-4))
    ccc = fit(CCC, DOW29[:, 1:4])
    @test dof(ccc) == 16
    @test ccc.spec.R[1, 2] ≈ 0.37095654552885643
    @test isapprox(stderror(ccc)[1], 0.06298215515406534, rtol=1e-3)
    cccs = simulate(ccc, T; rng=StableRNG(1))
    @test  cccs.data[end, 1] ≈ -0.8530862593689736
    @test coefnames(ccc) == ["ω₁", "β₁₁", "α₁₁", "μ₁", "ω₂", "β₁₂", "α₁₂", "μ₂", "ω₃", "β₁₃", "α₁₃", "μ₃", "ω₄", "β₁₄", "α₁₄", "μ₄"]
    io = IOBuffer()
    str = sprint(io -> show(io, ccc))
    @test startswith(str, "\n4-dim")
    io = IOBuffer()
    str = sprint(io -> show(io, ccc.spec))
    @test startswith(str, "DCC{0, 0")
end


text:
module MimiDICE2016

using Mimi
using ExcelReaders

include("helpers.jl")
include("parameters.jl")

include("marginaldamage.jl")

include("components/totalfactorproductivity_component.jl")
include("components/grosseconomy_component.jl")
include("components/emissions_component.jl")
include("components/co2cycle_component.jl")
include("components/radiativeforcing_component.jl")
include("components/climatedynamics_component.jl")
include("components/damages_component.jl")
include("components/neteconomy_component.jl")
include("components/welfare_component.jl")

export constructdice, getdiceexcel, getdicegams

const model_years = 2015:5:2510

function constructdice(p)

    m = Model()
    set_dimension!(m, :time, model_years)

    add_comp!(m, totalfactorproductivity, :totalfactorproductivity)
    add_comp!(m, grosseconomy, :grosseconomy)
    add_comp!(m, emissions, :emissions)
    add_comp!(m, co2cycle, :co2cycle)
    add_comp!(m, radiativeforcing, :radiativeforcing)
    add_comp!(m, climatedynamics, :climatedynamics)
    add_comp!(m, damages, :damages)
    add_comp!(m, neteconomy, :neteconomy)
    add_comp!(m, welfare, :welfare)

    # TFP COMPONENT
    set_param!(m, :totalfactorproductivity, :a0, p[:a0])
	set_param!(m, :totalfactorproductivity, :ga0, p[:ga0])
    set_param!(m, :totalfactorproductivity, :dela, p[:dela])

    # GROSS ECONOMY COMPONENT
    set_param!(m, :grosseconomy, :l, p[:l])
    set_param!(m, :grosseconomy, :gama, p[:gama])
    set_param!(m, :grosseconomy, :dk, p[:dk])
    set_param!(m, :grosseconomy, :k0, p[:k0])

    connect_param!(m, :grosseconomy, :AL, :totalfactorproductivity, :AL)
    connect_param!(m, :grosseconomy, :I, :neteconomy, :I)

    # EMISSIONS COMPONENT
    set_param!(m, :emissions, :gsigma1, p[:gsigma1])
    set_param!(m, :emissions, :dsig, p[:dsig])
	set_param!(m, :emissions, :eland0, p[:eland0])
	set_param!(m, :emissions, :deland, p[:deland])
	set_param!(m, :emissions, :e0, p[:e0])
    set_param!(m, :emissions, :MIU, p[:MIU])
    set_param!(m, :emissions, :cca0, p[:cca0])
	set_param!(m, :emissions, :cumetree0, p[:cumetree0])
    connect_param!(m, :emissions, :YGROSS, :grosseconomy, :YGROSS)

    # CO2 CYCLE COMPONENT
    set_param!(m, :co2cycle, :mat0, p[:mat0])
    set_param!(m, :co2cycle, :mu0, p[:mu0])
    set_param!(m, :co2cycle, :ml0, p[:ml0])
    set_param!(m, :co2cycle, :b12, p[:b12])
    set_param!(m, :co2cycle, :b23, p[:b23])
	set_param!(m, :co2cycle, :mateq, p[:mateq])
	set_param!(m, :co2cycle, :mleq, p[:mleq])
	set_param!(m, :co2cycle, :mueq, p[:mueq])
    connect_param!(m, :co2cycle, :E, :emissions, :E)

    # RADIATIVE FORCING COMPONENT
    set_param!(m, :radiativeforcing, :fex0, p[:fex0])
    set_param!(m, :radiativeforcing, :fex1, p[:fex1])
    set_param!(m, :radiativeforcing, :fco22x, p[:fco22x])
    set_param!(m, :radiativeforcing, :eqmat, p[:eqmat])
    connect_param!(m, :radiativeforcing, :MAT, :co2cycle, :MAT)

    # CLIMATE DYNAMICS COMPONENT
    set_param!(m, :climatedynamics, :fco22x, p[:fco22x])
    set_param!(m, :climatedynamics, :t2xco2, p[:t2xco2])
    set_param!(m, :climatedynamics, :tatm0, p[:tatm0])
    set_param!(m, :climatedynamics, :tocean0, p[:tocean0])
    set_param!(m, :climatedynamics, :c1, p[:c1])
    set_param!(m, :climatedynamics, :c3, p[:c3])
    set_param!(m, :climatedynamics, :c4, p[:c4])
    connect_param!(m, :climatedynamics, :FORC, :radiativeforcing, :FORC)

    # DAMAGES COMPONENT
    set_param!(m, :damages, :a1, p[:a1])
    set_param!(m, :damages, :a2, p[:a2])
    set_param!(m, :damages, :a3, p[:a3])
    connect_param!(m, :damages, :TATM, :climatedynamics, :TATM)
    connect_param!(m, :damages, :YGROSS, :grosseconomy, :YGROSS)

    # NET ECONOMY COMPONENT
    set_param!(m, :neteconomy, :MIU, p[:MIU])
    set_param!(m, :neteconomy, :expcost2, p[:expcost2])
    set_param!(m, :neteconomy, :pback, p[:pback])
	set_param!(m, :neteconomy, :gback, p[:gback])
    set_param!(m, :neteconomy, :S, p[:S])
    set_param!(m, :neteconomy, :l, p[:l])
    connect_param!(m, :neteconomy, :YGROSS, :grosseconomy, :YGROSS)
    connect_param!(m, :neteconomy, :DAMAGES, :damages, :DAMAGES)
	connect_param!(m, :neteconomy, :SIGMA, :emissions, :SIGMA)

    # WELFARE COMPONENT
    set_param!(m, :welfare, :l, p[:l])
    set_param!(m, :welfare, :elasmu, p[:elasmu])
    set_param!(m, :welfare, :rr, p[:rr])
    set_param!(m, :welfare, :scale1, p[:scale1])
    set_param!(m, :welfare, :scale2, p[:scale2])
    connect_param!(m, :welfare, :CPC, :neteconomy, :CPC)

    return m

end

function getdiceexcel(;datafile = joinpath(dirname(@__FILE__), "..", "data", "DICE2016R-090916ap-v2.xlsm"))
    params = getdice2016excelparameters(datafile)

    m = constructdice(params)

    return m
end

function getdicegams(;datafile = joinpath(dirname(@__FILE__), "..", "data", "DICE2016_IAMF_Parameters.xlsx"))
    params = getdice2016gamsparameters(datafile)

    m = constructdice(params)

    return m
end

# get_model function for standard Mimi API: use the Excel version
get_model = getdiceexcel

end # module


text:
#
# Optimal layout generator.
#

# Layout constraints.

struct Spec
    line_width::Int
    break_cost::Int
    spill_cost::Int
end

# Maps a starting column to the optimal layout and its cost.

mutable struct Rope <: AbstractVector{Tuple{Int,Int,Int,Int,Layout}}
    len::Int
    cols::Vector{Int}
    spans::Vector{Int}
    icepts::Vector{Int}
    slopes::Vector{Int}
    lts::Vector{Layout}

    cache::Dict{Spec, Dict{Layout, Rope}}

    Rope() = new(0, Int[], Int[], Int[], Int[], Layout[], Dict())
    Rope(::Nothing) = new(1, [0], [0], [0], [0], [literal("")], Dict())
end

size(rope::Rope) = (rope.len,)

IndexStyle(::Type{Rope}) = IndexLinear()

@inline function getindex(rope::Rope, k::Int)
    @boundscheck checkbounds(rope, k)
    @inbounds t = (rope.cols[k], rope.spans[k], rope.icepts[k], rope.slopes[k], rope.lts[k])
    t
end

best(rope::Rope) = rope.lts[1]

function extend!(rope::Rope, col::Int, span::Int, icept::Int, slope::Int, lt::Layout)
    @assert rope.len == 0 && col == 0 || rope.len > 0 && col > rope.cols[end]
    if rope.len > 0
        col′, span′, icept′, slope′, lt′ = rope[end]
        lt.blk == lt′.blk && lt.args == lt′.args &&
        span == span′ &&
        icept == icept′ + (col - col′) * slope′ &&
        slope == slope′ && return
    end
    rope.len += 1
    push!(rope.cols, col)
    push!(rope.spans, span)
    push!(rope.icepts, icept)
    push!(rope.slopes, slope)
    push!(rope.lts, lt)
    nothing
end

# Generator state.

mutable struct Formatter
    spec::Spec
    empty::Rope

    Formatter(spec::Spec, empty::Rope) = new(spec, empty)

    Formatter(line_width::Int = DEFAULT_LINE_WIDTH,
              break_cost::Int = DEFAULT_BREAK_COST,
              spill_cost::Int = DEFAULT_SPILL_COST) =
        new(Spec(line_width, break_cost, spill_cost), Rope(nothing))

    Formatter(io::IO,
              break_cost::Int = DEFAULT_BREAK_COST,
              spill_cost::Int = DEFAULT_SPILL_COST) =
        let (h, w) = displaysize(io)
            Formatter(w-1, break_cost, spill_cost)
        end
end

best_fit(args...) = best(fit(args...))

fit(lt::Layout) = fit(Formatter(), lt)

fit(io::IO, lt::Layout) = fit(Formatter(io), lt)

fit(fmt::Formatter, lt::Layout) =
    fit(fmt, lt, fmt.empty)

function fit(fmt::Formatter, lt::Layout, tail::Rope)
    spec_cache = get!(tail.cache, fmt.spec) do
        Dict{Layout, Rope}()
    end
    get!(spec_cache, lt) do
        fit(fmt, lt.blk, lt.args, tail)
    end
end

function fit(fmt::Formatter, blk::LiteralBlock, args::Vector{Layout}, tail::Rope)
    line_width = fmt.spec.line_width
    spill_cost = fmt.spec.spill_cost
    w = blk.len
    lt = Layout(blk, args)
    rope = Rope()
    if w < line_width
        extend!(rope, 0, w, 0, 0, lt)
        extend!(rope, line_width - w, w, 0, spill_cost, lt)
    else
        extend!(rope, 0, w, (w - line_width) * spill_cost, spill_cost, lt)
    end
    return horizontal(fmt, rope, tail)
end

function fit(fmt::Formatter, ::HorizontalBlock, args::Vector{Layout}, tail::Rope)
    rope = tail
    k = length(args)
    while k > 0
        rope = fit(fmt, args[k], rope)
        k -= 1
    end
    return rope
end

function fit(fmt::Formatter, ::VerticalBlock, args::Vector{Layout}, tail::Rope)
    rope = tail
    k = length(args)
    while k > 0
        if k == length(args)
            rope = fit(fmt, args[k], rope)
        else
            rope = vertical(fmt, fit(fmt, args[k]), rope)
        end
        k -= 1
    end
    return rope
end

function fit(fmt::Formatter, ::ChoiceBlock, args::Vector{Layout}, tail::Rope)
    if isempty(args)
        return tail
    end
    rope = fit(fmt, args[1], tail)
    k = 2
    while k <= length(args)
        rope = choice(fmt, rope, fit(fmt, args[k], tail))
        k += 1
    end
    return rope
end

function fit(fmt::Formatter, blk::PenaltyBlock, args::Vector{Layout}, tail::Rope)
    spec = fmt.spec
    spec′ = Spec(spec.line_width, spec.break_cost * blk.break_factor, spec.spill_cost * blk.spill_factor)
    fmt.spec = spec′
    try
        rope′ = fit(fmt, ChoiceBlock(), args, tail)
        if blk.cost == 0
            return rope′
        end
        rope = Rope()
        for (col, span, icept, slope, lt) in rope′
            extend!(rope, col, span, icept + blk.cost, slope, lt)
        end
        return rope
    finally
        fmt.spec = spec
    end
end

function horizontal(fmt::Formatter, rope1::Rope, rope2::Rope)
    if rope2 === fmt.empty
        return rope1
    end
    line_width = fmt.spec.line_width
    spill_cost = fmt.spec.spill_cost
    rope = Rope()
    col = 0
    i1 = 1
    i2 = searchsorted(rope2.cols, rope1.spans[i1]).stop
    while i1 <= length(rope1) && i2 <= length(rope2)
        col1, span1, icept1, slope1, lt1 = rope1[i1]
        col2, span2, icept2, slope2, lt2 = rope2[i2]
        span = span1 + span2
        icept = icept1 + (col - col1) * slope1 +
                icept2 + (col + span1 - col2) * slope2 +
                (col + span1 >= line_width ?
                    -(col + span1 - line_width) * spill_cost : 0)
        slope = slope1 + slope2 +
            (col + span1 >= line_width ? -spill_cost : 0)
        lt = lt1 * lt2
        extend!(rope, col, span, icept, slope, lt)
        if i1 < length(rope1) && i2 < length(rope2)
            col1 = rope1.cols[i1+1]
            col2 = rope2.cols[i2+1] - rope1.spans[i1]
            if col1 <= col2
                col = col1
                i1 += 1
                i2 = searchsorted(rope2.cols, col + rope1.spans[i1]).stop
            else
                col = col2
                i2 += 1
            end
        elseif i1 < length(rope1)
            col = rope1.cols[i1+1]
            i1 += 1
            i2 = searchsorted(rope2.cols, col + rope1.spans[i1]).stop
        elseif i2 < length(rope2)
            col = rope2.cols[i2+1] - rope1.spans[i1]
            i2 += 1
        else
            i1 += 1
        end
    end
    return rope
end

function vertical(fmt::Formatter, rope1::Rope, rope2::Rope)
    break_cost = fmt.spec.break_cost
    rope = Rope()
    col = 0
    i1 = 1
    i2 = 1
    while i1 <= length(rope1) && i2 <= length(rope2)
        col1, span1, icept1, slope1, lt1 = rope1[i1]
        col2, span2, icept2, slope2, lt2 = rope2[i2]
        span = span2
        icept = icept1 + (col - col1) * slope1 +
                icept2 + (col - col2) * slope2 +
                break_cost
        slope = slope1 + slope2
        lt = lt1 / lt2
        extend!(rope, col, span, icept, slope, lt)
        if i1 < length(rope1) && i2 < length(rope2)
            col1 = rope1.cols[i1+1]
            col2 = rope2.cols[i2+1]
            if col1 <= col2
                col = col1
                i1 += 1
            end
            if col1 >= col2
                col = col2
                i2 += 1
            end
        elseif i1 < length(rope1)
            col = rope1.cols[i1+1]
            i1 += 1
        elseif i2 < length(rope2)
            col = rope2.cols[i2+1]
            i2 += 1
        else
            i1 += 1
        end
    end
    return rope
end

function choice(fmt::Formatter, rope1::Rope, rope2::Rope)
    rope = Rope()
    col = 0
    i1 = 1
    i2 = 1
    while i1 <= length(rope1) && i2 <= length(rope2)
        col1, span1, icept1, slope1, lt1 = rope1[i1]
        col2, span2, icept2, slope2, lt2 = rope2[i2]
        cost1 = icept1 + (col - col1) * slope1
        cost2 = icept2 + (col - col2) * slope2
        if cost1 < cost2 || cost1 == cost2 && slope1 <= slope2
            span = span1
            icept = cost1
            slope = slope1
            lt = lt1
        else
            span = span2
            icept = cost2
            slope = slope2
            lt = lt2
        end
        extend!(rope, col, span, icept, slope, lt)
        xcol = col
        if slope1 != slope2
            xcol = ceil(Int, 1.0 * (col1 * slope1 - icept1 - col2 * slope2 + icept2) / (slope1 - slope2))
        end
        if i1 < length(rope1) && i2 < length(rope2)
            col1 = rope1.cols[i1+1]
            col2 = rope2.cols[i2+1]
            if xcol > col && xcol < col1 && xcol < col2
                col = xcol
            else
                if col1 <= col2
                    col = col1
                    i1 += 1
                end
                if col1 >= col2
                    col = col2
                    i2 += 1
                end
            end
        elseif i1 < length(rope1)
            if xcol > col && xcol < rope1.cols[i1+1]
                col = xcol
            else
                col = rope1.cols[i1+1]
                i1 += 1
            end
        elseif i2 < length(rope2)
            if xcol > col && xcol < rope2.cols[i2+1]
                col = xcol
            else
                col = rope2.cols[i2+1]
                i2 += 1
            end
        else
            if xcol > col
                col = xcol
            else
                i1 += 1
            end
        end
    end
    return rope
end


text:
# Helper function for sums excluding the ts's
function lloop(f::Function, ts::Vector{Int64}, Tend::Int64, K::Int64)
    ldenom = 0
    t = ts[2]
    k = 1
    for i=2:Tend
        if i==t
            k+=1
            if k < K
                t=ts[k+1]
            else
                t=Inf
            end
        else
            ldenom += f(i, k)
        end
    end
    return ldenom
end

function ntl_llikelihood(params::Vector{Float64}, ds::Vector{Int64},
        dcounts::Vector{Int64}, ts::Vector{Int64}, K::Int64, Tend::Int64)
    """
    Computes the log-likelihood of the data conditional on arrival times (Tj).

    # Arguments
    - `params::Vector{Float64}`: a length 1 array containing log(1-alpha),
      a transform on the NTL alpha parameter
    - `ds::Vector{Int64}`: a vector of the unique observed cluster sizes
      (multiplicity should be stored in `dcounts`)
    - `dcounts::Vector{Int64}`: corresponding to `ds`, the number of times
      each cluster size occured in data. The orders of `dcounts` and `ds` must
      match
    - `ts::Vector{Int64}`: The vector of observed arrival times T1, ..., TK
       The first element will be 1.
    - `K::Int64`: the total number of clusters, K = length(ts)
    - `Tend::Int64`: the termination time of the data, Tend = ds'*dcounts
    """
    # For optim
    alpha = 1 - exp(params[1])
    if alpha >= 1
        return -Inf
    end
    if alpha < -1e10
        return -lloop((i,k)->log(k), ts, Tend, K)
    end

    lnum = -K*lgamma(1-alpha) + dcounts' * lgamma.(ds - alpha)
    lden = lloop(((i, k)->log(i-1-alpha*k)), ts, Tend, K)

    return lnum - lden
end

function neg_grad_ntl_llikelihood!(storage::Vector{Float64},
        params::Vector{Float64}, ds::Vector{Int64},
        dcounts::Vector{Int64}, ts::Vector{Int64}, K::Int64, Tend::Int64)
    # For optim
    alpha = 1 - exp(params[1])
    trans_correct = -exp(params[1])
    if alpha >= 1
        return Inf
    end
    if alpha < -1e10
        return 0
    end

    glnum = K*digamma(1-alpha) - dcounts' * digamma.(ds - alpha)
    gldenom = lloop((i,k)->-k/(i-1-alpha*k), ts, Tend, K)

    # negative
    storage[1] = -trans_correct * (glnum - gldenom)
end

function geom_llikelihood(g::Vector{Float64}, K::Int64, Tend::Int64)
    # for optim
    g = exp(g[1])/(1+exp(g[1]))
    return (K-1)*log(g) + (Tend-K)*log(1-g)
end

function pyp_arr_llikelihood(params::Vector{Float64}, ds::Vector{Int64},
        dcounts::Vector{Int64}, ts::Vector{Int64}, K::Int64, Tend::Int64)
    # For optim
    tau = exp(params[1])/(1+exp(params[1]))
    theta = params[2]
    if theta <= -tau
        return -Inf
    end
    lnum = lgamma(theta + 1) + sum( log.(theta .+ tau.*collect(1:(K-1))) ) + lloop( (i,k)->log(i-1-k*tau), ts, Tend, K )
    ldenom =  lgamma(theta + Tend)
    return lnum - ldenom
end

function neg_grad_pyp_arr_llikelihood!(storage::Vector{Float64},
        params::Vector{Float64}, ds::Vector{Int64},
        dcounts::Vector{Int64}, ts::Vector{Int64}, K::Int64, Tend::Int64)
    # For optim
    tau = exp(params[1])/(1+exp(params[1]))
    theta = params[2]
    trans_correct_tau = exp(params[1])/(1+exp(params[1]))^2

    if theta <= -tau
        storage[1] = storage[2] = Inf
        return
    end

    # tau negative gradient
    glnum_tau = sum( k/(theta + tau*k ) for k=1:(K-1) ) + lloop((i,k)->-k/(i-1-k*tau), ts, Tend, K)
    gldenom_tau = 0
    storage[1] = -trans_correct_tau*(glnum_tau - gldenom_tau)
    # theta negative gradient
    glnum_theta = digamma(theta + 1) + sum( 1./(theta .+ tau.*collect(1:(K-1))) )
    gldenom_theta = digamma(theta + Tend)
    storage[2] = -(glnum_theta - gldenom_theta)

end

function pyp_llikelihood(params::Vector{Float64}, ds::Vector{Int64},
        dcounts::Vector{Int64}, ts::Vector{Int64}, K::Int64, Tend::Int64)
    # For optim
    tau = exp(params[1])/(1+exp(params[1]))
    theta = params[2]
    if theta <= -tau
        return -Inf
    end
    lnum = -K*lgamma(1-tau) + dcounts' * lgamma.(ds - tau) + sum(log.(theta+tau*collect(1:(K-1))))
    ldenom = sum( log.( (1:Tend-1) + theta ) )

    return lnum - ldenom
end

function neg_grad_pyp_llikelihood!(storage::Vector{Float64}, params::Vector{Float64},
        ds::Vector{Int64}, dcounts::Vector{Int64}, ts::Vector{Int64}, K::Int64,
        Tend::Int64)
    # For optim
    tau = exp(params[1])/(1+exp(params[1]))
    trans_correct = exp(params[1])/(1+exp(params[1]))^2
    theta = params[2]
    if theta <= -tau
        storage[1] = storage[2] = Inf
        return
    end

    storage[1] = -trans_correct * ( K*digamma(1-tau) - dcounts' * digamma.(ds - tau) + sum( k/(theta+tau*k) for k=1:(K-1) ) )
    storage[2] = -( sum( 1./(theta .+ tau*collect(1:(K-1))) ) - sum( 1./(theta .+ (1:(Tend-1))) ) )
end


text:
"""
    Magnetization(x, y, z)
    Magnetization{T}()
    Magnetization()

Create a mutable `Magnetization` object representing a 3D magnetization vector.

# Properties
- `x::Real`: x component of magnetization vector
- `y::Real`: y component of magnetization vector
- `z::Real`: z component of magnetization vector

# Examples
```jldoctest
julia> M = Magnetization()
Magnetization vector with eltype Float64:
 Mx = 0.0
 My = 0.0
 Mz = 0.0

julia> M.x = 1; M.y = 2; M.z = 3; M
Magnetization vector with eltype Float64:
 Mx = 1.0
 My = 2.0
 Mz = 3.0
```
"""
mutable struct Magnetization{T<:Real}
    x::T
    y::T
    z::T
end

Magnetization(x, y, z) = Magnetization(promote(x, y, z)...)
Magnetization{T}() where {T} = Magnetization(zero(T), zero(T), zero(T))
Magnetization() = Magnetization(0.0, 0.0, 0.0)

Base.show(io::IO, M::Magnetization) = print(io, "Magnetization(", M.x, ", ", M.y, ", ", M.z, ")")
Base.show(io::IO, ::MIME"text/plain", M::Magnetization{T}) where {T} =
    print(io, "Magnetization vector with eltype $T:\n Mx = ", M.x, "\n My = ", M.y, "\n Mz = ", M.z)

Base.copy(M::Magnetization) = Magnetization(M.x, M.y, M.z)

function Base.copyto!(dst::Magnetization, src::Magnetization)

    dst.x = src.x
    dst.y = src.y
    dst.z = src.z
    return nothing

end

Base.eltype(::Magnetization{T}) where {T} = T
Base.convert(::Type{Magnetization{T}}, M::Magnetization) where {T} = Magnetization(T(M.x), T(M.y), T(M.z))

Base.Vector(M::Magnetization) = [M.x, M.y, M.z]

function Base.copyto!(dst::Magnetization, src::AbstractVector)

    dst.x = src[1]
    dst.y = src[2]
    dst.z = src[3]
    return nothing

end

function Base.copyto!(dst::AbstractVector, src::Magnetization)

    dst[1] = src.x
    dst[2] = src.y
    dst[3] = src.z
    return nothing

end

Base.:(==)(M1::Magnetization, M2::Magnetization) = M1.x == M2.x && M1.y == M2.y && M1.z == M2.z
Base.isapprox(M1::Magnetization, M2::Magnetization; kwargs...) =
    isapprox(Vector(M1), Vector(M2); kwargs...)

"""
    signal(spin)
    signal(M)

Return the signal detected from the given spin or magnetization vector.

# Examples
```jldoctest
julia> signal(Spin(Magnetization(1, 2, 3), 1, 1000, 100, 0))
1.0 + 2.0im

julia> signal(MagnetizationMC((1, 2, 3), (1, 1, 1)))
2 + 3im
```
"""
signal(M::Magnetization) = complex(M.x, M.y)

"""
    MagnetizationMC((x1, y1, z1), ...)
    MagnetizationMC{T}(N)
    MagnetizationMC(N)

Create a `MagnetizationMC` object representing `N` 3D magnetization vectors in
the same physical location.

One can access the ith component magnetization vector by indexing the
`MagnetizationMC` object. Furthermore, iterating the `MagnetizationMC` object
iterates through each of the component magnetization vectors.

# Properties
- `M::NTuple{N,Magnetization{<:Real}}`: List of component magnetization vectors

# Examples
```jldoctest
julia> M = MagnetizationMC((1, 2, 3), (4, 5, 6))
2-compartment Magnetization vector with eltype Int64:
 Compartment 1:
  Mx = 1
  My = 2
  Mz = 3
 Compartment 2:
  Mx = 4
  My = 5
  Mz = 6

julia> M[2]
Magnetization vector with eltype Int64:
 Mx = 4
 My = 5
 Mz = 6

julia> foreach(m -> (m.x = 0; m.y = 1; m.z = 2), M)

julia> M
2-compartment Magnetization vector with eltype Int64:
 Compartment 1:
  Mx = 0
  My = 1
  Mz = 2
 Compartment 2:
  Mx = 0
  My = 1
  Mz = 2
```
"""
struct MagnetizationMC{T<:Real,N}
    M::NTuple{N,Magnetization{T}}

    function MagnetizationMC(M::Magnetization...)

        N = length(M)
        N > 1 || error(sprint(print, "MagnetizationMC expects 2 or more compartments, got ", N))
        T = promote_type(eltype.(M)...)
        args = convert.(Magnetization{T}, M)
        new{T,N}(args)

    end
end

MagnetizationMC(M::NTuple{3,Real}...) = MagnetizationMC((Magnetization(Mi...) for Mi in M)...)
MagnetizationMC{T}(N) where {T} = MagnetizationMC(ntuple(i -> Magnetization{T}(), N)...)
MagnetizationMC(N) = MagnetizationMC(ntuple(i -> Magnetization(), N)...)

function Base.show(io::IO, M::MagnetizationMC{T,N}) where {T,N}

    print(io, "MagnetizationMC((", M[1].x, ", ", M[1].y, ", ", M[1].z, "), (")
    for i = 2:N-1
        print(io, M[i].x, ", ", M[i].y, ", ", M[i].z, "), (")
    end
    print(io, M[N].x, ", ", M[N].y, ", ", M[N].z, "))")

end

function Base.show(io::IO, ::MIME"text/plain", M::MagnetizationMC{T,N}) where {T,N}

    print(io, "$N-compartment Magnetization vector with eltype $T:")
    for i = 1:N
        print(io, "\n Compartment ", i, ":")
        print(io, "\n  Mx = ", M[i].x)
        print(io, "\n  My = ", M[i].y)
        print(io, "\n  Mz = ", M[i].z)
    end

end

Base.copy(M::MagnetizationMC{T,N}) where {T,N} = MagnetizationMC(ntuple(i -> copy(M[i]), N)...)
Base.copyto!(dst::MagnetizationMC{T1,N}, src::MagnetizationMC{T2,N}) where {T1,T2,N} =
    foreach(i -> copyto!(dst[i], src[i]), 1:N)
Base.copyto!(dst::MagnetizationMC{T,N}, src::AbstractVector) where {T,N} =
    foreach(i -> copyto!(dst[i], view(src, 3i-2:3i)), 1:N)
Base.copyto!(dst::AbstractVector, src::MagnetizationMC{T,N}) where {T,N} =
    foreach(i -> copyto!(view(dst, 3i-2:3i), src[i]), 1:N)

Base.eltype(::MagnetizationMC{T,N}) where {T,N} = T
Base.getindex(M::MagnetizationMC, i) = M.M[i]
Base.iterate(M::MagnetizationMC{T,N}, i = 1) where {T,N} =  i > N ? nothing : (M.M[i], i + 1)
Base.convert(::Type{MagnetizationMC{T1,N}}, M::MagnetizationMC{T2,N}) where {T1,T2,N} = MagnetizationMC((convert(Magnetization{T1}, Mi) for Mi in M)...)
# This next definition of convert prevents StackOverflowErrors
Base.convert(::Type{MagnetizationMC{T,N}}, M::MagnetizationMC{T,N}) where {T,N} = M

function Base.Vector(M::MagnetizationMC{T,N}) where {T,N}

    v = Vector{T}(undef, 3N)
    for i = 1:N
        v[3i-2] = M[i].x
        v[3i-1] = M[i].y
        v[3i]   = M[i].z
    end
    return v

end

Base.:(==)(M1::MagnetizationMC{T1,N}, M2::MagnetizationMC{T2,N}) where {T1,T2,N} = all(M1[i] == M2[i] for i = 1:N)
Base.isapprox(M1::MagnetizationMC{T1,N}, M2::MagnetizationMC{T2,N}; kwargs...) where {T1,T2,N} = all(isapprox(M1[i], M2[i]; kwargs...) for i = 1:N)

signal(M::MagnetizationMC{T,N}) where {T,N} = sum(signal(M) for M in M)


text:
# Load data from text files

function load_data(folder, file_q_obs = "Q_obs.txt", file_tair = "Tair.txt",
                   file_prec = "Prec.txt", file_frac = "Frac.txt")

  # Read air temperature data

  str   = readline("$folder/$file_tair")
  nsep  = length(matchall(r";", str))
  tmp   = CSV.read("$folder/$file_tair", delim = ";", header = false,
                   dateformat="yyyy-mm-dd HH:MM", nullable = false, types = vcat(DateTime, repmat([Float64], nsep)))
  tair  = Array(tmp[:, 2:end])
  tair  = transpose(tair)

  # Read precipitation data

  str   = readline("$folder/$file_tair")
  nsep  = length(matchall(r";", str))
  tmp   = CSV.read("$folder/$file_prec", delim = ";", header = false,
                  dateformat="yyyy-mm-dd HH:MM", nullable = false, types = vcat(DateTime, repmat([Float64], nsep)))
  prec  = Array(tmp[:, 2:end])
  prec  = transpose(prec)

  # Read runoff data

  tmp   = CSV.read("$folder/$file_q_obs", delim = ";", header = false,
                   dateformat="yyyy-mm-dd HH:MM", nullable = false, types = [DateTime, Float64])
  q_obs = Array(tmp[:, 2])

  q_obs[q_obs .== -999.0] = NaN

  # Read elevation band data

  frac = readdlm("$folder/$file_frac")
  frac = squeeze(frac,2)

  # Get time data

  date = Array(tmp[:, 1])

  # Return data

  return date, tair, prec, q_obs, frac

end

# Load operational data from text files

function load_operational(folder, file_q_obs = "Q_obs.txt", file_tair = "Tair.txt",
                   file_prec = "Prec.txt", file_metadata = "metadata.txt")

  # Read air temperature data

  str   = readline("$folder/$file_tair")
  nsep  = length(matchall(r";", str))
  tmp   = CSV.read("$folder/$file_tair", delim = ";", header = false,
                   dateformat="yyyy-mm-dd HH:MM", nullable = false, types = vcat(DateTime, repmat([Float64], nsep)))
  tair  = Array(tmp[:, 2:end])
  tair  = transpose(tair)

  # Read precipitation data

  str   = readline("$folder/$file_tair")
  nsep  = length(matchall(r";", str))
  tmp   = CSV.read("$folder/$file_prec", delim = ";", header = false,
                  dateformat="yyyy-mm-dd HH:MM", nullable = false, types = vcat(DateTime, repmat([Float64], nsep)))
  prec  = Array(tmp[:, 2:end])
  prec  = transpose(prec)

  # Read runoff data

  tmp   = CSV.read("$folder/$file_q_obs", delim = ";", header = false,
                   dateformat="yyyy-mm-dd HH:MM", nullable = false, types = [DateTime, Float64])
  q_obs = Array(tmp[:, 2])

  q_obs[q_obs .== -999.0] = NaN

  # Read elevation band data

  metadata = readtable("$folder/$file_metadata", separator = ';')
  frac = convert(Array{Float64,1}, metadata[:area] / sum(metadata[:area]))

  # Get time data

  date = Array(tmp[:, 1])

  # Return data

  return date, tair, prec, q_obs, frac

end

# Crop data from start to stop date

function crop_data(date, tair, prec, q_obs, date_start, date_stop)

  # Find indicies

  istart = find(date .== date_start)
  istop = find(date .== date_stop)

  # Test if ranges are valid

  if isempty(istart) | isempty(istop)
    error("Cropping data outside range")
  end

  # Crop data

  date  = date[istart[1]:istop[1]]
  tair  = tair[:, istart[1]:istop[1]]
  prec  = prec[:, istart[1]:istop[1]]
  q_obs = q_obs[istart[1]:istop[1]]

  return date, tair, prec, q_obs

end

# Crop data before start date

function crop_data(date, tair, prec, q_obs, date_start)

  # Find indicies

  istart = find(date .== date_start)

  # Test if ranges are valid

  if isempty(istart)
    error("Cropping data outside range")
  end

  # Crop data

  date  = date[istart[1]:end]
  tair  = tair[:, istart[1]:end]
  prec  = prec[:, istart[1]:end]
  q_obs = q_obs[istart[1]:end]

  return date, tair, prec, q_obs

end

# Assign input data to snow model

function get_input(st_snow::Snow, prec, tair, itime)

  for izone in eachindex(st_snow.prec)

    st_snow.prec[izone] = prec[izone, itime]
    st_snow.tair[izone] = tair[izone, itime]

  end

end

# Assign input data to hydrological response model

function get_input(st_hydro::Hydro, prec, epot, itime)

  st_hydro.prec = prec[itime]
  st_hydro.epot = epot[itime]

end

function get_input(st_snow::Snow, st_hydro::Hydro, epot, itime)

  st_hydro.prec = sum(st_snow.frac .* st_snow.q_sim)
  st_hydro.epot = epot[itime]

end

max_stars_repo_name: jmgnve/Vann

max_forks_repo_name: jmgnve/Vann


text:
struct CoupledControlParams{T}
    V_min::T
    V_max::T

    k_V::T    # used only for initial seeding of linearization nodes
    k_s::T

    δ̇_max::T

    Q_Δs::T
    Q_Δψ::T
    Q_e::T
    W_β::T
    W_r::T
    W_HJI::T
    N_HJI::Int

    R_δ::T
    R_Δδ::T
    R_Fx::T
    R_ΔFx::T
end
function CoupledControlParams(;V_min=1.0,
                               V_max=15.0,
                               k_V=10/4/100,
                               k_s=10/4/10000,
                               δ̇_max=0.344,
                               Q_Δs=1.0,
                               Q_Δψ=1.0,
                               Q_e=1.0,
                               W_β=50/(10*π/180),
                               W_r=50.0,
                               W_HJI=500.0,
                               N_HJI=3,
                               R_δ=0.0,
                               R_Δδ=0.1,
                               R_Fx=0.0,
                               R_ΔFx=0.5)
    CoupledControlParams(V_min, V_max, k_V, k_s, δ̇_max, Q_Δs, Q_Δψ, Q_e, W_β, W_r, W_HJI, N_HJI, R_δ, R_Δδ, R_Fx, R_ΔFx)
end

function CoupledTrajectoryTrackingMPC(vehicle::Dict{Symbol,T}, trajectory::TrajectoryTube{T}; control_params=CoupledControlParams(),
                                      N_short=10, N_long=20, dt_short=T(0.01), dt_long=T(0.2), use_correction_step=true) where {T}
    dynamics = VehicleModel(vehicle)
    current_state = zeros(BicycleState{T})
    current_control = zeros(BicycleControl{T})
    time_steps = MPCTimeSteps(N_short, N_long, dt_short, dt_long, use_correction_step)

    N = 1 + N_short + N_long
    qs = rand(TrackingBicycleState{T}, N)    # not zeros so that construct_lateral_tracking_QP below works
    us = rand(BicycleControl2{T}, N)
    ps = zeros(TrackingBicycleParams{T}, N)
    tracking_dynamics = VehicleModel(vehicle, TrackingBicycleModel(vehicle))
    model, variables, parameters = construct_coupled_tracking_QP(tracking_dynamics, control_params, time_steps, qs, us, ps)
    TrajectoryTrackingMPC(vehicle, trajectory, dynamics, control_params,
                          current_state, current_control, 0, NaN,
                          time_steps,
                          qs, us, ps,
                          tracking_dynamics, model, variables, parameters)
end

function compute_linearization_nodes!(mpc::TrajectoryTrackingMPC{T},
                                      qs::Vector{TrackingBicycleState{T}},
                                      us::Vector{BicycleControl2{T}},
                                      ps::Vector{TrackingBicycleParams{T}}) where {T}
    traj = mpc.trajectory
    X  = mpc.tracking_dynamics
    U  = mpc.control_params
    q0 = mpc.current_state
    u0 = mpc.current_control
    TS = mpc.time_steps
    N_short, N_long = TS.N_short, TS.N_long
    ts, dt, prev_ts = TS.ts, TS.dt, TS.prev_ts

    s0, e0, t0 = path_coordinates(traj, q0)
    tj = traj[s0]
    Δs = s0 - traj(ts[1]).s
    Δψ = adiff(q0.ψ, tj.ψ)
    q  = TrackingBicycleState(Δs, q0.Ux, q0.Uy, q0.r, Δψ, e0)
    u  = BicycleControl2(u0)
    p  = TrackingBicycleParams(tj.V, tj.κ, T(0), T(0))
    if mpc.solved
        qs[1] = q
        us[1] = u
        ps[1] = p
        update_interpolations!(mpc.variables, mpc.model, prev_ts)
        for i in 2:N_short+N_long+1
            t = ts[i]
            if t < prev_ts[end]
                q = mpc.variables.q_interp(t)
                u = mpc.variables.u_interp(t) .* mpc.variables.u_normalization
            else
                q = mpc.variables.q_interp(prev_ts[end])
                u = mpc.variables.u_interp(prev_ts[end]) .* mpc.variables.u_normalization
            end
            s = traj(t).s + q.Δs
            tj = traj[s]
            p = TrackingBicycleParams(tj.V, tj.κ, T(0), T(0))
            qs[i] = q
            us[i] = u
            ps[i] = p
        end
    else    # steady-state estimates more-or-less the same as in decoupled MPC
        s  = s0
        sΔψ, cΔψ = sincos(Δψ)
        V  = q0.Ux*cΔψ - q0.Uy*sΔψ    # hypot(q0.Ux, q0.Uy)
        β0 = atan(q0.Uy, q0.Ux)
        r0 = q0.r
        δ0 = u0.δ
        Fyf0, Fyr0 = lateral_tire_forces(X.bicycle_model, q0, u0)
        for i in 1:N_short+N_long+1
            τ = (i == N_short+N_long+1 ? dt[i-1] : dt[i])
            tj = traj[s]
            Δs = s - traj(ts[i]).s
            A_des = tj.A + U.k_V*(tj.V - V)/τ + (isnan(mpc.time_offset) ? 0 : -U.k_s*Δs/τ/τ)
            A_des = min(max(A_des, (U.V_min - V)/τ), (U.V_max - V)/τ)
            if i == 1
                q̇ = mpc.dynamics(q0, BicycleControl2(u0), LocalRoadGeometry(tj))
                A = (q̇[4] - q0.r*q0.Uy)*cΔψ - (q̇[5] + q0.r*q0.Ux)*sΔψ
                # A = (q̇[4] - q0.r*q0.Uy)*cos(β0) + (q̇[5] + q0.r*q0.Ux)*sin(β0)
            elseif i <= N_short+1
                est = steady_state_estimates(X, V, A_des, tj.κ, num_iters=1, r=r0, β0=β0, δ0=δ0, Fyf0=Fyf0)
                q = TrackingBicycleState(Δs, q0.Ux, q0.Uy, q0.r, adiff(q0.ψ, tj.ψ), e0)
                u = BicycleControl2(est.δ, est.Fxf + est.Fxr)
                p = TrackingBicycleParams(tj.V, tj.κ, T(0), T(0))
                A = est.A
            else
                est = steady_state_estimates(X, V, A_des, tj.κ)
                q = TrackingBicycleState(Δs, est.Ux, est.Uy, est.r, -est.β, T(0))
                u = BicycleControl2(est.δ, est.Fxf + est.Fxr)
                p = TrackingBicycleParams(tj.V, tj.κ, T(0), T(0))
                A = est.A
            end
            qs[i] = q
            us[i] = u
            ps[i] = p
            i == N_short+N_long+1 && break
            V = V + A*τ
            s = s + V*τ + A*τ*τ/2
        end
    end
end

struct TrackingQPParams{T}
    Q_Δs  ::Parameter{Diagonal{T,Vector{T}},typeof(identity),true}
    Q_Δψ  ::Parameter{Diagonal{T,Vector{T}},typeof(identity),true}
    Q_e   ::Parameter{Diagonal{T,Vector{T}},typeof(identity),true}
    R_δ   ::Parameter{Diagonal{T,Vector{T}},typeof(identity),true}
    R_Δδ  ::Parameter{Diagonal{T,Vector{T}},typeof(identity),true}
    R_Fx  ::Parameter{Diagonal{T,Vector{T}},typeof(identity),true}
    R_ΔFx ::Parameter{Diagonal{T,Vector{T}},typeof(identity),true}
    W_β   ::Parameter{Vector{T},typeof(identity),true}
    W_r   ::Parameter{Vector{T},typeof(identity),true}
    W_HJI ::Parameter{Vector{T},typeof(identity),true}
    q_curr::Parameter{Vector{T},typeof(identity),true}
    u_curr::Parameter{Vector{T},typeof(identity),true}
    M_HJI ::Parameter{Matrix{T},typeof(identity),true}
    b_HJI ::Parameter{Matrix{T},typeof(identity),true}
    A     ::Vector{Parameter{Matrix{T},typeof(identity),true}}
    B     ::Vector{Parameter{Matrix{T},typeof(identity),true}}
    B0    ::Vector{Parameter{Matrix{T},typeof(identity),true}}
    Bf    ::Vector{Parameter{Matrix{T},typeof(identity),true}}
    c     ::Vector{Parameter{Vector{T},typeof(identity),true}}
    H     ::Vector{Parameter{Matrix{T},typeof(identity),true}}
    G     ::Vector{Parameter{Vector{T},typeof(identity),true}}
    δ_min ::Vector{Parameter{Vector{T},typeof(identity),true}}
    δ_max ::Vector{Parameter{Vector{T},typeof(identity),true}}
    Fx_max::Vector{Parameter{Vector{T},typeof(identity),true}}
    Δδ_min::Vector{Parameter{Vector{T},typeof(identity),true}}
    Δδ_max::Vector{Parameter{Vector{T},typeof(identity),true}}
end

struct TrackingQPVariables{T}
    q::Matrix{Variable}
    u::Matrix{Variable}
    σ::Matrix{Variable}
    σ_HJI::Vector{Variable}
    u_normalization::SVector{2,T}
    q_interp::GriddedInterpolation{TrackingBicycleState{T},1,TrackingBicycleState{T},Gridded{Linear},Tuple{Vector{T}}}
    u_interp::GriddedInterpolation{BicycleControl2{T},1,BicycleControl2{T},Gridded{Linear},Tuple{Vector{T}}}
end
function TrackingQPVariables(q::Matrix{Variable}, u::Matrix{Variable}, σ::Matrix{Variable}, σ_HJI::Vector{Variable}, u_normalization::SVector{2,T}) where {T}
    N = size(q, 2)
    q_interp = interpolate(T, TrackingBicycleState{T}, (T.(1:N),), zeros(TrackingBicycleState{T}, N), Gridded(Linear()))
    u_interp = interpolate(T, BicycleControl2{T}, (T.(1:N),), zeros(BicycleControl2{T}, N), Gridded(Linear()))
    TrackingQPVariables(q, u, σ, σ_HJI, u_normalization, q_interp, u_interp)
end

function update_interpolations!(variables::TrackingQPVariables, model::Model{T}, prev_ts) where {T}
    q_interp, u_interp = variables.q_interp, variables.u_interp
    q_interp.knots[1] .= prev_ts
    u_interp.knots[1] .= prev_ts
    reshape(reinterpret(T, q_interp.coefs), (length(eltype(q_interp)), length(q_interp))) .= value.(Ref(model), variables.q)
    reshape(reinterpret(T, u_interp.coefs), (length(eltype(u_interp)), length(u_interp))) .= value.(Ref(model), variables.u)
end

function construct_coupled_tracking_QP(dynamics::VehicleModel{T}, control_params, time_steps, qs, us, ps) where {T}
    N_short, N_long, dt = time_steps.N_short, time_steps.N_long, time_steps.dt
    u_normalization = SVector{2,T}(dynamics.control_limits.δ_max, max(-dynamics.control_limits.Fx_min, dynamics.control_limits.Fx_max))

    optimizer = OSQPOptimizer()
    MOI.set!(optimizer, OSQPSettings.Verbose(), false)
    MOI.set!(optimizer, OSQPSettings.WarmStart(), true)
    m = Model{T}(optimizer)

    Q_Δs   = Parameter(Diagonal(control_params.Q_Δs .* dt), m)
    Q_Δψ   = Parameter(Diagonal(control_params.Q_Δψ .* dt), m)
    Q_e    = Parameter(Diagonal(control_params.Q_e  .* dt), m)
    R_δ    = Parameter(Diagonal(control_params.R_δ  .* dt), m)
    R_Δδ   = Parameter(Diagonal(control_params.R_Δδ ./ dt), m)
    R_Fx   = Parameter(Diagonal(control_params.R_Fx  .* dt), m)
    R_ΔFx  = Parameter(Diagonal(control_params.R_ΔFx ./ dt), m)
    W_β    = Parameter(control_params.W_β .* dt, m)
    W_r    = Parameter(control_params.W_r .* dt, m)
    W_HJI  = Parameter(control_params.W_HJI .* ones(N_short), m)
    q_curr = Parameter(Array(qs[1]), m)
    u_curr = Parameter(Array(us[1] ./ u_normalization), m)
    M_HJI  = Parameter(zeros(T, 1, 2), m)
    b_HJI  = Parameter(ones(T, 1, N_short), m)
    A      = Parameter{Matrix{T},typeof(identity),true}[]
    B      = Parameter{Matrix{T},typeof(identity),true}[]
    B0     = Parameter{Matrix{T},typeof(identity),true}[]
    Bf     = Parameter{Matrix{T},typeof(identity),true}[]
    c      = Parameter{Vector{T},typeof(identity),true}[]
    H      = Parameter{Matrix{T},typeof(identity),true}[]
    G      = Parameter{Vector{T},typeof(identity),true}[]
    δ_min  = Parameter{Vector{T},typeof(identity),true}[]
    δ_max  = Parameter{Vector{T},typeof(identity),true}[]
    Fx_max = Parameter{Vector{T},typeof(identity),true}[]
    Δδ_min = Parameter{Vector{T},typeof(identity),true}[]
    Δδ_max = Parameter{Vector{T},typeof(identity),true}[]

    q     = [Variable(m) for i in 1:6, t in 1:N_short+N_long+1]    # (Δs, Ux, Uy, r, Δψ, e)
    u     = [Variable(m) for i in 1:2, t in 1:N_short+N_long+1]    # (δ, Fx)
    σ     = [Variable(m) for i in 1:2, t in 1:N_short+N_long]
    σ_HJI = [Variable(m) for t in 1:N_short]
    Δδ    = [Variable(m) for t in 1:N_short+N_long]
    ΔFx   = [Variable(m) for t in 1:N_short+N_long]
    δ     = u[1,:]
    Fx    = u[2,:]
    Ux    = q[2,:]
    @constraint(m, vec(σ) >= fill(T(0), 2*(N_short+N_long)))
    @constraint(m, σ_HJI  >= fill(T(0), N_short))
    @constraint(m, diff(δ) == Δδ)
    @constraint(m, diff(Fx) == ΔFx)
    @constraint(m, Ux >= fill(control_params.V_min, N_short+N_long+1))    # may need slacks
    @constraint(m, Ux <= fill(control_params.V_max, N_short+N_long+1))
    @constraint(m, Fx >= fill(dynamics.control_limits.Fx_min / u_normalization[2], N_short+N_long+1))

    @constraint(m, q[:,1] == q_curr)
    @constraint(m, u[:,1] == u_curr)
    for t in 1:N_short
        ZOHt = linearize(dynamics, qs[t], StepControl(dt[t], [us[t]; ps[t]]), keep_control_dims=SVector(1,2))
        At = push!(A, Parameter(Array(ZOHt.A), m))[end]
        Bt = push!(B, Parameter(Array(ZOHt.B .* u_normalization'), m))[end]
        ct = push!(c, Parameter(Array(ZOHt.c), m))[end]
        @constraint(m, At*q[:,t] + Bt*u[:,t] + ct == q[:,t+1])
    end
    @constraint(m, vec(M_HJI*u[:,1:N_short] + b_HJI) >= -σ_HJI)

    for t in N_short+1:N_short+N_long
        FOHt = linearize(dynamics, qs[t], RampControl(dt[t], [us[t]; ps[t]], [us[t+1]; ps[t+1]]), keep_control_dims=SVector(1,2))
        At  = push!(A,  Parameter(Array(FOHt.A),  m))[end]
        B0t = push!(B0, Parameter(Array(FOHt.B0 .* u_normalization'), m))[end]
        Bft = push!(Bf, Parameter(Array(FOHt.Bf .* u_normalization'), m))[end]
        ct  = push!(c,  Parameter(Array(FOHt.c),  m))[end]
        @constraint(m, At*q[:,t] + B0t*u[:,t] + Bft*u[:,t+1] + ct == q[:,t+1])
    end

    δ_max_, δ̇_max   = dynamics.control_limits.δ_max, control_params.δ̇_max
    Fx_max_, Px_max = dynamics.control_limits.Fx_max, dynamics.control_limits.Px_max
    for t in 1:N_short+N_long
        Uxt = qs[t+1].Ux
        Fxft, Fxrt = longitudinal_tire_forces(dynamics.longitudinal_params, us[t+1].Fx)
        δ_mint, δ_maxt, Ht, Gt = stable_limits(dynamics.bicycle_model, Uxt, Fxft, Fxrt)
        Ht = push!(H, Parameter(Array(Ht), m))[end]
        Gt = push!(G, Parameter(Array(Gt), m))[end]
        δ_mint  = push!(δ_min,  Parameter([max(δ_mint, -δ_max_) / u_normalization[1]], m))[end]
        δ_maxt  = push!(δ_max,  Parameter([min(δ_maxt, δ_max_)  / u_normalization[1]], m))[end]
        Fx_maxt = push!(Fx_max, Parameter([min(Px_max/Uxt, Fx_max_) / u_normalization[2]], m))[end]
        Uy_r = q[3:4,t+1]
        σt   = [σ[1,t],σ[1,t],σ[2,t],σ[2,t]]
        Δδt  = [Δδ[t]]
        Δδ_mint = push!(Δδ_min, Parameter([-δ̇_max*dt[t] / u_normalization[1]], m))[end]
        Δδ_maxt = push!(Δδ_max, Parameter([ δ̇_max*dt[t] / u_normalization[1]], m))[end]
        @constraint(m, [δ[t+1]] <= δ_maxt)
        @constraint(m, [δ[t+1]] >= δ_mint)
        @constraint(m, [Fx[t+1]] <= Fx_maxt)
        @constraint(m, Ht*Uy_r - Gt <= σt)
        @constraint(m, Δδt <= Δδ_maxt)
        @constraint(m, Δδt >= Δδ_mint)
    end

    Δs   = q[1,2:end]
    Δψ   = q[5,2:end]
    e    = q[6,2:end]
    δ⁺   = δ[2:end]
    Fx⁺  = Fx[2:end]
    σ1   = σ[1,:]
    σ2   = σ[2,:]
    obj  = @expression transpose(Δs)*Q_Δs*Δs +
                       transpose(Δψ)*Q_Δψ*Δψ +
                       transpose(e)*Q_e*e +
                       transpose(δ⁺)*R_δ*δ⁺ +
                       transpose(Δδ)*R_Δδ*Δδ +
                       transpose(Fx⁺)*R_Fx*Fx⁺ +
                       transpose(ΔFx)*R_ΔFx*ΔFx +
                       W_β⋅σ1 + W_r⋅σ2 + W_HJI⋅σ_HJI
    @objective(m, Minimize, obj)
    m, TrackingQPVariables(q, u, σ, σ_HJI, u_normalization), TrackingQPParams(Q_Δs, Q_Δψ, Q_e, R_δ, R_Δδ, R_Fx, R_ΔFx, W_β, W_r, W_HJI,
                                                                              q_curr, u_curr, M_HJI, b_HJI, A, B, B0, Bf, c, H, G,
                                                                              δ_min, δ_max, Fx_max, Δδ_min, Δδ_max)
end

function update_QP!(mpc::TrajectoryTrackingMPC, QPP::TrackingQPParams)
    dynamics = mpc.tracking_dynamics
    control_params  = mpc.control_params
    time_steps = mpc.time_steps
    N_short, N_long, dt = time_steps.N_short, time_steps.N_long, time_steps.dt
    qs, us, ps = mpc.qs, mpc.us, mpc.ps
    u_normalization = mpc.variables.u_normalization

    QPP.Q_Δs().diag  .= control_params.Q_Δs .* dt
    QPP.Q_Δψ().diag  .= control_params.Q_Δψ .* dt
    QPP.Q_e().diag   .= control_params.Q_e  .* dt
    QPP.R_δ().diag   .= control_params.R_δ  .* dt
    QPP.R_Δδ().diag  .= control_params.R_Δδ ./ dt
    QPP.R_Fx().diag  .= control_params.R_Fx  .* dt
    QPP.R_ΔFx().diag .= control_params.R_ΔFx ./ dt
    QPP.W_β()    .= control_params.W_β .* dt
    QPP.W_r()    .= control_params.W_r .* dt
    QPP.q_curr() .= qs[1]
    QPP.u_curr() .= us[1] ./ u_normalization

    for t in 1:N_short
        ZOHt = linearize(dynamics, qs[t], StepControl(dt[t], [us[t]; ps[t]]), keep_control_dims=SVector(1,2))
        QPP.A[t]() .= ZOHt.A
        QPP.B[t]() .= (ZOHt.B .* u_normalization')
        QPP.c[t]() .= ZOHt.c
    end
    relative_state = HJIRelativeState(mpc.current_state, mpc.other_car_state)
    M, b = compute_reachability_constraint(mpc.dynamics, mpc.HJI_cache, relative_state, mpc.HJI_ϵ, BicycleControl2(mpc.current_control))
    # QPP.W_HJI() .= control_params.W_HJI .* ones(N_short)
    QPP.W_HJI() .= control_params.W_HJI .* [ones(control_params.N_HJI); zeros(N_short - control_params.N_HJI)]
    QPP.M_HJI() .= (M .* u_normalization)'
    QPP.b_HJI() .= b
    for t in N_short+1:N_short+N_long
        FOHt = linearize(dynamics, qs[t], RampControl(dt[t], [us[t]; ps[t]], [us[t+1]; ps[t+1]]), keep_control_dims=SVector(1,2))
        QPP.A[t]() .= FOHt.A
        QPP.B0[t-N_short]() .= (FOHt.B0 .* u_normalization')
        QPP.Bf[t-N_short]() .= (FOHt.Bf .* u_normalization')
        QPP.c[t]() .= FOHt.c
    end
    δ_max, δ̇_max   = dynamics.control_limits.δ_max, control_params.δ̇_max
    Fx_max, Px_max = dynamics.control_limits.Fx_max, dynamics.control_limits.Px_max
    for t in 1:N_short+N_long
        Uxt = qs[t+1].Ux
        Fxft, Fxrt = longitudinal_tire_forces(dynamics.longitudinal_params, us[t+1].Fx)
        δ_mint, δ_maxt, Ht, Gt = stable_limits(dynamics.bicycle_model, Uxt, Fxft, Fxrt)
        QPP.H[t]() .= Ht
        QPP.G[t]() .= Gt
        QPP.δ_min[t]()  .= max(δ_mint, -δ_max) / u_normalization[1]
        QPP.δ_max[t]()  .= min(δ_maxt,  δ_max) / u_normalization[1]
        QPP.Fx_max[t]() .= min(Px_max/Uxt, Fx_max) / u_normalization[2]
        QPP.Δδ_min[t]() .= -δ̇_max*dt[t] / u_normalization[1]
        QPP.Δδ_max[t]() .=  δ̇_max*dt[t] / u_normalization[1]
    end
end

function get_next_control(mpc::TrajectoryTrackingMPC, variables::TrackingQPVariables)
    u_normalization = variables.u_normalization
    δ, Fx = value.(Ref(mpc.model), SVector(variables.u[1,2], variables.u[2,2]))
    BicycleControl(mpc.dynamics.longitudinal_params, BicycleControl2(δ, Fx) .* u_normalization)
end

max_stars_repo_name: Michelle-NYX/CollisionAvoidance

max_issues_repo_name: Michelle-NYX/CollisionAvoidance

max_forks_repo_name: Michelle-NYX/CollisionAvoidance


text:
isdefined(Base, :__precompile__) && __precompile__()

module VehicleModels

using Parameters
using JuMP
using Interpolations
using OrdinaryDiffEq
using DiffEqBase
using NLOptControl  # to use newConstraint!() and interpolateLagrange!() for checkCrash()

# funcitons in the VehicleModels.jl package
include("parameters.jl")
include("Three_DOF/Three_DOF.jl")
include("KinematicBicycle/KinematicBicycle.jl")
include("utils.jl")

export
  #########
  # Objects
  #########
  Vpara,

  ###########
  # Functions
  ###########
  checkCrash,
  FNames,

  # Three DOF
  ThreeDOFv1,
  ThreeDOFv2,
  ThreeDOFv2_expr,
  ThreeDOFv3,
  ThreeDOFv3_expr,

  # KinematicBicycle
  KinematicBicycle,
  KinematicBicycle_expr,
  KinematicBicycle2,
  KinematicBicycle_expr2,

  ###############################
  # Macros and support functions
  ###############################
  # Three DOF
  @F_YF,
  @F_YR,
  @FZ_RL,
  @FZ_RR,
  @FZ_FL,
  @FZ_FR,
  @Ax_min,
  @Ax_max,

  # parameters
  @unpack_Vpara,
  @pack_Vpara,

  # Parameters.jl
  @unpack,
  @pack

end # module


text:
using ArgParse

s=ArgParseSettings()
@add_arg_table s begin
    "--part2"
        help="Do part 2"
        action=:store_true
end
parsed_args = parse_args(ARGS, s)

function calculate_gamma(numbers, bitlen)
    counts = zeros(bitlen)
    for n in numbers
        for idx in 1:bitlen
            mask = 0x1000 >>> idx
            if (n & mask) != 0
                counts[idx] += 1
            end
        end
    end

    l = div(length(numbers), 2)
    gamma::UInt32 = 0
    for idx in 1:bitlen
        gamma <<= 1
        if counts[idx] >= l
            gamma |= 1
        end
    end
    gamma
end

function part1(bits)
    gamma = calculate_gamma(bits, 12)
    epsilon = xor(gamma, 0x0fff)
    epsilon * gamma
end

function filter_values(values, use_first)
    bit = 0x0800
    while bit != 0 && length(values) > 1
        set = count(v -> (v & bit) != 0 , values)
        if use_first
            good = set >= (length(values) - set) ? bit : 0
        else
            good = set < (length(values)-set) ? bit : 0
        end
        values = filter(b -> (b & bit) == good, values)
        bit >>>= 1
    end
    values[1]
end

function part2(values)
    o2_rating = filter_values(values, true)
    co2_rating = filter_values(values, false)
    o2_rating * co2_rating
end

function parseLine(line)
    parse(UInt32, line, base=2)
end

bits = [parseLine(line) for line in eachline("input/day_3.txt")]
if parsed_args["part2"]
    answer = part2(bits)
else
    answer = part1(bits)
end
println("The answer is: $answer")


text:
export Node, rx, ry, rz, xyz, txyz, scalexyz

"""
    Node([name],x,y,z)
    Node(xyz::Array{Float64,1})

`Node` objects `show` a FastHenry node command.
"""
struct Node <: Element
  name :: Symbol
  xyz1  :: Array{Float64,1}
  function Node(n,xyz::Array{Float64,1})
    xyz1 = Array{Float64}(undef,4)
    xyz1[1:3] = xyz[1:3]
    xyz1[4] = 1.0
    new(Symbol(n),xyz1)
  end
end
Node(xyz::Array{Float64,1}) = Node(Symbol(""),xyz)
Node(name,x,y,z) = Node(name,[x,y,z,1.0])
Node(x,y,z) = Node(Symbol(""),x,y,z)
Node(;name = Symbol(""), x=0, y=0, z=0) = Node(name,x,y,z)

"""
    xyz(n::Node)

Return the coordinate of a node.
"""
xyz(n::Node) = n.xyz1[1:3]

"""
    rx(α) = [1      0       0     0;
             0    cos(α)  sin(α)  0;
             0   -sin(α)  cos(α)  0;
             0      0       0     1]

[Rotation matrix](http://mathworld.wolfram.com/RotationMatrix.html) for angle α around x axis.
"""
rx(α) = [1.0      0.0       0.0     0.0;
         0.0     cos(α)    sin(α)   0.0;
         0.0    -sin(α)    cos(α)   0.0;
         0.0      0.0       0.0     1.0]

"""
    ry(β) = [cos(β)   0    -sin(β)  0;
               0      1       0     0;
             sin(β)   0     cos(β)  0;
               0      0       0     1]

[Rotation matrix](http://mathworld.wolfram.com/RotationMatrix.html) for angle β around y axis.
"""
ry(β) = [cos(β)   0.0    -sin(β)    0.0;
          0.0     1.0      0.0      0.0;
         sin(β)   0.0     cos(β)    0.0;
          0.0     0.0      0.0      1.0]

"""
    rz(γ) = [cos(γ) sin(γ) 0  0;
            -sin(γ) cos(γ) 0  0;
              0       0    1  0;
              0       0    0  1]

[Rotation matrix](http://mathworld.wolfram.com/RotationMatrix.html) for angle γ around z axis.
"""
rz(γ) = [cos(γ) sin(γ) 0.0  0.0;
        -sin(γ) cos(γ) 0.0  0.0;
          0.0    0.0   1.0  0.0;
          0.0    0.0   0.0  1.0]

"""
    txyz(p,q,r) = [1.0 0.0 0.0  p;
                   0.0 1.0 0.0  q;
                   0.0 0.0 1.0  r;
                   0.0 0.0 0.0 1.0]

[Translation matrix](http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/geometry/geo-tran.html).  Translate a coordinate by adding [p,q,r]
"""
txyz(p,q,r) = [1.0 0.0 0.0  p;
               0.0 1.0 0.0  q;
               0.0 0.0 1.0  r;
               0.0 0.0 0.0 1.0]
txyz(a) = txyz(a[1],a[2],a[3])

"""
    scalexyz(p,q,r) = [ p  0.0 0.0 0.0;
                       0.0  q  0.0 0.0;
                       0.0 0.0  r  0.0;
                       0.0 0.0 0.0 1.0]

[Scale matrix](http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/geometry/geo-tran.html).  Scale a coordinate [x, y, z] by [p, q, r].
"""
scalexyz(p,q,r) = [ p  0.0 0.0 0.0;
                   0.0  q  0.0 0.0;
                   0.0 0.0  r  0.0;
                   0.0 0.0 0.0 1.0]


text:
# Eric S. Tellez <eric.tellez@infotec.mx>

#Pkg.add("GZip")
#Pkg.add("Glob")
#Pkg.add("JSON")

using ArgParse
import JSON
import Glob

include("textmodel.jl")
include("io.jl")

function normpatterns(patlist::Vector{Vector{String}}, config::TextConfig)
    A = Vector{String}[]
    for patterns in patlist
        B = String[]
        for pat in patterns
            T = tokenize(pat, config)
            length(T) > 0 && push!(B, join(T, ' '))
        end
        length(B) > 0 && push!(A, B)
    end

    Set(A)
end

function normpatternsset(patlist::Vector{Vector{String}}, config::TextConfig)
    A = Vector{Set{String}}[]

    for patterns in patlist
        B = Set{String}[]
        for pat in patterns
            T = tokenize(pat, config)
            length(T) > 0 && push!(B, Set(T))
        end
        length(B) > 0 && push!(A, B)
    end

    Set(A)
end

function neardup(fun, input, config::TextConfig, key)
    # config.stem = false
    H = Set{String}()
    lineno = 0
    uniques = 0

    itertweets(input) do tweet
        lineno += 1
        text = normtext(getkeypath(tweet, key), config)

        if !(text in H)
            push!(H, text)
            fun(tweet)
            uniques += 1
        end

        if lineno % 100000 == 0
            info("Advance $uniques from $lineno")
        end
    end
end

function main()
    config = TextConfig()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--config"
        help = "a json-formatted dictionary with the TextModel configuration"
        arg_type = String
        default = "{}"
        "--stopwords"
        help = "the file containing the list of stopwords, one per line"
        arg_type = String
        default = ""
        "--key"
        help = "specifies the key pointing to the text"
        arg_type = String
        default = "text"
    end

    args = parse_args(ARGS, s)
    config = TextConfig()
    for (k, v) in JSON.parse(args["config"])
        setfield!(config, Symbol(k), v)
    end
    load(config, args["stopwords"])

    neardup(STDIN, config, args["key"]) do tweet
        write(STDOUT, JSON.json(tweet), '\n')
    end
end

main()

max_stars_repo_name: INGEOTEC/PatternAnalysisTools

max_issues_repo_name: INGEOTEC/PatternAnalysisTools

max_forks_repo_name: INGEOTEC/PatternAnalysisTools


text:
module Paxos

  include("./utils.jl")
  include("./nodes.jl")
  include("./configurations.jl")
  include("./ballots.jl")
  include("./transports.jl")
  include("./ledgers.jl")
  include("./protocols.jl")
  include("./leaders.jl")
  include("./clients.jl")

end # module


text:
# ============================================================================
# Thermostat model
# See https://fenix.tecnico.ulisboa.pt/downloadFile/3779579688470/lygeros.pdf,
# Section 1.3.4.
# ============================================================================

function thermostat()
    c_a = 0.1

    # automaton structure
    automaton = LightAutomaton(2)

    # mode on
    A = hcat(c_a)
    inv = HalfSpace([1.0], 22.0)  # x ≤ 22
    # @vars x v
    # @set x ≤ 22
    m_on = ConstrainedLinearContinuousSystem(A, inv)

    # mode off
    A = hcat(-c_a)
    inv = HalfSpace([-1.0], -18.0)  # x ≥ 18
    m_off = ConstrainedLinearContinuousSystem(A, inv)

    # modes
    modes = [m_on, m_off]

    # transition from on to off
    add_transition!(automaton, 1, 2, 1)
    guard = HalfSpace([-1.0], -21.0)  # x ≥ 21
    t_on2off = ConstrainedIdentityMap(2, guard)

    # transition from off to on
    add_transition!(automaton, 2, 1, 2)
    guard = HalfSpace([1.0], 19.0)  # x ≤ 19
    t_off2on = ConstrainedIdentityMap(2, guard)

    # transition annotations
    resetmaps = [t_on2off, t_off2on]

    # switching
    switchings = [AutonomousSwitching()]

    ℋ = HybridSystem(automaton, modes, resetmaps, switchings)

    # initial condition in mode on
    X0 = Singleton([18.0])
    initial_condition = [(1, X0)]

    return InitialValueProblem(ℋ, initial_condition)
end


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()

module Lumberjack

using Compat
import Base.show, Base.info, Base.log
#import Mocking: @mendable # TODO - figure out how to use Mocking on 0.5

if !isdefined(Base, :StackTraces)
    import StackTraces
end

# To avoid warnings, intentionally do not import:
# Base.error, Base.warn, Base.info

export log,
       debug, info, warn, error,
       add_saw, remove_saw, remove_saws,
       add_truck, remove_truck, remove_trucks,
       configure,

       TimberTruck,
       LumberjackTruck,
       CommonLogTruck,
       JsonTruck,
       Saw,
       FileRoller,
       Syslog

# -------

global _lumber_mill

include("saws.jl")
include("timbertruck.jl")
include("lumbermill.jl")
include("FileRoller.jl")
include("Syslog.jl")

# -------

function __init__()
    global _lumber_mill = LumberMill()
end

# -------

end


text:
"""
    theme(s::Symbol)

Specify the colour theme for plots.
"""
function theme(s::Symbol; kw...)
    # reset?
    if s == :none || s == :default
        PlotUtils.clibrary(:Plots)
        PlotUtils.default_cgrad(default = :sequential, sequential = :inferno)
        default(;
            bg        = :white,
            bglegend  = :match,
            bginside  = :match,
            bgoutside = :match,
            fg        = :auto,
            fglegend  = :match,
            fggrid    = :match,
            fgaxis    = :match,
            fgtext    = :match,
            fgborder  = :match,
            fgguide   = :match,
            palette   = :auto
        )
        return
    end

    # update the default gradient and other defaults
    thm = PlotThemes._themes[s]
    if thm.gradient != nothing
        PlotUtils.clibrary(:misc)
        PlotUtils.default_cgrad(default = :sequential, sequential = PlotThemes.gradient_name(s))
    end
    default(;
        bg       = thm.bg_secondary,
        bginside = thm.bg_primary,
        fg       = thm.lines,
        fgtext   = thm.text,
        fgguide  = thm.text,
        fglegend = thm.text,
        palette  = thm.palette,
        kw...
    )
end

@deprecate set_theme(s) theme(s)