text:
### Experimental ###
#################################

@params struct StressTemp{T}
    VTu::AbstractVector{T}
    Tu::AbstractVector{T}
    Te::AbstractMatrix{T}
    global_dofs::AbstractVector{Int}
end
function StressTemp(solver)
    @unpack u, problem = solver
    @unpack dh = problem.ch
    T = eltype(u)
    dim = TopOptProblems.getdim(problem)
    k = dim == 2 ? 3 : 6
    VTu = zero(MVector{k,T})
    Tu = similar(VTu)
    n_basefuncs = getnbasefunctions(solver.elementinfo.cellvalues)
    Te = zero(MMatrix{k,dim * n_basefuncs,T})
    global_dofs = zeros(Int, ndofs_per_cell(dh))

    return StressTemp(VTu, Tu, Te, global_dofs)
end
Zygote.@nograd StressTemp

@params mutable struct MacroVonMisesStress{T} <: AbstractFunction{T}
    utMu::AbstractVector{T}
    Mu::AbstractArray
    sigma_vm::AbstractVector{T}
    solver::Any
    global_dofs::AbstractVector{<:Integer}
    stress_temp::StressTemp
    fevals::Int
    reuse::Bool
    maxfevals::Int
end
function MacroVonMisesStress(solver; reuse=false, maxfevals=10^8)
    T = eltype(solver.u)
    dh = solver.problem.ch.dh
    k = ndofs_per_cell(dh)
    N = getncells(dh.grid)
    global_dofs = zeros(Int, k)
    Mu = zeros(SVector{k,T}, N)
    utMu = zeros(T, N)
    stress_temp = StressTemp(solver)
    sigma_vm = similar(utMu)

    return MacroVonMisesStress(
        utMu, Mu, sigma_vm, solver, global_dofs, stress_temp, 0, reuse, maxfevals
    )
end
function (ls::MacroVonMisesStress{T})(x) where {T}
    @unpack sigma_vm, Mu, utMu, stress_temp = ls
    @unpack reuse, solver, global_dofs = ls
    @unpack penalty, problem, xmin = solver
    E0 = problem.E
    ls.fevals += 1

    @assert length(global_dofs) == ndofs_per_cell(solver.problem.ch.dh)
    solver.vars .= x
    if !reuse
        solver()
        fill_Mu_utMu!(Mu, utMu, solver, stress_temp)
    end
    sigma_vm .= get_sigma_vm.(get_ρ.(x, Ref(penalty), xmin), utMu)
    return copy(sigma_vm)
end
function ChainRulesCore.rrule(vonmises::MacroVonMisesStress, x)
    @unpack Mu, utMu, stress_temp = vonmises
    @unpack reuse, solver, global_dofs = vonmises
    @unpack penalty, problem, u, xmin = solver
    dh = getdh(problem)
    @unpack Kes = solver.elementinfo
    E0 = problem.E
    vonmises.fevals += 1

    sigma_vm = vonmises(x)
    return sigma_vm,
    Δ -> (
        nothing,
        begin
            for e in 1:length(Mu)
                ρe = get_ρ(x[e], penalty, xmin)
                Mu[e] *= Δ[e] * ρe^2 / sigma_vm[e]
            end
            lhs = backsolve!(solver, Mu, global_dofs)
            map(1:length(Δ)) do e
                ρe, dρe = get_ρ_dρ(x[e], penalty, xmin)
                celldofs!(global_dofs, dh, e)
                t1 = Δ[e] * sigma_vm[e] / ρe * dρe
                @views t2 =
                    -dot(lhs[global_dofs], bcmatrix(Kes[e]) * u[global_dofs]) * dρe
                return t1 + t2
            end
        end,
    )
end
#getdim(f::MacroVonMisesStress) = length(f.sigma_vm)

###################################

@params struct MicroVonMisesStress{T} <: AbstractFunction{T}
    vonmises::MacroVonMisesStress{T}
end
function MicroVonMisesStress(args...; kwargs...)
    return MicroVonMisesStress(MacroVonMisesStress(args...; kwargs...))
end
function (f::MicroVonMisesStress)(x)
    @unpack vonmises = f
    @unpack sigma_vm, Mu, utMu, stress_temp = vonmises
    @unpack reuse, solver, global_dofs = vonmises
    @unpack penalty, problem, xmin = solver
    vonmises.fevals += 1
    @assert length(global_dofs) == ndofs_per_cell(solver.problem.ch.dh)
    solver.vars .= x
    if !reuse
        solver()
        fill_Mu_utMu!(Mu, utMu, solver, stress_temp)
    end
    out = sqrt.(utMu)
    sigma_vm .= get_ρ.(x, Ref(penalty), xmin) .* out
    return out
end
function ChainRulesCore.rrule(f::MicroVonMisesStress, x)
    @unpack vonmises = f
    @unpack Mu, utMu, stress_temp, sigma_vm = vonmises
    @unpack reuse, solver, global_dofs = vonmises
    @unpack penalty, problem, u, xmin = solver
    dh = getdh(problem)
    @unpack Kes = solver.elementinfo
    E0 = problem.E
    vonmises.fevals += 1

    out = f(x)
    return out,
    Δ -> (
        nothing,
        begin
            for e in 1:length(Mu)
                ρe = get_ρ(x[e], penalty, xmin)
                Mu[e] *= Δ[e] * ρe / sigma_vm[e]
            end
            lhs = backsolve!(solver, Mu, global_dofs)
            map(1:length(Δ)) do e
                ρe, dρe = get_ρ_dρ(x[e], penalty, xmin)
                celldofs!(global_dofs, dh, e)
                return -dot(lhs[global_dofs], bcmatrix(Kes[e]) * u[global_dofs]) * dρe
            end
        end,
    )
end

#########################################
# stress utils

function backsolve!(solver, Mu, global_dofs)
    dh = getdh(solver.problem)
    solver.rhs .= 0
    for i in 1:length(Mu)
        celldofs!(global_dofs, dh, i)
        solver.rhs[global_dofs] .+= Mu[i]
    end
    solver(; assemble_f=false)
    return solver.lhs
end

get_sigma_vm(ρ_e, utMu_e) = ρ_e * sqrt(utMu_e)

@inline function get_ϵ(u, ∇ϕ, i, j)
    return 1 / 2 * (u[i] * ∇ϕ[j] + u[j] * ∇ϕ[i])
end
@inline function apply_T!(Tu, u, dh, cellidx, global_dofs, cellvalues, ν, ::Val{2})
    # assumes cellvalues is initialized before passing
    Tu[:] .= 0
    # Generalize to higher order field basis functions
    q_point = 1
    n_basefuncs = getnbasefunctions(cellvalues)
    dim = 2
    for a in 1:n_basefuncs
        ∇ϕ = shape_gradient(cellvalues, q_point, a)
        _u = @view u[(@view global_dofs[(dim * (a - 1) + 1):(a * dim)])]
        ϵ_11 = get_ϵ(_u, ∇ϕ, 1, 1)
        ϵ_22 = get_ϵ(_u, ∇ϕ, 2, 2)
        ϵ_12 = get_ϵ(_u, ∇ϕ, 1, 2)

        ϵ_sum = ϵ_11 + ϵ_22

        temp1 = ν / (1 - ν^2)
        temp2 = ν * (1 + ν)

        Tu[1] += temp1 * ϵ_sum + temp2 * ϵ_11 # σ[1,1] / E
        Tu[2] += temp1 * ϵ_sum + temp2 * ϵ_22 # σ[2,2] / E
        Tu[3] += temp2 * ϵ_12 # σ[1,2] / E
    end
    return Tu
end
@inline function apply_T!(Tu, u, dh, global_dofs, cellvalues, ν, ::Val{3})
    # assumes cellvalues is initialized before passing
    Tu[:] .= 0
    q_point = 1
    n_basefuncs = getnbasefunctions(cellvalues)
    dim = 3
    for a in 1:n_basefuncs
        ∇ϕ = shape_gradient(cellvalues, q_point, a)
        _u = @view u[(@view global_dofs[(dim * (a - 1) + 1):(a * dim)])]
        ϵ_11 = get_ϵ(_u, ∇ϕ, 1, 1)
        ϵ_22 = get_ϵ(_u, ∇ϕ, 2, 2)
        ϵ_33 = get_ϵ(_u, ∇ϕ, 3, 3)
        ϵ_12 = get_ϵ(_u, ∇ϕ, 1, 2)
        ϵ_23 = get_ϵ(_u, ∇ϕ, 2, 3)
        ϵ_31 = get_ϵ(_u, ∇ϕ, 3, 1)

        ϵ_sum = ϵ_11 + ϵ_22 + ϵ_33

        temp1 = ν / (1 - ν^2)
        temp2 = ν * (1 + ν)

        Tu[1] += temp1 * ϵ_sum + temp2 * ϵ_11 # σ[1,1] / E
        Tu[2] += temp1 * ϵ_sum + temp2 * ϵ_22 # σ[2,2] / E
        Tu[3] += temp1 * ϵ_sum + temp2 * ϵ_33 # σ[3,3] / E
        Tu[4] += temp2 * ϵ_12 # σ[1,2] / E
        Tu[5] += temp2 * ϵ_23 # σ[2,3] / E
        Tu[6] += temp2 * ϵ_31 # σ[3,1] / E
    end
    return Tu
end

@inline function fill_T!(T, ::Val{3}, cellvalues, E0, ν)
    # assumes cellvalues is initialized before passing
    dim = 3
    temp1 = E0 * ν / (1 - ν^2)
    temp2 = E0 * ν * (1 + ν)
    q_point = 1
    n_basefuncs = size(T, 2) ÷ dim
    @assert size(T, 1) == 6
    for a in 1:n_basefuncs
        ∇ϕ = shape_gradient(cellvalues, q_point, a)
        cols = (dim * (a - 1) + 1):(dim * a)
        T[1, cols[1]] = (temp1 + temp2) * ∇ϕ[1]
        T[2, cols[1]] = temp1 * ∇ϕ[1]
        T[3, cols[1]] = temp1 * ∇ϕ[1]
        T[4, cols[1]] = temp2 * ∇ϕ[2] / 2
        T[5, cols[1]] = 0
        T[6, cols[1]] = temp2 * ∇ϕ[3] / 2

        T[1, cols[2]] = temp1 * ∇ϕ[2]
        T[2, cols[2]] = (temp1 + temp2) * ∇ϕ[2]
        T[3, cols[2]] = temp1 * ∇ϕ[2]
        T[4, cols[2]] = temp2 * ∇ϕ[1] / 2
        T[5, cols[2]] = temp2 * ∇ϕ[3] / 2
        T[6, cols[2]] = 0

        T[1, cols[3]] = temp1 * ∇ϕ[3]
        T[2, cols[3]] = temp1 * ∇ϕ[3]
        T[3, cols[3]] = (temp1 + temp2) * ∇ϕ[3]
        T[4, cols[3]] = 0
        T[5, cols[3]] = temp2 * ∇ϕ[2] / 2
        T[6, cols[3]] = temp2 * ∇ϕ[1] / 2
    end
    return T
end

@inline function fill_T!(T, ::Val{2}, cellvalues, E0, ν)
    # assumes cellvalues is initialized before passing
    dim = 2
    temp1 = E0 * ν / (1 - ν^2)
    temp2 = E0 * ν * (1 + ν)
    q_point = 1
    n_basefuncs = size(T, 2) ÷ dim
    @assert size(T, 1) == 3
    for a in 1:n_basefuncs
        ∇ϕ = shape_gradient(cellvalues, q_point, a)
        cols = (dim * (a - 1) + 1):(dim * a)
        T[1, cols[1]] = (temp1 + temp2) * ∇ϕ[1]
        T[2, cols[1]] = temp1 * ∇ϕ[1]
        T[3, cols[1]] = temp2 * ∇ϕ[2] / 2

        T[1, cols[2]] = temp1 * ∇ϕ[2]
        T[2, cols[2]] = (temp1 + temp2) * ∇ϕ[2]
        T[3, cols[2]] = temp2 * ∇ϕ[1] / 2
    end
    return T
end

function fill_Mu_utMu!(Mu, utMu, solver, stress_temp::StressTemp)
    @unpack problem, elementinfo, u = solver
    @unpack ch, ν, E = problem
    @unpack dh = ch
    @unpack VTu, Tu, Te, global_dofs = stress_temp

    _fill_Mu_utMu!(Mu, utMu, dh, elementinfo, u, E, ν, global_dofs, Tu, VTu, Te)
    return Mu, utMu
end

@inline function _fill_Mu_utMu!(
    Mu::AbstractVector,
    utMu::AbstractVector{T},
    dh::DofHandler{3},
    elementinfo,
    u,
    E0,
    ν,
    global_dofs=zeros(Int, ndofs_per_cell(dh)),
    Tu=zeros(T, 6),
    VTu=zeros(T, 6),
    Te=zeros(T, 6, ndofs_per_cell(dh)),
) where {T}
    dim = 3
    @unpack cellvalues = elementinfo
    for (cellidx, cell) in enumerate(CellIterator(dh))
        reinit!(cellvalues, cell)
        # Same for all elements
        fill_T!(Te, Val(3), cellvalues, E0, ν)
        celldofs!(global_dofs, dh, cellidx)
        @views mul!(Tu, Te, u[global_dofs])

        VTu[1] = Tu[1] - Tu[2] / 2 - Tu[3] / 2
        VTu[2] = -Tu[1] / 2 + Tu[2] - Tu[3] / 2
        VTu[3] = -Tu[1] / 2 - Tu[2] / 2 + Tu[3]
        VTu[4] = 3 * Tu[4]
        VTu[5] = 3 * Tu[5]
        VTu[6] = 3 * Tu[6]

        utMu_e = dot(Tu, VTu)
        @assert utMu_e >= 0
        utMu[cellidx] = utMu_e
        Mu[cellidx] = Te' * VTu
    end
    return Mu, utMu
end

@inline function _fill_Mu_utMu!(
    Mu::AbstractVector,
    utMu::AbstractVector{T},
    dh::DofHandler{2},
    elementinfo,
    u,
    E0,
    ν,
    global_dofs=zeros(Int, ndofs_per_cell(dh)),
    Tu=zeros(T, 3),
    VTu=zeros(T, 3),
    Te=zeros(T, 3, ndofs_per_cell(dh)),
) where {T}
    dim = 2
    @unpack cellvalues = elementinfo
    for (cellidx, cell) in enumerate(CellIterator(dh))
        reinit!(cellvalues, cell)
        fill_T!(Te, Val(2), cellvalues, E0, ν)
        celldofs!(global_dofs, dh, cellidx)
        @views mul!(Tu, Te, u[global_dofs])

        VTu[1] = Tu[1] - Tu[2] / 2
        VTu[2] = -Tu[1] / 2 + Tu[2]
        VTu[3] = 3 * Tu[3]

        utMu_e = dot(Tu, VTu)
        @assert utMu_e >= 0
        utMu[cellidx] = utMu_e
        Mu[cellidx] = Te' * VTu
    end
    return Mu, utMu
