text:
module SimMdlPrices

using Parameters
using LinearAlgebra
using QuantEcon
using SparseArrays
using StaticArrays
using Statistics
using Infiltrator
using StatsFuns
using BenchmarkTools
using InteractiveUtils

function simulate_ms_var_1_cond_shocks!(x_mat,s_mat,x0,s0,μ,Φ,Σ,Π,u,ϵ,T_start, g_start)

     simulate_markov_switch_init_cond_shock!(s_mat, s0, Π, u, T_start, g_start)
     simulate_msvar_cond_regime_path_shock!(x_mat, x0, s_mat, μ, Φ, Σ, ϵ,T_start, g_start)

end

function simulate_markov_switch_init_cond_shock!(s_mat,s0, Π, u,T_start,g_start)

     T,G = size(u)
     S = size(Π)[1]

     #s_mat = Array{Int64}(undef,T+1,G)
     if T_start == 1
          s_mat[1,:] .= s0
     end

     #Propogate
     Π₀ = cumsum(Π, dims = 2)
     cnt_g = 1
     for g ∈ g_start:G+g_start-1
          cnt_t  = 1
          for t ∈ T_start:T_start+T-1
               for s ∈ 1:S
                    if u[cnt_t, cnt_g] <= Π₀[s_mat[t, g], s]
                         s_mat[t+1, g] = s
                         break
                    end
               end
               cnt_t += 1
          end
          cnt_g += 1
     end

     return s_mat

end

function simulate_msvar_cond_regime_path_shock!(x_mat,x0,s_mat,μ,Φ,Σ,ϵ,T_start, g_start)

     #Unload and Preallocate
     N, N2 = size(x0)
     T = size(ϵ)[2]
     G = size(ϵ)[3]
     @assert N == size(ϵ)[1]

     if N2==1
          x0 = repeat(x0, 1, G)
      end

     #Propogate
     cnt_g = 1
     for g ∈ g_start:G+g_start-1
          cnt_t = 1
          if T_start==1
               x_mat[:, 1, g] = @view(x0[:, cnt_g])
          end

          for t ∈ T_start:T_start+T-1
               s = s_mat[t+1, g]
               x_mat[:, t+1, g] = μ[s] .+ Φ[s] *  @view(x_mat[:,t,g])  .+ Σ[s] * @view(ϵ[:, cnt_t, cnt_g])
               cnt_t += 1
          end

          cnt_g += 1

     end

end

function  construct_structural_matrices_macro(ϕ,δ,ρ₁,ρ₂,β₁,β₂,α₁,α₂,μ_g,μ_pi,m_g,m_r₁,m_r₂,σ_g,σ_pi,σ_r₁,σ_r₂)

     b0₁,b0₂ = construct_b0_macro_array(ϕ,δ,ρ₁,ρ₂,β₁,β₂)
     bm1₁,bm1₂ = construct_bm1_macro_array(μ_g,μ_pi,ρ₁,ρ₂);
     b1₁,b1₂ = construct_b1_macro_array(μ_g,μ_pi,ϕ,ρ₁,ρ₂,α₁,α₂)
     m0₁,m0₂ = construct_m0_macro_array(m_g,m_r₁,m_r₂)
     Γ₁,Γ₂ = construct_gamma_macro_array(σ_g,σ_pi,σ_r₁,σ_r₂)

     return b0₁,b0₂,bm1₁,bm1₂,b1₁,b1₂,m0₁,m0₂,Γ₁,Γ₂
end

function construct_b0_macro_array(ϕ,δ,ρ₁,ρ₂,β₁,β₂)
     b0₁ = @SMatrix [1. 0. ϕ; -δ  1. 0.; -(1. - ρ₁)*β₁ 0. 1.]
     b0₂ = @SMatrix [1. 0. ϕ; -δ  1. 0.; -(1. - ρ₂)*β₂ 0. 1.]

     return b0₁,b0₂
end

function construct_bm1_macro_array(μ_g, μ_pi, ρ₁, ρ₂)
     bm1₁ = @SMatrix [(1. - μ_g) 0. 0.; 0. (1. - μ_pi) 0.; 0. 0. ρ₁]
     bm1₂ = @SMatrix [(1. - μ_g) 0. 0.; 0. (1. - μ_pi) 0.; 0. 0. ρ₂]

     return bm1₁,bm1₂
end

function construct_b1_macro_array(μ_g,μ_pi,ϕ,ρ₁,ρ₂,α₁,α₂)

     b1₁ = @SMatrix [μ_g ϕ 0.; 0. μ_pi 0.; 0. (1. - ρ₁)*α₁ 0.]
     b1₂ = @SMatrix [μ_g ϕ 0.; 0. μ_pi 0.; 0. (1. - ρ₂)*α₂ 0.]

     return b1₁,b1₂
end

function construct_m0_macro_array(m_g,m_r₁,m_r₂)

     m0₁ = SMatrix{3,1,Float64}([m_g 0. m_r₁])
     m0₂ = SMatrix{3,1,Float64}([m_g 0. m_r₂]);

     return m0₁,m0₂
end

function construct_gamma_macro_array(σ_g, σ_pi, σ_r₁, σ_r₂)

     Γ₁ = @SMatrix [σ_g  0.  0.; 0. σ_pi  0.; 0. 0. σ_r₁]
     Γ₂ = @SMatrix [σ_g  0.  0.; 0. σ_pi  0.; 0. 0. σ_r₂]

     return Γ₁,Γ₂
end

function transform_struct_to_rf(b0₁, b0₂, bm1₁, bm1₂, b1₁, b1₂, m0₁, m0₂, Γ₁, Γ₂, π_m, π_d, Π_0, Π_X)

     maxₖ = 1000 #Maximum iterations for REE

     # Compute composite Markov transition matrix and its ergodic distribution
     Π = kron(π_m, π_d)
     mc = MarkovChain(Π)
     q = stationary_distributions(mc)[1]

     #Compute forward-looking RE equilibrium
     F1₁ = b0₁ \ bm1₁
     F1₂ = b0₂ \ bm1₂
     A1₁ = b0₁ \ b1₁
     A1₂ = b0₂ \ b1₂

     Σ1_11 = b0₁ \ Γ₁
     Σ1_12 = b0₁ \ Γ₂
     Σ1_21 = b0₂ \ Γ₁
     Σ1_22 = b0₂ \ Γ₂

     m1₁ = b0₁ \ m0₁
     m1₂ = b0₂ \ m0₂

     #Compute forward solution
     Φ, μ, Σ, det_dum = get_forward_solution_msre(π_m, A1₁, A1₂, F1₁, F1₂, m1₁, m1₂, Σ1_11, Σ1_12, Σ1_21,
                                                                                Σ1_22, maxₖ)

     #Define covariance "cell"
     cov = [Σ[1] * Σ[1]', Σ[2] * Σ[2]', Σ[3] * Σ[3]', Σ[4] * Σ[4]']

     #Compute risk-neautral dynamic
     μ_Q = [ μ[1] - cov[1] * Π_0,
                   μ[2] - cov[2] * Π_0,
                   μ[3] - cov[3] * Π_0,
                   μ[4] - cov[4] * Π_0 ]

     Φ_Q = [ Φ[1] - cov[1] * Π_X,
                    Φ[2] - cov[2] * Π_X,
                    Φ[3] - cov[3] * Π_X,
                    Φ[4] - cov[4] * Π_X ]

     return μ, Φ, μ_Q, Φ_Q, Σ, cov, Π, q, det_dum
end

function get_forward_solution_msre(π_m, A1₁, A1₂, F1₁, F1₂, m1₁, m1₂, Σ1_11, Σ1_12, Σ1_21, Σ1_22, maxₖ)

     Φ, k_term = get_rf_Φ(π_m, F1₁, F1₂, A1₁, A1₂, maxₖ = maxₖ)

     μ, Σ =  compute_μ_Σ_fmmsre(Φ, A1₁, A1₂, m1₁, m1₂, Σ1_11, Σ1_12, Σ1_21, Σ1_22, π_m)

     if k_term == maxₖ
          det_dum = 0
     else
          det_dum = 1
     end

     Φ_out = [Φ[1], Φ[1], Φ[2], Φ[2]]
     μ_out = [μ[1], μ[1], μ[2], μ[2]]
     return Φ_out, μ_out, Σ, det_dum
end

function get_rf_Φ(π_m, F1₁, F1₂, A1₁, A1₂; maxₖ = 1000, tolK = 0.000001)

     B₁ = F1₁
     B₂ = F1₂
     Ωkm1₁ = B₁
     Ωkm1₂ = B₂
     Ωkm1 = [Ωkm1₁,Ωkm1₂]
     global k_term = 0
     global Ωk₁ = similar(Ωkm1₁)
     global Ωk₂ = similar(Ωkm1₁)
     function f(Ωkm1,maxₖ,A1₁, A1₂,π_m,B₁,B₂,tolK)
               for i ∈ 1:maxₖ
                    #update Φₖ
                    Ωk₁ = (I - A1₁ * (π_m[1,1] .* Ωkm1[1] + π_m[1,2] .* Ωkm1[2] ) ) \ B₁
                    Ωk₂ = (I - A1₂ * (π_m[2,1] .* Ωkm1[1] + π_m[2,2] .* Ωkm1[2] ) ) \ B₂

                    #Check for convergence
                    dist = max( maximum(abs.( Ωk₁ - Ωkm1[1] )), maximum(abs.( Ωk₂ - Ωkm1[2] )))

                    #update values or break
                    if (dist <= tolK) || (i==maxₖ)
                         k_term = i
                         break
                    else
                         Ωkm1[1] = Ωk₁
                         Ωkm1[2] = Ωk₂
                    end

               end
          return Ωk₁, Ωk₂, k_term
     end

     Ωk₁,Ωk₂,k_term = f(Ωkm1,maxₖ,A1₁, A1₂,π_m,B₁,B₂,tolK)
     Φ = [ Ωk₁, Ωk₂ ]

     return Φ, k_term

end

function  compute_μ_Σ_fmmsre(Φ, A1₁, A1₂, m1₁, m1₂, Σ1_11, Σ1_12, Σ1_21, Σ1_22, π_m)

     N = size(A1₁)[1]

     m1 = [m1₁; m1₂]
     A₁_star = ( I - A1₁ * ( π_m[1,1]  .*  Φ[1] + π_m[1,2] .* Φ[2] ) )
     A₂_star = ( I - A1₂ * ( π_m[2,1]  .*  Φ[1] + π_m[2,2] .* Φ[2] ) )

     A_star = SMatrix{6,6}([(A₁_star-π_m[1,1].*A1₁) -π_m[1,2].*A1₁;
                     -π_m[2,1].*A1₂ (A₂_star - π_m[2,2].*A1₂) ])

     μ_vec = A_star \ m1
     μ₁ =SMatrix{3,1}(μ_vec[1:N])
     μ₂ = SMatrix{3,1}(μ_vec[N+1:end])

     μ = [μ₁, μ₂]

     Σ_11 = A₁_star \ Σ1_11
     Σ_12 = A₁_star \ Σ1_12
     Σ_21 = A₂_star \ Σ1_21
     Σ_22 = A₂_star \ Σ1_22

     Σ = [Σ_11, Σ_12, Σ_21, Σ_22]

      return μ , Σ
end

@with_kw struct ModelStructuralParams

     Π_0::SMatrix{3,1,Float64,3}
     Π_X::SMatrix{3,3,Float64,9}
     m_g::Float64
     m_r₁::Float64
     m_r₂::Float64
     ϕ::Float64
     δ::Float64
     μ_g::Float64
     μ_pi::Float64
     ρ₁::Float64
     ρ₂::Float64
     β₁::Float64
     β₂::Float64
     α₁::Float64
     α₂::Float64
     σ_g::Float64
     σ_pi::Float64
     σ_r₁::Float64
     σ_r₂::Float64
     σ_y::Float64
     γ_gₐ::Float64
     γ_piₐ::Float64
     γ_cₐ::Float64
     γ_gᵢ::Float64
     γ_piᵢ::Float64
     γ_cᵢ::Float64
     γ_gₒ::Float64
     γ_piₒ::Float64
     γ_cₒ::Float64
     ρₐ::Float64
     ρᵢ::Float64
     ρₒ::Float64
     σ_wₐ::Float64
     σ_zₐ::Float64
     σ_νₐ::Float64
     σ_qₐ::Float64
     σ_wᵢ::Float64
     σ_zᵢ::Float64
     σ_νᵢ::Float64
     σ_qᵢ::Float64
     σ_wₒ::Float64
     σ_zₒ::Float64
     σ_νₒ::Float64
     σ_qₒ::Float64
     λ::Float64
     π_m::SMatrix{2, 2, Float64, 4}
     π_d::SMatrix{2, 2, Float64, 4}

end

@with_kw struct ModelReducedFormParams

     μ::Vector{SMatrix{3,1,Float64,3}}
     μ_Q::Vector{SMatrix{3,1,Float64,3}}
     μ_ν::Vector{SMatrix{6,1,Float64,6}}
     μ_Q_ν::Vector{SMatrix{6,1,Float64,6}}
     Φ::Vector{SMatrix{3,3,Float64,9}}
     Φ_Q::Vector{SMatrix{3,3,Float64,9}}
     Φ_ν::Vector{SMatrix{6,6,Float64,36}}
     Φ_Q_ν::Vector{SMatrix{6,6,Float64,36}}
     Σ::Vector{SMatrix{3,3,Float64,9}}
     Σ_ν::Vector{SMatrix{6,6,Float64,36}}
     cov::Vector{SMatrix{3,3,Float64,9}}
     cov_ν::Vector{SMatrix{6,6,Float64,36}}
     Π::SMatrix{4,4,Float64,16}
     q::Vector{Float64}
     Δ::SMatrix{3, 7, Float64, 21}
     Δ_Q::SMatrix{3, 7, Float64, 21}
     σ_w::SMatrix{3, 1, Float64, 3}
     σ_z::SMatrix{3, 1, Float64, 3}
     σ_q::SMatrix{3, 1, Float64, 3}
     σ_ν::SMatrix{3, 1, Float64, 3}
     det_dum::Int64

end

