text:
module SSNVariability
using LinearAlgebra, Statistics
using Distributions

using OrdinaryDiffEq, StochasticDiffEq
import DiffEqNoiseProcess: CorrelatedWienerProcess

"""
    random_covariance_matrix(dims::Integer,diag_val::Real,k_dims::Integer=5)

Returns a random covariance matrix that is positive definite
and has off-diagonal elements.
# Arguments
- `d`: dimensions
- `diag_val`: scaling of the diagonal
- `k-dims`: to regulate off-diagonal elements
"""
function random_covariance_matrix(dims::Integer,diag_val::Real,k_dims::Integer=5)
  W = randn(dims,k_dims)
  S = W*W'+ Diagonal(rand(dims))
  temp_diag = Diagonal(inv.(sqrt.(diag(S))))
  S = temp_diag * S * temp_diag
  S .*= diag_val
  ret = Symmetric(S)
  return ret
end

abstract type  IOFunction{R} end
Base.Broadcast.broadcastable(x::IOFunction)=Ref(x)

function iofun!(dest::Vector{R},source::Vector{R},
    io::IOFunction{R}) where R
  @inbounds @simd for i in eachindex(dest)
    dest[i] = io(source[i])
  end
  return dest
end

function iofun(source::Vector{R},io::IOFunction{R}) where R
  return iofun!(similar(source),source,io)
end

struct IOIdentity{R} <: IOFunction{R}
  function IOIdentity()
    return new{Float64}()
  end
end
function (::IOIdentity{R})(x::R) where R
  return x
end
ioinv(x::R,::IOIdentity{R}) where R<:Real = x
ioprime(x::R,::IOIdentity{R}) where R<:Real = 1.0

struct ReLu{R} <: IOFunction{R}
  α::R
end

function (g::ReLu)(x::Real)
  re = x < 0 ? zero(x) : g.α*x
  return re
end
ioinv(x::R,g::ReLu{R}) where R<:Real = x < zero(R) ? zero(R) : x/g.α
ioprime(x::R,g::ReLu{R}) where R<:Real = x < zero(R) ? zero(R) : g.α

struct ReQuad{R} <:IOFunction{R}
  α::R
end

function (g::ReQuad{R})(x::R) where R
  ret = x < zero(R) ? zero(R) : g.α*x*x
  return ret
end
ioinv(x::R,g::ReQuad{R}) where R =  x < zero(R) ? zero(R) : sqrt(x/g.α)
ioprime(x::R,g::ReQuad{R}) where R = x < zero(R) ? zero(R) : 2.0*g.α*x

struct RecurrentNeuralNetwork{R}
  iofun::IOFunction{R}
  weight_matrix::Matrix{R}
  time_membrane::Vector{R}
  input::Vector{R}
  sigma_noise::Symmetric{R,Matrix{R}}
end
Base.Broadcast.broadcastable(x::RecurrentNeuralNetwork)=Ref(x)
function n_neurons(rnn::RecurrentNeuralNetwork)
  return size(rnn.weight_matrix,1)
end
function Base.ndims(rnn::RecurrentNeuralNetwork)
  return size(rnn.weight_matrix,1)
end

function hasnoise(rn::RecurrentNeuralNetwork{R}) where R
  return tr(rn.sigma_noise) > 100.0*eps(R)
end

# generate weight matrix
"""
        diagtozero!(M::AbstractMatrix{T}) where T
Replaces the matrix diagonal with zeros
"""
function diagtozero!(M::AbstractMatrix{T}) where T
    ms = min(size(M)...)
    zz = zero(T)
    @inbounds @simd for i in 1:ms
        M[i,i] = zz
    end
    return nothing
end

"""
        norm_sum_rows!(mat)
Rescales the matrix by row so that the sum of each row is 1.0
"""
function norm_sum_rows!(mat)
    normf=abs.(inv.(sum(mat,dims=2)))
    return broadcast!(*,mat,mat,normf)
end

const wmat2D_default=[ 2.5  -1.3
                 2.4   -1.  ]

input_base_default(n) =fill(7.,n)
time_membrane_default(ne,ni)=vcat(fill(0.02,ne),fill(0.01,ni))

 # two populations means Dale!
function make_wmat(ne::I,ni::I,w2D::Matrix{<:Real} ; noautapses=true) where I<: Integer
    @assert all(size(w2D) .== 2)
    if ne+ni==2 #nothing to do
        return w2D
    end
    # else, expand the 2D
    Wμ = w2D
    d = Exponential(1.0)
    W=rand(d,ne+ni,ne+ni)
    noautapses && diagtozero!(W)
    # now normalize by block
    pee = 1:ne
    pii = (1:ni) .+ ne
    let wee =  view(W,pee,pee)
        norm_sum_rows!(wee)
        wee .*= Wμ[1,1]
    end
    let wii =  view(W,pii,pii)
        norm_sum_rows!(wii)
        wii .*= Wμ[2,2]
    end
    let wei =  view(W,pee,pii) # inhibitory to excitatory
        norm_sum_rows!(wei)
        wei .*= Wμ[1,2]
    end
    let wie =  view(W,pii,pee) # excitatory to inhibitory
        norm_sum_rows!(wie)
        wie .*= Wμ[2,1]
    end
    return W
end

function make_wmat(ne::I,ni::I; noautapses=true) where I<: Integer
    return make_wmat(ne,ni,wmat2D_default;noautapses=noautapses)
end

# constructor with default values
function RecurrentNeuralNetwork(ne::Integer,ni::Integer;
      sigma_noise::Union{Nothing,AbstractMatrix{<:Real}}=nothing,
      iofunction::IOFunction= ReQuad(0.02))
    wmat = make_wmat(ne,ni)
    ntot=ne+ni
    input=input_base_default(ntot)
    taus=time_membrane_default(ne,ni)
    sigma_noise=something(sigma_noise,Symmetric(zeros(ntot,ntot)))
    return RecurrentNeuralNetwork(iofunction,wmat,taus,input,sigma_noise)
end

function scalebytime!(v::Vector{R},rnn::RecurrentNeuralNetwork{R}) where R
  v ./= rnn.time_membrane
  return v
end
function iofun!(dest::Vector{R},source::Vector{R},
    ntw::RecurrentNeuralNetwork{R}) where R
  return iofun!(dest,source,ntw.iofun)
end
function iofun(x::Real,ntw::RecurrentNeuralNetwork)
  return ntw.iofun(x)
end
function ioinv!(dest::Vector{R},source::Vector{R},
    ntw::RecurrentNeuralNetwork{R}) where R
  for i in eachindex(dest)
    @inbounds dest[i] = ioinv(source[i],ntw.iofun)
  end
  return dest
end
function ioinv(x::Real,ntw::RecurrentNeuralNetwork)
  return ioinv(x,ntw.iofun)
end
function ioprime(x::Real,ntw::RecurrentNeuralNetwork)
  return ioprime(x,ntw.iofun)
end

# Jacobian!

struct JGradPars{R}
    weights::Matrix{R}
    u::Matrix{R}
    ddgu_alloc::Vector{R}
    function JGradPars(ntw::RecurrentNeuralNetwork{R}) where R
        w = similar(ntw.weight_matrix)
        u = similar(w)
        ddgu_alloc = similar(ntw.time_membrane)
        new{R}(w,u,ddgu_alloc)
    end
end

function jacobian!(J,gradpars::Union{Nothing,JGradPars},
            u::Vector{R},dgu::Vector{R},ntw::RecurrentNeuralNetwork{R}) where R
    broadcast!(*,J, ntw.weight_matrix, transpose(dgu)) # multiply columnwise
    J -=I
    #normalize by taus, rowwise
    broadcast!(/,J,J,ntw.time_membrane)
    isnothing(gradpars) && return J
    @error "This part is not ready, yet!"
    # GRADIENTS !  W first
    #gradpars.weights .= gradpars.inv_taus * transpose(dgu)
    #broadcast!(/,gradpars.weights,transpose(dgu),ntw.membrane_taus)
    # now u
    #ddg!(gradpars.ddgu_alloc,u,ntw.gain_function)
    #broadcast!(*,gradpars.u,gradpars.inv_taus,
    #    ntw.weights, transpose(gradpars.ddgu_alloc))
    return J
end

"""
        jacobian(u,rn::RecurrentNetwork)
Jacobian matrix of the network dynamics at point u
"""
@inline function jacobian(u::Vector{R},rn::RecurrentNeuralNetwork{R}) where R
    J = similar(rn.weight_matrix)
    dgu = ioprime.(u,rn)
    return jacobian!(J,nothing,u,dgu,rn::RecurrentNeuralNetwork)
end

function spectral_abscissa(u::Vector{R},rn::RecurrentNeuralNetwork{R}) where R
    J=jacobian(u,rn)
    return maximum(real.(eigvals(J)))
end

### Dynamics

function du_nonoise(u::V,ntw::RecurrentNeuralNetwork{R}) where {R<:Real,V<:Vector{R}}
  return du_nonoise!(similar(u),u,ntw)
end
function du_nonoise!(dest::V,u::V,
        ntw::RecurrentNeuralNetwork{R}) where {R<:Real,V<:Vector{R}}
  r = ntw.iofun.(u)
  dest .= ntw.input .- u
  BLAS.gemv!('N',1.,ntw.weight_matrix,r,1.,dest)
  return scalebytime!(dest,ntw)
end

function run_network_nonoise(ntw::RecurrentNeuralNetwork{R},r_start::Vector{R},
    t_end::Real; verbose::Bool=false,stepsize=0.05) where R<:Real
  u0=ioinv.(r_start,ntw)
  f(du,u,p,t) = du_nonoise!(du,u,ntw)
  prob = ODEProblem(f,u0,(0.,t_end))
  solv = solve(prob,Tsit5();verbose=verbose,saveat=stepsize)
  ret_u = hcat(solv.u...)
  ret_r = iofun.(ret_u,ntw)
  return solv.t,ret_u,ret_r
end

function run_network_noise_simple(ntw::RecurrentNeuralNetwork{R},
    r_start::Vector{R},noiselevel::Real,t_end::Real;
    verbose::Bool=false,stepsize=0.05) where R<:Real
  u0=ioinv.(r_start,ntw)
  f(du,u,p,t) = du_nonoise!(du,u,ntw)
  σ_f(du,u,p,t) = fill!(du,noiselevel)
  prob = SDEProblem(f,σ_f,u0,(0.,t_end))
  solv =  solve(prob,EM();verbose=verbose,dt=stepsize)
  ret_u = hcat(solv.u...)
  ret_r = iofun.(ret_u,ntw)
  return solv.t,ret_u,ret_r
end

function run_network_noise(ntw::RecurrentNeuralNetwork{R},r_start::Vector{R},t_end::Real;
    stepsize=0.05) where R<:Real
  u0=ioinv.(r_start,ntw)
  chol = let c=cholesky(ntw.sigma_noise)
    Matrix(c.L)
  end
  f(du,u,p,t) = du_nonoise!(du,u,ntw)
  g(du,u,p,t) = (@. du = chol)
  prob = SDEProblem(f,g,u0,(0.,t_end);noise_rate_prototype=similar(chol))
  solv = solve(prob,EM();dt=stepsize)
  ret_u = hcat(solv.u...)
  ret_r = iofun.(ret_u,ntw)
  return solv.t,ret_u,ret_r
end

"""
        run_network_to_convergence(u0, rn::RecurrentNeuralNetwork ;
                t_end=80. , veltol=1E-4)
Runs the network as described in [`run_network_nonoise`](@ref), but stops as soon as
`norm(v) / n < veltol` where `v` is the velocity at time `t`.
If this condition is not satisfied (no convergence to attractor), it runs until `t_end` and prints a warning.
# Arguments
- `rn::RecurrentNeuralNetwork`
- `r_start::Vector` : initial conditions
- `t_end::Real`: the maximum time considered
- `veltol::Real` : the norm (divided by num. dimensions) for velocity at convergence
# Outputs
- `u_end`::Vector : the final state at convergence
- `r_end`::Vector : the final state at convergence as rate
"""
function run_network_to_convergence(ntw::RecurrentNeuralNetwork{R},
     r_start::Vector{R};t_end::Real=50. , veltol::Real=1E-1) where R
    function  condition(u,t,integrator)
        v = get_du(integrator)
        return norm(v) / length(v) < veltol
    end
    function affect!(integrator)
        savevalues!(integrator)
        return terminate!(integrator)
    end
    u0=ioinv.(r_start,ntw)
    cb=DiscreteCallback(condition,affect!)
    ode_solver = Tsit5()
    gu_alloc = similar(u0)
    f(du,u,p,t) = du_nonoise!(du,u,ntw)
    prob = ODEProblem(f,u0,(0.,t_end))
    out = solve(prob,Tsit5();verbose=false,callback=cb)
    u_out = out.u[end]
    t_out = out.t[end]
    if isapprox(t_out,t_end; atol=0.05)
        vel = du_nonoise(u_out,ntw)
        @warn "no convergence after max time $t_end"
        @info "the norm (divided by n) of the velocity is $(norm(vel)/length(vel)) "
    end
    return u_out,ntw.iofun.(u_out)
end

################################################
# Here is the algorithm adapted from Hennequin & Lengyel 2016
##################################################

function nu_fun1(mu::Vector{<:Real},sigma_diag::Vector{<:Real},α::Real)
  nr=Normal()
  return  @. α*(
    mu*cdf(nr,mu/sqrt(sigma_diag)) +
    sqrt(sigma_diag)*pdf(nr,mu/sqrt(sigma_diag)))
