text:
module Base

using JuMP
using DataFrames
using Interpolations
using CSV
using Printf

# These functions are required for MPC_Module.jl and PrettyPlots.jl (resultsDir!)
export
  evalConstraints!,
  opt2dfs!,
  dvs2dfs,
  plant2dfs!,
  postProcess!,
  optimize!,
  scale_tau,
  try_import,
  intervals,
  interpolateLagrange!,
  interpolateLinear!,
  interpolate_lagrange,
  State,
  Control,
  Constraint,
  Results,
  Settings,
  _Ipopt_defaults,
  _Ipopt_MPC,
  # _KNITRO_defaults,
  # _KNITRO_MPC,
  simulationModes,
  resultsDir!

################################################################################
# Constants
################################################################################
#####################################
# settings for both KNITRO and IPOPT
####################################
#outlev # (c.m.solver==:Ipopt) ? :print_level : :outlev # print level
# (KNITRO) absolute stopping tolerance for the feasibility error
# (Ipopt) Absolute tolerance on the constraint violation.
#feastol_abs # (c.m.solver==:Ipopt) ? :constr_viol_tol : :feastol_abs
# (KNITRO) maximum number of iterations before termination
# (Ipopt) Maximum number of iterations.
#maxit # (c.m.solver==:Ipopt) ? :max_iter : :maxit
# (KNITRO) in seconds, the maximum allowable CPU time before termination
# (Ipopt) A limit on CPU seconds that Ipopt can use to solve one problem
#maxtime_cpu # (c.m.solver==:Ipopt) ? :max_cpu_time : :maxtime_cpu
const _Ipopt_defaults=Dict(
  :print_level                =>0,
  :warm_start_init_point      =>"yes",
  :tol                        =>1e-8,
  :max_iter                   =>3000,
  :max_cpu_time               =>1e6,
  :dual_inf_tol               =>1.,
  :constr_viol_tol            =>0.0001,
  :compl_inf_tol              =>0.0001,
  :acceptable_tol             =>1e-6,
  :acceptable_constr_viol_tol =>0.01,
  :acceptable_dual_inf_tol    =>1e-10,
  :acceptable_compl_inf_tol   =>0.01,
  :acceptable_obj_change_tol  =>1e20,
  :diverging_iterates_tol     =>1e20
)

const _Ipopt_MPC=Dict(
  :print_level                =>0,
  :warm_start_init_point      =>"yes",
  :tol                        =>5e-1,
  :max_iter                   =>500,
  :max_cpu_time               =>0.47,
  :dual_inf_tol               =>5.,
  :constr_viol_tol            =>1e-1,
  :compl_inf_tol              =>1e-1,
  :acceptable_tol             =>1e-2,
  :acceptable_constr_viol_tol =>0.01,
  :acceptable_dual_inf_tol    =>1e10,
  :acceptable_compl_inf_tol   =>0.01,
  :acceptable_obj_change_tol  =>1e20,
  :diverging_iterates_tol     =>1e20
)

# const _KNITRO_defaults=Dict(
#  :outlev                       =>1,
#  :feastol                      =>1.0e-6,
#  :feastol_abs                  =>1e-3,
#  :ftol                         =>1e-15,
#  :ftol_iters                   =>5,
#  :infeastol                    =>1.0e-8,
#  :maxfevals                    =>-1,
#  :maxit                        =>0,
#  :maxtime_cpu                  =>1e8,
#  :maxtime_real                 =>1e8,
#  :opttol                       =>1e-6,
#  :opttol_abs                   =>1e-3,
#  :xtol                         =>1e-12,
#  :xtol_iters                   =>0,
#  :algorithm                    =>0,
#  :bar_initpt                   =>0,
#  :bar_murule                   =>0,
#  :bar_penaltycons              =>0,
#  :bar_penaltyrule              =>0,
#  :bar_switchrule               =>0,
#  :linesearch                   =>0,
#  :linsolver                    =>0,
#  :cg_pmem                      =>10,
#  :bar_initpt                   =>0,
#  :bar_penaltycons              =>0,
#  :bar_penaltyrule              =>0,
#  :bar_switchrule               =>0,
#  :linesearch                   =>0,
#  :tuner                        =>0
# )
#
# const _KNITRO_MPC=Dict(
#  :outlev                       =>0,
#  :feastol                      =>1.0e20,
#  :feastol_abs                  =>7e-2,
#  :ftol                         =>1e-15,
#  :ftol_iters                   =>5,
#  :infeastol                    =>1e-2,
#  :maxfevals                    =>-1,
#  :maxit                        =>0,
#  :maxtime_cpu                  =>1e8,
#  :maxtime_real                 =>0.47,
#  :opttol                       =>1.0e20,
#  :opttol_abs                   =>5e-1,
#  :xtol                         =>1e-12,
#  :xtol_iters                   =>0,
#  :algorithm                    =>1,
#  :bar_initpt                   =>0,
#  :bar_murule                   =>1,
#  :bar_penaltycons              =>0,
#  :bar_penaltyrule              =>0,
#  :bar_switchrule               =>0,
#  :linesearch                   =>0,
#  :linsolver                    =>4,
#  :cg_pmem                      =>0,
#  :bar_initpt                   => 3,
#  :bar_penaltycons              => 1,
#  :bar_penaltyrule              => 2,
#  :bar_switchrule               => 2,
#  :linesearch                   => 1,
#  :tuner                        =>0
# )

const simulationModes = [:OCP,:IP,:IPEP,:EP]

################################################################################
# structs
################################################################################
############################### control ########################################
mutable struct Control
  name::Vector{Any}
  description::Vector{Any}
  num::Int64
  pts::Int64
end
function Control()
  Control([],
          [],
          0,
          0)
end
# ############################# state  ##########################################
mutable struct State
  name::Vector{Any}
  description::Vector{Any}
  num::Int64
  pts::Int64
  model
end
function State()
  State([],
        [],
        0,
        0,
        Any)
end

############################## constraint ######################################
mutable struct Constraint
  name::Vector{Any}
  handle::Vector{Any}
  value::Vector{Any}
  nums  # range of indecies in g(x)
end
function Constraint()
  Constraint([],
             [],
             [],
             []);
end

############################### solver  ########################################
mutable struct Solver
    name
    settings
end

function Solver()
       Solver(:Ipopt,
              _Ipopt_defaults);
end

mutable struct PlantResults
  plant
  X0p
  X0a  # closest one in time to X0e
  X0e
  e
  dfsplant
  dfsplantPts::DataFrames.DataFrame # plant data extracted into a single DataFrame
  dfsX0p
  dfsX0a
  dfsX0e
  dfse
end

function PlantResults()
 PlantResults(
    DataFrame(), # plantdata
  Vector{Any}[],
  Vector{Any}[],
  Vector{Any}[],
  Matrix{Any}[], # e
           [],   # vector of plantdata dfs
  DataFrame(),
  DataFrame(),
  DataFrame(),
  DataFrame(),
  DataFrame()
      )
end

mutable struct OCPResults
  tctr                       # time vector for control
  tst                        # time vector for state
  x                          # JuMP states
  u                          # JuMP controls
  X                          # states
  U                          # controls
  X0                         # initial states for OCP
  CS                         # costates
  tpolyPts                   # time sample points for polynomials  (NOTE these interpolated solutions were developed for calulaing error, between them and a known Optimal solution)
  XpolyPts                   # state evaluated using Lagrange/Linear polynomial
  CSpolyPts                  # costate evaluated using Lagrange/Linear polynomial
  UpolyPts                   # control evaluated using Lagrane/Linear polynomial
  AlltpolyPts                   # time sample points for polynomials
  AllXpolyPts                   # state evaluated using Lagrange/Linear polynomial
  AllCSpolyPts                  # costate evaluated using Lagrange/Linear polynomial
  AllUpolyPts                   # control evaluated using Lagrane/Linear polynomial
  tpts                       # vector time sample points
  Xpts                       # vector state sample points
  Upts                       # vector control sample points
  CSpts                      # vector costate sample points
  x0Con                      # handle for initial state constraints
  x0sCon
  xfCon                      # handle for final state constraints
  xfsCon
  dynCon                     # dynamics constraints
  constraint::Constraint     # constraint handles and data
  evalNum                    # number of times optimization has been run
  iterNum                    # mics. data, perhaps an iteration number for a higher level algorithm
  status                     # optimization status
  tSolve                     # solve time for optimization
  objVal                     # objective function value
  dfs                        # results in DataFrame for plotting
  dfsOpt                     # optimization information in DataFrame for plotting
  dfsCon                     # constraint data
end

# Default Constructor
function OCPResults()
OCPResults( Vector{Any}[],# time vector for control
        Vector{Any}[],    # time vector for state
        Matrix{Any}[],    # JuMP states
        Matrix{Any}[],    # JuMP controls
        Matrix{Any}[],    # states
        Matrix{Any}[],    # controls
        Vector{Any}[],    # initial states for OCP
        [],               # costates
        [],               # time sample points for polynomials
        [],               # state evaluated using Lagrange polynomial
        [],               # costate evaluated using Lagrange polynomial
        [],               # control evaluated using Lagrange polynomial
        [],
        [],
        [],
        [],
        Vector{Any}[],    # vector time sample points
        Vector{Any}[],    # vector state sample points
        Vector{Any}[],    # vector control sample points
        Vector{Any}[],    # vector costate sample points
        nothing,          # handle for initial state constraints
        nothing,
        nothing,          # handle for final state constraints
        nothing,
        nothing,          # dynamics constraint
        Constraint(),     # constraint data
        1,                # current evaluation number
        [],               # mics. data, perhaps an iteration number for a higher level algorithm
        Symbol,           # optimization status
        Float64,          # solve time for optimization
        Float64,          # objective function value
        [],               # results in DataFrame for plotting
        DataFrame(),      # optimization information in DataFrame for plotting
        []                # constraint data
      )
end

mutable struct Results
  ocp::OCPResults
  ip::PlantResults
  ep::PlantResults
  resultsDir                 # string that defines results folder
  mainDir                    # string that defines main folder
end

function Results()
  Results(
  OCPResults(),
  PlantResults(),
  PlantResults(),
  string(pwd(),"/results/"),  # string that defines results folder
  pwd()                       # string that defines main folder
  )
end

# MPC Settings Class
mutable struct MPCSettings
 on::Bool
 mode::Symbol
 predictX0::Bool
 fixedTp::Bool
 IPKnown::Bool
 saveMode::Symbol
 maxSim::Int64                # maximum number of total MPC updates
 shiftX0::Bool                # a bool to indicate that a check on the feasibility of the linear constraints for the initial conditions should be performed
 lastOptimal::Bool            # a bool to indicate that the previous solution should be used if the current solution in not optimal
 printLevel::Int64
 expandGoal::Bool             # bool to indicate if the goal needs to be expanded
 enlargeGoalTolFactor::Int64  # scaling factor to enlare the goal
 onlyOptimal::Bool
end

function MPCSettings()
      MPCSettings(
      false,
      :OCP,
      false,
      true,
      true,
      :all,
      100,
      true,
      true,
      2,
      true,
      2,
      false
      )
end

# OCP Settings Class
mutable struct OCPSettings
  solver::Solver                # solver information
  finalTimeDV::Bool
  integrationMethod::Symbol
  integrationScheme::Symbol
  save::Bool                    # bool for saving data
  reset::Bool                   # bool for reseting data
  evalConstraints::Bool         # bool for evaluating duals of the constraints
  evalCostates::Bool            # bool for evaluating costates
  tfMin::Any                    # minimum final time
  tfMax::Any                    # maximum final time
  tfOptimal::Any                # known optimal final time
  numInterpPts::Int64           # number of points to sample polynomial running through collocation points
  cacheOnly::Bool               # bool for only caching the results when using optimize!()
  linearInterpolation::Bool     # bool for using linear interpolation even if integrationMethod ==:ps
  interpolationOn::Bool         # bool to indicate if user wants solution interpolated for them
  x0slackVariables::Bool
  xFslackVariables::Bool
end

# Default Constructor NOTE currently not using these, they get overwritten
function OCPSettings()
      OCPSettings(
         Solver(),           # default solver
         false,              # finalTimeDV
         :ts,                # integrationMethod
         :bkwEuler,          # integrationScheme
         true,               # bool for saving data
         false,              # bool for reseting data
         false,              # bool for evaluating duals of the constraints
         false,              # bool for evaluating costates
         0.001,              # minimum final time
         400.0,              # maximum final time
         false,              # known optimal final time
         250,                # number of points to sample polynomial running through collocation points
         false,              # bool for only caching the results when using optimize!()
         false,              # bool for using linear interpolation even if integrationMethod ==:ps
         false,              # bool to indicate if user wants solution interpolated for them
         false,
         false
                )
end

mutable struct Settings
 ocp::OCPSettings
 mpc::MPCSettings
end

function Settings()
 Settings(
  OCPSettings(),
  MPCSettings()
       )
end

"""
L = lagrange_basis_poly(x,x_data,Nc,j)
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 12/26/2016, Last Modified: 1/25/2016
Citations: This function was influenced by the lagrange() function [located here](https://github.com/pjabardo/Jacobi.jl/blob/master/src/gauss_quad.jl)
--------------------------------------------------------------------------------------\n
# Input Arguments
* `x`: point to approximate function value at
* `x_data`: x data to used calculate basis polynomials
* `Nc`: order of Lagrange interpolating polynomial
* `j`: index of interest

# Output Arguments
* `L`: Lagrange basis polynomials

A basic description of Lagrange interpolating polynomials is provided [here](http://127.0.0.1:8000/lagrange_poly.html#lagrange-poly)

"""
function lagrange_basis_poly(x,x_data,Nc,j)
    L = 1;
    for idx in 1:Nc+1 # use all of the data
      if idx!=j
        L = L*(x - x_data[idx])/(x_data[j]-x_data[idx]);
      end
    end
  return L