function construct_structural_parameters(params)

     PiX = @SMatrix [ params[4] 0. 0.; params[7] params[5] 0.; params[8] 0. params[6] ]
     Pi0 = SMatrix{3,1}([params[1] params[2] params[3]])
     beta1 = (1-params[16]) \ params[18]
     beta2 = (1-params[17]) \ params[19]
     alpha1 = (1-params[16])\params[20]
     alpha2 = (1-params[17])\params[21]
     pi_m = @SMatrix [ params[52] 1-params[52]; 1-params[53] params[53] ]
     pi_d = @SMatrix [ params[54] 1-params[54]; 1-params[55] params[55] ]

     msp = ModelStructuralParams(Π_0 = Pi0, Π_X = PiX, m_g = params[9], m_r₁ = params[10],
     m_r₂ = params[11], ϕ = params[12], δ = params[13], μ_g = params[14], μ_pi = params[15], ρ₁ = params[16],
     ρ₂ = params[17], β₁ = beta1, β₂ = beta2, α₁ = alpha1, α₂ = alpha2, σ_g = params[22], σ_pi = params[23],
     σ_r₁ = params[24], σ_r₂ = params[25], σ_y = params[26], γ_gₐ = params[27], γ_piₐ = params[28],
     γ_cₐ = params[29], γ_gᵢ = params[30], γ_piᵢ = params[31], γ_cᵢ = params[32], γ_gₒ = params[33],
     γ_piₒ = params[34], γ_cₒ = params[35], ρₐ = params[36], ρᵢ = params[37], ρₒ = params[38], σ_wₐ = params[39],
     σ_zₐ = params[40], σ_νₐ = params[41], σ_qₐ = params[42], σ_wᵢ = params[43], σ_zᵢ = params[44],
     σ_νᵢ = params[45], σ_qᵢ = params[46], σ_wₒ = params[47], σ_zₒ = params[48], σ_νₒ = params[49],
     σ_qₒ = params[50], λ = params[51], π_m = pi_m, π_d = pi_d)

     return msp

end

function params_to_rf(params)

     msp = construct_structural_parameters(params)

     @unpack Π_0, Π_X, m_g, m_r₁, m_r₂, ϕ, δ, μ_g, μ_pi, ρ₁, ρ₂, β₁, β₂, α₁, α₂, σ_g, σ_pi, σ_r₁, σ_r₂, σ_y, γ_gₐ,
     γ_piₐ, γ_cₐ, γ_gᵢ, γ_piᵢ, γ_cᵢ, γ_gₒ, γ_piₒ, γ_cₒ, ρₐ, ρᵢ, ρₒ, σ_wₐ, σ_zₐ, σ_νₐ, σ_qₐ, σ_wᵢ, σ_zᵢ, σ_νᵢ, σ_qᵢ, σ_wₒ, σ_zₒ,
     σ_νₒ, σ_qₒ, λ, π_m, π_d = msp

     b0₁,b0₂,bm1₁,bm1₂,b1₁,b1₂,m0₁,m0₂,Γ₁,Γ₂ =
          construct_structural_matrices_macro(ϕ, δ, ρ₁, ρ₂, β₁, β₂, α₁, α₂, μ_g, μ_pi, m_g, m_r₁, m_r₂, σ_g, σ_pi, σ_r₁,
                                                                       σ_r₂)

     μ, Φ, μ_Q, Φ_Q, Σ, cov, Π, q, det_dum =
          transform_struct_to_rf(b0₁, b0₂, bm1₁, bm1₂, b1₁, b1₂, m0₁, m0₂, Γ₁, Γ₂, π_m, π_d, Π_0, Π_X)

     μ_ν , Φ_ν, μ_Q_ν, Φ_Q_ν, Σ_ν, cov_ν, Δ, Δ_Q, σ_w, σ_z, σ_q, σ_ν =
          augment_macro_fsmsre_nu(μ, Φ, μ_Q, Φ_Q, Σ, γ_gₐ, γ_gᵢ, γ_gₒ, γ_piₐ, γ_piᵢ, γ_piₒ, γ_cₐ, γ_cᵢ, γ_cₒ, ρₐ, ρᵢ, ρₒ,
                                                        λ, σ_wₐ, σ_wᵢ, σ_wₒ, σ_zₐ, σ_zᵢ, σ_zₒ,σ_qₐ, σ_qᵢ, σ_qₒ, σ_νₐ, σ_νᵢ, σ_νₒ)

     mrfp = ModelReducedFormParams(μ = μ, μ_ν = μ_ν, μ_Q = μ_Q, μ_Q_ν = μ_Q_ν, Φ = Φ, Φ_ν = Φ_ν,
                                                               Φ_Q = Φ_Q, Φ_Q_ν = Φ_Q_ν, Σ = Σ, Σ_ν = Σ_ν, cov = cov,
                                                               cov_ν = cov_ν, Π = Π, q = q, det_dum = det_dum, Δ = Δ, Δ_Q = Δ_Q,
                                                               σ_w = σ_w, σ_z = σ_z, σ_q = σ_q, σ_ν = σ_ν)

     return mrfp
end

function augment_macro_fsmsre_nu(μ, Φ, μ_Q, Φ_Q, Σ, γ_gₐ, γ_gᵢ, γ_gₒ, γ_piₐ, γ_piᵢ, γ_piₒ, γ_cₐ, γ_cᵢ, γ_cₒ, ρₐ, ρᵢ,
                                                            ρₒ, λ, σ_wₐ, σ_wᵢ, σ_wₒ, σ_zₐ, σ_zᵢ, σ_zₒ, σ_qₐ, σ_qᵢ, σ_qₒ, σ_νₐ, σ_νᵢ, σ_νₒ)

     S = size(μ)[1]
     N = size(μ[1])[1]
     Δ = @SMatrix [γ_gₐ γ_piₐ 0.; γ_gᵢ γ_piᵢ 0.; γ_gₒ γ_piₒ 0. ]
     ρ = @SMatrix [ρₐ 0. 0.; 0. ρᵢ 0.; 0. 0. ρₒ]
     C = SMatrix{3,1}([γ_cₐ, γ_cᵢ, γ_cₒ])
     σ_z = SMatrix{3,1}([σ_zₐ, σ_zᵢ, σ_zₒ])
     C_Q = C - λ.*σ_z
     σ2_ν = @SMatrix [(σ_wₐ^2 + σ_zₐ^2) σ_zₐ*σ_zᵢ σ_zₐ*σ_zₒ;
                                    σ_zₐ*σ_zᵢ (σ_wᵢ^2 + σ_zᵢ^2) σ_zᵢ*σ_zₒ;
                                    σ_zₐ*σ_zₒ  σ_zᵢ*σ_zₒ (σ_wₒ^2 + σ_zₒ^2)]
     σ_ν = cholesky(σ2_ν).L

     μ_ν = [ [μ[1]; Δ*μ[1]+C], [μ[2]; Δ*μ[2]+C], [μ[3]; Δ*μ[3]+C], [μ[4]; Δ*μ[4]+C] ]

     μ_Q_ν = [ [μ_Q[1]; Δ*μ_Q[1]+C_Q], [μ_Q[2]; Δ*μ_Q[2]+C_Q], [μ_Q[3]; Δ*μ_Q[3]+C_Q],
                       [μ_Q[4];Δ*μ_Q[4]+C_Q] ]

     Φ_ν = [ SMatrix{6,6,Float64,36}([Φ[1] zeros(N,N); Δ*Φ[1] ρ]),
                  SMatrix{6,6,Float64,36}([Φ[2] zeros(N,N); Δ*Φ[2] ρ]),
                  SMatrix{6,6,Float64,36}( [Φ[3] zeros(N,N); Δ*Φ[3] ρ]),
                  SMatrix{6,6,Float64,36}([Φ[4] zeros(N,N); Δ*Φ[4] ρ])]

     Φ_Q_ν = [ SMatrix{6,6,Float64,36}([Φ_Q[1] zeros(N,N); Δ*Φ_Q[1] ρ]),
                       SMatrix{6,6,Float64,36}([Φ_Q[2] zeros(N,N); Δ*Φ_Q[2] ρ]),
                       SMatrix{6,6,Float64,36}([Φ_Q[3] zeros(N,N); Δ*Φ_Q[3] ρ]),
                       SMatrix{6,6,Float64,36}([Φ_Q[4] zeros(N,N); Δ*Φ_Q[4] ρ])]

     Σ_ν = [ SMatrix{6,6,Float64,36}([Σ[1] zeros(N,N); Δ*Σ[1] σ_ν]),
                 SMatrix{6,6,Float64,36}( [Σ[2] zeros(N,N); Δ*Σ[2] σ_ν]),
                 SMatrix{6,6,Float64,36}([Σ[3] zeros(N,N); Δ*Σ[3] σ_ν]),
                 SMatrix{6,6,Float64,36}([Σ[4] zeros(N,N); Δ*Σ[4] σ_ν])]

     cov_ν = [Σ_ν[1]*(Σ_ν[1]'),Σ_ν[2]*(Σ_ν[2]'),Σ_ν[3]*(Σ_ν[3]'),Σ_ν[4]*(Σ_ν[4]')]

     Δ_tmp = @SMatrix [γ_gₐ γ_piₐ 0. ρₐ 0. 0. γ_cₐ;
                                       γ_gᵢ γ_piᵢ 0. 0. ρᵢ 0. γ_cᵢ;
                                       γ_gₒ γ_piₒ 0. 0. 0. ρₒ γ_cₒ];

     Δ_Q = @SMatrix [γ_gₐ γ_piₐ 0. ρₐ 0. 0. (γ_cₐ-σ_zₐ*λ);
                                   γ_gᵢ γ_piᵢ 0. 0. ρᵢ 0. (γ_cᵢ-σ_zᵢ*λ);
                                   γ_gₒ γ_piₒ 0. 0. 0. ρₒ (γ_cₒ-σ_zₒ*λ) ];

    σ_w = SMatrix{3, 1, Float64, 3}([σ_wₐ; σ_wᵢ; σ_wₒ])
    σ_z = SMatrix{3, 1, Float64, 3}([σ_zₐ; σ_zᵢ; σ_zₒ])
    σ_q = SMatrix{3, 1, Float64, 3}([σ_qₐ; σ_qᵢ; σ_qₒ])
    σ_ν = SMatrix{3, 1, Float64, 3}([σ_νₐ; σ_νᵢ; σ_νₒ])

     return μ_ν , Φ_ν, μ_Q_ν, Φ_Q_ν, Σ_ν, cov_ν, Δ_tmp, Δ_Q, σ_w, σ_z, σ_q, σ_ν
end

function split_shock_mat(u_mat, ϵ_mat, w_mat, z_mat, g)

     gk = size(u_mat)[2]
     #T = size(u_mat)[1]
     k = Int(gk / g)

     u = Vector{Array{Float64,2}}(undef, g)
     ϵ = Vector{Array{Float64,3}}(undef, g)
     w = Vector{Array{Float64,3}}(undef, g)
     z = Vector{Array{Float64,2}}(undef, g)

     for gg ∈ 1:g
          u[gg] = u_mat[:,1+(gg-1)*k:gg*k]
          ϵ[gg] = ϵ_mat[:,:,1+(gg-1)*k:gg*k]
          w[gg] = w_mat[:,:,1+(gg-1)*k:gg*k]
          z[gg] = z_mat[:,1+(gg-1)*k:gg*k]
     end

     return u, ϵ, w, z
end