end


text:
import MitosisStochasticDiffEq as MSDE
using Mitosis
using Test, Random
using LinearAlgebra
using OrdinaryDiffEq
Random.seed!(12345)
@testset "backward Lyapunov filtering tests" begin
  # define SDE function
  f(u,p,t) = p[1]*u .+ p[2]
  g(u,p,t) = p[3]

  # time span
  tstart = 0.0
  tend = 1.0
  dt = 0.001
  trange = tstart:dt:tend

  # set of linear parameters Eq.~(2.2)
  B, β, σ̃ = -0.1, 0.2, 1.3
  plin = [B, β, σ̃]

  # initial values for ODE
  mynames = (:logscale, :μ, :Σ);
  c = randn()
  ν = randn()
  P = 10*rand()
  myvalues = [c, ν, P]
  NT = NamedTuple{mynames}(myvalues)

  # covariance filter
  kernel = MSDE.SDEKernel(Mitosis.AffineMap(B, β), Mitosis.ConstantMap(σ̃), trange, plin)
  message, solend = MSDE.backwardfilter(kernel, NT, alg=OrdinaryDiffEq.Tsit5())

  # Lyapunov filter
  message2, solend2 = MSDE.backwardfilter(kernel, NT,
    filter=MSDE.LyapunovFilter())

  @test message.ts == message2.ts
  @test P == message2.sol.PT
  @test ν == message2.sol.νT
  @test c == message2.sol.cT

  @test getindex.(message.soldis,1) ≈ getindex.(message2.soldis,1) rtol=1e-10
  @test getindex.(message.soldis,2) ≈ getindex.(message2.soldis,2) rtol=1e-10
  @test getindex.(message.soldis,3) ≈ getindex.(message2.soldis,3) rtol=1e-10
  @test isapprox(solend, solend2, rtol=1e-10)

  dim = 1
  m = 1
  B, β, σ̃ = fill(B, 1, 1), fill(β, 1), fill(σ̃, 1, 1)
  kernel = MSDE.SDEKernel(
      Mitosis.AffineMap(B, β), Mitosis.ConstantMap(σ̃), trange, [B, β, σ̃]
  )

  ν = fill(ν, 1)
  P = fill(P, 1, 1)
  NT = NamedTuple{mynames}([c, ν, P])

  message, solend  = MSDE.backwardfilter(kernel, NT, alg=OrdinaryDiffEq.Tsit5())
  # Lyapunov filter
  message2, solend2 = MSDE.backwardfilter(kernel, NT,
    filter=MSDE.LyapunovFilter())

  @test getindex.(message.soldis,1) ≈ getindex.(message2.soldis,1) rtol=1e-10
  @test getindex.(message.soldis,2) ≈ getindex.(message2.soldis,2) rtol=1e-10
  @test getindex.(message.soldis,3) ≈ getindex.(message2.soldis,3) rtol=1e-10
  @test isapprox(solend, solend2, rtol=1e-10)

  dim = 4
  m = 3
  B, β, σ̃ = [randn(dim,dim)-I, randn(dim), randn(dim,m)] # B, β, σtil
  kernel = MSDE.SDEKernel(
      Mitosis.AffineMap(B, β), Mitosis.ConstantMap(σ̃), trange, [B, β, σ̃]
  )

  c = randn()
  ν = randn(dim)
  P = randn(dim,dim)
  P = I + 0.1*P*P'
  NT = NamedTuple{mynames}([c, ν, P])

  message, solend  = MSDE.backwardfilter(kernel, NT,  alg=OrdinaryDiffEq.Tsit5())
  # Lyapunov filter
  message2, solend2 = MSDE.backwardfilter(kernel, NT,
    filter=MSDE.LyapunovFilter())

  @test getindex.(message.soldis,1) ≈ getindex.(message2.soldis,1) rtol=1e-9
  @test getindex.(message.soldis,2) ≈ getindex.(message2.soldis,2) rtol=1e-9
  @test getindex.(message.soldis,3) ≈ getindex.(message2.soldis,3) rtol=1e-9
  @test solend2 ≈ solend rtol=1e-9
end

ext: jl


text:
using Random

const ZOB_KEYS = begin
    local rng = MersenneTwister(1685)
    [rand(rng, UInt64) for i ∈ 1:64, j ∈ 1:14]
end

function zobrist(p::Piece, s::Square)::UInt64
    @inbounds ZOB_KEYS[s.val, p.val]
end

function zobep(s::Square)::UInt64
    @inbounds ZOB_KEYS[s.val, 7]
end

function zobcastle(castlerights::UInt8)
    @inbounds ZOB_KEYS[castlerights+1, 8]
end

function zobsidetomove()
    @inbounds ZOB_KEYS[64, 8]
end


text:
using GenomicAnnotations
using GenomicMaps
using ColorTypes

# You can add any kind of annotation that you want to display.
# Here, I add COG annotation:
function addcogs!(chr, filename)
    cogs = split.(readlines(filename), Ref('\t'))
    i = 1
    for gene in @genes(chr, :feature == "CDS")
        if gene.locus_tag == cogs[i][1]
            if occursin(r"\w", cogs[i][2])
                gene.cog = cogs[i][2]
            end
            i += 1
        end
    end
end

# Colour scheme for COG categories:
cogcolours = Dict("B"=>RGB{Float64}(1.0,0.630714,0.576563),
     "M"=>RGB{Float64}(0.756869,0.916499,0.965176),
     "I"=>RGB{Float64}(0.187839,0.54561,0.252343),
     "X"=>RGB{Float64}(0.540006,0.493982,0.813567),
     "Y"=>RGB{Float64}(0.0973617,0.285282,0.5329),
     "Z"=>RGB{Float64}(0.0418427,0.156645,0.341597),
     "L"=>RGB{Float64}(0.426131,0.0441442,0.0465628),
     "O"=>RGB{Float64}(0.518954,0.802339,0.930272),
     "F"=>RGB{Float64}(0.587882,0.865532,0.51112),
     "Q"=>RGB{Float64}(0.0,0.225356,0.101282),
     "D"=>RGB{Float64}(0.862653,0.958477,0.981395),
     "V"=>RGB{Float64}(0.188382,0.529206,0.795898),
     "U"=>RGB{Float64}(0.277786,0.635283,0.863472),
     "E"=>RGB{Float64}(0.711814,0.932724,0.629136),
     "T"=>RGB{Float64}(0.394211,0.72627,0.90426),
     "H"=>RGB{Float64}(0.32729,0.673206,0.326717),
     "P"=>RGB{Float64}(0.0232916,0.395886,0.180144),
     "G"=>RGB{Float64}(0.459895,0.779462,0.41097),
     "N"=>RGB{Float64}(0.641543,0.865092,0.94902),
     "K"=>RGB{Float64}(0.825431,0.118066,0.106858),
     "C"=>RGB{Float64}(0.835916,0.980813,0.770886),
     "R"=>RGB{Float64}(0.807625,0.787968,0.949453),
     "W"=>RGB{Float64}(0.137797,0.411028,0.686187),
     "A"=>RGB{Float64}(1.0,0.808314,0.771835),
     "S"=>RGB{Float64}(0.236943,0.0166779,0.407047),
     "J"=>RGB{Float64}(1.0,0.389569,0.336934))

# First download annotations for E. coli:
download("ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/005/845/GCA_000005845.2_ASM584v2/GCA_000005845.2_ASM584v2_genomic.gbff.gz", "ecoli.gbk.gz")

# Then read the annotations and add COGs:
chr = readgbk("ecoli.gbk.gz")[1]
addcogs!(chr, "ecoli_cogs.tsv")

# The output can be customised, see src/initialise.jl for all options. Here I
# provide a function that will be run on each gene to determine its colour:
colourby_cog = g -> unique(String.(split(get(g, :cog, ""), "")))
drawgenome(chr;
    outfile = "ecoli.svg",
    colourmap = cogcolours,
    colourfunction = colourby_cog,
    annotate = true,
    nbreaks = 40)


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id$
#*--------------------------------------------------------------------
type Heatex_LMTD
	Heatex_LMTD()=begin
		new(
			Heatex_Basic(),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Type of Heat Exchanger",
				:Valid=>["Counter Flow","Cocurrent Flow", "Shell and Tube"],
				:Default=>"Cocurrent Flow"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"LMTD Correction Factor Model",
				:Valid=>["Bowmann","Fakheri"],
				:Default=>"Bowmann"
			)),
			LMTD_Basic(Dict{Symbol,Any}(
				:Brief=>"LMTD Method of Calculation",
				:Symbol=>" "
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Capacity Ratio for LMTD Correction Fator",
				:Lower=>1e-6,
				:Hidden=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Non - Dimensional Variable for LMTD Correction Fator ",
				:Lower=>1e-6,
				:Hidden=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Non - Dimensional Variable for LMTD Correction Fator in Fakheri Equation",
				:Lower=>1e-6,
				:Hidden=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Non - Dimensional Variable for LMTD Correction Fator in Fakheri Equation",
				:Lower=>1e-6,
				:Symbol=>"\\phi",
				:Hidden=>true
			)),
			[
				:(Q = U*A*Method.LMTD*Method.Fc),
				:(Method.DT0 = _base_1.InletHot.T - _base_1.InletCold.T),
				:(Method.DTL = _base_1.OutletHot.T - _base_1.OutletCold.T),
				:(R=1),
				:(P=1),
				:(Phi = 1),
				:(Rho = 1),
				:(Method.Fc = 1),
				:(Method.DT0 = _base_1.InletHot.T - _base_1.OutletCold.T),
				:(Method.DTL = _base_1.OutletHot.T - _base_1.InletCold.T),
				:(R=1),
				:(P=1),
				:(Phi = 1),
				:(Rho = 1),
				:(Method.Fc = 1),
				:(Method.DT0 = _base_1.InletHot.T - _base_1.OutletCold.T),
				:(Method.DTL = _base_1.OutletHot.T - _base_1.InletCold.T),
				:(Phi = 1),
				:(Rho = 1),
				:(R*(_base_1.OutletCold.T - _base_1.InletCold.T ) = (_base_1.InletHot.T-_base_1.OutletHot.T)),
				:(P*(_base_1.InletHot.T- _base_1.InletCold.T)= (_base_1.OutletCold.T-_base_1.InletCold.T)),
				:(Method.Fc = (sqrt(2)*P)/((1-P)*ln( abs( ( 2-P*0.585786)/( 2-P*3.414214))))),
				:(Method.Fc = sqrt(R*R+1)*ln(abs((1-P*R)/(1-P)))/((1-R)*ln( abs( ( 2-P*(R+1-sqrt(R*R+1)))/ ( 2-P*(R + 1 + sqrt(R*R+1))))))),
				:(R*(_base_1.OutletCold.T - _base_1.InletCold.T ) = (_base_1.InletHot.T-_base_1.OutletHot.T)),
				:(P*(_base_1.InletHot.T- _base_1.InletCold.T)= (_base_1.OutletCold.T-_base_1.InletCold.T)),
				:(Phi = (sqrt(((_base_1.InletHot.T- _base_1.OutletHot.T)*(_base_1.InletHot.T- _base_1.OutletHot.T))+((_base_1.OutletCold.T - _base_1.InletCold.T)*(_base_1.OutletCold.T - _base_1.InletCold.T))))/(2*((_base_1.InletHot.T+ _base_1.OutletHot.T)-(_base_1.InletCold.T+ _base_1.OutletCold.T)))),
				:(Rho*(1-P*R) = (1-P)),
				:(Method.Fc = (4*Phi)/(ln(abs((1+2*Phi)/(1-2*Phi))))),
				:(Method.Fc = (2*Phi*(Rho+1)*ln(abs(Rho)))/( ln(abs((1+2*Phi)/(1-2*Phi)))*(Rho-1))),
			],
			[
				"Duty","Temperature Difference at Inlet","Temperature Difference at Outlet","R: Capacity Ratio for LMTD Correction Fator","P: Non - Dimensional Variable for LMTD Correction Fator"," Variable useless with this model"," Variable useless with this model","LMTD Correction Factor in Cocurrent Flow","Temperature Difference at Inlet","Temperature Difference at Outlet","R: Capacity Ratio for LMTD Correction Fator","P: Non - Dimensional Variable for LMTD Correction Fator"," Variable useless with this model"," Variable useless with this model","LMTD Correction Factor in Counter Flow","Temperature Difference at Inlet","Temperature Difference at Outlet"," Variable not in use with Bowmann equation"," Variable not in use with Bowmann equation","R: Capacity Ratio for LMTD Correction Fator when Shell and Tube","P: Non - Dimensional Variable for LMTD Correction Fator when Shell and Tube","LMTD Correction Fator when 1 Pass Shell Side","LMTD Correction Fator when 1 Pass Shell Side","R: Capacity Ratio for LMTD Correction Fator when Shell and Tube","P: Non - Dimensional Variable for LMTD Correction Fator when Shell and Tube","Non Dimensional Variable for LMTD Correction Fator in Fakheri Equation ","Non Dimensional Variable for LMTD Correction Fator in Fakheri Equation","LMTD Correction Fator when 1 Pass Shell Side","LMTD Correction Fator when 1 Pass Shell Side",
			],
			[:ExchangerType,:LMTDcorrection,],
			[:Method,:R,:P,:Rho,:Phi,]
		)
	end
	_base_1::Heatex_Basic
	ExchangerType::DanaSwitcher
	LMTDcorrection::DanaSwitcher
	Method::LMTD_Basic
	R::positive
	P::positive
	Rho::positive
	Phi::positive
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Heatex_LMTD
function setEquationFlow(in::Heatex_LMTD)
	addEquation(1)
	let switch=ExchangerType
		if switch=="Cocurrent Flow"
			addEquation(2)
			addEquation(3)
			addEquation(4)
			addEquation(5)
			addEquation(6)
			addEquation(7)
			addEquation(8)
		elseif switch=="Counter Flow"
			addEquation(9)
			addEquation(10)
			addEquation(11)
			addEquation(12)
			addEquation(13)
			addEquation(14)
			addEquation(15)
		elseif switch=="Shell and Tube"
			addEquation(16)
			addEquation(17)
			let switch=LMTDcorrection
				if switch=="Bowmann"
					addEquation(18)
					addEquation(19)
					addEquation(20)
					addEquation(21)
					if R == 1
						addEquation(22)
					else
						addEquation(23)
					end
				elseif switch=="Fakheri"
					addEquation(24)
					addEquation(25)
					addEquation(26)
					addEquation(27)
					if Rho == 1
						addEquation(28)
					else
						addEquation(29)
					end
				end
			end
		end
	end
