text:
using ModelConstructors, Nullables, SMC, Test, Distributed, Dates, DataFrames, OrderedCollections, FileIO, DataStructures, LinearAlgebra, StatsBase, Random, CSV, StateSpaceRoutines, HDF5, JLD2
import ModelConstructors: @test_matrix_approx_eq, @test_matrix_approx_eq_eps
@everywhere using DSGE, JLD2, Printf, LinearAlgebra, ModelConstructors, SMC

my_tests = [
            "core",
            "parameters",
            "models/an_schorfheide/an_schorfheide",
            "models/smets_wouters/smets_wouters",
            "models/smets_wouters_orig/smets_wouters_orig",
            "models/m990/m990",
            "models/m1002/m1002",
            "models/m1010/m1010",
            "models/poolmodel/poolmodel",
            "data/fred_data",
            "data/load_data",
            "data/load_data_poolmodel",
            "data/misc",
            "data/reverse_transform",
            "data/simulate_data",
            "data/transformations",
            "data/transform_data",
            "data/util",
            "solve/gensys",
            "solve/solve",
            "solve/solve_poolmodel",
            "estimate/filter",
            "estimate/cat",
            "estimate/posterior",
            "estimate/filter_poolmodel",
            "estimate/posterior_poolmodel",
            "estimate/estimate_bma",
            "estimate/hessian",
            "estimate/csminwel",
            "estimate/optimize",

            "estimate/metropolis_hastings",

            # "estimate/smc/smc",
            # "estimate/smc/helpers",
            # "estimate/smc/initialization",
            # "estimate/smc/util",
            # "estimate/smc/mutation",
            # "estimate/smc/resample",

            "forecast/drivers",
            "forecast/smooth",
            "forecast/forecast",
            "forecast/shock_decompositions",
            "forecast/impulse_responses",
            "forecast/io",
            "forecast/forecast_one",
            "forecast/util",
            "analysis/compute_meansbands",
            "analysis/df_to_table",
            "analysis/io",
            "analysis/meansbands",
            "analysis/moments",
            "analysis/util",
            "altpolicy/altpolicy",
            "altpolicy/taylor93",
            "altpolicy/taylor99",
            "scenarios/scenario",
            "scenarios/forecast",
            "scenarios/switching",
            "scenarios/drivers",
            "decomp/decompose_forecast",
            "plot/plot",
  	        "plot/util"
            #"packet/packet"
            ]

for test in my_tests
    test_file = string("$test.jl")
    @printf " * %s\n" test_file
    include(test_file)
end


text:
# using ScikitLearn
# @sk_import covariance: MinCovDet

import ..Utilities
abstract type SdErgm end

Base.string(x::SdErgm) = DynNets.name(x)
export string

function get_option(model::SdErgm, optName)

    defaultOptions = typeof(model)(model.staticModel).options

    if haskey(model.options, optName)
        return getindex(model.options, optName)

    elseif haskey(model.defaultOptions, optName)
        return getindex(defaultOptions, optName)
    else
        error("No option found")
    end
end

identify(model::SdErgm,UnPar::Array{<:Real,1}, idType ) = StaticNets.identify(model.staticModel, UnPar; idType = idType)

number_ergm_par(model::T where T <:SdErgm) = model.staticModel.nErgmPar

type_of_obs(model::SdErgm) =  StaticNets.type_of_obs(model.staticModel)

stats_from_mat(model::SdErgm, A ::Matrix{<:Real}) = StaticNets.stats_from_mat(model.staticModel, A ::Matrix{<:Real})

"""
Model to use for utilities and eventually fast sampling. The main use is to point pmle models to their mle counterpart, which can be used for sampling or for comparison purposes
"""
reference_model(model::SdErgm) = model

sample_ergm_sequence(model::SdErgm, N, parVecSeq_T::Matrix, nSample) = StaticNets.sample_ergm_sequence(model.staticModel, N, parVecSeq_T, nSample)

function seq_of_obs_from_seq_of_mats(model::T where T <:SdErgm, AT_in)

    AT = convert_to_array_of_mats(AT_in)
    T = length(AT)
    obsT = Array{type_of_obs(model), 1}(undef, T)
    for t in 1:T
        obsT[t] = DynNets.stats_from_mat(model, AT[t])
    end
    return obsT
end

#region options and conversions of parameters for optimization
function setOptionsOptim(model::T where T<: SdErgm; show_trace = false)
    "Set the options for the optimization required in the estimation of the model.
    For the optimization use the Optim package."
    tol = eps()*10
    maxIter = 300
    call_back = tr -> isnan(tr[end].value) |(!isfinite(tr[end].value))
    opt = Optim.Options(  g_tol = 1e-8,
                     x_tol = tol,
                     x_abstol = tol,
                     x_reltol = tol,
                     f_tol = tol,
                     f_reltol = tol,
                     f_abstol = tol,
                     iterations = maxIter,
                     show_trace = show_trace ,#false,#
                     store_trace = true ,#false,#
                     callback = call_back,
                     show_every=5)

    algo = NewtonTrustRegion(; initial_delta = 0.1,
                    delta_hat = 0.2,
                    eta = 0.1,
                    rho_lower = 0.25,
                    rho_upper = 0.75)
    algo = Newton(; alphaguess = LineSearches.InitialHagerZhang(),
    linesearch = LineSearches.BackTracking())
      return opt, algo
end

"""
Convert vector conaining both SD and constant parameters into array with one vector of static parameters for each ergm parameter. The vector can have eithre 1 or 3 elements, depending on wether the element is static or dynamical
"""
function vec_2_array_all_par(model::SdErgm, VecGasPar::Array{<:Real,1}, indTvPar :: BitArray{1})

    Npar = length(indTvPar)
    arrayAllPar = [zeros(Real, 1 + tv*2) for tv in indTvPar]#Array{Array{<:Real,1},1}(undef, Npar)
    last_i = 1

    for i=1:Npar
        nStatPerDyn = 1 + 2*indTvPar[i]
        arrayAllPar[i][:] = VecGasPar[last_i: last_i+nStatPerDyn-1]
            last_i = last_i + nStatPerDyn
    end
    return arrayAllPar
end

"""
The inverse of vec_2_array_all_par
"""
function array_2_vec_all_par(model::SdErgm, ArrayGasPar, indTvPar :: BitArray{1})

        NTvPar = sum(indTvPar)
        Npar = length(indTvPar)
        NgasPar = 2*NTvPar + Npar

        VecGasPar = zeros(NgasPar)
        last = 1
        for i=1:Npar
            VecGasPar[last:last+length(ArrayGasPar[i])-1] = ArrayGasPar[i]
            last = last + length(ArrayGasPar[i])
        end
        return VecGasPar
end

"""
Given vecAllPar divide it into a vector of Score Driven parameters and one of costant parameters
"""
function divide_SD_par_from_const(model::T where T <:SdErgm,  vecAllPar::Array{<:Real,1}; indTvPar = model.indTvPar)

    nTvPar = sum(indTvPar)
    nErgmPar = length(indTvPar)

    vecSDParAll = zeros(Real,3nTvPar )
    vConstPar = zeros(Real,nErgmPar-nTvPar)

    lastInputInd = 0
    lastIndSD = 0
    lastConstInd = 0
    #extract the vector of static SD parameters
    for i=1:nErgmPar
        if indTvPar[i]
            vecSDParAll[lastIndSD+1] = vecAllPar[lastInputInd + 1]
            vecSDParAll[lastIndSD+2] = vecAllPar[lastInputInd + 2]
            vecSDParAll[lastIndSD+3] = vecAllPar[lastInputInd + 3]
            lastInputInd +=3
            lastIndSD +=3
        else
            vConstPar[lastConstInd+1] = vecAllPar[lastInputInd  + 1]
            lastInputInd +=1
            lastConstInd +=1
        end
    end
    return vecSDParAll, vConstPar
end

function merge_SD_par_and_const(model::T where T <:SdErgm, vecSDPar::Array{<:Real,1}, vConstPar, indTvPar=model.indTvPar)

    nTvPar = sum(indTvPar)
    nConstPar = sum(.!indTvPar)
    nConstPar == length(vConstPar) ? () : error()

    nErgmPar = length(indTvPar)
    nAllPar = 3*nTvPar + nConstPar

    vecAllPar = zeros(Real, nAllPar)

    lastInputInd = 0
    lastConstInd = 0
    lastIndAll = 0
    lastIndSD = 0
    for i=1:nErgmPar
        if indTvPar[i]

            vecAllPar[lastIndAll+1] = vecSDPar[lastIndSD + 1]
            vecAllPar[lastIndAll+2] = vecSDPar[lastIndSD + 2]
            vecAllPar[lastIndAll+3] = vecSDPar[lastIndSD + 3]

            lastIndAll +=3
            lastIndSD +=3
        else
            vecAllPar[lastIndAll+1] = vConstPar[lastConstInd + 1]

            lastIndAll +=1
            lastInputInd +=1
            lastConstInd +=1
        end
    end
    return vecAllPar
end

"""
Restrict the  Score Driven parameters  to appropriate link functions to ensure that they remain in the region where the SD dynamics is well specified (basically 0<=B<1  A>=0)
"""
function restrict_SD_static_par(model::T where T <:SdErgm, vecUnSDPar::Array{<:Real,1})

    nSDPar = length(vecUnSDPar)
    nTvPar, rem = divrem(nSDPar,3)

    rem == 0 ? () : error()

    arrayOfVecsReSd = [ [vecUnSDPar[i], link_R_in_0_1(vecUnSDPar[i+1]), link_R_in_R_pos(vecUnSDPar[i+2]) ] for i in 1:3:nSDPar]

    vecReSDPar = reduce(vcat, arrayOfVecsReSd)

    return vecReSDPar
end

"""
    inverse of restrict_SD_static_par
"""
function unrestrict_SD_static_par(model::T where T <:SdErgm, vecReSDPar::Array{<:Real,1})

    nSDPar = length(vecReSDPar)
    nTvPar, rem = divrem(nSDPar,3)

    rem == 0 ? () : error()

    arrayOfVecsUnSd = [ [vecReSDPar[i], inv_link_R_in_0_1(vecReSDPar[i+1]), inv_link_R_in_R_pos(vecReSDPar[i+2]) ] for i in 1:3:nSDPar]

    vecUnSDPar = reduce(vcat, arrayOfVecsUnSd)

    return vecUnSDPar
end

"""
    separate SD parameters from constant ones, unrestrict SD and merge them back
"""
function unrestrict_all_par(model::T where T <:SdErgm, vAllPar)

    vSDRe, vConst = divide_SD_par_from_const(model, vAllPar)

    vSDUn = unrestrict_SD_static_par(model, vSDRe)

    merge_SD_par_and_const(model, vSDUn, vConst)
end

"""
    separate SD parameters from constant ones, restrict SD and merge them back
"""
function restrict_all_par(model::T where T <:SdErgm, vAllPar)

    vSDUn, vConst = divide_SD_par_from_const(model, vAllPar)

    vSDRe = restrict_SD_static_par(model, vSDUn)

    merge_SD_par_and_const(model, vSDRe, vConst)
end

"""
Map subset of parameters for integrated SD model (only A parameters and constant ergm parameters) to a full vector of parameters containing all static ergm parameters, w, B, and A. Setting w.=0 and B .=1
"""
function par_integr_2_full_unrestr(model, vecUnParIntegratedAll)
    NTvPar = sum(model.indTvPar)
    vecUnSDParA = vecUnParIntegratedAll[model.indTvPar]
    vConstPar = vecUnParIntegratedAll[.!model.indTvPar]

    vecUnSDPar = zeros(Real, 3*NTvPar)

    vecUnSDPar[1:3:end] .= 0
    vecUnSDPar[2:3:end] .= Inf
    vecUnSDPar[3:3:end] .= vecUnSDParA
    vecUnPar = merge_SD_par_and_const(model, vecUnSDPar, vConstPar)

    return vecUnPar
end

"""
Given the flag of constant parameters, a starting value for their unconditional means (their constant value, for those constant), return a starting point for the optimization
"""
function starting_point_optim(model::T where T <:SdErgm,  UM; indTargPar =  falses(100))
    indTvPar=model.indTvPar
    nTvPar = sum(indTvPar)
    nErgmPar = length(indTvPar)

    # #set the starting points for the optimizations
    B0_Re  = 0.98; B0_Un = log(B0_Re ./ (1 .- B0_Re ))
    ARe_min =0.0000000001
    A0_Re  = 0.01 ; A0_Un = log(A0_Re  .-  ARe_min)

    # starting values for the vector of parameters that have to be optimized
    vParOptim_0 = zeros(nErgmPar + nTvPar*2 )
    last = 0
    for i=1:nErgmPar
        if indTvPar[i]
            if indTargPar[i]
                vParOptim_0[last+1:last+2] = [ B0_Un; A0_Un]
                last+=2
            else
                vParOptim_0[last+1:last+3] = [UM[i]*(1 .- B0_Re) ; B0_Un; A0_Un]
                last+=3
            end
        else
            vParOptim_0[last+1] = UM[i]
            last+=1
        end
    end
    @debug "[starting_point_optim][UM = $UM, vParOptim_0 = $vParOptim_0, ARe_min = $ARe_min]"
    return vParOptim_0, ARe_min
end

convert_to_array_of_mats(AT::Array{Array{<:Any, 2}, 1}) = AT
convert_to_array_of_mats(AT::Array{<:Any, 1}) = [m for m in AT]
convert_to_array_of_mats(AT::Array{<:Any, 3}) = [AT[:,:,t] for t in 1:size(AT)[3]]

#endregion

target_function_t(model::SdErgm, obs_t, N, f_t) = StaticNets.obj_fun(model.staticModel, obs_t, N, f_t)

function target_function_t_grad(model::T where T<: SdErgm, obs_t, N, f_t)

    target_fun_t(x) = target_function_t(model, obs_t, N, x)

    grad_tot_t = ForwardDiff.gradient(target_fun_t, f_t)

    return grad_tot_t
end

function target_function_t_hess(model::T where T<: SdErgm, obs_t, N, f_t)

    target_fun_t(x) = target_function_t(model, obs_t, N, x)

    hess_tot_t = ForwardDiff.hessian(target_fun_t, f_t)

    return hess_tot_t
end

function target_function_t_fisher(model::T where T<: SdErgm, obs_t, N, f_t)

    error("need to code the fisher manually, cannot be obtained by Automatic Differentiation in the general case, only when the fisher equality holds")
    return fisher_info_t
end

