text:
############################ Reads Data from Excel ########################
################################################################################
#**
function get_Cost_Data()
    ks=cost_ks()
    ks.FC_ac=5.8#fixed AC cost
    ks.FC_dc=29#fixed DC cost
    ks.dc=0.25#penalization factor for different than 2 xfrms
    ################################
    #NOTE OSS slope modified here
    #max20%:0.0504
    #max10%:0.0462
    #mean:0.042
    #min10%:0.0378
    #min20%:0.0336
    ks.f_ct=0.042#generating plant variable cost orig:0.0085
    #0.042 is final from data
    ################################
    ks.p_ct=ks.f_ct#substructure variable cost
    ks.c_ct=0.1276#hvdc converter variable cost
    ks.Qc_oss=0.029#M£/MVAr
    ks.Qc_pcc=0.0174#M£/MVAr
    ks.life=25#lifetime of wind farm
    ks.T_op=219150#Operational lifetime in hours
    ks.E_op=0.00005#Energy price E/Wh
    ks.cf=30#Capitalization factor
    ks.FC_bld=10#Build cost
    ks.p2e=1.16#pounds/euro exchange
    #the cost of laying cables varies greatly - the original number is directly from North grid, but entso-e's number is less than half this, DTU tech catalogue is closer to high end entso-e
    #used high end of entso-e
    #HVDC is per cable therefore is *2, while AC is for a 3 phase cable
    ks.ldc=0.32#Cost of laying a single core DC cable
    #####################################################
    #max20%:0.48
    #max10%:0.44
    #mean:0.4
    #min10%:0.36
    #min20%:0.32
    #ks.lac=0.32*2.5#400kV as is 3 conductors #Cost of laying a single 3 phase AC cable#orig ks.lac=0.805
    ks.lac=0.32#220kV HVAC
    #0.032/0.3 is final from data
    ks.lmac=0.32#Cost of laying a single 3 phase MVAC cable
    #####################################################
    #ks.ldc=1.775#Cost of laying a single core DC cable
    #ks.lac=1.905#Cost of laying a single 3 phase AC cable
    #ks.lmac=1.69#Cost of laying a single 3 phase AC cable
    #ks.lmac=1.905#Cost of laying a single 3 phase AC cable
    ks.npv=1.04#discount rate - not used calc is done inline change at npv_hours(), npv_years()
    ks.opx_c=0.025#OPEX percent
    ks.opx_pl=0.02#dc cable
    ks.opx_x=0.0015#Offshore transformer
    ks.opx_co=0.02#offshore converter
    ks.opx_cp=0.007#onshore converter
    ks.conv_c=0.0589#c value of converter cost
    ks.conv_d=54.895#fixed construction converter cost
    ks.pac_e=0.085#e value of ac platform cost
    #ks.pac_e=0.105
    #ks.pac_f=16#f value of ac platform cost
    ################################
    #NOTE OSS offset modified here
    #min:2.5
    #mean:5
    #max:7.5
    ks.pac_f=20#orig=25new value to reflect better the base construction cost - note the variable cost function is also changed for this value see
    #20 is final from data
    #################################
    ks.pdc_g=0.125#g value of dc platform cost
    ks.pdc_h=165#h value of dc platform cost
    return ks
end
###########################################################
######################### cables ##########################
###########################################################
#original data is in excel file but reading file and can access directly or to put in a struct togther use with get_Cable_Data_from_file()
#**
function get_66kV_cables()
    cbls66=[[66.0, 95.0, 170.0, 0.17, 300.0, 187.92, 0.44], [66.0, 120.0, 140.0, 0.18, 340.0, 199.52, 0.43], [66.0, 150.0, 110.0, 0.19, 375.0, 211.12, 0.41], [66.0, 185.0,
90.0, 0.2, 420.0, 227.36, 0.4], [66.0, 240.0, 70.0, 0.22, 480.0, 251.72, 0.38], [66.0, 300.0, 50.0, 0.24, 530.0, 274.92, 0.37], [66.0, 400.0, 40.0, 0.26, 590.0, 306.24, 0.35],
[66.0, 500.0, 30.0, 0.29, 655.0, 345.68, 0.34], [66.0, 630.0, 27.1, 0.32, 715.0, 387.44, 0.33], [66.0, 800.0, 22.5, 0.35, 775.0, 436.16, 0.32], [66.0, 1000.0, 19.8, 0.38, 825.0, 482.56, 0.31]]
    return cbls66
end
#**
function get_220kV_cables()
    cbls220=[[220.0, 400.0, 40.1, 122.0, 665.0, 496.48, 0.457], [220.0, 500.0, 30.0, 136.0, 732.0, 597.4, 0.437], [220.0, 630.0, 27.1, 151.0, 808.0, 638.0, 0.415], [220.0, 800.0, 22.5, 163.0, 879.0, 783.0, 0.4], [220.0, 1000.0, 19.8, 177.0, 942.0, 812.0, 0.386]]
    #cbls220=[[220.0, 400.0, 60.1, 122.0, 10000, 10000000, 0.457]]
    return cbls220
end

#**
function get_400kV_cables()
        cbls400=[[400.0, 500.0, 30.0, 117.0, 776.0, 828.24, 1.42], [400.0, 630.0, 27.1, 125.0, 824.0, 916.4, 1.4], [400.0, 800.0, 22.5, 130.0, 870.0, 997.6, 1.37], [400.0, 1000.0, 19.8, 140.0, 932.0, 1154.2, 1.35], [400.0, 1200.0, 17.3, 170.0, 986.0, 1310.8, 1.33], [400.0, 1400.0, 15.8, 180.0, 1015.0, 1467.4, 1.32], [400.0, 1600.0, 14.7, 190.0, 1036.0, 1624.0, 1.31], [400.0, 2000.0, 10.4, 200.0, 1078.0, 1780.6, 0.44]]
        #cbls400=[[400.0, 2000.0, 13.2, 200.0, Inf, Inf, Inf]]
        return cbls400
end

#**
function get_150kV_cables()
    cbls150=[[150.0, 1000.0, 22.4, 0.7, 1644.0, 179.8, 0.7], [150.0, 1200.0, 19.2, 0.7, 1791.0, 208.8, 0.7], [150.0, 1400.0, 16.5, 0.7, 1962.0, 234.9, 0.7], [150.0, 1600.0, 14.4, 0.7, 2123.0, 261.0, 0.7], [150.0, 2000.0, 11.5, 0.7, 2407.0, 290.0, 0.7]]
    #cbls150=[[150.0, 1000.0, 22.4, 0.7, 1644.0, 10000000, 0.7]]
            return cbls150
end

#**
function get_300kV_cables()
    cbls300=[[300.0, 1000.0, 22.4, 0.7, 1644.0, 263.9, 0.7], [300.0, 1200.0, 19.2, 0.7, 1791.0, 295.8, 0.7], [300.0, 1400.0, 16.5, 0.7, 1962.0, 333.5, 0.7], [300.0, 1600.0, 14.4, 0.7, 2123.0, 371.2, 0.7], [300.0, 2000.0, 11.5, 0.7, 2407.0, 411.8, 0.7]]
    #cbls300=[[300.0, 1000.0, 22.4, 0.7, 1644.0, 10000000, 0.7]]
    return cbls300
end

#Fills in the physical data of a cable into the cable structure
#common data between ac and dc cables - note in dc cables capacitive and inductive reactance is not appliccable
#**
function fillOut_cable_struct(cbl,cbl_data,km,num)
    cbl.elec.volt=cbl_data[1]
    cbl.size=cbl_data[2]
    cbl.elec.ohm=cbl_data[3]*10^-3
    cbl.elec.farrad=cbl_data[4]*10^-9
    cbl.elec.amp=cbl_data[5]
    cbl.costs.perkm_cpx=cbl_data[6]*10^-3
    cbl.length=km
    cbl.elec.henry=cbl_data[7]*10^-3
    cbl.elec.xl=2*pi*cbl.elec.freq*cbl.elec.henry#cable inductive reactance
    cbl.elec.yc=2*pi*cbl.elec.freq*cbl.elec.farrad#cable capacitive reactance
    cbl.num=num
    cbl=get_cbl_failure_data(cbl)#Set failure data
    return cbl
end

#Fills in the physical data of an ac cable into the cable structure **
#**
function fillOut_cable_struct_ac(cbl,cbl_data,km,num)
    cbl=fillOut_cable_struct(cbl,cbl_data,km,num)#common data between ac and dc cables
    cbl.elec.mva=get_newQ_Capacity(cbl.elec.freq,km,cbl.elec.volt,cbl.elec.farrad,cbl.elec.amp)
    return cbl
end

#Fills in the physical data of an hvdc cable into the cable structure **
#**
function fillOut_cable_struct_dc(cbl,cbl_data,km,num)
    cbl=fillOut_cable_struct(cbl,cbl_data,km,num)#common data between ac and dc cables
    cbl.elec.mva=cbl.elec.volt*cbl.elec.amp*10^-3#dc does not change with distance
    return cbl
end

#**
function get_cbl_failure_data(cbl)
    #cbl.relia.fr=(0.1114/100)*cbl.length#/yr/100km
    #cbl.relia.mttr=1.971#months
    cbl.relia.fr=(0.1/100)*cbl.length#/yr/100km
    cbl.relia.mttr=2#months
    cbl.relia.mc=0.56#No longer used
    return cbl
end

#Calculates the new hvac cable capacity after 50-50 compensation at distance km. **
#**
function get_newQ_Capacity(f,l,v,q,a)
#Calculates the square of new hvac cable capacity after 50-50 compensation at distance km.
    mva=(sqrt(3)*v*10^3*a/10^6)^2-((0.5*((v*10^3)^2*2*pi*f*l*q))/10^6)^2
#takes square root if negative returns zero if negative
    if mva>=0
        mva=sqrt(mva)
    else
        mva=0.0
    end
 return mva
end

###########################################################
################### transformers ##########################
###########################################################
#Sets all options for transformer sizes in 10MVA steps **
#Max size considered 5000MVA => although not realistic it's to test 10GW HVDC
#**
function  get_Xfo_Data()
    xfos=Array{Float32,1}()
    push!(xfos,10)
    for i=50:10:1000
        push!(xfos,i)
    end
    #push!(xfos,2500)
    #push!(xfos,5000)
    return xfos
end

#failure data for transformers **
#**
function get_offshore_xfo_failure_data(xfo)
    #failure rate is faily consistent among sources but MTTR varies alot
    #Reference from North Sea Grid Annex - 120hrs
    #6 months is most commnly assumed value but seems old, more recently it's dropped to 2 or 3
    #Energy Transmission and Grid Integration of AC Offshore Wind Farms
    #Integration of Large Scale Wind Energy with Electrical Power Systems in China
    #Reliability Study Analysis of Electrical Systems within Offshore Wind Parks Elforsk report 07:65
    xfo.relia.fr=0.02#/yr# from CIGRE
    xfo.relia.mttr=3#month - assumed as it's middle of the values
    #xfo.relia.mc=2.8#
    #xfo.relia.mc=0.05#Not directly sourced but derived from Operation and maintenance requirements (No longer used)
    return xfo.relia
end

#failure data for transformers **
#**
function get_onshore_xfo_failure_data(xfo)
    #failure rate is faily consistent among sources but MTTR varies alot
    #Reference from North Sea Grid Annex - 120hrs
    #6 months is most commnly assumed value but seems old, more recently it's dropped to 2 or 3
    #Energy Transmission and Grid Integration of AC Offshore Wind Farms
    #Integration of Large Scale Wind Energy with Electrical Power Systems in China
    #Reliability Study Analysis of Electrical Systems within Offshore Wind Parks Elforsk report 07:65
    xfo.relia.fr=0.02#/yr# from CIGRE
    xfo.relia.mttr=0.6#month - assumed as it's middle of the values
    #xfo.relia.mc=2.8#
    #xfo.relia.mc=0.05#Not directly sourced but derived from Operation and maintenance requirements (No longer used)
    return xfo.relia