end

"""
y=interpolate_lagrange(ts[int],ts[int],stateMatrix[int][:,st])
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Created: 1/2/2017, Last Modified: 9/19/2017 \n
--------------------------------------------------------------------------------------\n
"""
function interpolate_lagrange(x::AbstractArray{T},x_data,y_data) where {T<:Number}
  Nc = length(x_data) - 1

  if length(x_data)!=length(y_data)
      error(string("\n",
                    "-------------------------------------------------------", "\n",
                    "There is an error with the data vector lengths!!", "\n",
                    "-------------------------------------------------------", "\n",
                    "The following variables should be equal:", "\n",
                    "length(x_data) = ",length(x_data),"\n",
                    "length(y_data) = ",length(y_data),"\n"
                    )
            )
    end
    ns = length(x)
    L = zeros(Float64,Nc+1,ns)
    x = x[:]; x_data = x_data[:]; y_data = y_data[:]; # make sure data is in a column
    for idx in 1:Nc+1
      for j in 1:ns
          L[idx,j] = lagrange_basis_poly(x[j],x_data,Nc,idx);
      end
    end
    return y_data'*L
end

"""
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 6/29/2017, Last Modified: 12/06/2019 \n
--------------------------------------------------------------------------------------\n
"""
function intervals(n,int,x,u)

  if typeof(x[1,1]) == JuMP.Variable
    # states
    x_int = Array{JuMP.Variable}(length(n.ocp.Nck_full[int]+1:n.ocp.Nck_full[int+1]),n.ocp.state.num)
    for st in 1:n.ocp.state.num # +1 adds the DV in the next interval
      x_int[:,st] = x[n.ocp.Nck_cum[int]+1:n.ocp.Nck_cum[int+1]+1,st];
    end

    # controls
    if int!=n.ocp.Ni
      u_int = Matrix{JuMP.Variable}(undef, length(n.ocp.Nck_full[int]+1:n.ocp.Nck_full[int+1]),n.ocp.control.num)
    else                    # -1 -> removing control in last mesh interval
      u_int = Matrix{JuMP.Variable}(undef, length(n.ocp.Nck_full[int]+1:n.ocp.Nck_full[int+1]-1),n.ocp.control.num)
    end
    for ctr in 1:n.ocp.control.num
      if int!=n.ocp.Ni          # +1 adds the DV in the next interval
        u_int[:,ctr] = u[n.ocp.Nck_cum[int]+1:n.ocp.Nck_cum[int+1]+1,ctr];
      else
        u_int[:,ctr] = u[n.ocp.Nck_cum[int]+1:n.ocp.Nck_cum[int+1],ctr];
      end
    end
  else
    # states
    x_int=Matrix{Any}(undef, length(n.ocp.Nck_full[int]+1:n.ocp.Nck_full[int+1]),n.ocp.state.num);
    for st in 1:n.ocp.state.num # +1 adds the DV in the next interval
      x_int[:,st] = x[n.ocp.Nck_cum[int]+1:n.ocp.Nck_cum[int+1]+1,st];
    end

    # controls
    if int!=n.ocp.Ni
      u_int = Matrix{Any}(undef, length(n.ocp.Nck_full[int]+1:n.ocp.Nck_full[int+1]),n.ocp.control.num)
    else                    # -1 -> removing control in last mesh interval
      u_int = Matrix{Any}(undef, length(n.ocp.Nck_full[int]+1:n.ocp.Nck_full[int+1]-1),n.ocp.control.num)
    end
    for ctr in 1:n.ocp.control.num
      if int!=n.ocp.Ni          # +1 adds the DV in the next interval
        u_int[:,ctr] = u[n.ocp.Nck_cum[int]+1:n.ocp.Nck_cum[int+1]+1,ctr]
      else
        u_int[:,ctr] = u[n.ocp.Nck_cum[int]+1:n.ocp.Nck_cum[int+1],ctr]
      end
    end
  end

  return x_int,u_int
end

"""
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Created: 9/19/2017, Last Modified: 12/06/2019 \n
--------------------------------------------------------------------------------------\n
"""
function interpolateLagrange!(n; numPts::Int64=250, tfOptimal::Any=false)
  # TODO throw an error if tfOptimal does not make sense given current solution
  if isa(tfOptimal,Bool)
    if n.s.ocp.finalTimeDV
      tf = getvalue(n.ocp.tf)
    else n.s.ocp.finalTimeDV
      tf = n.ocp.tf
    end
  else  # if there is a known optimal final time, then it is useful to evaluate the Lagrange polynomial at particular points to determine the error in the solution
    tf = tfOptimal
  end

  if isnan(tf)
      @warn "tf is a NaN cannot use it in interpolateLagrange!().\n
            Exiting interpolateLagrange!() without an interpolated solution."
    return nothing
  end

  if tf < 0.01
    @warn "tf needs to be greater than 0.01 to interpolate over solution.\n
          Exiting interpolateLagrange!() without an interpolated solution."
     return nothing
  end

  # sample points
  n.r.ocp.tpolyPts = [range(tf/n.ocp.Ni*(int-1),tf/n.ocp.Ni*int;length=numPts) .+ n.r.ocp.tst[1] for int in 1:n.ocp.Ni]
  n.r.ocp.XpolyPts = [[zeros(numPts) for int in 1:n.ocp.Ni] for st in 1:n.ocp.state.num]
  n.r.ocp.UpolyPts = [[zeros(numPts) for int in 1:n.ocp.Ni] for ctr in 1:n.ocp.control.num]
  if n.s.ocp.evalCostates; n.r.ocp.CSpolyPts = [[zeros(numPts) for int in 1:n.ocp.Ni] for st in 1:n.ocp.state.num] end

  # time data points
  t_st_int = [n.r.ocp.tst[n.ocp.Nck_cum[int]+1:n.ocp.Nck_cum[int+1]+1] for int in 1:n.ocp.Ni]
  t_ctr_int = [n.r.ocp.tctr[n.ocp.Nck_cum[int]+1:n.ocp.Nck_cum[int+1]+1] for int in 1:n.ocp.Ni-1]
  t_ctr_int = push!(t_ctr_int, n.r.ocp.tst[n.ocp.Nck_cum[n.ocp.Ni]+1:n.ocp.Nck_cum[n.ocp.Ni+1]]) # -1 -> removing control in last mesh interval

  for int in 1:n.ocp.Ni
    # controls and states for this interval
    x_int, u_int = intervals(n, int, copy(n.r.ocp.X), copy(n.r.ocp.U))

     # sample polynomial in interval at n.r.ocp.tpolyPts
     for st in 1:n.ocp.state.num
      n.r.ocp.XpolyPts[st][int] = interpolate_lagrange(n.r.ocp.tpolyPts[int], t_st_int[int], x_int[:,st])'
     end

     for ctr in 1:n.ocp.control.num
      n.r.ocp.UpolyPts[ctr][int] = interpolate_lagrange(n.r.ocp.tpolyPts[int], t_ctr_int[int], u_int[:,ctr])'
     end

      # sample polynomial in interval at n.r.ocp.tpolyPts NOTE costate is missing the last point, that is the t_st_int[int][1:end-1]
      if n.s.ocp.evalCostates && n.s.ocp.evalConstraints
         for st in 1:n.ocp.state.num
          n.r.ocp.CSpolyPts[st][int] = interpolate_lagrange(n.r.ocp.tpolyPts[int], t_st_int[int][1:end-1], n.r.ocp.CS[st][int])'
         end
      end
  end

  # extract result into vectors
  temp = [n.r.ocp.tpolyPts[int][1:end] for int in 1:n.ocp.Ni] # time
  n.r.ocp.tpts = [idx for tempM in temp for idx=tempM]
  totalPts = length(n.r.ocp.tpts)

  n.r.ocp.Xpts = Matrix{Float64}(undef, totalPts, n.ocp.state.num)
  for st in 1:n.ocp.state.num # states
    temp = [n.r.ocp.XpolyPts[st][int][1:end,:] for int in 1:n.ocp.Ni]
    n.r.ocp.Xpts[:,st] = [idx for tempM in temp for idx=tempM]
  end

  n.r.ocp.Upts = Matrix{Float64}(undef, totalPts, n.ocp.control.num)
  for ctr in 1:n.ocp.control.num # controls
    temp = [n.r.ocp.UpolyPts[ctr][int][1:end,:] for int in 1:n.ocp.Ni]
    n.r.ocp.Upts[:,ctr] = [idx for tempM in temp for idx=tempM]
  end

  if n.s.ocp.evalCostates && n.s.ocp.evalConstraints
    n.r.ocp.CSpts = Matrix{Float64}(undef, totalPts, n.ocp.state.num)
    for st in 1:n.ocp.state.num # states
      temp = [n.r.ocp.CSpolyPts[st][int][1:end,:] for int in 1:n.ocp.Ni]
      n.r.ocp.CSpts[:,st] = [idx for tempM in temp for idx=tempM]
    end
  end

  return nothing
end

"""
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Created: 10/4/2017, Last Modified: 12/06/2019 \n
--------------------------------------------------------------------------------------\n
"""
function interpolateLinear!(n; numPts::Int64=250, tfOptimal::Any=false)
  # TODO throw an error if tfOptimal does not make sense given current solution
  if isa(tfOptimal,Bool)
    if n.s.ocp.finalTimeDV
      tf = getvalue(n.ocp.tf)
    else n.s.ocp.finalTimeDV
      tf = n.ocp.tf
    end
  else  # if there is a known optimal final time, then it is useful to evaluate the Lagrange polynomial at particular points to determine the error in the solution
    tf = tfOptimal
  end

  if isnan(tf)
      @warn "tf is a NaN cannot use it in interpolateLinear!().\n
            Exiting interpolateLinear!() without an interpolated solution."
      return nothing
  end

  if tf < 0.01
    @warn "tf needs to be greater than 0.01 to interpolate over solution.\n
          Exiting interpolateLinear!() without an interpolated solution."
    return nothing
  end

  # sample points
  t = range(0,tf,length=numPts) .+ n.r.ocp.tst[1] # NOTE not tested
  n.r.ocp.tpts = convert(Array{Float64,1},t)
  n.r.ocp.Xpts = Matrix{Float64}(undef, numPts, n.ocp.state.num)
  n.r.ocp.Upts = Matrix{Float64}(undef, numPts, n.ocp.control.num)
  knots = (n.r.ocp.tst,)
  for st in 1:n.ocp.state.num
    sp_st = interpolate(knots,n.r.ocp.X[:,st],Gridded(Linear()))
    n.r.ocp.Xpts[:,st] = sp_st[n.r.ocp.tpts]
  end

  for ctr in 1:n.ocp.control.num
    if isequal(n.s.ocp.integrationMethod,:ps)
      sp_ctr = interpolate(knots,[n.r.ocp.U[:,ctr];0],Gridded(Linear()))
    else
      sp_ctr = interpolate(knots,n.r.ocp.U[:,ctr],Gridded(Linear()))
    end
    n.r.ocp.Upts[:,ctr] = sp_ctr[n.r.ocp.tpts]
  end
  return nothing
end
"""
scale_tau(tau,ta,tb)
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 12/23/2017, Last Modified: 12/06/2019 \n
--------------------------------------------------------------------------------------\n
"""
function scale_tau(tau,ta,tb)
  (tb - ta)/2*tau .+ (ta + tb)/2
end

"""
plant2dfs!(n,sol)
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 2/14/2017, Last Modified: 12/06/2019 \n
--------------------------------------------------------------------------------------\n
"""
function plant2dfs!(n,sol,U)
  tSample = range(sol.t[1],sol.t[end],length=n.mpc.ip.state.pts)
  dfs = DataFrame()
  if isempty(n.r.ip.plant)
    dfs[:t] = tSample
    for st in 1:n.mpc.ip.state.num
      dfs[n.mpc.ip.state.name[st]] = [sol(t)[st] for t in tSample]
    end
    for ctr in 1:n.mpc.ip.control.num
      dfs[n.mpc.ip.control.name[ctr]] = [U[ctr][t] for t in tSample]
    end
  else
    dfs[:t] = [n.r.ip.plant[:t]; tSample]
    for st in 1:n.mpc.ip.state.num
      dfs[n.mpc.ip.state.name[st]] = [n.r.ip.plant[n.mpc.ip.state.name[st]]; [sol(t)[st] for t in tSample]]
    end
    for ctr in 1:n.mpc.ip.control.num
      dfs[n.mpc.ip.control.name[ctr]] = [n.r.ip.plant[n.mpc.ip.control.name[ctr]]; [U[ctr][t] for t in tSample] ]
    end
  end
  n.r.ip.plant = dfs

  # TODO only run this if saving, it is redundant data, or maybe put time stamps or something in the above data
  dfs = DataFrame()
  dfs[:t] = tSample
  for st in 1:n.mpc.ip.state.num
    dfs[n.mpc.ip.state.name[st]] = [sol(t)[st] for t in tSample]
  end
  for ctr in 1:n.mpc.ip.control.num
    dfs[n.mpc.ip.control.name[ctr]] = [U[ctr][t] for t in tSample]
  end
  push!(n.r.ip.dfsplant,dfs)
  return nothing
end