end
function nu_fun2(mu::Vector{<:Real},sigma_diag::Vector{<:Real},α::Real)
  #nu1=nu_fun1(mu,sigma_diag,α)
  #return @. mu * nu1 + α*sigma_diag*cdf(Normal(),mu/sqrt(sigma_diag))
  nr=Normal()
  return  @. α*(
    mu^2*cdf(nr,mu/sqrt(sigma_diag)) +
    mu*sqrt(sigma_diag)*pdf(nr,mu/sqrt(sigma_diag)) +
    α*sigma_diag*cdf(nr,mu/sqrt(sigma_diag))
    )
end

function nu_fun(mu::Vector{<:Real},sigma_diag::Vector{<:Real},io::ReLu)
  return nu_fun1(mu,sigma_diag,io.α)
end
function nu_fun(mu::Vector{<:Real},sigma_diag::Vector{<:Real},io::ReQuad)
  return nu_fun2(mu,sigma_diag,io.α)
end

function gamma_fun1(mu::Vector{<:Real},sigma_diag::Vector{<:Real},α::Real)
  return @.  α*cdf(Normal(),mu/sqrt(sigma_diag))
end
function gamma_fun2(mu::Vector{<:Real},sigma_diag::Vector{<:Real},α::Real)
  nr=Normal()
  return @. 2*α*(
    mu*cdf(nr,mu/sqrt(sigma_diag)) +
    sqrt(sigma_diag)*pdf(nr,mu/sqrt(sigma_diag)))
end
function gamma_fun(mu::Vector{<:Real},sigma_diag::Vector{<:Real},io::ReLu)
  return gamma_fun1(mu,sigma_diag,io.α)
end
function gamma_fun(mu::Vector{<:Real},sigma_diag::Vector{<:Real},io::ReQuad)
  return gamma_fun2(mu,sigma_diag,io.α)
end

function j_thingy_old(mu::Vector{R},sigma_diag::Vector{R},
    ntw::RecurrentNeuralNetwork{R}) where R
  γ = gamma_fun(mu,sigma_diag,ntw.iofun)
  ret = copy(ntw.weight_matrix)
  broadcast!(*,ret,ret,transpose(γ))
  ret -= I
  return broadcast!(/,ret,ret,ntw.time_membrane)
end
function j_thingy(mu::Vector{R},sigma_diag::Vector{R},
    ntw::RecurrentNeuralNetwork{R}) where R
  T = Diagonal(ntw.time_membrane)
  γ = Diagonal(gamma_fun(mu,sigma_diag,ntw.iofun))
  return T\(ntw.weight_matrix*γ - I)
end

function dmu!(ret::V,mu::V,sigma_diag::V,
    ntw::RecurrentNeuralNetwork{R}) where {R<:Real,V<:Vector{R}}
  ν = nu_fun(mu,sigma_diag,ntw.iofun)
  ret .= ntw.input .- mu
  BLAS.gemv!('N',1.,ntw.weight_matrix,ν,1.,ret)
  return scalebytime!(ret,ntw)
end