function updatedGasPar( model::T where T<: SdErgm, obs_t, N, ftot_t::Array{<:Real,1}, I_tm1::Array{<:Real,2}, indTvPar::BitArray{1}, Wgas::Array{<:Real,1}, Bgas::Array{<:Real,1}, Agas::Array{<:Real,1})

    #= likelihood and gradients depend on all the parameters (ftot_t), but
    only the time vaying ones (f_t) are to be updated=#

    target_fun_val_t = target_function_t(model, obs_t, N, ftot_t)

    grad_tot_t = target_function_t_grad(model, obs_t, N, ftot_t)

    I_reg = I(size(I_tm1)[1])

    if model.scoreScalingType =="HESS_D"
        hess_tot_t = target_function_t_hess(model, obs_t, N, ftot_t)
        I_updt = -hess_tot_t
        I_t = I_updt
        s_t = (grad_tot_t./I_t[I(size(I_t)[1])])[indTvPar]
    elseif model.scoreScalingType =="MAX_LINKS"
        hess_tot_t = target_function_t_hess(model, obs_t, N, ftot_t)
        I_t = ones(size(grad_tot_t)).*(N^2-N)
        s_t = grad_tot_t[indTvPar]./I_t
    elseif model.scoreScalingType =="HESS"
        hess_tot_t = target_function_t_hess(model, obs_t, N, ftot_t)
        I_updt = -hess_tot_t
        I_t = I_updt
        s_t = (I_t\grad_tot_t)[indTvPar]

    elseif model.scoreScalingType =="FISH_D"
        extremeScaledScore = 200
        fish_tot_t = target_function_t_fisher(model, obs_t, N, ftot_t)
        I_updt = fish_tot_t
        I_t = I_updt
        s_t = clamp.( (grad_tot_t./sqrt.((I_t[I(size(I_t)[1])] )))[indTvPar], -extremeScaledScore, extremeScaledScore)

    end

    f_t = ftot_t[indTvPar] #Time varying ergm parameters
    f_tp1 = Wgas .+ Bgas.* f_t .+ Agas.*s_t

    ftot_tp1 = copy(ftot_t)
    ftot_tp1[indTvPar] = f_tp1 #of all parameters udate dynamic ones with GAS

    #any(isnan.(ftot_tp1)) ? error() : ()

    # I_t is the inverse of the S_t matrix in typical notations

    return ftot_tp1, target_fun_val_t, I_t, grad_tot_t
end

function score_driven_filter( model::T where T<: SdErgm, N, obsT, vResGasPar::Array{<:Real,1}, indTvPar::BitArray{1}; vConstPar ::Array{<:Real,1} = zeros(Real,2), ftot_0::Array{<:Real,1} = zeros(Real,2))

    nErgmPar = number_ergm_par(model)
    NTvPar   = sum(indTvPar)
    T= length(obsT)

    # Organize parameters of the GAS update equation
    Wvec = vResGasPar[1:3:3*NTvPar]
    Bvec = vResGasPar[2:3:3*NTvPar]
    Avec = vResGasPar[3:3:3*NTvPar]

    #compute unconditional means according to SD dgp
    UMallPar = zeros(Real,nErgmPar)
    UMallPar[indTvPar] =  Wvec ./ (1 .- Bvec)
    if !all(indTvPar) # if not all parameters are time varying
        UMallPar[.!indTvPar] = vConstPar
    end

    sum(ftot_0)==0 ? ftot_0 = UMallPar : ()# identify(model,UMallNodesIO)

    fVecT = ones(Real,nErgmPar,T)
    sVecT = ones(Real,nErgmPar,T)
    invScalMatT = ones(Real,nErgmPar, nErgmPar,T)
    logLikeVecT = ones(Real,T)

    if NTvPar==0
        I_tm1 = ones(1,1)
    else
        I_tm1 = Float64.(Diagonal{Real}(I,NTvPar))
    end

    fVecT[:,1] = ftot_0
    if sum(.!indTvPar) >0
        fVecT[.!indTvPar,1] = vConstPar
    end

    ftot_tp1 = zeros(nErgmPar)
    for t=1:T

        N_t = Utilities.get_N_t(N, t)

        obs_t = obsT[t]

        #obj fun at time t is objFun(obs_t, f_t)

        # predictive step
        ftot_tp1, logLikeVecT[t], invScalMatT[:,:,t], sVecT[:,t] = updatedGasPar(model, obs_t, N_t, fVecT[:, t], I_tm1, indTvPar, Wvec, Bvec, Avec)

        I_tm1 = invScalMatT[:,:,t]

        # would the following be the update step instead ??
        #fVecT[:,t], loglike_t, I_tm1, grad_t = updatedGasPar(model, obs_t, N, fVecT[:,t-1], I_tm1, indTvPar, Wvec, Bvec, Avec)

        if t!=T
            fVecT[:,t+1] = ftot_tp1
        end
    end

    return fVecT::Array{<:Real, 2}, logLikeVecT::Vector{<:Real}, sVecT::Array{<:Real,2}, invScalMatT::Array{<:Real, 3}
end

function score_driven_dgp( model::T where T<: SdErgm, N, dgpNT, vResGasPar::Array{<:Real,1}, indTvPar::BitArray{1}; vConstPar ::Array{<:Real,1} = zeros(Real,2), ftot_0::Array{<:Real,1} = zeros(Real,2))

    nErgmPar = number_ergm_par(model)
    NTvPar   = sum(indTvPar)

    T = dgpNT[2]
    N_const = dgpNT[1]
    N_T = ones(Int, T) .* N_const
    N = N_const

    # Organize parameters of the GAS update equation
    Wvec = vResGasPar[1:3:3*NTvPar]
    Bvec = vResGasPar[2:3:3*NTvPar]
    Avec = vResGasPar[3:3:3*NTvPar]

    # start values equal the unconditional mean,and  constant ones remain equal to the unconditional mean, hence initialize as:
    UMallPar = zeros(Real,nErgmPar)
    UMallPar[indTvPar] =  Wvec ./ (1 .- Bvec)
    if !all(indTvPar) # if not all parameters are time varying
        UMallPar[.!indTvPar] = vConstPar
    end

    fVecT = ones(Real,nErgmPar,T)
    sVecT = ones(Real,nErgmPar,T)
    invScalMatT = ones(Real,nErgmPar, nErgmPar,T)
    logLikeVecT = ones(Real,T)

    sum(ftot_0)==0 ? ftot_0 = UMallPar : ()# identify(model,UMallNodesIO)

    if NTvPar==0
        I_tm1 = ones(1,1)
    else
        I_tm1 = Float64.(Diagonal{Real}(I,NTvPar))
    end

    loglike = 0

    A_T = zeros(Int8, N, N, T)

    fVecT[:,1] = ftot_0

    ftot_tp1 = zeros(nErgmPar)
    for t=1:T
    #    println(t)

        diadProb = StaticNets.diadProbFromPars(StaticNets.ErgmDirBin0Rec0(), fVecT[:, t] )

        A_t = StaticNets.samplSingMatCan(StaticNets.ErgmDirBin0Rec0(), diadProb, N)

        A_T[:, : , t] = A_t

        obs_t = stats_from_mat(model, A_t)

        #obj fun at time t is objFun(obs_t, f_t)

        # predictive step
        ftot_tp1, logLikeVecT[t], invScalMatT[:,:,t], sVecT[:,t] = updatedGasPar(model, obs_t, N, fVecT[:, t], I_tm1, indTvPar, Wvec, Bvec, Avec)

        I_tm1 = invScalMatT[:,:,t]

        # would the following be the update step instead ??
        #fVecT[:,t], loglike_t, I_tm1, grad_t = updatedGasPar(model, obs_t,  N, fVecT[:,t-1], I_tm1, indTvPar, Wvec, Bvec, Avec)

        if t!=T
            fVecT[:,t+1] = ftot_tp1
        end
    end

    return fVecT::Array{Float64, 2}, A_T, sVecT, invScalMatT

end

estimate_single_snap_sequence(model::T where T<: SdErgm, obsT, aggregate=0) = StaticNets.estimate_sequence(model.staticModel, obsT)

"""
Sequence of loglikelihoods from SD filter as a function of unrestricted parameters (vector of both static ergm and static SD) and vector of initial parameters
"""
function seq_loglike_sd_filter(model::SdErgm, N, obsT,  vecUnPar::Array{<:Real,1}, ftot_0::T where T<:AbstractArray)

        vecUnSDPar, vConstPar = divide_SD_par_from_const(model, vecUnPar)

        oneInADterms  = (StaticNets.maxLargeVal + vecUnPar[1])/StaticNets.maxLargeVal

        vecReSDPar = restrict_SD_static_par(model, vecUnSDPar)

        foo, logLikeVecT, foo1 = score_driven_filter( model, N, obsT, vecReSDPar, model.indTvPar;  vConstPar =  vConstPar, ftot_0 = ftot_0 .* oneInADterms)

        return logLikeVecT
end

function seq_loglike_sd_filter(model, N, obsT,  vecUnPar::Array{<:Real,1}, ftot_0Fun::Function)
    @debug "[seq_loglike_sd_filter][N=$N, obs_1_sample = $(obsT[1][1:5,:]), vecUnPar=$vecUnPar, ftot_0Fun = $ftot_0Fun ]"

        vecUnSDPar, vConstPar = divide_SD_par_from_const(model, vecUnPar)

        vecReSDPar = restrict_SD_static_par(model, vecUnSDPar)

        ftot_0 = ftot_0Fun(vecReSDPar, vConstPar)

        return seq_loglike_sd_filter(model, N, obsT, vecUnPar, ftot_0)
end

"""
Estimate the SD and ergm static parameters
"""
function estimate(model::T where T<: SdErgm, N, obsT; indTvPar::BitArray{1}=model.indTvPar, indTargPar::BitArray{1} = falses(length(model.indTvPar)), initMeth :: String = get_option(model, "initMeth"), ftot_0Fixed :: Array{<:Real,1} = zeros(2), nObsFirstEst=10,  vParOptim_0 =zeros(2), shuffleObsInds::Union{Nothing, Vector{<:Int}} = nothing, show_trace = false )
    @debug "[estimate][start][model=$model, indTvPar=$indTvPar, indTargPar=$indTargPar, initMeth=$initMeth, ftot_0Fixed = $ftot_0Fixed,  vParOptim_0 = $vParOptim_0, shuffleObsInds=$shuffleObsInds, model.options = $(model.options)]"

    T = length(obsT);
    nErgmPar = number_ergm_par(model)
    NTvPar = sum(indTvPar)
    Logging.@debug( "[estimate][Estimating N = $N , T=$T]")

    optims_opt, algo = setOptionsOptim(model; show_trace = show_trace )

    staticPars = static_estimate(model, obsT)

    if sum(vParOptim_0) == 0
        vParOptim_0_tmp, ARe_min = starting_point_optim(model, staticPars; indTargPar = indTargPar)
        vParOptim_0 = vParOptim_0_tmp
        @debug "[estimate][using standard starting point][$(vParOptim_0_tmp)]"
    end

    ftot_0Fun(vecReSDPar, vConstPar) = ftot_0Fixed
    @show ftot_0Fixed

    # How should we handle the initialization of SD iteration ?
    if (initMeth == "uncMean") & (get_option(model,  "integrated") )

        error("Integrated filters cannot use the unconditional mean, using $initMeth instead")
    end

    if (initMeth == "uncMean") & (!get_option(model,  "integrated"))

        function ftot_0Fun(vecReSDPar, vConstPar)
            w = vecReSDPar[1:3:end]
            B = vecReSDPar[2:3:end]

            # unconditional mean of SD parameters
            UmSDPar  = unconditional_mean_auto_regr(w, B)

            ftot_0 = zeros(Real, nErgmPar)

            ftot_0[indTvPar] = UmSDPar
            ftot_0[.!indTvPar] = vConstPar
            return ftot_0
        end
    else
        if initMeth == "estimateFirstObs"
            ftot_0Fixed =  static_estimate(model, obsT[1:nObsFirstEst])
            ftot_0Fun(vecReSDPar, vConstPar) = ftot_0Fixed

        elseif typeof(initMeth) <: AbstractArray

            if length(ftot_0Fixed) == nErgmPar

                ftot_0Fixed =  static_estimate(model, obsT[1:nObsFirstEst])
                ftot_0Fun(vecReSDPar, vConstPar) = ftot_0Fixed
            else
                error("wrong length of initial ergm parameters")
            end
        end
    end

    #define the objective function for the optimization

    if !isnothing(shuffleObsInds)
       # was observations shuffing required ?  (used for explorative analysis on non-parametric bootstrap)
        length(shuffleObsInds) == T ? () : error("wrong lenght of shuffling indices")

    else
        # No shuffling required
        shuffleObsInds = collect(1:T)
    end

    #shuffled observations loglikelihood
    function objfunSdOptShuffled(vecUnPar::Array{<:Real,1})

        logLikeVecT = seq_loglike_sd_filter(model, N, obsT, vecUnPar, ftot_0Fun)
        return - sum(logLikeVecT[shuffleObsInds])
    end

    ADobjfunSdOpt = TwiceDifferentiable(objfunSdOptShuffled, vParOptim_0; autodiff = :forward);

    if haskey(model.options, "Firth")
        if model.options["Firth"]
            # log likelihood with firth additional term
            function objfunSdOptFirth(vecUnPar::Array{<:Real,1})
                objfunSdOpt(vecUnPar) -  0.5 * LinearAlgebra.logabsdet( FiniteDiff.finite_difference_hessian(objfunSdOpt, vecUnPar))[1]
            end
            ADobjfunSdOpt = TwiceDifferentiable(objfunSdOptFirth, vParOptim_0; autodiff = :forward);
        end
    end

    if get_option(model,  "integrated")
        # should we estimate  an integrated version of the SD filter?
        vParOptim_0 = ones(Real, NTvPar) * 0.001 #vParOptim_0[3:3:end]

        objfunSdOptInt(vecUnParIntegratedAll::Array{<:Real,1}) = objfunSdOptShuffled(par_integr_2_full_unrestr(model, vecUnParIntegratedAll))

        ADobjfunSdOpt = TwiceDifferentiable(objfunSdOptInt, vParOptim_0; autodiff = :forward);
    end

    #Run the optimization

    Logging.@debug("[estimate][Starting point for Optim $vParOptim_0]")
    Logging.@debug("[estimate][Starting point for Optim $(restrict_all_par(model, vParOptim_0))]")
    optim_out2  = optimize(ADobjfunSdOpt, vParOptim_0, algo, optims_opt)
    outParAllUn = Optim.minimizer(optim_out2)
    if get_option(model,  "integrated")
        outParAllUnFull = zeros(3*NTvPar)
        outParAllUnFull[1:3:end] .= 0
        outParAllUnFull[2:3:end] .= Inf
        outParAllUnFull[3:3:end] .= outParAllUn
        outParAllUn = deepcopy(outParAllUnFull)
    end

    Logging.@debug(optim_out2)

    vecAllReParHat = restrict_all_par(model, outParAllUn)

    vecParGasHat, vecParConstHat = divide_SD_par_from_const(model, vecAllReParHat)

    Logging.@debug("[estimate][Final paramters SD: $vecParGasHat , constant: $vecParConstHat ]")

    arrayAllParHat = vec_2_array_all_par(model, vecAllReParHat, indTvPar)

    ftot_0Hat = ftot_0Fun(vecParGasHat, vecParConstHat)

    conv_flag =  Optim.converged(optim_out2)

    @debug "[estimate][end]"
    return  (; arrayAllParHat, conv_flag, UM =zeros(2), ftot_0Hat)