"""
dvs2dfs(n)
# funtionality to save state, costate, and control data from optimization
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 2/10/2017, Last Modified: 12/06/2019 \n
--------------------------------------------------------------------------------------\n
"""
function dvs2dfs(n)
  dfs = DataFrame()
  dfs[:t] = n.r.ocp.tst
  for st in 1:n.ocp.state.num
    dfs[n.ocp.state.name[st]] = n.r.ocp.X[:,st]
  end
  for ctr in 1:n.ocp.control.num
    if n.s.ocp.integrationMethod==:tm
      dfs[!,n.ocp.control.name[ctr]] = n.r.ocp.U[:,ctr]
    else
      dfs[!,n.ocp.control.name[ctr]] = [n.r.ocp.U[:,ctr];NaN]
    end
  end

  if n.s.ocp.evalCostates && n.s.ocp.integrationMethod == :ps && n.s.ocp.evalConstraints
    CS_vector = Matrix{Float64}(undef, n.ocp.state.pts, n.ocp.state.num)
    for st in 1:n.ocp.state.num # states
      temp = [n.r.ocp.CS[st][int][1:end,:] for int in 1:n.ocp.Ni]
      CS_vector[1:end-1,st] = [idx for tempM in temp for idx=tempM]
      CS_vector[end,st] = NaN
    end
    for st in 1:n.ocp.state.num # states
      dfs[Symbol(n.ocp.state.name[st],:cs)] = CS_vector[:,st]
    end
  end

  return dfs
end

"""
opt2dfs!(n)
# funtionality to save optimization data
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 2/10/2017, Last Modified: 4/13/2018 \n
--------------------------------------------------------------------------------------\n
"""
function opt2dfs!(n;kwargs...)
  kw = Dict(kwargs);

  if !haskey(kw,:statusUpdate); statusUpdate=false;
  else; statusUpdate=get(kw,:statusUpdate,0);
  end

  # make sure that the feildnames are initialized
  if isempty(n.r.ocp.dfsOpt)
      n.r.ocp.dfsOpt = DataFrame(tSolve = [], objVal = [], status = [], iterNum = [], evalNum = [])
  end

  if !statusUpdate
    push!(n.r.ocp.dfsOpt[:tSolve], n.r.ocp.tSolve)
    push!(n.r.ocp.dfsOpt[:objVal], n.r.ocp.objVal)
    push!(n.r.ocp.dfsOpt[:status], n.r.ocp.status)
  else  # TODO consider removing and cleaning this up
    push!(n.r.ocp.dfsOpt[:tSolve], NaN)
    push!(n.r.ocp.dfsOpt[:objVal], NaN)
    if statusUpdate && (typeof(n.r.ocp.status)==Symbol)
      push!(n.r.ocp.dfsOpt[:status], n.r.ocp.status)
    else
      push!(n.r.ocp.dfsOpt[:status], NaN)
    end
  end
  if n.s.mpc.on
    push!(n.r.ocp.dfsOpt[:iterNum], n.r.ocp.iterNum) # can set iter_nums for a higher/lower level algoritm
  end
  push!(n.r.ocp.dfsOpt[:evalNum], n.r.ocp.evalNum-1)

  return nothing
end

"""
con2dfs(n)
# funtionality to save constraint data
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 2/20/2017, Last Modified: 4/13/2018 \n
--------------------------------------------------------------------------------------\n
"""
function con2dfs(n)
  dfs_con=DataFrame()
  dfs_con[:conVal]=n.r.ocp.constraint.value;
  return dfs_con
end

"""
postProcess!(n)
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 1/27/2017, Last Modified: 12/06/2019 \n
--------------------------------------------------------------------------------------\n
"""
function postProcess!(n;kwargs...)
  kw = Dict(kwargs)
  # check to see if the user is initializing while compensating for control delay
  if !haskey(kw,:Init); Init = false;
  else; Init = get(kw,:Init,0);
  end

  if n.s.ocp.save
    opt2dfs!(n)
  end

  # even if n.r.ocp.status==:Infeasible try to get solution. For the case that user may want to look at results to see where constraints where violated
  # in this case set =>  n.s.ocp.evalConstraints = true
  # http://jump.readthedocs.io/en/latest/refmodel.html#solve-status
  if !Init #&& (n.s.ocp.evalConstraints || ((n.r.ocp.status==:Optimal) || (n.r.ocp.status==:UserLimit)))
    if n.s.ocp.integrationMethod==:ps
      if n.s.ocp.finalTimeDV
        t = [scale_tau(n.ocp.ts[int],0.0,getvalue(n.ocp.tf)) for int in 1:n.ocp.Ni]     # scale time from [-1,1] to [t0,tf]
      else
        t = [scale_tau(n.ocp.ts[int],0.0,n.ocp.tf) for int in 1:n.ocp.Ni]
      end
      n.r.ocp.tctr = [idx for tempM in t for idx = tempM[1:end-1]] .+ getvalue(n.ocp.t0)
      n.r.ocp.tst = [n.r.ocp.tctr; t[end][end] .+ getvalue(n.ocp.t0)]
      # TODO check the above line... is t0 getting added on twice?
    elseif n.s.ocp.integrationMethod==:tm
      if n.s.ocp.finalTimeDV
        n.r.ocp.tctr = append!([0.0],cumsum(getvalue(n.ocp.dt))) .+ getvalue(n.ocp.t0)
      else
        n.r.ocp.tctr = append!([0.0],cumsum(n.ocp.dt)) .+ getvalue(n.ocp.t0)
      end
      n.r.ocp.tst = n.r.ocp.tctr
    end

    stateDataExists = false
    if n.r.ocp.status==:Optimal || (!n.s.mpc.onlyOptimal && n.s.mpc.on && isequal(n.mpc.v.evalNum,1))
      stateDataExists = true
      if (!isequal(n.r.ocp.status,:Optimal) && n.s.mpc.on && isequal(n.mpc.v.evalNum,1))
        @warn "There is no previous :Optimal solution to use since isequal(n.mpc.v.evalNum,1). \n
            Attemting to extract: ",n.r.ocp.status," solution. \n
            Setting: n.f.mpc.simFailed = [true, n.r.ocp.status] "
         n.f.mpc.simFailed = [true, n.r.ocp.status]
      end
      n.r.ocp.X = zeros(Float64,n.ocp.state.pts,n.ocp.state.num)
      n.r.ocp.U = zeros(Float64,n.ocp.control.pts,n.ocp.control.num)
      for st in 1:n.ocp.state.num
        n.r.ocp.X[:,st] = getvalue(n.r.ocp.x[:,st])
      end
      for ctr in 1:n.ocp.control.num
        n.r.ocp.U[:,ctr] = getvalue(n.r.ocp.u[:,ctr])
      end

    elseif n.s.mpc.on && n.s.mpc.lastOptimal && !n.s.mpc.onlyOptimal
      if !n.s.ocp.save
        error("This functionality currently needs to have n.s.ocp.save==true")
      end
      optIdx = findall(n.r.ocp.dfsOpt[:status].==:Optimal)[end]  # use the last :Optimal solution
      @show n.r.ocp.dfsOpt
      @show optIdx
      @show n.r.ocp.dfs[optIdx]
      @show n.mpc.v.t
      if n.r.ocp.dfs[optIdx][:t][1] > n.mpc.v.t
        timeIdx = 1
      elseif n.r.ocp.dfs[optIdx][:t][end] < n.mpc.v.t
        @warn "the current time is past the final time in the last :Optimal soultion.\n
              Setting: n.f.mpc.simFailed = [true, n.r.ocp.status]"
        n.f.mpc.simFailed = [true, n.r.ocp.status]
        return nothing
      else
        timeIdx = findall(n.r.ocp.dfs[optIdx][:t] .- n.mpc.v.t .<= 0)[end]     # find the nearest index in time
      end
      # TODO make an error message or fix      ERROR: LoadError: BoundsError: attempt to access 0-element Array{Int64,1} at index [0]
      n.r.ocp.tst = n.r.ocp.dfs[optIdx][:t][timeIdx:end]
      n.r.ocp.X = zeros(Float64,length(n.r.ocp.dfs[optIdx][n.ocp.state.name[1]][timeIdx:end]),n.ocp.state.num)
      if n.s.ocp.integrationMethod==:tm  # TODO try to
        n.r.ocp.tctr = n.r.ocp.dfs[optIdx][:t][timeIdx:end]
        n.r.ocp.U = zeros(Float64,length(n.r.ocp.dfs[optIdx][:,n.ocp.control.name[1]][timeIdx:end]),n.ocp.control.num)
      else
        n.r.ocp.tctr = n.r.ocp.dfs[optIdx][:t][timeIdx:end-1]
        n.r.ocp.U = zeros(Float64,length(n.r.ocp.dfs[optIdx][:,n.ocp.control.name[1]][timeIdx:end-1]),n.ocp.control.num)
      end
      for st in 1:n.ocp.state.num
        n.r.ocp.X[:,st] = n.r.ocp.dfs[optIdx][n.ocp.state.name[st]][timeIdx:end]
      end
      for ctr in 1:n.ocp.control.num
        if n.s.ocp.integrationMethod==:tm
          n.r.ocp.U[:,ctr] = n.r.ocp.dfs[optIdx][:,n.ocp.control.name[ctr]][timeIdx:end]
        else
          n.r.ocp.U[:,ctr] = n.r.ocp.dfs[optIdx][:,n.ocp.control.name[ctr]][timeIdx:end-1]
        end
      end
    else
      @warn "The solution is not Optimal \n
          Setting: n.f.mpc.simFailed = [true, n.r.ocp.status] "
       n.f.mpc.simFailed = [true, n.r.ocp.status]
    end
    if n.s.ocp.evalConstraints && n.r.ocp.status!=:Error  # note may want to remove the && arg
      evalConstraints!(n)
      # TODO make a note that costates can only be evaluated if .....
      if n.s.ocp.evalCostates && n.s.ocp.integrationMethod == :ps
        L1 = 0       # find index where dynamics constraints start
        for i in 1:length(n.r.ocp.constraint.name)
          if n.r.ocp.constraint.name[i] == :dyn_con
            L1 = n.r.ocp.constraint.nums[i][end][1]
          end
        end
        mpb = JuMP.internalmodel(n.ocp.mdl)
        c = MathProgBase.getconstrduals(mpb)
        # NOTE for now since costates are not defined for :tm methods, n.r.ocp.CS is in a different format than n.r.ocp.X
        # in the future if costates are defined for :tm methods this can be changed
        n.r.ocp.CS = [[zeros(Float64,n.ocp.Nck[int]) for int in 1:n.ocp.Ni] for st in 1:n.ocp.state.num]
        for int in 1:n.ocp.Ni
          b = 0
          for st in 1:n.ocp.state.num
            a = L1 + n.ocp.Nck[int]*(st-1)  # n.ocp.Nck[int]*(st-1) adds on indices for each additional state within the interval
            b = a + n.ocp.Nck[int] - 1      # length of the state within this interval
            n.r.ocp.CS[st][int] = -c[a:b]./n.ocp.ws[int]
          end
          L1 = b + 1 # adds indicies due to a change in the interval
        end
      end
    end
    if n.s.ocp.save && stateDataExists
      push!(n.r.ocp.dfs,dvs2dfs(n))
      push!(n.r.ocp.dfsCon,con2dfs(n))
      if n.s.ocp.interpolationOn
        if (n.r.ocp.status != :Error) # (n.r.ocp.status != :Infeasible) &&
          if n.s.ocp.integrationMethod==:ps && (n.r.ocp.status != :Infeasible) && !n.s.ocp.linearInterpolation
            interpolateLagrange!(n;numPts = n.s.ocp.numInterpPts, tfOptimal = n.s.ocp.tfOptimal)
            push!(n.r.ocp.AlltpolyPts,n.r.ocp.tpolyPts)
            push!(n.r.ocp.AllXpolyPts,n.r.ocp.XpolyPts)
            push!(n.r.ocp.AllUpolyPts,n.r.ocp.UpolyPts)
            if n.s.ocp.evalCostates && n.s.ocp.evalConstraints
              push!(n.r.ocp.AllCSpolyPts,n.r.ocp.CSpolyPts)
            end
          else
            interpolateLinear!(n;numPts = n.s.ocp.numInterpPts, tfOptimal = n.s.ocp.tfOptimal)
            if n.s.ocp.integrationMethod==:ps && !n.s.ocp.linearInterpolation
              push!(n.r.ocp.AlltpolyPts,nothing)
              push!(n.r.ocp.AllXpolyPts,nothing)
              push!(n.r.ocp.AllUpolyPts,nothing)
              if n.s.ocp.evalCostates && n.s.ocp.evalConstraints
                push!(n.r.ocp.AllCSpolyPts,nothing)
              end
            end
          end
        end
      end  # TODO save [] of interpolated pts for :tm methods
    end
  end
  return nothing
end
# TODO add some "finalPostProcess" to sace time
"""
optimize!(n)
# solves JuMP model and saves optimization data
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 2/6/2017, Last Modified: 4/13/2018 \n
--------------------------------------------------------------------------------------\n
"""
function optimize!(n;Iter::Int64=0)
  t1 = time(); status = JuMP.solve(n.ocp.mdl); t2 = time();
  if !n.s.ocp.cacheOnly
    n.r.ocp.status = status
    n.r.ocp.tSolve = t2 - t1
    n.r.ocp.objVal = getobjectivevalue(n.ocp.mdl)
    n.r.ocp.iterNum = Iter    # possible iteration number for a higher level algorithm
    n.r.ocp.evalNum = n.r.ocp.evalNum + 1
    postProcess!(n)      # temporarily save data
  end
  return nothing