end

#for now both the same for onshore as offshore until etter data is obtained

#**
function get_onshore_conv_failure_data(conv)
    conv.relia.fr=0.12#/yr
    conv.relia.mttr=0.6#month
    conv.relia.mc=0.56#
    return conv.relia
end

#**
function get_offshore_conv_failure_data(conv)
    conv.relia.fr=0.12#/yr
    conv.relia.mttr=2#month
    conv.relia.mc=0.56#
    return conv.relia
end

#sizes and costs reactors - same reactors each side
#set of reactors
#All reactor sizes and cost from North Sea Grid Final report AnnexAnnex
#**
function get_reactors(Q)
    reactors=Array{Int32,1}()
    qc=0
    while (Q>275)
        push!(reactors,350)
        qc=qc+10.065
        Q=Q-350
    end
    while (Q>150)
        push!(reactors,200)
        qc=qc+5.75
        Q=Q-200
    end
    while (Q>125)
        push!(reactors,150)
        qc=qc+4.315
        Q=Q-150
    end
    while (Q>75)
        push!(reactors,125)
        qc=qc+3.595
        Q=Q-125
    end
    while (Q>50)
        push!(reactors,75)
        qc=qc+2.155
        Q=Q-75
    end
    while (Q>0)
        push!(reactors,50)
        qc=qc+1.437
        Q=Q-50
    end
    if ((length(reactors)==0) && (Q>0 && Q<50))
        push!(reactors,50)
        qc=qc+1.437
    elseif (Q>0)
        println("Get Reactor Function, Problem assigning reactors.")
    end
    return qc,reactors
end

#Gets the data for equipment to be used, not used unless explicitly called, only for updating to new data if desired
#**
function get_Cable_Data_from_file()
    df = DataFrame(XLSX.readtable("/Users/shardy/Documents/GitHub/hv_offshore_topology_optimization/v1.0/greedy_search/economics/input_data/cost_data.xlsx", "eqp_data")...)
    eqps=cbls_data()#create an instance of object ks
    for index=1:length(df[:, 1])
        #get 33kv cable data
        if (typeof(df.kv_33[index]) != Missing)
            cbd=Float32[]
            push!(cbd,df.kv_33[index])
            push!(cbd,df.mm2_33[index])
            push!(cbd,df.ohms_per_km_33[index])
            push!(cbd,df.nf_per_km_33[index])
            push!(cbd,df.Amps_33[index])
            push!(cbd,df.euro_per_m_33[index])
            push!(cbd,df.mh_per_km_33[index])
            push!(eqps.cbls33kV,deepcopy(cbd))
        end

        #get 66kv cable data
        if (typeof(df.kv_66[index]) != Missing)
            cbd=Float32[]
            push!(cbd,df.kv_66[index])
            push!(cbd,df.mm2_66[index])
            push!(cbd,df.ohms_per_km_66[index])
            push!(cbd,df.nf_per_km_66[index])
            push!(cbd,df.Amps_66[index])
            push!(cbd,df.euro_per_m_66[index])
            push!(cbd,df.mh_per_km_66[index])
            push!(eqps.cbls66kV,deepcopy(cbd))
        end

        #get 132kv cable data
        if (typeof(df.kv_132[index]) != Missing)
            cbd=Float32[]
            push!(cbd,df.kv_132[index])
            push!(cbd,df.mm2_132[index])
            push!(cbd,df.ohms_per_km_132[index])
            push!(cbd,df.nf_per_km_132[index])
            push!(cbd,df.Amps_132[index])
            push!(cbd,df.euro_per_m_132[index])
            push!(cbd,df.mh_per_km_132[index])
            push!(eqps.cbls132kV,deepcopy(cbd))
        end

        #get 220kv cable data
        if (typeof(df.kv_220[index]) != Missing)
            cbd=Float32[]
            push!(cbd,df.kv_220[index])
            push!(cbd,df.mm2_220[index])
            push!(cbd,df.ohms_per_km_220[index])
            push!(cbd,df.nf_per_km_220[index])
            push!(cbd,df.Amps_220[index])
            push!(cbd,df.euro_per_m_220[index])
            push!(cbd,df.mh_per_km_220[index])
            push!(eqps.cbls220kV,deepcopy(cbd))
        end

        #get 400kv cable data
        if (typeof(df.kv_400[index]) != Missing)
            cbd=Float32[]
            push!(cbd,df.kv_400[index])
            push!(cbd,df.mm2_400[index])
            push!(cbd,df.ohms_per_km_400[index])
            push!(cbd,df.nf_per_km_400[index])
            push!(cbd,df.Amps_400[index])
            push!(cbd,df.euro_per_m_400[index])
            push!(cbd,df.mh_per_km_400[index])
            push!(eqps.cbls400kV,deepcopy(cbd))
        end

        #get 150kv cable data
        if (typeof(df.kv_150[index]) != Missing)
            cbd=Float32[]
            push!(cbd,df.kv_150[index])
            push!(cbd,df.mm2_150[index])
            push!(cbd,df.ohms_per_km_150[index])
            push!(cbd,df.nf_per_km_150[index])
            push!(cbd,df.Amps_150[index])
            push!(cbd,df.euro_per_m_150[index])
            push!(cbd,df.mh_per_km_150[index])
            push!(eqps.cbls150kV,deepcopy(cbd))
        end

        #get 300kv cable data
        if (typeof(df.kv_300[index]) != Missing)
            cbd=Float32[]
            push!(cbd,df.kv_300[index])
            push!(cbd,df.mm2_300[index])
            push!(cbd,df.ohms_per_km_300[index])
            push!(cbd,df.nf_per_km_300[index])
            push!(cbd,df.Amps_300[index])
            push!(cbd,df.euro_per_m_300[index])
            push!(cbd,df.mh_per_km_300[index])
            push!(eqps.cbls300kV,deepcopy(cbd))
        end
    end
    return eqps
end


text:
function get_orbital(;
    timestep=0.01,
    gravity=-9.81,
    spring=0.0,
    damper=0.0,
    num_bodies=5,
    T=Float64)

    # Parameters
    ex = [0.0; 0.0; 1.0]
    h = 1.0
    r = 0.05
    vert11 = [0.0; 0.0; h / 2.0]
    vert12 = -vert11

    # Links
    origin = Origin{T}()

    bodies = [Box(r, r, h, h, color=RGBA(1.0, 0.0, 0.0)) for i = 1:num_bodies]

    # Constraints
    jointb1 = JointConstraint(Fixed(origin, bodies[1];
        child_vertex=vert11))
    if num_bodies > 1
        joints = [
            jointb1;
            [
                JointConstraint(Orbital(bodies[i - 1], bodies[i], ex;
                parent_vertex=vert12,
                child_vertex=vert11,
                spring=spring,
                damper=damper)) for i = 2:num_bodies]
            ]
    else
        joints = [jointb1]
    end

    mech = Mechanism(origin, bodies, joints,
        gravity=gravity,
        timestep=timestep,
        spring=spring,
        damper=damper)

    return mech
end

function initialize_orbital!(mechanism::Mechanism{T};
    orientation=[pi / 4.0, pi / 8.0]) where T

    pbody = mechanism.bodies[1]
    joint = mechanism.joints[1]
    vert11 = joint.translational.vertices[2]
    vert12 = -vert11

    # set position and velocities
    set_maximal_configurations!(mechanism.origin, pbody,
        child_vertex=vert11,
        Δq=RotX(0.0))

    set_minimal_coordinates!(mechanism, mechanism.joints[2], orientation)

    zero_velocity!(mechanism)
end

ext: jl


text:
using JSON
@pyimport parameters_julia
out_dir = parameters_julia.out_dir
files_ID = parameters_julia.files_ID
month_w = parameters_julia.month_w
#instance1 = readstring(out_dir * "instance_comm.txt")

link_label_dict, link_label_dict_, link_length_dict = furInfo()

#load OD pair labels
odPairLabel = readstring(out_dir * "od_pair_label_dict_refined.json")
odPairLabel = JSON.parse(odPairLabel)

odPairLabel_ = readstring(out_dir * "od_pair_label_dict__refined.json")
odPairLabel_ = JSON.parse(odPairLabel_)

#load node-link incidence
nodeLink = readstring(out_dir * "node_link_incidence.json");
nodeLink = JSON.parse(nodeLink);

function extract_demandDict(out_dir, files_ID, month_w ,day, instance1, numZones)
	demandsDict = Dict()
	# get ground trueth demands, indexed by 0
	demandsDict[0] = iniDemand(out_dir * "data_traffic_assignment_uni-class/" * files_ID * "_trips_" * month_w * "_" * string(day) * "_" * instance1 * ".txt", numZones)
	# get initial demands, indexed by 1
	demandsDict[1] = iniDemand(out_dir * "data_traffic_assignment_uni-class/" * files_ID * "_trips_" * month_w * "_" * string(day) * "_" * instance1 * ".txt", numZones, 1)

	return demandsDict
end

function extract_dataf(out_dir, files_ID, month_w ,day, instance1, demandsDict)
	numNodes, numLinks, numODpairs, capacity, free_flow_time, ta_data_Apr_PM = paraNetwork(out_dir, files_ID, month_w, day, instance1, demandsDict)
	start_node = ta_data_Apr_PM.start_node
	end_node = ta_data_Apr_PM.end_node

	return numNodes, numLinks, numODpairs, capacity, free_flow_time, ta_data_Apr_PM, start_node, end_node
end


text:
# Default implementations of `sample`.
const PROGRESS = Ref(true)

"""
    setprogress!(progress::Bool)

Enable progress logging globally if `progress` is `true`, and disable it otherwise.
"""
function setprogress!(progress::Bool)
    @info "progress logging is $(progress ? "enabled" : "disabled") globally"
    PROGRESS[] = progress
    return progress
end

function StatsBase.sample(
    model::AbstractModel,
    sampler::AbstractSampler,
    arg;
    kwargs...
)
    return StatsBase.sample(Random.GLOBAL_RNG, model, sampler, arg; kwargs...)
end

"""
    sample([rng, ]model, sampler, N; kwargs...)

Return `N` samples from the `model` with the Markov chain Monte Carlo `sampler`.
"""
function StatsBase.sample(
    rng::Random.AbstractRNG,
    model::AbstractModel,
    sampler::AbstractSampler,
    N::Integer;
    kwargs...
)
    return mcmcsample(rng, model, sampler, N; kwargs...)
end

"""
    sample([rng, ]model, sampler, isdone; kwargs...)

Sample from the `model` with the Markov chain Monte Carlo `sampler` until a
convergence criterion `isdone` returns `true`, and return the samples.

The function `isdone` has the signature
```julia
isdone(rng, model, sampler, samples, iteration; kwargs...)
```
and should return `true` when sampling should end, and `false` otherwise.
"""
function StatsBase.sample(
    rng::Random.AbstractRNG,
    model::AbstractModel,
    sampler::AbstractSampler,
    isdone;
    kwargs...
)
    return mcmcsample(rng, model, sampler, isdone; kwargs...)
end

function StatsBase.sample(
    model::AbstractModel,
    sampler::AbstractSampler,
    parallel::AbstractMCMCEnsemble,
    N::Integer,
    nchains::Integer;
    kwargs...
)
    return StatsBase.sample(Random.GLOBAL_RNG, model, sampler, parallel, N, nchains;
                            kwargs...)
end

"""
    sample([rng, ]model, sampler, parallel, N, nchains; kwargs...)

Sample `nchains` Monte Carlo Markov chains from the `model` with the `sampler` in parallel
using the `parallel` algorithm, and combine them into a single chain.
"""
function StatsBase.sample(
    rng::Random.AbstractRNG,
    model::AbstractModel,
    sampler::AbstractSampler,
    parallel::AbstractMCMCEnsemble,
    N::Integer,
    nchains::Integer;
    kwargs...
)
    return mcmcsample(rng, model, sampler, parallel, N, nchains; kwargs...)