end

ext: jl


text:
export ContinuousUnivariateLogPdf

using Distributions

import DomainSets
import DomainIntegrals
import Base: isapprox

"""
    ContinuousUnivariateLogPdf{ D <: DomainSets.Domain, F } <: ContinuousUnivariateDistribution

Generic continuous univariate distribution in a form of domain specification and logpdf function. Can be used in cases where no
known analytical distribution available.

# Arguments
- `domain`: domain specificatiom from `DomainSets.jl` package
- `logpdf`: callable object that represents a `logpdf` of a distribution. Does not necessarily normalised.

```julia
fdist = ContinuousUnivariateLogPdf(DomainSets.FullSpace(), (x) -> -x^2)
```
"""
struct ContinuousUnivariateLogPdf{ D <: DomainSets.Domain, F } <: ContinuousUnivariateDistribution
    domain :: D
    logpdf :: F
end

ContinuousUnivariateLogPdf(f::Function) = ContinuousUnivariateLogPdf(DomainSets.FullSpace(), f)

(dist::ContinuousUnivariateLogPdf)(x::Real)                      = logpdf(dist, x)
(dist::ContinuousUnivariateLogPdf)(x::AbstractVector{ <: Real }) = logpdf(dist, x)

Distributions.support(dist::ContinuousUnivariateLogPdf) = Distributions.RealInterval(DomainSets.infimum(dist.domain), DomainSets.supremum(dist.domain))

Distributions.mean(dist::ContinuousUnivariateLogPdf)    = error("mean() is not defined for `ContinuousUnivariateLogPdf`.")
Distributions.median(dist::ContinuousUnivariateLogPdf)  = error("median() is not defined for `ContinuousUnivariateLogPdf`.")
Distributions.mode(dist::ContinuousUnivariateLogPdf)    = error("mode() is not defined for `ContinuousUnivariateLogPdf`.")
Distributions.var(dist::ContinuousUnivariateLogPdf)     = error("var() is not defined for `ContinuousUnivariateLogPdf`.")
Distributions.std(dist::ContinuousUnivariateLogPdf)     = error("std() is not defined for `ContinuousUnivariateLogPdf`.")
Distributions.cov(dist::ContinuousUnivariateLogPdf)     = error("cov() is not defined for `ContinuousUnivariateLogPdf`.")
Distributions.invcov(dist::ContinuousUnivariateLogPdf)  = error("invcov() is not defined for `ContinuousUnivariateLogPdf`.")
Distributions.entropy(dist::ContinuousUnivariateLogPdf) = error("entropy() is not defined for `ContinuousUnivariateLogPdf`.")

# We don't expect neither `pdf` nor `logpdf` to be normalised
Distributions.pdf(dist::ContinuousUnivariateLogPdf, x::Real)    = exp(logpdf(dist, x))

function Distributions.logpdf(dist::ContinuousUnivariateLogPdf, x::Real)
    @assert x ∈ dist.domain "x = $(x) does not belong to the domain of $dist"
    return dist.logpdf(x)
end

# These are fallbacks for various optimisation packages which may pass arguments as vectors
function Distributions.pdf(dist::ContinuousUnivariateLogPdf, x::AbstractVector{ <: Real })
    @assert length(x) === 1 "`ContinuousUnivariateLogPdf` expects either float or a vector of a single float as an input for the `pdf` function."
    return exp(logpdf(dist, first(x)))
end

function Distributions.logpdf(dist::ContinuousUnivariateLogPdf, x::AbstractVector{ <: Real })
    @assert length(x) === 1 "`ContinuousUnivariateLogPdf` expects either float or a vector of a single float as an input for the `logpdf` function."
    return logpdf(dist, first(x))
end

Base.precision(dist::ContinuousUnivariateLogPdf) = error("precision() is not defined for `ContinuousUnivariateLogPdf`.")

Base.convert(::Type{ ContinuousUnivariateLogPdf }, domain::D, logpdf::F) where { D <: DomainSets.Domain, F } = ContinuousUnivariateLogPdf{D, F}(domain, logpdf)

convert_eltype(::Type{ ContinuousUnivariateLogPdf }, ::Type{ T }, dist::ContinuousUnivariateLogPdf) where { T <: Real } = convert(ContinuousUnivariateLogPdf, dist.domain, dist.logpdf)

vague(::Type{ <: ContinuousUnivariateLogPdf }) = ContinuousUnivariateLogPdf(DomainSets.FullSpace(), (x) -> 1.0)

prod_analytical_rule(::Type{ <: ContinuousUnivariateLogPdf }, ::Type{ <: ContinuousUnivariateLogPdf }) = ProdAnalyticalRuleAvailable()

function prod(::ProdAnalytical, left::ContinuousUnivariateLogPdf, right::ContinuousUnivariateLogPdf)
    @assert left.domain == right.domain "Different domain types in product of generic `ContinuousUnivariateLogPdf` distributions. Left domain is $(left.domain), right is $(right.domain)."
    plogpdf = let left = left, right = right
        (x) -> logpdf(left, x) + logpdf(right, x)
    end
    return ContinuousUnivariateLogPdf(left.domain, plogpdf)
end

# This method is inaccurate and relies on various approximation methods, which may fail in different scenarious
# Current implementation of `isapprox` method supports only FullSpace and HalfLine domains with limited accuracy
function Base.isapprox(left::ContinuousUnivariateLogPdf, right::ContinuousUnivariateLogPdf; kwargs...)
    if left.domain !== right.domain
        return false
    end
    return culogpdf__isapprox(left.domain, left, right; kwargs...)
end

# https://en.wikipedia.org/wiki/Gauss–Hermite_quadrature
function culogpdf__isapprox(domain::DomainSets.FullSpace, left::ContinuousUnivariateLogPdf, right::ContinuousUnivariateLogPdf; kwargs...)
    return isapprox(zero(eltype(domain)), DomainIntegrals.integral(Q_GaussHermite(32), (x) -> exp(x ^ 2) * abs(left(x) - right(x))); kwargs...)
end

# https://en.wikipedia.org/wiki/Gauss–Laguerre_quadrature
function culogpdf__isapprox(domain::DomainSets.HalfLine, left::ContinuousUnivariateLogPdf, right::ContinuousUnivariateLogPdf; kwargs...)
    return isapprox(zero(eltype(domain)), DomainIntegrals.integral(Q_GaussLaguerre(32), (x) -> exp(x) * abs(left(x) - right(x))); kwargs...)
end

# We do not check typeof of a different functions because in most of the cases lambdas have different types, but it does not really mean that objects are different
function is_typeof_equal(left::ContinuousUnivariateLogPdf{D, F1}, right::ContinuousUnivariateLogPdf{D, F2}) where { D, F1 <: Function, F2 <: Function }
    return true
end


text:
using Test
include("../src/classification.jl")
using Main.Classification

@testset "Construct NeuralNet" begin
    X = [1 2 3 4 5
         1 2 3 4 5
         0 1 1 2 5]
    y = [1 1 0 1 1]
    @testset "Weave" begin
        net = NeuralNet([6, 10, 1])
        @test length(net.layers) == 3
        weave!(net, X, y)
        @test length(net.layers) == 3
        @test size(net.layers[1].b) == (6,)
        @test size(net.layers[2].b) == (10,)
        @test size(net.layers[3].b) == (1,)
        @test size(net.layers[1].W) == (6, 3)
        @test size(net.layers[2].W) == (10, 6)
        @test size(net.layers[3].W) == (1, 10)
    end
    @testset "FeedForward" begin
        net = NeuralNet([5, 2, 1])
        weave!(net, X, y)
        out = feedforward(net, X)
        @test size(out) == size(y)
    end
    @testset "FeedForward 2" begin
        net = NeuralNet([5, 2, 3])
        weave!(net, X, y)
        out = feedforward(net, X)
        @test size(out) == (3, size(y, 2))
    end
    # I messed up the types, causing these to fail
    # @testset "Backward Propogation" begin
        # net = NeuralNet([8, 6, 1])
        # weave!(net, X, y)
        # backpropagation(net, X, y)
        # @test true
    # end
    # @testset "Gradient Descent" begin
        # net = NeuralNet([8, 6, 1])
        # gd = GradientDescent()
        # fit!(net, gd, X, y)
        # @test true
    # end
end


text:
importall Bukdu

type UserController <: ApplicationController
end

index(::UserController) = render(JSON, "hello")

Router() do
    scope("/api") do
        resources("/users", UserController)
    end
end

using Base.Test
conn = (Router)(get, "/api/users")
@test 200 == conn.status
@test "application/json" == conn.resp_header["Content-Type"]
@test """\"hello\"""" == conn.resp_body

import HttpCommon: Request, Response
req = Request()
req.method = "GET"
req.resource = "/api/users"
res = Bukdu.handler(req, Response())
@test 200 == res.status
@test "application/json" == res.headers["Content-Type"]
@test """\"hello\"""" == String(res.data)

logs = []
before(render, JSON) do t
    push!(logs, "b $t")
end
after(render, JSON) do t
    push!(logs, "a $t")
end

conn = (Router)(get, "/api/users")
@test ["b hello", "a hello"] == logs

layout(::Layout, body) = """[$body]"""
show(::UserController) = render(JSON/Layout, "hello")
before(render, JSON/Layout) do t
    push!(logs, "bl $t")
end
after(render, JSON/Layout) do t
    push!(logs, "al $t")
end

empty!(logs)
req.method = "GET"
req.resource = "/api/users/1"
res = Bukdu.handler(req, Response())
@test 200 == res.status
@test "application/json" == res.headers["Content-Type"]
@test "[\"hello\"]" == String(res.data)
@test ["bl hello","b hello","a hello","al hello"] == logs


text:
# SPDX-License-Identifier: MIT

using CUDA
using CUDA.CUSPARSE
using JLD2
using MarkovModels
using SparseArrays
using LinearAlgebra

function main(fsmfile, T, B, N; warmup)
    # Load the Finite State Machine (FSM).
    fsm = jldopen(fsmfile, "r") do f
       f["fsm"]
    end

    # Number states in the FSM.
    S = length(fsm)

    # Create a batch of the same FSM.
    fsm = union(repeat([fsm], B)...)

    # Number of pdfs associated with the FSM. This is specific to the FSM
    # and should not be changed unless you use your own FSM.
    npdfs = 84

    # Length of the sequences of the batch. Here we assume that all
    # the sequence have the same length.
    lengths = repeat([N], B)

    # Generate some pseudo-(log)likelihood
    lhs = ones(T, npdfs, N, B)

    # Computing precision
    precision = T == Float64 ? "double" : "single"

    # Dummy run to force compilation
    if warmup pdfposteriors(fsm, lhs, lengths) end

    etime = time_ns()
    pdfposteriors(fsm, lhs, lengths)
    etime = (time_ns() - etime) / 1e9
    println("julia\t$precision\t$N\t$S\t$B\tsemifield\tcpu\t$etime")

    # Move the data and the FSM to the GPU
    fsm = fsm |> MarkovModels.gpu
    lhs = CuArray(lhs)

    # Dummy run to force compilation.
    if warmup pdfposteriors(fsm, lhs, lengths) end

    etime = time_ns()
    pdfposteriors(fsm, lhs, lengths)
    etime = (time_ns() - etime) / 1e9
    println("julia\t$precision\t$N\t$S\t$B\tsemifield\tgpu\t$etime")
end

main(
    "denominator_fsm.jld2", # "numerator_fsm.jld2"
    Float32,    # computing precision
    128,        # batch size
    700;        # sequence length
    warmup=true # if true, make a dummy run first to avoid compilation time
)

max_stars_repo_name: BUTSpeechFIT/HiddenMarkovModel

max_issues_repo_name: BUTSpeechFIT/MarkovChains

max_forks_repo_name: BUTSpeechFIT/ProbabilisticFSMs


text:
"""
Polygon boundaries are usually stored in a format called shapefile. For example, the shapefile of the world will consist of country names and each country will have coordinates of the boundaries.
```julia
load_shapefile("assets/mumbai_map/mumbai_districts.shp")
```
"""
function load_shapefile(filepath)
    table           = Shapefile.Table(filepath)                        ## Reading the shapefile
    shapefile_df    = DataFrame(table)                             ## Converting to dataframe
    geoms                        = Shapefile.shapes(table)         ## Reading shapes from the shapefile
    shapefile_df[!, "geometry"]  = geoms
    return shapefile_df
end

function make_polygon(geometry::CoordinateSystem, coords)
    array           = Array{Float16}[]          ## Declaring an empty array of Float32
    array           = coords[1]                 ## Accesssing n(th) polygon from the multipolgyon.
    poly    = []             ## Declaring a Luxor.Point type variable
    for x in array
        push!(poly, Point(long_to_column(geometry, x[1]), lat_to_row(geometry, x[2]))) # make FLOAT16/INT16
            # We convert the long,lat to 'Luxor.Point' and store them in points_luxor
    end
    return convert(Array{Point, 1}, poly)
end

function make_polygons(geometry::CoordinateSystem, geoms)
    #Makes a list of Luxor polygons from list of list of coordinates
    coord       = GeoInterface.coordinates(geoms)
    polygons    = Any[]
    for x in coord
        push!(polygons, make_polygon(geometry, x))
    end
    return polygons
end

"""
The polygon of a shapefile row can be make into a mask. This means all the points inside the polygon will be marked as 1, while the points outside will be marked as 0.

```julia
mumbai_districts = load_shapefile("assets/mumbai_map/mumbai_districts.shp")
district1 = mumbai_dists[1,:] # Select the first district
district1_mask = polygon_mask(my_coordinate_system, district1)
```
"""
function polygon_mask(geometry::CoordinateSystem, shapefile_row)
    geoms = shapefile_row.geometry
    polygons = make_polygons(geometry,geoms)
    points = zeros(Int32, geometry.height, geometry.width)
    for poly in polygons
        Boundary    = BoundingBox(poly)
        iend        = Int16(round(Boundary[2].y))
        j           = Int16(round(Boundary[1].x))
        jend        = Int16(round(Boundary[2].x))
        i           = Int16(round(Boundary[1].y))
        for e in i:iend
            for f in j:jend
                if isinside(Point(f, e), poly, allowonedge = true) == true
                    points[e,f] = 1
                end
            end
        end
    end
    return sparse(points)
end