end
"""
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 2/7/2017, Last Modified: 4/13/2018 \n
--------------------------------------------------------------------------------------\n
"""
function evalConstraints!(n)
  n.r.ocp.constraint.value=[];   # reset values
  n.r.ocp.constraint.nums=[]; s=1;
  for i = 1:length(n.r.ocp.constraint.handle)
    if n.r.ocp.constraint.name[i]==:dyn_con  # state constraits
      dfs=Vector{DataFrame}(n.ocp.state.num);
      con=DataFrame(step=1);
      l=0;
      for st in 1:n.ocp.state.num
        if n.s.ocp.integrationMethod==:ps
          temp=[getdual(n.r.ocp.constraint.handle[i][int][:,st]) for int in 1:n.ocp.Ni];
          vals=[idx for tempM in temp for idx=tempM];
          dfs[st]=DataFrame(step=1:sum(n.ocp.Nck);Dict(n.ocp.state.name[st] => vals)...);
          l=l+length(vals);
        else
          dfs[st]=DataFrame(step=1:n.ocp.N;Dict(n.ocp.state.name[st] => getdual(n.r.ocp.constraint.handle[i][:,st]))...);
          l=l+length(n.r.ocp.constraint.handle[i][:,st]);
        end
        if st==1; con=dfs[st]; else; con=join(con,dfs[st],on=:step); end
      end
    else
      S=0;
      try
        S=JuMP.size(n.r.ocp.constraint.handle[i])
      catch
        error("\n For now, the constraints cannot be in this form: \n
        con=@NLconstraint(mdl,n.r.ocp.u[1,1]==param); \n
        Write it in array form: \n
          con=@NLconstraint(mdl,[i=1],n.r.ocp.u[i,1]==param); \n")
      end
      if length(S)==1
        con = DataFrame(step=1:length(n.r.ocp.constraint.handle[i]);Dict(n.r.ocp.constraint.name[i] => getdual(n.r.ocp.constraint.handle[i][:]))...);
        l=S[1];
      elseif length(S)==2
        dfs=Vector{DataFrame}(S[1]);
        con=DataFrame(step=1);
        for idx in 1:S[1]
          try
            dfs[idx] = DataFrame(step=1:S[2];Dict(n.r.ocp.constraint.name[i] => getdual(n.r.ocp.constraint.handle[i][idx,:]))...);
          catch
            dfs[idx] = DataFrame(step=1:S[2];Dict(n.r.ocp.constraint.name[i] => NaN)...); # fix for case where all of the states are not being constrainted, but some are within some XF_tol
          end
          if idx==1;
            con = dfs[idx]
          else
            con = join(con,dfs[idx],on=:step,makeunique=true)
          end
        end
        l = S[1]*S[2]
      end
    end
    f = s + l - 1
    num = (i,n.r.ocp.constraint.name[i],@sprintf("length = %0.0f",l),string("indecies in g(x) = "),(s,f))
    push!(n.r.ocp.constraint.nums,num)
    push!(n.r.ocp.constraint.value,con)
    s = f + 1
  end
  return nothing
end

"""
description = string(
" *  \n ")

resultsDir!(r;resultsName=resultsName,description=description)
# removes results folder and creates a new one
# TODO consider putting in a warning or some sort of an interaction with user
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 3/26/2017, Last Modified: 5/29/2017 \n
--------------------------------------------------------------------------------------\n
"""
function resultsDir!(n;resultsName::String = "",description::DataFrame = DataFrame())
 results_dir=string(n.r.mainDir,"/results/",resultsName)  # define directories
 n.r.resultsDir=results_dir;

 if isdir(n.r.resultsDir)
   rm(n.r.resultsDir; recursive=true)
   print("\n The old results have all been deleted! \n \n")
 end
 mkdir(n.r.resultsDir)# create directory

 cd(n.r.resultsDir)
   CSV.write("description.csv", description; quotechar = ' ')
 cd(n.r.mainDir)
 return nothing
end

end # module


text:
module inputs
export ρ0, dt, g0, tmax, t0, β, y0, p0, R, m,r, calibre, v0

const ρ0 = 1.225 #kg/m3
const g0 = 9.80665 #m/s2
const R = 287.05287
const tmax = 1.0 #s
const dt = 0.01 #s
const t0 = 288.15#K ou 15degC
const β = -6.5 #K/km
const y0 = 0.0 #m
const p0 = 101325#Pa
const m = 0.423#g
const r = 6356766# m
const calibre = 0.03 #m
const v0 = 900.0 #m/s

R  =
Ω =
L = 494.7 #mm
AZ =
S =
d = 104.8 # mm
#ρ =
ρ0 = 1.225 #kg/m3
m = 14.97 #kg
Iy = 0.23118 #kg/m2
Ix = 0.02326 #kg/m2
δF =
g = 32.174*(1-0.0026*cos(2*L))*0.3048
n = 38.96
CG = 183.4 #mm from base
T0 = 20 #degrees

end

ext: jl


text:
if record_reference_images
    cd(homedir()) do
        isdir(dirname(recordpath)) || run(`git clone git@github.com:SimonDanisch/ReferenceImages.git`)
        isdir(recordpath) && rm(recordpath, recursive = true, force = true)
        mkdir(recordpath)
    end
end


text:
# This is an example on training the TurbulenceConvection.jl implementation
# of the EDMF scheme with data generated using PyCLES (an LES solver) or
# TurbulenceConvection.jl (perfect model setting).
#
# This example parallelizes forward model evaluations using Julia's pmap() function.
# This parallelization strategy is efficient for small ensembles (N_ens < 20). For
# calibration processes involving a larger number of ensemble members, parallelization
# using HPC resources directly is advantageous.

# Import modules to all processes
using Distributed
@everywhere using Pkg
@everywhere Pkg.activate(dirname(dirname(dirname(@__DIR__))))
@everywhere begin
    using ArgParse
    using CalibrateEDMF
    using CalibrateEDMF.DistributionUtils
    using CalibrateEDMF.Pipeline
    src_dir = dirname(pathof(CalibrateEDMF))
    using CalibrateEDMF.HelperFuncs
    # Import EKP modules
    using EnsembleKalmanProcesses
    using EnsembleKalmanProcesses.ParameterDistributions
end
# include(joinpath(@__DIR__, "../../../src/viz/ekp_plots.jl"))
using JLD2
import Dates

s = ArgParseSettings()
@add_arg_table s begin
    "--config"
    help = "Inverse problem config file path"
    arg_type = String
end
parsed_args = parse_args(ARGS, s)
include(parsed_args["config"])
config = get_config()

# Initialize calibration process
outdir_path = init_calibration(config; config_path = parsed_args["config"], mode = "pmap")
priors = deserialize_prior(load(joinpath(outdir_path, "prior.jld2")))

# Dispatch SCM eval functions to workers
@everywhere begin
    scm_eval_train(version) = versioned_model_eval(version, $outdir_path, "train", $config)
    scm_eval_validation(version) = versioned_model_eval(version, $outdir_path, "validation", $config)
end

# Calibration process
N_iter = config["process"]["N_iter"]
@info "Running EK updates for $N_iter iterations"
for iteration in 1:N_iter
    @time begin
        @info "   iter = $iteration"
        versions = readlines(joinpath(outdir_path, "versions_$(iteration).txt"))
        ekp = load(ekobj_path(outdir_path, iteration))["ekp"]
        pmap(scm_eval_train, versions)
        pmap(scm_eval_validation, versions)
        ek_update(ekp, priors, iteration, config, versions, outdir_path)
    end
end
@info "Calibration completed. $(Dates.now())"


text:
export ft_setchars

function ft_setchars(ft_handle::Culong, eventch::UInt8, eventchen::Bool,
                                        errorch::UInt8, errorchen::Bool)
  if eventchen
    eventcharacterenable = 0x01
  else
    eventcharacterenable = 0x00
  end
  if errorchen
    errorcharacterenable = 0x01
  else
    errorcharacterenable = 0x00
  end
  ft_status = ccall((:FT_SetChars, d2xx),
                     Cuint,
                     (Culong, Cuchar, Cuchar, Cuchar, Cuchar),
                     ft_handle, eventch, eventcharacterenable, errorch, errorcharacterenable)
  checkstatus(ft_status)
  return nothing
end


text:
##
using GLMakie
using Random

data = rand(10,50000) # 10 channels, 50000 samples
nchan = size(data,1)
nsamp = size(data,2)
srate = 1000
#x_zoom = slider(LinRange(0, 100, 100), raw = true, camera = campixel!, start = 10)
#y_zoom = slider(LinRange(0, 100, 100), raw = true, camera = campixel!, start = 10)
#x_pos = slider(LinRange(0, nsamp, 1000), raw = true, camera = campixel!, start = 0)

##

fig = Figure(resolution = (2560, 1440))
axis = fig[1, 1] = Axis(fig, title = "EEG")

x_zoom = layout[2, 1] = Slider(scene, range = 1:0.1:100, startvalue = 100,tellheight=true);
x_pos  = layout[4, 1] = Slider(scene, range = 0:1:nsamp, startvalue = 1,tellheight = true);
y_zoom = layout[:, 2] = Slider(scene, range = 0:0.1:10, horizontal = false,
         tellwidth = true, height = nothing, width = Auto(),startvalue=3);

##

data_zoom = []
for ch = 1:size(data,1)
    ix_range = lift((from,delta)->range(Int(floor(from+1)),stop=Int(floor(1+from+delta))),x_pos.value,x_zoom.value)
    d = lift((ix,zoom)->ch .+ zoom.*data[ch,ix],ix_range,y_zoom.value)
    append!(data_zoom,[d])
    lines!(axis,@lift($ix_range.*1/srate),data_zoom[ch])
end

#curr_xlim =    xlims!(x_zoom.value,x_zoom.value:)
#ylims!([-y_zoom[end][:value].*y_zoom_fac,y_zoom[end][:value].*y_zoom_fac])

#p = heatmap()
scene
#final = hbox(p, vbox(x_pos,x_zoom,y_zoom))#, parent = Scene(resolution = (1500, 1500)))


text:
using CompScienceMeshes, BEAST
using LinearAlgebra, Pkg

Pkg.activate(@__DIR__)
# Γ = readmesh(joinpath(@__DIR__,"sphere2.in"))
Γ = readmesh(joinpath(dirname(pathof(BEAST)),"../examples/sphere2.in"))
# Γ = meshrectangle(1.0, 1.0, 0.05, 3)
X = lagrangec0d1(Γ)
@show numfunctions(X)

κ = 2π; γ = im*κ
a = -1Helmholtz3D.hypersingular(gamma=γ)
b = Helmholtz3D.doublelayer(gamma=γ) - 0.5Identity()

uⁱ = Helmholtz3D.planewave(wavenumber=κ, direction=ẑ)
f = strace(uⁱ,Γ)
g = ∂n(uⁱ)

@hilbertspace u
@hilbertspace v
eq1 = @discretise a[v,u] == g[v] u∈X v∈X
eq2 = @discretise b[v,u] == f[v] u∈X v∈X

x1 = gmres(eq1)
x2 = gmres(eq2)

# @assert norm(x1-x2)/norm(x1+x2) <  0.5e-2

fcr1, geo1 = facecurrents(x1, X)
fcr2, geo2 = facecurrents(x2, X)

using Plots
plot(title="Comparse 1st and 2nd kind eqs.")
plot!(norm.(fcr1),c=:blue,label="1st")
scatter!(norm.(fcr2),c=:red,label="2nd")

import Plotly
Plotly.plot(patch(Γ, norm.(fcr1)))
Plotly.plot(patch(Γ, norm.(fcr2)))

ys = range(-2,2,length=200)
zs = range(-2,2,length=200)
pts = [point(0.5, y, z) for y in ys, z in zs]

nf = BEAST.HH3DDoubleLayerNear(wavenumber=κ)
near = BEAST.potential(nf, pts, x1, X, type=ComplexF64)
inc = uⁱ.(pts)
tot = near + inc


text:
using Documenter, NbodyGradient

makedocs(sitename="NbodyGradient",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    )
)

deploydocs(
    repo = "github.com/ericagol/NbodyGradient.git"
)


text:
# Lower triangular matrices
module Lower

include("Random.jl")
using .Random

function Cholesky!(L,A,d)
    # Compute L lower triangular such that L*L' = A, given a d-by-d positive definite matrix A.
    # Only L is modified by this procedure.
    for i = 1:d
        s = 0.
        for j = 1:i-1
            r = 0.
            for k = 1:j-1; r += L[i,k]*L[j,k]; end
            L[i,j] = (A[i,j]-r)/L[j,j]
            s += L[i,j]*L[i,j]
        end
        L[i,i] = sqrt(A[i,i]-s)
        for j = i+1:d; L[i,j] = 0.; end
    end
end

function reverse_Cholesky!(L,A,d)
    # Compute L lower triangular such that L'*L = A, given a d-by-d positive definite matrix A.
    # Only L is modified by this procedure.
    for j = d:-1:1
        s = 0.
        for i = d:-1:j+1
            r = 0.
            for k = i+1:d; r += L[k,i]*L[k,j]; end
            L[i,j] = (A[i,j]-r)/L[i,i]
            s += L[i,j]*L[i,j]
        end
        L[j,j] = sqrt(A[j,j]-s)
        for i = 1:j-1; L[i,j] = 0.; end
    end
end

function Cholesky_inverse!(L,M,A,d)
    # Compute L and M lower triangular such that L*L' = inv(A) and M'*M = A,
    # given a d-by-d positive definite matrix A. Only L and M are modified by this procedure.
    reverse_Cholesky!(M,A,d) # M'*M = A
    inverse!(L,M,d) # L = inv(M)
end

function scale!(L,c,d)
    # Multiply L by c, for a scalar c and a d-by-d matrix L.
    for i = 1:d*d; L[i] *= c; end
end

function multiply!(L,M,N,d)
    # Compute L = M*N, for d-by-d lower triangular matrices M and N.
    # NOTE: This is carefully implemented so that L and N can occupy the same memory.
    for i = d:-1:1, j = i:-1:1
        r = 0.
        for k = j:i
            r += M[i,k]*N[k,j]
        end
        L[i,j] = r
    end
end

function multiplyMNt!(L,M,N,d)
    # Compute L = M*N', for d-by-d lower triangular matrices M and N.
    for i = 1:d, j = 1:d
        r = 0.
        for k = 1:min(i,j)
            r += M[i,k]*N[j,k]
        end
        L[i,j] = r
    end