end

# Default implementations of regular and parallel sampling.

function mcmcsample(
    rng::Random.AbstractRNG,
    model::AbstractModel,
    sampler::AbstractSampler,
    N::Integer;
    progress = PROGRESS[],
    progressname = "Sampling",
    callback = nothing,
    discard_initial = 0,
    thinning = 1,
    chain_type::Type=Any,
    kwargs...
)
    # Check the number of requested samples.
    N > 0 || error("the number of samples must be ≥ 1")
    Ntotal = thinning * (N - 1) + discard_initial + 1

    # Start the timer
    start = time()
    local state

    @ifwithprogresslogger progress name=progressname begin
        # Determine threshold values for progress logging
        # (one update per 0.5% of progress)
        if progress
            threshold = Ntotal ÷ 200
            next_update = threshold
        end

        # Obtain the initial sample and state.
        sample, state = step(rng, model, sampler; kwargs...)

        # Discard initial samples.
        for i in 1:(discard_initial - 1)
            # Update the progress bar.
            if progress && i >= next_update
                ProgressLogging.@logprogress i/Ntotal
                next_update = i + threshold
            end

            # Obtain the next sample and state.
            sample, state = step(rng, model, sampler, state; kwargs...)
        end

        # Run callback.
        callback === nothing || callback(rng, model, sampler, sample, state, 1; kwargs...)

        # Save the sample.
        samples = AbstractMCMC.samples(sample, model, sampler, N; kwargs...)
        samples = save!!(samples, sample, 1, model, sampler, N; kwargs...)

        # Update the progress bar.
        itotal = 1 + discard_initial
        if progress && itotal >= next_update
            ProgressLogging.@logprogress itotal / Ntotal
            next_update = itotal + threshold
        end

        # Step through the sampler.
        for i in 2:N
            # Discard thinned samples.
            for _ in 1:(thinning - 1)
                # Obtain the next sample and state.
                sample, state = step(rng, model, sampler, state; kwargs...)

                # Update progress bar.
                if progress && (itotal += 1) >= next_update
                    ProgressLogging.@logprogress itotal / Ntotal
                    next_update = itotal + threshold
                end
            end

            # Obtain the next sample and state.
            sample, state = step(rng, model, sampler, state; kwargs...)

            # Run callback.
            callback === nothing || callback(rng, model, sampler, sample, state, i; kwargs...)

            # Save the sample.
            samples = save!!(samples, sample, i, model, sampler, N; kwargs...)

            # Update the progress bar.
            if progress && (itotal += 1) >= next_update
                ProgressLogging.@logprogress itotal / Ntotal
                next_update = itotal + threshold
            end
        end
    end

    # Get the sample stop time.
    stop = time()
    duration = stop - start
    stats = SamplingStats(start, stop, duration)

    return bundle_samples(
        samples,
        model,
        sampler,
        state,
        chain_type;
        stats=stats,
        discard_initial=discard_initial,
        thinning=thinning,
        kwargs...
    )
end

function mcmcsample(
    rng::Random.AbstractRNG,
    model::AbstractModel,
    samplers::Vector{<:AbstractSampler},
    swap_every::Integer,
    N::Integer;
    progress = PROGRESS[],
    progressname = "Replica Exchange: Sampling ($(length(samplers)) replicas)",
    callback = nothing,
    discard_initial = 0,
    thinning = 1,
    chain_type::Type=Any,
    kwargs...
)
    """
    Replica exchange
    """
    println("Replica Exchange AbstractMCMC.mcmcsample")
    function swap_β(samplers::Vector{<:AbstractSampler}, states, start::Integer)
        L = length(samplers) - 1
        for sampler_id in start:2:L
            println("Attempting swap b/w $sampler_id, $(sampler_id+1)")
            T = typeof(states[sampler_id].z.ℓπ)
            for (name, typ) in zip(fieldnames(T), T.types)
                println("type of the fieldname $name is $typ")
            end
            println(states[sampler_id].z.ℓπ.value)
            logα = - (samplers[sampler_id].alg.β - samplers[sampler_id + 1].alg.β) * (states[sampler_id].z.ℓπ.value - states[sampler_id+1].z.ℓπ.value)
            if log(1-Random.rand(rng)) ≤ logα
                println("Successful swap b/w $sampler_id, $(sampler_id+1)")
                @set samplers[sampler_id].alg.β, samplers[sampler_id + 1].alg.β = samplers[sampler_id + 1].alg.β, samplers[sampler_id].alg.β
            else
                println("Failed swap b/w $sampler_id, $(sampler_id+1)")
            end
        end
    end
    # Check the number of requested samples.
    N > 0 || error("the number of samples must be ≥ 1")
    Ntotal = thinning * (N - 1) + discard_initial + 1

    # Start the timer
    start = time()
    local states, state, itotal
    states, samples_per_beta = [], Dict()
    βs = [sampler.alg.β for sampler in samplers]
    @ifwithprogresslogger progress name=progressname begin
        # Determine threshold values for progress logging
        # (one update per 0.5% of progress)
        if progress
            threshold = Ntotal ÷ 200
            next_update = threshold
        end

        for (sampler_id, sampler) in enumerate(samplers)
            # Obtain the initial sample and state.
            sample, state = step(rng, model, sampler; kwargs...)

            push!(states, state)

            # Discard initial samples.
            for i in 1:(discard_initial - 1)
                # Update the progress bar.
                if progress && i >= next_update && sampler_id == length(samplers)
                    ProgressLogging.@logprogress i/Ntotal
                    next_update = i + threshold
                end

                # Obtain the next sample and state.
                sample, states[sampler_id] = step(rng, model, sampler, states[sampler_id]; kwargs...)
            end

            # Run callback.
            callback === nothing || callback(rng, model, sampler, sample, states[sampler_id], 1; kwargs...)

            # Save the sample.
            samples = AbstractMCMC.samples(sample, model, sampler, N; kwargs...)
            samples = save!!(samples, sample, 1, model, sampler, N; kwargs...)
            samples_per_beta[sampler.alg.β] = samples
            # Update the progress bar if it is the last replica
            itotal = 1 + discard_initial
            println(itotal)
            if progress && itotal >= next_update && sampler_id == length(samplers)
                ProgressLogging.@logprogress itotal / Ntotal
                next_update = itotal + threshold
            end
        end
        println(itotal)
        # Step through the sampler.
        for i in 2:N
            # # Discard thinned samples.
            # for _ in 1:(thinning - 1)
            #     # Obtain the next sample and state.
            #     sample, state = step(rng, model, sampler, state; kwargs...)

            #     # Update progress bar.
            #     if progress && (itotal += 1) >= next_update
            #         ProgressLogging.@logprogress itotal / Ntotal
            #         next_update = itotal + threshold
            #     end
            # end

            if i % swap_every == 0
                if i % (2*swap_every) == 0
                    swap_β(samplers, states, 2) # swap even indices
                else
                    swap_β(samplers, states, 1) # swap odd indices
                end
            end
            for (sampler_id, sampler) in enumerate(samplers)
                # Load state of replica 'sampler_id'
                state = states[sampler_id]

                # Obtain the next sample and state for the replica
                sample, state = step(rng, model, sampler, state; kwargs...)

                # Run callback for the replica
                callback === nothing || callback(rng, model, sampler, sample, state, i; kwargs...)

                # Save state of replica 'sampler_id'
                states[sampler_id] = state

                # Save the sample for the replica
                samples_per_beta[sampler.alg.β] = save!!(samples_per_beta[sampler.alg.β], sample, i, model, sampler, N; kwargs...)

                # Update the progress bar if it is the last replica
                if progress && (itotal += 1) >= next_update && sampler_id == length(samplers)
                    ProgressLogging.@logprogress itotal / Ntotal
                    next_update = itotal + threshold
                end
            end
        end
    end

    # Get the sample stop time.
    stop = time()
    duration = stop - start
    stats = SamplingStats(start, stop, duration)

    # println(samples_per_beta)

    bundled_samples = Dict()
    for β in βs
        for (sampler_id, sampler) in enumerate(samplers)
            if sampler.alg.β == β
                state = states[sampler_id]
                samples = samples_per_beta[β]
                bundled_samples[β] =  bundle_samples(
                    samples,
                    model,
                    sampler,
                    state,
                    chain_type;
                    stats=stats,
                    discard_initial=discard_initial,
                    thinning=thinning,
                    kwargs...
                )
            end
        end
    end
    return bundled_samples
end

function mcmcsample(
    rng::Random.AbstractRNG,
    model::AbstractModel,
    sampler::AbstractSampler,
    isdone;
    chain_type::Type=Any,
    progress = PROGRESS[],
    progressname = "Convergence sampling",
    callback = nothing,
    discard_initial = 0,
    thinning = 1,
    kwargs...
)

    # Start the timer
    start = time()
    local state

    @ifwithprogresslogger progress name=progressname begin
        # Obtain the initial sample and state.
        sample, state = step(rng, model, sampler; kwargs...)

        # Discard initial samples.
        for _ in 2:discard_initial
            # Obtain the next sample and state.
            sample, state = step(rng, model, sampler, state; kwargs...)
        end

        # Run callback.
        callback === nothing || callback(rng, model, sampler, sample, state, 1; kwargs...)

        # Save the sample.
        samples = AbstractMCMC.samples(sample, model, sampler; kwargs...)
        samples = save!!(samples, sample, 1, model, sampler; kwargs...)

        # Step through the sampler until stopping.
        i = 2

        while !isdone(rng, model, sampler, samples, state, i; progress=progress, kwargs...)
            # Discard thinned samples.
            for _ in 1:(thinning - 1)
                # Obtain the next sample and state.
                sample, state = step(rng, model, sampler, state; kwargs...)
            end

            # Obtain the next sample and state.
            sample, state = step(rng, model, sampler, state; kwargs...)

            # Run callback.
            callback === nothing || callback(rng, model, sampler, sample, state, i; kwargs...)

            # Save the sample.
            samples = save!!(samples, sample, i, model, sampler; kwargs...)

            # Increment iteration counter.
            i += 1
        end
    end

    # Get the sample stop time.
    stop = time()
    duration = stop - start
    stats = SamplingStats(start, stop, duration)

    # Wrap the samples up.
    return bundle_samples(
        samples,
        model,
        sampler,
        state,
        chain_type;
        stats=stats,
        discard_initial=discard_initial,
        thinning=thinning,
        kwargs...
    )
end