text:
"""
   Define the partial derivative operator for v_z with PML boundary
"""
function Mvz(gammaz_stagger::Vector{Tv}, buoy::Matrix{Tv},
         dz::Tv, dt::Tv, fd::Vector{Tv}) where {Tv<:AbstractFloat}

    # size of the model
    (nz, nx) = size(buoy)
    if nz != length(gammaz_stagger)
       error("size mismatch")
    end

    # precision order
    order = length(fd)

    # data format
    data_format = eltype(buoy)
    double      = data_format(2.0)

    # the coefficients in front of v_z
    MvzBvz = ones(data_format, nz)
    for i = 1 : nz
        MvzBvz[i] = (double - dt * gammaz_stagger[i]) / (double + dt * gammaz_stagger[i])
    end

    # the coefficients in front of p
    MvzBp  = ones(data_format, nz*nx)
    idx = 0
    for j = 1 : nx
        for i = 1 : nz
            idx = idx + 1

            # averge of density
            if i < nz
               bavg = (buoy[i,j] + buoy[i+1,j]) / double
            else
               bavg = buoy[i,j]
            end

            MvzBp[idx] = (double * bavg * dt) / (double+dt*gammaz_stagger[i])
        end
    end

    # spatial partial derivative
    dpdz = zeros(data_format, nz, nz)
    for ir = 1 : nz
        lower = ir - order + 1 >= 1  ? ir-order+1 : 1
        upper = ir + order     <= nz ? ir+order   : nz

        for ic = lower : upper
            if ic-ir <= 0
               idx = ir-ic + 1
               dpdz[ir,ic] = - fd[idx] / dz
            elseif ic-ir >= 1
               idx = ic-ir
               dpdz[ir,ic] =   fd[idx] / dz
            end
        end
    end

    return MvzBvz, MvzBp, sparse(dpdz)
end

"""
   Define the partial derivative operator for v_z with rigid boundary
"""
function Rvz(buoy::Matrix{Tv}, dz::Tv, dt::Tv, fd::Vector{Tv}) where {Tv<:AbstractFloat}

    # size of the model
    (nz, nx) = size(buoy)

    # precision order
    order = length(fd)

    # data format
    data_format = eltype(buoy)
    double      = data_format(2.0)

    # the coefficients in front of p
    RvzBp  = ones(data_format, nz*nx)
    idx = 0
    for j = 1 : nx
        for i = 1 : nz
            idx = idx + 1

            # averge of density
            if i < nz
               bavg = (buoy[i,j] + buoy[i+1,j]) / double
            else
               bavg = buoy[i,j]
            end

            RvzBp[idx] = bavg * dt
        end
    end

    # spatial partial derivative
    rpdz = zeros(data_format, nz, nz)
    for ir = 1 : nz
        lower = ir - order + 1 >= 1  ? ir-order+1 : 1
        upper = ir + order     <= nz ? ir+order   : nz

        for ic = lower : upper
            if ic-ir <= 0
               idx = ir-ic + 1
               rpdz[ir,ic] = - fd[idx] / dz
            elseif ic-ir >= 1
               idx = ic-ir
               rpdz[ir,ic] =   fd[idx] / dz
            end
        end
    end

    return RvzBp, sparse(rpdz)
end

"""
   Define the partial derivative operator for v_x with PML boundary
"""
function Mvx(gammax_stagger::Vector{Tv}, buoy::Matrix{Tv},
         dx::Tv, dt::Tv, fd::Vector{Tv}) where {Tv<:AbstractFloat}

    # size of the model
    (nz, nx) = size(buoy)
    if nx != length(gammax_stagger)
       error("size mismatch")
    end

    # precision order
    order = length(fd)

    # data format
    data_format = eltype(buoy)
    double      = data_format(2.0)

    # the coefficients in front of v_x
    MvxBvx = ones(data_format, nx)
    for j = 1 : nx
        MvxBvx[j] = (double - dt * gammax_stagger[j]) / (double + dt * gammax_stagger[j])
    end

    # the coefficients in front of p
    MvxBp  = ones(data_format, nz*nx)
    idx = 0
    for j = 1 : nx
        for i = 1 : nz
            idx = idx + 1

            # averge of density
            if j < nx
               bavg = (buoy[i,j] + buoy[i,j+1]) / double
            else
               bavg = buoy[i,j]
            end

            MvxBp[idx] = (double * bavg * dt) / (double + dt * gammax_stagger[j])
        end
    end

    # spatial partial derivative dpdx
    dpdx = zeros(data_format, nx, nx)
    for ir = 1 : nx
        lower = ir - order + 1 >= 1  ? ir-order+1 : 1
        upper = ir + order     <= nx ? ir+order   : nx

        for ic = lower : upper
            if ic-ir <= 0
               idx = ir - ic + 1
               dpdx[ir,ic] = - fd[idx] / dx
            elseif ic-ir >= 1
               idx = ic - ir
               dpdx[ir,ic] =   fd[idx] / dx
            end
        end
    end

    return MvxBvx, MvxBp, sparse(dpdx)
end

"""
   Define the partial derivative operator for v_x with rigid boundary
"""
function Rvx(buoy::Matrix{Tv}, dx::Tv, dt::Tv, fd::Vector{Tv}) where {Tv<:AbstractFloat}

    # size of the model
    (nz, nx) = size(buoy)

    # precision order
    order = length(fd)

    # data format
    data_format = eltype(buoy)
    double      = data_format(2.0)

    # the coefficients in front of p
    RvxBp  = ones(data_format, nz*nx)
    idx = 0
    for j = 1 : nx
        for i = 1 : nz
            idx = idx + 1

            # averge of density
            if j < nx
               bavg = (buoy[i,j] + buoy[i,j+1]) / double
            else
               bavg = buoy[i,j]
            end

            RvxBp[idx] = bavg * dt
        end
    end

    # spatial partial derivative dpdx
    rpdx = zeros(data_format, nx, nx)
    for ir = 1 : nx
        lower = ir - order + 1 >= 1  ? ir-order+1 : 1
        upper = ir + order     <= nx ? ir+order   : nx

        for ic = lower : upper
            if ic-ir <= 0
               idx = ir - ic + 1
               rpdx[ir,ic] = - fd[idx] / dx
            elseif ic-ir >= 1
               idx = ic - ir
               rpdx[ir,ic] =   fd[idx] / dx
            end
        end
    end

    return RvxBp, sparse(rpdx)
end

"""
   Define the partial derivative operator for p_z with PML boundary
"""
function Mpz(gammaz_integer::Vector{Tv}, bulk::Matrix{Tv},
         dz::Tv, dt::Tv, fd::Vector) where {Tv<:AbstractFloat}

    # size of the model
    (nz, nx) = size(bulk)
    if nz != length(gammaz_integer)
       error("size mismatch")
    end

    # precision order
    order = length(fd)

    # data format
    data_format = eltype(bulk)
    double      = data_format(2.0)

    # the coefficients in front of p_z
    MpzBpz = ones(data_format, nz)
    for i = 1 : nz
        MpzBpz[i] = (double - dt * gammaz_integer[i]) / (double + dt * gammaz_integer[i])
    end

    # the coefficients in front of v_z
    MpzBvz = ones(data_format, nz*nx)
    idx = 0
    for j = 1 : nx
        for i = 1 : nz
            idx = idx + 1
            MpzBvz[idx] = (double * bulk[i,j] * dt) / (double + dt * gammaz_integer[i])
        end
    end

    # spatial partial derivative dvdz
    dvdz = zeros(data_format, nz, nz)
    for ir = 1 : nz
        lower = ir - order     >= 1  ? ir-order   : 1
        upper = ir + order - 1 <= nz ? ir+order-1 : nz

        for ic = lower : upper
            if ic-ir <= -1
               idx = ir-ic
               dvdz[ir,ic] = - fd[idx] / dz
            elseif ic-ir >= 0
               idx = ic-ir+1
               dvdz[ir,ic] =   fd[idx] / dz
            end
        end
    end

    return MpzBpz, MpzBvz, sparse(dvdz)
end

"""
   Define the partial derivative operator for p_x with PML boundary
"""
function Mpx(gammax_integer::Vector{Tv}, bulk::Matrix{Tv},
         dx::Tv, dt::Tv, fd::Vector{Tv}) where {Tv<:AbstractFloat}

    # size of the model
    (nz, nx) = size(bulk)
    if nx != length(gammax_integer)
       error("size mismatch")
    end

    # precision order
    order = length(fd)

    # data format
    data_format = eltype(bulk)
    double      = data_format(2.0)

    # the coefficients in front of p_x
    MpxBpx = ones(data_format, nx)
    for j = 1 : nx
        MpxBpx[j] = (double - dt * gammax_integer[j]) / (double + dt * gammax_integer[j])
    end

    # the coefficients in front of v_x
    MpxBvx = ones(data_format, nz*nx)
    idx = 0
    for j = 1 : nx
        for i = 1 : nz
            idx = idx + 1
            MpxBvx[idx] = (double * bulk[i,j] * dt) / (double + dt * gammax_integer[j])
        end
    end

    # spatial partial derivative dvdx
    dvdx = zeros(data_format, nx, nx)
    for ir = 1 : nx
        lower = ir - order     >= 1  ? ir-order   : 1
        upper = ir + order - 1 <= nx ? ir+order-1 : nx

        for ic = lower : upper
            if ic-ir <= -1
               idx = ir - ic
               dvdx[ir,ic] = - fd[idx] / dx
            elseif ic-ir >= 0
               idx = ic - ir + 1
               dvdx[ir,ic] =   fd[idx] / dx
            end
        end
    end

    return MpxBpx, MpxBvx, sparse(dvdx)
end

"""
   Define the partial derivative operator for p_z with rigid boundary
"""
function Rpzx(bulk::Matrix{Tv}, dz::Tv, dx::Tv, dt::Tv, fd::Vector) where {Tv<:AbstractFloat}

    # size of the model
    (nz, nx) = size(bulk)

    # precision order
    order = length(fd)

    # data format
    data_format = eltype(bulk)
    double      = data_format(2.0)

    # the coefficients in front of v_z
    RpBv = ones(data_format, nz*nx)
    idx = 0
    for j = 1 : nx
        for i = 1 : nz
            idx = idx + 1
            RpBv[idx] = bulk[i,j] * dt
        end
    end

    # spatial partial derivative dvdz
    rvdz = zeros(data_format, nz, nz)
    for ir = 1 : nz
        lower = ir - order     >= 1  ? ir-order   : 1
        upper = ir + order - 1 <= nz ? ir+order-1 : nz

        for ic = lower : upper
            if ic-ir <= -1
               idx = ir-ic
               rvdz[ir,ic] = - fd[idx] / dz
            elseif ic-ir >= 0
               idx = ic-ir+1
               rvdz[ir,ic] =   fd[idx] / dz
            end
        end
    end

    # spatial partial derivative dvdx
    rvdx = zeros(data_format, nx, nx)
    for ir = 1 : nx
        lower = ir - order     >= 1  ? ir-order   : 1
        upper = ir + order - 1 <= nx ? ir+order-1 : nx

        for ic = lower : upper
            if ic-ir <= -1
               idx = ir - ic
               rvdx[ir,ic] = - fd[idx] / dx
            elseif ic-ir >= 0
               idx = ic - ir + 1
               rvdx[ir,ic] =   fd[idx] / dx
            end
        end
    end

    return RpBv, sparse(rvdz), sparse(rvdx)
end

# """
#    sparse matrix for finite-difference grid with rigid boundary
# """
# struct RigidFDStencil{Ti<:Int64, Tv<:AbstractFloat}
#     MvzBp  :: SparseMatrixCSC{Tv,Ti}
#     MvxBp  :: SparseMatrixCSC{Tv,Ti}
#     MpzBvz :: SparseMatrixCSC{Tv,Ti}
#     MpxBvx :: SparseMatrixCSC{Tv,Ti}
# end
#
# """
#    constructor for finite difference stencil with rigid boundary condition
# """
# function RigidFDStencil(params::TdParams)
#
#     # expand physical model parameters
#     lambda = vp2lambda(params.vel, params.rho)
#
#     MvzBp  = Mvz(params.rho, params.dz, params.dt, params.fd_coefficients)
#     MvxBp  = Mvx(params.rho, params.dx, params.dt, params.fd_coefficients)
#     MpzBvz = Mpz(lambda    , params.dz, params.dt, params.fd_coefficients)
#     MpxBvx = Mpx(lambda    , params.dx, params.dt, params.fd_coefficients)
#
#     return RigidFDStencil(MvzBp, MvxBp, MpzBvz, MpxBvx)
# end