function dSigma!(dΣ::M,Σ::M,μ::V,
    ntw::RecurrentNeuralNetwork{R}) where {R<:Real,M<:Matrix{R},V<:Vector{R}}
  J = j_thingy(μ,diag(Σ),ntw)
  dΣ .= ntw.sigma_noise .+ (J*Σ) .+ (Σ*J')
  return dΣ
end

#=
function dmu_step!(ret::V, mu::V,Σ::Matrix{R},δ::Real,
    ntw::RecurrentNeuralNetwork{R}) where {R<:Real,V<:Vector{R}}
  dmu!(ret,mu,diag(Σ),ntw)
  ret .= ret.*δ .+ mu
  return ret
end

function dSigma_step_uniform_noise!(ret::M,Σ::M,μ::V,δ::Real,
    ntw::RecurrentNeuralNetwork{R}) where {R<:Real,V<:Vector{R},M<:Matrix{R}}
  Σdiag=diag(Σ)
  J=j_thingy(μ,Σdiag,ntw)
  lmul!(δ,J)
  J+=I
  copy!(ret,ntw.sigma_noise)
  return BLAS.gemm!('N','T',1.0,J*Σ,J,δ,ret)
end
function dmuSigma_step!(dest_mu::V,dest_sigma::M,μ::V,Σ::M,δ::Real,
    ntw::RecurrentNeuralNetwork{R}) where {R<:Real,V<:Vector{R},M<:Matrix{R}}
  dSigma_step!(dest_sigma,Σ,μ,δ,ntw)
  dmu_step!(dest_mu,μ,dest_sigma,δ,ntw)
  return nothing
end
=#

## initial conditions

function mustart_sigmastart(ntw::RecurrentNeuralNetwork{R}) where R
  mustart = copy(ntw.input)
  if hasnoise(ntw)
    Tm = diagm(-inv.(ntw.time_membrane ))
    sigma_start = lyap(Tm,Matrix(ntw.sigma_noise))
  else
    ntot = ndims(ntw)
    sigma_start = zeros(ntot,ntot)
  end
  return mustart,sigma_start
end

# noise in a linear system
# if S is noise covariance for dx = A x dt + noise
# then the covariance of x, called P , satisfies 0 = A P + P A' + S
function cov_linear(A::Matrix{Float64},cov_noise::AbstractMatrix{Float64})
  return lyap(A,cov_noise)
end

# dx/dt = A x + h + noise
# then  E[x(t)] = - A^-1 * h
function stable_point_linear(A::Matrix{Float64},h::Vector{Float64})
  return -A\h
end

function mean_cov_linear_ntw(ntw::RecurrentNeuralNetwork{R}) where R
  @assert typeof(ntw.iofun) <: IOIdentity
  T = Diagonal(ntw.time_membrane)
  A = T\(ntw.weight_matrix-I)
  h = T\ntw.input
  μ = stable_point_linear(A,h)
  Σ = cov_linear(A,Matrix(ntw.sigma_noise))
  return (μ,Σ)
end

end # module


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

using ACE.SphericalHarmonics: RSHBasis, index_y
using StaticArrays
using JuLIP: AbstractCalculator, Atoms, JVec
using JuLIP.Potentials: MSitePotential, SZList, ZList
using NeighbourLists: neigs

import JuLIP, JuLIP.MLIPs
import JuLIP: cutoff, alloc_temp, alloc_temp_d
import JuLIP.Potentials: evaluate!, evaluate_d!
import JuLIP.MLIPs: combine
import Base: Dict, convert, ==

export RSHIP

struct RSHIP{T, NZ, TJ} <: MSitePotential
   J::TJ
   SH::RSHBasis{T}
   zlist::SZList{NZ}
   # -------------- A and AA datastructures + coefficients
   alists::NTuple{NZ, AList}
   aalists::NTuple{NZ, AAList}
   coeffs::NTuple{NZ, Vector{T}}
end

cutoff(ship::RSHIP) = cutoff(ship.J)

==(S1::RSHIP, S2::RSHIP) =
      all( getfield(S1, i) == getfield(S2, i)
           for i = 1:fieldcount(RSHIP) )

Base.length(ship::RSHIP) = sum(length.(ship.aalists))

# BO + 1 because BO is the number of neighbours not the actual body-order
bodyorder(ship::RSHIP{BO}) where {BO} = BO + 1

# ------------------------------------------------------------
#   FIO code
# ------------------------------------------------------------

Dict(ship::RSHIP{T,NZ}) where {T, NZ} = Dict(
      "__id__" => "ACE_RSHIP",
      "J" => Dict(ship.J),
      "SH_maxL" => ship.SH.maxL,   # TODO: replace this with Dict(SH)
      "T" => string(fltype(ship.SH)),
      "zlist" => Dict(ship.zlist),
      "alists" => [Dict.(ship.alists)...],
      "aalists" => [Dict.(ship.aalists)...],
      "coeffs" => [ ship.coeffs[i] for i = 1:NZ  ]
   )

convert(::Val{:ACE_RSHIP}, D::Dict) = RSHIP(D)

# bodyorder - 1 is because BO is the number of neighbours
# not the actual body-order
function RSHIP(D::Dict)
   T = Meta.eval(Meta.parse(D["T"]))
   J = TransformedJacobi(D["J"])
   SH = SHBasis(D["SH_maxL"], T)
   zlist = read_dict(D["zlist"])
   NZ = length(zlist)
   alists = ntuple(i -> AList(D["alists"][i]), NZ)
   aalists = ntuple(i -> AAList(D["aalists"][i], alists[i]), NZ)
   coeffs = ntuple(i -> T.(D["coeffs"][i]), NZ)
   return  RSHIP(J, SH, zlist, alists, aalists, coeffs)
end

# ------------------------------------------------------------
#   Evaluation code
# ------------------------------------------------------------

alloc_temp(ship::RSHIP{T,NZ}, N::Integer) where {T, NZ} =
   (     J = alloc_B(ship.J),
         Y = alloc_B(ship.SH),
         A = [ zeros(T, length(ship.alists[iz])) for iz=1:NZ ],
      tmpJ = alloc_temp(ship.J),
      tmpY = alloc_temp(ship.SH),
         R = zeros(JVec{T}, N),
         Z = zeros(Int16, N),
         AA = [ zeros(T, length(ship.aalists[iz])) for iz=1:NZ ],
           )

# compute one site energy
function evaluate!(tmp, ship::RSHIP{T},
                   Rs::AbstractVector{JVec{T}},
                   Zs::AbstractVector{<:Integer},
                   z0::Integer) where {T}
   iz0 = z2i(ship, z0)
   precompute_A!(tmp.A[iz0], tmp, ship.alists[iz0], Rs, Zs, ship)
   aalist = ship.aalists[iz0]
   A = tmp.A[iz0]
   c = ship.coeffs[iz0]
   Es = zero(T)
   @inbounds for iAA = 1:length(aalist)
      Esi = c[iAA]
      for α = 1:aalist.len[iAA]
         Esi *= A[aalist.i2Aidx[iAA, α]]
      end
      Es += Esi
   end
   return Es
end

# compute one site energy
# experimental new code. Not very fast unfortunately ...
function evaluate_new!(tmp, ship::RSHIP{T},
                   Rs::AbstractVector{JVec{T}},
                   Zs::AbstractVector{<:Integer},
                   z0::Integer) where {T}
   iz0 = z2i(ship, z0)
   precompute_A!(tmp.A[iz0], tmp, ship.alists[iz0], Rs, Zs, ship)
   aalist = ship.aalists[iz0]
   A = tmp.A[iz0]
   c = ship.coeffs[iz0]
   AA = tmp.AA[iz0]
   copy!(AA, c)
   for α = 1:length(aalist)
      for iAA = 1:length(aalist)
         if α <= aalist.len[iAA]
            AA[iAA] *= A[aalist.i2Aidx[iAA, α]]
         end
      end
   end
   return sum(AA)
end

alloc_temp_d(ship::RSHIP{T, NZ}, N::Integer) where {T, NZ} =
      ( J = alloc_B(ship.J),
       dJ = alloc_dB(ship.J),
        Y = alloc_B(ship.SH),
       dY = alloc_dB(ship.SH),
        A = [ zeros(T, length(ship.alists[iz])) for iz=1:NZ ],
     dAco = [ zeros(T, length(ship.alists[iz])) for iz=1:NZ ],
     tmpJ = alloc_temp_d(ship.J),
     tmpY = alloc_temp_d(ship.SH),
       dV = zeros(JVec{T}, N),
        R = zeros(JVec{T}, N),
        Z = zeros(Int16, N)
      )

function evaluate_d!(dEs, tmp, ship::RSHIP{T},
                     Rs::AbstractVector{JVec{T}},
                     Zs::AbstractVector{<:Integer},
                     z0::Integer
                     ) where {T}
   iz0 = z2i(ship, z0)
   alist = ship.alists[iz0]
   aalist = ship.aalists[iz0]

   # stage 1: precompute all the A values
   precompute_A!(tmp.A[iz0], tmp, alist, Rs, Zs, ship)

   # stage 2: compute the coefficients for the ∇A_{klm} = ∇ϕ_{klm}
   A = tmp.A[iz0]
   dAco = tmp.dAco[iz0]
   _evaluate_d_stage2!(dAco, A, aalist, ship.coeffs[iz0], ship)

   # stage 3: get the gradients
   fill!(dEs, zero(JVec{T}))
   for (iR, (R, Z)) in enumerate(zip(Rs, Zs))
      R̂ = R / norm(R)
      evaluate_d!(tmp.J, tmp.dJ, tmp.tmpJ, ship.J, norm(R))
      evaluate_d!(tmp.Y, tmp.dY, tmp.tmpY, ship.SH, R)
      iz = z2i(ship, Z)
      for iA = alist.firstz[iz]:(alist.firstz[iz+1]-1)
         zklm = alist[iA]
         ik, iy = zklm.k+1, index_y(zklm.l, zklm.m)
         dEs[iR] += dAco[iA] * (
                         tmp.J[ik] * tmp.dY[iy] + tmp.dJ[ik] * tmp.Y[iy] * R̂)
      end
   end
   return dEs
end

function _evaluate_d_stage2!(dAco::AbstractVector{CT}, A, aalist, c, ship::RSHIP
                            ) where {CT}
   fill!(dAco, 0)
   for iAA = 1:length(aalist)
      for α = 1:aalist.len[iAA]
         CxA_α = c[iAA]
         for β = 1:aalist.len[iAA]
            if β != α
               iAβ = aalist.i2Aidx[iAA, β]
               CxA_α *= A[iAβ]
            end
         end
         iAα = aalist.i2Aidx[iAA, α]
         dAco[iAα] += CxA_α
      end
   end
end

ext: jl


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    Coordinates,
    cartesian2d, cartesian3d, internalcoords,
    point, segment,
    x, y, z, u, v, ux, uy, uz, vx, vy, vz,
    vector, distance,
    setcoord!,
    x_components, y_components, z_components,
    radiantophase

abstract type Coordinates end

# TODO: Matrix shape (node-wise or coordinate-wise?)

radiantophase(angle) = mod((angle + 2pi) / 2pi)


text:
export afun

function afun(x::Tx, params::Dict{String, Any}) where {ETx<:Number, Tx<:AbstractArray{ETx}}

if params["logical"] == 1
    xv = vec(x)
    Ind = params["Ind"]
    xv[.~Ind] = zero(ETx)
else
    xv = vec(x)
    Ind = params["Ind"]
    xv[Ind] .= zero(ETx)
end

return xv
end


text:
@doc raw"""
The linear trend model is defined by:
```math
\begin{gather*}
    \begin{aligned}
        y_{t} &=  \mu_{t} + \gamma_{t} + \varepsilon_{t} \quad &\varepsilon_{t} \sim \mathcal{N}(0, \sigma^2_{\varepsilon})\\
        \mu_{t+1} &= \mu_{t} + \nu_{t} + \xi_{t} \quad &\xi_{t} \sim \mathcal{N}(0, \sigma^2_{\xi})\\
        \nu_{t+1} &= \nu_{t} + \zeta_{t} \quad &\zeta_{t} \sim \mathcal{N}(0, \sigma^2_{\zeta})\\
    \end{aligned}
\end{gather*}
```

# Example
```jldoctest
julia> model = LocalLinearTrend(rand(100))
LocalLinearTrend model
```

See more on [Finland road traffic fatalities](@ref)

# References
 * Durbin, James, & Siem Jan Koopman. (2012). "Time Series Analysis by State Space Methods:
    Second Edition." Oxford University Press. pp. 44
"""
mutable struct LocalLinearTrend <: StateSpaceModel
    hyperparameters::HyperParameters
    system::LinearUnivariateTimeInvariant
    results::Results

    function LocalLinearTrend(y::Vector{Fl}) where Fl
        Z = Fl.([1.0; 0.0])
        T = Fl.([1 1; 0 1])
        R = Fl.([1 0; 0 1])
        d = zero(Fl)
        c = zeros(Fl, 2)
        H = one(Fl)
        Q = Fl.([1 0; 0 1])

        system = LinearUnivariateTimeInvariant{Fl}(y, Z, T, R, d, c, H, Q)

        names = ["sigma2_ε", "sigma2_ξ", "sigma2_ζ"]
        hyperparameters = HyperParameters{Fl}(names)

        return new(hyperparameters, system, Results{Fl}())
    end
end

function default_filter(model::LocalLinearTrend)
    Fl = typeof_model_elements(model)
    steadystate_tol = Fl(1e-5)
    a1 = zeros(Fl, 2)
    P1 = Fl(1e6) .* Matrix{Fl}(I, 2, 2)
    return UnivariateKalmanFilter(a1, P1, 2, steadystate_tol)
end

function initial_hyperparameters!(model::LocalLinearTrend)
    Fl = typeof_model_elements(model)
    initial_hyperparameters = Dict{String,Fl}(
        "sigma2_ε" => var(model.system.y),
        "sigma2_ξ" => var(model.system.y),
        "sigma2_ζ" => one(Fl),
    )
    set_initial_hyperparameters!(model, initial_hyperparameters)
    return model
end

function constrain_hyperparameters!(model::LocalLinearTrend)
    constrain_variance!(model, "sigma2_ε")
    constrain_variance!(model, "sigma2_ξ")
    constrain_variance!(model, "sigma2_ζ")
    return model
end

function unconstrain_hyperparameters!(model::LocalLinearTrend)
    unconstrain_variance!(model, "sigma2_ε")
    unconstrain_variance!(model, "sigma2_ξ")
    unconstrain_variance!(model, "sigma2_ζ")
    return model
end

function fill_model_system!(model::LocalLinearTrend)
    model.system.H = get_constrained_value(model, "sigma2_ε")
    model.system.Q[1] = get_constrained_value(model, "sigma2_ξ")
    model.system.Q[end] = get_constrained_value(model, "sigma2_ζ")
    return model
end

function reinstantiate(::LocalLinearTrend, y::Vector{Fl}) where Fl
    return LocalLinearTrend(y)
end

has_exogenous(::LocalLinearTrend) = false


text:
__precompile__()

module Colors

using FixedPointNumbers, ColorTypes, Reexport, Compat
@reexport using ColorTypes
# deprecated exports
export U8, U16

@compat AbstractGray{T} = Color{T,1}
using ColorTypes: TransparentGray
@compat AbstractAGray{C<:AbstractGray,T} = AlphaColor{C,T,2}
@compat AbstractGrayA{C<:AbstractGray,T} = ColorAlpha{C,T,2}
@compat Color3{T} = Color{T,3}
@compat Transparent4{C<:Color3,T} = TransparentColor{C,T,4}

import Base: ==, +, -, *, /
import Base: convert, eltype, hex, isless, linspace, show, typemin, typemax

# Additional exports, not exported by ColorTypes
export weighted_color_mean,
       hex, @colorant_str,
       protanopic, deuteranopic, tritanopic,
       distinguishable_colors, whitebalance,
       colordiff, DE_2000, DE_94, DE_JPC79, DE_CMC, DE_BFD, DE_AB, DE_DIN99, DE_DIN99d, DE_DIN99o,
       MSC, sequential_palette, diverging_palette, colormap,
       colormatch, CIE1931_CMF, CIE1964_CMF, CIE1931J_CMF, CIE1931JV_CMF

# Early utilities
include("utilities.jl")

# Include other module components
include("conversions.jl")
include("promotions.jl")
include("algorithms.jl")
include("parse.jl")
include("differences.jl")
include("colormaps.jl")
include("display.jl")
include("colormatch.jl")

end # module


text:
# Data forms are forms applied over vectors. The purpose is twofold: it's a
# slightly more convenient way of generating a bunch of similar forms, and
# secondly, it allows us to exploit some of the nicer features of d3.

# A DataProperty is a leaf node that may only be attached to a DataForm.
type DataProperty{T <: PropertyPrimitive}
    ds::Vector{AbstractArray}
end

primitive{T}(dp::DataProperty{T}) = T

fill(ds::AbstractArray) = DataProperty{Fill}(AbstractArray[ds])
stroke(ds::AbstractArray) = DataProperty{Stroke}(AbstractArray[ds])
svgclass(ds::AbstractArray) = DataProperty{SVGClass}(AbstractArray[ds])
d3embed(code::String) = DataProperty{D3Embed}(AbstractArray[[code]])

type DataForm{T <: FormPrimitive} <: Form
    # A matrix of parameters to T's constructor.
    ds::Vector{AbstractArray}

    # A list of DataProperty objects applied to the generated Forms in parallel.
    dataprops::List{DataProperty}

    # Scalar property applied to every form.
    property::Property

    # shallow copy constructor
    function DataForm(a::DataForm{T})
        new(a.ds, a.dataprops, a.property)
    end

    function DataForm(ds::Vector{AbstractArray})
        new(ds, ListNil{DataProperty}(), empty_property)
    end
end

DataForm{T}(a::DataForm{T}) = DataForm{T}(a)

primitive{T}(df::DataForm{T}) = T
removable(df::DataForm) = false
children(df::DataForm) = []

copy(df::DataForm) = DataForm(df)

function compose(df::DataForm, p::Property)
    df = copy(df)
    df.property = combine(df.property, p)
    df
end

# Parallel compose.
function compose(df::DataForm, dp::DataProperty)
    df = copy(df)
    df.dataprops = cons(dp, df.dataprops)
    df
end

# DataProperties with regular forms.
compose(f::FormTree, dp::DataProperty) = compose(f, convert(PropertySeq, dp))

# DataForm constructors
function ellipse(x::AbstractArray, y::AbstractArray,
                 rx::AbstractArray, ry::AbstractArray)
    DataForm{Ellipse}(AbstractArray[x, y, rx, ry])
end

function circle(x::AbstractArray, y::AbstractArray, r::AbstractArray)
    ellipse(x, y, r, r)
end

# Convert a DataProperty to a sequence of properties for processing by a non-d3
# backend.
function convert(::Type{PropertySeq}, dp::DataProperty)
    T = primitive(dp)
    PropertySeq(T([d[1] for d in dp.ds]...))
end

# Compose data properties with canvases
compose(a::CanvasTree, b::DataProperty) = compose(a, convert(PropertySeq, b))

# Default handling of data forms. Generate scalar forms and feed them into the
# backend.
function draw{T}(backend::Backend, t::NativeTransform, unit_box::BoundingBox,
                 box::NativeBoundingBox, form::DataForm{T})
    # generate properties
    n = max(map(length, form.ds))
    ps = Property[empty_property for _ in 1:n]
    for dataprop in form.dataprops
        T = primitive(dataprop)
        for (i, row) in enumerate(zip({take(cycle(d), n) for d in dataprop.ds}...))
            ps[i] = combine(ps[i], PropertySeq(T(row...)))
        end
    end

    # combine with non-data properties
    for i in 1:n
        ps[i] = combine(ps[i], form.property)
    end

    # generate scalar forms
    for (i, row) in enumerate(zip([take(cycle(d), n) for d in form.ds]...))
        push_property(backend, ps[i])
        draw(backend, t, unit_box, box, T(row...))
        pop_property(backend)
    end
end

# Generate a d3 '.attr(...)' expression for x-coordinate values.
#
# Args:
#   backend: D3 backend instance.
#   ds: Data vector.
#   unit_box: Parent canvases coordinate system.
#   box: Parent canvases absolute bounding box.
#   dataindexes: Dict mapping indexes in the backend's data matrix to an index
#               within the zipped row.
#
# Returns:
#   A string containing a d3.js ".attr" expression.
#
# Modifies:
#   dataindexes
#
function make_d3_x_attr{T, N}(backend::D3, ds::AbstractArray{T, N},
                              unit_box::BoundingBox, box::NativeBoundingBox,
                              dataindexes::Dict{Int, Int})
    mx = box.width.value / unit_box.width.value
    bx = box.x0.value - unit_box.x0.value * box.width.value / unit_box.width.value

    # For singleton values, we compute the transform in place.
    if length(ds) == 1
        x = native_measure(x_measure(ds[1]), unit_box, box, backend)
        svg_fmt_float(x.value + box.x0.value)
    else
        idx = data_idx(backend, ds)
        if !haskey(dataindexes, idx)
            dataindexes[idx] = length(dataindexes)
        end
        rowidx = dataindexes[idx]
        if T <: Real
            @sprintf("function(d) { return d[%d] * %e + %s; }",
                     rowidx, mx, svg_fmt_float(bx))
        else
            @sprintf("function(d) { return d[%d]; }", rowidx)
        end
    end
end

function make_d3_width_attr{T, N}(backend::D3, ds::AbstractArray{T, N},
                                  unit_box::BoundingBox, box::NativeBoundingBox,
                                  dataindexes::Dict{Int, Int})
    mx = box.width.value / unit_box.width.value

    if length(ds) == 1
        x = native_measure(x_measure(ds[1]), unit_box, box, backend)
        svg_fmt_float(x.value)
    else
        idx = data_idx(backend, ds)
        if !haskey(dataindexes, idx)
            dataindexes[idx] = length(dataindexes)
        end
        rowidx = dataindexes[idx]
        if T <: Real
            @sprintf("function(d) { return d[%d] * %e; }",
                     rowidx, mx)
        else
            @sprintf("function(d) { return d[%d]; }", rowidx)
        end
    end
end

# Generate a d3 '.attr(...)' expression for y-coordinate values.
#
# Args:
#   backend: D3 backend instance.
#   ds: Data vector.
#   unit_box: Parent canvases coordinate system.
#   box: Parent canvases absolute bounding box.
#   dataindexes: Dict mapping indexes in the backend's data matrix to an index
#               within the zipped row.
#
# Returns:
#   A string containing a d3.js ".attr" expression.

# Modifies:
#   dataindexes
#
function make_d3_y_attr{T, N}(backend::D3, ds::AbstractArray{T, N},
                              unit_box::BoundingBox, box::NativeBoundingBox,
                              dataindexes::Dict{Int, Int})
    my = box.height.value / unit_box.height.value
    by = box.y0.value - unit_box.y0.value * box.height.value / unit_box.height.value

    # For singleton values, we compute the transform in place.
    if length(ds) == 1
        y = native_measure(x_measure(ds[1]), unit_box, box, backend)
        svg_fmt_float(y.value + box.y0.value)
    else
        idx = data_idx(backend, ds)
        if !haskey(dataindexes, idx)
            dataindexes[idx] = length(dataindexes)
        end
        rowidx = dataindexes[idx]
        if T <: Real
            @sprintf("function(d) { return d[%d] * %e + %s; }",
                     rowidx, my, svg_fmt_float(by))
        else
            @sprintf("function(d) { return d[%d]; }", rowidx)
        end
    end
end

function make_d3_height_attr{T, N}(backend::D3, ds::AbstractArray{T, N},
                                   unit_box::BoundingBox, box::NativeBoundingBox,
                                   dataindexes::Dict{Int, Int})
    my = box.height.value / unit_box.height.value

    # For singleton values, we compute the transform in place.
    if length(ds) == 1
        y = native_measure(x_measure(ds[1]), unit_box, box, backend)
        svg_fmt_float(y.value)
    else
        idx = data_idx(backend, ds)
        if !haskey(dataindexes, idx)
            dataindexes[idx] = length(dataindexes)
        end
        rowidx = dataindexes[idx]
        if T <: Real
            @sprintf("function(d) { return d[%d] * %e; }",
                     rowidx, my)
        else
            @sprintf("function(d) { return d[%d]; }", rowidx)
        end
    end
end

# Build a series of ".attr" expressions representing properties applied to a
# data form.
function make_d3_data_property_expr(backend::D3,
                                    unit_box::BoundingBox,
                                    box::NativeBoundingBox,
                                    dataindexes::Dict{Int, Int},
                                    dps::List{DataProperty})
    expr = IOBuffer()
    for dp in dps
        print(expr, make_d3_data_property_expr(backend, unit_box,
                                               box, dataindexes, dp))
    end
    takebuf_string(expr)
end

function make_d3_data_property_expr(backend::D3,
                                    unit_box::BoundingBox,
                                    box::NativeBoundingBox,
                                    dataindexes::Dict{Int, Int},
                                    dp::DataProperty{Fill})
    colors = dp.ds[1]
    if length(colors) == 1
        @sprintf(".attr(\"fill\", %s)\n", json(color(colors[1])))
    else
        idx = data_idx(backend, colors)
        if !haskey(dataindexes, idx)
            dataindexes[idx] = length(dataindexes)
        end
        rowidx = dataindexes[idx]
        @sprintf(".attr(\"fill\", function(d) { return d[%d]; })\n", rowidx)
    end
end

function make_d3_data_property_expr(backend::D3,
                                    unit_box::BoundingBox,
                                    box::NativeBoundingBox,
                                    dataindexes::Dict{Int, Int},
                                    dp::DataProperty{Stroke})
    colors = dp.ds[1]
    if length(colors) == 1
        @sprintf(".attr(\"stroke\", %s)\n", json(color(colors[1])))
    else
        idx = data_idx(backend, colors)
        if !haskey(dataindexes, idx)
            dataindexes[idx] = length(dataindexes)
        end
        rowidx = dataindexes[idx]
        @sprintf(".attr(\"stroke\", function(d) { return d[%d]; })\n", rowidx)
    end
end

function make_d3_data_property_expr(backend::D3,
                                    unit_box::BoundingBox,
                                    box::NativeBoundingBox,
                                    dataindexes::Dict{Int, Int},
                                    dp::DataProperty{SVGClass})
    classes = dp.ds[1]
    if length(classes) == 1
        @sprintf(".attr(\"class\", \"%s\")\n", escape_string(classes[1]))
    else
        idx = data_idx(backend, classes)
        if !haskey(dataindexes, idx)
            dataindexes[idx] = length(dataindexes)
        end
        rowidx = dataindexes[idx]
        @sprintf(".attr(\"class\", function(d) { return d[%d]; })\n", rowidx)
    end
end

function make_d3_data_property_expr(backend::D3,
                                    unit_box::BoundingBox,
                                    box::NativeBoundingBox,
                                    dataindexes::Dict{Int, Int},
                                    dp::DataProperty{D3Embed})
    string(dp.ds[1][1], "\n")
end

# Build expression suitable for the argument of a '.data(...)' call.
#
# Args:
#   dataindexes: Dict mapping indexes in the backend's data matrix to an index
#               within the zipped row.
#
function make_d3_data_expr(dataindexes::Dict{Int, Int})
    idxpairs = Array(Tuple, length(dataindexes))
    for (i, (k, v)) in enumerate(dataindexes)
        idxpairs[i] = (v, k)
    end
    sort!(idxpairs)
    @sprintf("d3.zip(%s)",
             join(["data[$(idxpair[2])]" for idxpair in idxpairs], ","))
end

function draw(backend::D3, t::NativeTransform, unit_box::BoundingBox,
              box::NativeBoundingBox, form::DataForm{Ellipse})

    class = next_dataform_class(backend)
    x, y, rx, ry = form.ds

    dataindexes = Dict{Int, Int}()
    x_expr = make_d3_x_attr(backend, x, unit_box, box, dataindexes)
    y_expr = make_d3_y_attr(backend, y, unit_box, box, dataindexes)

    dp_expr = make_d3_data_property_expr(backend, unit_box, box,
                                         dataindexes, form.dataprops)
    push_property(backend, form.property)

    if rx == ry
        r_expr = make_d3_width_attr(backend, rx, unit_box, box, dataindexes)
        @printf(backend.out,
                "g.selectAll(\"%s\")
                  .data(%s)
                  .enter()
                  .append(\"circle\")\n",
                  class, make_d3_data_expr(dataindexes))
        @printf(backend.out, ".attr(\"cx\", %s)\n", x_expr)
        @printf(backend.out, ".attr(\"cy\", %s)\n", y_expr)
        @printf(backend.out, ".attr(\"r\", %s)\n", r_expr)
    else
        rx_expr = make_d3_width_attr(backend, rx, unit_box, box, dataindexes)
        ry_expr = make_d3_height_attr(backend, ry, unit_box, box, dataindexes)
        @printf(backend.out,
                "g.selectAll(\"%s\")
                  .data(%s)
                  .enter()
                  .append(\"ellipse\")\n",
                  class, make_d3_data_expr(dataindexes))
        @printf(backend.out, ".attr(\"cx\", %s)\n", x_expr)
        @printf(backend.out, ".attr(\"cy\", %s)\n", y_expr)
        @printf(backend.out, ".attr(\"rx\", %s)\n", rx_expr)
        @printf(backend.out, ".attr(\"ry\", %s)\n", ry_expr)
    end

    print(backend.out, dp_expr)
    print(backend.out, ";\n")
    pop_property(backend)
end

# What we need now: text and lines. Also, serialization of the Aesthetics object
# into JSON.


text:
module Tokenization

importall ..JukaiNLP
using Merlin

export train, Tokenizer

include("tagset.jl")
include("tokenizer.jl")
include("train.jl")

end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
function acqr_ksmoother(a,xtt,Ptt,xtt1,Ptt1)
	#
	# Kalman smoother for model
	#
	# x_{t+1} = ax_{t} + w_{t}
	# y_{t}   = cx_{t} + v_{t}
	#
	# cov(w_{t},v_{t}) = [q 0;0 r]
	#
	# javier.cara@upm.es, 2020-05
	#

	nt = length(xtt)

	# allocation
	xtN = zeros(nt+1)
	PtN = zeros(nt+1)
	Pt1tN = zeros(nt)

	# values for t = nt+1
	xtN[nt+1] = xtt1[nt+1]
	PtN[nt+1] = Ptt1[nt+1]

	# values for t=nt
	xtN[nt] = xtt[nt]
	PtN[nt] = Ptt[nt]
	Pt1tN[nt] = a*Ptt[nt]

	# smother
	for t = nt-1:-1:1
		# Kalman Smoother matrix J
		Jt = (Ptt[t]*a) / Ptt1[t+1]

		xtN[t] = xtt[t] + Jt*( xtN[t+1] - xtt1[t+1] )
		PtN[t] = Ptt[t] + Jt^2*( PtN[t+1] - Ptt1[t+1] )

		Pt1tN[t] = PtN[t+1]*Jt
	end

	return xtN,PtN,Pt1tN

end


text:
using DiffEqJump, DiffEqBase
using Test, HypergeometricFunctions

Nsims = 1e4

# ABC model A + B <--> C
reactstoch = [
    [1 => 1, 2 => 1],
    [3 => 1],
]
netstoch = [
    [1 => -1, 2 => -1, 3 => 1],
    [1 => 1, 2 => 1, 3 => -1]
]
rates = [0.1, 1.]
u0 = [500,500,0]
tspan = (0.0,5.0)
prob = DiscreteProblem([500,500,0],(0.0,2.0), rates)
majumps = MassActionJump(rates, reactstoch, netstoch)

function getmean(jprob,Nsims)
    Amean = 0
    for i = 1:Nsims
        sol = solve(jprob,SSAStepper())
        Amean += sol[1,end]
    end
    Amean /= Nsims
    Amean
end

K = 1/.1
function analyticmean(u, K)
    α = u[1]; β = u[2]; γ = u[3]
    @assert β ≥ α "A(0) must not exceed B(0)"
    K * (α+γ)/(β-α+1) * pFq([-α-γ+1], [β-α+2], -K) / pFq([-α-γ], [β-α+1], -K)
end
analytic_mean = analyticmean(u0, K)

algs = DiffEqJump.JUMP_AGGREGATORS
relative_tolerance = 0.01
for alg in algs
    local jprob = JumpProblem(prob,alg,majumps,save_positions=(false,false))
    local Amean = getmean(jprob, Nsims)
    @test abs(Amean - analytic_mean)/analytic_mean < relative_tolerance
end


text:
# Create a Tile object with position, number of mine neighbours, and states describing
# if tile: is flagged, is shown and if is a mine.
mutable struct Tile
    shown::Bool
    is_mine::Bool
    x::Int
    y::Int
    num_mines::Int
    is_flagged::Bool
end

# Return the string to show on screen based on current state of cell
function to_string(tile)
    # If tile is shown, return the corresponding character
    if tile.shown
        if tile.is_mine == true
            return "•"
        else
            return tile.num_mines
        end

    # If tile is not shown but is flagged, show a flag. Else, show a box with a question mark.
    elseif tile.is_flagged
        return "!"
    else
        return "�"
    end
end

# Count the number of mine neighbours for a tile.
function count_neighbours(board, tile)
    # If the tile itself is a mine, let the number of mine neighbours be -1
    if (tile.is_mine)
        tile.num_mines = -1
        return -1
    else
        # Iterate over neighbours in grid, and count number of mine neighbours
        total = 0
        for x_offset in -1:1
            for y_offset in -1:1
                i = tile.x + x_offset
                j = tile.y + y_offset
                if ((i in 1:length(board)) && (j in 1:length(board[1])) && (board[i][j].is_mine))
                    total += 1
                end
            end
        end

        # Return the number obtained
        tile.num_mines = total
        return total
    end
end

max_stars_repo_name: soumitradev/Minesweeper

max_issues_repo_name: soumitradev/Minesweeper

max_forks_repo_name: soumitradev/Minesweeper


text:
"""
```julia
gradient(O::Expression, vars::AbstractVector{<:Expression}; simplify = true)
```

A helper function for computing the gradient of an expression with respect to
an array of variable expressions.
"""
function gradient(O::Expression, vars::AbstractVector{<:Expression}; simplify = true)
    [expand_derivatives(Differential(v)(O),simplify) for v in vars]
end

"""
```julia
jacobian(ops::AbstractVector{<:Expression}, vars::AbstractVector{<:Expression}; simplify = true)
```

A helper function for computing the Jacobian of an array of expressions with respect to
an array of variable expressions.
"""
function jacobian(ops::AbstractVector{<:Expression}, vars::AbstractVector{<:Expression}; simplify = true)
    [expand_derivatives(Differential(v)(O),simplify) for O in ops, v in vars]
end

"""
```julia
hessian(O::Expression, vars::AbstractVector{<:Expression}; simplify = true)
```

A helper function for computing the Hessian of an expression with respect to
an array of variable expressions.
"""
function hessian(O::Expression, vars::AbstractVector{<:Expression}; simplify = true)
    [expand_derivatives(Differential(v2)(Differential(v1)(O)),simplify) for v1 in vars, v2 in vars]
end

function simplified_expr(O::Operation)
  if O isa Constant
    return O.value
  elseif isa(O.op, Differential)
    return :(derivative($(simplified_expr(O.args[1])),$(simplified_expr(O.op.x))))
  elseif isa(O.op, Variable)
    isempty(O.args) && return O.op.name
    return Expr(:call, Symbol(O.op), simplified_expr.(O.args)...)
  end
  return Expr(:call, Symbol(O.op), simplified_expr.(O.args)...)
end

simplified_expr(c::Constant) = c.value

function simplified_expr(eq::Equation)
    Expr(:(=), simplified_expr(eq.lhs), simplified_expr(eq.rhs))
end

simplified_expr(eq::AbstractArray) = simplified_expr.(eq)


text:
# test InertialPose3

using RoME
# , Distributions
using JLD

@load joinpath(dirname(@__FILE__),"preintstationarydata.jld")
# DATA[1]

global N = 100
global fg = initfg()

global initCov = 0.001*Matrix{Float64}(LinearAlgebra.I, 15,15)
[initCov[i,i] = 0.0001^2 for i in 4:6];
[initCov[i,i] = 0.0002^2 for i in 10:15];
global odoCov = deepcopy(initCov)

println("Adding PriorInertialPose3 to graph...")
global v1 = addVariable!(fg, :x1, InertialPose3, N=N) #0.1*randn(15,N)
global initPosePrior = PriorInertialPose3( MvNormal( zeros(15), initCov) )
global f1  = addFactor!(fg, [v1], initPosePrior)

global n = 1
global v2 = addVariable!(fg, :x2, InertialPose3, dims=15,  N=N)
global noise = MvNormal(zeros(15),(DATA[n][3]+DATA[n][3]')*0.5 )
global inerodo = InertialPose3(noise,DATA[n][1],DATA[n][2])
global f2  = addFactor!(fg, [v1;v2], inerodo )

initVariable!(fg, :x2, N=N)

# drawDensityMesh(fg, :x2)

plotKDE(fg, :x2, dims=[1;2]  , title="x,y")
plotKDE(fg, :x2, dims=[2;3]  , title="y,z")
plotKDE(fg, :x2, dims=[6;3]  , title="ψ,z")
plotKDE(fg, :x2, dims=[4;5]  , title=",")
plotKDE(fg, :x2, dims=[7;8]  , title=",")
plotKDE(fg, :x2, dims=[9;10] , title=",")
plotKDE(fg, :x2, dims=[11;12], title=",")
plotKDE(fg, :x2, dims=[13;14], title=",")
plotKDE(fg, :x2, dims=[15]   , title=",")

DATA[n][3][1,1]

global res = zeros(15)
global idx = 1
# pos, so3, vel, bw, ba
global meas = vectoarr2([0;0;9.81/2; 0;0;0; 0;0;9.81; 0;0;0; 0;0;0])
global wIPi = zeros(15,1)
global wIPj = zeros(15,1)

inerodo(res, nothing, idx, (meas,), wIPi, wIPj)

using Optim

ggo = (x) -> inerodo(res, nothing, idx, (meas,), wIPi, vectoarr2(x))

ggos = (x) -> ggo([x[1:5]...,0.0,x[6:end]...])

global ret = optimize(ggo, zeros(15))

using Gadfly

global obj = (x) -> ggo([0,0.0,0, 0,0,x, 0,0,0, 0,0,0, 0,0,0])
plot(obj, -2, 2)

global ran = range(-2,stop=2,length=100)
ggoxy = (x,y) -> ggo([0,0,0.0, 0,0,x, y,0,0, 0,0,0, 0,0,0])
plot(z=ggoxy, x=ran, y=ran, Geom.contour)

@show res
using NLsolve

gg = (res, x) -> inerodo(res, nothing, idx, (meas,), wIPi, vectoarr2(x))

global ret = nlsolve(gg, wIPj[:])

@show ret.zero

using KernelDensityEstimate
p = kde!( getSample(inerodo, 100)[1] )
plotKDE(p, dims=[1;2]  )
plotKDE(p, dims=[6;3]  )
plotKDE(p, dims=[4;5]  )
plotKDE(p, dims=[7;8]  )
plotKDE(p, dims=[9;10] )
plotKDE(p, dims=[11;12])
plotKDE(p, dims=[13;14])
plotKDE(p, dims=[15]   )

# Juno.breakpoint("/home/dehann/.julia/v0.5/IncrementalInference/src/ApproxConv.jl", 22)

using Gadfly

spy(DATA[n][3])

norm(DATA[n][3]-DATA[n][3]')

ensureAllInitialized!(fg)

tree, smt, hist = solveTree!(fg)
# global tree = wipeBuildNewTree!(fg)
# inferOverTree!(fg, tree)

plotKDE(fg, :x1, dims=[4;5])

writeGraphPdf(fg)
run(`evince fg.pdf`)
Base.rm("fg.pdf")

#


text:
export kfn

import ..KFNModel

using mlpack._Internal.io

import mlpack_jll
const kfnLibrary = mlpack_jll.libmlpack_julia_kfn

# Call the C binding of the mlpack kfn binding.
function kfn_mlpackMain()
  success = ccall((:kfn, kfnLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module kfn_internal
  import ..kfnLibrary

import ...KFNModel

# Get the value of a model pointer parameter of type KFNModel.
function IOGetParamKFNModel(paramName::String)::KFNModel
  KFNModel(ccall((:IO_GetParamKFNModelPtr, kfnLibrary), Ptr{Nothing}, (Cstring,), paramName))
end

# Set the value of a model pointer parameter of type KFNModel.
function IOSetParamKFNModel(paramName::String, model::KFNModel)
  ccall((:IO_SetParamKFNModelPtr, kfnLibrary), Nothing, (Cstring, Ptr{Nothing}), paramName, model.ptr)
end

# Serialize a model to the given stream.
function serializeKFNModel(stream::IO, model::KFNModel)
  buf_len = UInt[0]
  buf_ptr = ccall((:SerializeKFNModelPtr, kfnLibrary), Ptr{UInt8}, (Ptr{Nothing}, Ptr{UInt}), model.ptr, Base.pointer(buf_len))
  buf = Base.unsafe_wrap(Vector{UInt8}, buf_ptr, buf_len[1]; own=true)
  write(stream, buf)
end
# Deserialize a model from the given stream.
function deserializeKFNModel(stream::IO)::KFNModel
  buffer = read(stream)
  KFNModel(ccall((:DeserializeKFNModelPtr, kfnLibrary), Ptr{Nothing}, (Ptr{UInt8}, UInt), Base.pointer(buffer), length(buffer)))
end
end # module

"""
    kfn(; [algorithm, epsilon, input_model, k, leaf_size, percentage, query, random_basis, reference, seed, tree_type, true_distances, true_neighbors, verbose])

This program will calculate the k-furthest-neighbors of a set of points. You may
specify a separate set of reference points and query points, or just a reference
set which will be used as both the reference and query set.

For example, the following will calculate the 5 furthest neighbors of eachpoint
in `input` and store the distances in `distances` and the neighbors in
`neighbors`:

```julia
julia> using CSV
julia> input = CSV.read("input.csv")
julia> distances, neighbors, _ = kfn(k=5, reference=input)
```

The output files are organized such that row i and column j in the neighbors
output matrix corresponds to the index of the point in the reference set which
is the j'th furthest neighbor from the point in the query set with index i.  Row
i and column j in the distances output file corresponds to the distance between
those two points.

# Arguments

 - `algorithm::String`: Type of neighbor search: 'naive', 'single_tree',
      'dual_tree', 'greedy'.  Default value `dual_tree`.

 - `epsilon::Float64`: If specified, will do approximate furthest neighbor
      search with given relative error. Must be in the range [0,1).  Default
      value `0`.

 - `input_model::KFNModel`: Pre-trained kFN model.
 - `k::Int`: Number of furthest neighbors to find.  Default value `0`.

 - `leaf_size::Int`: Leaf size for tree building (used for kd-trees, vp
      trees, random projection trees, UB trees, R trees, R* trees, X trees,
      Hilbert R trees, R+ trees, R++ trees, and octrees).  Default value `20`.

 - `percentage::Float64`: If specified, will do approximate furthest
      neighbor search. Must be in the range (0,1] (decimal form). Resultant
      neighbors will be at least (p*100) % of the distance as the true furthest
      neighbor.  Default value `1`.

 - `query::Array{Float64, 2}`: Matrix containing query points (optional).
 - `random_basis::Bool`: Before tree-building, project the data onto a
      random orthogonal basis.  Default value `false`.

 - `reference::Array{Float64, 2}`: Matrix containing the reference
      dataset.
 - `seed::Int`: Random seed (if 0, std::time(NULL) is used).  Default
      value `0`.

 - `tree_type::String`: Type of tree to use: 'kd', 'vp', 'rp', 'max-rp',
      'ub', 'cover', 'r', 'r-star', 'x', 'ball', 'hilbert-r', 'r-plus',
      'r-plus-plus', 'oct'.  Default value `kd`.

 - `true_distances::Array{Float64, 2}`: Matrix of true distances to
      compute the effective error (average relative error) (it is printed when
      -v is specified).
 - `true_neighbors::Array{Int, 2}`: Matrix of true neighbors to compute
      the recall (it is printed when -v is specified).
 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `distances::Array{Float64, 2}`: Matrix to output distances into.
 - `neighbors::Array{Int, 2}`: Matrix to output neighbors into.
 - `output_model::KFNModel`: If specified, the kFN model will be output
      here.

"""
function kfn(;
             algorithm::Union{String, Missing} = missing,
             epsilon::Union{Float64, Missing} = missing,
             input_model::Union{KFNModel, Missing} = missing,
             k::Union{Int, Missing} = missing,
             leaf_size::Union{Int, Missing} = missing,
             percentage::Union{Float64, Missing} = missing,
             query = missing,
             random_basis::Union{Bool, Missing} = missing,
             reference = missing,
             seed::Union{Int, Missing} = missing,
             tree_type::Union{String, Missing} = missing,
             true_distances = missing,
             true_neighbors = missing,
             verbose::Union{Bool, Missing} = missing,
             points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, kfnLibrary), Nothing, ());

  IORestoreSettings("k-Furthest-Neighbors Search")

  # Process each input argument before calling mlpackMain().
  if !ismissing(algorithm)
    IOSetParam("algorithm", convert(String, algorithm))
  end
  if !ismissing(epsilon)
    IOSetParam("epsilon", convert(Float64, epsilon))
  end
  if !ismissing(input_model)
    kfn_internal.IOSetParamKFNModel("input_model", convert(KFNModel, input_model))
  end
  if !ismissing(k)
    IOSetParam("k", convert(Int, k))
  end
  if !ismissing(leaf_size)
    IOSetParam("leaf_size", convert(Int, leaf_size))
  end
  if !ismissing(percentage)
    IOSetParam("percentage", convert(Float64, percentage))
  end
  if !ismissing(query)
    IOSetParamMat("query", query, points_are_rows)
  end
  if !ismissing(random_basis)
    IOSetParam("random_basis", convert(Bool, random_basis))
  end
  if !ismissing(reference)
    IOSetParamMat("reference", reference, points_are_rows)
  end
  if !ismissing(seed)
    IOSetParam("seed", convert(Int, seed))
  end
  if !ismissing(tree_type)
    IOSetParam("tree_type", convert(String, tree_type))
  end
  if !ismissing(true_distances)
    IOSetParamMat("true_distances", true_distances, points_are_rows)
  end
  if !ismissing(true_neighbors)
    IOSetParamUMat("true_neighbors", true_neighbors, points_are_rows)
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("distances")
  IOSetPassed("neighbors")
  IOSetPassed("output_model")
  # Call the program.
  kfn_mlpackMain()

  return IOGetParamMat("distances", points_are_rows),
         IOGetParamUMat("neighbors", points_are_rows),
         kfn_internal.IOGetParamKFNModel("output_model")
end


text:
__precompile__(false)

module ReviseTest

square(x) = x^2

cube(x) = error("cube")

fourth(x) = x^4  # this is an addition to the file

module Internal

mult2(x) = error("mult2")
mult3(x) = 3*x

end  # Internal

end


text:
using BenchmarkTools
using ExtremeFloats
using Quadmath
using ArbNumerics
using MultiFloats

function log_sum_exp(exponents, ::Type{T}) where T
    sum(x -> exp(T(x)), exponents) |> log |> Float64
end

n = 1000
exponents = rand(n) .* 10.0 .^ rand(-10:10, n)

@info "Running log-sum-exp benchmarks..."
for T in (Float64, BigFloat, ArbFloat, Quadmath.Float128, MultiFloats.Float64x8, ExtremeFloat)
    @info T
    try
        result = @btime log_sum_exp($exponents, $T)
        @show result
    catch ex
        @info "Calculation failed."
    end
    println()
end


text:
module TestTraits

using MLJBase
using Test
using MLJIteration
using ..DummyModel

model = DummyIterativeModel()
imodel = IteratedModel(model=model, measure=mae)

@test is_wrapper(imodel)
@test !MLJBase.caches_data_by_default(imodel)
@test !supports_weights(imodel)
@test !supports_class_weights(imodel)
@test load_path(imodel) == "MLJIteration.DeterministicIteratedModel"
@test package_name(imodel) == "MLJIteration"
@test package_uuid(imodel) == "614be32b-d00c-4edb-bd02-1eb411ab5e55"
@test package_url(imodel) == "https://github.com/JuliaAI/MLJIteration.jl"
@test package_license(imodel) == "MIT"
@test is_pure_julia(imodel) == is_pure_julia(model)
@test input_scitype(imodel) == input_scitype(model)
@test output_scitype(imodel) == output_scitype(model)
@test target_scitype(imodel) == target_scitype(model)

end

true


text:
using Documenter, KDESieve

makedocs(;
    modules=[KDESieve],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/rjww/KDESieve.jl/blob/{commit}{path}#L{line}",
    sitename="KDESieve.jl",
    authors="Robert Woods",
    assets=String[],
)

deploydocs(;
    repo="github.com/rjww/KDESieve.jl",
)


text:
struct Gaussian{T}
    scaling::T
    width::T
    delay::T
end

Gaussian(;scaling=1.0, width, delay) = Gaussian(scaling, width, delay)

(g::Gaussian)(s::Real) = 4*g.scaling/(g.width*√π) * exp(-(4*(s-g.delay)/g.width)^2)

function creategaussian(width,s0,scaling=one(typeof(width)))
    #f(s) = 4*scaling/(width*sqrt(π)) * exp(-(4*(s-s0)/width)^2)
    Gaussian(scaling, width, s0)
    #f(s) = scaling * exp(-(4*(s-s0)/width)^2)
end

function  fouriertransform(g::Gaussian)
    scaling = g.scaling
    width = g.width
    s0 = g.delay
    ft(w) = scaling * exp(-im*w*s0 - (width*w/8)^2) / sqrt(2π)
end

function fouriertransform(a::Array, dt, t0, dim=1)
    n = size(a,dim)
    dω = 2π / (n*dt)
    b = fftshift(fft(a, dim), dim) * dt / sqrt(2π)
    ω0 = -dω * div(n,2)
    b, dω, ω0
end

fouriertransform(a::Array; stepsize, offset, dim=1) = fouriertransform(a, stepsize, offset, dim)

derive(g::Gaussian) =  s -> g(s) * (-8 * (s-g.delay)/g.width) * (4/g.width)

struct ErrorFunction{T}
    scaling::T
    width::T
    delay::T
end

function (f::ErrorFunction)(s)
    f.scaling * 0.5 * (1 + erf(4*(s-f.delay)/f.width))
end

function integrate(f::Gaussian)
    return ErrorFunction(f.scaling, f.width, f.delay)
end


text:
#md ---
#md layout: post
#md title: "kinetic models in real life"
#md date: 2019-01-10 00:00:00 +0000
#md categories: blog
#md mathjax: true
#md ---

# # kinetic models in real life

# >In this post I'll introduce **computer based kinetic models** and how to apply them to real-life problems.

# >To give an intuitive hint about *kinetic models* and their applications we might think about gossip spreading. There are more relevant examples, e.g. the epidemiologic study of virus diffusion, but for the moment I'd rather keep a general and *nice* example. Classial statistics woulg give an overall interpretation suggesting that 1% of the people are gossip aware. Classical statistics, at least in its simpler form, cannot identify who knows the gossip and who does not know. On the contrary, kinetic models try to fill this gap by identifying, exactly, who is gossip aware.

# >While statistics tell you that there will be a Martin Luther King every 100 years, kinetic modelling is trying to tell you who's gonna be.

# Instead of giving a formal mathematical definition, we might want to uvail this concept by constructing together an example. Instead of referring to some epidemiologic model, e.g. a hillness speading, despite pertinet sounds less dramatic and sillier to image how a given gossip can spread from person to person.
# Intuitively what we are planning to do is rather simple. One person, in a room, knows a gossip and meeting with a colleague share such a piece of gossip. Now te two people moving in room will eventually meet with other two colleagues and the gossip will be shared among four people. And clearly, so on.
# The room example is intuitive enought and allow for smooth programming without many technicalities, still it misses the most of the model: the gossip being broght outside of the room to a pub, to a different building, to people above some age since the too youngs cannot be interested, and so on.
# It is rather important that while the kinetic effects are important, we still want to measure some fluid (or aggregated) quantities, such as the total number of people aware of the gossip and
# the radius of influence of the gossip. How far is the gossip spreading? These models can answer this question.

# let's place 100 people in a room,
using Plots
n_people=1000
x, y, flag=rand(n_people), rand(n_people), zeros(n_people)
people = zeros(n_people,3)
flag[1]=1
[people[i,:]=[x[i],y[i],flag[i]] for i=1:n_people]
# and gather the information in a matrix *people* with 3 columns. The first column contains the person position along the *x* direction, the second along the *y* direction and the third a flag that track if the person is a ware of the gossip *flag=1* or not *flag=0*.
# We want to set up this function to plot people position,
function plot_(people)
    rectangle(x, y, w, h) = Shape(x .+ [0,w,w,0], y .+ [0,0,h,h])
    plot(rectangle(0,0,1,1), fill="white", xlim=(-.2,1.2), ylim=(-.2,1.2), background_color=RGB(0.2,0.2,0.2), label="room")
    scatter!( people[ findall(x->x==0., people[:,3]) ,1], people[ findall(x->x==0., people[:,3]) ,2], color="rgb(234, 153, 153)", shape=:circle, ms=8, alpha=1.0, label="not aware")
    scatter!( people[ findall(x->x==1., people[:,3]) ,1], people[ findall(x->x==1., people[:,3]) ,2], color="rgb(119, 221, 119)", shape=:square, ms=8, alpha=1.0, label="aware")
end
plot_(people)
savefig(".\\_img\\002_kinetic_models_in_real_life\\002_kinetic_models_in_real_life_01.svg")
# ![dots and squares](IMAGEFOLDER/002_kinetic_models_in_real_life_01.svg)

# We plot with a circle a person that is not aware of the gossip and with a *square* the persons that are aware.

# We now create a function *move_()* that advance people position.
function move_(people)
    for idx = 1:size(people)[1]
        people[idx,:]= [ people[idx,1]+rand()*0.01, people[idx,2]+rand()*0.01, people[idx,3] ]
        if people[idx,1] > 1.
            people[idx,1] = 1-people[idx,1]
        end
        if people[idx,1] < 0.
            people[idx,1] = -people[idx,1]
        end
        if people[idx,2] > 1.
            people[idx,2] = 1-people[idx,2]
        end
        if people[idx,2] < 0.
            people[idx,2] = -people[idx,2]
        end
    end
    return people
end
# Despite a person trajectory depends on what the person is doing, the person is going to a colleague or to the printer, so shoud have a given trajectory depending on the schope we simplify by assuming that
# a person new position only depends on the present position and his advancment is normal distributed both in *x* and *y*. This assumption, goes, alongside with the assumption (for the moment) that points are volumeless and do not hit one another.
# These simplifications are operated to focus our attention on the *collision operator*, when two points are close enought the two people interact and if one of the persons is aware of the gossip shares it.
function collision_(people)
    new_flag = zeros( size(people)[1] )
    for idx = 1:size(people)[1]
        if people[idx,3] == 1.
            distances = ( (people[:,1].-people[idx,1]).^2 .+ (people[:,2].-people[idx,2]).^2 ).^0.5
            distances = (distances .< 0.05)
            new_flag[distances] .= 1
        end
    end
    people[:,3] = new_flag
    return people
end

# Let's adavance the system and let's check its evolution.
for time = 1:10
    move_(people)
    collision_(people)
end
plot_(people)
savefig(".\\_img\\002_kinetic_models_in_real_life\\002_kinetic_models_in_real_life_02.svg")
# ![time evolution](IMAGEFOLDER/002_kinetic_models_in_real_life_02.svg)

# # A brief disquisition on the underlying maths
# What I personally like about kinetic models is the possibility to explain them in simple words with down to heath example.
# I reckon the gossip example is rather simple, eventually funny, and intuitively simple. The underlying math follow a different formalism known as Bolzmann's statistics.
# A person, that we have represented with a circle or eventually a square, can be modelled mathematically as a Dirac Delta centered where the person is,
# ```math
# p(x_i) = \delta(x-x_i)  \:\:  :: i=1,...,N
# ```
# since each person p(x_i) is also moving is also characterised by a velocity, that using the same notation,
# ```math
# p(x_i) = \delta(x-x_i) \: \delta(v-v_i)  \:\:   :: i=1,...,N
# ```
# and to describe the entire population we sum all over the people,
# ```math
# P = \Sum_{i=0}_N \delta(x-x_i) \delta(v-v_i)
# ```
# The complexity of the problem scales linearly with the number of people *N*, and analytically we quickly encount a calculation obstaacle whether we wish to solve the system evolution by hand.
# To be able to treat analytically such a problem it is convenient to treat the ansemble behaviour of *P*, ansembled over infinitesimal volumes and infinitesimal time-intervals.

# -----------------
cd("C:\\Users\\a.marocchino\\Documents\\codes\\albz.github.io")                                                           #src
using Literate                                                                                #src
# preprocess for notebooks                                                                                #src
function setimagefolder(content)                                                                                #src
    content = replace(content, "IMAGEFOLDER" => "$IMAGEFOLDER")                                                           #src
    return content                                                                                #src
end                                                                                #src
# for Jupyter notebook, put images in subfolder                                                                           #src
IMAGEFOLDER = "_img\\002_kinetic_models_in_real_life"                                                                        #src
Literate.notebook("_source\\002_kinetic_models_in_real_life.jl", "notebooks", preprocess = setimagefolder)                   #src
# for Markdown/Jekyll notebook, put images in "/images"                                                                   #src
IMAGEFOLDER = "/_img/002_kinetic_models_in_real_life"                                                                         #src
Literate.markdown("_source\\002_kinetic_models_in_real_life.jl", ".", name="_posts\\002_kinetic_models_in_real_life",               #src
 preprocess = setimagefolder,                                                                                #src
 documenter=false)                                                                                #src


text:
# Investor's Risk Preference
abstract type RiskPreference end

"""
Risk Neutrality
"""
struct RiskNeutral <: RiskPreference end

"""
This struct contains the parameters for modeling Risk Aversion,
    where a project utility is calculated as follows:
    U = constant - multiplier * e ^ (- risk_coefficient ^ NPV)
"""
struct RiskAverse <: RiskPreference
    constant::Float64
    multiplier::Float64
    risk_coefficient::Float64
end

get_risk_coefficient(preference::RiskAverse) = preference.risk_coefficient
get_constant(preference::RiskAverse) = preference.constant
get_multiplier(preference::RiskAverse) = preference.multiplier


text:
using Flux, Zygote
abstract type CurvatureInterface end

"""
    jacobians(curvature::CurvatureInterface, X::AbstractArray)

Computes the Jacobian `∇f(x;θ)` where `f: ℝᴰ ↦ ℝᴷ`.
"""
function jacobians(curvature::CurvatureInterface, X::AbstractArray)
    nn = curvature.model
    # Output:
    ŷ = nn(X)
    # Jacobian:
    𝐉 = jacobian(() -> nn(X),Flux.params(nn))
    𝐉 = reduce(hcat,[𝐉[θ] for θ ∈ curvature.params])
    return 𝐉, ŷ
end

"""
    gradients(curvature::CurvatureInterface, X::AbstractArray, y::Number)

Compute the gradients with respect to the loss function: `∇ℓ(f(x;θ),y)` where `f: ℝᴰ ↦ ℝᴷ`.
"""
function gradients(curvature::CurvatureInterface, X::AbstractArray, y::Number)
    model = curvature.model
    𝐠 = gradient(() -> curvature.loss(X,y),Flux.params(model))
    return 𝐠
end
struct EmpiricalFisher <: CurvatureInterface
    model::Any
    loss::Function
    params::AbstractArray
end

"""
    full(curvature::EmpiricalFisher, d::Union{Tuple,NamedTuple})

Compute the full empirical Fisher.
"""
function full(curvature::EmpiricalFisher, d::Tuple)
    x, y = d
    𝐠 = gradients(curvature, x, y)
    𝐠 = reduce(vcat,[vec(𝐠[θ]) for θ ∈ curvature.params])

    # Empirical Fisher:
    H = 𝐠 * 𝐠'
    return H

end


text:
is_valid(x, ::Val{:one_hot}) = sum(x) == 1

function binarize(x, d::D, ::Val{:one_hot}) where {T <: Number, D <: DiscreteDomain{T}}
    y = zeros(T, length(d))
    is_in = false
    for (i, v) in enumerate(get_domain(d))
        if x == v
            y[i] = 1
            is_in = true
            break
        end
    end
    @assert is_in "The value $x is not in the domain $d"
    return y
end

function debinarize(x, d::D, ::Val{:one_hot}) where {T <: Number, D <: DiscreteDomain{T}}
    valid = false
    val = typemax(T)
    for (b, v) in Iterators.zip(x, get_domain(d))
        if b == 1
            valid && break
            valid = true
            val = v
        end
    end
    return val
end

function debinarize(x, domains, ::Val{:one_hot})
    start = 0
    stop = 0
    function aux(d)
        start = stop + 1
        stop += length(d)
        return debinarize(@view(x[start:stop]), d, Val(:one_hot))
    end
    return map(aux, domains)
end


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, ["libddierckx"], :libddierckx),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/kbarbary/DierckxBuilder/releases/download/v1.0.0-1"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/libddierckx.v1.0.0.arm-linux-gnueabihf-gcc4.tar.gz", "c83224796f4a8fdbf3d74ba91f2a430680611be3a1e7c8ef1517377dad9bf606"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/libddierckx.v1.0.0.arm-linux-gnueabihf-gcc7.tar.gz", "2246e492f55df8759bf957e8ffa562bfd241206470fe1d7c78c67d060e86e69a"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/libddierckx.v1.0.0.arm-linux-gnueabihf-gcc8.tar.gz", "f096dbacfd47aec9033fac72be9e4a3922d53a275c051c42c053c3e385c03c28"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/libddierckx.v1.0.0.i686-linux-gnu-gcc4.tar.gz", "2a52576cd189121dd584d7d9ff3f00b1a7ed469081217033e4b793e4f29ae85b"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/libddierckx.v1.0.0.i686-linux-gnu-gcc7.tar.gz", "a1859e7cbd1c4d7b6fab82cd583dba8ac0f3562cb049c427c55d96c4d36d6eed"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/libddierckx.v1.0.0.i686-linux-gnu-gcc8.tar.gz", "dfa3de97899cff88dba969c9951fce210f01fb34da89a64d0fd78887e563628a"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/libddierckx.v1.0.0.i686-w64-mingw32-gcc4.tar.gz", "919eceff111d9cbc6431f239bf6f23097656e93937e66ffdcab67a47c172e941"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/libddierckx.v1.0.0.i686-w64-mingw32-gcc7.tar.gz", "38d1dc9f63ad9a8465c6b49ea4336d93d9d92487b84809d77d78717022406e44"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/libddierckx.v1.0.0.i686-w64-mingw32-gcc8.tar.gz", "e0e7fa9534449f6c397efd195de2af54a9e0b010c603adb22119ec6f1f49feef"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/libddierckx.v1.0.0.x86_64-apple-darwin14-gcc7.tar.gz", "ff0e49b74a05012cfc1198da460b0e9d36dca90adf53c3964836f86a2ad7275c"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/libddierckx.v1.0.0.x86_64-apple-darwin14-gcc8.tar.gz", "809a205cd7a8df3d3723e9c013b855c9cb034ff881d3eda820f26809d1efad8a"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/libddierckx.v1.0.0.x86_64-linux-gnu-gcc4.tar.gz", "2b4061d85a1f5b40f37880e27be00b1b3dc7c628d88e1a024c64492aed067339"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/libddierckx.v1.0.0.x86_64-linux-gnu-gcc7.tar.gz", "5b6c3b511ca64e59054a3c196c7e840d3b8d2d4ad04beefb2bd246573f67e936"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/libddierckx.v1.0.0.x86_64-linux-gnu-gcc8.tar.gz", "f75ae37c668c276922254ef6320e2ff0824bae1a88e4f554396a873df7885f29"),
    FreeBSD(:x86_64, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/libddierckx.v1.0.0.x86_64-unknown-freebsd11.1-gcc4.tar.gz", "dcbd52553fe730852ee31c1fc0f2b0f6b5aa6db5f820505aa9ba77b73ebcc4d7"),
    FreeBSD(:x86_64, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/libddierckx.v1.0.0.x86_64-unknown-freebsd11.1-gcc7.tar.gz", "22ed4639567896a469c00dd7e5cd4ac97038f4763d055698c1073a0932c13b32"),
    FreeBSD(:x86_64, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/libddierckx.v1.0.0.x86_64-unknown-freebsd11.1-gcc8.tar.gz", "f550f78201b27c82e41103b6d5db482c266281eb01968a33b98b82b9fb40b171"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/libddierckx.v1.0.0.x86_64-w64-mingw32-gcc4.tar.gz", "d67dcd3521c973e3deb0efecdf192871ba5c22c994426394d474fa6c7ab7a9a7"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/libddierckx.v1.0.0.x86_64-w64-mingw32-gcc7.tar.gz", "98ebdea7468d0fcd0c5ec266766c0cdb4e9f9c5cb89dbd53f5102f646a533385"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/libddierckx.v1.0.0.x86_64-w64-mingw32-gcc8.tar.gz", "eb54c673b95cae7b0743858a789aa770422b823f97bfde2efca16b13800c3818"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
dl_info = choose_download(download_info, platform_key_abi())
if dl_info === nothing && unsatisfied
    # If we don't have a compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\") is not supported by this package!")
end

# If we have a download, and we are unsatisfied (or the version we're
# trying to install is not itself installed) then load it up!
if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
    # Download and install binaries
    install(dl_info...; prefix=prefix, force=true, verbose=verbose)
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)


text:
# a support the Julia-ish interface when it maps cleanly onto a distributed
# vector

export set_values1!, get_values1!, assembly_begin, assembly_end, size_local,
       size_global, length_local, length_global, local_indices, fill_zero!,
       diagonal_shift!

"""
  Create a vector of a given size.  Users can specify either the global
  dimension or the local dimension
"""
function PetscVec(mglobal::Integer, comm::MPI_Comm; mlocal=PETSC_DECIDE)
  vec = PetscVec(comm)
  VecSetSizes(vec::PetscVec, mlocal, mglobal)

  return vec
end

"""
  Create a PetscVec, setting both the type and the format.  Users can specify
  either the local or global dimensions
"""
function PetscVec(mglobal::Integer, format, comm::MPI_Comm; mlocal=PETSC_DECIDE)
  vec = PetscVec(mglobal, comm, mlocal=mlocal)
  VecSetType(vec, format)
  return vec
end

"""
  PetscDestroy for AbstractVector.  No-op
"""
function PetscDestroy(vec::AbstractVector)
end

# 1-based indexing that unifies regular matrices and Petsc matrices
# TODO: generalize to AbstractArray when possible
#----------------------------------------------------------------------------

"""
  1-based indexing for both regular vectors and Petsc vector

  **Inputs**

   * vals: values to add/insert into the vector, must be length(idx)
   * flag: INSERT_VALUES or ADD_VALUES

  **Inputs/Outputs**

  * vec: the vector, can be a Petsc vector or a julia vector
  * idx: (global) indices to add/insert vals into

  idx is listed as input/output because it may be modified during the function.
  It will be returned to its original values when the function exits.
  This is necessary to accomodate Petscs zero-based indexing interface

"""
function set_values1!(vec::PetscVec, idx::Array{PetscInt}, vals::Array{PetscScalar}, flag::Integer=INSERT_VALUES)
  for i=1:length(idx)
    idx[i] -= 1
  end

  err = VecSetValues(vec, idx, vals, flag)

  for i=1:length(idx)
    idx[i] += 1
  end

  return err
end

"""
  Method for AbstractVector
"""
function set_values1!{T}(vec::AbstractVector, idx::Array{PetscInt}, vals::Array{T}, flag::Integer=INSERT_VALUES)
  if flag == INSERT_VALUES
    for i in idx
      vec[i] = vals[i]
    end
  elseif flag == ADD_VALUES
    for i in idx
      vec[i] += vals[i]
    end
  end

  return PetscErrorCode(0)
end

"""
  Like [`set_values1!`](@ref) but for retrieving values.  Note that Petsc
  only supports retrieving values from the local part of the vector

  **Inputs**

   * vec: a vector, can be a julia vector or a Petsc vector.

  **Inputs/Outputs**

   * idx: indices to retrieve
   * vals: array to put the values into

"""
function get_values1!(vec::PetscVec, idx::Array{PetscInt}, vals::Array{PetscScalar})
  for i=1:length(idx)
    idx[i] -= 1
  end

  err = VecGetValues(vec, idx, vals)

  for i=1:length(idx)
    idx[i] += 1
  end

  return err
end

"""
  Method for AbstractVector
"""
function get_values1!(vec::AbstractVector, idx::Array{PetscInt}, vals::Array)
  for i=1:length(idx)
    vals[i] = vec[idx[i]]
  end

  return PetscErrorCode(0)
end

"""
  Calls [`VecAssemblyBegin`](@ref).  No-op for Julia vectors

  **Inputs**

   * vec: AbstractVector
"""
function assembly_begin(vec::PetscVec)
  VecAssemblyBegin(vec)
end

function assembly_begin(vec::AbstractVector)
end

"""
  Calls [`VecAssemblyEnd`](@ref).  No-op for Julia vectors

  **Inputs**

   * vec: AbstractVector
"""
function assembly_end(vec::PetscVec)
  VecAssemblyEnd(vec)
end

function assembly_end(vec::AbstractVector)
end

"""
  Print non-Petsc vector to a given IO (a Julia IO, not a Petsc IO).  Defaults
  to printing to STDOUT.

  **Inputs**

   * b: AbstractVector
   * f: IO, defaults to STDOUT
"""
function PetscView(b::AbstractVector, f::IO=STDOUT)
  println(f, "b = \n", a)
end

"""
  Size of local part of vector

  **Inputs**

   * A: AbstractVector
"""
function size_local(A::PetscVec)
  return (VecGetLocalSize(A), )
end

function size_local(A::AbstractVector)
  return size(A)
end

"""
  Size of global vector

  **Inputs**

   * A: AbstractVector
"""
function size_global(A::PetscVec)
  return (VecGetSize(A), )
end

function size_global(A::AbstractVector)
  return size(A)
end

"""
  Length of local part of vector

  **Inputs**

   * A: AbstractVector
"""
function length_local(A::AllVectors)
  return size_local(A)[1]
end

"""
  Length of global vector

  **Inputs**

   * A: AbstractVector

"""
function length_global(A::AllVectors)
  return size_local(A)[1]
end

"""
  Returns a UnitRange containing the (1-based) global indices owned by this
  process.

  **Inputs**

   * A: AbstractVector

  **Outputs**

   * rng: UnitRange
"""
function local_indices(vec::PetscVec)

  low, high = VecGetOwnershipRange(vec)
  return Int(low+1):Int(high)
end

function local_indices(vec::AbstractVector)
  return 1:length(vec)  #TODO: update in future versions of Julia
end

"""
  Fill a vector with zeros

  **Inputs**

   * A: AbstractVector
"""
function fill_zero!(A::PetscVec)
  VecSet(A, PetscScalar(0.0))
end

function fill_zero!(A::AbstractVector)
  fill!(A, 0.0)
end

import Base.fill!

"""
  Base.fill! for Petsc vectors
"""
function fill!(A::PetscVec, a)
  VecSet(A, PetscScalar(a))
end

import Base.maximum

"""
  Base.maximum
"""
function maximum(x::PetscVec)
  mval, idx = VecMax(x)
  return mval
end

import Base.minimum

"""
  Base.minimum
"""
function minimum(x::PetscVec)
  mval, idx = VecMin(x)
  return mval
end

import Base.copy

"""
  Base.copy
"""
function copy(vec::PetscVec)
  b2 = VecDuplicate(vec)
  VecCopy(vec, b2)

  return b2
end

import Base.copy!
"""
  Base.copy!
"""
function copy!(dest::PetscVec, src::PetscVec)
  VecCopy(src, dest)
end

#TODO:
#      +, .*, -, (check for in-place versions),
#

#TODO: extend blas routines

import Base.scale!
"""
  Base.scale!
"""
function scale!(vec::PetscVec, a::Number)
  _a = PetscScalar(a)
  VecScale(vec, _a)
end

"""
  Add a given value to all elements of the vector

  **Inputs**

   * A: AbstractVector
   * a: number to shift by
"""
function diagonal_shift!(A::PetscVec, a::Number)
  VecShift(A, PetscScalar(a))
end

function diagonal_shift!(A::AbstractVector, a::Number)
  @simd for i=1:length(A)
    A[i] += a
  end
end

import Base.norm
"""
  Base.norm
"""
function norm(A::PetscVec, p=2)
  if p == 1
    _p = NORM_1
  elseif p == 2
    _p = NORM_2
  elseif p == Inf
    _p = NORM_INFINITY
  else
    error("Petsc Vectors do not support norm p = $p")
  end

  VecNorm(A, _p)
end

import Base.dot
"""
  Base.dot

  Dot product where the *first* vector is conjugated.  This is is the reverse
  of VecDot, where the *second* vector is conjugated
"""
function dot(x::PetscVec, y::PetscVec)
  VecDot(y, x)
end

import Base.sum
"""
  Base.sum
"""
function sum(x::PetscVec)
  VecSum(x)
end


text:
"""
    translate_symlinks(root::AbstractString; verbose::Bool=false)

Walks through the root directory given within `root`, finding all symlinks that
point to an absolute path within `root`, and rewriting them to be a relative
symlink instead, increasing relocatability.
"""
function translate_symlinks(root::AbstractString; verbose::Bool=false)
    for f in collect_files(root, islink)
        link_target = readlink(f)
        if isabspath(link_target) && startswith(link_target, "/workspace")
            new_link_target = relpath(link_target, replace(dirname(f), root => "/workspace/destdir"))
            if verbose
                @info("Translating $f to point to $(new_link_target)")
            end
            rm(f; force=true)
            symlink(new_link_target, f)
        end
    end
end

"""
    warn_deadlinks(root::AbstractString)

Walks through the given `root` directory, finding broken symlinks and warning
the user about them.  This is used to catch instances such as a build recipe
copying a symlink that points to a dependency; by doing so, it implicitly
breaks relocatability.
"""
function warn_deadlinks(root::AbstractString)
    for f in collect_files(root, islink; exclude_externalities=false)
        link_target = readlink(f)
        if !startswith(link_target, "/")
            # If the link is relative, prepend the dirname of `f` to
            # `link_target`, otherwise `isfile(link_target)` will always be
            # false, as the test isn't performed in dirname(f).  Why Not using
            # `isabspath`?  Because that command is platform-dependent, but we
            # always build in a Unix-like environment.
            link_target = joinpath(dirname(f), link_target)
        end
        if !ispath(link_target)
            @warn("Broken symlink: $(relpath(f, root))")
        end
    end
end


text:
using Distributions
using Base.Test

# Test the consistency between the recursive and nonrecursive computation of the pdf
# of the Binomial distribution
srand(1234)
for (p, n) in [(0.6, 10), (0.8, 6), (0.5, 40), (0.04, 20), (1., 100), (0., 10), (0.999999, 1000), (1e-7, 1000)]

    d = Binomial(n, p)

    a = pdf(d, 0:n)
    for t=0:n
        @test pdf(d, t) ≈ a[1+t]
    end

    li = rand(0:n, 2)
    rng = minimum(li):maximum(li)
    b = pdf(d, rng)
    for t in rng
        @test pdf(d, t) ≈ b[t - first(rng) + 1]
    end

end


text:
using Javis

function ground(args...)
    background("white")
    sethue("white")
end

video = Video(400, 450)

U = 50
SHIFT = Point(0, -1.7U)

function tangram(obj)
    if obj == :triangle1
        sethue("#C9E891")
        [O, Point(2U, 0), Point(0, 2U)]
    elseif obj == :triangle2
        sethue("#69C0D2")
        [O, Point(-2U, 0), Point(0, 2U)]
    elseif obj == :triangle3
        sethue("#F8CB2D")
        [O, Point(-U, U), Point(U, U)]
    elseif obj == :triangle4
        sethue("#F195C8")
        [O, Point(-U, 0), Point(0, -U)]
    elseif obj == :triangle5
        sethue("#F9F224")
        [O, Point(U, 0), Point(0, -U)]
    elseif obj == :box
        sethue("#A796C2")
        [O, Point(0, -U), Point(U, -U), Point(U, 0)]
    elseif obj == :parallelogram
        sethue("#EF3E62")
        [O, Point(U, -U), Point(0, -U), Point(-U, 0)]
    end
end

javis(
    video,
    [
        BackgroundAction(1:400, ground),
        Action(
            1:400,
            (args...) -> poly(tangram(:triangle1), :fill, close = true);
            subactions = [
                SubAction(1:1, Translation(O, SHIFT)),
                SubAction(50:60, Translation(O, Point(0, 2U))),
                SubAction(60:70, Rotation(0.0, -π / 2)),
                SubAction(70:80, Translation(O, Point(-2U, 0))),
                SubAction(80:90, Translation(O, Point(U, -U))),
            ],
        ),
        Action(
            1:400,
            (args...) -> poly(tangram(:triangle2), :fill, close = true);
            subactions = [SubAction(1:1, Translation(O, SHIFT))],
        ),
        Action(
            1:400,
            (args...) -> poly(tangram(:triangle3), :fill, close = true);
            subactions = [
                SubAction(1:1, Translation(O, SHIFT)),
                SubAction(1:1, Translation(O, Point(0, -2U))),
                SubAction(100:110, Translation(O, Point(-3U, 0))),
                SubAction(110:120, Rotation(0.0, -π / 2)),
                SubAction(120:130, Translation(O, Point(-5.5U, 0))),
                SubAction(130:140, Translation(O, Point(0, U))),
            ],
        ),
        Action(
            1:400,
            (args...) -> poly(tangram(:triangle5), :fill, close = true);
            subactions = [
                SubAction(1:1, Translation(O, SHIFT)),
                SubAction(1:1, Translation(O, Point(U, 0))),
                SubAction(150:160, Translation(O, Point(0, 4.5 * U))),
                SubAction(160:170, Rotation(0.0, -3π / 4)),
                SubAction(
                    170:180,
                    Translation(
                        O,
                        Point(sqrt((2 - √2 / 2)^2 / 2) * U, -sqrt((2 - √2 / 2)^2 / 2) * U),
                    ),
                ),
            ],
        ),
        Action(
            1:400,
            (args...) -> poly(tangram(:parallelogram), :fill, close = true);
            subactions = [
                SubAction(1:1, Translation(O, SHIFT)),
                SubAction(1:1, Translation(O, Point(-U, 0))),
                SubAction(200:210, Translation(O, Point(U, -U))),
                SubAction(210:220, Translation(O, Point(U, 0))),
                SubAction(220:230, Rotation(0.0, π / 2)),
                SubAction(230:240, Translation(O, Point(2U, 0))),
                SubAction(240:250, Translation(O, Point(0, U))),
            ],
        ),
        Action(
            1:400,
            (args...) -> poly(tangram(:triangle4), :fill, close = true);
            subactions = [
                SubAction(1:1, Translation(O, SHIFT)),
                SubAction(250:260, Translation(O, Point(0, -U))),
                SubAction(260:270, Translation(O, Point(2U, 0))),
                SubAction(270:280, Rotation(0.0, -π / 2)),
                SubAction(280:290, Translation(O, Point(-4U, 0))),
                SubAction(290:300, Translation(O, Point(0.5 * U, -0.5 * U))),
            ],
        ),
        Action(
            1:400,
            (args...) -> poly(tangram(:box), :fill, close = true);
            subactions = [
                SubAction(1:1, Translation(O, SHIFT)),
                SubAction(300:310, Translation(O, Point(-U, 0))),
                SubAction(310:320, Rotation(0.0, -π / 6)),
            ],
        ),
    ];
    pathname = "logo.gif",
)


text:
include("rachford_rice.jl")
include("bubble_pressure.jl")
include("LLE_pressure.jl")
include("three_phase.jl")
include("crit_mix.jl")
include("PT_flash.jl")
export bubble_pressure, LLE_pressure, three_phase, crit_mix


text:
export GradientOp

"""
    gradOp(T::Type, shape::NTuple{1,Int64})

1d gradient operator for an array of size `shape`
"""
GradientOp(T::Type, shape::NTuple{1,Int64}) = GradientOp(T,shape,1)

"""
    gradOp(T::Type, shape::NTuple{2,Int64})

2d gradient operator for an array of size `shape`
"""
function GradientOp(T::Type, shape::NTuple{2,Int64})
  return vcat( GradientOp(T,shape,1), GradientOp(T,shape,2) )
end

"""
    gradOp(T::Type, shape::NTuple{3,Int64})

3d gradient operator for an array of size `shape`
"""
function GradientOp(T::Type, shape::NTuple{3,Int64})
  return vcat( GradientOp(T,shape,1), GradientOp(T,shape,2), GradientOp(T,shape,3) )
end

"""
    gradOp(T::Type, shape::NTuple{N,Int64}, dim::Int64) where N

directional gradient operator along the dimension `dim`
for an array of size `shape`
"""
function GradientOp(T::Type, shape::NTuple{N,Int64}, dim::Int64) where N
  nrow = div( (shape[dim]-1)*prod(shape), shape[dim] )
  ncol = prod(shape)
  return LinearOperator{T}(nrow, ncol, false, false,
                    x-> grad(x,shape,dim),
                    x-> grad_t(x,shape,dim),
                    x-> grad_t(x,shape,dim) )
end

# directional gradients
function grad(img::T, shape::NTuple{1,Int64}, dim::Int64) where T<:AbstractVector
  return img[1:end-1].-img[2:end]
end

function grad(img::T, shape::NTuple{2,Int64}, dim::Int64) where T<:AbstractVector
  img = reshape(img,shape)

  if dim==1
    return vec(img[1:end-1,:].-img[2:end,:])
  end

  return vec(img[:,1:end-1].-img[:,2:end])
end

function grad(img::T, shape::NTuple{3,Int64}, dim::Int64) where T<:AbstractVector
  img = reshape(img,shape)

  if dim==1
    return vec(img[1:end-1,:,:].-img[2:end,:,:])
  elseif dim==2
    return vec(img[:,:,1:end-1].-img[:,:,2:end])
  end

  return vec(img[:,:,1:end-1].-img[:,:,2:end])
end

# adjoint of directional gradients
function grad_t(g::T, shape::NTuple{1,Int64}, dim::Int64) where T<:AbstractVector
  nx = shape[1]
  img = similar(g,nx)
  img .= zero(eltype(g))

  img[1:nx-1] .= g
  img[2:nx] .-= g

  return img
end

function grad_t(g::T, shape::NTuple{2,Int64}, dim::Int64) where T<:AbstractVector
  nx,ny = shape
  img = similar(g,nx,ny)
  img .= zero(eltype(g))

  if dim==1
    g = reshape(g,nx-1,ny)
    img[1:nx-1,:] .= g
    img[2:nx,:] .-= g
  else
    g = reshape(g,nx,ny-1)
    img[:,1:ny-1] .= g
    img[:,2:ny] .-= g
  end
  return vec(img)
end

function grad_t(g::T, shape::NTuple{3,Int64}, dim::Int64) where T<:AbstractVector
  nx,ny,nz = shape
  img = similar(g,nx,ny,nz)
  img .= zero(eltype(g))

  if dim==1
    g = reshape(g,nx-1,nx,nz)
    img[1:nx-1,:,:] .= g
    img[2:nx,:,:] .-= g
  elseif dim==2
    g = reshape(g,nx,ny-1,nz)
    img[:,1:ny-1,:] .= g
    img[:,2:ny,:] .-= g
  else
    g = reshape(g,nx,nx,nz-1)
    img[:,:,1:nz-1] .= g
    img[:,:,2:nz] .-= g
  end
  return vec(img)
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: mobile
using AWS.Compat
using AWS.UUIDs

"""
    create_project()
    create_project(params::Dict{String,<:Any})

 Creates an AWS Mobile Hub project.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"contents"`:  ZIP or YAML file which contains configuration settings to be used when
  creating the project. This may be the contents of the file downloaded from the URL provided
  in an export project operation.
- `"name"`:  Name of the project.
- `"region"`:  Default region where project resources should be created.
- `"snapshotId"`:  Unique identifier for an exported snapshot of project configuration.
  This snapshot identifier is included in the share URL when a project is exported.
"""
create_project(; aws_config::AbstractAWSConfig=global_aws_config()) = mobile("POST", "/projects"; aws_config=aws_config)
create_project(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mobile("POST", "/projects", params; aws_config=aws_config)

"""
    delete_project(project_id)
    delete_project(project_id, params::Dict{String,<:Any})

 Delets a project in AWS Mobile Hub.

# Arguments
- `project_id`:  Unique project identifier.

"""
delete_project(projectId; aws_config::AbstractAWSConfig=global_aws_config()) = mobile("DELETE", "/projects/$(projectId)"; aws_config=aws_config)
delete_project(projectId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mobile("DELETE", "/projects/$(projectId)", params; aws_config=aws_config)

"""
    describe_bundle(bundle_id)
    describe_bundle(bundle_id, params::Dict{String,<:Any})

 Get the bundle details for the requested bundle id.

# Arguments
- `bundle_id`:  Unique bundle identifier.

"""
describe_bundle(bundleId; aws_config::AbstractAWSConfig=global_aws_config()) = mobile("GET", "/bundles/$(bundleId)"; aws_config=aws_config)
describe_bundle(bundleId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mobile("GET", "/bundles/$(bundleId)", params; aws_config=aws_config)

"""
    describe_project(project_id)
    describe_project(project_id, params::Dict{String,<:Any})

 Gets details about a project in AWS Mobile Hub.

# Arguments
- `project_id`:  Unique project identifier.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"syncFromResources"`:  If set to true, causes AWS Mobile Hub to synchronize information
  from other services, e.g., update state of AWS CloudFormation stacks in the AWS Mobile Hub
  project.
"""
describe_project(projectId; aws_config::AbstractAWSConfig=global_aws_config()) = mobile("GET", "/project", Dict{String, Any}("projectId"=>projectId); aws_config=aws_config)
describe_project(projectId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mobile("GET", "/project", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("projectId"=>projectId), params)); aws_config=aws_config)

"""
    export_bundle(bundle_id)
    export_bundle(bundle_id, params::Dict{String,<:Any})

 Generates customized software development kit (SDK) and or tool packages used to integrate
mobile web or mobile app clients with backend AWS resources.

# Arguments
- `bundle_id`:  Unique bundle identifier.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"platform"`:  Developer desktop or target application platform.
- `"projectId"`:  Unique project identifier.
"""
export_bundle(bundleId; aws_config::AbstractAWSConfig=global_aws_config()) = mobile("POST", "/bundles/$(bundleId)"; aws_config=aws_config)
export_bundle(bundleId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mobile("POST", "/bundles/$(bundleId)", params; aws_config=aws_config)

"""
    export_project(project_id)
    export_project(project_id, params::Dict{String,<:Any})

 Exports project configuration to a snapshot which can be downloaded and shared. Note that
mobile app push credentials are encrypted in exported projects, so they can only be shared
successfully within the same AWS account.

# Arguments
- `project_id`:  Unique project identifier.

"""
export_project(projectId; aws_config::AbstractAWSConfig=global_aws_config()) = mobile("POST", "/exports/$(projectId)"; aws_config=aws_config)
export_project(projectId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mobile("POST", "/exports/$(projectId)", params; aws_config=aws_config)

"""
    list_bundles()
    list_bundles(params::Dict{String,<:Any})

 List all available bundles.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`:  Maximum number of records to list in a single response.
- `"nextToken"`:  Pagination token. Set to null to start listing bundles from start. If
  non-null pagination token is returned in a result, then pass its value in here in another
  request to list more bundles.
"""
list_bundles(; aws_config::AbstractAWSConfig=global_aws_config()) = mobile("GET", "/bundles"; aws_config=aws_config)
list_bundles(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mobile("GET", "/bundles", params; aws_config=aws_config)

"""
    list_projects()
    list_projects(params::Dict{String,<:Any})

 Lists projects in AWS Mobile Hub.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`:  Maximum number of records to list in a single response.
- `"nextToken"`:  Pagination token. Set to null to start listing projects from start. If
  non-null pagination token is returned in a result, then pass its value in here in another
  request to list more projects.
"""
list_projects(; aws_config::AbstractAWSConfig=global_aws_config()) = mobile("GET", "/projects"; aws_config=aws_config)
list_projects(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mobile("GET", "/projects", params; aws_config=aws_config)

"""
    update_project(project_id)
    update_project(project_id, params::Dict{String,<:Any})

 Update an existing project.

# Arguments
- `project_id`:  Unique project identifier.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"contents"`:  ZIP or YAML file which contains project configuration to be updated. This
  should be the contents of the file downloaded from the URL provided in an export project
  operation.
"""
update_project(projectId; aws_config::AbstractAWSConfig=global_aws_config()) = mobile("POST", "/update", Dict{String, Any}("projectId"=>projectId); aws_config=aws_config)
update_project(projectId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mobile("POST", "/update", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("projectId"=>projectId), params)); aws_config=aws_config)


text:
module MKL

using MKL_jll

JULIA_VER_NEEDED = v"1.7.0-DEV.641"
VERSION > JULIA_VER_NEEDED && using LinearAlgebra

if Base.USE_BLAS64
    const MKLBlasInt = Int64
else
    const MKLBlasInt = Int32
end

@enum Threading begin
    THREADING_INTEL
    THREADING_SEQUENTIAL
    THREADING_PGI
    THREADING_GNU
    THREADING_TBB
end

@enum Interface begin
    INTERFACE_LP64
    INTERFACE_ILP64
    INTERFACE_GNU
end

function set_threading_layer(layer::Threading = THREADING_INTEL)
    err = ccall((:MKL_Set_Threading_Layer, libmkl_rt), Cint, (Cint,), layer)
    err == -1 && throw(ErrorException("return value was -1"))
    return nothing
end

function set_interface_layer(interface = Base.USE_BLAS64 ? INTERFACE_ILP64 : INTERFACE_LP64)
    err = ccall((:MKL_Set_Interface_Layer, libmkl_rt), Cint, (Cint,), interface)
    err == -1 && throw(ErrorException("return value was -1"))
    return nothing
end

function __init__()
    if MKL_jll.is_available()
        set_threading_layer()
        set_interface_layer()
        VERSION > JULIA_VER_NEEDED && BLAS.lbt_forward(libmkl_rt, clear=true)
    end
end

function mklnorm(x::Vector{Float64})
    ccall((:dnrm2_, libmkl_rt), Float64,
          (Ref{MKLBlasInt}, Ptr{Float64}, Ref{MKLBlasInt}),
          length(x), x, 1)
end

VERSION > JULIA_VER_NEEDED && include("install.jl")

end # module