function mcmcsample(
    rng::Random.AbstractRNG,
    model::AbstractModel,
    sampler::AbstractSampler,
    ::MCMCThreads,
    N::Integer,
    nchains::Integer;
    progress = PROGRESS[],
    progressname = "Sampling ($(min(nchains, Threads.nthreads())) threads)",
    kwargs...
)
    # Check if actually multiple threads are used.
    if Threads.nthreads() == 1
        @warn "Only a single thread available: MCMC chains are not sampled in parallel"
    end

    # Check if the number of chains is larger than the number of samples
    if nchains > N
        @warn "Number of chains ($nchains) is greater than number of samples per chain ($N)"
    end

    # Copy the random number generator, model, and sample for each thread
    # NOTE: As of May 17, 2020, this relies on Julia's thread scheduling functionality
    #       that distributes a for loop into equal-sized blocks and allocates them
    #       to each thread. If this changes, we may need to rethink things here.
    interval = 1:min(nchains, Threads.nthreads())
    rngs = [deepcopy(rng) for _ in interval]
    models = [deepcopy(model) for _ in interval]
    samplers = [deepcopy(sampler) for _ in interval]

    # Create a seed for each chain using the provided random number generator.
    seeds = rand(rng, UInt, nchains)

    # Set up a chains vector.
    chains = Vector{Any}(undef, nchains)

    @ifwithprogresslogger progress name=progressname begin
        # Create a channel for progress logging.
        if progress
            channel = Channel{Bool}(length(interval))
        end

        Distributed.@sync begin
            if progress
                # Update the progress bar.
                Distributed.@async begin
                    # Determine threshold values for progress logging
                    # (one update per 0.5% of progress)
                    threshold = nchains ÷ 200
                    nextprogresschains = threshold

                    progresschains = 0
                    while take!(channel)
                        progresschains += 1
                        if progresschains >= nextprogresschains
                            ProgressLogging.@logprogress progresschains/nchains
                            nextprogresschains = progresschains + threshold
                        end
                    end
                end
            end

            Distributed.@async begin
                try
                    Threads.@threads for i in 1:nchains
                        # Obtain the ID of the current thread.
                        id = Threads.threadid()

                        # Seed the thread-specific random number generator with the pre-made seed.
                        subrng = rngs[id]
                        Random.seed!(subrng, seeds[i])

                        # Sample a chain and save it to the vector.
                        chains[i] = StatsBase.sample(subrng, models[id], samplers[id], N;
                                                     progress = false, kwargs...)

                        # Update the progress bar.
                        progress && put!(channel, true)
                    end
                finally
                    # Stop updating the progress bar.
                    progress && put!(channel, false)
                end
            end
        end
    end

    # Concatenate the chains together.
    return chainsstack(tighten_eltype(chains))
end

function mcmcsample(
    rng::Random.AbstractRNG,
    model::AbstractModel,
    sampler::AbstractSampler,
    ::MCMCDistributed,
    N::Integer,
    nchains::Integer;
    progress = PROGRESS[],
    progressname = "Sampling ($(Distributed.nworkers()) processes)",
    kwargs...
)
    # Check if actually multiple processes are used.
    if Distributed.nworkers() == 1
        @warn "Only a single process available: MCMC chains are not sampled in parallel"
    end

    # Check if the number of chains is larger than the number of samples
    if nchains > N
        @warn "Number of chains ($nchains) is greater than number of samples per chain ($N)"
    end

    # Create a seed for each chain using the provided random number generator.
    seeds = rand(rng, UInt, nchains)

    # Set up worker pool.
    pool = Distributed.CachingPool(Distributed.workers())

    local chains
    @ifwithprogresslogger progress name=progressname begin
        # Create a channel for progress logging.
        if progress
            channel = Distributed.RemoteChannel(() -> Channel{Bool}(Distributed.nworkers()))
        end

        Distributed.@sync begin
            if progress
                # Update the progress bar.
                Distributed.@async begin
                    # Determine threshold values for progress logging
                    # (one update per 0.5% of progress)
                    threshold = nchains ÷ 200
                    nextprogresschains = threshold

                    progresschains = 0
                    while take!(channel)
                        progresschains += 1
                        if progresschains >= nextprogresschains
                            ProgressLogging.@logprogress progresschains/nchains
                            nextprogresschains = progresschains + threshold
                        end
                    end
                end
            end

            Distributed.@async begin
                try
                    chains = Distributed.pmap(pool, seeds) do seed
                        # Seed a new random number generator with the pre-made seed.
                        Random.seed!(rng, seed)

                        # Sample a chain.
                        chain = StatsBase.sample(rng, model, sampler, N;
                                                 progress = false, kwargs...)

                        # Update the progress bar.
                        progress && put!(channel, true)

                        # Return the new chain.
                        return chain
                    end
                finally
                    # Stop updating the progress bar.
                    progress && put!(channel, false)
                end
            end
        end
    end

    # Concatenate the chains together.
    return chainsstack(tighten_eltype(chains))
end

function mcmcsample(
    rng::Random.AbstractRNG,
    model::AbstractModel,
    sampler::AbstractSampler,
    ::MCMCSerial,
    N::Integer,
    nchains::Integer;
    progressname = "Sampling",
    kwargs...
)
    # Check if the number of chains is larger than the number of samples
    if nchains > N
        @warn "Number of chains ($nchains) is greater than number of samples per chain ($N)"
    end

    # Sample the chains.
    chains = map(
        i -> StatsBase.sample(rng, model, sampler, N; progressname = string(progressname, " (Chain ", i, " of ", nchains, ")"),
        kwargs...),
        1:nchains
    )

    # Concatenate the chains together.
    return chainsstack(tighten_eltype(chains))
end

tighten_eltype(x) = x
tighten_eltype(x::Vector{Any}) = map(identity, x)


text:
"""
    increase_time!(status::KalmanStatus)
    increase_time!(status::SizedKalmanStatus)

Increase time by one unit.
"""
increase_time!(status::KalmanStatus) = status.t += 1;
increase_time!(status::SizedKalmanStatus) = increase_time!(status.online_status);