# """
#    sparse matrix for finite-difference grid with PML boundary
# """
# struct ObsorbFDStencil{Ti<:Int64, Tv<:AbstractFloat}
#     MvzBvz :: Vector{Tv}
#     MvzBp  :: SparseMatrixCSC{Tv,Ti}
#     MvxBvx :: Vector{Tv}
#     MvxBp  :: SparseMatrixCSC{Tv,Ti}
#     MpzBpz :: Vector{Tv}
#     MpzBvz :: SparseMatrixCSC{Tv,Ti}
#     MpxBpx :: Vector{Tv}
#     MpxBvx :: SparseMatrixCSC{Tv,Ti}
# end
#
# """
#    sparse matrix for finite-difference grid with rigid boundary
# """
# struct RigidFDStencil{Ti<:Int64, Tv<:AbstractFloat}
#     MvzBp  :: SparseMatrixCSC{Tv,Ti}
#     MvxBp  :: SparseMatrixCSC{Tv,Ti}
#     MpzBvz :: SparseMatrixCSC{Tv,Ti}
#     MpxBvx :: SparseMatrixCSC{Tv,Ti}
# end
#
# """
#    compute lambda from velocity and density
# """
# function vp2lambda(vp::Matrix{Tv}, rho::Matrix{Tv}) where {Tv<:AbstractFloat}
#
#     (m , n ) = size(vp)
#     (m1, n1) = size(rho)
#     if m1 != m || n1 != n
#        error("size mismatch")
#     end
#
#     lambda = zeros(eltype(vp), m, n)
#     for j = 1 : n
#         for i = 1 : m
#             lambda[i,j] = rho[i,j] * vp[i,j]^2
#         end
#     end
#     return lambda
# end
#
# """
#    Define the partial derivative operator for v_z with PML boundary
# """
# function Mvz(vz_pml::SparseMatrixCSC{Tv,Ti}, rho::Matrix{Tv},
#          dz::Tv, dt::Tv, fd::Vector{Tv}) where {Ti<:Int64, Tv<:AbstractFloat}
#
#     # data format
#     data_format = eltype(rho)
#     double_dt   = data_format(2.0 * dt)
#
#     # finite-difference coefficients
#     order = length(fd)
#     c     = fd / dz
#
#     # model size
#     (m, n) = size(vz_pml)
#     c_diag = zeros(data_format, m*n)
#     denum  = zeros(data_format, m*n)
#     for ix = 1 : n
#         for iz = 1 : m
#
#             if iz < m
#                a = (rho[iz+1,ix]+rho[iz,ix]) / double_dt
#             else
#                a = rho[iz,ix] / dt
#             end
#             b = vz_pml[iz,ix] / 2
#             denum[(ix-1)*m+iz] = 1.0   / (a+b)
#             c_diag[(ix-1)*m+iz]= (a-b) / (a+b)
#         end
#     end
#
#     # spatial partial derivative
#     tmp = zeros(data_format, m, m)
#     for iz = 1 : m
#         lower = iz - order + 1 >= 1 ? iz-order+1 : 1
#         upper = iz + order     <= m ? iz+order   : m
#
#         for ix = lower : upper
#             if ix-iz <= 0
#                idx = iz-ix + 1
#                tmp[iz,ix] = - c[idx]
#             elseif ix-iz >= 1
#                idx = ix-iz
#                tmp[iz,ix] =   c[idx]
#             end
#         end
#     end
#
#     MvzBvz = c_diag
#     MvzBp  = spdiagm(0 => denum) * kron(spdiagm(0=>ones(data_format,n)), sparse(tmp))
#
#     return MvzBvz, MvzBp
#
# end
#
# """
#    Define the partial derivative operator for v_z with rigid boundary
# """
# function Mvz(rho::Matrix{Tv}, dz::Tv, dt::Tv, fd::Vector{Tv}) where {Tv<:AbstractFloat}
#
#     # data format
#     data_format = eltype(rho)
#     double_dt   = data_format(2.0 * dt)
#
#     # finite-difference coefficients
#     order = length(fd)
#     c     = fd / dz
#
#     # model size
#     (m, n) = size(rho)
#     denum  = zeros(data_format, m*n)
#     for ix = 1 : n
#         for iz = 1 : m
#
#             if iz < m
#                a = (rho[iz+1,ix]+rho[iz,ix]) / double_dt
#             else
#                a = rho[iz,ix] / dt
#             end
#             denum[(ix-1)*m+iz] = 1.0 / a
#         end
#     end
#
#     # spatial partial derivative
#     tmp = zeros(data_format, m, m)
#     for iz = 1 : m
#         lower = iz - order + 1 >= 1 ? iz-order+1 : 1
#         upper = iz + order     <= m ? iz+order   : m
#
#         for ix = lower : upper
#             if ix-iz <= 0
#                idx = iz-ix + 1
#                tmp[iz,ix] = - c[idx]
#             elseif ix-iz >= 1
#                idx = ix-iz
#                tmp[iz,ix] =   c[idx]
#             end
#         end
#     end
#
#     MvzBp  = spdiagm(0 => denum) * kron(spdiagm(0=>ones(data_format,n)), sparse(tmp))
#     return MvzBp
# end
#
# """
#    Define the partial derivative operator for v_x with PML boundary
# """
# function Mvx(vx_pml::SparseMatrixCSC{Tv,Ti}, rho::Matrix{Tv},
#          dx::Tv, dt::Tv, fd::Vector{Tv}) where {Ti<:Int64, Tv<:AbstractFloat}
#
#     # data format
#     data_format  = eltype(rho)
#     double_dt    = data_format(2.0 * dt)
#
#     # finite-difference coefficients
#     order = length(fd)
#     c     = fd / dx
#
#     # model size
#     (m, n) = size(vx_pml)
#     c_diag = zeros(data_format, m*n)
#     denum  = zeros(data_format, m*n)
#     for ix = 1 : n
#         for iz = 1 : m
#
#             if ix < n
#                a = (rho[iz,ix+1]+rho[iz,ix]) / double_dt
#             else
#                a = rho[iz,ix] / dt
#             end
#             b = vx_pml[iz,ix] / 2
#             denum[(ix-1)*m+iz] = 1     / (a+b)
#             c_diag[(ix-1)*m+iz]= (a-b) / (a+b)
#         end
#     end
#
#     # spatial partial derivative
#     tmp = zeros(data_format, n, n)
#     for iz = 1 : n
#         lower = iz - order + 1 >= 1 ? iz-order+1 : 1
#         upper = iz + order     <= n ? iz+order   : n
#
#         for ix = lower : upper
#             if ix-iz <= 0
#                idx = iz-ix + 1
#                tmp[iz,ix] = - c[idx]
#             elseif ix-iz >= 1
#                idx = ix-iz
#                tmp[iz,ix] =   c[idx]
#             end
#         end
#     end
#
#     MvxBvx = c_diag
#     MvxBp  = spdiagm(0 => denum) * kron(sparse(tmp), spdiagm(0=>ones(data_format,m)))
#
#     return MvxBvx, MvxBp
#
# end
#
# """
#    Define the partial derivative operator for v_x with rigid boundary
# """
# function Mvx(rho::Matrix{Tv}, dx::Tv, dt::Tv, fd::Vector{Tv}) where {Tv<:AbstractFloat}
#
#     # data format
#     data_format  = eltype(rho)
#     double_dt = data_format(2.0 * dt)
#
#     # finite-difference coefficients
#     order = length(fd)
#     c     = fd / dx
#
#     # model size
#     (m, n) = size(rho)
#     denum  = zeros(data_format, m*n)
#     for ix = 1 : n
#         for iz = 1 : m
#             # average density
#             if ix < n
#                a = (rho[iz,ix+1]+rho[iz,ix]) / double_dt
#             else
#                a = rho[iz,ix] / dt
#             end
#             denum[(ix-1)*m+iz] = 1.0 / a
#         end
#     end
#
#     # spatial partial derivative
#     tmp = zeros(data_format, n, n)
#     for iz = 1 : n
#         lower = iz - order + 1 >= 1 ? iz-order+1 : 1
#         upper = iz + order     <= n ? iz+order   : n
#
#         for ix = lower : upper
#             if ix-iz <= 0
#                idx = iz-ix + 1
#                tmp[iz,ix] = - c[idx]
#             elseif ix-iz >= 1
#                idx = ix-iz
#                tmp[iz,ix] =   c[idx]
#             end
#         end
#     end
#
#     MvxBp  = spdiagm(0 => denum) * kron(sparse(tmp), spdiagm(0=>ones(data_format,m)))
#     return MvxBp
# end
#
# """
#    Define the partial derivative operator for p_z with PML boundary
# """
# function Mpz(pz_pml::SparseMatrixCSC{Tv,Ti}, lambda::Matrix{Tv},
#          dz::Tv, dt::Tv, fd::Vector{Tv}) where {Ti<:Int64, Tv<:AbstractFloat}
#
#     # data format
#     data_format = eltype(lambda)
#
#     # finite-difference coefficients
#     order = length(fd)
#     c     = fd / dz
#
#     # model size
#     (m, n) = size(pz_pml)
#     c_diag = zeros(data_format, m*n)
#     denum  = zeros(data_format, m*n)
#     for ix = 1 : n
#         for iz = 1 : m
#             a = 1/dt + pz_pml[iz,ix]/2
#             b = 1/dt - pz_pml[iz,ix]/2
#             denum[(ix-1)*m+iz] = lambda[iz,ix] / a
#             c_diag[(ix-1)*m+iz]= b / a
#         end
#     end
#
#     # spatial partial derivative
#     tmp = zeros(data_format, m, m)
#     for iz = 1 : m
#         lower = iz - order     >= 1 ? iz-order   : 1
#         upper = iz + order - 1 <= m ? iz+order-1 : m
#
#         for ix = lower : upper
#             if ix-iz <= -1
#                idx = iz-ix
#                tmp[iz,ix] = - c[idx]
#             elseif ix-iz >= 0
#                idx = ix-iz+1
#                tmp[iz,ix] =   c[idx]
#             end
#         end
#     end
#
#     MpzBpz = c_diag
#     MpzBvz = spdiagm(0 => denum) * kron(spdiagm(0=>ones(data_format,n)), sparse(tmp))
#
#     return MpzBpz, MpzBvz
# end
#
# """
#    Define the partial derivative operator for p_z with rigid boundary
# """
# function Mpz(lambda::Matrix{Tv}, dz::Tv, dt::Tv, fd::Vector{Tv}) where {Tv<:AbstractFloat}
#
#     # data format
#     data_format = eltype(lambda)
#     c           = fd / dz
#
#     # model size
#     (m, n) = size(lambda)
#     denum  = zeros(data_format, m*n)
#     for ix = 1 : n
#         for iz = 1 : m
#             denum[(ix-1)*m+iz] = lambda[iz,ix] * dt
#         end
#     end
#
#     # spatial partial derivative
#     order = length(fd)
#     tmp   = zeros(data_format, m, m)
#     for iz = 1 : m
#         lower = iz - order     >= 1 ? iz-order   : 1
#         upper = iz + order - 1 <= m ? iz+order-1 : m
#
#         for ix = lower : upper
#             if ix-iz <= -1
#                idx = iz-ix
#                tmp[iz,ix] = - c[idx]
#             elseif ix-iz >= 0
#                idx = ix-iz+1
#                tmp[iz,ix] =   c[idx]
#             end
#         end
#     end
#
#     MpzBvz = spdiagm(0 => denum) * kron(spdiagm(0=>ones(data_format,n)), sparse(tmp))
#     return MpzBvz
#
# end
#
# """
#    Define the partial derivative operator for p_x with PML boundary
# """
# function Mpx(px_pml::SparseMatrixCSC{Tv,Ti}, lambda::Matrix{Tv},
#          dx::Tv, dt::Tv, fd::Vector{Tv}) where {Tv<:AbstractFloat, Ti<:Int64}
#
#     # data format
#     data_format = eltype(lambda)
#
#     # finite-difference coefficients
#     order = length(fd)
#     c     = fd / dx
#
#     (m, n) = size(px_pml)
#     c_diag = zeros(data_format, m*n)
#     denum  = zeros(data_format, m*n)
#     for ix = 1 : n
#         for iz = 1 : m
#             a = 1/dt + px_pml[iz,ix]/2
#             b = 1/dt - px_pml[iz,ix]/2
#             denum[(ix-1)*m+iz] = lambda[iz,ix] / a
#             c_diag[(ix-1)*m+iz]= b / a
#         end
#     end
#
#     # spatial partial derivative
#     tmp = zeros(data_format, n, n)
#     for iz = 1 : n
#         lower = iz - order     >= 1 ? iz-order   : 1
#         upper = iz + order - 1 <= n ? iz+order-1 : n
#
#         for ix = lower : upper
#             if ix-iz <= -1
#                idx = iz-ix
#                tmp[iz,ix] = - c[idx]
#             elseif ix-iz >= 0
#                idx = ix-iz+1
#                tmp[iz,ix] =   c[idx]
#             end
#         end
#     end
#
#     MpxBpx = c_diag
#     MpxBvx = spdiagm(0 => denum) * kron(sparse(tmp), spdiagm(0=>ones(data_format,m)))
#
#     return  MpxBpx, MpxBvx
# end
#
# """
#    Define the partial derivative operator for p_x with rigid boundary
# """
# function Mpx(lambda::Matrix{Tv}, dx::Tv, dt::Tv, fd::Vector{Tv}) where {Tv<:AbstractFloat}
#
#     # data format
#     data_format = eltype(lambda)
#
#     # finite-difference coefficients
#     order = length(fd)
#     c     = fd / dx
#
#     (m, n) = size(lambda)
#     denum  = zeros(data_format, m*n)
#     for ix = 1 : n
#         for iz = 1 : m
#             denum[(ix-1)*m+iz] = lambda[iz,ix] * dt
#         end
#     end
#
#     # spatial partial derivative
#     tmp = zeros(data_format, n, n)
#     for iz = 1 : n
#         lower = iz - order     >= 1 ? iz-order   : 1
#         upper = iz + order - 1 <= n ? iz+order-1 : n
#
#         for ix = lower : upper
#             if ix-iz <= -1
#                idx = iz-ix
#                tmp[iz,ix] = - c[idx]
#             elseif ix-iz >= 0
#                idx = ix-iz+1
#                tmp[iz,ix] =   c[idx]
#             end
#         end
#     end
#
#     MpxBvx = spdiagm(0 => denum) * kron(sparse(tmp), spdiagm(0=>ones(data_format,m)))
#
#     return MpxBvx
# end
#
# """
#    constructor for obsorbing finite difference stencil
# """
# function ObsorbFDStencil(params::TdParams)
#
#     # pml coefficient
#     pml = PMLCoefficients(params)
#
#     # padding physical model parameters
#     lambda = vp2lambda(params.vel, params.rho)
#     lambda = model_padding(lambda    , params.npml, params.free_surface)
#     rho    = model_padding(params.rho, params.npml, params.free_surface)
#
#     (MvzBvz, MvzBp ) = Mvz(pml.vz, rho   , params.dz, params.dt, params.fd_coefficients)
#     (MvxBvx, MvxBp ) = Mvx(pml.vx, rho   , params.dx, params.dt, params.fd_coefficients)
#     (MpzBpz, MpzBvz) = Mpz(pml.pz, lambda, params.dz, params.dt, params.fd_coefficients)
#     (MpxBpx, MpxBvx) = Mpx(pml.px, lambda, params.dx, params.dt, params.fd_coefficients)
#
#     return ObsorbFDStencil(MvzBvz, MvzBp , MvxBvx, MvxBp,
#                            MpzBpz, MpzBvz, MpxBpx, MpxBvx)
# end
#
# """
#    constructor for finite difference stencil with rigid boundary condition
# """
# function RigidFDStencil(params::TdParams)
#
#     # expand physical model parameters
#     lambda = vp2lambda(params.vel, params.rho)
#
#     MvzBp  = Mvz(params.rho, params.dz, params.dt, params.fd_coefficients)
#     MvxBp  = Mvx(params.rho, params.dx, params.dt, params.fd_coefficients)
#     MpzBvz = Mpz(lambda    , params.dz, params.dt, params.fd_coefficients)
#     MpxBvx = Mpx(lambda    , params.dx, params.dt, params.fd_coefficients)
#
#     return RigidFDStencil(MvzBp, MvxBp, MpzBvz, MpxBvx)
# end

# old version
# function compute_relative_error(fd::Vector{Tv}) where {Tv<:AbstractFloat}
#
#     M = length(fd)
#
#     interval = 0.001 * pi
#     beta = collect(interval : interval : pi)
#     N    = length(beta)
#     err  = zeros(N)
#
#     for i = 1 : N
#         tmp = 0.0
#         for m = 1 : M
#             tmp = tmp + fd[m] * sin((m-0.5)*beta[i])
#         end
#         err[i] = 2.0 * tmp / beta[i] - 1.0
#     end
#
#     return beta, err
# end