function simulate_nu_cond_x_i_shock_rn!(ν_mat,x_mat,init_ν,mrfp,w,z,T_start, g_start)

     N = size(x_mat)[1]
     n_ts = size(w)[2]
     G = size(w)[3]
     N_re, N2 = size(init_ν)
     @unpack σ_w, σ_z, Δ_Q = mrfp

     if N2==1
          init_ν = repeat(init_ν,1,G)
     end

     ## Scale systematic shocks by std. dev. loading
     w_mat = repeat(σ_w, 1, n_ts, G) .* w

     tmp_z_mat = Array{Float64}(undef, N_re, n_ts, G)
     for g ∈ 1:G
          tmp_z_mat[:, :, g] = repeat(z[:,g]', N_re)
     end

     z_mat = repeat(σ_z, 1, n_ts, G) .* tmp_z_mat
     wz_mat = w_mat + z_mat

     cnt_g = 1
     for g=g_start:G+g_start-1
          if T_start == 1
               ν_mat[:,1,g] = @view(init_ν[:, cnt_g])
          end
          wz_path = @view(wz_mat[:, :, cnt_g])
          ##Create X_mat = [x_t; nu_{t-1}; 1]'; an (N + N_re + 1) x n_ts Array
          X_mat = Array{Float64}(undef, N+N_re+1, n_ts+1)
          X_mat[1:N, 1:end-1] = @view(x_mat[:, T_start+1:T_start+n_ts, g])
          X_mat[end,:] = ones(1, n_ts+1)
          X_mat[N+1:N+N_re, 1] = @view(ν_mat[:, 1, g])

          ##Compute nu_t iteratively
          cnt = 1
          for i=T_start:n_ts+T_start-1

               X_mat[N+1:N+N_re, cnt+1] = Δ_Q * @view(X_mat[:, cnt]) + @view(wz_path[: ,cnt])
               cnt +=1

          end

          ν_mat[:, T_start+1:T_start+n_ts, g] = @view(X_mat[N+1:N+N_re, 2:end])
          cnt_g += 1

     end

     return ν_mat
end

function compute_mc_real_estate_Q_cond_x_i_nofull!( Q_mc, η_mc, Q_std_mc, η_std_mc,
                                                                                Q_std_prep, η_std_prep,x_mat, ν_mat,T_bar,G)
     T = size(x_mat)[2] - 1
     n_η = T_bar
     N_re = size(ν_mat)[1]

     # Preallocate
     Δ_ν = [0. 0. -1. 1. 0. 0.;
                 0. 0. -1. 0. 1. 0.;
                 0. 0. -1. 0. 0. 1.]

     m1_η_sim = Array{Float64}(undef, N_re, T_bar, G)
     x_path_mat = [x_mat; ν_mat]

     for g=1:G

          x_path = @view(x_path_mat[:,1:T_bar,g])
          m1_η_sim[:, :, g] = cumsum(Δ_ν * x_path, dims = 2)

     end
     η_sim = exp.(m1_η_sim)
     η_mc[:,1:T] .= dropdims(mean(η_sim[:,:,1:G], dims = 3), dims = 3)
     Q_mc .= sum(η_mc[:,1:T], dims = 2)
     #Get std mat
     N_std = Int64(floor(G / 30))
     for i=1:30

          η_std_prep[:, 1:T_bar, i] =  dropdims(mean(@view(η_sim[:, :, 1+(i-1)*N_std:i*N_std]), dims = 3), dims = 3)
          Q_std_prep[:, i] = sum(@view(η_std_prep[:, 1:n_η, i]), dims = 2)

     end

     η_std_mc .= dropdims(std(η_std_prep, dims = 3), dims = 3) ./ sqrt(30)
     Q_std_mc .= std(Q_std_prep, dims = 2) ./ sqrt(30)

end

function compute_mc_term_structure_cond_x(x_mat, maturity_mat)
     G = size(x_mat)[3]
     N_terms = size(maturity_mat)[2]

     δ = [0. 0. -1.]
     ts_sim = Array{Float64}(undef, N_terms, G)
     m1_ts_sim = Array{Float64}(undef, N_terms,G)
     m2_ts_sim = Array{Float64}(undef, N_terms,G)

     for g=1:G
          x_path = @view(x_mat[:, :, g])
          for n=1:N_terms
             m1_ts_sim[n,g] = sum(δ *@view(x_path[:, 1:maturity_mat[n]]), dims = 2)[1]
          end
     end
     m2_ts_sim = m1_ts_sim.^2
     ts_sim = exp.(m1_ts_sim)

     ts_mc = -log.(mean(ts_sim, dims = 2)) ./ maturity_mat'
     m1_ts_mc = mean(m1_ts_sim, dims = 2);
     m2_ts_mc = mean(m2_ts_sim, dims = 2);

     N_std = Int64(floor(G/30))
     ts_std_prep = Array{Float64}(undef, N_terms, 30)
     for i=1:30
          ts_std_prep[:, i] =  -log.(mean(@view(ts_sim[:, 1+(i-1)*N_std:i*N_std]), dims = 2)) ./ maturity_mat'
     end
     ts_std_mc = std(ts_std_prep, dims = 2) ./ sqrt(30)

     return ts_mc, ts_std_mc, m1_ts_mc, m2_ts_mc, ts_std_prep
end

function   get_total_mean_and_std!(Q_mc_tot, η_mc_tot, Q_std_mc_tot, Q_mc_cell, eta_mc_cell,
                                                            Q_std_prep_cell)

     Q_mc_tot .= compute_mean_over_subsamples(Q_mc_cell)
     η_mc_tot .= compute_mean_over_subsamples(eta_mc_cell)

     Q_std_mc_tot .= compute_std_over_subsamples(Q_std_prep_cell)
     #η_std_mc_tot .= compute_std_over_subsamples(η_std_prep_cell)

end

function compute_std_over_subsamples(in_cell_std)
     L = size(in_cell_std)[1]
     N1, N2 = size(in_cell_std[1])

     if L>1
          in_std_mat = dropdims(sum(reshape(reduce(hcat, in_cell_std), N1, N2, L), dims = 3) ./ L, dims = 3)
      else
          in_std_mat = in_cell_std[1]
      end

     b = size(in_std_mat)[2]
     out_std_mat = std(in_std_mat, dims = 2) ./ sqrt(b)

     return out_std_mat

end

function compute_mean_over_subsamples(in_cell)

     n,m = size(in_cell[1])
     L = size(in_cell)[1]

     if m==1
         mean_mat = mean(reduce(hcat,in_cell), dims = 2)
     else
         mean_mat =  dropdims(sum(reshape(reduce(hcat, in_cell), n, m, L), dims = 3) ./ L, dims = 3)
     end

     return mean_mat

end

function simulate_model_prices_cond_shock_acc_ts(x_init, ν_init, mrfp, T_sim, u, ϵ, w, z, tol, del, η_tol, t, s,
                                                                                n_grps, L)
     # Unload inputs
     N_macro = 3
     N_re = 3
     M = size(u)[1]
     G = size(ϵ)[3]
     @assert floor(G / L)== (G/L)
     gg = Int64(G / L)

     @unpack μ_Q, Φ_Q, Σ, Π = mrfp

     Q_mc = Array{Float64}(undef, N_re, 1)
     Q_std_mc = Array{Float64}(undef, N_re, 1)
     #η_std_mc = Array{Float64}(undef, N_re, T_sim)
     η_mc= Array{Float64}(undef, N_re, T_sim)
     x_mat = Array{Float64}(undef, N_macro, T_sim+1, G)
     ν_mat = Array{Float64}(undef, N_re, T_sim + 1, G)
     s_mat = Array{Int64}(undef, T_sim+1, G)

     u_cell, ϵ_cell, w_cell, z_cell  = split_shock_mat(u, ϵ, w, z, L)

     X = Array(Array(x_init[:, t]')')
     ν = Array(Array(ν_init[:, t]')')
     T_bar = [0]

     outer_while!(Q_mc, η_mc, Q_std_mc, x_mat, ν_mat, s_mat, T_bar, T_sim, n_grps, gg, X, ν, s, μ_Q,
                                    Φ_Q, Σ, Π, u_cell, ϵ_cell, w_cell, z_cell, mrfp, η_tol, N_re, del, L, tol)

     return Q_mc, η_mc, Q_std_mc, s_mat, x_mat, ν_mat, T_bar

end

function inner_while!(x_mat,s_mat, ν_mat, T_bar, Q_mc_tmp, η_mc_tmp, Q_std_mc_tmp, η_std_mc_tmp,
                                   Q_std_prep, η_std_prep, l, n_grps, T_tmp, μ_Q, Φ_Q, Σ, Π, u_cell, ϵ_cell, w_cell,z_cell,
                                   mrfp, η_tol,gg)

     trig_2 = 0
     cnt_2 = 1
     while (trig_2==0)&&(cnt_2<=n_grps)

          tmp_ndx = (1+(cnt_2-1)*T_tmp):(cnt_2*T_tmp)
          g_start = 1 + (l-1)*gg
          #Simulate model for length T_tmp
          simulate_ms_var_1_cond_shocks!(x_mat, s_mat, x_mat[:, tmp_ndx[1], :],
               s_mat[tmp_ndx[1], :], μ_Q, Φ_Q, Σ, Π, u_cell[l][tmp_ndx, :], ϵ_cell[l][: ,tmp_ndx, :],tmp_ndx[1],g_start)

          #Simulate \nu
          simulate_nu_cond_x_i_shock_rn!(ν_mat, x_mat, ν_mat[:, tmp_ndx[1], :], mrfp, w_cell[l][:, tmp_ndx, :],
                                                                 z_cell[l][tmp_ndx, :],tmp_ndx[1],g_start)

          #Simulate Q
          compute_mc_real_estate_Q_cond_x_i_nofull!(Q_mc_tmp, η_mc_tmp, Q_std_mc_tmp, η_std_mc_tmp,
                                                                                Q_std_prep, η_std_prep , x_mat[:, 1:tmp_ndx[end]+1, :],
                                                                                ν_mat[:, 1:tmp_ndx[end]+1, :],tmp_ndx[end],l*gg)
          #Check that all η_{T_sim,x_t)< cutoff for all t
          if (~any(η_mc_tmp[:, tmp_ndx[end]] .>= η_tol)) || (cnt_2==n_grps)

               #Set loop trigger and T_bar
               T_bar[1]  = cnt_2 * T_tmp
               trig_2 = 1

          else

               cnt_2 +=  1

          end
     end

     return

end

function outer_while!(Q_mc, η_mc, Q_std_mc, x_mat, ν_mat, s_mat, T_bar, T_sim, n_grps, gg, X, ν, s,
                                        μ_Q, Φ_Q, Σ, Π, u_cell, ϵ_cell, w_cell, z_cell, mrfp, η_tol, N_re, del, L, tol)

     T_tmp = Int64(T_sim ./ n_grps)
     G = gg * L
     x_mat[:, 1, :] = repeat(X, 1, 1, G)
     ν_mat[:, 1, :] = repeat(ν, 1, 1, G)
     s_mat[1, :] .= s

     trig = 0
     l = 1

     η_std_prep = Array{Float64}(undef, N_re, T_sim, 30)
     Q_std_prep = Array{Float64}(undef, N_re, 30)
     Q_mc_tmp = Array{Float64}(undef, N_re, 1)
     η_mc_tmp = Array{Float64}(undef, N_re, T_sim)
     Q_std_mc_tmp = Array{Float64}(undef, N_re, 1)
     η_std_mc_tmp = Array{Float64}(undef, N_re, T_sim)
     Q_mc_tot = Array{Float64}(undef, N_re,1)
     η_mc_tot = Array{Float64}(undef, N_re, T_sim)
     Q_std_mc_tot = Array{Float64}(undef, N_re, 1)
     #η_std_mc_tot = Array{Float64}(undef, N_re, T_sim)
     Q_mc_cell = Vector{Array{Float64,2}}(undef, L)
     Q_std_mc_cell = Vector{Array{Float64,2}}(undef, L)
     η_mc_cell = Vector{Array{Float64,2}}(undef, L)
     η_std_mc_cell = Vector{Array{Float64,2}}(undef, L)
     Q_std_prep_cell = Vector{Array{Float64,2}}(undef, L)
     η_std_prep_cell = Vector{Array{Float64,3}}(undef, L)

     while ((trig==0) && (l<=L)) #early simulation check

          if l==1

               #Detect early truncation point, T_bar
               inner_while!(x_mat, s_mat, ν_mat, T_bar, Q_mc_tmp, η_mc_tmp, Q_std_mc_tmp, η_std_mc_tmp,
                                        Q_std_prep, η_std_prep, l, n_grps, T_tmp, μ_Q, Φ_Q, Σ, Π, u_cell, ϵ_cell, w_cell, z_cell,
                                         mrfp, η_tol,gg)
               TT_bar = T_bar[1]

          else
               TT_bar = T_bar[1]
               g_start = 1 + (l-1)*gg
               #Run simulation for T_bar periods
               # Simulate Macro model
               simulate_ms_var_1_cond_shocks!(x_mat, s_mat,X, s, μ_Q, Φ_Q, Σ,Π,u_cell[l][1:TT_bar, :],
                                                                                ϵ_cell[l][:, 1:TT_bar, :],1,g_start)

               #Simulate \nu
               simulate_nu_cond_x_i_shock_rn!(ν_mat, x_mat, ν, mrfp ,w_cell[l][:, 1:TT_bar, :], z_cell[l][1:TT_bar, :],
                                                                      1, g_start)

               #Approximate real estate Q
               compute_mc_real_estate_Q_cond_x_i_nofull!(Q_mc_tmp, η_mc_tmp, Q_std_mc_tmp, η_std_mc_tmp,
                                                                                Q_std_prep, η_std_prep, x_mat[:,1:TT_bar+1,:],
                                                                                ν_mat[:,1:TT_bar+1,:],TT_bar, l * gg)
          end

          Q_mc_cell[l] = copy(Q_mc_tmp)
          η_mc_cell[l] = copy(η_mc_tmp)
          Q_std_mc_cell[l] = copy(Q_std_mc_tmp)
          η_std_mc_cell[l] = copy(η_std_mc_tmp)
          Q_std_prep_cell[l] = copy(Q_std_prep)
          η_std_prep_cell[l] = copy(η_std_prep)

          #Compute total mean and variance
          get_total_mean_and_std!(Q_mc_tot, η_mc_tot, Q_std_mc_tot, Q_mc_cell[1:l], η_mc_cell[1:l],
                                                    Q_std_prep_cell[1:l])

          #Check condition
          if (!any((2 .* (1. .- normcdf.((1.2 .* tol .* Q_mc_tot) ./ Q_std_mc_tot)) .> del)) || (l==L) )

               trig = 1;
               println("l: $l")
               println("T_bar: $T_bar")

               Q_mc .= Q_mc_tot
               η_mc .= η_mc_tot
               Q_std_mc .= Q_std_mc_tot
               #η_std_mc[:, 1:T_bar] = η_std_mc_tot

          else

              l = l+1;

          end
     end

     return

end

export simulate_markov_switch_init_cond_shock, simulate_ms_var_1_cond_shocks,
     simulate_msvar_cond_regime_path_shock, construct_gamma_macro_array, construct_m0_macro_array,
     construct_b1_macro_array, construct_bm1_macro_array, construct_b0_macro_array,
     construct_structural_matrices_macro, transform_struct_to_rf, get_forward_solution_msre, get_rf_Φ,
     compute_drift_sigma_fmmsre, compute_μ_Σ_fmmsre, construct_structural_parameters, params_to_rf,
     augment_macro_fsmsre_nu, simulate_nu_cond_x_i_shock_rn, compute_mc_real_estate_Q_cond_x_i_nofull,
     compute_mc_term_structure_cond_x, compute_std_over_subsamples, compute_mean_over_subsamples,
     simulate_model_prices_cond_shock_acc_ts, split_shock_mat

end


text:
#
# DataStream() should just convert things to appropriate
# type of DataStream, each of which implements the AbstractDataStream
# protocol. At present this can be done using a file on disk through
# the FileDataStream type or using an in-memory DataFrame through the
# the DataFrameDataStream type.
#

abstract AbstractDataStream

##############################################################################
#
# FileDataStream
#
##############################################################################

# TODO: Decide whether to keep column types. Reader no longer supports
# explicit a priori type information.
type FileDataStream <: AbstractDataStream
  filename::String
  stream::IO # Any
  separator::Char
  quotation_character::Char
  missingness_indicators::Vector
  header::Bool
  column_names::Vector
  column_types::Vector
  minibatch_size::Int
end

function FileDataStream{T <: String}(filename::T, minibatch_size::Int)
  stream = OUTPUT_STREAM
  separator = determine_separator(filename)
  quotation_character = DEFAULT_QUOTATION_CHARACTER
  missingness_indicators = DEFAULT_MISSINGNESS_INDICATORS
  header = true
  short_circuit = true
  io = open(filename, "r")
  column_names = determine_column_names(io,
                                        separator,
                                        quotation_character,
                                        header)
  close(io)
  column_types = {Any for i in 1:length(column_names)}
  FileDataStream(filename, stream, separator,
                 quotation_character, missingness_indicators,
                 header, column_names, column_types,
                 minibatch_size)
end

function FileDataStream{T <: String}(filename::T)
  FileDataStream(filename, 1)
end

function DataStream{T <: String}(filename::T, minibatch_size::Int)
  FileDataStream(filename, minibatch_size)
end

function DataStream{T <: String}(filename::T)
  FileDataStream(filename, 1)
end

function start(ds::FileDataStream)
  ds.stream = open(ds.filename, "r")

  # Read one line to remove header in advance
  if ds.header
    readline(ds.stream)
  end

  next_df = read_minibatch(ds.stream,
                           ds.separator,
                           ds.quotation_character,
                           ds.missingness_indicators,
                           colnames(ds),
                           ds.minibatch_size)

  return next_df
end

function next(ds::FileDataStream, df::DataFrame)
  next_df = read_minibatch(ds.stream,
                           ds.separator,
                           ds.quotation_character,
                           ds.missingness_indicators,
                           colnames(ds),
                           ds.minibatch_size)
  (df, next_df)
end

function done(ds::FileDataStream, df::DataFrame)
  if nrow(df) == 0
    close(ds.stream)
    return true
  else
    return false
  end
end

coltypes(ds::FileDataStream) = ds.column_types
colnames(ds::FileDataStream) = ds.column_names

##############################################################################
#
# IODataStream
#
##############################################################################

type IODataStream <: AbstractDataStream
  stream::IO
  separator::Char
  quotation_character::Char
  missingness_indicators::Vector
  header::Bool
  column_names::Vector
  minibatch_size::Int
end

function IODataStream(io::IO, minibatch_size::Int)
  stream = io
  separator = DEFAULT_SEPARATOR
  quotation_character = DEFAULT_QUOTATION_CHARACTER
  missingness_indicators = DEFAULT_MISSINGNESS_INDICATORS
  header = true
  column_names = ["FILLER"]
  short_circuit = true
  IODataStream(stream, separator,
               quotation_character, missingness_indicators,
               header, column_names, minibatch_size)
end

function IODataStream(io::IO)
  IODataStream(io, 1)
end

function DataStream(io::IO, minibatch_size::Int)
  IODataStream(io, minibatch_size)
end

function DataStream(io::IO)
  IODataStream(io, 1)
end

function start(ds::IODataStream)
  # Seek to the start of the stream
  # May be impossible for some streams like STDIN
  try
    seek(ds.stream, 0)
  catch
    warn("Could not seek to the start of the input stream")
  end

  # Read one line to remove header in advance
  if ds.header
    line = chomp(readline(ds.stream))
    if strlen(line) == 0
      error("Empty header line")
    end
    items = Array(UTF8String, strlen(line))
    current_field = Array(Char, strlen(line))
    ds.column_names = split_separated_line(line, ds.separator, ds.quotation_character, items, current_field)
  else
    error("Currently only IODataStream's with headers are supported")
  end

  df = read_minibatch(ds.stream,
                      ds.separator,
                      ds.quotation_character,
                      ds.missingness_indicators,
                      colnames(ds),
                      ds.minibatch_size)

  return df
end

function next(ds::IODataStream, df::DataFrame)
  next_df = read_minibatch(ds.stream,
                           ds.separator,
                           ds.quotation_character,
                           ds.missingness_indicators,
                           colnames(ds),
                           ds.minibatch_size)
  (df, next_df)
end

function done(ds::IODataStream, df::DataFrame)
  if nrow(df) == 0
    return true
  else
    return false
  end
end

coltypes(ds::IODataStream) = ds.column_types
colnames(ds::IODataStream) = ds.column_names

##############################################################################
#
# DataFrameDataStream
#
##############################################################################

type DataFrameDataStream <: AbstractDataStream
  df::DataFrame
  minibatch_size::Int
end

function start(ds::DataFrameDataStream)
  return 1
end

function next(ds::DataFrameDataStream, i::Int)
  (ds.df[i:(i + ds.minibatch_size - 1), :], i + ds.minibatch_size)
end

function done(ds::DataFrameDataStream, i::Int)
  return i > nrow(ds.df)
end

DataStream(df::DataFrame, minibatch_size::Int) = DataFrameDataStream(df, minibatch_size)

DataStream(df::DataFrame) = DataFrameDataStream(df, 1)

coltypes(ds::DataFrameDataStream) = coltypes(ds.df)
colnames(ds::DataFrameDataStream) = colnames(ds.df)

##############################################################################
#
# MatrixDataStream
#
##############################################################################

type MatrixDataStream{T} <: AbstractDataStream
  m::Matrix{T}
  minibatch_size::Int
end

function start(ds::MatrixDataStream)
  return 1
end

function next(ds::MatrixDataStream, i::Int)
  (DataFrame(ds.m[i:(i + ds.minibatch_size - 1), :]), i + ds.minibatch_size)
end

function done(ds::MatrixDataStream, i::Int)
  return i > size(ds.m, 1)
end

DataStream{T}(m::Matrix{T}, minibatch_size::Int) = MatrixDataStream{T}(m, minibatch_size)

DataStream{T}(m::Matrix{T}) = MatrixDataStream(m, 1)

coltypes{T}(ds::MatrixDataStream{T}) = {T for i in 1:size(ds.m, 2)}
colnames(ds::MatrixDataStream) = generate_column_names(size(ds.m, 2))

##############################################################################
#
# Streaming data functions
#
##############################################################################

function colsums(ds::AbstractDataStream)
  p = length(colnames(ds))
  sums = zeros(p)
  ns = zeros(Int, p)

  for minibatch in ds
    for row_index in 1:nrow(minibatch)
      for column_index in 1:p
        if coltypes(minibatch)[column_index] <: Real && !isna(minibatch[row_index, column_index])
          sums[column_index] += minibatch[row_index, column_index]
          ns[column_index] += 1
        end
      end
    end
  end

  result_types = {Float64 for i in 1:p}
  results = DataFrame(result_types, colnames(ds), 1)

  for column_index in 1:p
    if ns[column_index] != 0
      results[1, column_index] = sums[column_index]
    end
  end

  return results
end

function colprods(ds::AbstractDataStream)
  p = length(colnames(ds))
  prods = zeros(p)
  ns = ones(Int, p)

  for minibatch in ds
    for row_index in 1:nrow(minibatch)
      for column_index in 1:p
        if coltypes(minibatch)[column_index] <: Real && !isna(minibatch[row_index, column_index])
          prods[column_index] *= minibatch[row_index, column_index]
          ns[column_index] += 1
        end
      end
    end
  end

  result_types = {Float64 for i in 1:p}
  results = DataFrame(result_types, colnames(ds), 1)

  for column_index in 1:p
    if ns[column_index] != 0
      results[1, column_index] = prods[column_index]
    end
  end

  return results
end

function colmeans(ds::AbstractDataStream)
  p = length(colnames(ds))
  sums = zeros(p)
  ns = zeros(Int, p)

  for minibatch in ds
    for row_index in 1:nrow(minibatch)
      for column_index in 1:p
        if coltypes(minibatch)[column_index] <: Real && !isna(minibatch[row_index, column_index])
          sums[column_index] += minibatch[row_index, column_index]
          ns[column_index] += 1
        end
      end
    end
  end

  result_types = {Float64 for i in 1:p}
  results = DataFrame(result_types, colnames(ds), 1)

  for column_index in 1:p
    if ns[column_index] != 0
      results[1, column_index] = sums[column_index] / ns[column_index]
    end
  end

  return results
end

function colvars(ds::AbstractDataStream)
  p = length(colnames(ds))
  means = zeros(p)
  deltas = zeros(p)
  m2s = zeros(p)
  vars = zeros(p)
  ns = zeros(Int, p)

  for minibatch in ds
    for row_index in 1:nrow(minibatch)
      for column_index in 1:p
        if coltypes(minibatch)[column_index] <: Real && !isna(minibatch[row_index, column_index])
          ns[column_index] += 1
          deltas[column_index] = minibatch[row_index, column_index] - means[column_index]
          means[column_index] += deltas[column_index] / ns[column_index]
          m2s[column_index] = m2s[column_index] + deltas[column_index] * (minibatch[row_index, column_index] - means[column_index])
          vars[column_index] = m2s[column_index] / (ns[column_index] - 1)
        end
      end
    end
  end

  result_types = {Float64 for i in 1:p}
  results = DataFrame(result_types, colnames(ds), 1)

  for column_index in 1:p
    if ns[column_index] != 0
      results[1, column_index] = vars[column_index]
    end
  end

  return results
end

function colstds(ds::AbstractDataStream)
  vars = colvars(ds)
  stds = deepcopy(vars)
  column_types = coltypes(vars)
  for j in 1:length(column_types)
    if column_types[j] <: Real
      stds[1, j] = sqrt(vars[1, j])
    end
  end
  return stds
end

function colmins(ds::AbstractDataStream)
  p = length(colnames(ds))
  mins = [Inf for i in 1:p]
  ns = zeros(Int, p)

  for minibatch in ds
    for row_index in 1:nrow(minibatch)
      for column_index in 1:p
        if coltypes(minibatch)[column_index] <: Real && !isna(minibatch[row_index, column_index])
          if minibatch[row_index, column_index] < mins[column_index]
            mins[column_index] = minibatch[row_index, column_index]
            ns[column_index] += 1
          end
        end
      end
    end
  end

  result_types = {Float64 for i in 1:p}
  df = DataFrame(result_types, colnames(ds), 1)

  for column_index in 1:p
    if ns[column_index] != 0
      df[1, column_index] = mins[column_index]
    end
  end

  return df
end

function colmaxs(ds::AbstractDataStream)
  p = length(colnames(ds))
  maxs = [-Inf for i in 1:p]
  ns = zeros(Int, p)

  for minibatch in ds
    for row_index in 1:nrow(minibatch)
      for column_index in 1:p
        if coltypes(minibatch)[column_index] <: Real && !isna(minibatch[row_index, column_index])
          if minibatch[row_index, column_index] > maxs[column_index]
            maxs[column_index] = minibatch[row_index, column_index]
            ns[column_index] += 1
          end
        end
      end
    end
  end

  result_types = {Float64 for i in 1:p}
  df = DataFrame(result_types, colnames(ds), 1)

  for column_index in 1:p
    if ns[column_index] != 0
      df[1, column_index] = maxs[column_index]
    end
  end

  return df
end

function colranges(ds::AbstractDataStream)
  p = length(colnames(ds))
  mins = [Inf for i in 1:p]
  maxs = [-Inf for i in 1:p]
  ns = zeros(Int, p)

  for minibatch in ds
    for row_index in 1:nrow(minibatch)
      for column_index in 1:p
        if coltypes(minibatch)[column_index] <: Real && !isna(minibatch[row_index, column_index])
          ns[column_index] += 1
          if minibatch[row_index, column_index] < mins[column_index]
            mins[column_index] = minibatch[row_index, column_index]
          end
          if minibatch[row_index, column_index] > maxs[column_index]
            maxs[column_index] = minibatch[row_index, column_index]
          end
        end
      end
    end
  end

  result_types = {Float64 for i in 1:p}
  df_mins = DataFrame(result_types, colnames(ds), 1)
  df_maxs = DataFrame(result_types, colnames(ds), 1)

  for column_index in 1:p
    if ns[column_index] != 0
      df_mins[1, column_index] = mins[column_index]
      df_maxs[1, column_index] = maxs[column_index]
    end
  end

  return (df_mins, df_maxs)
end

# Two-pass algorithm for covariance and correlation
function cov_pearson(ds::AbstractDataStream)
  p = length(colnames(ds))

  # Make one pass to compute means
  means = colmeans(ds)

  # Now compute covariances during second pass
  ns = zeros(Int, p, p)
  covariances = datazeros(p, p)

  for minibatch in ds
    for row_index in 1:nrow(minibatch)
      for column_index in 1:p
        for alt_column_index in 1:p
          if coltypes(minibatch)[column_index] <: Real &&
                !isna(minibatch[row_index, column_index]) &&
                coltypes(minibatch)[alt_column_index] <: Real &&
                !isna(minibatch[row_index, alt_column_index])
            ns[column_index, alt_column_index] += 1
            n = ns[column_index, alt_column_index]
            a = minibatch[row_index, column_index] - means[1, column_index]
            b = minibatch[row_index, alt_column_index] - means[1, alt_column_index]
            covariances[column_index, alt_column_index] = ((n - 1) / n) * covariances[column_index, alt_column_index] + (a * b) / n
          end
        end
      end
    end
  end

  # Scale estimates by (n / (n - 1))
  for i in 1:p
    for j in 1:p
      if ns[i, j] <= 2
        covariances[i, j] = NA
      else
        n = ns[i, j]
        covariances[i, j] *= (n / (n - 1))
      end
    end
  end

  return covariances
end

function cor_pearson(ds::AbstractDataStream)
  covariances = cov(ds)
  correlations = deepcopy(covariances)
  p = nrow(correlations)
  for i in 1:p
    for j in 1:p
      correlations[i, j] = covariances[i, j] / sqrt(covariances[i, i] * covariances[j, j])
    end
  end
  return correlations
end

function ref(ds::AbstractDataStream, i::Int)
  cur_i = 0
  for df in ds
    if cur_i + nrow(df) > i
      return df[i - cur_i, :]
    end
  end
  error("Did not find requested row")
end

# TODO: Stop returning empty DataFrame at the end of a stream
#       (NOTE: Probably not possible because we don't know nrows.)
# TODO: Implement
#        * colentropys
#        * colcardinalities
#        * colmedians
#        * colffts
#        * colnorms

ext: jl


text:
using PCLCommon
using PCLIO
using PCLVisualization
using Cxx
using OpenNI2

const ni2 = OpenNI2

const w = 640
const h = 480
const pxfmt = ni2.PIXEL_FORMAT_DEPTH_1_MM

function getPointXYZ(stream, deptharr)
    w, h = size(deptharr)
    depthvec = vec(deptharr)
    depthp = pointer(depthvec)
    cloud = PointCloud{PointXYZ}(w, h)
    icxx"$(cloud.handle)->is_dense = false;"
    pointsptr = icxx"&$(cloud.handle)->points[0];"
    icxx"""
    const float millimeterToMeter = 0.001f;
    for (int ri = 0; ri < $h; ++ri) {
        for (int ci = 0; ci < $w; ++ci) {
            auto p = $(pointsptr) + $w * ri + ci;
            auto z = $(depthp) + $w * ri + ci;
            openni::CoordinateConverter::convertDepthToWorld(
                *$(stream.handle), ci, ri, *z, &p->x, &p->y, &p->z);
            p->x *= millimeterToMeter;
            p->y *= millimeterToMeter;
            p->z *= millimeterToMeter;
        }
    }
    """
    cloud
end

ni2.initialize()
device = ni2.DevicePtr()
ni2.open(device)
ni2.setDepthColorSyncEnabled(device, false)

di = ni2.getDeviceInfo(device)
@show ni2.getName(di)
@show ni2.getVendor(di)

depth = ni2.VideoStreamPtr()
ni2.create(depth, device, ni2.SENSOR_DEPTH)
ni2.setMirroringEnabled(depth, true)

modes = ni2.getSupportedVideoModes(ni2.getSensorInfo(device, ni2.SENSOR_DEPTH))
for mode in modes
    if Int(ni2.getResolutionX(mode)) == w &&
            Int(ni2.getResolutionY(mode)) == h &&
            ni2.getPixelFormat(mode) == pxfmt
        println("$w x $h, PIXEL_FORMAT_DEPTH_1_MM: video mode found")
        ni2.setVideoMode(depth, mode)
        break
    end
end

ni2.start(depth)

frame = ni2.VideoFrameRef()

info("Prepare PCL visualizer...")
global viewer = PCLVisualizer("pcl visualizer")

while true
    @assert ni2.waitForAnyStream([depth]) == 0
    ni2.readFrame(depth, frame)
    arr = convert(Array{Cushort,2}, frame)
    @assert size(arr) == (w, h)

    cloud = getPointXYZ(depth, arr)
    if !updatePointCloud(viewer, cloud, id="depth")
        addPointCloud(viewer, cloud, id="depth")
    end
    spinOnce(viewer, 1)

    rand() > 0.98 && gc(false)
end

ni2.stop(depth)
ni2.destroy(depth)
ni2.close(device)
ni2.shutdown()

close(viewer)


text:
#!/usr/local/bin/julia
# ARGS=["753TD","753ND","tmp3.tsv","753_pair"]

using DataFrames
using CSV
using DelimitedFiles
using Statistics

tumor_id=ARGS[1]
normal_id=ARGS[2]
file1=ARGS[3]
pair_id=ARGS[4]
file1a=pair_id*".raw.tsv"

isString(x::Number)=false
isString(x::Missing)=false
isString(x::AbstractString)=true

df=CSV.File(file1,delim='\t') |> DataFrame

# size(df)
# describe(df)
# print(df)
FIELDS=names(df)
FIELDSK=[:CHRO, :POS, :REF, :ALT, :FILTER, :POPAF, :NLOD, :PON, :TLOD, :NORMAL_AD_REF, :NORMAL_AD_ALT, :TUMOR_AD_REF, :TUMOR_AD_ALT]

for i=1:length(FIELDS)
    if ~(FIELDS[i] in FIELDSK)
        deletecols!(df, FIELDS[i])
    end
end

df[:VCF_REF]=map(x->x,df[:REF])
df[:VCF_ALT]=map(x->x,df[:ALT])
df[:VCF_POS]=map(x->x,df[:POS])

for c in names(df)
    println(c)
    k=ismissing(df[c])
    if k == false
        continue
    end
    if Statistics.mean(k)==1.0
        deletecols!(df, c)
        continue
    end
    if isa(df[c],Array{String,1})
        df[k,c] = ""
    end
end
a= df[:CHRO]
if !isa(a[1],Int)
    a=map(x -> replace(x,r"[X]"=> "23"), a)
    a=map(x -> replace(x,r"[Y]"=> "24"), a)
    a=map(x -> replace(x,r"[MT]"=> "25"), a)
    a=map(x -> replace(x,r"[M]"=> "25"), a)
    a=map(x -> parse(Int32,x), a)
end
df[:a] = a
sort!(df, [:a, :POS])
deletecols!(df, [:a])

println("")
println(file1a)
println("")

open(file1a, "w") do f
    writedlm(f, reshape(names(df), 1, length(names(df))), '\t')
    writedlm(f, convert(Matrix,df), '\t')
end

a=df[:ALT]
r=df[:REF]
n=length(a)
ar=fill("",n)
for i in 1:n
    #println(i)
    ar[i]=string(a[i],":",r[i])
end

l1=map(x -> length(x), ar)
df[:INDEL]=map(x -> x>3,l1)
df[:Variant_Type]=fill("SNV",n)
df[df[:INDEL],:Variant_Type]="INDEL"
dfindel=df[df[:INDEL],:]
dfsnv=df[.~df[:INDEL],:]
df2=df

# set ALT0 to original ALT,shorten ALT to at most 5 bases, then copy back after merge

df2[:ALT0]=df2[:ALT]
alt0=df2[:ALT0]
alt=fill("",n)
n=length(alt0)
l0=map(x -> length(x), alt0)
for i in 1:n
    l1=min(l0[i],5)
    alt1=string(alt0[i])
    alt[i]=alt1[1:l1]
end
df2[:ALT]=alt

deletecols!(df, [:INDEL, :Variant_Type,:ALT0])

# label M1 with M1 flag
a= df[:CHRO]
if !isa(a[1],Int)
    a=map(x -> replace(x,r"[X]"=> "23"), a)
    a=map(x -> replace(x,r"[Y]"=> "24"), a)
    a=map(x -> replace(x,r"[MT]"=> "25"), a)
    a=map(x -> replace(x,r"[M]"=> "25"), a)
    a=map(x -> parse(Int32,x), a)
end
df[:a] = a
sort!(df, [:a, :POS])
deletecols!(df, [:a])

for c in names(df)
    println(c)
    #k=findall(map(x->ismissing(x),df[c]))
    k1=map(x->ismissing(x)*1,df[c])
    k=findall(k1.>0)
    k2=findall(k1.==0)
        #k=fomismissing(df[c])
    if Statistics.mean(k1)==1.0
        deletecols!(df, c)
        continue
    end
    println(typeof(df[k2[1],c]))
    if isString(df[k2[1],c])
        df[k,c] = ""
    end
end

# maflite fields
# SNP: build    chr start   end ref_allele  tum_allele1 tum_allele2 tumor_barcode   normal_barcode  tumor_f init_t_lod  t_lod_fstar t_alt_count t_ref_count judgement
# INDEL: build  contig  start_position  end_position    ref_allele  alt_allele1 alt_allele2 tumor_name  normal_name tumor_f n_ref_count n_alt_count t_ref_count t_alt_count init_n_lod  init_t_lod  judgement

df[:end]=df[:POS]
# multiple alt alleles
df[:multiALT]=map(x->1*(length(split(x,":"))>1),df[:ALT])
k=findall(df[:multiALT].>0)
df[:ALT]=map(x->split(x,":")[1],df[:ALT])
print(join(names(df),"\n"))
df[:POPAF]=map(x->split(x,":")[1],df[:POPAF])
df[:NLOD]=map(x->split(x,":")[1],df[:NLOD])
df[:TLOD]=map(x->split(x,":")[1],df[:TLOD])

# MNPs
df[:MNP]=1 .*[(length(r)>1)&(length(r)==length(a))  for (r,a) = zip(df[:REF], df[:ALT])]
LA=map(x->length(x),df[:ALT])
LR=map(x->length(x),df[:REF])
k=findall(df[:MNP].>0)
df[k,:end]=df[k,:POS]+LA[k].-1

# start for indels is POS+1
df[:DEL]=1 .*[(length(r)>1)&(length(r)>length(a))  for (r,a) = zip(df[:REF], df[:ALT])]
p=df[:POS]
p=p+df[:DEL]
df[:POS]=p
k = findall(df[:DEL].>0)
ref=df[k,:REF]
ref=map(x->x[2:end],ref)
alt=df[k,:ALT]
alt=map(x->x[2:end],alt)
k1 = findall(map(x-> length(x)<1, alt))
alt[k1]=map(x-> "-", alt[k1])
df[k,:REF]=ref
df[k,:ALT]=alt
dref = map(x-> length(x), ref)
dalt = map(x-> length(x), alt)
ddel = dref
df[k,:end]=df[k,:POS]+ddel

# calc end for INS
df[:INS]=1 .*[(length(a)>1)&(length(a)>length(r))  for (r,a) = zip(df[:REF], df[:ALT])]
k = findall(df[:INS].>0)
ref=df[k,:REF]
ref=map(x->x[2:end],ref)
k1 = findall(map(x-> length(x)<1, ref))
ref[k1]=map(x-> "-", ref[k1])
alt=df[k,:ALT]
alt=map(x->x[2:end],alt)
df[k,:REF]=ref
df[k,:ALT]=alt
df[k,:end]=df[k,:POS].+1

# add build column
df[:build]=fill("37",size(df,1))
df[:tumor_barcode]=fill(tumor_id,size(df,1))
df[:normal_barcode]=fill(normal_id,size(df,1))
df[:judgement]=fill("KEEP",size(df,1))

# M1 counts
df[:n_alt_count]=df[:NORMAL_AD_ALT]
df[:n_ref_count]=df[:NORMAL_AD_REF]
df[:t_lod_fstar]=df[:TLOD]
df[:t_alt_count]=df[:TUMOR_AD_ALT]
df[:t_ref_count]=df[:TUMOR_AD_REF]
#rename!(df, [:CHRO,:POS, :REF, :ALT], [:chr, :start, :ref_allele,:alt_allele])#from = [:DB, :HCNT, :MAX_ED,:RPA,:STR, :RU, :NLOD,:TUMOR_ALT_F1R2,:MIN_ED,:ECNT, :TLOD, :TUMOR_PID,:TUMOR_AF,:TUMOR_FOXOG, :TUMOR_QSS, :TUMOR_ALT_F2R1,:TUMOR_AD_REF,:TUMOR_AD_ALT, :TUMOR_REF_F2R1, :TUMOR_REF_F1R2,:NORMAL_AF,:NORMAL_AD_REF,:NORMAL_AD_ALT]
rename!(df, [f => t  for (f,t) = zip([:CHRO,:POS, :REF, :ALT], [:chr, :start, :ref_allele,:alt_allele])])

maf=df

permutecols!(maf, [:chr,:start,:end,:ref_allele,:alt_allele,:tumor_barcode,:normal_barcode,:POPAF,:NLOD,:PON,:TLOD,:build,:judgement,:n_alt_count,:n_ref_count,:t_lod_fstar,:t_alt_count,:t_ref_count,:NORMAL_AD_REF,:NORMAL_AD_ALT,:TUMOR_AD_REF,:TUMOR_AD_ALT,:FILTER,:multiALT, :MNP, :DEL, :INS,:VCF_REF,:VCF_ALT,:VCF_POS])

for c in names(maf)
    #println(c)
    if ~isString(maf[1,c])
        maf[c] = map(x -> string(x),maf[c])
    end
end

for c in names(maf)
    println(c)
    k1=map(x->ismissing(x)*1,maf[c])
    k=findall(k1.>0)
    k2=findall(k1.==0)
    if mean(k1)==1.0
        deletecols!(maf, c)
        continue
    end
    k=findall(map(x->x=="NA",maf[c]))
    maf[k,c] = ""
    k=findall(map(x->uppercase(x)=="TRUE",maf[c]))
    maf[k,c] = "1"
    k=findall(map(x->uppercase(x)=="FALSE",maf[c]))
    maf[k,c] = "0"
end

maflite_all=pair_id*".m2.all.maflite.tsv"
open(maflite_all, "w") do f
    writedlm(f, reshape(names(maf), 1, length(names(maf))), '\t')
    #writedlm(f, convert(Array,maf), '\t')
    writedlm(f, convert(Matrix,maf), '\t')
end

kpass = findall(map(x-> uppercase(x)=="PASS", maf[:FILTER]))
maf=maf[kpass,:]
deletecols!(maf,[:FILTER])

maflite_pass=pair_id*".m2.pass.maflite.tsv"
open(maflite_pass, "w") do f
    writedlm(f, reshape(names(maf), 1, length(names(maf))), '\t')
    writedlm(f, convert(Matrix,maf), '\t')
end


text:
# Autogenerated wrapper script for EarCut_jll for i686-w64-mingw32
export libearcut

JLLWrappers.@generate_wrapper_header("EarCut")
JLLWrappers.@declare_library_product(libearcut, "libearcut.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libearcut,
        "bin\\libearcut.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using JuLIP, Test
using JuLIP.Testing
using JuLIP.FIO

# comparison of Atoms.data dictionaries
import Base
function Base.isapprox(d1::Dict{Any, JuLIP.JData{T}}, d2::Dict{Any, JuLIP.JData{T}}; tol = sqrt(eps(T))) where T
    for (k1, v1) in d1
       k1 ∈ keys(d2) || (@error "key $k1 not in d2"; return false)
       if v1.data isa AbstractArray || v1.data isa AbstractFloat
           isapprox(v1.data, d2[k1].data; atol=tol)  || (@error "key $k1 value mismatch $(v1.data) !≈ $(d2[v2].data)"; return false)
       else
           v1.data == d2[k1].data || (@error "key $k1 value mismatch $(v1.data) != $(d2[v2].data)"; return false)
       end
   end
   return true
end

h3("Testing single `Atoms` <-> `Dict`")
at = bulk(:Cu, cubic=true) * 3
set_pbc!(at, (true, false, true))
rattle!(at, 0.1)
D = write_dict(at)
at1 = Atoms(D)
println(@test(at == at1))
at2 = read_dict(D)
println(@test(at == at2))

h3("Test JSON fio")
fn = tempname()
save_dict(fn, D)
D1 = load_dict(fn)
# D1 == D  => this will be false so don't test it!
at3 = read_dict(D1)
println(@test at3 == at1)

h3("Test array of Atoms <-> Dict")
ats = [ (bulk(:Cu) * rand(2:4)) for n = 1:5 ]
Ds = Dict("ats" => write_dict.(ats))
ats1 = read_dict.(Ds["ats"])
println(@test ats1 == ats)

h3("Test JSON fio for array")
fn = tempname()
save_dict(fn, Ds)
Ds1 = load_dict(fn)
ats2 = read_dict.(Ds1["ats"])
println(@test ats == ats2)

h3("Test ExtXYZ fio for Atoms")
@testset "extxyz" begin
    filename = tempname() * ".xyz"
    seq0 = [bulk(:Si) * 3 for i=1:10]
    for atoms in seq0
        rattle!(atoms, 0.1)
        set_calculator!(atoms, StillingerWeber())
        set_data!(atoms, "energy", energy(atoms))
        set_data!(atoms, "stress", stress(atoms))
        set_data!(atoms, "forces", forces(atoms))
    end
    write_extxyz(filename, seq0)
    seq1 = read_extxyz(filename)
    @test all(seq0 .≈ seq1)

    data0 = [atoms.data for atoms in seq0]
    data1 = [atoms.data for atoms in seq1]
    @test all(isapprox.(data1, data0; tol=1e-6))

    seq2 = read_extxyz(filename, 4:10)
    frame = read_extxyz(filename, 4)
    @test all(seq1[4:10] .≈ seq2)

    data1 = [atoms.data for atoms in seq1[4:10]]
    data2 = [atoms.data for atoms in seq2]
    @test all(isapprox.(data1, data2; tol=1e-6))

    f = open(filename, "r")
    seq3 = read_extxyz(f)
    close(f)

    @test all(seq1 .≈ seq3)
    data1 = [atoms.data for atoms in seq1]
    data3 = [atoms.data for atoms in seq3]
    @test all(isapprox.(data1, data3; tol=1e-6))

    at5 = read_extxyz(filename, 5)
    @test at5[1] ≈ seq1[5]
    data5 = [atoms.data for atoms in at5]
    data1 = [atoms.data for atoms in seq1[5:5]]
    @test all(isapprox.(data1, data5; tol=1e-6))
end


text:
module Repo # Octo

using ..DBMS
using ..Backends
using ..AdapterBase
using ..Queryable: Structured, SubQuery, FromItem
using ..Octo: Raw
using ..Schema # Schema.validates
using ..Pretty: show

struct NeedsConnectError <: Exception
    msg
end

struct ExecuteResult
end

@enum RepoLogLevel::Int32 begin
    LogLevelDebugSQL = -1
    LogLevelInfo = 0
end

const current = Dict{Symbol, Union{Nothing, Module, RepoLogLevel}}(
    :adapter => nothing,
    :loader => nothing,
    :log_level => LogLevelInfo
)

const color_params = :yellow

current_loader() = current[:loader]
function current_adapter() # throw Repo.NeedsConnectError
    if current[:adapter] isa Nothing
        throw(NeedsConnectError("Needs to Repo.connect"))
    else
        current[:adapter]
    end
end

function set_log_level(level::RepoLogLevel)
    current[:log_level] = level
end

"""
    Repo.debug_sql(debug::Bool = true)
"""
function debug_sql(debug::Bool = true)
    current[:log_level] = debug ? LogLevelDebugSQL : LogLevelInfo
end

function print_debug_sql_params(io, params)
    printstyled(io, "   ")
    for (idx, x) in enumerate(params)
        printstyled(io, repr(x), color=color_params)
        length(params) != idx && printstyled(io, ", ")
    end
end

function print_debug_sql(stmt::Structured, params = nothing)
    if current[:log_level] <= LogLevelDebugSQL
        buf = IOBuffer()
        show(IOContext(buf, :color=>true), MIME"text/plain"(), stmt)
        !(params isa Nothing) && print_debug_sql_params(IOContext(buf, :color=>true), params)
        @info String(take!(buf))
    end
end

# Repo.connect
"""
    Repo.connect(; adapter::Module, database::Union{Nothing,Type{D} where {D <: DBMS.AbstractDatabase}}=nothing, kwargs...)
"""
function connect(; adapter::Module, database::Union{Nothing,Type{D} where {D <: DBMS.AbstractDatabase}}=nothing, kwargs...)
    if database !== nothing
        adapter.Database[:ID] = database
    end
    AdapterBase.current[:database] = adapter.DatabaseID()

    current[:adapter] = adapter
    loader = Backends.backend(adapter)
    current[:loader] = loader

    Base.invokelatest(loader.db_connect; kwargs...)
end

# Repo.disconnect
"""
    Repo.disconnect()
"""
function disconnect()
    loader = current_loader()
    disconnected = loader.db_disconnect()
    current[:adapter] = nothing
    current[:loader] = nothing
    disconnected
end

function _get_primary_key(M) # throw Schema.PrimaryKeyError
    Tname = Base.typename(M)
    info = Schema.tables[Tname]
    if haskey(info, :primary_key)
        a = current_adapter()
        table = a.from(M)
        primary_key = Symbol(info[:primary_key])
        a.Field(table, primary_key)
    else
        throw(Schema.PrimaryKeyError(""))
    end
end

function _get_primary_key_with(M, nt::NamedTuple) # throw Schema.PrimaryKeyError
    Tname = Base.typename(M)
    info = Schema.tables[Tname]
    if haskey(info, :primary_key)
        key = _get_primary_key(M)
        primary_key = Symbol(info[:primary_key])
        pk = getfield(nt, primary_key)
        (key, pk)
     else
        throw(Schema.PrimaryKeyError(""))
    end
end

# Repo.query

"""
    Repo.query(stmt::Structured)
"""
function query(stmt::Structured)
    a = current_adapter()
    sql = a.to_sql(stmt)
    print_debug_sql(stmt)
    loader = current_loader()
    loader.query(sql)
end

"""
    Repo.query(M::Type)
"""
function query(M::Type)
    a = current_adapter()
    table = a.from(M)
    query([a.SELECT * a.FROM table])
end

"""
    Repo.query(from::FromItem)
"""
function query(from::FromItem)
    query(from.__octo_model)
end

"""
    Repo.query(subquery::SubQuery)
"""
function query(subquery::SubQuery)
    query(subquery.__octo_query)
end

"""
    Repo.query(rawquery::Octo.Raw)
"""
function query(rawquery::Raw)
    query([rawquery])
end

### Repo.query - vals::Vector
"""
    Repo.query(stmt::Structured, vals::Vector)
"""
function query(stmt::Structured, vals::Vector) # throw Backends.UnsupportedError
    a = current_adapter()
    prepared = a.to_sql(stmt)
    print_debug_sql(stmt, vals)
    loader = current_loader()
    loader.query(prepared, vals)
end

### Repo.query - pk
"""
    Repo.query(M::Type, pk::Union{Int, String})
"""
function query(M::Type, pk::Union{Int, String}) # throw Schema.PrimaryKeyError
    a = current_adapter()
    table = a.from(M)
    key = _get_primary_key(M)
    query([a.SELECT * a.FROM table a.WHERE key == pk])
end

"""
    Repo.query(from::FromItem, pk::Union{Int, String})
"""
function query(from::FromItem, pk::Union{Int, String}) # throw Schema.PrimaryKeyError
    query(from.__octo_model, pk)
end

### Repo.query - pk_range
"""
    Repo.query(M::Type, pk_range::UnitRange{Int64})
"""
function query(M::Type, pk_range::UnitRange{Int64}) # throw Schema.PrimaryKeyError
    a = current_adapter()
    table = a.from(M)
    key = _get_primary_key(M)
    query([a.SELECT * a.FROM table a.WHERE key a.BETWEEN pk_range.start a.AND pk_range.stop])
end

"""
    Repo.query(from::FromItem, pk_range::UnitRange{Int64}
"""
function query(from::FromItem, pk_range::UnitRange{Int64}) # throw Schema.PrimaryKeyError
    query(from.__octo_model, pk_range)
end

### Repo.query - nt::NamedTuple
"""
    Repo.query(stmt::Structured, nt::NamedTuple)
"""
function query(stmt::Structured, nt::NamedTuple)
    a = current_adapter()
    v = vec(stmt)
    for (idx, kv) in enumerate(pairs(nt))
        key = a.Field(nothing, kv.first)
        push!(v, key == kv.second)
        idx != length(nt) && push!(v, a.AND)
    end
    query(v)
end

"""
    Repo.query(M::Type, nt::NamedTuple)
"""
function query(M::Type, nt::NamedTuple)
    a = current_adapter()
    table = a.from(M)
    v = [a.SELECT, *, a.FROM, table, a.WHERE]
    for (idx, kv) in enumerate(pairs(nt))
        key = Base.getproperty(table, kv.first)
        push!(v, key == kv.second)
        idx != length(nt) && push!(v, a.AND)
    end
    query(v)
end

"""
    Repo.query(from::FromItem, nt::NamedTuple)
"""
function query(from::FromItem, nt::NamedTuple)
    query(from.__octo_model, nt)
end

"""
    Repo.query(subquery::SubQuery, nt::NamedTuple)
"""
function query(subquery::SubQuery, nt::NamedTuple)
    query(subquery.__octo_query, nt)
end

"""
    Repo.query(rawquery::Octo.Raw, nt::NamedTuple)
"""
function query(rawquery::Raw, nt::NamedTuple)
    query([rawquery], nt)
end

# Repo.get
"""
    Repo.get(M::Type, pk::Union{Int, String})
"""
function get(M::Type, pk::Union{Int, String}) # throw Schema.PrimaryKeyError
    query(M, pk)
end

"""
    Repo.get(M::Type, pk_range::UnitRange{Int64})
"""
function get(M::Type, pk_range::UnitRange{Int64}) # throw Schema.PrimaryKeyError
    query(M, pk_range)
end

"""
    Repo.get(M::Type, nt::NamedTuple)
"""
function get(M::Type, nt::NamedTuple)
    query(M, nt)
end

# Repo.execute
"""
    Repo.execute(stmt::Structured)
"""
function execute(stmt::Structured)
    a = current_adapter()
    sql = a.to_sql(stmt)
    print_debug_sql(stmt)
    loader = current_loader()
    loader.execute(sql)
    nothing
end

"""
    Repo.execute(stmt::Structured, vals::Vector)
"""
function execute(stmt::Structured, vals::Vector)
    a = current_adapter()
    prepared = a.to_sql(stmt)
    print_debug_sql(stmt, vals)
    loader = current_loader()
    loader.execute(prepared, vals)
    nothing
end

"""
    Repo.execute(stmt::Structured, nts::Vector{<:NamedTuple})
"""
function execute(stmt::Structured, nts::Vector{<:NamedTuple})
    a = current_adapter()
    prepared = a.to_sql(stmt)
    print_debug_sql(stmt, nts)
    loader = current_loader()
    loader.execute(prepared, nts)
    nothing
end

execute(raw::AdapterBase.Raw)                            = execute([raw])
execute(raw::AdapterBase.Raw, nt::NamedTuple)            = execute([raw], [nt])
execute(raw::AdapterBase.Raw, nts::Vector{<:NamedTuple}) = execute([raw], nts)
execute(raw::AdapterBase.Raw, vals::Vector)              = execute([raw], vals)

# Repo.insert!
"""
    Repo.insert!(M::Type, nts::Vector{<:NamedTuple})
"""
function insert!(M::Type, nts::Vector{<:NamedTuple})
    if !isempty(nts)
        Schema.validates.(M, nts) # throw InvalidChangesetError
        a = current_adapter()
        table = a.from(M)
        nt = first(nts)
        fieldnames = a.Enclosed(collect(keys(nt)))
        values = a.placeholders(length(nt))
        execute([a.INSERT a.INTO table fieldnames a.VALUES values], nts)
   end
   nothing
end

"""
    Repo.insert!(M, nt::NamedTuple)
"""
function insert!(M, nt::NamedTuple)
    insert!(M, [nt])
end

"""
    Repo.insert!(M::Type, vals::Vector{<:Tuple})
"""
function insert!(M::Type, vals::Vector{<:Tuple})
    if !isempty(vals)
        a = current_adapter()
        table = a.from(M)
        execute([a.INSERT a.INTO table a.VALUES a.VectorOfTuples(vals)])
   end
   nothing
end

# Repo.update!
"""
    Repo.update!(M::Type, nt::NamedTuple)
"""
function update!(M::Type, nt::NamedTuple) # throw Schema.PrimaryKeyError
    Schema.validates(M, nt) # throw InvalidChangesetError
    a = current_adapter()
    (key, pk) = _get_primary_key_with(M, nt)
    table = a.from(M)
    rest = filter(kv -> kv.first != key.name, pairs(nt))
    v = Any[a.UPDATE, table, a.SET]

    # tup[1] : idx
    # tup[2] : kv
    push!(v, tuple(map(tup -> a.Field(table, tup[2].first) == a.placeholder(tup[1]), enumerate(rest))...))

    push!(v, a.WHERE)
    push!(v, key == pk)
    execute(v, collect(values(rest)))
    nothing
end

# Repo.delete!
"""
    Repo.delete!(M::Type, nt::NamedTuple)
"""
function delete!(M::Type, nt::NamedTuple) # throw Schema.PrimaryKeyError
    a = current_adapter()
    (key, pk) = _get_primary_key_with(M, nt)
    table = a.from(M)
    execute([a.DELETE a.FROM table a.WHERE key == pk])
    nothing
end

"""
    Repo.delete!(M::Type, pk_range::UnitRange{Int64})
"""
function delete!(M::Type, pk_range::UnitRange{Int64}) # throw Schema.PrimaryKeyError
    a = current_adapter()
    table = a.from(M)
    key = _get_primary_key(M)
    execute([a.DELETE a.FROM table a.WHERE key a.BETWEEN pk_range.start a.AND pk_range.stop])
    nothing
end

end # module Octo.Repo


text:
# This file is part of the TaylorSeries.jl Julia package, MIT license
#
# Luis Benet & David P. Sanders
# UNAM
#
# MIT Expat license
#

## Differentiating ##
"""
    derivative(a)

Return the `Taylor1` polynomial of the differential of `a::Taylor1`.
The last coefficient is set to zero.

The function `differentiate` is an exact synonym of `derivative`.
"""
function derivative(a::Taylor1)
    res = zero(a)
    @inbounds for ord in eachindex(a)
        derivative!(res, a, ord)
        # ord == a.order-1 && break
    end
    return res
end

"""
    differentiate

An exact synonym of [`derivative`](@ref).
"""
const differentiate = derivative

"""
    derivative!(res, a) --> nothing

In-place version of `derivative`. Compute the `Taylor1` polynomial of the
differential of `a::Taylor1` and save it into `res`. The last coefficient is
set to zero.
"""
function derivative!(res::Taylor1, a::Taylor1)
    @inbounds for ord in eachindex(a)
        derivative!(res, a, ord)
        # ord == a.order-1 && break
    end
    res[a.order] = zero(a[0])
    nothing
end

"""
    derivative!(p, a, k) --> nothing

Update in-place the `k-th` expansion coefficient `p[k]` of `p = derivative(a)`
for both `p` and `a` `Taylor1`.

The coefficients are given by

```math
p_k = (k+1)a_{k+1}.
```

"""
derivative!(p::Taylor1, a::Taylor1, k::Int) = k < a.order ? p[k] = (k+1)*a[k+1] : nothing

"""
    derivative(a, n)

Compute recursively the `Taylor1` polynomial of the n-th derivative of
`a::Taylor1`.
"""
function derivative(a::Taylor1{T}, n::Int) where {T<:Number}
    @assert a.order ≥ n ≥ 0
    if n==0
        return a
    else
        res = deepcopy(a)
        for i in 1:n
            derivative!(res, res)
        end
        return res
    end
end

"""
    derivative(n, a)

Return the value of the `n`-th derivative of the polynomial `a`.
"""
function derivative(n::Int, a::Taylor1{T}) where {T<:Number}
    @assert a.order ≥ n ≥ 0
    factorial( widen(n) ) * a[n] :: T
end

## Integrating ##
"""
    integrate(a, [x])

Return the integral of `a::Taylor1`. The constant of integration
(0-th order coefficient) is set to `x`, which is zero if ommitted.
"""
function integrate(a::Taylor1{T}, x::S) where {T<:Number,S<:Number}
    order = get_order(a)
    aa = a[0]/1 + zero(x)
    R = typeof(aa)
    coeffs = Array{typeof(aa)}(undef, order+1)
    fill!(coeffs, zero(aa))
    @inbounds for i = 1:order
        coeffs[i+1] = a[i-1] / i
    end
    @inbounds coeffs[1] = convert(R, x)
    return Taylor1(coeffs, a.order)
end
integrate(a::Taylor1{T}) where {T<:Number} = integrate(a, zero(a[0]))

## Differentiation ##
"""
    derivative(a, r)

Partial differentiation of `a::HomogeneousPolynomial` series with respect
to the `r`-th variable.
"""
function derivative(a::HomogeneousPolynomial, r::Int)
    @assert 1 ≤ r ≤ get_numvars()
    # @show zero(a[1]) zero(eltype(a))
    T = eltype(a)
    a.order == 0 && return HomogeneousPolynomial(zero(a[1]), 0)
    @inbounds num_coeffs = size_table[a.order]
    coeffs = zeros(T,num_coeffs)
    @inbounds posTb = pos_table[a.order]
    @inbounds num_coeffs = size_table[a.order+1]

    @inbounds for i = 1:num_coeffs
        iind = coeff_table[a.order+1][i]
        n = iind[r]
        n == 0 && continue
        iind[r] -= 1
        kdic = in_base(get_order(),iind)
        pos = posTb[kdic]
        coeffs[pos] = n * a[i]
        iind[r] += 1
    end

    return HomogeneousPolynomial{T}(coeffs, a.order-1)
end
derivative(a::HomogeneousPolynomial, s::Symbol) = derivative(a, lookupvar(s))

"""
    derivative(a, r)

Partial differentiation of `a::TaylorN` series with respect
to the `r`-th variable. The `r`-th variable may be also
specified through its symbol.
"""
function derivative(a::TaylorN, r=1::Int)
    T = eltype(a)
    coeffs = Array{HomogeneousPolynomial{T}}(undef, a.order)

    @inbounds for ord = 1:a.order
        coeffs[ord] = derivative( a[ord], r)
    end
    return TaylorN{T}( coeffs, a.order )
end
derivative(a::TaylorN, s::Symbol) = derivative(a, lookupvar(s))

"""
    derivative(a::TaylorN{T}, ntup::NTuple{N,Int})

Return a `TaylorN` with the partial derivative of `a` defined
by `ntup::NTuple{N,Int}`, where the first entry is the number
of derivatives with respect to the first variable, the second is
the number of derivatives with respect to the second, and so on.
"""
function derivative(a::TaylorN, ntup::NTuple{N,Int}) where {N}

    @assert N == get_numvars() && all(ntup .>= 0)

    sum(ntup) > a.order && return zero(a)
    sum(ntup) == 0 && return copy(a)

    aa = copy(a)
    for nvar in 1:get_numvars()
        for numder in 1:ntup[nvar]
            aa = derivative(aa, nvar)
        end
    end

    return aa
end

"""
    derivative(ntup::NTuple{N,Int}, a::TaylorN{T})

Returns the value of the coefficient of `a` specified by
`ntup::NTuple{N,Int}`, multiplied by the corresponding
factorials.
"""
function derivative(ntup::NTuple{N,Int}, a::TaylorN) where {N}

    @assert N == get_numvars() && all(ntup .>= 0)

    c = getcoeff(a, [ntup...])
    for ind = 1:get_numvars()
        c *= factorial(ntup[ind])
    end

    return c
end

## Gradient, jacobian and hessian
"""
```
    gradient(f)
    ∇(f)
```

Compute the gradient of the polynomial `f::TaylorN`.
"""
function gradient(f::TaylorN)
    T = eltype(f)
    numVars = get_numvars()
    grad = Array{TaylorN{T}}(undef, numVars)
    @inbounds for nv = 1:numVars
        grad[nv] = derivative(f, nv)
    end
    return grad
end
const ∇ = TaylorSeries.gradient

"""
```
    jacobian(vf)
    jacobian(vf, [vals])
```

Compute the jacobian matrix of `vf`, a vector of `TaylorN` polynomials,
evaluated at the vector `vals`. If `vals` is ommited, it is evaluated at zero.
"""
function jacobian(vf::Array{TaylorN{T},1}) where {T<:Number}
    numVars = get_numvars()
    @assert length(vf) == numVars
    jac = Array{T}(undef, numVars, numVars)

    @inbounds for comp = 1:numVars
        jac[:,comp] = vf[comp][1][1:end]
    end

    return transpose(jac)
end
function jacobian(vf::Array{TaylorN{T},1}, vals::Array{S,1}) where {T<:Number,S<:Number}
    R = promote_type(T,S)
    numVars = get_numvars()
    @assert length(vf) == numVars == length(vals)
    jac = Array{R}(undef, numVars, numVars)

    for comp = 1:numVars
        @inbounds grad = gradient( vf[comp] )
        @inbounds for nv = 1:numVars
            jac[nv,comp] = evaluate(grad[nv], vals)
        end
    end

    return transpose(jac)
end

function jacobian(vf::Array{Taylor1{TaylorN{T}},1}) where {T<:Number}
    vv = convert(Array{TaylorN{Taylor1{T}},1}, vf)
    jacobian(vv)
end

"""
```
    jacobian!(jac, vf)
    jacobian!(jac, vf, [vals])
```

Compute the jacobian matrix of `vf`, a vector of `TaylorN` polynomials
evaluated at the vector `vals`, and write results to `jac`. If `vals` is ommited,
it is evaluated at zero.
"""
function jacobian!(jac::Array{T,2}, vf::Array{TaylorN{T},1}) where {T<:Number}
    numVars = get_numvars()
    @assert length(vf) == numVars
    @assert (numVars, numVars) == size(jac)
    for comp2 = 1:numVars
        for comp1 = 1:numVars
            @inbounds jac[comp1,comp2] = vf[comp1][1][comp2]
        end
    end
    nothing
end
function jacobian!(jac::Array{T,2}, vf::Array{TaylorN{T},1},
        vals::Array{T,1}) where {T<:Number}

    numVars = get_numvars()
    @assert length(vf) == numVars == length(vals)
    @assert (numVars, numVars) == size(jac)
    for comp = 1:numVars
        @inbounds for nv = 1:numVars
            jac[nv,comp] = evaluate(derivative(vf[nv], comp), vals)
        end
    end
    nothing
end

"""
```
    hessian(f)
    hessian(f, [vals])
```

Return the hessian matrix (jacobian of the gradient) of `f::TaylorN`,
evaluated at the vector `vals`. If `vals` is ommited, it is evaluated at
zero.
"""
hessian(f::TaylorN{T}, vals::Array{S,1}) where {T<:Number,S<:Number} =
    (R = promote_type(T,S); jacobian( gradient(f), vals::Array{R,1}) )

hessian(f::TaylorN{T}) where {T<:Number} = hessian( f, zeros(T, get_numvars()) )

"""
```
    hessian!(hes, f)
    hessian!(hes, f, [vals])
```

Return the hessian matrix (jacobian of the gradient) of `f::TaylorN`,
evaluated at the vector `vals`, and write results to `hes`. If `vals` is
ommited, it is evaluated at zero.
"""
hessian!(hes::Array{T,2}, f::TaylorN{T}, vals::Array{T,1}) where {T<:Number} =
    jacobian!(hes, gradient(f), vals)

hessian!(hes::Array{T,2}, f::TaylorN{T}) where {T<:Number} =
    jacobian!(hes, gradient(f))

##Integration
"""
    integrate(a, r)

Integrate the `a::HomogeneousPolynomial` with respect to the `r`-th
variable. The returned `HomogeneousPolynomial` has no added constant of
integration. If the order of a corresponds to `get_order()`, a zero
`HomogeneousPolynomial` of 0-th order is returned.

"""
function integrate(a::HomogeneousPolynomial, r::Int)
    @assert 1 ≤ r ≤ get_numvars()

    order_max = get_order()
    a.order == order_max && return HomogeneousPolynomial(zero(a[1]/1), 0)

    @inbounds posTb = pos_table[a.order+2]
    @inbounds num_coeffs = size_table[a.order+1]

    T = promote_type(eltype(a), eltype(a[1]/1))
    coeffs = zeros(T, size_table[a.order+2])

    @inbounds for i = 1:num_coeffs
        iind = coeff_table[a.order+1][i]
        n = iind[r]
        n == order_max && continue
        iind[r] += 1
        kdic = in_base(get_order(), iind)
        pos = posTb[kdic]
        coeffs[pos] = a[i] / (n+1)
        iind[r] -= 1
    end

    return HomogeneousPolynomial(coeffs, a.order+1)
end
integrate(a::HomogeneousPolynomial, s::Symbol) = integrate(a, lookupvar(s))

"""
    integrate(a, r, [x0])

Integrate the `a::TaylorN` series with respect to the `r`-th variable,
where `x0` the integration constant and must be independent
of the `r`-th variable; if `x0` is ommitted, it is taken as zero.
"""
function integrate(a::TaylorN, r::Int)
    T = promote_type(eltype(a), eltype(a[0]/1))
    order_max = min(get_order(), a.order+1)
    coeffs = zeros(HomogeneousPolynomial{T}, order_max)

    @inbounds for ord = 0:order_max-1
        coeffs[ord+1] = integrate( a[ord], r)
    end

    return TaylorN(coeffs)
end
function integrate(a::TaylorN, r::Int, x0::TaylorN)
    # Check constant of integration is independent of re
    @assert derivative(x0, r) == 0.0 """
    The integration constant ($x0) must be independent of the
    $(_params_TaylorN_.variable_names[r]) variable"""

    res = integrate(a, r)
    return x0+res
end
integrate(a::TaylorN, r::Int, x0) =
    integrate(a,r,TaylorN(HomogeneousPolynomial([convert(eltype(a),x0)], 0)))

integrate(a::TaylorN, s::Symbol) = integrate(a, lookupvar(s))
integrate(a::TaylorN, s::Symbol, x0::TaylorN) = integrate(a, lookupvar(s), x0)
integrate(a::TaylorN, s::Symbol, x0) = integrate(a, lookupvar(s), x0)


text:
import Entropy
import Base.Test

#Testing nsb_entropy

function init()
	srand(1234)
	x = rand(1:100,50)
	n = Entropy.counts(x)
	S_nsb = Entropy.NSBEntropy.NSBEntropy(n,100)
	return S_nsb
end

function test_B_xiK()
	S_nsb = init()
	B =  Entropy..B_xiK(1.2, S_nsb)
	return Base.Test.@test_approx_eq B 1.3853866920075761
end

function test_find_nsb_entropy()
	S_nsb = init()
	Entropy.find_nsb_entropy(S_nsb,1e-5)
	return Base.Test.@test_approx_eq S_nsb.S_nsb 4.484256684719789
end

function test_mlog_evidence()
	S_nsb = init()
	nsb_mlog = Entropy.mlog_evidence(200*S_nsb.K,S_nsb)
	return Base.Test.@test_approx_eq nsb_mlog 46.05806318661056
end

S_nsb = init()
ms = Entropy.meanS2(370.0850,[26,9,2], [1,2,3], 100)
Base.Test.@test_approx_eq ms 20.0490
B =  Entropy.B_xiK(1.2, S_nsb)
Base.Test.@test_approx_eq B 1.3853866920075761

nsb_mlog = Entropy.mlog_evidence(200*S_nsb.K,S_nsb)
Base.Test.@test_approx_eq nsb_mlog 59.85364833069893

Entropy.find_nsb_entropy(S_nsb,1e-5)
Base.Test.@test_approx_eq S_nsb.S_nsb 4.3952671710926685


text:
const CACHE_DIR = "../_cache/"

function getcache(dur)
  days = numofdays(dur)
  res = Result[]

  for day in days
      p = daycachefilename(day)
      if ispath(p) == false
          continue
      end
      println("reading " * p)
      raw = read(p, String)
      c = Result(raw)
      push!(res, c)
  end

  return res
end

function savecache(res)
  if ispath(CACHE_DIR) == false
      mkdir(CACHE_DIR)
  end
  p = rescachefilename(res)
  if ispath(p) == false
      savetofile(res.raw, p)
  end
end

function daycachefilename(day)
  return CACHE_DIR * Dates.format(day, "yyyy_mm_dd") * ".csv"
end

function rescachefilename(res)
  return CACHE_DIR * Dates.format(res.header.start_time, "yyyy_mm_dd") * ".csv"
end

function findday(results, day)
  for res in results
      d = Date(Dates.year(res.header.start_time),
      Dates.month(res.header.start_time),
      Dates.day(res.header.start_time))
      if d == day
          return res
      end
  end

  return false
end


text:
"""
    nout(v)

Return the number of output neurons of vertex `v`.

This is typically the number of rows/columns of a parameter `Matrix` for fully connected layers or the number of output channels
in a convolutional layer.

Note that NaiveNASlib does not have the capability to figure this out by itself so this must be provided by the function wrapped
in the vertex. For [`SizeTransparent`](@ref) vertices, NaiveNASlib will default to computing `nout(v)` from the sizes of the inputs.

There are three ways to implement `nin` for a function/callable `f` of type `F`:

    1. nout(f::F)
    2. nout(st::ST, f) and st = shapetrait(f)
    3. nout(f::F, t::MutationTrait, v::AbstractVertex)

While `1` is the most straight forward, `2` can be useful of there are many different `f`s which happen to share a common method
for determining the size.
Option `3` is when the implementation might want to use other information from `v` or `t` and is left as an escape hatch.
"""
function nout end

"""
    nin(v)

Return the number of input neurons of vertex `v`.

This is typically the number of rows/columns of a parameter `Matrix` for fully connected layers or the number of input channels
in a convolutional layer.

Note that NaiveNASlib does not have the capability to figure this out by itself so this must be provided by the function wrapped
in the vertex. For [`SizeTransparent`](@ref) vertices, NaiveNASlib will default to computing `nin(v)` from the sizes of the inputs.

There are three ways to implement `nin` for a function/callable `f` of type `F`:

    1. nin(f::F)
    2. nin(st::ST, f) and st = shapetrait(f)
    3. nin(f::F, t::MutationTrait, v::AbstractVertex)

While `1` is the most straight forward, `2` can be useful of there are many different `f`s which happen to share a common method
for determining the size.
Option `3` is when the implementation might want to use other information from `v` or `t` and is left as an escape hatch.
"""
function nin end

"""
    Δsize!([utilityfun], [s::AbstractΔSizeStrategy], vs::AbstractVector{<:AbstractVertex})

Change size of (potentially) all vertices in `vs` according to the provided [`AbstractΔSizeStrategy`](@ref) (default
[`UpperInsertBound`](@ref)).

Argument `utilityfun` provides a vector `utility = utilityfun(vx)` for any vertex `vx` in the same graph as `v` where
`utility[i] > utility[j]` indicates that output neuron index `i` shall be preferred over `j` for vertex `vx`. It may also
provide a scalar which will be used as utility of all neurons of `vx`. If not provided, [`defaultutility(vx)`](@ref)
will be used.
"""
Δsize!(vs::AbstractVector{<:AbstractVertex}) = Δsize!(defaultutility, vs)
Δsize!(utilityfun, vs::AbstractVector{<:AbstractVertex}) = Δsize!(utilityfun, UpperInsertBound(0), vs)
Δsize!(s::AbstractΔSizeStrategy, vs::AbstractVector{<:AbstractVertex}) = Δsize!(Δsizetype(vs), s, vs)
Δsize!(utilityfun, s::AbstractΔSizeStrategy, vs::AbstractVector{<:AbstractVertex}) = Δsize!(utilityfun, Δsizetype(vs), s, vs)

"""
    Δsize!(s::AbstractΔSizeStrategy)
    Δsize!(utilityfun, s::AbstractΔSizeStrategy)

Change size of (potentially) all vertices which `s` has a chance to impact the size of.

Argument `utilityfun` provides a vector `utility = utilityfun(vx)` for any vertex `vx` in the same graph as `v` where
`utility[i] > utility[j]` indicates that output neuron index `i` shall be preferred over `j` for vertex `vx`. It may
also provide a scalar which will be used as utility of all neurons of `vx`. If not provided, [`defaultutility(vx)`](@ref)
will be used.
"""
Δsize!(s::AbstractΔSizeStrategy) = Δsize!(s, add_participants!(s))
Δsize!(utilityfun, s::AbstractΔSizeStrategy) = Δsize!(utilityfun, s, add_participants!(s))

"""
    Δsize!([utilityfun], [s::AbstractΔSizeStrategy], g::CompGraph)
    Δsize!([utilityfun], [s::AbstractΔSizeStrategy], v::AbstractVertex)

Change size of (potentially) all vertices of graph `g` (or graph to which `v` is connected) according to the provided
[`AbstractΔSizeStrategy`](@ref) (default [`UpperInsertBound`](@ref)).

Return true of operation was successful, false otherwise.

Argument `utilityfun` provides a vector `utility = utilityfun(vx)` for any vertex `vx` in the same graph as `v` where
`utility[i] > utility[j]` indicates that output neuron index `i` shall be preferred over `j` for vertex `vx`. It may
also provide a scalar which will be used as utility of all neurons of `vx`. If not provided, [`defaultutility(vx)`](@ref)
will be used.
"""
Δsize!(g::CompGraph) = Δsize!(defaultutility, g)
Δsize!(utilityfun, g::CompGraph) = Δsize!(utilityfun, vertices(g))
Δsize!(s::AbstractΔSizeStrategy, g::CompGraph) = Δsize!(defaultutility, s, g)
Δsize!(utilityfun, s::AbstractΔSizeStrategy, g::CompGraph) = Δsize!(utilityfun, s, vertices(g))

Δsize!(v::AbstractVertex) = Δsize!(defaultutility, v)
Δsize!(utilityfun, v::AbstractVertex) = Δsize!(utilityfun, all_in_graph(v))
Δsize!(s::AbstractΔSizeStrategy, v::AbstractVertex) =Δsize!(defaultutility, s::AbstractΔSizeStrategy, v)
Δsize!(utilityfun, s::AbstractΔSizeStrategy, v::AbstractVertex) = Δsize!(utilityfun, s, all_in_graph(v))

"""
    Δnin!([utilityfun], v, Δ...)
    Δnin!([utilityfun], v1 => Δ1, v2 => Δ2,...)
    Δnin!([utilityfun], Δs::AbstractDict)

Change input size of all provided vertices with the associated `Δ` and make the appropriate changes to other vertices
so that the graph is aligned w.r.t activations. Return `true` if successful (`false` if not successful).

For `Δ`s provided as integers it must be possible to change the size by exactly `Δ` or else the attempt will be considered failed.
A failed attempt will be retried immediately in relaxed form where the wanted size changes are moved to the objective.
The relaxation means that input size might not change by exactly `Δ`. Use [`relaxed(Δ)`](@ref) to indicate that a size change is
relaxed in the initial attempt.

For vertices with more than one input, the size change must be expressed as a tuple with one element per input.
Use `missing` to indicate that no special treatment is needed for an input. Both `missing` and [`relaxed`](@ref) can be
mixed freely inside and outside the tuples (see examples).

Note that the above constrain makes `Δnin!` much more cumbersome to use compared to [`Δnout!`](@ref) and in most cases
there are no direct advantages of using `Δnin!` over [`Δnout!`](@ref) as they both boil down to the same thing.

Argument `utilityfun` provides a vector `utility = utilityfun(vx)` for any vertex `vx` in the same graph as `v` where
`utility[i] > utility[j]` indicates that output neuron index `i` shall be preferred over `j` for vertex `vx`. It may also provide a scalar which will be used as utility of all neurons of `vx`. If not provided, [`defaultutility(vx)`](@ref) will be used.

Note that `Δnin!([utilityfun], args...)` is equivalent to  `Δsize!([utilityfun], ΔNin(args...))`.

$(generic_Δnin_docstring_examples("Δnin!"; header="```jldoctest", footer=""))
julia> Δnin!(v3, relaxed(3), missing, 2);

julia> Δnin!(v3 => (relaxed(3), missing, 2), v4 => relaxed(-2, 0)) do v
           randn(nout(v))
       end;
```
"""
function Δnin! end

"""
    Δnout!([utilityfun], v, Δ...)
    Δnout!([utilityfun], v1 => Δ1, v2 => Δ2,...)
    Δnout!([utilityfun], Δs::AbstractDict)

Change output size of all provided vertices with the associated `Δ` and make the appropriate changes to other vertices
so that the graph is aligned w.r.t activations. Return `true` if successful (`false` if not successful).

For `Δ`s provided as integers it must be possible to change the size by exactly `Δ` or else the attempt will be considered failed.
A failed attempt will be retried immediately in relaxed form where the wanted size changes are moved to the objective.
The relaxation means that output size might not change by exactly `Δ`. Use [`relaxed(Δ)`](@ref) to indicate that a size change is
relaxed in the initial attempt.

Argument `utilityfun` provides a vector `utility = utilityfun(vx)` for any vertex `vx` in the same graph as `v` where
`utility[i] > utility[j]` indicates that output neuron index `i` shall be preferred over `j` for vertex `vx`. It may also provide a scalar which will be used as utility of all neurons of `vx`. If not provided, [`defaultutility(vx)`](@ref) will be used.

Note that `Δnout!([utilityfun], args...)` is equivalent to `Δsize!([utilityfun], ΔNout(args...))`.

$(generic_Δnout_docstring_examples("Δnout!"; header="```jldoctest", footer=""))
julia> Δnout!(v1, relaxed(2));

julia> Δnout!(v1 => relaxed(2), v2 => -1) do v
           randn(nout(v))
       end;
```
"""
function Δnout! end

# Just another name for Δsize with the corresponding direction
Δnin!(args::Union{Pair{<:AbstractVertex}, AbstractDict{<:AbstractVertex}}...) = Δsize!(ΔNin(args...))
Δnout!(args::Union{Pair{<:AbstractVertex}, AbstractDict{<:AbstractVertex}}...) = Δsize!(ΔNout(args...))
Δnin!(v::AbstractVertex, Δs...) = Δnin!(v => Δs)
Δnout!(v::AbstractVertex, Δ) = Δnout!(v => Δ)

Δnin!(f, args::Union{Pair{<:AbstractVertex}, AbstractDict{<:AbstractVertex}}...) = Δsize!(f, ΔNin(args...))
Δnout!(f, args::Union{Pair{<:AbstractVertex}, AbstractDict{<:AbstractVertex}}...) = Δsize!(f, ΔNout(args...))
Δnin!(f, v::AbstractVertex, Δs...) = Δnin!(f, v => Δs)
Δnout!(f, v::AbstractVertex, Δ) = Δnout!(f, v => Δ)


text:
interval_indicators = ["VWAP", "AD", "OBV",
                       "TRANGE", "STOCH", "STOCHF", "BOP",
                       "ULTOSC", "SAR", "ADOSC"]

for func in interval_indicators
    fname = Symbol(func)
    @eval begin
        function ($fname)(symbol::String, interval::String; client = GLOBAL[], datatype::String="json", parser = "default", kwargs...)
            @argcheck in(interval, ["1min", "5min", "15min", "30min", "60min", "daily", "weekly", "monthly"])

            requiredParams = "&symbol=$symbol&interval=$interval"
            optionalParams = _parse_params(kwargs)

            uri = _form_uri_head(client, uppercase($func)) * requiredParams * optionalParams * _form_uri_tail(client, nothing, datatype)
            data = retry(_get_request, delays=Base.ExponentialBackOff(n=3, first_delay=5, max_delay=1000))(uri)
            p = _parser(parser, datatype)
            return p(data)
        end
    export $fname
    end
end

interval_seriestype_indicators = ["HT_TRENDLINE", "HT_SINE", "HT_TRENDMODE",
                                  "HT_DCPERIOD", "HT_DCPHASE", "HT_PHASOR",
                                  "MACD", "MAMA", "MACDEXT", "APO", "PPO"]

for func in interval_seriestype_indicators
    fname = Symbol(func)
    @eval begin
        function ($fname)(symbol::String, interval::String, series_type::String; client = GLOBAL[], datatype::String="json", parser = "default", kwargs...)
            @argcheck in(interval, ["1min", "5min", "15min", "30min", "60min", "daily", "weekly", "monthly"])
            @argcheck in(series_type, ["open", "high", "low", "close"])

            requiredParams = "&symbol=$symbol&interval=$interval&series_type=$series_type"
            optionalParams = _parse_params(kwargs)

            uri = _form_uri_head(client, uppercase($func)) * requiredParams * optionalParams * _form_uri_tail(client, nothing, datatype)
            data = retry(_get_request, delays=Base.ExponentialBackOff(n=3, first_delay=5, max_delay=1000))(uri)
            p = _parser(parser, datatype)
            return p(data)
        end
    export $fname
    end
end

interval_timeperiod_indicators = ["ADX", "ADXR",
                                  "AROON", "NATR", "WILLR",
                                  "CCI", "AROONOSC", "MFI",
                                  "DX", "MINUS_DI", "PLUS_DI",
                                  "MINUS_DM", "PLUS_DM", "ATR"]

for func in interval_timeperiod_indicators
    fname = Symbol(func)
    @eval begin
        function ($fname)(symbol::String, interval::String, time_period::Int64; client = GLOBAL[], datatype::String="json", parser = "default", kwargs...)
            @argcheck in(interval, ["1min", "5min", "15min", "30min", "60min", "daily", "weekly", "monthly"])
            @argcheck time_period > 0

            requiredParams = "&symbol=$symbol&interval=$interval&time_period=$time_period"
            optionalParams = _parse_params(kwargs)

            uri = _form_uri_head(client, uppercase($func)) * requiredParams * optionalParams * _form_uri_tail(client, nothing, datatype)
            data = retry(_get_request, delays=Base.ExponentialBackOff(n=3, first_delay=5, max_delay=1000))(uri)
            p = _parser(parser, datatype)
            return p(data)
        end
        export $fname
    end
end

interval_timeperiod_seriestype_indicators = ["EMA", "SMA", "WMA",
                                             "DEMA", "TEMA", "TRIMA",
                                             "KAMA", "RSI", "T3",
                                             "STOCHRSI", "MOM", "CMO",
                                             "ROC", "ROCR", "TRIX",
                                             "BBANDS", "MIDPOINT", "MIDPRICE"]

for func in interval_timeperiod_seriestype_indicators
    fname = Symbol(func)
    @eval begin
        function ($fname)(symbol::String, interval::String, time_period::Int64, series_type::String; client = GLOBAL[], datatype::String="json", parser = "default", kwargs...)
            @argcheck in(interval, ["1min", "5min", "15min", "30min", "60min", "daily", "weekly", "monthly"])
            @argcheck in(series_type, ["open", "high", "low", "close"])
            @argcheck time_period > 0

            requiredParams = "&symbol=$symbol&interval=$interval&time_period=$time_period&series_type=$series_type&datatype=$datatype&"
            optionalParams = _parse_params(kwargs)
            uri = _form_uri_head(client, uppercase($func)) * requiredParams * optionalParams * _form_uri_tail(client, nothing, datatype)
            data = retry(_get_request, delays=Base.ExponentialBackOff(n=3, first_delay=5, max_delay=1000))(uri)
            p = _parser(parser, datatype)
            return p(data)
        end
    export $fname
    end
end


text:
"""
    ExponentialLikelihood(l::AbstractLink=ExpLink())

Exponential likelihood with scale given by `l(f)`.

```math
    p(y|f) = Exponential(y | l(f))
```
"""
struct ExponentialLikelihood{Tl<:AbstractLink}
    invlink::Tl
end

ExponentialLikelihood() = ExponentialLikelihood(ExpLink())

(l::ExponentialLikelihood)(f::Real) = Exponential(l.invlink(f))

function (l::ExponentialLikelihood)(fs::AbstractVector{<:Real})
    return Product(Exponential.(l.invlink.(fs)))
end