"""
    apriori(X::FloatVector, settings::KalmanSettings)

Kalman filter a-priori prediction for X.

# Arguments
- `X`: Last expected value of the states
- `settings`: KalmanSettings struct

    apriori(P::SymMatrix, settings::KalmanSettings)

Kalman filter a-priori prediction for P.

# Arguments
- `P`: Last conditional covariance the states
- `settings`: KalmanSettings struct
"""
apriori(X::FloatVector, settings::KalmanSettings) = settings.C * X;
apriori(P::SymMatrix, settings::KalmanSettings) = Symmetric(settings.C * P * settings.C' + settings.DQD)::SymMatrix;

"""
    initialise_apriori!(settings::KalmanSettings, status::KalmanStatus)

Initialise X_prior, P_prior and loglik.
"""
function initialise_apriori!(settings::KalmanSettings, status::KalmanStatus)

    # First a-priori prediction for X and P
    status.X_prior = apriori(settings.X0, settings);
    status.P_prior = apriori(settings.P0, settings);

    # Initialise buffers
    status.buffer_J1 = similar(status.X_prior);
    status.buffer_J2 = similar(status.P_prior);
    status.buffer_m_m = similar(status.P_prior);

    # Initialise loglikelihood
    if settings.compute_loglik == true
        status.loglik = 0.0;
    end
end

"""
    initialise_status_history!(settings::KalmanSettings, status::OnlineKalmanStatus)
    initialise_status_history!(settings::KalmanSettings, status::DynamicKalmanStatus)
    initialise_status_history!(settings::KalmanSettings, status::SizedKalmanStatus)

Initialise the `status.history_*` arrays when needed.
"""
function initialise_status_history!(settings::KalmanSettings, status::OnlineKalmanStatus)
    if settings.store_history == true
        error("History must not be stored with an `OnlineKalmanStatus` struct. Try with an `DynamicKalmanStatus` or `SizedKalmanStatus` struct.")
    end
end

function initialise_status_history!(settings::KalmanSettings, status::DynamicKalmanStatus)
    if settings.store_history == true
        status.history_X_prior = Array{FloatVector,1}();
        status.history_X_post = Array{FloatVector,1}();
        status.history_P_prior = Array{SymMatrix,1}();
        status.history_P_post = Array{SymMatrix,1}();
        status.history_e = Array{FloatVector,1}();
        status.history_inv_F = Array{SymMatrix,1}();
        status.history_L = Array{FloatMatrix,1}();
    end
end

function initialise_status_history!(settings::KalmanSettings, status::SizedKalmanStatus)
    if settings.store_history == false
        error("History must be stored with a `SizedKalmanStatus` struct. Try with an `OnlineKalmanStatus` or `DynamicKalmanStatus` struct.")
    end
end

"""
    apriori!(old_X_prior::Nothing, settings::KalmanSettings, status::KalmanStatus)
    apriori!(old_X_prior::Nothing, settings::KalmanSettings, status::SizedKalmanStatus)

Kalman filter a-priori prediction for t==1.

# Arguments
- `old_X_prior`: latest a-priori prediction (i.e., `nothing` for the first point in time)
- `settings`: KalmanSettings struct
- `status`: KalmanStatus struct

    apriori!(old_X_prior::FloatVector, settings::KalmanSettings, status::KalmanStatus)
    apriori!(old_X_prior::FloatVector, settings::KalmanSettings, status::SizedKalmanStatus)

Kalman filter a-priori prediction.

# Arguments
- `old_X_prior`: latest a-priori prediction
- `settings`: KalmanSettings struct
- `status`: KalmanStatus struct
"""
function apriori!(old_X_prior::Nothing, settings::KalmanSettings, status::KalmanStatus)
    initialise_apriori!(settings, status);
    initialise_status_history!(settings, status);
end

function apriori!(old_X_prior::Nothing, settings::KalmanSettings, status::SizedKalmanStatus)
    initialise_apriori!(settings, status.online_status);
    initialise_status_history!(settings, status);
end

function apriori!(old_X_prior::FloatVector, settings::KalmanSettings, status::KalmanStatus)
    status.X_prior = apriori(status.X_post, settings);
    status.P_prior = apriori(status.P_post, settings);
end

apriori!(old_X_prior::FloatVector, settings::KalmanSettings, status::SizedKalmanStatus) = apriori!(old_X_prior, settings, status.online_status);

"""
    call_apriori!(settings::KalmanSettings, status::KalmanStatus)
    call_apriori!(settings::KalmanSettings, status::SizedKalmanStatus)

API to call `apriori!`.
"""
call_apriori!(settings::KalmanSettings, status::KalmanStatus) = apriori!(status.X_prior, settings, status);
call_apriori!(settings::KalmanSettings, status::SizedKalmanStatus) = apriori!(status.online_status.X_prior, settings, status);

"""
    find_observed_data(settings::KalmanSettings, status::KalmanStatus)
    find_observed_data(settings::KalmanSettings, status::SizedKalmanStatus)

Return position of the observed measurements at time status.t.

# Arguments
- `settings`: KalmanSettings struct
- `status`: KalmanStatus struct

    find_observed_data(settings::KalmanSettings, t::Int64)

Return position of the observed measurements at time t.

# Arguments
- `settings`: KalmanSettings struct
- `status`: KalmanStatus struct
"""
function find_observed_data(settings::KalmanSettings, status::KalmanStatus)
    if status.t <= settings.Y.T
        Y_t_all = @view settings.Y.data[:, status.t];
        ind_not_missings = findall(ismissing.(Y_t_all) .== false);
        if length(ind_not_missings) > 0
            return ind_not_missings;
        end
    end
end

find_observed_data(settings::KalmanSettings, status::SizedKalmanStatus) = find_observed_data(settings, status.online_status);

function find_observed_data(settings::KalmanSettings, t::Int64)
    if t <= settings.Y.T
        Y_t_all = @view settings.Y.data[:, t];
        ind_not_missings = findall(ismissing.(Y_t_all) .== false);
        if length(ind_not_missings) > 0
            return ind_not_missings;
        end
    end
end

"""
    update_P_post!(P_post_old::SymMatrix, status::KalmanStatus, K_t::FloatMatrix, R_t::SubArray{Float64})
    update_P_post!(P_post_old::Nothing, status::KalmanStatus, K_t::FloatMatrix, R_t::SubArray{Float64})

Update status.P_post.
"""
function update_P_post!(P_post_old::SymMatrix, status::KalmanStatus, K_t::FloatMatrix, R_t::SubArray{Float64})
    mul!(status.buffer_m_m, status.L, status.P_prior);
    mul!(status.P_post.data, status.buffer_m_m, status.L');
    mul!(status.buffer_m_n_obs, K_t, R_t);
    mul!(status.P_post.data, status.buffer_m_n_obs, K_t', 1.0, 1.0);
end

function update_P_post!(P_post_old::Nothing, status::KalmanStatus, K_t::FloatMatrix, R_t::SubArray{Float64})
    mul!(status.buffer_m_m, status.L, status.P_prior);
    status.P_post = Symmetric(status.buffer_m_m*status.L');
    mul!(status.buffer_m_n_obs, K_t, R_t);
    mul!(status.P_post.data, status.buffer_m_n_obs, K_t', 1.0, 1.0);
end

"""
    update_loglik!(status::KalmanStatus)

Update status.loglik.

# Arguments
- `status`: KalmanStatus struct
"""
function update_loglik!(status::KalmanStatus)
    status.loglik -= 0.5*(-logdet(status.inv_F) + status.e'*status.inv_F*status.e);
end

"""
    aposteriori!(settings::KalmanSettings, status::KalmanStatus, ind_not_missings::IntVector)
    aposteriori!(settings::KalmanSettings, status::SizedKalmanStatus, ind_not_missings::IntVector)

Kalman filter a-posteriori update. Measurements are observed (or partially observed) at time t.

The update for the covariance matrix is implemented by using the Joseph's stabilised form (Bucy and Joseph, 1968).

# Arguments
- `settings`: KalmanSettings struct
- `status`: KalmanStatus struct
- `ind_not_missings`: Position of the observed measurements

    aposteriori!(settings::KalmanSettings, status::KalmanStatus, ind_not_missings::Nothing)
    aposteriori!(settings::KalmanSettings, status::SizedKalmanStatus, ind_not_missings::Nothing)

Kalman filter a-posteriori update. All measurements are not observed at time t.

# Arguments
- `settings`: KalmanSettings struct
- `status`: KalmanStatus struct
- `ind_not_missings`: Empty array
"""
function aposteriori!(settings::KalmanSettings, status::KalmanStatus, ind_not_missings::IntVector)

    # Convenient views
    B_t = @view settings.B[ind_not_missings, :];
    R_t = @view settings.R[ind_not_missings, ind_not_missings];

    # Forecast error
    status.e = settings.Y.data[ind_not_missings, status.t];
    mul!(status.e, B_t, status.X_prior, -1.0, 1.0);

    # Convenient shortcut for the forecast error covariance matrix and Kalman gain
    # The line below initialises `status.buffer_m_n_obs` for the current point in time
    status.buffer_m_n_obs = status.P_prior*B_t';

    # Inverse of the forecast error covariance matrix
    F_t = settings.R[ind_not_missings, ind_not_missings];
    mul!(F_t, status.buffer_m_n_obs', B_t', 1.0, 1.0);
    status.inv_F = inv(Symmetric(F_t));

    # Kalman gain
    K_t = status.buffer_m_n_obs*status.inv_F;

    # Convenient shortcut for the Joseph form and needed statistics for the Kalman smoother
    status.L = Matrix(1.0I, settings.m, settings.m);
    mul!(status.L, K_t, B_t, -1.0, 1.0);

    # A posteriori estimates: X_post
    status.X_post = copy(status.X_prior);
    mul!(status.X_post, K_t, status.e, 1.0, 1.0);

    # A posteriori estimates: P_post (P_post is updated using the Joseph form)
    update_P_post!(status.P_post, status, K_t, R_t);

    # Update log likelihood
    if settings.compute_loglik == true
        update_loglik!(status);
    end
end

aposteriori!(settings::KalmanSettings, status::SizedKalmanStatus, ind_not_missings::IntVector) = aposteriori!(settings, status.online_status, ind_not_missings);

function aposteriori!(settings::KalmanSettings, status::KalmanStatus, ind_not_missings::Nothing)
    status.X_post = copy(status.X_prior);
    status.P_post = copy(status.P_prior);
    status.e = zeros(1);
    status.inv_F = Symmetric(zeros(1,1));
    status.L = Matrix(1.0I, settings.m, settings.m);
end

aposteriori!(settings::KalmanSettings, status::SizedKalmanStatus, ind_not_missings::Nothing) = aposteriori!(settings, status.online_status, ind_not_missings);

"""
    update_status_history!(settings::KalmanSettings, status::OnlineKalmanStatus)
    update_status_history!(settings::KalmanSettings, status::DynamicKalmanStatus)
    update_status_history!(settings::KalmanSettings, status::SizedKalmanStatus)

Update the `status.history_*` arrays with a new entry, when needed.
"""
update_status_history!(settings::KalmanSettings, status::OnlineKalmanStatus) = nothing;

function update_status_history!(settings::KalmanSettings, status::DynamicKalmanStatus)
    if settings.store_history == true
        push!(status.history_X_prior, status.X_prior);
        push!(status.history_X_post, status.X_post);
        push!(status.history_P_prior, status.P_prior);
        push!(status.history_P_post, status.P_post);
        push!(status.history_e, status.e);
        push!(status.history_inv_F, status.inv_F);
        push!(status.history_L, status.L);
    end
end

function update_status_history!(settings::KalmanSettings, status::SizedKalmanStatus)
    status.history_X_prior[status.online_status.t] = status.online_status.X_prior;
    status.history_X_post[status.online_status.t] = status.online_status.X_post;
    status.history_P_prior[status.online_status.t] = status.online_status.P_prior;
    status.history_P_post[status.online_status.t] = status.online_status.P_post;
    status.history_e[status.online_status.t] = status.online_status.e;
    status.history_inv_F[status.online_status.t] = status.online_status.inv_F;
    status.history_L[status.online_status.t] = status.online_status.L;
end

"""
    kfilter!(settings::KalmanSettings, status::KalmanStatus)

Kalman filter: a-priori prediction and a-posteriori update.

# Arguments
- `settings`: KalmanSettings struct
- `status`: KalmanStatus struct
"""
function kfilter!(settings::KalmanSettings, status::KalmanStatus)

    # Update status.t
    increase_time!(status);

    # A-priori prediction
    call_apriori!(settings, status);

    # Handle missing observations
    ind_not_missings = find_observed_data(settings, status);

    # Ex-post update
    aposteriori!(settings, status, ind_not_missings);

    # Update `status.history_*`
    update_status_history!(settings, status);
end

"""
    reset_kalman_status!(status::KalmanStatus)

Reset essential entries in `status` to their original state (i.e., when t==0)
"""
function reset_kalman_status!(status::KalmanStatus)
    status.t = 0;
    status.X_prior = nothing;
end

"""
    kfilter_full_sample(settings::KalmanSettings)

Run Kalman filter from t=1 to T.
"""
function kfilter_full_sample(settings::KalmanSettings, status::KalmanStatus=DynamicKalmanStatus())
    for t=1:settings.Y.T
        kfilter!(settings, status);
    end

    return status;
end

"""
    kfilter_full_sample!(settings::KalmanSettings, status::SizedKalmanStatus)

Run Kalman filter from t=1 to `history_length` and update `status` in-place.
"""
function kfilter_full_sample!(settings::KalmanSettings, status::SizedKalmanStatus)
    reset_kalman_status!(status.online_status);
    for t=1:status.history_length
        kfilter!(settings, status);
    end
end

"""
    kforecast(settings::KalmanSettings, X::Union{FloatVector, Nothing}, h::Int64)

Forecast X up to h-step ahead.

# Arguments
- `settings`: KalmanSettings struct
- `X`: State vector
- `h`: Forecast horizon

    kforecast(settings::KalmanSettings, X::Union{FloatVector, Nothing}, P::Union{SymMatrix, Nothing}, h::Int64)

Forecast X and P up to h-step ahead.

# Arguments
- `settings`: KalmanSettings struct
- `X`: State vector
- `P`: Covariance matrix of the states
- `h`: Forecast horizon
"""
function kforecast(settings::KalmanSettings, Xt::Union{FloatVector, Nothing}, h::Int64)

    # Initialise forecast history
    history_X = Array{FloatVector,1}();

    X = copy(Xt);

    # Loop over forecast horizons
    for horizon=1:h
        X = apriori(X, settings);
        push!(history_X, X);
    end

    # Return output
    return history_X;
end

function kforecast(settings::KalmanSettings, Xt::Union{FloatVector, Nothing}, Pt::Union{SymMatrix, Nothing}, h::Int64)

    # Initialise forecast history
    history_X = Array{FloatVector,1}();
    history_P = Array{SymMatrix,1}();

    X = copy(Xt);
    P = copy(Pt);

    # Loop over forecast horizons
    for horizon=1:h
        X = apriori(X, settings);
        P = apriori(P, settings);
        push!(history_X, X);
        push!(history_P, P);
    end

    # Return output
    return history_X, history_P;
end

"""
    retrieve_status_t(status::KalmanStatus)
    retrieve_status_t(status::SizedKalmanStatus)

Retrieve last filtering point in time.
"""
retrieve_status_t(status::KalmanStatus) = status.t;
retrieve_status_t(status::SizedKalmanStatus) = status.online_status.t;

"""
    update_smoothing_factors!(settings::KalmanSettings, status::KalmanStatus, ind_not_missings::IntVector, J1::FloatVector, J2::SymMatrix, e::FloatVector, inv_F::SymMatrix, L::FloatMatrix)
    update_smoothing_factors!(settings::KalmanSettings, status::SizedKalmanStatus, ind_not_missings::IntVector, J1::FloatVector, J2::SymMatrix, e::FloatVector, inv_F::SymMatrix, L::FloatMatrix)

Update J1 and J2 with a-posteriori recursion.

    update_smoothing_factors!(settings::KalmanSettings, status::KalmanStatus, ind_not_missings::Nothing, J1::FloatVector, J2::SymMatrix)
    update_smoothing_factors!(settings::KalmanSettings, status::SizedKalmanStatus, ind_not_missings::Nothing, J1::FloatVector, J2::SymMatrix)
    update_smoothing_factors!(settings::KalmanSettings, status::KalmanStatus, ind_not_missings::Nothing, J1::FloatVector, J2::SymMatrix, e::FloatVector, inv_F::SymMatrix, L::FloatMatrix)

Update J1 and J2 with a-priori recursion when all series are missing.
"""
function update_smoothing_factors!(settings::KalmanSettings, status::KalmanStatus, ind_not_missings::IntVector, J1::FloatVector, J2::SymMatrix, e::FloatVector, inv_F::SymMatrix, L::FloatMatrix)

    # Retrieve coefficients
    B_t = @view settings.B[ind_not_missings, :];
    B_inv_F = B_t'*inv_F;
    mul!(status.buffer_m_m, L', settings.C');

    # Compute J1
    mul!(status.buffer_J1, status.buffer_m_m, J1);
    mul!(J1, B_inv_F, e);
    J1 .+= status.buffer_J1;

    # Compute J2
    mul!(status.buffer_J2, status.buffer_m_m, J2);
    mul!(J2.data, status.buffer_J2, status.buffer_m_m');
    mul!(J2.data, B_inv_F, B_t, 1.0, 1.0);
end

update_smoothing_factors!(settings::KalmanSettings, status::SizedKalmanStatus, ind_not_missings::IntVector, J1::FloatVector, J2::SymMatrix, e::FloatVector, inv_F::SymMatrix, L::FloatMatrix) = update_smoothing_factors!(settings, status.online_status, ind_not_missings, J1, J2, e, inv_F, L);

function update_smoothing_factors!(settings::KalmanSettings, status::KalmanStatus, ind_not_missings::Nothing, J1::FloatVector, J2::SymMatrix)
    copyto!(J1, settings.C'*J1);
    mul!(status.buffer_J2, settings.C', J2);
    mul!(J2.data, status.buffer_J2, settings.C);
end

update_smoothing_factors!(settings::KalmanSettings, status::SizedKalmanStatus, ind_not_missings::Nothing, J1::FloatVector, J2::SymMatrix) = update_smoothing_factors!(settings, status.online_status, ind_not_missings, J1, J2);
update_smoothing_factors!(settings::KalmanSettings, status::KalmanStatus, ind_not_missings::Nothing, J1::FloatVector, J2::SymMatrix, e::FloatVector, inv_F::SymMatrix, L::FloatMatrix) = update_smoothing_factors!(settings, status, ind_not_missings, J1, J2);

"""
    backwards_pass(Xp::FloatVector, Pp::SymMatrix, J1::FloatVector)

Backward pass for the state vector.

    backwards_pass(Pp::SymMatrix, J2::SymMatrix)

Backward pass for the covariance of the states.
"""
backwards_pass(Xp::FloatVector, Pp::SymMatrix, J1::FloatVector) = Xp + Pp*J1;
backwards_pass(Pp::SymMatrix, J2::SymMatrix) = Symmetric(Pp - Pp*J2*Pp);

"""
    ksmoother(settings::KalmanSettings, status::KalmanStatus)

Kalman smoother: RTS smoother from the last evaluated time period in status to t==0.

The smoother is implemented following the approach proposed in Durbin and Koopman (2012).

# Arguments
- `settings`: KalmanSettings struct
- `status`: KalmanStatus struct
"""
function ksmoother(settings::KalmanSettings, status::KalmanStatus)

    # Initialise smoother history
    history_X = Array{FloatVector,1}();
    history_P = Array{SymMatrix,1}();

    J1 = zeros(settings.m);
    J2 = Symmetric(zeros(settings.m, settings.m));

    # Loop over t
    for t=retrieve_status_t(status):-1:1

        # Pointers
        Xp = status.history_X_prior[t];
        Pp = status.history_P_prior[t];
        e = status.history_e[t];
        inv_F = status.history_inv_F[t];
        L = status.history_L[t];

        # Handle missing observations
        ind_not_missings = find_observed_data(settings, t);

        # Smoothed estimates for t
        update_smoothing_factors!(settings, status, ind_not_missings, J1, J2, e, inv_F, L);
        pushfirst!(history_X, backwards_pass(Xp, Pp, J1));
        pushfirst!(history_P, backwards_pass(Pp, J2));
    end

    # Compute smoothed estimates for t==0
    update_smoothing_factors!(settings, status, nothing, J1, J2);
    X0 = backwards_pass(settings.X0, settings.P0, J1);
    P0 = backwards_pass(settings.P0, J2);

    # Return output
    return history_X, history_P, X0, P0;
end


text:
function assert_current_branch_is(branch_name)
    _current_branch = current_branch()
    if _current_branch != branch_name
        throw(ErrorException("Expected current branch to be \"$(branch_name)\", but the current branch is actually \"$(_current_branch)\""))
    end
    return nothing
end

function assert_file_contains(filename::AbstractString, pattern)
    filecontents = read(filename, String)
    if !(occursin(pattern, filecontents))
        throw(ErrorException("The file at \"$(filename)\" does not contain the pattern \"$(pattern)\""))
    end
    return nothing
end

function assert_string_startswith(s::AbstractString, prefix::AbstractString)
    if !(startswith(s, prefix))
        throw(ErrorException("\"$(s)\" does not start with \"$(prefix)\""))
    end
    return nothing
end


text:
@testset "Parity" begin
    ft = job_fairtensor()
    M = [true_positive_rate]
    df = disparity(m, ft; refGrp="Board")
    df = parity(df, 0.2)
    @test names(df)[3] == :true_positive_rate_parity
    arr = df[:, :true_positive_rate_parity]
    @test arr[1] && isnan(arr[2]) && !arr[3]
end


text:
import MacroTools

"""
"""
function SimplePipeline(
        objectsvector::T;
        name::S = "",
        ) where
        S <: AbstractString where
        T <: AbstractVector{F} where
        F <: AbstractFittable
    result = SimplePipeline{S,T}(name,
                                 objectsvector)
    return result
end

SimplePipeline(x::T) where T = SimplePipeline(T[x])
SimplePipeline(p::SimplePipeline) = p

Base.IteratorEltype(itertype::Type{SimplePipeline{S, T}}) where S where T =
    Base.IteratorEltype(T)
Base.IteratorSize(itertype::Type{SimplePipeline{S, T}}) where S where T =
    Base.IteratorSize(T)
Base.IndexStyle(itertype::Type{SimplePipeline{S, T}}) where S where T =
    Base.IndexStyle(T)

MacroTools.@forward SimplePipeline.objectsvector Base.axes
MacroTools.@forward SimplePipeline.objectsvector Base.eachindex
MacroTools.@forward SimplePipeline.objectsvector Base.eltype
MacroTools.@forward SimplePipeline.objectsvector Base.empty!
MacroTools.@forward SimplePipeline.objectsvector Base.findall
MacroTools.@forward SimplePipeline.objectsvector Base.findfirst
MacroTools.@forward SimplePipeline.objectsvector Base.findlast
MacroTools.@forward SimplePipeline.objectsvector Base.findnext
MacroTools.@forward SimplePipeline.objectsvector Base.findprev
MacroTools.@forward SimplePipeline.objectsvector Base.firstindex
MacroTools.@forward SimplePipeline.objectsvector Base.getindex
MacroTools.@forward SimplePipeline.objectsvector Base.isassigned
MacroTools.@forward SimplePipeline.objectsvector Base.isempty
MacroTools.@forward SimplePipeline.objectsvector Base.iterate
MacroTools.@forward SimplePipeline.objectsvector Base.lastindex
MacroTools.@forward SimplePipeline.objectsvector Base.length
MacroTools.@forward SimplePipeline.objectsvector Base.ndims
MacroTools.@forward SimplePipeline.objectsvector Base.setindex!
MacroTools.@forward SimplePipeline.objectsvector Base.size
MacroTools.@forward SimplePipeline.objectsvector Base.vec
MacroTools.@forward SimplePipeline.objectsvector Base.view

Base.:|>(x::AbstractFittable, y::AbstractFittable) = SimplePipeline(x) |>
                                                     SimplePipeline(y)

# alternatively, we could replace `Union{F1, F2}` with `typejoin(F1, F2)`
function Base.:|>(p1::SimplePipeline{S1, T1},
                  p2::SimplePipeline{S2, T2}) where
                  S1 where T1 <:AbstractVector{F1} where F1 where
                  S2 where T2 <:AbstractVector{F2} where F2
    length_1 = length(p1)
    length_2 = length(p2)
    new_objectsvector = Vector{Union{F1, F2}}(undef, length_1 + length_2)
    for i = 1:length_1
        new_objectsvector[i] = p1[i]
    end
    for i = 1:length_2
        new_objectsvector[length_1 + i] = p2[i]
    end
    result = SimplePipeline(new_objectsvector;
                            name = string(p1.name, p2.name))
    return result
end

"""
"""
function ispipeline end

ispipeline(::Any) = false
ispipeline(::AbstractPipeline) = true

"""
"""
function isflat end

isflat(x::SimplePipeline) = !any(ispipeline.(x.objectsvector))

"""
"""
function flatten(::Type{SimplePipeline}, p::SimplePipeline)
    result = _flatten(SimplePipeline, p)
    return result
end

flatten(p::SimplePipeline) = flatten(SimplePipeline, p::SimplePipeline)

# alternatively, we could replace `Union{typeof.(temp_objects)...}` with
# `typejoin(typeof.(temp_objects)...)`
function _flatten(::Type{SimplePipeline}, p::SimplePipeline)
    temp_names = Vector{Any}(undef, 0)
    temp_objects = Vector{Any}(undef, 0)
    for i = 1:length(p)
        object = p[i]
        push!(temp_names, _flatten_name(SimplePipeline, object))
        append!(temp_objects, _flatten_objects(SimplePipeline, object))
    end
    new_F = Union{typeof.(temp_objects)...}
    new_objects::Vector{new_F} = Vector{new_F}(undef,length(temp_objects))
    for j = 1:length(temp_objects)
        new_objects[j] = temp_objects[j]
    end
    new_name::String = string(p.name, join(temp_names, ""))
    new_pipeline = SimplePipeline(new_objects; name = new_name)
    return new_pipeline
end

_flatten_name(::Type{SimplePipeline}, x) = ""
_flatten_name(::Type{SimplePipeline}, p::SimplePipeline) =
    _flatten(SimplePipeline, p).name

_flatten_objects(::Type{SimplePipeline}, x) = [x]
_flatten_objects(::Type{SimplePipeline}, p::SimplePipeline) =
    _flatten(SimplePipeline, p).objectsvector

"""
"""
function set_max_epochs!(
        p::SimplePipeline,
        new_max_epochs::Integer,
        )
    for i = 1:length(p)
        set_max_epochs!(
            p[i],
            new_max_epochs,
            )
    end
    return nothing
end

"""
"""
function set_feature_contrasts!(
        p::SimplePipeline,
        feature_contrasts::AbstractFeatureContrasts,
        )
    for i = 1:length(p)
        set_feature_contrasts!(p[i], feature_contrasts)
    end
    return nothing
end

"""
"""
function get_underlying(
        p::SimplePipeline;
        saving::Bool = false,
        loading::Bool = false,
        )
    underlying = [
        get_underlying(
            o;
            saving=saving,
            loading=loading,
            ) for o in p.objectsvector
        ]
    if saving || loading
    else
        delete_nothings!(underlying)
        if length(underlying) == 0
            underlying = nothing
        elseif length(underlying) == 1
            underlying = underlying[1]
        end
    end
    return underlying
end

"""
"""
function get_history(
        p::SimplePipeline;
        saving::Bool = false,
        loading::Bool = false,
        )
    history = [
        get_history(
            o;
            saving = saving,
            loading = loading,
            ) for o in p.objectsvector
        ]
    if saving || loading
    else
        delete_nothings!(history)
        if length(history) == 0
            history = nothing
        elseif length(history) == 1
            history = history[1]
        end
    end
    return history
end

"""
"""
function parse_functions!(simplelinearpipeline::SimplePipeline)
    for i = 1:length(simplelinearpipeline)
        parse_functions!(simplelinearpipeline[i])
    end
    return nothing
end

"""
"""
function fit!(
        simplelinearpipeline::SimplePipeline,
        varargs...;
        kwargs...
        )
    output = fit!(
        simplelinearpipeline[1],
        varargs...;
        kwargs...
        )
    for i = 2:length(simplelinearpipeline)
        input = tuplify(output)
        output = fit!(
            simplelinearpipeline[i],
            input...;
            kwargs...
            )
    end
    return output
end

"""
"""
function predict(
        simplelinearpipeline::SimplePipeline,
        varargs...;
        kwargs...
        )
    output = predict(
        simplelinearpipeline[1],
        varargs...;
        kwargs...
        )
    for i = 2:length(simplelinearpipeline)
        input = tuplify(output)
        output = predict(
            simplelinearpipeline[i],
            input...;
            kwargs...
            )
    end
    return output
end

"""
"""
function predict_proba(
        simplelinearpipeline::SimplePipeline,
        varargs...;
        kwargs...
        )
    output = predict_proba(
        simplelinearpipeline[1],
        varargs...
        )
    for i = 2:length(simplelinearpipeline)
        input = tuplify(output)
        output = predict_proba(
            simplelinearpipeline[i],
            input...;
            kwargs...
            )
    end
    return output
end


text:
using AES, BenchmarkTools, Test

# AES128
k=[0x54,0x68,0x61,0x74,0x73,0x20,0x6D,0x79,0x20,0x4B,0x75,0x6E,0x67,0x20,0x46,0x75]
p=[0x54,0x77,0x6F,0x20,0x4F,0x6E,0x65,0x20,0x4E,0x69,0x6E,0x65,0x20,0x54,0x77,0x6F]
c=[0x29,0xC3,0x50,0x5F,0x57,0x14,0x20,0xF6,0x40,0x22,0x99,0xB3,0x1A,0x02,0xD7,0x3A]
key = AES128Key(k)

@test c == AES.AESEncryptBlock(p, key)
@test p == AES.AESDecryptBlock(c, key)
@test p == AES.AESDecryptBlock(AES.AESEncryptBlock(p,key),key)

# Performance Checks
cache = AES.gen_cache(key)
@btime AES.AESEncryptBlock!(p, p, key, cache)
@btime AES.AESDecryptBlock!(c, c, key, cache)

# AES192
k = [0x00, 0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08, 0x09, 0x0A, 0x0B, 0x0C, 0x0D, 0x0E, 0x0F, 0x10, 0x11, 0x12, 0x13, 0x14, 0x15, 0x16, 0x17]
p = [0x00, 0x11, 0x22, 0x33, 0x44, 0x55, 0x66, 0x77, 0x88, 0x99, 0xAA, 0xBB, 0xCC, 0xDD, 0xEE, 0xFF]
c = [0xdd, 0xa9, 0x7c, 0xa4, 0x86, 0x4c, 0xdf, 0xe0, 0x6e, 0xaf, 0x70, 0xa0, 0xec, 0x0d, 0x71, 0x91]
key = AES192Key(k)

@test c == AES.AESEncryptBlock(p, key)
@test p == AES.AESDecryptBlock(c, key)
@test p == AES.AESDecryptBlock(AES.AESEncryptBlock(p,key),key)

# Performance Checks
cache = AES.gen_cache(key)
@btime AES.AESEncryptBlock!(p, p, key, cache)
@btime AES.AESDecryptBlock!(c, c, key, cache)


text:
# echarts built-in for average, min, max
function yline!(ec::EChart, value::String; series::Int = 1, kwargs...)

	#Defensive; empty markLine already added by default to every chart
    ec.series[series].markLine == nothing ? ec.series[series].markLine = deepcopy(MarkLine()) : nothing

	#push data onto []
    push!(ec.series[series].markLine.data, Dict{Any, Any}("type" => value, "lineStyle" => LineStyle()
		)
	)

	for (k, v) in kwargs
		setfield!(ec.series[series].markLine.data[end]["lineStyle"].normal, k, v)
	end

    return ec

end

#arbitrary number line
function yline!(ec::EChart, value::Real; series::Int = 1, kwargs...)

	#Defensive; empty markLine already added by default to every chart
    ec.series[series].markLine == nothing ? ec.series[series].markLine = deepcopy(MarkLine()) : nothing

	#push data onto []
    push!(ec.series[series].markLine.data, Dict("yAxis" => value, "lineStyle" => LineStyle()
		)
	)

	for (k, v) in kwargs
		setfield!(ec.series[series].markLine.data[end]["lineStyle"], k, v)
	end

    return ec

end

#arbitrary line: need to confirm that the xAxis line can only be xAxis => value, not min/max/average
function xline!(ec::EChart, value::Any; series::Int = 1, kwargs...)

	#Defensive; empty markLine already added by default to every chart
    ec.series[series].markLine == nothing ? ec.series[series].markLine = deepcopy(MarkLine()) : nothing

	#push data onto []
    push!(ec.series[series].markLine.data, Dict("xAxis" => value, "lineStyle" => LineStyle()))

	for (k, v) in kwargs
		setfield!(ec.series[series].markLine.data[end]["lineStyle"], k, v)
	end

    return ec

end


text:
# define an abstrat tree node type - concrete types are TreeSplit and TreeLeaf
abstract type Node{T<:AbstractFloat} end

# store perf info of each variable
mutable struct SplitInfo{L, T<:AbstractFloat, S<:Int}
    gain::T
    ∑δL::SVector{L,T}
    ∑δ²L::SVector{L,T}
    ∑𝑤L::SVector{1,T}
    ∑δR::SVector{L,T}
    ∑δ²R::SVector{L,T}
    ∑𝑤R::SVector{1,T}
    gainL::T
    gainR::T
    𝑖::S
    feat::S
    cond::T
end

struct TreeNode{L, T<:AbstractFloat, S<:Integer, B<:Bool}
    left::S
    right::S
    feat::S
    cond::T
    gain::T
    pred::SVector{L,T}
    split::B
end

TreeNode(left::S, right::S, feat::S, cond::T, gain::T, L::S) where {T<:AbstractFloat, S<:Integer} = TreeNode{L,T,S,Bool}(left, right, feat, cond, gain, zeros(SVector{L,T}), true)
TreeNode(pred::SVector{L,T}) where {L,T} = TreeNode(0, 0, 0, zero(T), zero(T), pred, false)

# single tree is made of a root node that containes nested nodes and leafs
struct TrainNode{L, T<:AbstractFloat, S<:Integer}
    parent::S
    depth::S
    ∑δ::SVector{L,T}
    ∑δ²::SVector{L,T}
    ∑𝑤::SVector{1,T}
    gain::T
    𝑖::Vector{S}
    𝑗::Vector{S}
end

# single tree is made of a root node that containes nested nodes and leafs
struct Tree{L, T<:AbstractFloat, S<:Integer}
    nodes::Vector{TreeNode{L,T,S,Bool}}
end

# eval metric tracking
mutable struct Metric
    iter::Int
    metric::Float32
end
Metric() = Metric(0, Inf)

# gradient-boosted tree is formed by a vector of trees
struct GBTree{L, T<:AbstractFloat, S<:Integer}
    trees::Vector{Tree{L,T,S}}
    params::EvoTypes
    metric::Metric
    K::Int
    levels
end


text:
@testset "Natural numbers (example)" begin

# Test the natural numbers and addition
clauses = @julog [
    nat(0) <<= true,
    nat(s(N)) <<= nat(N),
    add(0, Y, Y) <<= true,
    add(s(X), Y, s(Z)) <<= add(X, Y, Z)
]

# Is 1 a natural number?
sat, subst = resolve(@julog(nat(s(0))), clauses)
@test sat == true

# Is 5 a natural number?
sat, subst = resolve(@julog(nat(s(s(s(s(s(0))))))), clauses)
@test sat == true

# Is 1 + 1 = 2?
sat, subst = resolve(@julog(add(s(0), s(0), s(s(0)))), clauses)
@test sat == true

# What are all the ways to add up to 3?
sat, subst = resolve(@julog(add(A, B, s(s(s(0))))), clauses)
subst = Set(subst)
@test @varsub({A => 0, B => s(s(s(0)))}) in subst
@test @varsub({A => s(0), B => s(s(0))}) in subst
@test @varsub({A => s(s(0)), B => s(0)}) in subst
@test @varsub({A => s(s(s(0))), B => 0}) in subst

end


text:
# type RichMethod
#     static_parameters
#     variables
#     arguments
#     macro_form
#     lowered_form
#     typed_form
# end

function slot_mapping(lam_typed)
    if isa(lam_typed.slotnames, Void) || isa(lam_typed.slottypes, Void)
        return [(SlotNumber(-1), ("", Void))]
    end
    slotnames = copy(lam_typed.slotnames)
    slottypes = copy(lam_typed.slottypes)

    slots = [(SlotNumber(i), (string(name), slottypes[i])) for (i, name) in enumerate(slotnames)]
    ssaslot = [(SSAValue(i-1), ("ssa_$(i-1)", t)) for (i,t) in enumerate(lam_typed.ssavaluetypes)]
    vcat(slots, ssaslot)
end

function get_lambda(pass, f, types)
    lambda = pass(f, types)
    if isa(lambda, Vector)
        if isempty(lambda)
            args = join(map(t->"::$t", types), ", ")
            error("$f($args) couldn't be found")
        end
        return first(lambda)
    else
        isa(lambda, LambdaInfo) && return lambda
        error("Not sure what's up with returntype of $pass. Returned: $lambda")
    end
end

function get_method(f, types)
    if !all(isleaftype, types)
        error("Not all types are concrete: $types")
    end
    if !applicable(f, types)
        error("Method $f with signature $types not defined")
    end
    precompile(f, types) # make sure there is a specialization
    x = methods(f, types)
    if length(x) != 1
        error("
            More than one method found for signature $f $types.
            Please use more specific types!
        ")
    end
    first(x)
end

function macro_form(f, types)
    m = get_method(f, types)
    file = string(m.file)
    linestart = m.line
    @show file linestart
    code, str = open(file) do io
        line = ""
        for i=1:linestart-1
            line = readline(io)
        end
        @show line
        try # lines can be one off, which will result in a parse error
            parse(line)
        catch e
            line = readline(io)
        end
        @show line
        while !eof(io)
            line = line*readline(io)
            e = Base.parse_input_line(line; filename=file)
            if !(isa(e,Expr) && e.head === :incomplete)
                return e, line
            end
        end
    end
    code, str
end

# deal with all variances in base that should really be tuples but are something else
_tuple(x) = (x,)
_tuple(x::Core.SimpleVector) = tuple(x...)
_tuple(x::Tuple) = x
_tuple{T<:Tuple}(x::Type{T}) = tuple(x.parameters...)

function get_static_parameters(f, types)
    m = get_method(f, types)
    mi = m.specializations.func
    spnames = map(x->x.name, _tuple(m.tvars))
    sptypes = _tuple(mi.sparam_vals)
    spnames, sptypes
end

function get_ast(li::LambdaInfo)
    ast = li.code
    if isa(ast, Vector{UInt8})
        return Base.uncompressed_ast(li)
    end
    ast
end

function filter_expr(keep, ast)
    replace_or_drop(x->(false, x), x->!keep(x), identity, Any[ast])[1]
end
function replace_expr(f, ast)
    replace_or_drop(f, x->false, Any[ast])[1]#ever drop
end
function replace_or_drop(f, drop, ast::Vector, result=[])
    for elem in ast
        replace_or_drop(f, drop, elem, result)
    end
    result
end

function replace_or_drop(f, drop, ast, result=[])
    drop(ast) && return result
    replace, replacement = f(ast)
    if replace
        push!(result, replacement)
    else
        expr = if isa(ast, Expr)
            nexpr = Expr(ast.head)
            @show nexpr
            replace_or_drop(f, drop, ast.args, nexpr.args)
            nexpr
        else
            ast
        end
        push!(result, expr)
    end
    result
end

function remove_static_params(ast, static_params)
    sp_names, sp_types = static_params
    replace_expr(ast) do expr
        if isa(expr, Symbol)
            idx = findfirst(sp_names, expr)
            if idx != 0
                return true, sp_types[idx]
            end
        end
        false, expr
    end
end

_typeof{T}(x::Type{T}) = Type{T}
_typeof{T}(x::T) = T
function extract_type(x::Symbol, slots)
    m = current_module()
    if isdefined(m, x)
        _typeof(getfield(m, x))
    else
        Symbol
    end
end
extract_type{T}(x::T, slots) = T
function extract_type(x::Expr, slots)
    if x.head == :(::)
        x.args[2]
    else
        extract_type(add_typing(x, slots), slots)
    end
end

function get_func(x::Expr)
    x.head == :curly && return get_func(x.args[1])
end
function get_func(x::Symbol)
    getfield(current_module(), x)
end
function get_func(x::GlobalRef)
    getfield(x.mod, x.name)
end
function extract_func(x::Expr, slots)
    @assert x.head == :call
    f = get_func(x.args[1])
    args = x.args[2:end]
    typed = add_typing(args, slots)
    _typed = !isa(typed, Vector) ? Any[typed] :typed
    f, _typed
end

function return_type(f, types)
    x = Base.return_types(f, types)
    if length(x) == 1
        x[1]
    else
        Union{x...}
    end
end

function add_typing(ast, slot_dict)
    replace_expr(ast) do x
        if isa(x, Expr) && x.head == :-> # cant type dem lambdas
            return true, x
        end
        if isa(x, Expr) && x.head == :(::) # already typed
            return true, x
        end
        if haskey(slot_dict, x)
            return true, Expr(:(::), x, slot_dict[x])
        end
        if is_call(x)
            f, typed_args = extract_func(x, slot_dict)
            types = tuple(map(x->extract_type(x, slot_dict), typed_args)...)
            T = return_type(f, types)
            x.args[2:end] = typed_args # type args
            return true, Expr(:(::), x, T)
        end
        return false, x
    end
end

is_linenumber(x) = false
is_linenumber(x::LineNumberNode) = true
function is_linenumber(x::Expr)
    x.head == :line
end

is_call(x) = false
function is_call(x::Expr)
    x.head == :call
end

using Plots; glvisualize()
using GLVisualize
plot(rand(100))
w = GLVisualize.current_screen()

function clean_form(f, types)
    ast, str = macro_form(f, types)
    static_params = get_static_parameters(f, types)
    ast = filter_expr(x->!is_linenumber(x), ast)
    remove_static_params(ast, static_params)
end
function clean_typed(f, types)
    ast = clean_form(f, types)
    slots = slot_mapping(get_lambda(code_typed, f, types))
    slot_dict = Dict()
    for (k, (name, T)) in slots
        if !isa(k, SSAValue)
            slot_dict[k] = T
        end
        slot_dict[Symbol(name)] = T
    end
    add_typing(ast, slot_dict)
end


text:
abstract type AbstractNumDisplay end

### low level methods

"""
    scan_rate(d::AbstractNumDisplay)
Getter to display actual refresh rate.

## Return

- `Int` value in Hz
"""
function scan_rate(d::AbstractNumDisplay)
    throw("Method is not applicable for $(typeof(d)) type")
end

"""
    size(d::AbstractNumDisplay)
Getter to show the maximal number of digits available for the display.
This is equal to the number of `digits_pins`.
"""
function size(d::AbstractNumDisplay)
    throw("Method is not applicable for $(typeof(d)) type")
end

"""
    shutdown_mode_off(d::AbstractNumDisplay)

Sets the normal mode, i.e. display digits.
It should be run after initialization of device.
"""
function shutdown_mode_off(d::AbstractNumDisplay)
    throw("Method is not applicable for $(typeof(d)) type")
end

"""
    shutdown_mode_on(d::AbstractNumDisplay)

In shutdown mode you can update digits values and set different settings but nothing will be shown.
Shutdown can be used to save power or as an alarm to flash the display by
successively entering and leaving shutdown mode.
"""
function shutdown_mode_on(d::AbstractNumDisplay)
    throw("Method is not applicable for $(typeof(d)) type")
end

"""
    function test_mode_off(d::AbstractNumDisplay)

Sets the normal mode, i.e. display digits.
"""
function test_mode_off(d::AbstractNumDisplay)
    throw("Method is not applicable for $(typeof(d)) type")
end

"""
    function test_mode_on(d::AbstractNumDisplay)

Display-test mode turns all LEDs on.
Test mode overrides shutdown mode. Display remain in display-test mode
until the display is reconfigured for normal operation.
"""
function test_mode_on(d::AbstractNumDisplay)
    throw("Method is not applicable for $(typeof(d)) type")
end

"""
    function set_limit(d::AbstractNumDisplay, limit::Int = size(d))

The function sets how many digits are displayed.
Since the number of scanned digits affects
the display brightness, the scan-limit register should
not be used to blank portions of the display.

The default value equal to the maximal value.

## Arguments

- `d` : device instance
- `limit` : number of active digits, from 1 to size(d)
    1 means only digit 1 is shown. 4 meand 8 digits will be shown.
"""
function set_limit(d::AbstractNumDisplay, limit::Int = size(d))
    throw("Method is not applicable for $(typeof(d)) type")
end

"""
    function decode_mode(d::AbstractNumDisplay, decode::UInt8 = 0b1111_1111)

The decode-mode register sets BCD code (0-9, E, H,
L, P, and -) or no-decode operation for each digit.

When the code B decode mode is used, the decoder
looks only at the lower nibble of the data in the digit
registers (D4–D1), disregarding bits D5–D7.
D8 sets the decimal point (DP).

When no-decode is selected, data bits D8–D1
correspond to the segment lines

all digits in decode mode => 0xff = 0b1111_1111
all digits in none-decode mode => 0x00 = 0b0000_0000
first decode other none mode => 0x01 = 0b0000_0001

## Arguments

- `d` : device instance
- `decode` : Each bit in the register corresponds to one digit. A logic high
    selects code decoding while logic low bypasses the decoder.
"""
function decode_mode(d::AbstractNumDisplay, decode::UInt8 = 0b1111_1111)
    throw("Method is not applicable for $(typeof(d)) type")
end

"""
    function set_intensity(d::AbstractNumDisplay, intensity::Int = 16)

Digital control of display brightness is provided by an
internal pulse-width modulator.

## Arguments

- `d` : device instance
- `intensity` : the value from 1 to 16 controlling the brighness.
    The modulator scales the average segment current in 16 steps.
"""
function set_intensity(d::AbstractNumDisplay, intensity::Int = 16)
    throw("Method is not applicable for $(typeof(d)) type")
end

"""
    write_digit(
        d::AbstractNumDisplay,
        value::UInt8,
        position::Int
    )

Writes a digit to the position. The result of the execution is changing one digit.

The result will depend on decode mode, see `decode_mode()` function.
The most significant bit sets the DP (dot) state.

In non-decode mode the bits 7-0 manage the sectors a, b, c, d, e, f, g.

In decode mode only bits 3-0 play a role. The sequence of numbers is the following: 0,1,2,3,4,5,6,7,8,9,-,E,H,L,P,blank

## Arguments

- `d` : object representing display device

- `value` : `UInt8` value updating the buffer value.

- `position` : number of digit to write starting from 1 which mean least significant digit.
    The maximal value depends on available digits, so it should be `<= length(indicator.digit_pins)`
"""
function write_digit(
    d::AbstractNumDisplay,
    value::UInt8,
    position::Int # starting from less significant
)
    throw("Method is not applicable for $(typeof(d)) type")
end

### high level methods

"""
    write_number(
        indicator::AbstractNumDisplay,
        number::Int,
        dp_position::Union{Int, Nothing} = nothing
    )

Writes the decimal value to the display.

## Arguments

- `indicator` : object representing display device

- `number` : decimal number. The maximal possible value depends on number of digits and will be checked.

- `dp_position` : position of digit starting from less significant digit or nothing.

## Example

```
# d is 4-digit display
write_number(d, 123)
# the result is _123

write_number(d, 1234)
# the result is 1234

write_number(d, 12345)
# throws an error
```
"""
function write_number(
    d::AbstractNumDisplay,
    number::Int,
    dp_position::Union{Int, Nothing} = nothing
)
    @assert 0 <= number < 10^size(d) "number $number cannot be displayed on $(size(d))-digits indicator"

    write_number(d, digits(number), dp_position)
end

"""
Writes several digits
Runs display in decode mode and display several digits
element in numer_vector must be Integers
"""
function write_number(
    d::AbstractNumDisplay,
    number_vector::AbstractArray{D},
    dp_position::Union{Int, Nothing} = nothing
) where D <: Union{Int, Nothing}
    l = length(number_vector)

    @assert 0 <= l <= size(d) "number_vector of length $l cannot be displayed on $(size(d))-digits indicator"
    set_limit(d, l)
    decode_mode(d)

    for i in 1:l
        num_i = number_vector[i] !== nothing ? number_vector[i] : 15
        is_dot = dp_position == i
        @assert 0 <= num_i <= 15 "number in vector must be between 0 and 15, got $num_i"

        digit_i = is_dot << 7 | num_i # dot coded by most significant bit
        write_digit(d, UInt8(digit_i), i)
    end
end

"""
works in segment-wise mode
set active segments DP A B C D E F G
"""
function write_symbols(
    d::AbstractNumDisplay,
    byte_vector::AbstractArray{UInt8},
)
    l = length(byte_vector)
    set_limit(d, l)
    decode_mode(d, 0x0)

    for i in 1:l
        write_digit(d, byte_vector[i], i)
    end
end

function write_symbols(
    d::AbstractNumDisplay,
    symbol_vector::AbstractArray{Char},
)
    # CHAR_TRANSLATOR
    byte_vector = [get(CHAR_TRANSLATOR, char, 0b00000000) for char in symbol_vector]
    write_symbols(d, byte_vector)
end

function write_symbols(
    d::AbstractNumDisplay,
    s::String
)
    symbol_vector = collect(s)
    write_symbols(d, reverse(symbol_vector))
end


text:
"""
    LazyArray(gen::Base.Generator)

Wrapper around Base.Generator that also indexes like an array. This is needed to make ComponentArrays
that hold arrays of ComponentArrays
"""
struct LazyArray{T,N,G} <: AbstractArray{T,N}
    gen::G
    LazyArray{T}(gen) where T = new{T, ndims(gen), typeof(gen)}(gen)
    LazyArray(gen::Base.Generator{A,F}) where {A,F} = new{eltype(A), ndims(gen), typeof(gen)}(gen)
end

const LazyVector{T,G} = LazyArray{T,1,G}
const LazyMatrix{T,G} = LazyArray{T,2,G}

Base.getindex(a::LazyArray, i...) =  _un_iter(getfield(a, :gen), i)

_un_iter(iter, idxs) = _un_iter(iter.f, iter.iter, idxs)
_un_iter(f, iter::Base.Generator, idxs) = f(_un_iter(iter.f, iter.iter, idxs))
_un_iter(f, iter::Base.Iterators.ProductIterator, idxs) = f(getindex.(iter.iterators, idxs))
_un_iter(f, iter, idxs) = f(iter[idxs...])

Base.getproperty(a::LazyArray, s::Symbol) = LazyArray(getproperty(item, s) for item in a)

Base.iterate(a::LazyArray) = iterate(getfield(a, :gen))
Base.iterate(a::LazyArray, state...) = iterate(getfield(a, :gen), state...)

Base.collect(a::LazyArray) = collect(getfield(a, :gen))

Base.length(a::LazyArray) = length(getfield(a, :gen))

Base.size(a::LazyArray) = size(getfield(a, :gen))

Base.eltype(::LazyArray{T,N,G}) where {T,N,G} = T

Base.show(io::IO, a::LazyArray) = show(io, collect(a))
function Base.show(io::IO, mime::MIME"text/plain", a::LazyArray)
    arr = collect(a)
    rep = repr(mime, arr)
    return print(replace(rep, r"(\d+-element )?((Vector|Array){(.+)?})" => s"\1LazyArray{\4}"; count=1))
end


text:
using Test: @test, @testset
using LogicCircuits
using ProbabilisticCircuits

@testset "Probabilistic circuits learner tests" begin
    data = dataset(twenty_datasets("nltcs"); do_shuffle=false, batch_size=-1)
    train_x = train(data)
    (pc, vtree) = learn_struct_prob_circuit(train_x; pseudocount = 1.0, algo = "chow-liu", algo_kwargs=(α=1.0, clt_root="graph_center"),
            vtree = "chow-liu", vtree_kwargs=(vtree_mode="balanced",))

    # simple test
    @test pc isa ProbΔ
    @test vtree isa PlainVtree
    @test num_variables(vtree) == num_features(data)
    @test check_parameter_integrity(pc)
    @test num_parameters(pc) == 74

    # test below has started to fail -- unclear whether that is a bug or randomness...?
    # @test pc[28].log_thetas[1] ≈ -1.1870882896239272 atol=1.0e-7

    # is structured decomposable
    for (n, vars) in variable_scopes(pc)
        @test vars == BitSet(variables(origin(n).vtree))
    end

    # all evidence sums to 1
    N = num_features(train_x);
    data_all = XData(generate_data_all(N))
    fc = FlowΔ(pc, max_batch_size(train_x), Bool, opts_accumulate_flows)
    calc_prob_all = log_likelihood_per_instance(fc, data_all)
    calc_prob_all = exp.(calc_prob_all)
    sum_prob_all = sum(calc_prob_all)
    @test sum_prob_all ≈ 1 atol = 1.0e-7;
end