# the coefficients obtained by minimize absolute error
# fdc = Vector{Vector{Float64}}(9)
# fdc[1] = [0.1129136e1, -0.4304542e-1]
# fdc[2] = [0.1186247e1, -0.7266808e-1, 0.6351497e-2]
# fdc[3] = [0.1218159e1, -0.9397218e-1, 0.1519043e-1, -0.1742128e-2]
# fdc[4] = [0.1236425e1, -0.1081130   , 0.2339911e-1, -0.5061550e-2, 0.7054313e-3]
# fdc[5] = [0.1247576e1, -0.1174969   , 0.2997288e-1, -0.8741572e-2, 0.2262285e-2, -0.3745306e-3]
# fdc[6] = [0.1254380e1, -0.1235307   , 0.3467231e-1, -0.1192915e-1, 0.4057090e-2, -0.1191005e-2, 0.2263204e-3]
# fdc[7] = [0.1259012e1, -0.1277647   , 0.3820715e-1, -0.1458251e-1, 0.5845385e-2, -0.2213861e-2, 0.7243880e-3, -0.1566173e-3]
# fdc[8] = [0.1262147e1, -0.1306967   , 0.4075792e-1, -0.1665221e-1, 0.7377057e-2, -0.3258150e-2, 0.1336259e-2, -0.4775830e-3, 0.1151664e-3]
# fdc[9] = [0.1264362e1, -0.1327958   , 0.4264687e-1, -0.1824918e-1, 0.8656223e-2, -0.4200034e-2, 0.1989180e-2, -0.8686637e-3, 0.3342741e-3, -0.8854090e-4]

# M = 5
#
# (x_axis, err ) = compute_relative_error(fdc[M-1])
# (x_axis, err1) = compute_relative_error(rfdc[M-1])
#
# cfdc = finite_difference_coefficient(M)
# (x_axis, err2) = compute_relative_error(cfdc)
#
# plot(x_axis, err, label="new_a")
# plot(x_axis, err1, label="new_r")
# plot(x_axis, err2, label="old"); legend();

# """
#    finite difference stencil for updating vz
# """
# function Mvz(vz_pml::SparseMatrixCSC{Tv,Ti}, rho::Matrix{Tv},
#          dz::Tv, dt::Tv) where {Ti<:Int64, Tv<:AbstractFloat}
#
#     data_format  = eltype(rho)
#     (m, n) = size(vz_pml)
#     double_dt = data_format(2.0 * dt)
#
#     a1 = data_format( 9 / 8 )
#     a2 = data_format(-1 / 24)
#     c1 = a1 / dz
#     c2 = a2 / dz
#
#     c3    = zeros(data_format, m*n)
#     denum = zeros(data_format, m*n)
#
#     for ix = 1 : n
#         for iz = 1 : m
#             if iz < m
#                a = (rho[iz+1,ix]+rho[iz,ix]) / double_dt
#             else
#                a = rho[iz,ix] / dt
#             end
#             b = vz_pml[iz,ix] / 2
#
#             denum[(ix-1)*m+iz] = 1.0   / (a+b)
#             c3[(ix-1)*m+iz]    = (a-b) / (a+b)
#         end
#     end
#
#     tmp = zeros(data_format, m, m)
#     tmp[1,1] = -c1; tmp[1,2] = c1; tmp[1,3] = c2;
#
#     for iz = 2: m-2
#         tmp[iz,iz  ] = -c1; tmp[iz, iz+1] = c1;
#         tmp[iz,iz-1] = -c2; tmp[iz, iz+2] = c2;
#     end
#
#     tmp[m-1,m-1] = -c1; tmp[m-1, m  ] = c1; tmp[m-1, m-2] = -c2;
#     tmp[m  ,m  ] = -c1; tmp[m  , m-1] =-c2;
#
#     MvzBvz = c3
#     MvzBp  = spdiagm(denum) * kron(speye(data_format, n), sparse(tmp))
#
#     return MvzBvz, MvzBp
#
# end
#
# """
#    rigid boundary
# """
# function Mvz(rho::Matrix{Tv}, dz::Tv, dt::Tv) where {Tv<:AbstractFloat}
#
#     data_format  = eltype(rho)
#     (m, n) = size(rho)
#     double_dt = data_format(2.0 * dt)
#
#     a1 = data_format( 9 / 8 )
#     a2 = data_format(-1 / 24)
#     c1 = a1 / dz
#     c2 = a2 / dz
#
#     denum = zeros(data_format, m*n)
#
#     for ix = 1 : n
#         for iz = 1 : m
#             if iz < m
#                a = (rho[iz+1,ix]+rho[iz,ix]) / double_dt
#             else
#                a = rho[iz,ix] / dt
#             end
#             denum[(ix-1)*m+iz] = 1 / a
#         end
#     end
#
#     tmp = zeros(data_format, m, m)
#     tmp[1,1] = -c1; tmp[1,2] = c1; tmp[1,3] = c2;
#
#     for iz = 2: m-2
#         tmp[iz,iz  ] = -c1; tmp[iz, iz+1] = c1;
#         tmp[iz,iz-1] = -c2; tmp[iz, iz+2] = c2;
#     end
#
#     tmp[m-1,m-1] = -c1; tmp[m-1, m  ] = c1; tmp[m-1, m-2] = -c2;
#     tmp[m  ,m  ] = -c1; tmp[m  , m-1] =-c2;
#
#     MvzBp  = spdiagm(denum) * kron(speye(data_format, n), sparse(tmp))
#
#     return MvzBp
#
# end

# """
#    finite difference stencil for updating vx
# """
# function Mvx(vx_pml::SparseMatrixCSC{Tv,Ti}, rho::Matrix{Tv},
#          dx::Tv, dt::Tv) where {Ti<:Int64, Tv<:AbstractFloat}
#
#     data_format  = eltype(rho)
#     (m, n) = size(vx_pml)
#     double_dt = data_format(2.0 * dt)
#
#     a1 = data_format(9/8  )
#     a2 = data_format(-1/24)
#     c1 = a1 / dx
#     c2 = a2 / dx
#
#     c3    = zeros(data_format, m*n)
#     denum = zeros(data_format, m*n)
#
#     for ix = 1 : n
#         for iz = 1 : m
#             if ix < n
#                a = (rho[iz,ix+1]+rho[iz,ix]) / double_dt
#             else
#                a = rho[iz,ix] / dt
#             end
#             b = vx_pml[iz,ix] / 2
#
#             denum[(ix-1)*m+iz] = 1     / (a+b)
#             c3[(ix-1)*m+iz]    = (a-b) / (a+b)
#         end
#     end
#
#     tmp = zeros(data_format, n, n)
#     tmp[1,1] = -c1; tmp[1,2] = c1; tmp[1,3] = c2;
#
#     for ix = 2: n-2
#         tmp[ix,ix  ] = -c1; tmp[ix, ix+1] = c1;
#         tmp[ix,ix-1] = -c2; tmp[ix, ix+2] = c2;
#     end
#
#     tmp[n-1,n-1] = -c1; tmp[n-1, n  ] = c1; tmp[n-1,n-2] = -c2;
#     tmp[n  ,n  ] = -c1; tmp[n  , n-1] =-c2;
#
#     MvxBvx = c3
#     MvxBp = spdiagm(denum) * kron(sparse(tmp), speye(data_format, m))
#
#     return MvxBvx, MvxBp
#
# end
#
# """
#    rigid boundary Mvx
# """
# function Mvx(rho::Matrix{Tv}, dx::Tv, dt::Tv) where {Tv<:AbstractFloat}
#
#     data_format  = eltype(rho)
#     (m, n) = size(rho)
#     double_dt = data_format(2.0 * dt)
#
#     a1 = data_format(9/8  )
#     a2 = data_format(-1/24)
#     c1 = a1 / dx
#     c2 = a2 / dx
#     denum = zeros(data_format, m*n)
#
#     for ix = 1 : n
#         for iz = 1 : m
#             if ix < n
#                a = (rho[iz,ix+1]+rho[iz,ix]) / double_dt
#             else
#                a = rho[iz,ix] / dt
#             end
#             denum[(ix-1)*m+iz] = 1 / a
#         end
#     end
#
#     tmp = zeros(data_format, n, n)
#     tmp[1,1] = -c1; tmp[1,2] = c1; tmp[1,3] = c2;
#
#     for ix = 2: n-2
#         tmp[ix,ix  ] = -c1; tmp[ix, ix+1] = c1;
#         tmp[ix,ix-1] = -c2; tmp[ix, ix+2] = c2;
#     end
#
#     tmp[n-1,n-1] = -c1; tmp[n-1, n  ] = c1; tmp[n-1,n-2] = -c2;
#     tmp[n  ,n  ] = -c1; tmp[n  , n-1] =-c2;
#
#     MvxBp = spdiagm(denum) * kron(sparse(tmp), speye(data_format, m))
#
#     return MvxBp
#
# end

# """
#    finite difference stencil for updating pz
# """
# function Mpz(pz_pml::SparseMatrixCSC{Tv,Ti}, lambda::Matrix{Tv},
#          dz::Tv, dt::Tv) where {Ti<:Int64, Tv<:AbstractFloat}
#
#     data_format  = eltype(lambda)
#     (m, n) = size(pz_pml)
#
#     a1 = data_format( 9/8 )
#     a2 = data_format(-1/24)
#     c1 = a1 / dz
#     c2 = a2 / dz
#
#     c3    = zeros(data_format, m*n)
#     denum = zeros(data_format, m*n)
#
#     for ix = 1 : n
#         for iz = 1 : m
#             a = 1/dt + pz_pml[iz,ix]/2
#             b = 1/dt - pz_pml[iz,ix]/2
#             denum[(ix-1)*m+iz] = lambda[iz,ix] / a
#             c3[(ix-1)*m+iz]= b / a
#         end
#     end
#
#     tmp = zeros(data_format, m, m)
#     tmp[1,1] =  c1; tmp[1,2] = c2;
#     tmp[2,1] = -c1; tmp[2,2] = c1; tmp[2,3] = c2;
#
#     for iz = 3: m-1
#         tmp[iz,iz-1] = -c1; tmp[iz, iz  ] = c1;
#         tmp[iz,iz-2] = -c2; tmp[iz, iz+1] = c2;
#     end
#
#     tmp[m,m-1] = -c1; tmp[m,m] = c1; tmp[m, m-2] = -c2
#
#     MpzBpz = c3
#     MpzBvz = spdiagm(denum) * kron(speye(data_format, n), sparse(tmp))
#
#     return MpzBpz, MpzBvz
#
# end

# """
#    rigid boundary Mpz with high-order finite difference approximation
# """
# function Mpz(lambda::Matrix{Tv}, dz::Tv, dt::Tv) where {Tv<:AbstractFloat}
#
#     data_format  = eltype(lambda)
#     (m, n) = size(lambda)
#
#     a1 = data_format( 9/8 )
#     a2 = data_format(-1/24)
#     c1 = a1 / dz
#     c2 = a2 / dz
#
#     c3    = zeros(data_format, m*n)
#     denum = zeros(data_format, m*n)
#
#     for ix = 1 : n
#         for iz = 1 : m
#             denum[(ix-1)*m+iz] = lambda[iz,ix] * dt
#         end
#     end
#
#     tmp = zeros(data_format, m, m)
#     tmp[1,1] =  c1; tmp[1,2] = c2;
#     tmp[2,1] = -c1; tmp[2,2] = c1; tmp[2,3] = c2;
#
#     for iz = 3: m-1
#         tmp[iz,iz-1] = -c1; tmp[iz, iz  ] = c1;
#         tmp[iz,iz-2] = -c2; tmp[iz, iz+1] = c2;
#     end
#
#     tmp[m,m-1] = -c1; tmp[m,m] = c1; tmp[m, m-2] = -c2
#
#     MpzBvz = spdiagm(denum) * kron(speye(data_format, n), sparse(tmp))
#
#     return MpzBvz
# end

# """
#    finite difference stencil for updating px
# """
# function Mpx(px_pml::SparseMatrixCSC{Tv,Ti}, lambda::Matrix{Tv},
#          dx::Tv, dt::Tv) where {Tv<:AbstractFloat, Ti<:Int64}
#
#     data_format  = eltype(lambda)
#     (m, n) = size(px_pml)
#
#     a1 = data_format( 9/8 )
#     a2 = data_format(-1/24)
#     c1 = a1 / dx
#     c2 = a2 / dx
#
#     c3    = zeros(data_format, m*n)
#     denum = zeros(data_format, m*n)
#
#     for ix = 1 : n
#         for iz = 1 : m
#             a = 1/dt + px_pml[iz,ix]/2
#             b = 1/dt - px_pml[iz,ix]/2
#             denum[(ix-1)*m+iz] = lambda[iz,ix] / a
#             c3[(ix-1)*m+iz]= b / a
#         end
#     end
#
#     tmp = zeros(data_format, n, n)
#     tmp[1,1] =  c1; tmp[1,2] = c2;
#     tmp[2,1] = -c1; tmp[2,2] = c1; tmp[2,3] = c2;
#
#     for ix = 3: n-1
#         tmp[ix,ix-1] = -c1; tmp[ix, ix  ] = c1;
#         tmp[ix,ix-2] = -c2; tmp[ix, ix+1] = c2;
#     end
#
#     tmp[n,n-1] = -c1; tmp[n,n] = c1; tmp[n,n-2] = -c2;
#
#     MpxBpx = c3
#     MpxBvx = spdiagm(denum) * kron(sparse(tmp), speye(data_format, m))
#
#     return  MpxBpx, MpxBvx
# end
#
#
# """
#    rigid boundary Mpx
# """
# function Mpx(lambda::Matrix{Tv}, dx::Tv, dt::Tv) where {Tv<:AbstractFloat}
#
#     data_format  = eltype(lambda)
#     (m, n) = size(lambda)
#
#     a1 = data_format( 9/8 )
#     a2 = data_format(-1/24)
#     c1 = a1 / dx
#     c2 = a2 / dx
#
#     denum = zeros(data_format, m*n)
#
#     for ix = 1 : n
#         for iz = 1 : m
#             denum[(ix-1)*m+iz] = lambda[iz,ix] * dt
#         end
#     end
#
#     tmp = zeros(data_format, n, n)
#     tmp[1,1] =  c1; tmp[1,2] = c2;
#     tmp[2,1] = -c1; tmp[2,2] = c1; tmp[2,3] = c2;
#
#     for ix = 3: n-1
#         tmp[ix,ix-1] = -c1; tmp[ix, ix  ] = c1;
#         tmp[ix,ix-2] = -c2; tmp[ix, ix+1] = c2;
#     end
#
#     tmp[n,n-1] = -c1; tmp[n,n] = c1; tmp[n,n-2] = -c2;
#
#     MpxBvx = spdiagm(denum) * kron(sparse(tmp), speye(data_format, m))
#
#     return MpxBvx
# end