end

function multiplyMtN!(L,M,N,d)
    # Compute L = M'*N, for d-by-d lower triangular matrices M and N.
    # NOTE: This is carefully implemented so that L and N can occupy the same memory.
    for i = 1:d
        for j = 1:d
            r = 0.
            for k = max(i,j):d
                r += M[k,i]*N[k,j]
            end
            L[i,j] = r
        end
    end
end

function logdetsq(L,d)
    # Return log(det(L*L'))=log(det(L'*L))=log(det(L)^2) for a d-by-d triangular matrix L.
    r = 0.
    for i = 1:d; r += log(L[i,i]); end
    return 2.0*r
end

function quadratic(x,y,L,d)
    # Return (x-y)'*L*L'*(x-y) for vectors x,y and a d-by-d lower triangular matrix L.
    r = 0.
    for j = 1:d
        s = 0.
        for i = j:d; s += (x[i]-y[i])*L[i,j]; end
        r += s*s
    end
    return r
end

function sample!(L,M,nu,d)
    # Fill L with a randomly generated lower triangular matrix such that L*L' ~ Wishart(M*M',nu),
    # that is, L*L' is Wishart with scale matrix V=M*M' and nu degrees of freedom.
    # L and M are d-by-d lower triangular matrices, and nu is a real number greater than d-1.
    #
    # This uses the Bartlett decomposition.
    # The same approach is used by Matlab's wishrnd function, so it should be valid
    # for any nu > d-1, including noninteger nu.
    @assert(nu>d-1)
    for i = 1:d
        for j = 1:i-1; L[i,j] = randn(); end
        L[i,i] = sqrt(Random.chi_square(nu-i+1))
        for j = i+1:d; L[i,j] = 0.; end
    end
    multiply!(L,M,L,d)
    return L
end

function inverse!(L,M,d)
    # Compute L = inv(M), for an invertible d-by-d lower triangular matrix M.
    # NOTE: This is faster than inv(M) for d up to around 100, but slower after that (as of Dec 18, 2014).
    for i = 1:d
        for j = i+1:d; L[i,j] = 0.; end
        @assert(M[i,i] != 0.)
        L[i,i] = 1.0/M[i,i]
        for j = i-1:-1:1
            r = 0.
            for k = j+1:i
                r += L[i,k]*M[k,j]
            end
            L[i,j] = -r*L[j,j]  # Multiplication is faster than division
        end
    end
end

function solve_Lx_eq_y!(L,x,y,d)
    # Solve L*x = y for x, for a d-by-d lower triangular matrix L. Only x is modified by this procedure.
    # This uses the forward-substitution algorithm, avoiding the need to invert L.
    # The algorithm is valid if x and y occupy the same memory.
    for i = 1:d
        s = 0.
        for j = 1:i-1
            s += L[i,j]*x[j]
        end
        x[i] = (y[i]-s)/L[i,i]
    end
end

function solve_Ltx_eq_y!(L,x,y,d)
    # Solve L'*x = y for x, for a d-by-d lower triangular matrix L. Only x is modified by this procedure.
    # This uses the back-substitution algorithm, avoiding the need to invert L'.
    # The algorithm is valid if x and y occupy the same memory.
    for i = d:-1:1
        s = 0.
        for j = i+1:d
            s += L[j,i]*x[j]
        end
        x[i] = (y[i]-s)/L[i,i]
    end
end

function sample_Normal!(x,m,L,d)
    # Sample x ~ Normal(m,inv(L*L')), where L is an invertible d-by-d lower triangle matrix.
    # This is done by solving L'*(x-m) = z for x, where z ~ N(0,I). Only x is modified by this procedure.
    for i = d:-1:1
        s = 0.
        for j = i+1:d
            s += L[j,i]*x[j]
        end
        x[i] = (randn()-s)/L[i,i]
    end
    for i = 1:d; x[i] += m[i]; end
    return x
end

end # module Lower


text:
#####################################################################################
# Projected API
#####################################################################################
"""
    projected_integrator(ds::DynamicalSystem, projection, complete_state; kwargs...) → integ

Return an integrator that produces iterations of the dynamical system `ds` on a
projected state space. See [Integrator API](@ref) for handling integrators.

The `projection` defines the projected space. If `projection isa AbstractVector{Int}`,
then the projected space is simply the variable indices that `projection` contains.
Otherwise, `projection` can be an arbitrary function that given the state of the
original system, returns the state in the projected space. In this case the projected
space can be equal, or even higher-dimensional than the original.

`complete_state` produces the state for the original system from the projected state.
`complete_state` can always be a function that given the projected state returns the full
state. However, if `projection isa AbstractVector{Int}`, then `complete_state` can
also be a vector that contains the values of the _remaining_ variables of the system,
i.e., those _not_ contained in the projected space. Obviously in this case
the projected space needs to be lower-dimensional than the original.

Notice that `projected_integrator` does not require an invertible projection,
`complete_state` is only used during [`reinit!`](@ref). Of course the internal
integrator operates in the full state space, where the dynamic rule has been defined.
The projection always happens as a last step, e.g., during [`get_state`](@ref).

## Keyword Arguments
* `u0`: initial state
* `diffeq` is a `NamedTuple` (or `Dict`) of keyword arguments propagated into
  `init` of DifferentialEquations.jl.

## Examples
Case 1: project 5-dimensional system to its last two dimensions.
```julia
ds = Systems.lorenz96(5)
projection = [4, 5]
complete_state = [0.0, 0.0, 0.0] # completed state just in the plane of last two dimensions
pinteg = projected_integrator(ds, projection, complete_state)
reinit!(pinteg, [0.2, 0.4])
step!(pinteg)
get_state(pinteg)
```
Case 2: custom projection to general functions of state.
```julia
ds = Systems.lorenz96(5)
projection(u) = [sum(u), sqrt(u[1]^2 + u[2]^2)]
complete_state(y) = repeat(y[1]/5, 5)
pinteg = # same as in above example...
````
"""
function projected_integrator(ds::GeneralizedDynamicalSystem, projection, complete_state;
        u0 = get_state(ds), diffeq = NamedTuple()
	)
    if projection isa AbstractVector{Int}
        @assert all(1 .≤ projection .≤ dimension(ds))
        projection = SVector(projection...)
        y = u0[projection]
    else
        @assert projection(u0) isa AbstractVector
        y = projection(u0)
    end
    if complete_state isa AbstractVector
        @assert projection isa AbstractVector{Int}
        @assert length(complete_state) + length(projection) == dimension(ds)
        remidxs = setdiff(1:dimension(ds), projection)
        @assert !isempty(remidxs)
    else
        @assert length(complete_state(y)) == dimension(ds)
        remidxs = nothing
    end
    integ = integrator(ds, u0; diffeq)
    u = zeros(dimension(ds))
	return ProjectedIntegrator{
        typeof(projection), length(y), typeof(complete_state),
        typeof(remidxs), typeof(integ)}(projection, complete_state, u, remidxs, integ)
end

struct ProjectedIntegrator{P, PD, C, R, I} <: GeneralizedDynamicalSystem
    projection::P
    complete_state::C
    u::Vector{Float64} # dummy variable for a state in full state space
    remidxs::R
	integ::I
end
isdiscretetime(p::ProjectedIntegrator) = isdiscretetime(p.integ)
DelayEmbeddings.dimension(::ProjectedIntegrator{P, PD}) where {P, PD} = PD

integrator(pinteg::ProjectedIntegrator, args...; kwargs...) = pinteg
#integrator(p::ProjectedIntegrator) = p
get_state(pinteg::ProjectedIntegrator{<:Function}) =
    pinteg.projection(get_state(pinteg.integ))
get_state(pinteg::ProjectedIntegrator{<:SVector}) =
    get_state(pinteg.integ)[pinteg.projection]

function SciMLBase.step!(pinteg::ProjectedIntegrator, args...)
	step!(pinteg.integ, args...)
	return
end

function Base.show(io::IO, pinteg::ProjectedIntegrator)
    println(io, "Integrator of a projected system")
    println(io,  rpad(" rule f: ", 14), get_rule_for_print(pinteg.integ))
    println(io,  rpad(" projection: ", 14), pinteg.projection)
    println(io,  rpad(" dimension: ", 14), dimension(pinteg))
    println(io,  rpad(" complete state: ", 14), pinteg.complete_state)
end

function SciMLBase.reinit!(pinteg::ProjectedIntegrator{P, D, <:AbstractVector}, y) where {P, D}
    u = pinteg.u
    u[pinteg.projection] .= y
    u[pinteg.remidxs] .= pinteg.complete_state
    reinit!(pinteg.integ, u)
end

function SciMLBase.reinit!(pinteg::ProjectedIntegrator{P, D, <:Function}, y) where {P, D}
    reinit!(pinteg.integ, pinteg.complete_state(y))
end

current_time(pinteg::ProjectedIntegrator) = current_time(pinteg.integ)
function (pinteg::ProjectedIntegrator{P})(t)  where {P}
    u = pinteg.integ(t)
    if P <: Function
        return pinteg.projection(u)
    elseif P <: SVector
        return u[pinteg.projection]
    end
end


text:
const JD2000 = 2451544.5
const JD1950 = 2433282.5
const MJD = 2400000.5

function julian(dt::DateTime, base::String="2000")
    return julian(year(dt), month(dt), day(dt), hour(dt), minute(dt), second(dt), base)
end

function julian(year::Int, month::Int, day::Int, hour::Int, minute::Int, second::Int, base::String="2000")
    return julian(year, month, day, hour, minute, float(second), base)
end

function julian(year::Int, month::Int, day::Int, hour::Int, minute::Int, second::FloatingPoint, base::String="2000")
    jd = (367.0 * year - floor((7.0 * (year + floor((month + 9.0) / 12.0))) * 0.25) + floor(275.0 * month / 9.0) + day + 1721013.5 + ((second / 60.0 + minute) / 60.0 + hour) / 24.0)
    if base == "2000"
        jd -= JD2000
    elseif base == "1950"
        jd -= JD1950
    elseif lowercase(base) == "mjd"
        jd -= MJD
    elseif lowercase(base) == "jd"
    else
        error("Unknown base epoch.")
    end
    return jd
end

function gregorian(jd::FloatingPoint, base::String="2000")
    if base == "2000"
        jd += JD2000
    elseif base == "1950"
        jd += JD1950
    elseif lowercase(base) == "mjd"
        jd += MJD
    elseif lowercase(base) == "jd"
        #pass
    else
        error("Unknown base epoch.")
    end
    jd += .5
    Z = itrunc(jd)
    F = jd - Z
    if Z < 2299161
        A = Z
    else
        alpha = int((Z - 1867216.25)/36524.25)
        A = Z + 1 + alpha - int(alpha/4)
    end
    B = A + 1524
    C = itrunc((B - 122.1)/365.25)
    D = itrunc(365.25*C)
    E = itrunc((B - D)/30.6001)
    day = ifloor(B - D - itrunc(30.6001*E) + F)
    month = E < 14 ? E-1 : E-13
    year = month > 2 ? C-4716 : C-4715
    t = F*86400
    hmod = mod(t, 3600)
    second = mod(hmod, 60)
    hour = itrunc((t - hmod)/3600)
    minute = itrunc((hmod - second)/60)
    return year, month, day, hour, minute, second
end

function gregorian(jd::Int, base::String="2000")
    return gregorian(float(jd))
end


text:
"""
    add_scalar_constraint(model::MOI.ModelLike,
                          func::MOI.AbstractScalarFunction,
                          set::MOI.AbstractScalarSet;
                          allow_modify_function::Bool=false)

Adds the scalar constraint obtained by moving the constant term in `func` to
the set in `model`. If `allow_modify_function` is `true` then the function
`func`, can be modified.
"""
function add_scalar_constraint end

function add_scalar_constraint(model::MOI.ModelLike, func::MOI.SingleVariable,
                               set::MOI.AbstractScalarSet;
                               allow_modify_function::Bool=false)
    return MOI.add_constraint(model, func, set)
end
function add_scalar_constraint(model::MOI.ModelLike,
                               func::Union{MOI.ScalarAffineFunction{T},
                                           MOI.ScalarQuadraticFunction{T}},
                               set::MOI.AbstractScalarSet;
                               allow_modify_function::Bool=false) where T
    set = shift_constant(set, -func.constant)
    if !allow_modify_function
        func = copy(func)
    end
    func.constant = zero(T)
    return MOI.add_constraint(model, func, set)
end


text:
function disconnect!(window::AbstractScreen, signal)
    disconnect!(to_native(window), signal)
end
window_area(scene, native_window) = not_implemented_for(native_window)
window_open(scene, native_window) = not_implemented_for(native_window)
mouse_buttons(scene, native_window) = not_implemented_for(native_window)
mouse_position(scene, native_window) = not_implemented_for(native_window)
scroll(scene, native_window) = not_implemented_for(native_window)
keyboard_buttons(scene, native_window) = not_implemented_for(native_window)
unicode_input(scene, native_window) = not_implemented_for(native_window)
dropped_files(scene, native_window) = not_implemented_for(native_window)
hasfocus(scene, native_window) = not_implemented_for(native_window)
entered_window(scene, native_window) = not_implemented_for(native_window)

function register_callbacks(scene::Scene, native_window)

    window_area(scene, native_window)
    window_open(scene, native_window)
    mouse_buttons(scene, native_window)
    mouse_position(scene, native_window)
    scroll(scene, native_window)
    keyboard_buttons(scene, native_window)
    unicode_input(scene, native_window)
    dropped_files(scene, native_window)
    hasfocus(scene, native_window)
    entered_window(scene, native_window)

end