end
function atributes(in::Heatex_LMTD,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/Heatex"
	fields[:Brief]="Simplified model for Heat Exchangers"
	fields[:Info]="This model perform material and heat balance using the Log Mean Temperature Difference Approach.
This shortcut calculation does not require exchanger configuration or geometry data.

== Assumptions ==
* Steady-State operation;
* No heat loss to the surroundings.

== Specify ==
* The Inlet streams: Hot and Cold.

== References ==
[1] E.A.D. Saunders, Heat Exchangers: Selection, Design and
 Construction, Longman, Harlow, 1988.

[2] Taborek, J., Shell-and-tube heat exchangers, in Heat Exchanger Design Handbook, Vol. 3
 Hemisphere Publishing Corp., New York, 1988.

[3] Fakheri, A. , Alternative approach for determining log mean temperature difference correction factor
 and number of shells of shell and tube heat exchangers, Journal of Enhanced Heat Transfer, v. 10, p. 407- 420, 2003.
"
	drive!(fields,_)
	return fields
end
Heatex_LMTD(_::Dict{Symbol,Any})=begin
	newModel=Heatex_LMTD()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using LogisticRegression
using Test

@testset "LogisticRegression.jl" begin
    include("logistic_regression.jl")
end


text:
export FFNNLayer, activate, update!

"""
`type FFNNLayer`

Type representing a Neural Network 1-D layer with `N` neurons and, eventually, a bias unit.

### Fields
* `neurons` (`Vector{Float64}`): Vector with each neuron and, if `bias = true`, a bias unit (index 1)
* `activation` (`Function`): Activation function
* `bias` (`Bool`): True if there is a bias unit in this layer
"""
type FFNNLayer
    neurons::Vector{Float64}
    activation::Function
    bias::Bool
end

#############################
#        CONSTRUCTORS       #
#############################

"""
`FFNNLayer(n::Integer; bias::Bool = true)`

Construct a 1-D layer of a Neural Network with `n` neurons and, eventually, a bias unit. The layer has `tanh` as activation function

### Arguments
* `n` (`Int`): Number of neurons in this layer (not counting the eventual bias unit)

### Keyword Arguments
* `bias` (`Bool`, `true` by default): True if there is a bias unit in this layer
"""
function FFNNLayer(n::Integer; bias::Bool = true)
    if bias
        return FFNNLayer(vcat([1.0], zeros(n)), tanh, bias)
    else
        return FFNNLayer(zeros(n), tanh, bias)
    end
end

"""
`FFNNLayer(n::Integer; bias::Bool = true)`

Construct a 1-D layer of a Neural Network with `n` neurons, `f` as activation function and, eventually, a bias unit.

### Arguments
* `n` (`Int`): Number of neurons in this layer (not counting the eventual bias unit)
* `f` (`Function`): Activation function of this layer

### Keyword Arguments
* `bias` (`Bool`, `true` by default): True if there is a bias unit in this layer
"""
function FFNNLayer(n::Integer, f::Function; bias::Bool = true)
    if bias
        return FFNNLayer(vcat([1.0], zeros(n)), f, bias)
    else
        return FFNNLayer(zeros(n), f, bias)
    end
end

############################
#      BASIC FUNCTIONS     #
############################
length(l::FFNNLayer) = l.bias ? length(l.neurons)-1 : length(l.neurons)
size(l::FFNNLayer)   = length(l)
endof(l::FFNNLayer)  = length(l)

function show(io::IO, l::FFNNLayer)
    N = length(l)
    print(io, "Feedforward Neural Net Layer: ", N, " neurons")
    if l.bias
        print(io, " + bias unit")
    end
    print(io, ", ", string(l.activation))
end

getindex(l::FFNNLayer, i) = l.bias ? l.neurons[i+1] : l.neurons[i]

function setindex!(l::FFNNLayer, x, i)
    if l.bias
        l.neurons[i+1] = x
    else
        l.neurons[i]   = x
    end
end

start(l::FFNNLayer) = 1
done(l::FFNNLayer, s) = s > length(l)
next(l::FFNNLayer, s) = (l[s], s+1)

function ==(l1::FFNNLayer, l2::FFNNLayer)
    (size(l1.neurons) == size(l2.neurons)) &&
    (l1.activation == l2.activation) &&
    (l1.bias == l2.bias)
end

############################
#      OTHER FUNCTIONS     #
############################
"Activation of a Neural Network layer"
activate(l::FFNNLayer) = l.activation(l)

"Update the internal values of the neurons of a layer "
update!(l::FFNNLayer, x::Vector{Float64}) = l[1:end] = x


text:
test_gage_principal_axes(
    [0. 1. 1.;
     0. 0. 0.;],
    [0.  1.  0.;
     sqrt(2)/2. 0. -sqrt(2)/2.;
     sqrt(2)/2. 0. sqrt(2)/2. ], "vector in the Y-Z plane")


text:
@testset "cl" begin
    for n in vcat(collect(1:16), [1000, 1001, 1_000_000])
        data = open(readdlm, joinpath(@__DIR__, "data", "Cl$(n).txt"))

        for r in 1:size(data, 1)
            row      = data[r, :]
            x        = row[1]
            expected = row[2]
            sgn      = -(-1)^n

            cl  = ClausenFunctions.cl(n, x)
            clm = sgn*ClausenFunctions.cl(n, -x)

            @test cl == clm
            @test cl ≈ expected atol=1e-14
        end

    end

    @test_throws DomainError ClausenFunctions.cl(0, 1.0)
    @test_throws DomainError ClausenFunctions.cl(-1, 1.0)
    @test_throws DomainError ClausenFunctions.cl(-2, 1.0)
end


text:
# In this simulation the synapse is the smallest component.
# A synapse is feed by a StreamMerger's output

const INITIAL_PRE_T = 0.0 # -1000000000.0

mutable struct Synapse <: AbstractSynapse
    soma::AbstractSoma
    dendrite::AbstractDendrite
    compartment::AbstractCompartment
   	model::Model.ModelData

    id::Int64

    # true = excititory, false = inhibitory
    excititory::Bool

    # The weight of the synapse
    w::Float64

    wMax::Float64
   	wMin::Float64

    # The stream (aka Merger) that feeds into this synapse
    stream::AbstractBitStream

	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	# Surge
	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	# Surge base value
   	amb::Float64

	# Surge peak
   	ama::Float64

	# Surge window
   	tsw::Float64

	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	# new surge ion concentration
	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	# concentration base. We should always have a minimum concentration
	# as a result of a spike
	# Surge is calculated at the arrival of a spike
	# surge = amb - ama*e^(-psp/tsw) == rising curve
   	surge::Float64

	# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

	# The time-mark at which a spike arrived at a synapse
   	preT::Float64

	# The current ion concentration
    psp::Float64

	# =============================================================
	# Learning rules:
	# =============================================================
	#
	# Depression pair-STDP, Potentiation is triplet.
	# "tao"s control the rate of decay. Larger values means a slower decay.
	# Smaller values equals a sharper decay.
	# -----------------------------------

	# denominator, positive window time decay
   	taoP::Float64

	# denominator, negative window time decay
   	taoN::Float64

	# Ratio of mRate/taoX
   	tao::Float64

	# -----------------------------------
	# Weight dependence
	# -----------------------------------
	# F-(w) = λ⍺w^µ, F+(w) = λ(1-w)^µ
   	mu::Float64 # µ
   	lambda::Float64 # λ
   	alpha::Float64 # ⍺

	# -----------------------------------
	# Suppression
	# -----------------------------------
   	taoI::Float64
   	prevEffTrace::Float64

   	learningRateSlow::Float64 # Unused
   	learningRateFast::Float64 # Unused

	# -----------------------------------
	# Fall off
	# -----------------------------------
	distanceEfficacy::Float64

	# The farther the synapse is from the Soma the less of an influence
	# this synapse has. The function can either be linear or non-linear.
	# The default is linear.
	# Note: no matter how far a synapse is it will still have an influence
	# otherwise it's useless.
	# distance's value = 1.0 for synapses closest to soma. The farther out
	# the value reaches a minimum of around ~0.25.
	# distance is multiplied into soma.psp.
	distance::Float64

    function Synapse(soma::AbstractSoma, dendrite::AbstractDendrite, compartment::AbstractCompartment, model::Model.ModelData)
        o = new()
        o.model = model
        o.soma = soma
        o.dendrite = dendrite
        o.compartment = compartment
        o.excititory = true # default to excite type
        o.preT = INITIAL_PRE_T
        o.id = 0
        add_synapse!(compartment, o)

        o
    end
end

function set_stream!(syn::Synapse, stream::AbstractBitStream)
   	syn.stream = stream
end

function set_as_inhibit!(syn::Synapse)
    syn.excititory = false
end

function initialize!(syn::Synapse)
	# Focus the model on the correct synapse.
   	Model.set_active_synapse!(syn.model, syn.id)

	# Set properties based on model. These drive the other properties.
   	syn.taoP = Model.taoP(syn.model)
   	syn.taoN = Model.taoN(syn.model)
   	syn.mu = Model.mu(syn.model)
   	syn.distance = Model.distance(syn.model)
   	syn.lambda = Model.lambda(syn.model)
   	syn.amb = Model.amb(syn.model)
   	syn.w = Model.weight(syn.model)
   	syn.alpha = Model.alpha(syn.model)
   	syn.learningRateFast = Model.learning_rate_fast(syn.model)
   	syn.learningRateSlow = Model.learning_rate_slow(syn.model)
   	syn.taoI = Model.taoI(syn.model)
   	syn.ama = Model.ama(syn.model)

	# Calc this synapses's reaction to the AP based on its
	# distance from the soma.
   	syn.distanceEfficacy = AP_efficacy(syn.compartment.dendrite, syn.distance)

   	syn.wMax = Model.w_max(syn.model)
   	syn.wMin = Model.w_min(syn.model)

   	# println("___ Synapse properties ___")
   	# println("| taoI: ", syn.taoI)
   	# println("| taoP: ", syn.taoP)
   	# println("| taoN: ", syn.taoN)
   	# println("| mu: ", syn.mu)
   	# println("| distance: ", syn.distance)
   	# println("| lambda: ", syn.lambda)
   	# println("| ama: ", syn.ama)
   	# println("| amb: ", syn.amb)
   	# println("| w: ", syn.w)
   	# println("| alpha: ", syn.alpha)
   	# println("| learningRateFast: ", syn.learningRateFast)
   	# println("| learningRateSlow: ", syn.learningRateSlow)
   	# println("| distanceEfficacy: ", syn.distanceEfficacy)
   	# println("| wMax: ", syn.wMax)
   	# println("| wMin: ", syn.wMin)
   	# println("---------------------------")

end

# Reset properties to default values from model.
function reset!(syn::Synapse)
   	syn.prevEffTrace = 1.0
   	syn.surge = 0.0
   	syn.psp = 0.0
    syn.preT = 0.0

	# Reset weights back to best guess values.
    syn.wMax = syn.compartment.weight_max
   	syn.w = syn.wMax / syn.compartment.weight_divisor
end

# This is considered the 1st pass of the simulation per time step.
# Each "pass" is broken up into several steps.
# Internal values are 'moved' to the outputs.
# Learning rules are applied.
# Prepare for next integration step
function pre_integrate!(syn::Synapse)
end

# The main integration. Returns a float value
function integrate!(syn::Synapse, span_t::Int64, t::Int64)
    triplet_integration(syn, span_t, t)
end

function post_integrate!(syn::Synapse)
end

# =============================================================
# Triplet:
# =============================================================
# Pre trace, Post slow and fast traces.
#
# Depression: fast post trace with at pre spike
# Potentiation: slow post trace at post spike
function triplet_integration(syn::Synapse, span_t::Int64, t::Int64)
	# Calc psp based on current dynamics: (t - preT). As dt increases
	# psp will decrease asymtotically to zero.
    dt = Float64(t) - syn.preT

   	dwD = 0.0
   	dwP = 0.0
   	updateWeight = false

	# The output of the stream is the input to this synapse.
	# The stream is almost always a StreamMerger
	   if output(syn.stream) == 1
		# A spike has arrived on the input to this synapse.
       	# println("(", t, ") syn: ", syn.id)
       	syn.surge = if syn.excititory
            syn.psp + syn.ama * exp(-syn.psp / syn.taoP)
       	else
            syn.psp + syn.ama * exp(-syn.psp / syn.taoN)
        end

		# #######################################
		# Depression LTD
		# #######################################
		# Read post trace and adjust weight accordingly.
        dwD = syn.prevEffTrace * weight_factor(syn, false) * syn.soma.apFast

        syn.prevEffTrace = efficacy(syn, dt)

        syn.preT = t
        dt = 0.0

        updateWeight = true
    end

   	syn.psp = if syn.excititory
       	syn.surge * exp(-dt / syn.taoP)
   	else
       	syn.surge * exp(-dt / syn.taoN)
   	end

	# If an AP occurred (from the soma) we read the current psp value and add it to the "w"
   	if output(syn.soma) == 1.0
		# #######################################
		# Potentiation LTP
		# #######################################
		# Read pre trace (aka psp) and slow AP trace for adjusting weight accordingly.
		#     Post efficacy                                          weight dependence                 triplet sum
        dwP = efficacy_trace(syn.soma) * syn.distanceEfficacy * weight_factor(syn, true) * (syn.psp + syn.soma.apSlowPrior)
        updateWeight = true
    end

	# Finally update the weight.
   	if updateWeight
        syn.w = max(min(syn.w + dwP - dwD, syn.wMax), syn.wMin)
    end

	# Return the "value" of this synapse for this "t"
   	value = if syn.excititory
       	syn.psp * syn.w * syn.distance
   	else
       	-syn.psp * syn.w * syn.distance # is inhibitory
   	end

	# Collecting is centralized in streams.jl for consistency.
   	collect_synapse!(syn.soma.samples, syn, span_t)

   	(value, syn.w)
end

# Each spike of pre-synaptic neuron j sets the presynaptic spike
# efficacy j to 0
# whereafter it recovers exponentially to 1 with a time constant
# τj = toaJ
# In other words, the efficacy of a spike is suppressed by
# the proximity of a previous spike.
function efficacy(syn::Synapse, dt)
   	1.0 - exp(-dt / syn.taoI)
end

# mu = 0.0 = additive, mu = 1.0 = multiplicative
function weight_factor(syn::Synapse, potentiation::Bool)
    if potentiation
        return syn.lambda * ((1.0 - syn.w / syn.wMax)^syn.mu)
   	end

   	syn.lambda * syn.alpha * ((syn.w / syn.wMax)^syn.mu)
end


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
Pkg.instantiate()

# ## Setup
using Dojo
using Random
using LinearAlgebra
using JLD2
include(joinpath(@__DIR__, "algorithms/ars.jl")) # augmented random search

# # Training

# ## Environment
env = get_environment(:halfcheetah, timestep=0.05)
obs = reset(env);

# ## Random policy
Random.seed!(38)
hp = HyperParameters(
    main_loop_size=10,
    horizon=80,
    n_directions=6,
    b=6,
    step_size=0.02)
input_size = length(obs)
output_size = length(env.input_previous)
policy = Policy(input_size, output_size, hp)
normalizer = Normalizer(input_size);

# ## Train policy
train_time = @elapsed train(env, policy, normalizer, hp)

# ## Evaluate policy
reward = rollout_policy(policy.θ, env, normalizer, hp, reset=env->reset(env,reset_noise_scale=0.0))

# @save joinpath(@__DIR__, "results/ars_demonstration_dynamic.jld2") train_time reward policy normalizer
@load joinpath(@__DIR__, "results/ars_demonstration_dynamic.jld2") train_time reward policy normalizer

# ## rollout trained policy
traj = display_policy(env,
    policy,
    normalizer,
    hp,
    reset=env->reset(env,reset_noise_scale=0.0))

for t = 1:length(traj)
    traj[t][2] += -0.50
end

# ## Open Visualizer
open(env.vis)

# ## Visualize Policy
visualize(env, traj)
set_camera!(env.vis,
    cam_pos=[0.0, -35.0, 0.0],
    zoom=15)

# set_camera!(env.vis,
#     cam_pos=[3.0, -6.0, 2.0],
#     zoom=2)
# render_static(env.vis)
# open(joinpath(@__DIR__, "results/ars_demonstration_dynamic.html"), "w") do file
#     write(file, static_html(env.vis))
# end
include("../generate_notebooks.jl")


text:
# TO-DO :
# -tenir compte du nombre de parametres pour les initialvalues dans la fct fit
# -vérifier la convergence pour l'optimisation dans la fct fit

"""
    nonstatEGPpower(σ, ξ₀, ξ₁, κ₀, κ₁)

Extension non stationnaire du modèle 1 de Naveau. À détailler...
"""
struct nonstatEGPpower{T<:Real} <: ContinuousUnivariateDistribution
    σ₀::T    # scale parameter
    σ₁::T    # scale parameter
    ξ₀::T    # rate of upper tail decay
    ξ₁::T    # rate of upper tail decay
    κ₀::T    # shape of the lower tail
    κ₁::T    # shape of the lower tail

    function nonstatEGPpower{T}(σ₀::T, σ₁::T, ξ₀::T, ξ₁::T, κ₀::T, κ₁::T) where {T <: Real}
        new{T}(σ₀, σ₁, ξ₀, ξ₁, κ₀, κ₁)
    end
end

function nonstatEGPpower(σ₀::T, σ₁::T, ξ₀::T, ξ₁::T, κ₀::T, κ₁::T; check_args=true) where {T <: Real}
    #check_args && @check_args(EGPpower, σ > zero(σ) && κ > zero(κ))
    return nonstatEGPpower{T}(σ₀, σ₁, ξ₀, ξ₁, κ₀, κ₁)
end

nonstatEGPpower(σ₀::Real, σ₁::Real, ξ₀::Real, ξ₁::Real, κ₀::Real, κ₁::Real) = nonstatEGPpower(promote(σ₀, σ₁, ξ₀, ξ₁, κ₀, κ₁)...)
nonstatEGPpower(σ₀::Integer, σ₁::Integer, ξ₀::Integer, ξ₁::Integer, κ₀::Integer, κ₁::Integer) = nonstatEGPpower(float(σ₀), float(σ₁), float(ξ₀), float(ξ₁), float(κ₀), float(κ₁))

params(d::nonstatEGPpower) = (d.σ₀, d.σ₁, d.ξ₀, d.ξ₁, d.κ₀, d.κ₁)

"""
    EGPpower(fm::nonstatEGPpower{T}, covariate::Real)

Pour obtenir le modèle stationnaire à partir de la valeur de la variable explicative.
"""
function EGPpower(fm::nonstatEGPpower{T}, covariate::Real) where T<:Real
    σ = fm.σ₀ + (fm.σ₁ * covariate)
    ξ = fm.ξ₀ + (fm.ξ₁ * covariate)
    κ = fm.κ₀ + (fm.κ₁ * covariate)
    return EGPpower(σ, ξ, κ)
end

#### Evaluation

function logpdf(d::nonstatEGPpower{T}, x::Real, covariate::Real) where T<:Real
    #TODO : assert σ not equals to zero
    μ = 0
    σ = d.σ₀ + d.σ₁*covariate
    ξ = d.ξ₀ + d.ξ₁*covariate
    κ = d.κ₀ + d.κ₁*covariate

    pd = GeneralizedPareto(μ, 1, ξ)

    lg(v::Real) = log(κ) + (κ-1)*log(v)

    p = -log(σ) + Distributions.logpdf(pd, x/σ) + lg(Distributions.cdf(pd, x/σ))

    return p
end

pdf(d::nonstatEGPpower, x::Real, covariate::Real) = exp(logpdf(d, x, covariate))

function logcdf(d::nonstatEGPpower{T}, x::Real, covariate::Real) where T<:Real
    #TODO : assert σ not equals to zero
    μ = 0
    σ = d.σ₀ + d.σ₁*covariate
    ξ = d.ξ₀ + d.ξ₁*covariate
    κ = d.κ₀ + d.κ₁*covariate

    pd = GeneralizedPareto(μ, 1, ξ)

    #G(v::Real, κ::Real) = v^κ
    #lG(v::Real, κ::Real) = κ*log(v)

    p = κ*Distributions.logcdf(pd, x/σ)

    return p
end

cdf(d::nonstatEGPpower, x::Real, covariate::Real) = exp(logcdf(d, x, covariate))

function quantile(d::nonstatEGPpower{T}, p::Real, covariate::Real) where T<:Real
    @assert zero(p)<p<one(p) "the quantile level should be between 0 and 1."

    σ = d.σ₀ + d.σ₁*covariate
    ξ = d.ξ₀ + d.ξ₁*covariate
    κ = d.κ₀ + d.κ₁*covariate

    invG(p::Real, κ::Real) = p^(1/κ)

    x = (σ/ξ)*(((1 - invG(p,κ))^(-ξ)) - 1)

    return x
end

#### Parameter estimation

function EGPnonstatpowerfit(data::Array{<:Real,1},
    covariate::Array{<:Real,1};  # Vecteur de la variable explicative
    initialvalues::Vector{<:Real}=Float64[],
    scalecov::Bool=true, # variable explicative sur σ ?
    lowertailcov::Bool=false,  # variable explicative sur κ ?
    uppertailcov::Bool=false,  # variable explicative sur ξ ?
    censoring::Real=0)

    # TO-DO :
    # -tenir compte du nombre de parametres pour les initialvalues
    # -voir pour ajouter la censure

    if isempty(initialvalues)
        σ₀₀, σ₁₀, ξ₀₀, ξ₁₀, κ₀₀, κ₁₀ = [1, 1, 0.15, 0, 1, 1]
    else
        σ₀₀, σ₁₀, ξ₀₀, ξ₁₀, κ₀₀, κ₁₀ = initialvalues
    end

    if !lowertailcov & !uppertailcov & !scalecov
        return EGPpowerfit(data, initialvalues = [σ₀₀, ξ₀₀, κ₀₀], censoring = censoring)
    else
        r_data = data[data .>= censoring]
        l_data = fill(censoring, count(data .< censoring))

        r_cov = covariate[data .>= censoring]
        l_cov = covariate[data .< censoring]

        # Model 1 : covariate on σ only
        function loglike_sigma(σ₀, σ₁, ξ₀, κ₀)
            if σ₀*σ₁ <= 0 || κ₀ <= 0
                return -Inf
            else
                pd = nonstatEGPpower(σ₀, σ₁, ξ₀, 0, κ₀, 0)
                #ll = sum(logpdf.(pd, data, covariate))
                ll = sum(logcdf.(pd, l_data, l_cov)) + sum(logpdf.(pd, r_data, r_cov))
                return ll
            end
        end
        # Model 2 : covariate on σ and ξ
        function loglike_sigma_xi(σ₀, σ₁, ξ₀, ξ₁, κ₀)
            if σ₀*σ₁ <= 0 || κ₀ <= 0
                return -Inf
            else
                pd = nonstatEGPpower(σ₀, σ₁, ξ₀, ξ₁, κ₀, 0)
                #ll = sum(logpdf.(pd, data, covariate))
                ll = sum(logcdf.(pd, l_data, l_cov)) + sum(logpdf.(pd, r_data, r_cov))
                return ll
            end
        end
        # Model 3 : covariate on σ, ξ and κ
        function loglike_all(σ₀, σ₁, ξ₀, ξ₁, κ₀, κ₁)
            if σ₀*σ₁ <= 0 || κ₀*κ₁ <= 0
                return -Inf
            else
                pd = nonstatEGPpower(σ₀, σ₁, ξ₀, ξ₁, κ₀, κ₁)
                #ll = sum(logpdf.(pd, data, covariate))
                ll = sum(logcdf.(pd, l_data, l_cov)) + sum(logpdf.(pd, r_data, r_cov))
                return ll
            end
        end
        # Model 4 : covariate on σ and κ
        function loglike_sigma_kappa(σ₀, σ₁, ξ₀, κ₀, κ₁)
            if σ₀*σ₁ <= 0 || κ₀*κ₁ <= 0
                return -Inf
            else
                pd = nonstatEGPpower(σ₀, σ₁, ξ₀, 0, κ₀, κ₁)
                #ll = sum(logpdf.(pd, data, covariate))
                ll = sum(logcdf.(pd, l_data, l_cov)) + sum(logpdf.(pd, r_data, r_cov))
                return ll
            end
        end
        # Model 5 : covariate on κ only
        function loglike_kappa(σ₀, ξ₀, κ₀, κ₁)
            if σ₀ <= 0 || κ₀*κ₁ <= 0
                return -Inf
            else
                pd = nonstatEGPpower(σ₀, 0, ξ₀, 0, κ₀, κ₁)
                #ll = sum(logpdf.(pd, data, covariate))
                ll = sum(logcdf.(pd, l_data, l_cov)) + sum(logpdf.(pd, r_data, r_cov))
                return ll
            end
        end
        # Model 6 : covariate on ξ and κ
        function loglike_xi_kappa(σ₀, ξ₀, ξ₁, κ₀, κ₁)
            if σ₀ <= 0 || κ₀*κ₁ <= 0
                return -Inf
            else
                pd = nonstatEGPpower(σ₀, 0, ξ₀, ξ₁, κ₀, κ₁)
                #ll = sum(logpdf.(pd, data, covariate))
                ll = sum(logcdf.(pd, l_data, l_cov)) + sum(logpdf.(pd, r_data, r_cov))
                return ll
            end
        end
        # Model 7 : covariate on ξ only
        function loglike_xi(σ₀, ξ₀, ξ₁, κ₀)
            if σ₀ <= 0 || κ₀ <= 0
                return -Inf
            else
                pd = nonstatEGPpower(σ₀, 0, ξ₀, ξ₁, κ₀, 0)
                #ll = sum(logpdf.(pd, data, covariate))
                ll = sum(logcdf.(pd, l_data, l_cov)) + sum(logpdf.(pd, r_data, r_cov))
                return ll
            end
        end

        # Fonctions objectives
        fobj_1(θ) = -loglike_sigma(θ...)
        fobj_2(θ) = -loglike_sigma_xi(θ...)
        fobj_3(θ) = -loglike_all(θ...)
        fobj_4(θ) = -loglike_sigma_kappa(θ...)
        fobj_5(θ) = -loglike_kappa(θ...)
        fobj_6(θ) = -loglike_xi_kappa(θ...)
        fobj_7(θ) = -loglike_xi(θ...)

        if !lowertailcov & !uppertailcov & scalecov
            # Model 1 : covariate on σ only
            res = optimize(fobj_1, [σ₀₀, σ₁₀ , ξ₀₀, κ₀₀])
            fittedmodel = nonstatEGPpower(Optim.minimizer(res)[1], Optim.minimizer(res)[2], Optim.minimizer(res)[3], 0, Optim.minimizer(res)[4], 0)

        elseif !lowertailcov & uppertailcov & scalecov
            # Model 2 : covariate on σ and ξ
            res = optimize(fobj_2, [σ₀₀, σ₁₀ , ξ₀₀, ξ₁₀, κ₀₀])
            fittedmodel = nonstatEGPpower(Optim.minimizer(res)[1], Optim.minimizer(res)[2], Optim.minimizer(res)[3], Optim.minimizer(res)[4], Optim.minimizer(res)[5], 0)

        elseif lowertailcov & uppertailcov & scalecov
            # Model 3 : covariate on σ, ξ and κ
            res = optimize(fobj_3, [σ₀₀, σ₁₀, ξ₀₀, ξ₁₀, κ₀₀, κ₁₀])
            fittedmodel = nonstatEGPpower(Optim.minimizer(res)[1], Optim.minimizer(res)[2], Optim.minimizer(res)[3], Optim.minimizer(res)[4], Optim.minimizer(res)[5], Optim.minimizer(res)[6])

        elseif lowertailcov & !uppertailcov & scalecov
            # Model 4 : covariate on σ and κ
            res = optimize(fobj_4, [σ₀₀, σ₁₀, ξ₀₀, κ₀₀, κ₁₀])
            fittedmodel = nonstatEGPpower(Optim.minimizer(res)[1], Optim.minimizer(res)[2], Optim.minimizer(res)[3], 0, Optim.minimizer(res)[4], Optim.minimizer(res)[5])

        elseif lowertailcov & !uppertailcov & !scalecov
            # Model 5 : covariate on κ only
            res = optimize(fobj_5, [σ₀₀, ξ₀₀, κ₀₀, κ₁₀])
            fittedmodel = nonstatEGPpower(Optim.minimizer(res)[1], 0, Optim.minimizer(res)[2], 0, Optim.minimizer(res)[3], Optim.minimizer(res)[4])

        elseif lowertailcov & uppertailcov & !scalecov
            # Model 6 : covariate on ξ and κ
            res = optimize(fobj_6, [σ₀₀, ξ₀₀, ξ₁₀, κ₀₀, κ₁₀])
            fittedmodel = nonstatEGPpower(Optim.minimizer(res)[1], 0, Optim.minimizer(res)[2], Optim.minimizer(res)[3], Optim.minimizer(res)[4], Optim.minimizer(res)[5])

        elseif !lowertailcov & uppertailcov & !scalecov
            # Model 7 : covariate on ξ only
            res = optimize(fobj_7, [σ₀₀, ξ₀₀, ξ₁₀, κ₀₀])
            fittedmodel = nonstatEGPpower(Optim.minimizer(res)[1], 0, Optim.minimizer(res)[2], Optim.minimizer(res)[3], Optim.minimizer(res)[4], 0)
        end

        return fittedmodel, EGPpower.(fittedmodel, covariate)

        # TO-DO : - vérifier la convergence, ex :
            #if Optim.converged(res)
            #    σ̂, ξ̂, κ̂ = [Optim.minimizer(res)[1], Optim.minimizer(res)[2], Optim.minimizer(res)[3]]
            #else
            #    @warn "The maximum likelihood algorithm did not find a solution. Maybe try with different initial values or with another method. The returned values are the initial values."
            #    σ̂, ξ̂, κ̂  = [initialvalues[1], initialvalues[2], initialvalues[3]]
            #end
    end
end

"""
    BIC(fm, data, covariate)

"""
function BIC(fm::nonstatEGPpower{T}, data::Array{T,1}, covariate::Array{T,1}; censoring::Real=0) where T<:Real
    n = size(data,1)
    k = size(params(fm),1)
    fm.σ₁ == 0 ? k = k-1 : nothing
    fm.ξ₁ == 0 ? k = k-1 : nothing
    fm.κ₁ == 0 ? k = k-1 : nothing

    r_data = data[data .>= censoring]
    l_data = fill(censoring, count(data .< censoring))

    r_cov = covariate[data .>= censoring]
    l_cov = covariate[data .< censoring]

    ll = sum(logcdf.(fm, l_data, l_cov)) + sum(logpdf.(fm, r_data, r_cov))

    return ll - k/2*log(n)
end


text:
using Statistics
using StatsBase:sample
using Base.Threads:@threads
using BenchmarkTools
using EvoTrees
using SIMD
using LoopVectorization

n = Int(1e6)
nvars = 100
nbins = 64
𝑖 = collect(1:n)
𝑗 = collect(1:nvars)
δ = rand(n)
δ² = rand(n)
𝑤 = rand(n)

hist_δ = zeros(nbins, nvars)
hist_δ² = zeros(nbins, nvars)
hist_𝑤 = zeros(nbins, nvars)
X_bin = reshape(sample(UInt8.(1:nbins), n * nvars), n, nvars)

function iter_1(X_bin, hist_δ, δ, 𝑖, 𝑗)
    # hist_δ .*= 0.0
    @inbounds @threads for j in 𝑗
        @inbounds for i in 𝑖
            hist_δ[X_bin[i,j], j] += δ[i]
        end
    end
end

@time iter_1(X_bin, hist_δ, δ, 𝑖, 𝑗)
@btime iter_1(X_bin, hist_δ, δ, 𝑖, 𝑗)

function iter_1A(X_bin, hist_δ, δ, 𝑖, 𝑗)
    # hist_δ .*= 0.0
    @inbounds @threads for j in 𝑗
        @inbounds for i in 𝑖
            hist_δ[X_bin[i,j], j] += δ[i]
        end
    end
end

@time iter_1A(X_bin, hist_δ, δ, 𝑖, 𝑗)
@btime iter_1A(X_bin, hist_δ, δ, 𝑖, 𝑗)

### 3 features in seperate hists
function iter_1B(X_bin, hist_δ, hist_δ², hist_𝑤, δ, δ², 𝑤, 𝑖, 𝑗)
    # hist_δ .= 0.0
    # hist_δ² .= 0.0
    # hist_𝑤 .= 0.0
    @inbounds @threads for j in 𝑗
        @inbounds for i in 𝑖
            hist_δ[X_bin[i,j], j] += δ[i]
            hist_δ²[X_bin[i,j], j] += δ²[i]
            hist_𝑤[X_bin[i,j], j] += 𝑤[i]
        end
    end
end

@time iter_1B(X_bin, hist_δ, hist_δ², hist_𝑤, δ, δ², 𝑤, 𝑖, 𝑗)
@btime iter_1B($X_bin, $hist_δ, $hist_δ², $hist_𝑤, $δ, $δ², $𝑤, $𝑖, $𝑗)

𝑖2 = sample(𝑖, 500000, replace=false, ordered=true)
𝑗2 = sample(𝑗, 50, replace=false, ordered=true)
@time iter_1B(X_bin, hist_δ, hist_δ², hist_𝑤, δ, δ², 𝑤, 𝑖2, 𝑗2)
@btime iter_1B($X_bin, $hist_δ, $hist_δ², $hist_𝑤, $δ, $δ², $𝑤, $𝑖2, $𝑗2)

### 3 features in common hists
hist_δ𝑤 = zeros(3, nbins, nvars)
function iter_2(X_bin, hist_δ𝑤, δ, δ², 𝑤, 𝑖, 𝑗)
    # hist_δ𝑤 .= 0.0
    @inbounds @threads for j in 𝑗
        @inbounds @simd for i in 𝑖
            hist_δ𝑤[1, X_bin[i,j], j] += δ[i]
            hist_δ𝑤[2, X_bin[i,j], j] += δ²[i]
            hist_δ𝑤[3, X_bin[i,j], j] += 𝑤[i]
        end
    end
end

@time iter_2(X_bin, hist_δ𝑤, δ, δ², 𝑤, 𝑖, 𝑗)
@btime iter_2($X_bin, $hist_δ𝑤, $δ, $δ², $𝑤, $𝑖, $𝑗)

### 3 features in common hists - gradients/weight in single matrix
hist_δ𝑤 = zeros(3, nbins, nvars)
δ𝑤 = rand(3, n)

function iter_3(X_bin, hist_δ𝑤, δ𝑤, 𝑖, 𝑗)
    # hist_δ𝑤 .= 0.0
    @inbounds @threads for j in 𝑗
        @inbounds for i in 𝑖
            hist_δ𝑤[1, X_bin[i,j], j] += δ𝑤[1,i]
            hist_δ𝑤[2, X_bin[i,j], j] += δ𝑤[2,i]
            hist_δ𝑤[3, X_bin[i,j], j] += δ𝑤[3,i]
        end
    end
end

@time iter_3(X_bin, hist_δ𝑤, δ𝑤, 𝑖, 𝑗)
@btime iter_3($X_bin, $hist_δ𝑤, $δ𝑤, $𝑖, $𝑗)

# 39.849 ms (81 allocations: 8.22 KiB)
function iter_3B(X_bin, hist_δ𝑤, δ𝑤, 𝑖, 𝑗)
    @inbounds @threads for j in 𝑗
        @inbounds @simd for i in 𝑖
            hist_δ𝑤[1, X_bin[i,j], j] += δ𝑤[1,i]
            hist_δ𝑤[2, X_bin[i,j], j] += δ𝑤[2,i]
            hist_δ𝑤[3, X_bin[i,j], j] += δ𝑤[3,i]
        end
    end
end

@time iter_3B(X_bin, hist_δ𝑤, δ𝑤, 𝑖, 𝑗)
@btime iter_3B($X_bin, $hist_δ𝑤, $δ𝑤, $𝑖, $𝑗)

# 38.027 ms (81 allocations: 8.22 KiB)
function iter_3C(X_bin, hist_δ𝑤, δ𝑤, 𝑖, 𝑗)
    @inbounds @threads for j in 𝑗
        @inbounds for chunk in 1:10000
            @inbounds for iter in 1:100
                i = iter + 100 * (chunk - 1)
                hist_δ𝑤[1, X_bin[i,j], j] += δ𝑤[1,i]
                hist_δ𝑤[2, X_bin[i,j], j] += δ𝑤[2,i]
                hist_δ𝑤[3, X_bin[i,j], j] += δ𝑤[3,i]
            end
        end
    end
end

@time iter_3C(X_bin, hist_δ𝑤, δ𝑤, 𝑖, 𝑗)
@btime iter_3C($X_bin, $hist_δ𝑤, $δ𝑤, $𝑖, $𝑗)

### 3 features in common hists - vector of matrix hists - gradients/weight in single matrix
hist_δ𝑤_vec = [zeros(3, nbins) for j in 1:nvars]
δ𝑤 = rand(3, n)

# 33.762 ms (81 allocations: 8.22 KiB)
function iter_4(X_bin, hist_δ𝑤_vec, δ𝑤, 𝑖, 𝑗)
    # [hist_δ𝑤_vec[j] .= 0.0 for j in 𝑗]
    @inbounds @threads for j in 𝑗
        @inbounds @simd for i in 𝑖
            hist_δ𝑤_vec[j][1, X_bin[i,j]] += δ𝑤[1,i]
            hist_δ𝑤_vec[j][2, X_bin[i,j]] += δ𝑤[2,i]
            hist_δ𝑤_vec[j][3, X_bin[i,j]] += δ𝑤[3,i]
        end
    end
end

@time iter_4(X_bin, hist_δ𝑤_vec, δ𝑤, 𝑖, 𝑗)
@btime iter_4($X_bin, $hist_δ𝑤_vec, $δ𝑤, $𝑖, $𝑗)

### 3 features in common hists - vector of matrix hists - gradients/weight in single matrix
hist_δ𝑤_vec = [zeros(3, nbins) for j in 1:nvars]
δ𝑤 = rand(n, 3)

# 40.005 ms (81 allocations: 8.22 KiB)
function iter_4B(X_bin, hist_δ𝑤_vec, δ𝑤, 𝑖, 𝑗)
    # [hist_δ𝑤_vec[j] .= 0.0 for j in 𝑗]
    @inbounds @threads for j in 𝑗
        @inbounds @simd for i in 𝑖
            hist_δ𝑤_vec[j][1, X_bin[i,j]] += δ𝑤[i,1]
            hist_δ𝑤_vec[j][2, X_bin[i,j]] += δ𝑤[i,2]
            hist_δ𝑤_vec[j][3, X_bin[i,j]] += δ𝑤[i,3]
        end
    end
end

@time iter_4B(X_bin, hist_δ𝑤_vec, δ𝑤, 𝑖, 𝑗);
@btime iter_4B($X_bin, $hist_δ𝑤_vec, $δ𝑤, $𝑖, $𝑗);

### 3 features in common hists - vector of vec hists - gradients/weight in single vector
hist_δ𝑤_vec = [zeros(3 * nbins) for j in 1:nvars]
δ𝑤 = rand(3 * n)

function iter_5(X_bin, hist_δ𝑤_vec, δ𝑤, 𝑖, 𝑗)
    # [hist_δ𝑤_vec[j] .= 0.0 for j in 𝑗]
    @inbounds @threads for j in 𝑗
        @inbounds @simd for i in 𝑖
            # @inbounds if mask[i]
            id = 3 * i - 2
            hid = 3 * X_bin[i,j] - 2
            hist_δ𝑤_vec[j][hid] += δ𝑤[id]
            hist_δ𝑤_vec[j][hid + 1] += δ𝑤[id + 1]
            hist_δ𝑤_vec[j][hid + 2] += δ𝑤[id + 2]
            # end
        end
    end
end

@time iter_5(X_bin, hist_δ𝑤_vec, δ𝑤, 𝑖, 𝑗)
@btime iter_5($X_bin, $hist_δ𝑤_vec, $δ𝑤, $𝑖, $𝑗)

# 𝑖2 = sample(𝑖, 900000, replace=false, ordered=true)
# 𝑖3 = view(𝑖2, 100001:650000)
using Random
𝑖2 = sample(𝑖, 500000, replace=false, ordered=true)
𝑗2 = sample(𝑗, 50, replace=false, ordered=true)
@time iter_5(X_bin, hist_δ𝑤_vec, δ𝑤, 𝑖2, 𝑗2)
@btime iter_5($X_bin, $hist_δ𝑤_vec, $δ𝑤, $𝑖2, $𝑗2)

### 3 features in common hists - vector of vec hists - gradients/weight in single vector - explicit loop
hist_δ𝑤_vec = [zeros(3 * nbins) for j in 1:nvars]
δ𝑤 = rand(3, n)

function iter_5B(X_bin, hist_δ𝑤_vec, δ𝑤, 𝑖, 𝑗)
    @inbounds @threads for j in 𝑗
        @inbounds @simd for i in 𝑖
            hid = 3 * X_bin[i,j] - 2
            hist_δ𝑤_vec[j][hid] += δ𝑤[1, i]
            hist_δ𝑤_vec[j][hid + 1] += δ𝑤[2, i]
            hist_δ𝑤_vec[j][hid + 2] += δ𝑤[3, i]
        end
    end
end

@time iter_5B(X_bin, hist_δ𝑤_vec, δ𝑤, 𝑖, 𝑗)
@btime iter_5B($X_bin, $hist_δ𝑤_vec, $δ𝑤, $𝑖, $𝑗)

### 3 features in common hists - vector of vec hists - gradients/weight in single vector - explicit loop
hist_δ𝑤_vec = [zeros(3 * nbins) for j in 1:nvars]
δ𝑤 = rand(n, 3)

function iter_5C(X_bin, hist_δ𝑤_vec, δ𝑤, 𝑖, 𝑗)
    @inbounds @threads for j in 𝑗
        @inbounds @simd for i in 𝑖
            hid = 3 * X_bin[i,j] - 2
            hist_δ𝑤_vec[j][hid] += δ𝑤[i,1]
            hist_δ𝑤_vec[j][hid + 1] += δ𝑤[i,2]
            hist_δ𝑤_vec[j][hid + 2] += δ𝑤[i,3]
        end
    end
end

@time iter_5C(X_bin, hist_δ𝑤_vec, δ𝑤, 𝑖, 𝑗)
@btime iter_5C($X_bin, $hist_δ𝑤_vec, $δ𝑤, $𝑖, $𝑗)

function iter_6(X_bin, hist_δ𝑤_vec, δ𝑤, 𝑖, 𝑗, K)
    # [hist_δ𝑤_vec[j] .= 0.0 for j in 𝑗]
    @inbounds @threads for j in 𝑗
        @inbounds @simd for i in 𝑖
            id = 3 * i - 2
            hid = 3 * X_bin[i,j] - 3
            for k in 1:K
                hist_δ𝑤_vec[j][hid + k] += δ𝑤[id + k]
            end
        end
    end
end

K = 3
@time iter_6(X_bin, hist_δ𝑤_vec, δ𝑤, 𝑖, 𝑗, K)
@btime iter_6($X_bin, $hist_δ𝑤_vec, $δ𝑤, $𝑖, $𝑗, $K)

### 3 features in common hists - vector of vec hists - gradients/weight in single vector - with node assignations
hist_δ𝑤_vec = [[zeros(3 * nbins) for n in 1:16] for j in 1:nvars]
δ𝑤 = rand(3 * n)
𝑛 = sample(1:16, n)

function iter_7(X_bin, hist_δ𝑤_vec::Vector{Vector{Vector{T}}}, δ𝑤::Vector{T}, 𝑖, 𝑗, 𝑛) where T
    # [hist_δ𝑤_vec[j][n] .= 0.0 for n in 𝑛]
    @inbounds @threads for j in 𝑗
        @inbounds @simd for i in 𝑖
            id = 3 * i - 2
            hid = 3 * X_bin[i,j] - 2
            n = 𝑛[i]

            hist_δ𝑤_vec[j][n][hid] += δ𝑤[id]
            hist_δ𝑤_vec[j][n][hid + 1] += δ𝑤[id + 1]
            hist_δ𝑤_vec[j][n][hid + 2] += δ𝑤[id + 2]
        end
    end
end

@time iter_7(X_bin, hist_δ𝑤_vec, δ𝑤, 𝑖, 𝑗, 𝑛)
@btime iter_7($X_bin, $hist_δ𝑤_vec, $δ𝑤, $𝑖, $𝑗, $𝑛)

using StatsBase:sample
using BenchmarkTools

n_obs = Int(1e6)
n_vars = 100
n_bins = 64
K = 3
𝑖 = collect(1:n_obs)
δ = rand(K, n_obs)
hist = zeros(K, n_bins, n_vars);
X_bin = sample(UInt8.(1:n_bins), n_obs * n_vars);
X_bin = reshape(X_bin, n_obs, n_vars);
𝑖_sample = sample(𝑖, Int(n_obs / 2), ordered=true)

function iter_1(X_bin, hist, δ, 𝑖)
    hist .= 0.0
    @inbounds @simd for i in 𝑖
        @inbounds for k in 1:3
            hist[k, X_bin[i,1], 1] += δ[k, i]
        end
    end
end

@btime iter_1($X_bin, $hist, $δ, $𝑖_sample)


text:
export accuracy, bce, logistic, nll, zeroone

"""
    nll(scores, labels; dims=1, average=true)

Return the negative log likelihood for a single batch of data given an unnormalized `scores`
matrix and an `Integer` array of correct `labels`. The `scores` matrix should have size
`(classes,instances)` if `dims=1` or `(instances,classes)` if `dims=2`. `labels[i]` should be
in `1:classes` to indicate the correct class for instance i, or 0 to skip instance i.

The return value is `(total/count)` if `average=true` and `(total,count)` if `average=false`
where `count` is the number of instances not skipped (i.e. `label != 0`) and `total` is their
total negative log likelihood.

## Example

Let's assume that there are three classes (cat, dog, ostrich) and just 2 instances with
the unnormalized score `scores[:,1]` and `scores[:,2]` respectively. The first instance
is actually a cat and the second instance a dog:

    scores = [12.2    0.3;
               2.0   21.5;
               0.0  -21.0]
    labels = [1, 2]
    nll(scores,labels)
    # returns 2.1657e-5

The probabilites are derived from the scores and the negative log-probabilities corresponding
to the labels are averaged:

    probabilites = exp.(scores) ./ sum(exp.(scores),dims=1)
    -(log(probabilites[labels[1],1]) + log(probabilites[labels[2],2]))/2
    # returns 2.1657e-5

"""
function nll(scores,labels::AbstractArray{<:Integer}; dims=1, average=true)
    indices = findindices(scores,labels,dims=dims)
    lp = logsoftmax(scores,dims=dims)[indices]
    average ? (-sum(lp) / length(lp)) : (-sum(lp), length(lp))
end

"""
    accuracy(scores, labels; dims=1, average=true)

Given an unnormalized `scores` matrix and an `Integer` array of correct `labels`, return the
ratio of instances where the correct label has the maximum score. `dims=1` means instances are
in columns, `dims=2` means instances are in rows. Use `average=false` to return the pair
(ncorrect,count) instead of the ratio (ncorrect/count). The valid labels should be integers in
the range `1:numclasses`, if `labels[i] == 0`, instance i is skipped.

"""
function accuracy(scores,labels::AbstractArray{<:Integer}; dims=1, average=true)
    indices = findindices(scores,labels,dims=dims)
    ycpu = convert(Array,value(scores))
    (maxval,maxind) = findmax(ycpu,dims=dims)
    maxind = LinearIndices(ycpu)[maxind]
    maxind = vec(maxind)[vec(labels) .!= 0]
    correct = (maxind .== indices)
    average ? (sum(correct) / length(correct)) : (sum(correct), length(correct))
end

"""
    logistic(scores, labels; average=true)

Computes logistic loss given predicted unnormalized scores and answer labels for a binary
prediction task.

    log.(1 .+ exp.(-labels .* scores))

Label values should be {-1,1}. Scores are unrestricted.  The return value is `(total/count)`
if `average=true` and `(total,count)` if `average=false` where `count` is the number of
instances and `total` is their total loss.

See also `bce` which computes the same loss with {0,1} labels.

Reference: https://towardsdatascience.com/nothing-but-numpy-understanding-creating-binary-classification-neural-networks-with-e746423c8d5c
"""
function logistic(scores, labels::AbstractVector{<:Integer}; average=true)
    labels = ((labels .+ 1) .÷ 2) # (-1,1)->(0,1)
    bce(scores,labels; average=average)
end

"""
    bce(scores, labels; average=true)

Computes binary cross entropy loss given predicted unnormalized scores and answer labels for a
binary prediction task. Label values should be in {0,1}. Scores are unrestricted and will be
converted to probabilities using

    probs = 1 ./ (1 .+ exp.(-scores))

The loss calculated is

    -(labels .* log.(probs) .+ (1 .- labels) .* log.(1 .- probs))

The return value is `(total/count)` if `average=true` and `(total,count)` if `average=false`
where `count` is the number of instances and `total` is their total loss.

See also `logistic` which computes the same loss with {-1,1} labels.

Reference: https://towardsdatascience.com/nothing-but-numpy-understanding-creating-binary-classification-neural-networks-with-e746423c8d5c
"""
function bce(scores, labels::AbstractVector{<:Integer}; average=true)
    labels = oftype(scores,labels)
    l = max.(0, scores) .- labels .* scores .+ log.(1 .+ exp.(-abs.(scores)))
    average ? (sum(l)/length(l)) : (sum(l),length(l))
end

# Indexing help for nll and accuracy

function findindices(scores, labels::AbstractArray{<:Integer}; dims=1)
    ninstances = length(labels)
    nindices = 0
    indices = Vector{Int}(undef,ninstances)
    if dims == 1                   # instances in first dimension
        y1 = size(scores,1)
        y2 = div(length(scores),y1)
        if ninstances != y2; throw(DimensionMismatch()); end
        @inbounds for j=1:ninstances
            if labels[j] == 0; continue; end
            indices[nindices+=1] = (j-1)*y1 + labels[j]
        end
    elseif dims == 2               # instances in last dimension
        y2 = size(scores,ndims(scores))
        y1 = div(length(scores),y2)
        if ninstances != y1; throw(DimensionMismatch()); end
        @inbounds for j=1:ninstances
            if labels[j] == 0; continue; end
            indices[nindices+=1] = (labels[j]-1)*y1 + j
        end
    else
        error("findindices only supports dims = 1 or 2")
    end
    return (nindices == ninstances ? indices : view(indices,1:nindices))
end

"""
    nll(model; data, dims=1, average=true, o...)

Compute the negative log likelihood for a model over a dataset:

    nll(model(inputs; kwargs...), labels; dims) for (inputs,labels) in data

and return `(total/count)` if `average=true` or `(total,count)` if `average=false` where
`count` is the number of instances not skipped (instances with `label==0` are skipped) and
`total` is their total negative log likelihood.

The `model` should be a function returning scores given inputs, and data should be an iterable
of `(inputs,labels)` pairs. The valid labels should be integers in the range `1:numclasses`,
if `labels[i] == 0`, instance i is skipped.

"""
function nll(model; data, dims=1, average=true, o...)
    sum = cnt = 0
    for (x,y) in data
        (z,n) = nll(model(x; o...), y; dims=dims, average=false)
        sum += z; cnt += n
    end
    average ? sum / cnt : (sum, cnt)
end

"""
    accuracy(model; data, dims=1, average=true, o...)

Compute the number of correct predictions of a model over a dataset:

    accuracy(model(inputs; kwargs...), labels; dims) for (inputs,labels) in data

and return `(ncorrect/count)` if `average=true` or `(ncorrect,count)` if `average=false` where
`count` is the number instances not skipped (instances with `label==0` are skipped) and
`ncorrect` is the number of them correctly labeled by the model.

The `model` should be a function returning scores given inputs, and data should be an iterable
of `(inputs,labels)` pairs. The valid labels should be integers in the range `1:numclasses`,
if `labels[i] == 0`, instance i is skipped.

"""
function accuracy(model; data, dims=1, average=true, o...)
    sum = cnt = 0
    for (x,y) in data
        (z,n) = accuracy(model(x; o...), y; dims=dims, average=false)
        sum += z
        cnt += n
    end
    average ? sum / cnt : (sum, cnt)
end

# DEPRECATE:

# The two-arg model,data calls cause mix-ups with the scores,labels calls:
function nll(model,data; o...)
    @warn "nll(model,data; o...) is deprecated, please use nll(model; data=data, o...)" maxlog=1
    nll(model; data=data, o...)
end

function accuracy(model,data; o...)
    @warn "accuracy(model,data; o...) is deprecated, please use accuracy(model; data=data, o...)" maxlog=1
    accuracy(model; data=data, o...)
end

# We need the (model,x,y) interface to implement regularization:
function nll(f, x, y; dims=1, average=true, o...)
    @warn "nll(f,x,y; o...) is deprecated, please use nll(f(x),y; o...) instead." maxlog=1
    nll(f(x; o...), y; dims=dims, average=average)
end

function accuracy(f, x, y; dims=1, average=true, o...)
    @warn "accuracy(f,x,y; o...) is deprecated, please use accuracy(f(x),y; o...) instead." maxlog=1
    accuracy(f(x; o...), y; dims=dims, average=average)
end

# We need the (weights,data,predict) interface to support the old interface:
function nll(w, data, f::Function; dims=1, average=true, o...)
    @warn "nll(weights,data,func; o...) is deprecated, please use nll(x->func(weights,x); data=data, o...) instead." maxlog=1
    nll(x->f(w,x;o...), data; dims=dims, average=average)
end

function accuracy(w, data, f::Function; dims=1, average=true, o...)
    @warn "accuracy(weights,data,func; o...) is deprecated, please use accuracy(x->func(weights,x); data=data, o...) instead." maxlog=1
    accuracy(x->f(w,x;o...), data; dims=dims, average=average)
end

"zeroone loss is equal to 1 - accuracy"
function zeroone(x...; o...)
    @warn "zeroone() is deprecated, please use 1-accuracy()" maxlog=1
    1 - accuracy(x...; o...)
end


text:
typeof_model_elements(model::StateSpaceModel) = eltype(model.system.y)
num_states(model::StateSpaceModel) = num_states(system(model))
system(model::StateSpaceModel) = model.system
isunivariate(model::StateSpaceModel) = isa(model.system.y, Vector)


text:
using Test
using SafeTestsets

@safetestset "Legendre tests" begin include("legendre.jl") end
@safetestset "Bessel tests" begin include("bessel.jl") end
@safetestset "Mie series tests" begin include("scattering.jl") end

max_stars_repo_name: Riarrieta/MieSeries

max_issues_repo_name: Riarrieta/MieSeries

max_forks_repo_name: Riarrieta/MieSeries


text:
__precompile__()
module Multilane

import StatsBase: WeightVec, sample

using POMDPs
import POMDPs: actions, discount, isterminal, iterator
import POMDPs: rand, reward
import POMDPs: solve, action
import POMDPs: update, initialize_belief
#using GenerativeModels
#import GenerativeModels: generate_s, generate_sr, initial_state, generate_o, generate_sor

import Distributions: Dirichlet, Exponential, Gamma, rand

import Iterators

import Base: ==, hash, length, vec, +, -, *, .*, ^, .^, .-, /, sqrt, zero, abs, max

# import POMDPToolbox: Particle, ParticleBelief

using DataFrames
using DataArrays
using ProgressMeter
using PmapProgressMeter
using POMDPToolbox
using ParticleFilters

import MCTS # so that we can define node_tag, etc.
# using RobustMCTS # for RobustMDP
import BasicPOMCP # for particle filter

# using Reel
# using AutomotiveDrivingModels
# using AutoViz

import Mustache
import JLD

using Parameters

include("debug.jl")

# package code goes here
export
    PhysicalParam,
    CarState,
    MLState,
    MLAction,
    MLObs,
    BehaviorModel,
    MLMDP,
    MLPOMDP,
    OriginalMDP,
    OriginalRewardModel,
    IDMMOBILModel,
    IDMParam,
    MOBILParam,
    IDMMOBILBehavior,
    MLPhysicalState,
    CarPhysicalState,
    RobustNoCrashMDP,
    FixedBehaviorNoCrashMDP,
    StochasticBehaviorNoCrashMDP,
    RobustMLSolver,
    RobustMLPolicy,
    MLPOMDPSolver,
    MLPOMDPAgent,
    DiscreteBehaviorSet,
    DiscreteBehaviorBelief,
    WeightUpdateParams,
    UniformIDMMOBIL,
    standard_uniform,
    MLMPCSolver,
    MLPOMDPSolver

export
    NoCrashRewardModel,
    NoCrashIDMMOBILModel,
    NoCrashMDP,
    NoCrashPOMDP,
    Simple,
    SimpleSolver,
    BehaviorSolver,
    IDMLaneSeekingSolver

export
    SingleBehaviorSolver,
    SingleBehaviorPolicy,
    single_behavior_state

export
    get_neighborhood, #testing VVV
    get_dv_ds,
    is_lanechange_dangerous,
    get_rear_accel, #testing ^^^
    get_idm_dv,
    get_mobil_lane_change,
    is_crash,
    detect_braking,
    max_braking,
    braking_ids

export #data structure stuff
    ste,
    test_run,
    evaluate,
    merge_results!,
    rerun,
    assign_keys,
    test_run_return_policy,
    lambda,
    TestSet,
    gen_initials,
    run_simulations,
    fill_stats!,
    sbatch_spawn,
    gather_results,
    relaxed_initial_state,
    nan

export # POMDP belief stuff
    ParticleUpdater,
    create_belief,
    update,
    rand,
    sample,
    initialize_belief,
    BehaviorRootUpdater,
    BehaviorRootUpdaterStub,
    AggressivenessBelief,
    AggressivenessUpdater,
    agg_means,
    agg_stds,
    aggressiveness,
    BehaviorParticleUpdater,
    BehaviorParticleBelief,
    ParticleGenerator,
    param_means,
    param_stds

export
    MaxBrakeMetric,
    NumBehaviorBrakesMetric

export
    GaussianCopula

export
    include_visualization,
    #generate_sr,
    create_state,
    NoCrashProblem,

    HardBrakeOverlay,
    InfoOverlay,
    CarIDOverlay,
    CarVelOverlay,
    gen_initials_cz,
    relaxed_initial_state_cz,
    relaxed_initial_state_cz2,
    test_run_cz,
    simulate_cz,
    MLMPCAgent,
    BehaviorGenerator,
    standard_uniform,
    relaxed_initial_state_TRI,
    relaxed_initial_state_cz_control,
    get_traffic_density,
    relaxed_initial_state_cz_emergency,
    EmergencyRewardModel,
    judge_collision

include("triangular.jl")
include("physical.jl")
include("MDP_types.jl")
include("crash.jl")
include("IDM.jl")
include("MOBIL.jl")
include("behavior.jl")
include("copula.jl")
include("behavior_gen.jl")
include("success_model.jl")
include("no_crash_model.jl")

include("metrics.jl")
include("evaluation.jl")
include("test_sets.jl")
include("relax.jl")
include("robust_mdp.jl")
include("pomdp_glue.jl")
include("most_likely_mpc.jl")
include("single_behavior.jl")

#include("heuristics.jl")

include("beliefs.jl")
include("aggressiveness_particle_filter.jl")
include("Traffic_Simulation_TRI.jl")
#=
include("uniform_particle_filter.jl")
include("heuristics.jl")
include("tree_vis.jl")
include("sherlock.jl")
=#
include_visualization() = include(joinpath(Pkg.dir("Multilane"),"src","visualization.jl"))

#if gethostname() == "Theresa"
kkkk = 1
if kkkk == 1
    println("Automatically loading visualization components.")
    include("visualization.jl")
    export
        display_sim,
        show_state,
        show_sim,
        display_sim,
        save_frame,
        visualize,
        visualize_cz
end

end # module

max_stars_repo_name: zhcao/HighwayExit

max_issues_repo_name: zhcao/HighwayExit

max_forks_repo_name: zhcao/HighwayExit


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Author: Rafael de P. Soares and Paula B. Staudt
#* $Id$
#*--------------------------------------------------------------------
type pfr
	pfr()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of components"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of reactions"
			)),
			fill(DanaReal(Dict{Symbol,Any}(
				:Brief=>"Stoichiometric Matrix"
			)),(NComp, NReac)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of points of discretization",
				:Default=>10
			)),
			fill(molweight(Dict{Symbol,Any}(
				:Brief=>"Component Mol Weight"
			)),(NComp)),
			length(Dict{Symbol,Any}(
				:Brief=>"Reactor Length"
			)),
			area(Dict{Symbol,Any}(
				:Brief=>"Cross section area"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet Stream",
				:PosX=>0,
				:PosY=>0.5076,
				:Symbol=>"_{in}"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Outlet Stream",
				:PosX=>1,
				:PosY=>0.5236,
				:Symbol=>"_{out}"
			)),
			fill(vapour_stream()),
			fill(vol_mol()),
			fill(dens_mass()),
			fill(heat_rate()),
			fill(mol(Dict{Symbol,Any}(
				:Brief=>"Molar holdup"
			)),(NComp, NDisc)),
			fill(mol(Dict{Symbol,Any}(
				:Brief=>"Molar holdup"
			)),(NDisc)),
			fill(conc_mol(Dict{Symbol,Any}(
				:Brief=>"Components concentration",
				:Lower=>-1e-6
			)),(NComp, NDisc)),
			fill(energy(Dict{Symbol,Any}(
				:Brief=>"Total Energy Holdup on element"
			)),(NDisc)),
			fill(reaction_mol()),
			fill(heat_reaction()),
			[
				:(str(1).F = Inlet.F),
				:(str(1).T = Inlet.T),
				:(str(1).P = Inlet.P),
				:(str(1).z = Inlet.z),
				:(Outlet.F = str(NDisc+1).F),
				:(Outlet.T = str(NDisc+1).T),
				:(Outlet.P = str(NDisc+1).P),
				:(Outlet.z = str(NDisc+1).z),
				:(Outlet.h = str(NDisc+1).h),
				:(Outlet.v = str(NDisc+1).v),
				:(diff(M([1:NComp],[1:NDisc])) = (str([1:NDisc]).F*str([1:NDisc]).[1:NDisc]([1:NComp]) - str([1:NDisc]+1).F*str([1:NDisc]+1).[1:NDisc]([1:NComp])) + sum(stoic([1:NComp],:)*r(:, [1:NDisc])) * Across*L/NDisc),
				:(diff(E([1:NDisc])) = str([1:NDisc]).F*str([1:NDisc]).h - str([1:NDisc]+1).F*str([1:NDisc]+1).h + sum(Hr(:,[1:NDisc])*r(:,[1:NDisc])) * Across*L/NDisc - q([1:NDisc])),
				:(E([1:NDisc]) = Mt([1:NDisc])*str([1:NDisc]+1).h - str([1:NDisc]+1).P*Across*L/NDisc),
				:(str([1:NDisc]+1).F*vol([1:NDisc]+1)*rho([1:NDisc]+1) = str([1:NDisc]).F*vol([1:NDisc])*rho([1:NDisc])),
				:(C(:,[1:NDisc]) * Across*L/NDisc = M(:,[1:NDisc])),
				:(Across*L/NDisc = Mt([1:NDisc]) * vol([1:NDisc])),
				:(str([1:NDisc]+1).[1:NDisc] * Mt([1:NDisc]) = M(:,[1:NDisc])),
				:(vol([1:NDisc+1]) = PP.VapourVolume(str([1:NDisc+1]).T, str([1:NDisc+1]).P, str([1:NDisc+1]).[1:NDisc+1])),
				:(rho([1:NDisc+1]) = PP.VapourDensity(str([1:NDisc+1]).T, str([1:NDisc+1]).P, str([1:NDisc+1]).[1:NDisc+1])),
			],
			[
				"Inlet boundary","","","","Outlet boundary","","","","","","Component Molar Balance","Energy Balance","Energy Holdup","mass flow is considered constant","Molar concentration","Geometrical constraint","Molar fraction","Specific Volume","Specific Mass",
			],
			[:PP,:NComp,:NReac,:stoic,:NDisc,:Mw,:L,:Across,],
			[:Inlet,:Outlet,:str,:vol,:rho,:q,:M,:Mt,:C,:E,:r,:Hr,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	NReac::DanaInteger
	stoic::Array{DanaReal}
	NDisc::DanaInteger
	Mw::Array{molweight}
	L::length
	Across::area
	Inlet::stream
	Outlet::stream
	str::Array{vapour_stream}
	vol::Array{vol_mol}
	rho::Array{dens_mass}
	q::Array{heat_rate}
	M::Array{mol}
	Mt::Array{mol}
	C::Array{conc_mol}
	E::Array{energy}
	r::Array{reaction_mol}
	Hr::Array{heat_reaction}
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export pfr
function set(in::pfr)
	Mw = PP.MolecularWeight()

end
function setEquationFlow(in::pfr)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
	addEquation(11)
	addEquation(12)
	addEquation(13)
	addEquation(14)
	addEquation(15)
	addEquation(16)
	addEquation(17)
	#Hf(:,z) = PP.IdealGasEnthalpyOfFormation(str(z+1).T);
	#Hr(:,z) = -sum(stoic*Hf(:, z));
	addEquation(18)
	addEquation(19)
end
function atributes(in::pfr,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Brief]="Model of a Generic PFR with constant mass holdup"
	fields[:Icon]="icon/pfr"
	fields[:Info]="== Requires the information of ==
* Reaction values
* Heat of reaction
* Pressure profile
"
	drive!(fields,_)
	return fields
end
pfr(_::Dict{Symbol,Any})=begin
	newModel=pfr()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
include("map.jl")

@valueequals mutable struct Side
    map::Map
    data::Dict{String, Any}

    Side(map::Map, data::Dict{String, Any}=Dict{String, Any}()) = new(map, data)
end

const decoderBlacklist = String[
    "Filler", "Style", "levels"
]

function Base.Dict(s::Side)
    data = Dict(s.map)
    encodeMetadata!(data, s.data)

    return data
end

function loadDataAsDict!(target::Dict{String, Any}, data::DecodedElement)
    name = data.name
    target[name] = Dict{String, Any}()

    for (attr, value) in attributes(data)
        target[name][attr] = value
    end

    for child in children(data)
        loadDataAsDict!(target[name], child)
    end
end

function dictAsEncodedData(name::String, data::Dict{String, Any})
    res = Dict{String, Any}(
        "__name" => name,
        "__children" => Dict{String, Any}[]
    )

    for (k, v) in data
        if isa(v, Dict)
            push!(res["__children"], dictAsEncodedData(k, v))

        else
            res[k] = v
        end
    end

    if isempty(res["__children"])
        delete!(res, "__children")
    end

    return res
end

function loadMetadata(data::DecodedElement)
    res = Dict{String, Any}()

    for child in children(data)
        if !(child.name in decoderBlacklist)
            loadDataAsDict!(res, child)
        end
    end

    return res
end

loadSide(fn::String) = loadSide(decodeMap(fn))
function loadSide(data::DecodedElement)
    return Side(loadMap(data), loadMetadata(data))
end

function encodeMetadata!(target::Dict{String, Any}, data::Dict{String, Any})
    for (k, v) in data
        push!(target["__children"], dictAsEncodedData(k, v))
    end
end

function encodeSide(side::Side, outfile::String)
    return encodeMap(Dict(side), outfile)
end

function getSideName(side::Side)
    if haskey(side.data, "meta") && haskey(side.data["meta"], "Name") && !isempty(side.data["meta"]["Name"])
        return side.data["meta"]["Name"]

    else
        return side.map.package
    end
end

max_stars_repo_name: CelestialCartographers/Maple

max_issues_repo_name: CelestialCartographers/Maple

max_forks_repo_name: CelestialCartographers/Maple


text:
# An implementation of CGNE for the solution of the consistent
# (under/over-determined or square) linear system
#
#  Ax = b.
#
# The method seeks to solve the minimum-norm problem
#
#  min ‖x‖²  s.t. Ax = b,
#
# and is equivalent to applying the conjugate gradient method
# to the linear system
#
#  AA'y = b.
#
# This method is also known as Craig's method, CGME, and other
# names, and is described in
#
# J. E. Craig. The N-step iteration procedures.
# Journal of Mathematics and Physics, 34(1):64--73, 1955.
#
# which is based on Craig's thesis from MIT:
#
# J. E. Craig. Iterations Procedures for Simultaneous Equations.
# Ph.D. Thesis, Department of Electrical Engineering, MIT, 1954.
#
# Dominique Orban, <dominique.orban@gerad.ca>
# Montréal, QC, April 2015.

export cgne

"""Solve the consistent linear system

  Ax + √λs = b

using the Conjugate Gradient (CG) method, where λ ≥ 0 is a regularization
parameter. This method is equivalent to applying CG to the normal equations
of the second kind

  (AA' + λI) y = b

but is more stable. When λ = 0, this method solves the minimum-norm problem

  min ‖x‖₂  s.t. Ax = b.

When λ > 0, it solves the problem

  min ‖(x,s)‖₂  s.t. Ax + √λs = b.

CGNE produces monotonic errors ‖x-x*‖₂ but not residuals ‖r‖₂.
It is formally equivalent to CRAIG, though can be slightly less accurate,
but simpler to implement. Only the x-part of the solution is returned.

A preconditioner M may be provided in the form of a linear operator.
"""
function cgne(A :: AbstractLinearOperator, b :: AbstractVector{T};
              M :: AbstractLinearOperator=opEye(size(A,1)), λ :: Float64=0.0,
              atol :: Float64=1.0e-8, rtol :: Float64=1.0e-6, itmax :: Int=0,
              verbose :: Bool=false) where T <: Number

  m, n = size(A);
  size(b, 1) == m || error("Inconsistent problem size");
  verbose && @printf("CGNE: system of %d equations in %d variables\n", m, n);

  x = zeros(T, n);
  r = copy(b)
  z = M * r
  rNorm = @knrm2(m, r)   # Marginally faster than norm(r)
  rNorm ≈ 0 && return x, SimpleStats(true, false, [rNorm], [], "x = 0 is a zero-residual solution");
  λ > 0 && (s = copy(r));

  # The following vector copy takes care of the case where A is a LinearOperator
  # with preallocation, so as to avoid overwriting vectors used later. In other
  # case, this should only add minimum overhead.
  p = copy(A' * z);

  # Use ‖p‖ to detect inconsistent system.
  # An inconsistent system will necessarily have AA' singular.
  # Because CGNE is equivalent to CG applied to AA'y = b, there will be a
  # conjugate direction u such that u'AA'u = 0, i.e., A'u = 0. In this
  # implementation, p is a substitute for A'u.
  pNorm = @knrm2(n, p)

  γ = @kdot(m, r, z)  # Faster than γ = dot(r, z);
  iter = 0;
  itmax == 0 && (itmax = m + n);

  rNorms = [rNorm;];
  ɛ_c = atol + rtol * rNorm;  # Stopping tolerance for consistent systems.
  ɛ_i = atol + rtol * pNorm;  # Stopping tolerance for inconsistent systems.
  verbose && @printf("%5s  %8s\n", "Aprod", "‖r‖")
  verbose && @printf("%5d  %8.2e\n", 1, rNorm);

  status = "unknown";
  solved = rNorm <= ɛ_c;
  inconsistent = (rNorm > 1.0e+2 * ɛ_c) && (pNorm <= ɛ_i);
  tired = iter >= itmax;

  while ! (solved || inconsistent || tired)
    q = A * p;
    λ > 0 && @kaxpy!(m, λ, s, q)
    δ = @kdot(n, p, p)   # Faster than dot(p, p);
    λ > 0 && (δ += λ * @kdot(m, s, s))
    α = γ / δ;
    @kaxpy!(n,  α, p, x)     # Faster than x = x + α * p;
    @kaxpy!(m, -α, q, r)     # Faster than r = r - α * q;
    z = M * r
    γ_next = @kdot(m, r, z)  # Faster than γ_next = dot(r, z);
    β = γ_next / γ;
    @kscal!(n, β, p)
    @kaxpy!(n, 1.0, A' * z, p)   # Faster than p = A' * z + β * p;
    pNorm = @knrm2(n, p)
    if λ > 0
      @kscal!(m, β, s)
      @kaxpy!(m, 1.0, r, s)   # s = r + β * s;
    end
    γ = γ_next;
    rNorm = sqrt(γ_next);
    push!(rNorms, rNorm);
    iter = iter + 1;
    verbose && @printf("%5d  %8.2e\n", 1 + 2 * iter, rNorm);
    solved = rNorm <= ɛ_c;
    inconsistent = (rNorm > 1.0e+2 * ɛ_c) && (pNorm <= ɛ_i);
    tired = iter >= itmax;
  end

  status = tired ? "maximum number of iterations exceeded" : (inconsistent ? "system probably inconsistent" : "solution good enough given atol and rtol")
  stats = SimpleStats(solved, inconsistent, rNorms, [], status);
  return (x, stats);
end


text:
"""
    struct YenState{T, U}

Designed for yen k-shortest-paths calculations.
"""
struct YenState{T,U <: Integer} <: AbstractPathState
    dists::Vector{T}
    paths::Vector{Vector{U}}
end

"""
    yen_k_shortest_paths(g, source, target, distmx=weights(g), K=1; maxdist=Inf);

Perform [Yen's algorithm](http://en.wikipedia.org/wiki/Yen%27s_algorithm)
on a graph, computing k-shortest distances between `source` and `target` other vertices.
Return a [`YenState`](@ref) that contains distances and paths.
"""
function yen_k_shortest_paths(g::AbstractGraph,
    source::U,
    target::U,
    distmx::AbstractMatrix{T}=weights(g),
    K::Int=1;
    maxdist=Inf) where T <: Real where U <: Integer

    source == target && return YenState{T,U}([U(0)], [[source]])

    dj = dijkstra_shortest_paths(g, source, distmx)
    path = enumerate_paths(dj)[target]
    isempty(path) && return YenState{T,U}(Vector{T}(), Vector{Vector{U}}())

    dists = Array{T,1}()
    push!(dists, dj.dists[target])
    A = [path]
    B = PriorityQueue()
    gcopy = deepcopy(g)

    for k = 1:(K - 1)
        for j = 1:length(A[k])
            # Spur node is retrieved from the previous k-shortest path, k − 1
            spurnode = A[k][j]
            #  The sequence of nodes from the source to the spur node of the previous k-shortest path
            rootpath = A[k][1:j]

            # Store the removed edges
            edgesremoved = Array{Tuple{Int,Int},1}()
            # Remove the links of the previous shortest paths which share the same root path
            for ppath in A
                if length(ppath) > j && rootpath == ppath[1:j]
                    u = ppath[j]
                    v = ppath[j + 1]
                    if has_edge(gcopy, u, v)
                        rem_edge!(gcopy, u, v)
                        push!(edgesremoved, (u, v))
                    end
                end
            end

            # Remove node of root path and calculate dist of it
            distrootpath = 0.
            for n = 1:(length(rootpath) - 1)
                u = rootpath[n]
                nei = copy(neighbors(gcopy, u))
                for v in nei
                    rem_edge!(gcopy, u, v)
                    push!(edgesremoved, (u, v))
                end

                # Evaluate distante of root path
                v = rootpath[n + 1]
                distrootpath += distmx[u, v]
            end

            # Calculate the spur path from the spur node to the sink
            djspur = dijkstra_shortest_paths(gcopy, spurnode, distmx)
            spurpath = enumerate_paths(djspur)[target]
            if !isempty(spurpath)
                # Entire path is made up of the root path and spur path
                pathtotal = [rootpath[1:(end - 1)]; spurpath]
                distpath  = distrootpath + djspur.dists[target]
                # Add the potential k-shortest path to the heap
                if !haskey(B, pathtotal)
                    enqueue!(B, pathtotal, distpath)
                end
            end

            for (u, v) in edgesremoved
                add_edge!(gcopy, u, v)
            end
        end

        # No more paths in B
        isempty(B) && break
        mindistB = peek(B)[2]
        # The path with minimum distance in B is higher than maxdist
        mindistB > maxdist && break
        push!(dists, peek(B)[2])
        push!(A, dequeue!(B))
    end

    return YenState{T,U}(dists, A)
end


text:
using Flux
using Zygote: Buffer

struct HoltWinters
    input_size
    output_size
    freq

    level_smooth::Vector{Float32}
    trend_smooth::Vector{Float32}
    season_smooth::Vector{Float32}

    level_init::Vector{Float32}
    trend_init::Vector{Float32}
    season_init::Vector{Float32}
end

function HoltWinters(input_size, output_size, freq)
    level_smooth = [0.5]
    trend_smooth = [0.5]
    season_smooth = [0.5]

    level_init = [0.5]
    trend_init = [0.5]
    season_init = repeat([0.5], freq)

    HoltWinters(
        input_size, output_size, freq,
        level_smooth, trend_smooth, season_smooth,
        level_init, trend_init, season_init)
end

function (model::HoltWinters)(input)
    level_smooth = model.level_smooth[1]
    trend_smooth = model.level_smooth[1]
    season_smooth = model.season_smooth[1]

    level = Buffer([], Float64, model.input_size)
    trend = Buffer([], Float64, model.input_size)
    season = Buffer([], Float64, model.input_size)

    for t in 1:model.input_size
        level_prev = t == 1 ? model.level_init[1] : level[t - 1]
        trend_prev = t == 1 ? model.trend_init[1] : trend[t - 1]
        season_prev = t <= model.freq ? model.season_init[t] : season[t - model.freq]

        level[t] = (level_smooth * (input[t] - season_prev)
            + (1 - level_smooth) * (level_prev + trend_prev))

        trend[t] = (trend_smooth * (level[t] - level_prev)
            + (1 - trend_smooth) * trend_prev)

        season[t] = (season_smooth * (input[t] - level_prev - trend_prev)
            + (1 - season_smooth) * season_prev)
    end

    pred = Buffer([], Float64, model.output_size)

    for h in 1:model.output_size
        k = floor(Int64, (h - 1) / model.freq)

        pred[h] = (level[model.input_size]
            + h * trend[model.input_size]
            + season[model.input_size + h - model.freq * (k + 1)])
    end

    copy(pred)
end

#model = HoltWinters(input_size, output_size, info.freq)
#model(windows[2][1])


text:
module TestPITMatrix

using Copulas
using Base.Test

## single observation
kk = PITMatrix([0.1, 0.3])

kk2 = [0.1 NaN;
       NaN 0.3]
expOut = PITMatrix(reshape(kk2, 2, 2, 1))

@test isequal(kk, expOut)

## multiple observations
kk = PITMatrix([0.1 0.3;
                0.4 0.7;
                0.6 0.4])

@test size(dg(kk)) == (2, 2, 3)
@test kk[1, 1] == [0.1, 0.4, 0.6]
@test isnan(dg(kk)[1, 2])

## throw error during construction
@test_throws Exception PITMatrix([0.1 0.3;
                                  0.4 1.2;
                                  0.6 0.4])

@test_throws Exception PITMatrix([0.1 0.3;
                                  0.4 -0.2;
                                  0.6 0.4])

@test_throws Exception PITMatrix([0.1 0.3;
                                  0.4 -0.2;
                                  0.6 0.4])

kk = reshape([0.1 0.2;
              0.4 0.3;
              0.5 0.1], 3, 2, 1)

@test_throws Exception PITMatrix(kk)

end


text:
using Test
using JuMP
using MIPVerify: Linear, check_size
@isdefined(TestHelpers) || include("../../TestHelpers.jl")

@testset "linear.jl" begin

    @testset "Linear" begin
        @testset "With Bias" begin
            @testset "Matched Size" begin
                height = 10
                matrix = ones(2, height)
                bias = ones(height)
                p = Linear(matrix, bias)
                @test p.matrix == matrix
                @test p.bias == bias
                @testset "check_size" begin
                    @test check_size(p, (2, 10)) === nothing
                    @test_throws AssertionError check_size(p, (2, 9))
                end
            end
            @testset "Unmatched Size" begin
                matrix_height = 10
                bias_height = 5
                matrix = ones(2, matrix_height)
                bias = ones(bias_height)
                @test_throws AssertionError Linear(matrix, bias)
            end
            @testset "Multiplying by >1-dimensional array" begin
                height = 15
                matrix = ones(2, height)
                bias = ones(height)
                p = Linear(matrix, bias)
                @test_throws ArgumentError p(ones(3, 5))
            end
        end

        @testset "Base.show" begin
            height = 13
            matrix = ones(37, height)
            bias = ones(height)
            p = Linear(matrix, bias)
            io = IOBuffer()
            Base.show(io, p)
            @test String(take!(io)) == "Linear(37 -> 13)"
        end

        @testset "matmul" begin
            matrix = [[1, 2] [3, 4]]
            bias = [5, 6]
            p = Linear(matrix, bias)
            @testset "Real * Real" begin
                @test p([7, 8]) == [28, 59]
            end
            @testset "JuMP.AffExpr * Real" begin
                m = TestHelpers.get_new_model()
                x = @variable(m, start = 7)
                y = @variable(m, start = 8)
                @test getvalue(p([x, y])) == [28, 59]
            end
        end

    end

end


text:
struct IndexablePoint2D <: VoronoiDelaunay.AbstractPoint2D
    _x::Float64
    _y::Float64
    _index::Int64
end

IndexablePoint2D(x::Float64, y::Float64) = IndexablePoint2D(x, y, -1)
getx(p::IndexablePoint2D) = p._x
gety(p::IndexablePoint2D) = p._y
Base.getindex(p::IndexablePoint2D) = p._index
Base.getindex(::VoronoiDelaunay.Point2D) = -1

for op in [:+, :-]
	@eval begin
		Base.$op(p::VoronoiDelaunay.AbstractPoint2D, q::VoronoiDelaunay.AbstractPoint2D) =
        VoronoiDelaunay.Point2D($op(getx(p), getx(q)), $op(gety(p), gety(q)))
	end
end

Base.:*(a::Float64, p::VoronoiDelaunay.AbstractPoint2D) = VoronoiDelaunay.Point2D(a*getx(p), a*gety(p))

function Base.isapprox(p::VoronoiDelaunay.AbstractPoint2D, q::VoronoiDelaunay.AbstractPoint2D)
    isapprox(getx(p), getx(q)) && isapprox(gety(p), gety(q))
end

function Base.abs2(A, B)
    abs2(getx(A) - getx(B)) + abs2(gety(A) - gety(B))
end

function GeometryBasics.Point2(p::VoronoiDelaunay.AbstractPoint2D)
    GeometryBasics.Point(getx(p), gety(p))
end

getx(p::GeometryBasics.Point2) = p[1]
gety(p::GeometryBasics.Point2) = p[2]


text:
# This file was generated, do not modify it. # hide
@show estimates.measurement[1]
@show mean(estimates.per_fold[1])


text:
using SQLite

db = SQLite.DB("quotes.db")   # SQLiteDB() etc., are not exported

SQLite.tables(db)

SQLite.columns(db,"quotes")

# New syntax, using dataframes

df = DataFrame(SQLite.Query(db, "select count(*) as N from quotes"))

df[1]

df = DataFrame(SQLite.Query(db, "select * from quotes limit 10"))

sql =  "select q.quoname, q.author, c.catname from quotes q ";
sql *= "join categories c on q.cid = c.id limit 5";

df = DataFrame(SQLite.Query(db,sql))

# The feather package can be used to store snapshots of queries
using Feather

Feather.write("QuoView01.feather", df)

# Retrieving fro a feather file is "lazy"
# That is the metadata is read but records are fetched when they are referenced
# So feather files can hold very large datasets with little memory overhead.

dfx = Feather.read("QuoView01.feather");
size(dfx)

dfx[1,1]   # Get first quote

sql =  "select q.quoname from quotes q ";
sql *= " where q.author = 'Oscar Wilde'";
SQLite.Query(db,sql) |> DataFrame

sql =  "select q.quoname, q.author, c.catname from quotes q ";
sql *= "join categories c on q.cid = c.id";

df = DataFrame(SQLite.Query(db,sql));
nrows = size(df)[1];
println("Number of rows: $nrows" )

# Assumes Postgres is running and the Chinook dataset loaded into a database 'chinnok'
using LibPQ, DataStreams
conn = LibPQ.Connection("dbname=chinook");

res = execute(conn, "SELECT count(*) FROM \"Album\"");  # Note escaping the quotes
Data.stream!(res, NamedTuple)

qry = fetch!(NamedTuple, execute(conn, "SELECT count(*) as NA FROM \"Album\""));
qry[1][1]

qry = fetch!(NamedTuple, execute(conn, "SELECT * FROM \"MediaType\""))

qry[:Name]

qry.FirstName[1]

sqlx = """select e."FirstName", e."LastName", count(i."InvoiceId") as "Sales"
 from "Employee" as e
 join "Customer" as c on e."EmployeeId" = c."SupportRepId"
 join "Invoice" as i on i."CustomerId" = c."CustomerId"
 group by e."EmployeeId"
""";

qry = fetch!(NamedTuple, execute(conn, sqlx))

using Printf
for i in 1:length(qry)
    @printf("%s %s has %d sales\n",
        qry.FirstName[i],qry.LastName[i],qry.Sales[i])
end


text:
################################################################################

"Reset hidden Markov model object"
function reset!(self::HMM)
  self.tbM = [x .= [0; repeat(-1:-1, length(self.tbM[1]) - 1)] for x in self.tbM]
  return self.tbM
end

################################################################################

max_stars_repo_name: grabherr/HiddenMarkovModelReaders

max_issues_repo_name: grabherr/HiddenMarkovModelReaders

max_forks_repo_name: grabherr/HiddenMarkovModelReaders


text:
"""
Struct for Geometric Brownian Motion

		gbmProcess=GeometricBrownianMotion(σ::num1,μ::num2) where {num1,num2 <: Number}

Where:\n
		σ	=	volatility of the process.
		μ	=	drift of the process.
		x0	=	initial value.
"""
mutable struct GeometricBrownianMotion{num <: Number, num1 <: Number, num2 <: Number, numtype <: Number} <: ItoProcess{numtype}
    σ::num
    μ::num1
    x0::num2
    function GeometricBrownianMotion(σ::num, μ::num1, x0::num2) where {num <: Number, num1 <: Number, num2 <: Number}
        if σ <= 0
            error("Volatility must be positive")
        else
            zero_typed = zero(num) + zero(num1)
            return new{num, num1, num2, typeof(zero_typed)}(σ, μ, x0)
        end
    end
end

export GeometricBrownianMotion;

function simulate!(X, mcProcess::GeometricBrownianMotion, mcBaseData::AbstractMonteCarloConfiguration, T::Number)
    @assert T > 0
    σ_gbm = mcProcess.σ
    mu_gbm = mcProcess.μ
    μ_bm = mu_gbm - σ_gbm^2 / 2
    simulate!(X, BrownianMotion(σ_gbm, μ_bm), mcBaseData, T)
    S0 = mcProcess.x0
    # @. X=S0*exp(X);
    broadcast!(x -> S0 * exp(x), X, X)

    nothing
end