# """
#    constructor for finite difference sparse matrix
# """
# function FiniteDiffMatrix(par::PhysicalModel)
#
#     # pml coefficient
#     pml = PMLCoefficient(par)
#
#     # expand physical model parameters
#     lambda = vp2lambda(par.vel, par.rho)
#     lambda = model_expand(lambda , par.npml, par.free_surface)
#     rho    = model_expand(par.rho, par.npml, par.free_surface)
#
#     (MvzBvz, MvzBp ) = Mvz(pml.vz, rho   , par.dz, par.dt)
#     (MvxBvx, MvxBp ) = Mvx(pml.vx, rho   , par.dx, par.dt)
#     (MpzBpz, MpzBvz) = Mpz(pml.pz, lambda, par.dz, par.dt)
#     (MpxBpx, MpxBvx) = Mpx(pml.px, lambda, par.dx, par.dt)
#
#     return FiniteDiffMatrix(MvzBvz, MvzBp , MvxBvx, MvxBp,
#                             MpzBpz, MpzBvz, MpxBpx, MpxBvx)
#
# end
#
# """
#    finite difference sparse matrix with rigit boundary condition
# """
# function RigidFiniteDiffMatrix(par::PhysicalModel)
#
#     # expand physical model parameters
#     lambda = vp2lambda(par.vel, par.rho)
#
#     MvzBp  = Mvz(par.rho, par.dz, par.dt)
#     MvxBp  = Mvx(par.rho, par.dx, par.dt)
#     MpzBvz = Mpz(lambda , par.dz, par.dt)
#     MpxBvx = Mpx(lambda , par.dx, par.dt)
#
#     return RigidFiniteDiffMatrix(MvzBp , MvxBp,
#                                  MpzBvz, MpxBvx)
#
# end


text:
export QLearner,
       GradientBanditLearner,
       MonteCarloLearner, MonteCarloExploringStartLearner, OffPolicyMonteCarloLearner,
       TDLearner, OffPolicyTDLearner, DoubleLearner, QLearner, DifferentialTDLearner, TDλReturnLearner,
       ReinforceLearner,  ReinforceBaselineLearner,
       DQN, DoubleDQN

include("gradient_bandit_learner.jl")
include("monte_carlo_learner.jl")
include("temporal_difference_learner.jl")
include("reinforce_learner.jl")
include("dqn.jl")


text:
module PointGroups

using MLStyle: @match

using ..Crystallography: Triclinic,
    Monoclinic,
    Orthorhombic,
    Tetragonal,
    Cubic,
    Trigonal,
    Hexagonal

export PointGroup
export pointgroups, hermann_mauguin, international, schönflies, schoenflies, orderof, laueclasses

abstract type PointGroup end
abstract type CyclicGroup <: PointGroup end
struct C1 <: CyclicGroup end
struct C2 <: CyclicGroup end
struct C3 <: CyclicGroup end
struct C4 <: CyclicGroup end
struct C6 <: CyclicGroup end
struct C2v <: CyclicGroup end
struct C3v <: CyclicGroup end
struct C4v <: CyclicGroup end
struct C6v <: CyclicGroup end
struct C1h <: CyclicGroup end
struct C2h <: CyclicGroup end
struct C3h <: CyclicGroup end
struct C4h <: CyclicGroup end
struct C6h <: CyclicGroup end
abstract type SpiegelGroup <: PointGroup end
struct S2 <: SpiegelGroup end
struct S4 <: SpiegelGroup end
struct S6 <: SpiegelGroup end
abstract type DihedralGroup <: PointGroup end
struct D2 <: DihedralGroup end
struct D3 <: DihedralGroup end
struct D4 <: DihedralGroup end
struct D6 <: DihedralGroup end
struct D2h <: DihedralGroup end
struct D3h <: DihedralGroup end
struct D4h <: DihedralGroup end
struct D6h <: DihedralGroup end
struct D2d <: DihedralGroup end
struct D3d <: DihedralGroup end
abstract type TetrahedralGroup <: PointGroup end
struct T <: TetrahedralGroup end
struct Th <: TetrahedralGroup end
struct Td <: TetrahedralGroup end
abstract type OctahedralGroup <: PointGroup end
struct O <: OctahedralGroup end
struct Oh <: OctahedralGroup end
const C1v = C1h
const Cs = C1h
const Ci = S2
const C3i = S6
const D1 = C2
const V = D2
const D1h = C2v
const D1d = C2h
const V2 = D2d
const Vh = D2h
function PointGroup(str::AbstractString)
    str = lowercase(filter(!isspace, str))
    return if isnumeric(first(str))  # International notation
        @match str begin
            "1" => C1()
            "1̄" => S2()
            "2" => C2()
            "m" => C1h()
            "3" => C3()
            "4" => C4()
            "4̄" => S4()
            "2/m" => C2h()
            "222" => D2()
            "mm2" => C2v()
            "3̄" => S6()
            "6" => C6()
            "6̄" => C3h()
            "32" => D3()
            "3m" => C3v()
            "mmm" => D2h()
            "4/m" => C4h()
            "422" => D4()
            "4mm" => C4v()
            "4̄2m" => D2d()
            "6/m" => C6h()
            "23" => T()
            "3̄m" => D3d()
            "622" => D6()
            "6mm" => C6v()
            "6̄m2" => D3h()
            "4/mmm" => D4h()
            "6/mmm" => D6h()
            "m3̄" => Th()
            "432" => O()
            "4̄3m" => Td()
            "m3̄m" => Oh()
            _ => throw(ArgumentError("\"$str\" is not a valid point group symbol!"))
        end
    else  # Schoenflies notation
        @match str begin
            "c1" => C1()
            "c2" => C2()
            "c3" => C3()
            "c4" => C4()
            "c6" => C6()
            "c2v" => C2v()
            "c3v" => C3v()
            "c4v" => C4v()
            "c6v" => C6v()
            "c1h" => C1h()
            "c2h" => C2h()
            "c3h" => C3h()
            "c4h" => C4h()
            "c6h" => C6h()
            "s2" => S2()
            "s4" => S4()
            "s6" => S6()
            "d2" => D2()
            "d3" => D3()
            "d4" => D4()
            "d6" => D6()
            "d2h" => D2h()
            "d3h" => D3h()
            "d4h" => D4h()
            "d6h" => D6h()
            "d2d" => D2d()
            "d3d" => D3d()
            "t" => T()
            "th" => Th()
            "td" => Td()
            "o" => O()
            "oh" => Oh()
            "c1v" => C1h()
            "cs" => C1h()
            "ci" => S2()
            "c3i" => S6()
            "d1" => C2()
            "v" => D2()
            "d1h" => C2v()
            "d1d" => C2h()
            "v2" => D2d()
            "vh" => D2h()
            _ => throw(ArgumentError("\"$str\" is not a valid point group symbol!"))
        end
    end
end

struct LaueClass{T<:PointGroup} end

laueclasses(::Triclinic) = (LaueClass{Ci}(),)
laueclasses(::Monoclinic) = (LaueClass{C2h}(),)
laueclasses(::Orthorhombic) = (LaueClass{D2h}(),)
laueclasses(::Tetragonal) = (LaueClass{C4h}(), LaueClass{D4h}())
laueclasses(::Trigonal) = (LaueClass{C3i}(), LaueClass{D3d}())
laueclasses(::Hexagonal) = (LaueClass{C6h}(), LaueClass{D6h}())
laueclasses(::Cubic) = (LaueClass{Th}(), LaueClass{Oh}())

pointgroups(::Triclinic) = (C1(), S2())
pointgroups(::Monoclinic) = (C2(), C1h(), C2h())
pointgroups(::Orthorhombic) = (D2(), C2v(), D2h())
pointgroups(::Tetragonal) = (C4(), S4(), C4h(), D4(), C4v(), D2d(), D4h())
pointgroups(::Trigonal) = (C3(), S6(), D3(), C3v(), D3d())
pointgroups(::Hexagonal) = (C6(), C3h(), C6h(), D6(), C6v(), D3h(), D6h())
pointgroups(::Cubic) = (T(), Th(), O(), Td(), Oh())

hermann_mauguin(::C1) = "1"
hermann_mauguin(::S2) = "1̄"
hermann_mauguin(::C2) = "2"
hermann_mauguin(::C1h) = "m"
hermann_mauguin(::C3) = "3"
hermann_mauguin(::C4) = "4"
hermann_mauguin(::S4) = "4̄"
hermann_mauguin(::C2h) = "2/m"
hermann_mauguin(::D2) = "222"
hermann_mauguin(::C2v) = "mm2"
hermann_mauguin(::S6) = "3̄"
hermann_mauguin(::C6) = "6"
hermann_mauguin(::C3h) = "6̄"
hermann_mauguin(::D3) = "32"
hermann_mauguin(::C3v) = "3m"
hermann_mauguin(::D2h) = "mmm"
hermann_mauguin(::C4h) = "4/m"
hermann_mauguin(::D4) = "422"
hermann_mauguin(::C4v) = "4mm"
hermann_mauguin(::D2d) = "4̄2m"
hermann_mauguin(::C6h) = "6/m"
hermann_mauguin(::T) = "23"
hermann_mauguin(::D3d) = "3̄m"
hermann_mauguin(::D6) = "622"
hermann_mauguin(::C6v) = "6mm"
hermann_mauguin(::D3h) = "6̄m2"
hermann_mauguin(::D4h) = "4/mmm"
hermann_mauguin(::D6h) = "6/mmm"
hermann_mauguin(::Th) = "m3̄"
hermann_mauguin(::O) = "432"
hermann_mauguin(::Td) = "4̄3m"
hermann_mauguin(::Oh) = "m3̄m"
const international = hermann_mauguin

schönflies(g::PointGroup) = isconcretetype(typeof(g)) ? string(nameof(typeof(g))) : throw(ArgumentError("$g is not a specific point group!"))
const schoenflies = schönflies

orderof(::C1) = 1
orderof(::S2) = 2
orderof(::C2) = 2
orderof(::C1h) = 2
orderof(::C3) = 3
orderof(::C4) = 4
orderof(::S4) = 4
orderof(::C2h) = 4
orderof(::D2) = 4
orderof(::C2v) = 4
orderof(::S6) = 6
orderof(::C6) = 6
orderof(::C3h) = 6
orderof(::D3) = 6
orderof(::C3v) = 6
orderof(::D2h) = 8
orderof(::C4h) = 8
orderof(::D4) = 8
orderof(::C4v) = 8
orderof(::D2d) = 8
orderof(::C6h) = 12
orderof(::T) = 12
orderof(::D3d) = 12
orderof(::D6) = 12
orderof(::C6v) = 12
orderof(::D3h) = 12
orderof(::D4h) = 16
orderof(::D6h) = 24
orderof(::Th) = 24
orderof(::O) = 24
orderof(::Td) = 24
orderof(::Oh) = 48

end


text:
using MoreBasisFunctions: vandermonde_matrix, vandermonde_matrix_inverse

n = 11
x = collect(range(0., stop=1., length=n))
y = rand(n)

V = vandermonde_matrix(x)
W = vandermonde_matrix_inverse(x)

x = V\y
a = W*y

@test maximum(abs.((a .- x) ./ a)) < eps(Float32)


text:
tests = ["common", "env", "dbi", "cur"]

for t in tests
    fp = "$t.jl"
    println("* running $fp ...")
    include(fp)
end


text:
using CPIDataBase, Dates
using CPIDataBase.TestHelpers
using DataFrames
using Test

# Creation of types
include("create_types.jl")

# InflationCombination
include("InflationCombination.jl")

# Test some VarCPIBase and CountryStructure operations
include("operations.jl")

# Inflation tests with MixedCountryStructure
include("inflation.jl")

# Trees creation and operations
include("cpitree.jl")


text:
abstract BioUnit                <: Integer

bitstype 8 AminoAcid            <: BioUnit
bitstype 8 Nucleotide           <: BioUnit

## Conversions ##

convert{T<:BioUnit}(::Type{T}, value::UInt8) = box(T,unbox(UInt8,value));
convert{T<:BioUnit}(::Type{UInt8}, value::T) = box(UInt8,unbox(T,value));

const _types = (Int8, Int16, UInt16, Int32, UInt32, Char, Int64, UInt64, Float32, Float64, Bool)

for t in _types
    @eval convert{T<:BioUnit}(::Type{$t}, x::T) = convert($t,UInt8(x))
    @eval convert{T<:BioUnit}(::Type{T},  x::$t)= convert(T,UInt8(x))
end

## Conversions ##

aminoacid{T<:Number}(value::T)   = convert(AminoAcid, value)
nucleotide{T<:Number}(value::T)   = convert(Nucleotide,value)

aminoacid(value::Char)   = convert(AminoAcid, Int(value))
nucleotide(value::Char)   = convert(Nucleotide, Int(value))

aminoacid(vec::Vector)  = convert(Vector{AminoAcid}, vec)
nucleotide(vec::Vector)  = convert(Vector{Nucleotide},vec)

aminoacid(str::ASCIIString)  = convert(Vector{AminoAcid}, str)
nucleotide(str::ASCIIString)  = convert(Vector{Nucleotide},str)

aminoacid(mat::Matrix)  = convert(Matrix{AminoAcid}, mat)
nucleotide(mat::Matrix)  = convert(Matrix{Nucleotide},mat)

macro aa_str(s);  :(reinterpret(AminoAcid,  @b_str($s) )); end

macro nt_str(s);  :(reinterpret(Nucleotide, @b_str($s) )); end

## Promotion rules ##

promote_rule{Tb<:BioUnit,T<:Number}(::Type{Tb}, ::Type{T}) = T
promote_rule{T<:BioUnit}(::Type{T}, ::Type{Union{}}) = T

## Comparisons ##

=={T<:BioUnit}(x::T, y::T) = ==(UInt8(x),UInt8(y))
!={T<:BioUnit}(x::T, y::T) = !=(UInt8(x),UInt8(y))
<{T<:BioUnit}(x::T, y::T)  = <( UInt8(x),UInt8(y))
<={T<:BioUnit}(x::T, y::T) = <=(UInt8(x),UInt8(y))

## Bitwise Operations ##

~{T<:BioUnit}(x::T)            = convert(T,~UInt8(x))
(&){T<:BioUnit}(x::T, y::T)    = convert(T,(&)(UInt8(x),UInt8(y)))
|{T<:BioUnit}(x::T, y::T)      = convert(T,|(UInt8(x),UInt8(y)))
($){T<:BioUnit}(x::T, y::T)    = convert(T,($)(UInt8(x),UInt8(y)))
<<{T<:BioUnit}(x::T, y::Int32) = convert(T,<<(UInt8(x),y))
>>{T<:BioUnit}(x::T, y::Int32) = convert(T,>>(UInt8(x),y))
>>>{T<:BioUnit}(x::T, y::Int32)= convert(T,>>>(UInt8(x),y))

## Characteristics ##

typemin{T<:BioUnit}(::Type{T}) = convert(T, 0x00 )
typemax{T<:BioUnit}(::Type{T}) = convert(T, 0xff )

## Show as Character ##

show{T<:BioUnit}(io::IO,x::T) = (write(io,UInt8(x)); nothing)

## Bits representation ##

bits{T<:BioUnit}(x::T) = bin(reinterpret(UInt8,x),8)

## Array & Strings ##

convert{N}(::Type{Array{AminoAcid,N}}, x::Array{AminoAcid,N}) = x
convert{N}(::Type{Array{Nucleotide,N}}, x::Array{Nucleotide,N}) = x