button_key(x::Type{T}) where {T} = error("Must be a keyboard or mouse button. Found: $T")
button_key(x::Type{Keyboard.Button}) = :keyboardstate
button_key(x::Type{Mouse.Button}) = :mousebuttonstate
button_key(x::Set{T}) where {T} = button_key(T)
button_key(x::T) where {T} = button_key(T)

"""
    ispressed(scene, buttons)

Returns true if all `buttons` are pressed in the given `scene`. `buttons` can be
a `Vector` or `Tuple` of `Keyboard` buttons (e.g. `Keyboard.a`), `Mouse` buttons
(e.g. `Mouse.left`) and `nothing`.
"""
function ispressed(scene::SceneLike, button::Union{Vector, Tuple})
    all(x-> ispressed(scene, x), button)
end

# TODO this is a bit shady, but maybe a nice api!
# So you can use void whenever you don't care what is pressed
ispressed(scene::SceneLike, ::Nothing) = true

function ispressed(buttons::Set{T}, button::T) where T <: Union{Keyboard.Button, Mouse.Button}
    return button in buttons
end

function ispressed(buttons::Set{T}, button::Set{T}) where T <: Union{Keyboard.Button, Mouse.Button}
    return issubset(button, buttons)
end

"""
    ispressed(scene, button)

Returns true if `button` is pressed in the given `scene`. The `button` can be
a `Keyboard` button (e.g. `Keyboard.a`), a `Mouse` button (e.g. `Mouse.left`)
or `nothing`. In the latter case `true` is always returned.
"""
function ispressed(scene::SceneLike, button)
    buttons = getfield(events(scene), button_key(button))
    ispressed(buttons, button)
end

"""
Picks a mouse position.  Implemented by the backend.
"""
function pick end

"""
    onpick(func, plot)
Calls `func` if one clicks on `plot`.  Implemented by the backend.
"""
function onpick end


text:
function mrp_quaternion_map(mrp)
	n2 = mrp[1] * mrp[1] + mrp[2] * mrp[2] + mrp[3] * mrp[3]
    M = 2/(1+n2)
    return SVector{4}([(1-n2)/(1+n2), M*mrp[1], M*mrp[2], M*mrp[3]])
end

mrp_rotation_matrix(mrp) = quaternion_rotation_matrix(mrp_quaternion_map(mrp))

mrp = rand(3)

norm(MRP(mrp...) - mrp_rotation_matrix(mrp))

@variables mrp_sym[1:3]

MRP(mrp_sym...)

mrp_rotation_matrix(mrp_sym)


text:
module TestUtil

using FactCheck
importall GradientBoost.Util

facts("Util functions") do
  context("err_must_be_overriden throws an error") do
    @fact_throws err_must_be_overriden()
  end

  context("weighted_median works") do
    weights = [1, 1, 1, 1, 1]
    values = [1, 2, 3, 4, 5]
    @fact weighted_median(weights, values) => 3

    weights = [1, 1, 1, 1, 1]
    values = [5, 4, 3, 2, 1]
    @fact weighted_median(weights, values) => 3

    weights = [1, 1, 1, 1, 1]
    values = [4, 5, 2, 3, 1]
    @fact weighted_median(weights, values) => 3

    weights = [1, 3, 1, 1, 1]
    values = [1, 2, 3, 4, 5]
    @fact weighted_median(weights, values) => 2

    weights = [1, 1, 1, 1, 100]
    values = [1, 2, 3, 4, 5]
    @fact weighted_median(weights, values) => 5

    weights = [100, 1, 1, 1, 1]
    values = [1, 2, 3, 4, 5]
    @fact weighted_median(weights, values) => 1

    weights = [0.5, 0.5, 1, 1, 1]
    values = [1, 2, 3, 4, 5]
    @fact weighted_median(weights, values) => 3
  end

  context("holdout returns proportional partitions") do
    n = 10
    right_prop = 0.3
    (left, right) = holdout(n, right_prop)

    @fact size(left, 1) => n - (n * right_prop)
    @fact size(right, 1) => n * right_prop
    @fact intersect(left, right) => isempty
    @fact size(union(left, right), 1) => n
  end
end

end # module


text:
"Previous interpolation type"
abstract type Previous <: InterpMode end

function interp(xi::Vector{T}, yi::VecOrMat{T}, ::Type{Previous};
                copy=true, sorted=false) where T
    _xi, _yi = _prepare_input(xi, yi, copy, sorted)

    function f(x)
        x < _xi[1] && return nothing
        idx = searchsortedlast(_xi, x)
        idx>length(_xi) && (idx = length(_xi))
        return ndims(_yi) == 1 ? _yi[idx] : _yi[:, idx]
    end

    return f
end


text:
module ConstraintTransformsTests

using Celeste.Model: CanonicalParams, ids, NUM_COLOR_COMPONENTS

using Celeste.DeterministicVI.ElboMaximize: elbo_constraints

using Celeste.DeterministicVI.ConstraintTransforms: Constraint, BoxConstraint,
                                    SimplexConstraint, ParameterConstraint, ConstraintBatch,
                                    inv_logit, logit, to_free,
                                    to_free!, to_bound, to_bound!, simplexify, unsimplexify,
                                    enforce, enforce!, allocate_free_params,
                                    TransformDerivatives, differentiate!,
                                    propagate_derivatives!

import Celeste.Transform, Celeste.SensitiveFloats

include(joinpath(dirname(@__FILE__), "SampleData.jl"))

using Base.Test

srand(1)

####################
# Constraint Types #
####################

@test BoxConstraint <: Constraint
@test SimplexConstraint <: Constraint
@test ParameterConstraint <: Constraint

n = 3
lower = 0.01 / n
scale = rand()

box = BoxConstraint(lower, 2*lower, scale)
@test box.lower === lower
@test box.upper === 2*lower
@test box.scale === scale

@test_throws AssertionError BoxConstraint(lower, lower/2, scale)
@test_throws AssertionError BoxConstraint(-Inf, lower, scale)
@test_throws AssertionError BoxConstraint(lower, Inf, scale)
@test_throws AssertionError BoxConstraint(lower, 2*lower, 0.0)
@test_throws AssertionError BoxConstraint(lower, 2*lower, -scale)
@test_throws AssertionError BoxConstraint(lower, 2*lower, Inf)

simplex = SimplexConstraint(lower, scale, n)
@test simplex.lower === lower
@test simplex.scale === scale
@test simplex.n === n

@test_throws AssertionError SimplexConstraint(lower, scale, 0)
@test_throws AssertionError SimplexConstraint(lower, scale, 1)
@test_throws AssertionError SimplexConstraint(-lower, scale, 1)
@test_throws AssertionError SimplexConstraint(2*lower, scale, 1)
@test_throws AssertionError SimplexConstraint(lower, 0.0, 2)
@test_throws AssertionError SimplexConstraint(lower, 0.0, n)
@test_throws AssertionError SimplexConstraint(lower, -scale, n)
@test_throws AssertionError SimplexConstraint(lower, Inf, n)

c = ParameterConstraint(box, n)
@test c.constraint === box
@test c.indices === n:n

c = ParameterConstraint(box, 1:n)
@test c.constraint === box
@test c.indices === 1:n

c = ParameterConstraint(simplex, n)
@test c.constraint === simplex
@test c.indices === n:n

c = ParameterConstraint(simplex, 1:n)
@test c.constraint === simplex
@test c.indices === 1:n

#################################################
# `to_bound`/`to_bound!` & `to_free`/`to_free!` #
#################################################

function valid_boxed_parameters!(bound, indices, c::BoxConstraint)
    population = linspace(c.lower, c.upper)
    for src in eachindex(bound), i in indices
        bound[src][i] = rand(population)
    end
    return bound
end

function valid_simplexed_parameters!(bound, indices, c::SimplexConstraint)
    @assert length(indices) == c.n
    for src in eachindex(bound)
        params = view(bound[src], indices)
        for i in eachindex(params)
            params[i] = min(rand() + nextfloat(c.lower), 1.0)
        end
        normalize!(params, 1)
        @assert all(c.lower .< params .< 1.0)
        @assert isapprox(sum(params), 1.0)
    end
    return bound
end

# BoxConstraint #
#---------------#

x = rand()
@test inv_logit(x) === -log(1.0 / x - 1)
@test logit(x) === 1.0 / (1.0 + exp(-x))
@test isapprox(logit(inv_logit(x)), x)

lower, scale = rand(2)

c = BoxConstraint(lower, 2*lower, scale)
bound = 1.5 * c.lower

@test_throws AssertionError to_free(c.lower, c)
@test_throws AssertionError to_free(c.upper, c)
@test to_free(bound, c) === inv_logit((bound - c.lower) / (c.upper - c.lower)) * c.scale

free = to_free(bound, c)

@test_throws AssertionError to_bound(-Inf, c)
@test_throws AssertionError to_bound(Inf, c)
@test to_bound(free, c) === bound

# SimplexConstraint #
#-------------------#

x = rand()
c = SimplexConstraint(0.01/5, rand(), 5)

@test simplexify(x, c) === (1 - c.n * c.lower) * x + c.lower
@test unsimplexify(x, c) === (x - c.lower) / (1 - c.n * c.lower)
@test unsimplexify(simplexify(x, c), c) === x

function correct_to_free(bound, bound_range, c::SimplexConstraint)
    log_last = log(unsimplexify(bound[bound_range[end]], c))
    free = similar(bound, length(bound_range) - 1)
    for i in 1:length(free)
        free[i] = c.scale * (log(unsimplexify(bound[bound_range[i]], c)) - log_last)
    end
    return free
end

n_params = 10
n_sources = 3
bound_range = 1:5
free_range = 3:6
bound = collect(rand(n_params) for _ in 1:n_sources)
valid_simplexed_parameters!(bound, bound_range, c)
old_bound = deepcopy(bound)
free = collect(zeros(n_params - 1) for _ in 1:n_sources)

for src in 1:n_sources
    to_free!(free[src], free_range, bound[src], bound_range, c)
    @test free[src][free_range] == correct_to_free(bound[src], bound_range, c)
    for i in setdiff(1:n_sources, src), j in setdiff(1:length(free[src]), free_range)
        @test free[i][j] === 0.0
    end

    bound[src][bound_range] = 0.0
    to_bound!(bound[src], bound_range, free[src], free_range, c)
    @test isapprox(old_bound[src], bound[src])
    @test all(old_bound[i] == bound[i] for i in setdiff(1:n_sources, src))

    copy!(bound[src], old_bound[src])
    fill!(free[src], 0.0)
end

# ParameterConstraint #
#---------------------#

# BoxConstraint

n_params = 10
n_sources = 3
lower, scale = rand(2)
c = ParameterConstraint(BoxConstraint(lower, 2*lower, scale), 3:6)
bound = collect(rand(n_params) for _ in 1:n_sources)
valid_boxed_parameters!(bound, c.indices, c.constraint)
old_bound = deepcopy(bound)
free = collect(zeros(n_params - 1) for _ in 1:n_sources)

for src in 1:n_sources
    free_index = to_free!(free[src], bound[src], c, 1)
    @test free_index == (length(c.indices) + 1)
    for i in 1:length(c.indices)
        @test free[src][i] == to_free(bound[src][c.indices[i]], c.constraint)
    end

    free_index = to_bound!(bound[src], free[src], c, 1)
    @test free_index == (length(c.indices) + 1)
    @test isapprox(bound[src], old_bound[src])
    @test all(old_bound[i] == bound[i] for i in setdiff(1:n_sources, src))

    copy!(bound[src], old_bound[src])
    fill!(free[src], 0.0)
end

# SimplexConstraint

n_params = 10
n_sources = 3
c = ParameterConstraint(SimplexConstraint(0.01/5, rand(), 5), 3:7)
bound = collect(rand(n_params) for _ in 1:n_sources)
valid_simplexed_parameters!(bound, c.indices, c.constraint)
old_bound = deepcopy(bound)
free = collect(zeros(n_params - 1) for _ in 1:n_sources)

for src in 1:n_sources
    free_index = to_free!(free[src], bound[src], c, 1)
    @test free_index == length(c.indices)
    @test free[src][1:(free_index-1)] == correct_to_free(bound[src], c.indices, c.constraint)

    free_index = to_bound!(bound[src], free[src], c, 1)
    @test free_index == length(c.indices)
    @test isapprox(bound[src], old_bound[src])
    @test all(old_bound[i] == bound[i] for i in setdiff(1:n_sources, src))

    copy!(bound[src], old_bound[src])
    fill!(free[src], 0.0)
end

########################
# `enforce`/`enforce!` #
########################

# BoxConstraint #
#---------------#

lower, scale = rand(2)
c = BoxConstraint(lower, 2*lower, scale)

@test enforce(c.lower, c) === nextfloat(c.lower)
@test enforce(c.upper, c) === prevfloat(c.upper)
@test enforce(-Inf, c) === nextfloat(c.lower)
@test enforce(Inf, c) === prevfloat(c.upper)
@test enforce(1.5 * c.lower, c) === 1.5 * c.lower

# SimplexConstraint #
#-------------------#

n = 10
lower = 0.01 / n
c = SimplexConstraint(lower, rand(), n)
bound = valid_simplexed_parameters!([rand(n)], 1:n, c)[1]
old_bound = copy(bound)

enforce!(bound, c)
@test bound == old_bound

bound[2] = 1.0
bound[4] = c.lower
bound[6] = Inf
bound[7] = c.lower * 1.5

enforce!(bound, c)
@test all(c.lower .< bound .< 1.0)
@test isapprox(sum(bound), 1.0)

new_bound = copy(old_bound)
new_bound[2] = prevfloat(1.0)
new_bound[4] = nextfloat(c.lower)
new_bound[6] = prevfloat(1.0)
new_bound[7] = c.lower * 1.5
rescale = (1 - c.n * c.lower) / (sum(new_bound) - c.n * c.lower)
new_bound[:] = [nextfloat(c.lower) + rescale * (i - c.lower) for i in new_bound]
@test bound == new_bound

# ParameterConstraint #
#---------------------#

# BoxConstraint

n_params = 10
n_sources = 3
lower, scale = rand(2)
c = BoxConstraint(lower, 2*lower, scale)
pc = ParameterConstraint(c, 3:7)
bound = collect(rand(n_params) for _ in 1:n_sources)
valid_boxed_parameters!(bound, pc.indices, c)
old_bound = deepcopy(bound)

for src in 1:n_sources
    enforce!(bound[src], pc)

    for i in 1:n_sources
        @test old_bound[i] == bound[i]
        bound[i][3] = c.lower
        bound[i][4] = c.upper
        bound[i][5] = -Inf
        bound[i][6] = Inf
        bound[i][7] = 1.5 * c.lower
    end

    enforce!(bound[src], pc)

    @test bound[src][3] === nextfloat(c.lower)
    @test bound[src][4] === prevfloat(c.upper)
    @test bound[src][5] === nextfloat(c.lower)
    @test bound[src][6] === prevfloat(c.upper)
    @test bound[src][7] === 1.5 * c.lower

    for other_src in setdiff(1:n_sources, src)
        @test bound[other_src][3] === c.lower
        @test bound[other_src][4] === c.upper
        @test bound[other_src][5] === -Inf
        @test bound[other_src][6] === Inf
        @test bound[other_src][7] === 1.5 * c.lower
    end

    for i in 1:n_sources
        for other_param in (1, 2, 8, 9, 10)
            @test bound[i][other_param] == old_bound[i][other_param]
        end
        copy!(bound[i], old_bound[i])
    end
end

# SimplexConstraint

n_params = 10
n_sources = 3
n = 5
lower = 0.01 / n
c = SimplexConstraint(lower, rand(), n)
pc = ParameterConstraint(c, 3:7)
bound = collect(rand(n_params) for _ in 1:n_sources)
valid_simplexed_parameters!(bound, pc.indices, c)
old_bound = deepcopy(bound)

for src in 1:n_sources
    enforce!(bound[src], pc)

    for i in 1:n_sources
        @test old_bound[i] == bound[i]
        bound[i][3] = 1.0
        bound[i][4] = c.lower
        bound[i][6] = Inf
        bound[i][7] = c.lower * 1.5
    end

    new_bound = deepcopy(bound)
    enforce!(bound[src], pc)

    @test all(c.lower .< bound[src][pc.indices] .< 1.0)
    @test isapprox(sum(bound[src][pc.indices]), 1.0)

    new_bound[src][3] = prevfloat(1.0)
    new_bound[src][4] = nextfloat(c.lower)
    new_bound[src][6] = prevfloat(1.0)
    new_bound[src][7] = c.lower * 1.5
    rescale = (1 - c.n * c.lower) / (sum(new_bound[src][pc.indices]) - c.n * c.lower)
    for i in pc.indices
        new_bound[src][i] = nextfloat(c.lower) + rescale * (new_bound[src][i] - c.lower)
    end

    for i in 1:n_sources
        @test new_bound[i] == bound[i]
        copy!(bound[i], old_bound[i])
    end
end

###################
# ConstraintBatch #
###################

loc_width, loc_scale = rand(), rand()
n_sources = 4
bound = collect(rand(length(CanonicalParams)) for _ in 1:n_sources)
constraints = elbo_constraints(bound, loc_width, loc_scale)

for src in 1:n_sources
    u1_box = constraints.boxes[src][1].constraint
    u2_box = constraints.boxes[src][2].constraint
    u1, u2 = bound[src][ids.pos[1]], bound[src][ids.pos[2]]
    @test u1_box === BoxConstraint(u1 - loc_width, u1 + loc_width, loc_scale)
    @test u2_box === BoxConstraint(u2 - loc_width, u2 + loc_width, loc_scale)
end

free = allocate_free_params(bound, constraints)

@test length(free) == length(bound) == n_sources

free_length = length(CanonicalParams) - length(constraints.simplexes[1])

for src in 1:n_sources
    @test length(free[src]) == free_length
end

enforce!(bound, constraints)
old_bound = deepcopy(bound)
to_free!(free, bound, constraints)
old_free = deepcopy(free)
to_bound!(bound, free, constraints)

@test all(isapprox(old_bound[src], bound[src]) for src in 1:n_sources)

# This test currently fails, but that's probably inevitable as long as doing a
# round-trip unconstrain/reconstrain on a parameter set yields parameters that
# are only appromixately equal to the original parameters (which seems like a
# tolerable behavior)
#
# to_free!(free, bound, constraints)
# @test all(isapprox(old_free[src], free[src]) for src in 1:n_sources)

#######################################
# star-only parameter transformations #
#######################################

n_sources = 10

bound = collect(rand(length(CanonicalParams)) for _ in 1:n_sources)

function star_constraints(bound; loc_width = 1.5e-3, loc_scale = 1.0)
    n_sources = length(bound)
    boxes = Vector{Vector{ParameterConstraint{BoxConstraint}}}(n_sources)
    simplexes = Vector{Vector{ParameterConstraint{SimplexConstraint}}}(n_sources)
    for src in 1:n_sources
        i1, i2 = ids.pos[1], ids.pos[2]
        u1, u2 = bound[src][i1], bound[src][i2]
        boxes[src] = [
            ParameterConstraint(BoxConstraint(u1 - loc_width, u1 + loc_width, loc_scale), i1),
            ParameterConstraint(BoxConstraint(u2 - loc_width, u2 + loc_width, loc_scale), i2),
            ParameterConstraint(BoxConstraint(-1.0, 10.0, 1.0), ids.flux_loc[1]),
            ParameterConstraint(BoxConstraint(-10.0, 10.0, 1.0), ids.color_mean[:, 1]),
            ParameterConstraint(BoxConstraint(1e-4, 1.0, 1.0), ids.color_var[:, 1])
        ]
        simplexes[src] = [
            ParameterConstraint(SimplexConstraint(0.005, 1.0, 2), ids.is_star)
            ParameterConstraint(SimplexConstraint(0.01/NUM_COLOR_COMPONENTS, 1.0, NUM_COLOR_COMPONENTS), ids.k[:, 1])
        ]
    end
    return ConstraintBatch(boxes, simplexes)
end

constraints = star_constraints(bound)

free = allocate_free_params(bound, constraints)

@test length(free) == n_sources
@test all(length(free[src]) == 19 for src in 1:n_sources)

enforce!(bound, constraints)
old_bound = deepcopy(bound)
to_free!(free, bound, constraints)
old_free = deepcopy(free)
to_bound!(bound, free, constraints)
@test all(isapprox(old_bound[src], bound[src]) for src in 1:n_sources)

###########################################
# check derivatives against old Transform #
###########################################

# check derivatives #
#-------------------#

ea, vp, catalog = SampleData.gen_n_body_dataset(1)

bound = vp
constraints = elbo_constraints(bound)
free = allocate_free_params(bound, constraints)

enforce!(bound, constraints)
to_free!(free, bound, constraints)

derivs = TransformDerivatives(bound, free)

differentiate!((y, x) -> to_bound!(y, x, constraints, 1), derivs, free[1])

dt = Transform.get_mp_transform(vp, [1])
td = Transform.get_transform_derivatives(vp, [1], dt.bounds)

@test isapprox(td.dparam_dfree, derivs.jacobian)

for i in 1:size(derivs.hessian, 1)
    @test isapprox(td.d2param_dfree2[i], derivs.hessian[i, :, :])
end

# check SensitiveFloat propagation #
#----------------------------------#

ea, vp, catalog = SampleData.gen_n_body_dataset(1)

bound = vp
constraints = elbo_constraints(bound)
free = allocate_free_params(bound, constraints)

enforce!(bound, constraints)
to_free!(free, bound, constraints)

derivs = TransformDerivatives(bound, free)
dt = Transform.get_mp_transform(vp, [1])

sf_bound1 = SensitiveFloats.SensitiveFloat{Float64}(length(bound[1]), 1, true, true)
sf_bound1.v[] = rand()
rand!(sf_bound1.d)
rand!(sf_bound1.h)
sf_bound2 = SensitiveFloats.SensitiveFloat{Float64}(length(bound[1]), 1, true, true)
sf_bound2.v[] = sf_bound1.v[]
copy!(sf_bound2.d, sf_bound1.d)
copy!(sf_bound2.h, sf_bound1.h)
sf_free1 = SensitiveFloats.SensitiveFloat{Float64}(length(free[1]), 1, true, true)
sf_free2 = SensitiveFloats.SensitiveFloat{Float64}(length(free[1]), 1, true, true)

propagate_derivatives!(to_bound!, sf_bound1, sf_free1, free, constraints, derivs)
Transform.transform_sensitive_float!(dt, sf_free2, sf_bound2, vp, [1])

@test isapprox(sf_free1.v[], sf_free2.v[])
@test isapprox(sf_free1.d, sf_free2.d)
@test isapprox(sf_free1.h, sf_free2.h)

end # module


text:
"""
    similarstructure_from_indices

Return the structure of that similar to the input and should always be specified for
custom array/tensor types.
"""
function similarstructure_from_indices end

"""
    memsize(a::Any)

Return the size of the `a` in memory. This is a generic definition, not very efficient.
Provide more efficient version if possible
"""
memsize(a::Any) = Base.summarysize(a)

"""
    similar_from_structure(A, T, structure)

Return an object that is similar to `A` with `eltype = T` and the structure specified by
`structure`.

Generic definitions, should be overwritten if your array/tensor type does not support
`Base.similar(object, eltype, structure)`.
"""
function similar_from_structure(A, T, structure)
    if isbits(T)
        similar(A, T, structure)
    else
        fill!(similar(A, T, structure), zero(T)) # this fixes BigFloat issues
    end
end

"""
    similar_from_indices(T::Type, p1::IndexTuple, p2::IndexTuple, A, CA::Symbol)

Return the object that is similar to `A` which has an `eltype` given by `T` and whose
left indices correspond to the indices `p1` from `op(A)`, and its right indices correspond
to the indices `p2` from `op(A)`, where `op` is `conj` if `CA == :C` or does nothing if
`CA == :N` (default).
"""
function similar_from_indices(T::Type, p1::IndexTuple, p2::IndexTuple, A, CA::Symbol)
    structure = similarstructure_from_indices(T, p1, p2, A, CA)
    similar_from_structure(A, T, structure)
end

"""
    similar_from_indices(T::Type, poA::IndexTuple, poB::IndexTuple,
                p1::IndexTuple, p2::IndexTuple, A, B, CA::Symbol, CB::Symbol)

Return an object that is similar to the tensor map `A`, i.e., the HomSpace from
the domain to the codomain, where the codomain is the tensor product of spaces selected by
`p1` from the list of spaces that are selected by `poA` from `op(A)` and `poB` from `op(B)`
in sequence, and the domain is that by `p2`. Th `opA` is `conj` if `CA == :C` or does
nothing if `CA == :N` (default), and similarly for `opB`.
"""
function similar_from_indices(T::Type, poA::IndexTuple, poB::IndexTuple,
                p1::IndexTuple, p2::IndexTuple, A, B, CA::Symbol, CB::Symbol)
    structure = similarstructure_from_indices(T, poA, poB, p1, p2, A, B, CA, CB)
    similar_from_structure(A, T, structure)
end

"""
    similartype_from_indices(T::Type, p1, p2, A, CA)

Return the type of the object that obtained by `similar_from_indices`. This should work
generically but can be overwritten.
"""
function similartype_from_indices(T::Type, p1, p2, A, CA)
    Core.Compiler.return_type(similar_from_indices,
                                Tuple{Type{T}, typeof(p1), typeof(p2), typeof(A), Symbol})
end

"""
    similartype_from_indices(T::Type, poA, poB, p1, p2, A, B, CA, CB)

Return the type of the object that obtained by `similar_from_indices`. This should work
generically but can be overwritten.
"""
function similartype_from_indices(T::Type, poA, poB, p1, p2, A, B, CA, CB)
    Core.Compiler.return_type(similar_from_indices, Tuple{Type{T}, typeof(poA), typeof(poB),
                typeof(p1), typeof(p2), typeof(A), typeof(B), Symbol, Symbol})
end

"""
    cached_similar_from_indices(sym::Symbol, T::Type, p1::IndexTuple, p2::IndexTuple,
                                    A, CA::Symbol)

The version of `similar_from_indices` that uses cache is `use_cache() == true`. This is
generic, should probably not be overwritten.
"""
function cached_similar_from_indices(sym::Symbol, T::Type, p1::IndexTuple, p2::IndexTuple,
                                        A, CA::Symbol)
    if use_cache()
        structure = similarstructure_from_indices(T, p1, p2, A, CA)
        typ = similartype_from_indices(T, p1, p2, A, CA)
        key = (sym, taskid(), typ, structure)
        C::typ = get!(cache, key) do
            similar_from_indices(T, p1, p2, A, CA)
        end
        return C
    else
        return similar_from_indices(T, p1, p2, A, CA)
    end
end

"""
    cached_similar_from_indices(sym::Symbol, T::Type, poA::IndexTuple, poB::IndexTuple,
                        p1::IndexTuple, p2::IndexTuple, A, B, CA::Symbol, CB::Symbol)

The version of `similar_from_indices` that uses cache is `use_cache() == true`. This is
generic, should probably not be overwritten.
"""
function cached_similar_from_indices(sym::Symbol, T::Type, poA::IndexTuple, poB::IndexTuple,
                        p1::IndexTuple, p2::IndexTuple, A, B, CA::Symbol, CB::Symbol)

    if use_cache()
        structure = similarstructure_from_indices(T, poA, poB, p1, p2, A, B, CA, CB)
        typ = similartype_from_indices(T, poA, poB, p1, p2, A, B, CA, CB)
        key = (sym, taskid(), typ, structure)
        C::typ = get!(cache, key) do
            similar_from_indices(T, poA, poB, p1, p2, A, B, CA, CB)
        end
        return C
    else
        return similar_from_indices(T, poA, poB, p1, p2, A, B, CA, CB)
    end