convert{N}(::Type{Array{AminoAcid,N}}, x::Array{Nucleotide,N}) = reinterpret(AminoAcid,s)
convert{N}(::Type{Array{Nucleotide,N}}, x::Array{AminoAcid,N}) = reinterpret(Nucleotide,s)

convert{N}(::Type{Array{UInt8,N}},s::Array{AminoAcid,N}) = reinterpret(UInt8,s)
convert{N}(::Type{Array{AminoAcid,N}},s::Array{UInt8,N}) = reinterpret(AminoAcid,s)

convert{N}(::Type{Array{UInt8,N}},s::Array{Nucleotide,N}) = reinterpret(UInt8,s)
convert{N}(::Type{Array{Nucleotide,N}},s::Array{UInt8,N}) = reinterpret(Nucleotide,s)

convert{T<:BioUnit}(::Type{ASCIIString}, seq::Vector{T}) = ASCIIString(reinterpret(UInt8,copy(seq)))
convert{T<:BioUnit}(::Type{Vector{T}}, str::ASCIIString) = reinterpret(T,copy(str.data))

bytestring{T<:BioUnit}(seq::Vector{T}) = bytestring(reinterpret(UInt8,seq))

## HASH (for Julia 0.3)
hash{T<:BioUnit}(x::T) = hash(UInt8(x))
hash{T<:BioUnit}(seq::Vector{T}) = hash(bytestring(seq))


text:
# wrapper to convert Conic solver into LPQP solver

# To enable LPQP support from a Conic solver, define, e.g.,
# LinearQuadraticModel(s::ECOSSolver) = ConicToLPQPBridge(ConicModel(s))

mutable struct ConicToLPQPBridge <: AbstractLinearQuadraticModel
    m::AbstractConicModel
    A::SparseMatrixCSC{Float64,Int}
    collb::Vector{Float64}
    colub::Vector{Float64}
    obj::Vector{Float64}
    rowlb::Vector{Float64}
    rowub::Vector{Float64}
    sense::Symbol
    varboundmap::Vector{Int}
    SOCconstrs::Vector{Vector{Int}} # x'x <= y^2, y is first index in vector
    RSOCconstrs::Vector{Vector{Int}} # x'x <= y*z, y is first index in vector, z is second
    vartypes::Vector{Symbol}
end

ConicToLPQPBridge(s::AbstractConicModel) = ConicToLPQPBridge(s, sparse(Int[],Int[],Float64[]), Float64[], Float64[], Float64[], Float64[], Float64[], :Min, Int[], Array{Vector{Int}}(undef, 0),Array{Vector{Int}}(undef, 0), Symbol[])

export ConicToLPQPBridge

# Loads the provided problem data to set up the linear programming problem:
# min c'x
# st  lb <= Ax <= ub
#      l <=  x <= u
# where sense = :Min or :Max
function loadproblem!(wrap::ConicToLPQPBridge, A, collb, colub, obj, rowlb, rowub, sense)
    wrap.A = A
    wrap.collb = collb
    wrap.colub = colub
    wrap.obj = obj
    wrap.rowlb = rowlb
    wrap.rowub = rowub
    wrap.sense = sense
end

const notsoc_error = "For conic solvers, only quadratic constraints in second-order cone format (x'x <= y^2) or rotated second-order cone format (x'x <= yz) are supported"

function addquadconstr!(wrap::ConicToLPQPBridge, linearidx, linearval, quadrowidx, quadcolidx, quadval, sense, rhs)
    if length(linearidx) > 0 || length(linearval) > 0 || sense != '<' || rhs != 0 || count((quadval .== 0.0) .| (quadval .== 1.0)) != length(quadval) - 1 || count(quadval .== -1.0) != 1
        error(notsoc_error)
    end
    length(quadrowidx) == length(quadcolidx) == length(quadval) || error("Inconsistent dimensions")
    # check if SOC or RSOC or neither
    n_pos_on_diag = 0
    off_diag_idx  = 0
    neg_diag_idx  = 0
    n = length(quadval)
    nz = 0
    for i in 1:n
        quadval[i] == 0 && continue
        nz += 1
        if quadrowidx[i] == quadcolidx[i]
            if quadval[i] == 1
                n_pos_on_diag += 1
            elseif quadval[i] == -1
                neg_diag_idx = i
            else
                error(notsoc_error)
            end
        else
            if quadval[i] == -1
                if !(neg_diag_idx == off_diag_idx == 0)
                    error(notsoc_error)
                end
                off_diag_idx = i
            else
                error(notsoc_error)
            end
        end
    end
    if n_pos_on_diag == nz-1 && neg_diag_idx > 0
        # SOC
        SOCconstr = Int[]
        push!(SOCconstr, quadcolidx[neg_diag_idx])
        for i in 1:n
            if quadval[i] == 1
                push!(SOCconstr,quadcolidx[i])
            end
        end
        push!(wrap.SOCconstrs, SOCconstr)
    elseif n_pos_on_diag == nz-1 && off_diag_idx > 0
        # Rotated SOC
        RSOCconstr = Int[]
        push!(RSOCconstr, quadcolidx[off_diag_idx])
        push!(RSOCconstr, quadrowidx[off_diag_idx])
        for i in 1:n
            if quadval[i] == 1
                push!(RSOCconstr,quadcolidx[i])
            end
        end
        push!(wrap.RSOCconstrs, RSOCconstr)

    else
        error(notsoc_error)
    end
end

function optimize!(wrap::ConicToLPQPBridge)
    A = wrap.A
    collb = wrap.collb
    colub = wrap.colub
    obj = wrap.obj
    rowlb = wrap.rowlb
    rowub = wrap.rowub
    empty!(wrap.varboundmap)
    #@show full(A), collb, colub, obj, rowlb, rowub

    (nvar = length(collb)) == length(colub) || error("Unequal lengths for column bounds")
    (nrow = length(rowlb)) == length(rowub) || error("Unequal lengths for row bounds")

    constr_cones = []
    var_cones = [(:Free,1:nvar)]

    # for each variable bound, create a new constraint
    # x <= u  ==> u - x >= 0
    # x >= l  ==> l - x <= 0
    b = Float64[]
    I = Int[] # sparse tuples
    J = Int[]
    extrarows = 0
    for j in 1:nvar
        if collb[j] == colub[j] # Fixed variable
            extrarows += 1
            push!(I,extrarows)
            push!(J,j)
            push!(b,collb[j])
            push!(constr_cones, (:Zero,extrarows:extrarows))
            push!(wrap.varboundmap, j)
            continue
        end
        if collb[j] != -Inf
            # Variable has lower bound
            extrarows += 1
            push!(I,extrarows)
            push!(J,j)
            push!(b, collb[j])
            push!(constr_cones, (:NonPos,extrarows:extrarows))
            push!(wrap.varboundmap, j)
        end
        if colub[j] != Inf
            # Variable has upper bound
            extrarows += 1
            push!(I,extrarows)
            push!(J,j)
            push!(b, colub[j])
            push!(constr_cones, (:NonNeg,extrarows:extrarows))
            push!(wrap.varboundmap, j)
        end
    end
    #@show extrarows

    for it in 1:nrow
        # Equality constraint
        if rowlb[it] == rowub[it]
            # a'x = b ==> b - a'x = 0
            push!(b, rowlb[it])
            push!(constr_cones,(:Zero,it.+(extrarows:extrarows)))
        # Range constraint - not supported
        elseif rowlb[it] != -Inf && rowub[it] != Inf
            error("Ranged constraints unsupported!")
        # Less-than constraint
        elseif rowlb[it] == -Inf
            # a'x <= b ==> b - a'x >= 0
            push!(b, rowub[it])
            push!(constr_cones,(:NonNeg,it.+(extrarows:extrarows)))
        # Greater-than constraint
        else
            # a'x >= b ==> b - a'x <= 0
            push!(b, rowlb[it])
            push!(constr_cones,(:NonPos,it.+(extrarows:extrarows)))
        end
    end

    # Now SOC constraints.
    # We can't trivially apply these to variables because they might have overlapping indices.
    # Instead, append to A.
    ISOC = Int[]
    JSOC = Int[]
    VSOC = Float64[]
    SOCconstrs = wrap.SOCconstrs
    soc_row = 0
    for it in 1:length(SOCconstrs)
        n = length(SOCconstrs[it])
        if collb[SOCconstrs[it][1]] < 0
            error("Invalid second-order conic constraint: x'x <= y^2 requires y >= 0")
        end
        for k in 1:n
            push!(ISOC, soc_row+k)
            push!(JSOC, SOCconstrs[it][k])
            push!(VSOC, -1.0)
            push!(b, 0.0)
        end
        push!(constr_cones, (:SOC, (nrow+extrarows+soc_row+1):(nrow+extrarows+soc_row+n)))
        soc_row += n
    end
    RSOCconstrs = wrap.RSOCconstrs
    for it in 1:length(RSOCconstrs)
        n = length(RSOCconstrs[it])
        if collb[RSOCconstrs[it][1]] < 0 || collb[RSOCconstrs[it][2]] < 0
            error("Invalid rotated second-order conic constraint: x'x <= yz requires y >= 0 and z >= 0")
        end
        for k in 1:n
            push!(ISOC, soc_row+k)
            push!(JSOC, RSOCconstrs[it][k])
            if k <= 2
                push!(VSOC, -1.0/sqrt(2))
            else
                push!(VSOC, -1.0)
            end
            push!(b, 0.0)
        end
        push!(constr_cones, (:SOCRotated, (nrow+extrarows+soc_row+1):(nrow+extrarows+soc_row+n)))
        soc_row += n
    end

    A = vcat(sparse(I,J,ones(extrarows),extrarows, nvar),
             A,
             sparse(ISOC,JSOC,VSOC,soc_row,nvar))

    if wrap.sense == :Max
        obj = -obj
    end

    #@show obj, full(A), b, constr_cones, var_cones
    loadproblem!(wrap.m, obj, A, b, constr_cones, var_cones)
    if !all(t -> t == :Cont, wrap.vartypes)
        setvartype!(wrap.m, wrap.vartypes)
    end
    optimize!(wrap.m)

end

getsolution(wrap::ConicToLPQPBridge) = getsolution(wrap.m)
function getconstrsolution(wrap::ConicToLPQPBridge)
    wrap.A * getsolution(wrap.m)
end
status(wrap::ConicToLPQPBridge) = status(wrap.m)
for f in [:getobjval, :getobjbound]
    @eval function ($f)(wrap::ConicToLPQPBridge)
        if wrap.sense == :Max
            return -$f(wrap.m)
        else
            return $f(wrap.m)
        end
    end
end

function getreducedcosts(wrap::ConicToLPQPBridge)
    redcost = zeros(length(wrap.collb))
    conedual = getdual(wrap.m)
    for i in 1:length(wrap.varboundmap)
        varidx = wrap.varboundmap[i]
        redcost[varidx] -= conedual[i]
    end
    if wrap.sense == :Max
        rmul!(redcost, -1.0)
    end
    return redcost
end

function getconstrduals(wrap::ConicToLPQPBridge)
    constrduals = zeros(length(wrap.rowlb))
    offset = length(wrap.varboundmap)
    conedual = getdual(wrap.m)
    for i in (offset+1):(offset+length(wrap.rowlb))
        constrduals[i-offset] -= conedual[i]
    end
    if wrap.sense == :Max
        rmul!(constrduals, -1.0)
    end
    return constrduals
end

numconstr(wrap::ConicToLPQPBridge) = size(wrap.A, 1)
numvar(wrap::ConicToLPQPBridge) = size(wrap.A, 2)
getvarLB(wrap::ConicToLPQPBridge) = wrap.collb
getvarUB(wrap::ConicToLPQPBridge) = wrap.colub
getconstrLB(wrap::ConicToLPQPBridge) = wrap.rowlb
getconstrUB(wrap::ConicToLPQPBridge) = wrap.rowub
getobj(wrap::ConicToLPQPBridge) = wrap.obj
getsense(wrap::ConicToLPQPBridge) = wrap.sense
function setvarLB!(wrap::ConicToLPQPBridge, l)
    wrap.collb = l
end
function setvarUB!(wrap::ConicToLPQPBridge, u)
    wrap.colub = u
end
function setconstrLB!(wrap::ConicToLPQPBridge, lb)
    wrap.rowlb = lb
end
function setconstrUB!(wrap::ConicToLPQPBridge, ub)
    wrap.rowub = ub
end
function setobj!(wrap::ConicToLPQPBridge, obj)
    wrap.obj = obj
end
function setsense!(wrap::ConicToLPQPBridge, sense)
    wrap.sense = sense
end
function addvar!(wrap::ConicToLPQPBridge, constridx::AbstractArray{T}, constrcoef, l, u, objcoef) where T<:Integer
    wrap.A = [wrap.A sparsevec(constridx, constrcoef, size(wrap.A, 1))]
    push!(wrap.collb, l)
    push!(wrap.colub, u)
    push!(wrap.obj, objcoef)
end
function addconstr!(wrap::ConicToLPQPBridge, varidx::AbstractArray{T}, coef, lb, ub) where T<:Integer
    wrap.A = [wrap.A; sparsevec(varidx, coef, size(wrap.A, 2))']
    push!(wrap.rowlb, lb)
    push!(wrap.rowub, ub)
end
function setvartype!(wrap::ConicToLPQPBridge, v)
    wrap.vartypes = copy(v)
end


text:
#=
acausal_structures:
- Julia version:
- Author: marcin
- Date: 2018-05-10
=#

using BayesNets.CPDs.DiscreteCPD
using BayesNets.DAG
using BayesNets.DiscreteBayesNet

include("discrete_mcpd.jl")
include("hermitian_matrix_distribution.jl")
include("acausal_structure.jl")

export
    AcausalStructure,
    DiscreteQCPD,
    HermitianMatrix,
    event


text:
__precompile__()

module ImageTracking

using Images
using ImageFiltering
using Interpolations
using StaticArrays
using LinearAlgebra
using CoordinateTransformations
using Statistics

abstract type VisualizationMethod end
struct ColorBased <: VisualizationMethod end

abstract type AbstractFlowError end
struct EndpointError <: AbstractFlowError end
struct AngularError <: AbstractFlowError end

abstract type AbstractCoordinateConvention end
struct RasterConvention <: AbstractCoordinateConvention end
struct CartesianConvention <: AbstractCoordinateConvention end

include("core.jl")
include("optical_flow.jl")
include("haar.jl")
include("utility.jl")

export

	# main functions
    optical_flow,
    optical_flow!,

	# other functions
	haar_coordinates,
	haar_features,

	# other functions
	ColorBased,
	EndpointError,
	AngularError,
	polynomial_expansion,
	RasterConvention,
	CartesianConvention,
	visualize_flow,
	read_flow_file,
	evaluate_flow_error,
	calculate_statistics,

	# optical flow algorithms
	LucasKanade,
	Farneback,

    # types that select implementation
    ConvolutionImplementation,
    MatrixImplementation

end