end


text:
### A Pluto.jl notebook ###
# v0.12.7

using Markdown
using InteractiveUtils

# ╔═╡ 8e1bdc16-2342-11eb-22a6-cb3487fd6818
using LinearAlgebra

# ╔═╡ c810852c-2342-11eb-3c4c-b35d8ca9bb77
using Plots

# ╔═╡ e2bdc99a-2340-11eb-36d1-5b26c095e69d
md"""
# Ma première cellule `Markdown`

Une première équation ``\LaTeX``
```math
\overline{\rho}_i \equiv \frac{n_i}{V}.
```

Pour exécuter une cellule, deux options :

1. Cliquer sur le bouton *play* en bas à droite de la cellule,
1. Presser *Shift* puis *Enter* séquentiellement en maintenant *Shift* appuyé.

"""

# ╔═╡ 7e43d9a4-2341-11eb-2a9b-6f2fab9f6803
function hello(name = "Louis")
	println("Hello, $name !")
end

# ╔═╡ aab3a9b0-2341-11eb-37a1-476a3acf758b
hello("Martin")

# ╔═╡ bccedea6-2341-11eb-3a3c-65363b8441ec
begin
	a = 3
	b = 2
end

# ╔═╡ d52c2852-2341-11eb-1ff2-393759544efe
a, b

# ╔═╡ df4ac668-2341-11eb-0203-2f2059b0c3e9
md"""
Voici grossièrement comment un système de notebook fonctionne.

* Jupyter,
* Pluto.

Le système Pluto.jl est ce qu'on appelle un environnement réactif.

"""

# ╔═╡ 3512729e-2342-11eb-161c-f31631663ff0
c = [1., 2., 4., 10.]

# ╔═╡ 646d1382-2342-11eb-1e6a-c9f8f95c9d0d
typeof(c)

# ╔═╡ 68594ec0-2342-11eb-1b54-8133d0d00cc1
A = [1. 0. 3. 4.
	5. 6. 7. 8.
	1. 0. 3. 4.
	5. 6. 7. 8.]

# ╔═╡ 85dcfd16-2342-11eb-1759-bf7897e56472
typeof(A)

# ╔═╡ 8aed5238-2342-11eb-0bc9-8fbe5e1929fb
A * c

# ╔═╡ 962ebdbc-2342-11eb-1a08-85f6c5fb0acd
eigvals(A)

# ╔═╡ d12ab39e-2342-11eb-181b-895534c67768
plot(exp, lw = 2, xlim = (-2, 3), label = "exp")

# ╔═╡ a2a81b28-2343-11eb-2190-63139c0ce9af
md"""
Le package `DifferentialEquations.jl` permet de résoudre : EDO, DAE, DDE, stochastique

* Problèmes de Cauchy,
* Problèmes aux limites,
* ...

"""

# ╔═╡ Cell order:
# ╟─e2bdc99a-2340-11eb-36d1-5b26c095e69d
# ╠═7e43d9a4-2341-11eb-2a9b-6f2fab9f6803
# ╠═aab3a9b0-2341-11eb-37a1-476a3acf758b
# ╠═bccedea6-2341-11eb-3a3c-65363b8441ec
# ╠═d52c2852-2341-11eb-1ff2-393759544efe
# ╟─df4ac668-2341-11eb-0203-2f2059b0c3e9
# ╠═3512729e-2342-11eb-161c-f31631663ff0
# ╠═646d1382-2342-11eb-1e6a-c9f8f95c9d0d
# ╠═68594ec0-2342-11eb-1b54-8133d0d00cc1
# ╠═85dcfd16-2342-11eb-1759-bf7897e56472
# ╠═8aed5238-2342-11eb-0bc9-8fbe5e1929fb
# ╠═8e1bdc16-2342-11eb-22a6-cb3487fd6818
# ╠═962ebdbc-2342-11eb-1a08-85f6c5fb0acd
# ╠═c810852c-2342-11eb-3c4c-b35d8ca9bb77
# ╠═d12ab39e-2342-11eb-181b-895534c67768
# ╟─a2a81b28-2343-11eb-2190-63139c0ce9af


text:
using LinearAlgebra
using CVChannel
using Test
using DelimitedFiles

"""
In this script we look at the PPT relaxation of the communication value of the
Werner-Holevo channel for multiple copies and its relation to the parameters
λ and the number of copies.
"""

@testset "PPTcv in dimension" begin
    d_list = append!(append!([3],[10:10:100;]),[100:25:750;])
    n_list = [1:4;]
    λ_list = [0,0.25,0.5]
    println("Here we are going to generate data about the PPT relaxation of cv.")
    println("Our interest is largely in how it scales in d.")
    println("We scan over for n ≤ 4, λ in ", λ_list," and d in ", d_list)
    data_table = zeros(length(n_list),length(d_list),length(λ_list))
    λ_ctr = 1
    for λ_val in λ_list
        println("Now evaluating for λ = ", λ_val)
        for n_val in n_list
            println("Now evaluating for n = ", n_val)
            d_ctr = 1
            for d_val in d_list
                A,B,g,a = generalWHLPConstraints(n_val,d_val,λ_val*ones(n_val))
                cv_ppt, opt_var = wernerHolevoCVPPT(n_val,d_val,A,B,g,a)
                data_table[n_val,d_ctr,λ_ctr] = cv_ppt
                d_ctr += 1
            end
        end
        λ_ctr +=1
    end

    header = hcat(["n\\dim:"], d_list')
    data_table_0 = vcat(header,hcat(n_list,data_table[:,:,1]))
    data_table_1 = vcat(header,hcat(n_list,data_table[:,:,2]))
    data_table_2 = vcat(header,hcat(n_list,data_table[:,:,3]))
    println("Here is a susbset of the results.")
    println("Note that each converges to a value as d increases and that n=4 is numerically unstable.")
    println("\nλ=", λ_list[1],":")
    show(stdout, "text/plain", data_table_0[:,1:12])
    println("\nλ=", λ_list[2],":")
    show(stdout, "text/plain", data_table_1[:,1:3:end])
    println("\nλ=", λ_list[3],":")
    show(stdout, "text/plain", data_table_2[:,1:3:end])

    ell = size(data_table_0)[2]
    data_to_save = vcat(data_table_0,zeros(ell)')
    data_to_save = vcat(data_to_save,vcat(data_table_1,zeros(ell)'))
    data_to_save = vcat(data_to_save,data_table_2)
    info_vec = Vector{Union{Nothing,String}}(nothing, size(data_to_save)[1])
    info_vec[1] = "INFO:"
    info_vec[2] = "Generated by multi-copy-werner-holevo-dimension-scaling.jl"
    info_vec[3] = "Tables descend as lambda increases"
    info_vec[4] = "d_list" * string(d_list)
    info_vec[5] = "n_list" * string(n_list)
    info_vec[6] = "lambda_list" * string(λ_list)
    data_to_save = hcat(data_to_save,info_vec)

    print("\nPlease name the file you'd like to write the results to: \n")
    file_name = readline()
    file_to_open = string(file_name,".csv")
    writedlm(file_to_open, data_to_save, ',')

    print("\nWe verify that the values converge for all lambda for")
    print(" n<4, because for n=4, the LP doesn't behave well (see printed data). \n")
    @test isapprox(data_table_0[2:4,38],data_table_0[2:4,39],atol=1e-2)
    @test isapprox(data_table_1[2:4,38],data_table_1[2:4,39],atol=1e-2)
    @test isapprox(data_table_2[2:4,38],data_table_2[2:4,39],atol=1e-2)
end

max_forks_repo_name: ChitambarLab/cv-channel


text:
export verify_host

"""
    verify_host(url::AbstractString, [transport::AbstractString]) :: Bool

The `verify_host` function tells the caller whether the identity of a host
should be verified when communicating over secure transports like TLS or SSH.
The `url` argument may be:

1. a proper URL staring with `proto://`
2. an `ssh`-style bare host name or host name prefixed with `user@`
3. an `scp`-style host as above, followed by `:` and a path location

In each case the host name part is parsed out and the decision about whether to
verify or not is made based solely on the host name, not anything else about the
input URL. In particular, the protocol of the URL does not matter (more below).

The `transport` argument indicates the kind of transport that the query is
about. The currently known values are `SSL` (alias `TLS`) and `SSH`. If the
transport is ommitted, the query will return `true` only if the host name should
not be verified regardless of transport.

The host name is matched against the host patterns in the relavent environment
variables depending on whether `transport` is supplied and what its value is:

- `JULIA_NO_VERIFY_HOSTS` — hosts that should not be verified for any transport
- `JULIA_SSL_NO_VERIFY_HOSTS` — hosts that should not be verified for SSL/TLS
- `JULIA_SSH_NO_VERIFY_HOSTS` — hosts that should not be verified for SSH
- `JULIA_ALWAYS_VERIFY_HOSTS` — hosts that should always be verified

The values of each of these variables is a comma-separated list of host name
patterns with the following syntax — each pattern is split on `.` into parts and
each part must one of:

1. A literal domain name component consisting of one or more ASCII letter,
   digit, hyphen or underscore (technically not part of a legal host name, but
   sometimes used). A literal domain name component matches only itself.
2. A `**`, which matches zero or more domain name components.
3. A `*`, which match any one domain name component.

When matching a host name against a pattern list in one of these variables, the
host name is split on `.` into components and that sequence of words is matched
against the pattern: a literal pattern matches exactly one host name component
with that value; a `*` pattern matches exactly one host name component with any
value; a `**` pattern matches any number of host name components. For example:

- `**` matches any host name
- `**.org` matches any host name in the `.org` top-level domain
- `example.com` matches only the exact host name `example.com`
- `*.example.com` matches `api.example.com` but not `example.com` or
  `v1.api.example.com`
- `**.example.com` matches any domain under `example.com`, including
  `example.com` itself, `api.example.com` and `v1.api.example.com`
"""
function verify_host(
    url :: AbstractString,
    transport :: Union{AbstractString, Nothing} = nothing,
)
    host = url_host(url)
    env_host_pattern_match("JULIA_ALWAYS_VERIFY_HOSTS", host) && return true
    env_host_pattern_match("JULIA_NO_VERIFY_HOSTS", host) && return false
    transport = transport === nothing ? nothing : uppercase(transport)
    return if transport in ("SSL", "TLS")
        !env_host_pattern_match("JULIA_SSL_NO_VERIFY_HOSTS", host)
    elseif transport == "SSH"
        !env_host_pattern_match("JULIA_SSH_NO_VERIFY_HOSTS", host)
    else
        true # do verify
    end
end

function url_host(url::AbstractString)
    m = match(r"^(?:[a-z]+)://(?:[^@/]+@)?([-\w\.]+)"ai, url)
    m !== nothing && return m.captures[1]
    m = match(r"^(?:[-\w\.]+@)?([-\w\.]+)(?:$|:)"a, url)
    m !== nothing && return m.captures[1]
    return nothing # couldn't parse
end

const MATCH_ANY_RE = r""
const MATCH_NONE_RE = r"$.^"

env_host_pattern_match(var::AbstractString, host::AbstractString) =
    occursin(env_host_pattern_regex(var), host)
env_host_pattern_match(var::AbstractString, host::Nothing) =
    env_host_pattern_regex(var) === MATCH_ANY_RE

const ENV_HOST_PATTERN_LOCK = ReentrantLock()
const ENV_HOST_PATTERN_CACHE = Dict{String,Tuple{String,Regex}}()

function env_host_pattern_regex(var::AbstractString)
    lock(ENV_HOST_PATTERN_LOCK) do
        value = get(ENV, var, nothing)
        if value === nothing
            delete!(ENV_HOST_PATTERN_CACHE, var)
            return MATCH_NONE_RE
        end
        old_value, regex = get(ENV_HOST_PATTERN_CACHE, var, (nothing, nothing))
        old_value == value && return regex
        regex = host_pattern_regex(value, var)
        ENV_HOST_PATTERN_CACHE[var] = (value, regex)
        return regex
    end
end

if !@isdefined(contains)
    contains(needle) = haystack -> occursin(needle, haystack)
end

function host_pattern_regex(value::AbstractString, var::AbstractString="")
    match_any = false
    patterns = Vector{String}[]
    for pattern in split(value, r"\s*,\s*", keepempty=false)
        match_any |= pattern == "**"
        parts = split(pattern, '.')
        # emit warning but ignore any pattern we don't recognize;
        # this allows adding syntax without breaking old versions
        if !all(contains(r"^([-\w]+|\*\*?)$"a), parts)
            in = isempty(var) ? "" : " in ENV[$(repr(var))]"
            @warn("bad host pattern$in: $(repr(pattern))")
            continue
        end
        push!(patterns, parts)
    end
    match_any && return MATCH_ANY_RE
    isempty(patterns) && return MATCH_NONE_RE
    regex = ""
    for parts in patterns
        re = ""
        for (i, part) in enumerate(parts)
            re *= if i < length(parts)
                part == "*"  ? "[-\\w]+\\." :
                part == "**" ? "(?:[-\\w]+\\.)*" : "$part\\."
            else
                part == "*"  ? "[-\\w]+" :
                part == "**" ? "(?:[-\\w]+\\.)*[-\\w]+" : part
            end
        end
        regex = isempty(regex) ? re : "$regex|$re"
    end
    return Regex("^(?:$regex)\$", "ai")
end