text:
function interpolate_parameters_slrdamages(p, v, d, t)
    # interpolation of parameters, see notes in run_timestep for why these parameters
    if is_first(t)
        for annual_year = 2015:(gettime(t))
            yr = annual_year - 2015 + 1
            frac = annual_year - 2015
            fraction_timestep = frac / (gettime(t) - 2015)
            for r in d.region

                if use_linear || use_logburke || use_logpopulation
                    # linear
                    v.s_sealevel_ann[yr] = p.s_sealevel[t] * fraction_timestep + p.s0_initialSL * (1 - fraction_timestep)
                    v.tct_per_cap_totalcostspercap_ann[yr, r] = p.tct_per_cap_totalcostspercap[t, r] * fraction_timestep
                elseif use_logwherepossible
                    # log
                    v.s_sealevel_ann[yr] = p.s_sealevel[t]^fraction_timestep * p.s0_initialSL^(1 - fraction_timestep)
                    # linear (because tct is zero at the start of the model run)
                    v.tct_per_cap_totalcostspercap_ann[yr, r] = p.tct_per_cap_totalcostspercap[t, r] * fraction_timestep
                end

            end
        end
    else
        for annual_year = (gettime(t - 1) + 1):(gettime(t))
            yr = annual_year - 2015 + 1
            frac = annual_year - gettime(t - 1)
            fraction_timestep = frac / ((gettime(t)) - (gettime(t - 1)))

            for r in d.region
                if use_linear || use_logburke || use_logpopulation
                    # linear
                    v.s_sealevel_ann[yr] = p.s_sealevel[t] * fraction_timestep + p.s_sealevel[t - 1] * (1 - fraction_timestep)
                    v.tct_per_cap_totalcostspercap_ann[yr, r] = p.tct_per_cap_totalcostspercap[t, r] * fraction_timestep + p.tct_per_cap_totalcostspercap[t - 1, r] * (1 - fraction_timestep)
                elseif use_logwherepossible
                    # log
                    v.s_sealevel_ann[yr] = p.s_sealevel[t]^fraction_timestep * p.s_sealevel[t - 1]^(1 - fraction_timestep)
                    # linear (because tct can turn negative)
                    v.tct_per_cap_totalcostspercap_ann[yr, r] = p.tct_per_cap_totalcostspercap[t, r] * fraction_timestep + p.tct_per_cap_totalcostspercap[t - 1, r] * (1 - fraction_timestep)
                else
                    error("NO INTERPOLATION METHOD SELECTED! Specify linear or logarithmic interpolation.")
                end

            end
        end
    end
end

function calc_SLRDamages(p, v, d, t, annual_year)
    # setting the year for entry in lists.
    yr = annual_year - 2015 + 1 # + 1 because of 1-based indexing in Julia

    for r in d.region
        v.cons_percap_aftercosts_ann[yr, r] = p.cons_percap_consumption_ann[yr, r] - v.tct_per_cap_totalcostspercap_ann[yr, r] - p.act_percap_adaptationcosts_ann[yr, r]

        # apply the boundary condition if consumption minus adaptation and mitigation costs drops below the threshold
        if p.use_convergence == 1. && (v.cons_percap_aftercosts_ann[yr, r] < p.cbabsn_pcconsumptionbound_neighbourhood)
            v.cons_percap_aftercosts_ann[yr, r] = p.cbabsn_pcconsumptionbound_neighbourhood - 0.5 * p.cbaux1_pcconsumptionbound_auxiliary1 +
                                p.cbaux1_pcconsumptionbound_auxiliary1 * exp(p.cbaux2_pcconsumptionbound_auxiliary2 *
                                                                        (p.cons_percap_consumption_noconvergence_ann[yr,r] - v.tct_per_cap_totalcostspercap_ann[yr, r] - p.act_percap_adaptationcosts_ann[yr, r] - p.cbabsn_pcconsumptionbound_neighbourhood)) /
                                                                    (1 + exp(p.cbaux2_pcconsumptionbound_auxiliary2 *
                                                                            (p.cons_percap_consumption_noconvergence_ann[yr,r] - v.tct_per_cap_totalcostspercap_ann[yr, r] - p.act_percap_adaptationcosts_ann[yr, r] - p.cbabsn_pcconsumptionbound_neighbourhood)))
        elseif p.use_convergence != 1. && (v.cons_percap_aftercosts_ann[yr, r] < p.cbabs_pcconsumptionbound)
            v.cons_percap_aftercosts_ann[yr, r] = p.cbabs_pcconsumptionbound
        end
	if v.cons_percap_aftercosts_ann[yr, r] < 0
	    v.cons_percap_aftercosts_ann[yr, r] = 0
	end

        v.gdp_percap_aftercosts_ann[yr,r] = v.cons_percap_aftercosts_ann[yr, r] / (1 - p.save_savingsrate / 100)

        if (v.s_sealevel_ann[yr] - p.atl_adjustedtolerablelevelofsealevelrise_ann[yr,r]) < 0
            v.i_regionalimpactSLR_ann[yr,r] = 0
        else
            v.i_regionalimpactSLR_ann[yr,r] = v.s_sealevel_ann[yr] - p.atl_adjustedtolerablelevelofsealevelrise_ann[yr,r]
        end

        v.iref_ImpactatReferenceGDPperCapSLR_ann[yr,r] = p.wincf_weightsfactor_sea[r] * ((p.W_SatCalibrationSLR + p.iben_SLRInitialBenefit * p.scal_calibrationSLR) *
            (v.i_regionalimpactSLR_ann[yr,r] / p.scal_calibrationSLR)^p.pow_SLRImpactFxnExponent - v.i_regionalimpactSLR_ann[yr,r] * p.iben_SLRInitialBenefit)

        v.igdp_ImpactatActualGDPperCapSLR_ann[yr,r] = v.iref_ImpactatReferenceGDPperCapSLR_ann[yr,r] *
                (v.gdp_percap_aftercosts_ann[yr,r] / p.GDP_per_cap_focus_0_FocusRegionEU)^p.ipow_SLRIncomeFxnExponent

        if v.igdp_ImpactatActualGDPperCapSLR_ann[yr,r] < p.isatg_impactfxnsaturation
            v.isat_ImpactinclSaturationandAdaptationSLR_ann[yr,r] = v.igdp_ImpactatActualGDPperCapSLR_ann[yr,r]
        else
            v.isat_ImpactinclSaturationandAdaptationSLR_ann[yr,r] = p.isatg_impactfxnsaturation +
                ((100 - p.save_savingsrate) - p.isatg_impactfxnsaturation) *
                ((v.igdp_ImpactatActualGDPperCapSLR_ann[yr,r] - p.isatg_impactfxnsaturation) /
                (((100 - p.save_savingsrate) - p.isatg_impactfxnsaturation) +
                (v.igdp_ImpactatActualGDPperCapSLR_ann[yr,r] - p.isatg_impactfxnsaturation)))
        end
        if v.i_regionalimpactSLR_ann[yr,r] < p.impmax_maxSLRforadaptpolicySLR[r]
            v.isat_ImpactinclSaturationandAdaptationSLR_ann[yr,r] = v.isat_ImpactinclSaturationandAdaptationSLR_ann[yr,r] * (1 - p.imp_actualreductionSLR_ann[yr,r] / 100)
        else
            v.isat_ImpactinclSaturationandAdaptationSLR_ann[yr,r] = v.isat_ImpactinclSaturationandAdaptationSLR_ann[yr,r] * (1 - (p.imp_actualreductionSLR_ann[yr,r] / 100) * p.impmax_maxSLRforadaptpolicySLR[r] /
                v.i_regionalimpactSLR_ann[yr,r])
        end

        v.isat_per_cap_SLRImpactperCapinclSaturationandAdaptation_ann[yr,r] = (v.isat_ImpactinclSaturationandAdaptationSLR_ann[yr,r] / 100) * v.gdp_percap_aftercosts_ann[yr,r]
        v.rcons_per_cap_SLRRemainConsumption_ann[yr,r] = v.cons_percap_aftercosts_ann[yr,r] - v.isat_per_cap_SLRImpactperCapinclSaturationandAdaptation_ann[yr,r]
        v.rgdp_per_cap_SLRRemainGDP_ann[yr,r] = v.rcons_per_cap_SLRRemainConsumption_ann[yr,r] / (1 - p.save_savingsrate / 100)

    end

end

@defcomp SLRDamages begin
    region = Index()
    year = Index()

    y_year = Parameter(index=[time], unit="year")
    y_year_ann = Parameter(index=[year], unit="year")
    y_year_0 = Parameter(unit="year")

    # incoming parameters from SeaLevelRise
    s_sealevel = Parameter(index=[time], unit="m")
    s0_initialSL = Parameter(unit="m", default=0.18999999999999997)               # mode initial sea level, from PAGE-ICE
    s_sealevel_ann = Variable(index=[year], unit="m")                                                                                # interpolate (SeaLevelRise)
    # incoming parameters to calculate consumption per capita after Costs
    cons_percap_consumption = Parameter(index=[time, region], unit="\$/person")
    cons_percap_consumption_ann = Parameter(index=[year, region], unit="\$/person")
    cons_percap_consumption_0 = Parameter(index=[region], unit="\$/person")
    tct_per_cap_totalcostspercap = Parameter(index=[time,region], unit="\$/person")
    tct_per_cap_totalcostspercap_ann = Variable(index=[year,region], unit="\$/person")                                                 # interpolate (TotalAbatement)
    act_percap_adaptationcosts = Parameter(index=[time, region], unit="\$/person")
    act_percap_adaptationcosts_ann = Parameter(index=[year, region], unit="\$/person")

    # component parameters
    impmax_maxSLRforadaptpolicySLR = Parameter(index=[region], unit="m")

    save_savingsrate = Parameter(unit="%", default=15.00) # pp33 PAGE09 documentation, "savings rate".
    wincf_weightsfactor_sea = Parameter(index=[region], unit="")
    W_SatCalibrationSLR = Parameter(default=1.0) # pp33 PAGE09 documentation, "Sea level impact at calibration sea level rise"
    ipow_SLRIncomeFxnExponent = Parameter(default=-0.30)
    pow_SLRImpactFxnExponent = Parameter(default=0.7333333333333334)
    iben_SLRInitialBenefit = Parameter(default=0.00)
    scal_calibrationSLR = Parameter(default=0.5)
    GDP_per_cap_focus_0_FocusRegionEU = Parameter(unit="\$/person", default=34298.93698672955)

    # component variables
    cons_percap_aftercosts = Variable(index=[time, region], unit="\$/person")
    cons_percap_aftercosts_ann = Variable(index=[year, region], unit="\$/person")
    gdp_percap_aftercosts = Variable(index=[time, region], unit="\$/person")
    gdp_percap_aftercosts_ann = Variable(index=[year, region], unit="\$/person")

    atl_adjustedtolerablelevelofsealevelrise = Parameter(index=[time,region], unit="m") # meter
    atl_adjustedtolerablelevelofsealevelrise_ann = Parameter(index=[year,region], unit="m") # meter
    imp_actualreductionSLR = Parameter(index=[time, region], unit="%")
    imp_actualreductionSLR_ann = Parameter(index=[year, region], unit="%")
    i_regionalimpactSLR = Variable(index=[time, region], unit="m")
    i_regionalimpactSLR_ann = Variable(index=[year, region], unit="m")

    iref_ImpactatReferenceGDPperCapSLR = Variable(index=[time, region])
    iref_ImpactatReferenceGDPperCapSLR_ann = Variable(index=[year, region])
    igdp_ImpactatActualGDPperCapSLR = Variable(index=[time, region])
    igdp_ImpactatActualGDPperCapSLR_ann = Variable(index=[year, region])

    isatg_impactfxnsaturation = Parameter(unit="unitless")
    isat_ImpactinclSaturationandAdaptationSLR = Variable(index=[time,region])
    isat_ImpactinclSaturationandAdaptationSLR_ann = Variable(index=[year,region])
    isat_per_cap_SLRImpactperCapinclSaturationandAdaptation = Variable(index=[time, region], unit="\$/person")
    isat_per_cap_SLRImpactperCapinclSaturationandAdaptation_ann = Variable(index=[year, region], unit="\$/person")

    rcons_per_cap_SLRRemainConsumption = Variable(index=[time, region], unit="\$/person") # include?
    rcons_per_cap_SLRRemainConsumption_ann = Variable(index=[year, region], unit="\$/person") # include?
    rgdp_per_cap_SLRRemainGDP = Variable(index=[time, region], unit="\$/person")
    rgdp_per_cap_SLRRemainGDP_ann = Variable(index=[year, region], unit="\$/person")

    ###############################################
    # Growth Effects - additional variables and parameters
    ###############################################
    # new parameters for convergence boundary system
    use_convergence =                           Parameter(unit="none", default=1.)
    cbabs_pcconsumptionbound =                  Parameter(unit="\$/person", default=740.65)
    cbabsn_pcconsumptionbound_neighbourhood =   Parameter(unit="\$/person")
    cbaux1_pcconsumptionbound_auxiliary1 =      Parameter(unit="none")
    cbaux2_pcconsumptionbound_auxiliary2 =      Parameter(unit="none")
    cons_percap_consumption_noconvergence =     Parameter(index=[time, region], unit="\$/person")
    cons_percap_consumption_noconvergence_ann = Parameter(index=[year, region], unit="\$/person")
    ###############################################

    function run_timestep(p, v, d, t)

        interpolate_parameters_slrdamages(p, v, d, t)

        for r in d.region
            v.cons_percap_aftercosts[t, r] = p.cons_percap_consumption[t, r] - p.tct_per_cap_totalcostspercap[t, r] - p.act_percap_adaptationcosts[t, r]

            # apply the boundary condition if consumption minus adaptation and mitigation costs drops below the threshold
            if p.use_convergence == 1. && (v.cons_percap_aftercosts[t, r] < p.cbabsn_pcconsumptionbound_neighbourhood)
                v.cons_percap_aftercosts[t, r] = p.cbabsn_pcconsumptionbound_neighbourhood - 0.5 * p.cbaux1_pcconsumptionbound_auxiliary1 +
                                    p.cbaux1_pcconsumptionbound_auxiliary1 * exp(p.cbaux2_pcconsumptionbound_auxiliary2 *
                                                                            (p.cons_percap_consumption_noconvergence[t,r] - p.tct_per_cap_totalcostspercap[t, r] - p.act_percap_adaptationcosts[t, r] - p.cbabsn_pcconsumptionbound_neighbourhood)) /
                                                                        (1 + exp(p.cbaux2_pcconsumptionbound_auxiliary2 *
                                                                                (p.cons_percap_consumption_noconvergence[t,r] - p.tct_per_cap_totalcostspercap[t, r] - p.act_percap_adaptationcosts[t, r] - p.cbabsn_pcconsumptionbound_neighbourhood)))
            elseif p.use_convergence != 1. && (v.cons_percap_aftercosts[t, r] < p.cbabs_pcconsumptionbound)
                v.cons_percap_aftercosts[t, r] = p.cbabs_pcconsumptionbound
            end

            v.gdp_percap_aftercosts[t,r] = v.cons_percap_aftercosts[t, r] / (1 - p.save_savingsrate / 100)

            if (p.s_sealevel[t] - p.atl_adjustedtolerablelevelofsealevelrise[t,r]) < 0
                v.i_regionalimpactSLR[t,r] = 0
            else
                v.i_regionalimpactSLR[t,r] = p.s_sealevel[t] - p.atl_adjustedtolerablelevelofsealevelrise[t,r]
            end

            v.iref_ImpactatReferenceGDPperCapSLR[t,r] = p.wincf_weightsfactor_sea[r] * ((p.W_SatCalibrationSLR + p.iben_SLRInitialBenefit * p.scal_calibrationSLR) *
                (v.i_regionalimpactSLR[t,r] / p.scal_calibrationSLR)^p.pow_SLRImpactFxnExponent - v.i_regionalimpactSLR[t,r] * p.iben_SLRInitialBenefit)

            v.igdp_ImpactatActualGDPperCapSLR[t,r] = v.iref_ImpactatReferenceGDPperCapSLR[t,r] *
                    (v.gdp_percap_aftercosts[t,r] / p.GDP_per_cap_focus_0_FocusRegionEU)^p.ipow_SLRIncomeFxnExponent

            if v.igdp_ImpactatActualGDPperCapSLR[t,r] < p.isatg_impactfxnsaturation
                v.isat_ImpactinclSaturationandAdaptationSLR[t,r] = v.igdp_ImpactatActualGDPperCapSLR[t,r]
            else
                v.isat_ImpactinclSaturationandAdaptationSLR[t,r] = p.isatg_impactfxnsaturation +
                    ((100 - p.save_savingsrate) - p.isatg_impactfxnsaturation) *
                    ((v.igdp_ImpactatActualGDPperCapSLR[t,r] - p.isatg_impactfxnsaturation) /
                    (((100 - p.save_savingsrate) - p.isatg_impactfxnsaturation) +
                    (v.igdp_ImpactatActualGDPperCapSLR[t,r] - p.isatg_impactfxnsaturation)))
            end
            if v.i_regionalimpactSLR[t,r] < p.impmax_maxSLRforadaptpolicySLR[r]
                v.isat_ImpactinclSaturationandAdaptationSLR[t,r] = v.isat_ImpactinclSaturationandAdaptationSLR[t,r] * (1 - p.imp_actualreductionSLR[t,r] / 100)
            else
                v.isat_ImpactinclSaturationandAdaptationSLR[t,r] = v.isat_ImpactinclSaturationandAdaptationSLR[t,r] * (1 - (p.imp_actualreductionSLR[t,r] / 100) * p.impmax_maxSLRforadaptpolicySLR[r] /
                    v.i_regionalimpactSLR[t,r])
            end

            v.isat_per_cap_SLRImpactperCapinclSaturationandAdaptation[t,r] = (v.isat_ImpactinclSaturationandAdaptationSLR[t,r] / 100) * v.gdp_percap_aftercosts[t,r]
            v.rcons_per_cap_SLRRemainConsumption[t,r] = v.cons_percap_aftercosts[t,r] - v.isat_per_cap_SLRImpactperCapinclSaturationandAdaptation[t,r]
            v.rgdp_per_cap_SLRRemainGDP[t,r] = v.rcons_per_cap_SLRRemainConsumption[t,r] / (1 - p.save_savingsrate / 100)

        end

        # calculate  for this specific year
        if is_first(t)
            for annual_year = 2015:(gettime(t))
                calc_SLRDamages(p, v, d, t, annual_year)
            end
        else
            for annual_year = (gettime(t - 1) + 1):(gettime(t))
                calc_SLRDamages(p, v, d, t, annual_year)
            end
        end

    end
end

# Still need this function in order to set the parameters than depend on
# readpagedata, which takes model as an input. These cannot be set using
# the default keyword arg for now.
function addslrdamages(model::Model)
    SLRDamagescomp = add_comp!(model, SLRDamages)

    SLRDamagescomp[:wincf_weightsfactor_sea] = readpagedata(model, "data/wincf_weightsfactor_sea.csv")
    SLRDamagescomp[:impmax_maxSLRforadaptpolicySLR] = readpagedata(model, "data/impmax_sealevel.csv")

    return SLRDamagescomp
end


text:
# Fallback functions for conjugates

posterior_canon(pri, G::IncompleteFormulation, x) = posterior_canon(pri, suffstats(G, x))
posterior_canon(pri, G::IncompleteFormulation, x, w) = posterior_canon(pri, suffstats(G, x, w))

posterior{P<:Distribution}(pri::P, ss::SufficientStats) = convert(P, posterior_canon(pri, ss))
posterior{P<:Distribution}(pri::P, G::IncompleteFormulation, x) = convert(P, posterior_canon(pri, G, x))
posterior{P<:Distribution}(pri::P, G::IncompleteFormulation, x, w) = convert(P, posterior_canon(pri, G, x, w))

posterior_rand(pri, ss::SufficientStats) = rand(posterior_canon(pri, ss))
posterior_rand(pri, G::IncompleteFormulation, x) = rand(posterior_canon(pri, G, x))
posterior_rand(pri, G::IncompleteFormulation, x, w) = rand(posterior_canon(pri, G, x, w))

posterior_rand!(r::Array, pri, ss::SufficientStats) = rand!(posterior_canon(pri, ss), r)
posterior_rand!(r::Array, pri, G::IncompleteFormulation, x) = rand!(posterior_canon(pri, G, x), r)
posterior_rand!(r::Array, pri, G::IncompleteFormulation, x, w) = rand!(posterior_canon(pri, G, x, w), r)

posterior_mode(pri, ss::SufficientStats) = mode(posterior_canon(pri, ss))
posterior_mode(pri, G::IncompleteFormulation, x) = mode(posterior_canon(pri, G, x))
posterior_mode(pri, G::IncompleteFormulation, x, w) = mode(posterior_canon(pri, G, x, w))

fit_map(pri, G::IncompleteFormulation, x) = complete(G, pri, posterior_mode(pri, G, x))
fit_map(pri, G::IncompleteFormulation, x, w) = complete(G, pri, posterior_mode(pri, G, x, w))

posterior_randmodel(pri, G::IncompleteFormulation, x) = complete(G, pri, posterior_rand(pri, G, x))
posterior_randmodel(pri, G::IncompleteFormulation, x, w) = complete(G, pri, posterior_rand(pri, G, x, w))

ext: jl


text:
using Metal

@show devices()
dev = MtlDevice(1)

#lib = MTL.LibraryWithFile(d, "default.metallib")

src = """
#include <metal_stdlib>
using namespace metal;

/// This is a Metal Shading Language (MSL) function equivalent to the add_arrays() C function, used to perform the calculation on a GPU.
kernel void add_arrays(device const float* inA,
                       device const float* inB,
                       device float* result,
                       uint index [[thread_position_in_grid]])
{
    // the for-loop is replaced with a collection of threads, each of which
    // calls this function.
    result[index] = inA[index] + inB[index];
}
"""

bufferSize = 128
bufferA = MtlBuffer(Float32, dev, bufferSize, storage=Shared)
bufferB = MtlBuffer(Float32, dev, bufferSize, storage=Shared)
bufferC = MtlBuffer(Float32, dev, bufferSize, storage=Shared)

#vecA = unsafe_wrap(Vector{Float32}, convert(Ptr{Float32}, content(bufferA)), bufferSize)
vecA = unsafe_wrap(Vector{Float32}, bufferA, tuple(bufferSize))
vecB = unsafe_wrap(Vector{Float32}, bufferB, tuple(bufferSize))
vecC = unsafe_wrap(Vector{Float32}, bufferC, tuple(bufferSize))

using Random
rand!.([vecA, vecB])

## Setup
opts = MtlCompileOptions()
lib = MtlLibrary(dev, src, opts)

fun = MtlFunction(lib, "add_arrays")
pip_addfun = MtlComputePipelineState(dev, fun)
queue = global_queue(dev) #MtlCommandQueue(dev)

##
cmd = MTL.commit!(queue) do cmdbuf
    MtlComputeCommandEncoder(cmdbuf) do enc
        MTL.set_function!(enc, pip_addfun)
        MTL.set_buffers!(enc,
                                [bufferA, bufferB, bufferC],
                                [0,0,0], 1:3)
        gridSize = MtSize(length(vecA), 1, 1)
        threadGroupSize = min(length(vecA), pip_addfun.maxTotalThreadsPerThreadgroup)
        threadGroupSize = MTL.MtSize(threadGroupSize, 1, 1)
        MTL.append_current_function!(enc, gridSize, threadGroupSize)
    end
end

# Execute
wait(cmd)

@show vecC


text:
struct TPUSession
    sess::TensorFlow.Session
    job_name::String
    tpu_devices::Vector{DeviceIndex}
    topology::TopologyProto
end
function Base.show(io::IO, sess::TPUSession)
    print(io, "TPUSession(<$(length(sess.tpu_devices)) TPU chips in $(topo_string(sess.topology)) topology>)")
end

Base.run(sess::TPUSession, args...; kwargs...) = Base.run(tf_session(sess), args...; kwargs...)

const tf = TensorFlow
function TPUSession(grpc_endpoints::Vector{String})
    if length(grpc_endpoints) == 1
        job_name = "tpu_worker"
        sess = Session(Graph(); target="grpc://$(grpc_endpoints[1])")
    else
        job_name = "tpu_job"
        job = tf.tensorflow.JobDef(
            name = job_name,
            tasks = Dict(x[1]-1=>x[2] for x in pairs(grpc_endpoints)))

        cluster = tf.tensorflow.ClusterDef(job=[job])

        config = tf.tensorflow.ConfigProto()
        config.cluster_def = cluster

        sess = Session(Graph(), config; target="grpc://$(grpc_endpoints[1])")
    end
    devices = collect(TensorFlow.DeviceList(sess))
    tpu_devices = DeviceIndex[]
    for d in devices
        if d.device_type == "TPU"
            dev = TensorFlow.Device(d.name)
            push!(tpu_devices, DeviceIndex(dev.parts[3].index, dev.parts[4].index))
        end
    end
    topology = initialize_tpu!(sess)
    TPUSession(sess, job_name, tpu_devices, topology)
end
TPUSession(endpoint::String) = TPUSession([endpoint])

all_tpu_devices(sess::TPUSession) = sess.tpu_devices

function make_xrt_execute_on(device::DeviceIndex, com::XRTCompilation, inputs...)
    make_run_op(com, inputs...; device=tf_tpu_device(com.sess, device))
end

function make_infeed_on(sess, device::DeviceIndex, args...)
    _make_infeed_op(sess,
        args...;
        device = tf_host_cpu_device(sess, device),
        device_ordinal = device.index)
end

function make_outfeed_on(sess, device::DeviceIndex, args...)
    make_outfeed_op(sess,
        args...;
        device = tf_host_cpu_device(sess, device),
        device_ordinal = device.index)
end

tf_session(sess) = sess
tf_session(sess::TPUSession) = sess.sess
tf_graph(sess) = tf_session(sess).graph

job_name(sess::TPUSession) = sess.job_name

# Global session management
let sess = nothing
    global global_session, initialize!, set_global_session!
    function try_initialize_from_env!()
        if haskey(ENV, "COLAB_TPU_ADDR")
            set_global_session!(TPUSession(ENV["COLAB_TPU_ADDR"]))
        elseif haskey(ENV, "COLAB_GPU")
            @warn "Running on Colab, but TPU address is not set. Did you select the TPU accelerator option?"
        end
    end
    function global_session()
        if sess === nothing
            try_initialize_from_env!()
        end
        sess === nothing && error("Global TPU session not initialized. Call XLA.initialize!(endpoint) or pass an active session directly.")
        sess
    end
    function set_global_session!(new_sess::TPUSession)
        sess = new_sess
    end
    function initialize!(endpoint; reset=false)
        (!reset && sess !== nothing) &&
            error("Cowardly refusing to override existing session. Use reset=true to force.")
        set_global_session!(TPUSession(endpoint))
    end
end


text:
using Documenter, SpreadboostDummy

makedocs(
	format = :html,
	sitename = "Distributed Boosting",
	pages = ["index.md",
			 "installation.md",
			 "usage.md",
			 "binaries.md"]
	)
deploydocs(
	repo= "",
	target = "build",
	deps = nothing,
	make = nothing
	)

max_stars_repo_name: danielazoeller/DistributedBoosting

max_issues_repo_name: danielazoeller/DistributedBoosting

max_forks_repo_name: danielazoeller/DistributedBoosting


text:
# generate LLVM wrappers

using Clang

function wrap(config, destdir)
    @info("Wrapping LLVM C API at $destdir")
    if !isdir(destdir)
        mkdir(destdir)
    end

    includedir = readchomp(`$config --includedir`)
    cppflags = split(readchomp(`$config --cppflags`))

    # TODO: this should really be taken care of by Clang.jl...
    if haskey(ENV, "LLVM_CONFIG")
        clang_config = ENV["LLVM_CONFIG"]
    else
        clang_config = joinpath(dirname(pathof(Clang)), "..", "deps", "usr", "tools", "llvm-config")
    end
    clang_libdir = readchomp(`$clang_config --libdir`)
    clang_version = readchomp(`$clang_config --version`)

    # Set-up arguments to clang
    clang_includes  = map(x->x[3:end], filter( x->startswith(x,"-I"), cppflags))
    push!(clang_includes, "$clang_libdir/clang/$clang_version/include")
    clang_extraargs =                  filter(x->!startswith(x,"-I"), cppflags)

    # Recursively discover LLVM C API headers (files ending in .h)
    header_dirs = String[joinpath(includedir, "llvm-c")]
    header_files = String[]
    while !isempty(header_dirs)
        parent = pop!(header_dirs)
        children = readdir(parent)
        for child in children
            path = joinpath(parent, child)
            if isdir(path)
                push!(header_dirs, path)
            elseif isfile(path) && endswith(path, ".h")
                push!(header_files, path)
            end
        end
    end

    context = init(;
                    headers = header_files,
                    output_file = "$destdir/libLLVM_h.jl",
                    common_file = "$destdir/libLLVM_common.jl",
                    clang_includes = convert(Vector{String}, clang_includes),
                    clang_args = convert(Vector{String}, clang_extraargs),
                    header_library = x->"libllvm",
                    header_wrapped = (top,cursor)->occursin("include/llvm", cursor)
                  )

    run(context)
end

length(ARGS) == 2 || error("Usage: wrap.jl /path/to/llvm-config target")
config = ARGS[1]
ispath(config) || error("llvm-config at $config is't a valid path")

# Manual clean-up:
# - remove build-host details (HAVE_INTTYPES_H, LLVM_DEFAULT_TARGET_TRIPLE etc) in libLLVM_common.jl
# - remove LLVMInitializeAll and LLVMInitializeNative wrappers (these are macros)
# - remove "# Skipping ..." comments by Clang.jl
# - replace `ccall\(\((:.+), libllvm\), (.*)\)` with `@apicall($1, $2)`
# - replace `const (LLVMOpaque.*) = Cvoid` with `struct $1 end`
# - use `gawk -i inplace '/^[[:blank:]]*$/ { print; next; }; {cur = seen[$0]; if(!seen[$0]++ || (/^end$/ && !prev) || /^.*Clang.*$/) print $0; prev=cur}' libLLVM_h.jl` to remove duplicates
# - use `cat -s` to remove duplicate empty lines

target = ARGS[2]
wrap(config, target)


text:
abstract type AbstractGene end

"""
Struct for storing information on genomic locations. `strand` can be '+', '-',
or '.' when the strand is irrelevant. `order` is used to store discontiguous
sequences, indicated in the GenBank file with the order() and join() operators.
"""
struct Locus
    position::UnitRange{Int}
    strand::Char
    complete_left::Bool
    complete_right::Bool
    order::Vector{UnitRange{Int}}
    join::Bool
end

"""
    Locus()
    Locus(position::UnitRange{Int})
    Locus(position::UnitRange{Int}, strand::Char)

"""
Locus() = Locus(1:1, '.', true, true, UnitRange{Int}[], false)
Locus(position::UnitRange{Int}) = Locus(position, '.', true, true, UnitRange{Int}[], false)
Locus(position::UnitRange{Int}, strand::Char) = Locus(position, strand, true, true, UnitRange{Int}[], false)
Locus(position::UnitRange{Int}, strand::Char, complete_left, complete_right) = Locus(position, strand, complete_left, complete_right, UnitRange{Int}[], false)

Base.convert(::Type{Locus}, x::UnitRange{Int}) = Locus(x)
function Base.convert(::Type{Locus}, x::StepRange{Int, Int})
    if x.step == -1
        return Locus(x.stop:x.start, '-')
    elseif x.step == 1
        return Locus(x.start:x.stop)
    end
    throw(DomainError(x, "`x` must have a step of 1 or -1"))
end

"""
Struct for storing annotations for a single chromosome, plasmid, contig, etc.
Contains five fields: `name`, `sequence`, `header`, `genes`, and `genedata`.
Annotations are stored as a `DataFrame` in `genedata`, but can be accessed
more easily through `genes` using the API provided in this module.
"""
mutable struct Chromosome{G <: AbstractGene}
    name::String
    sequence::LongDNASeq
    header::String
    genes::Vector{G}
    genedata::DataFrame
    function Chromosome{G}() where G
        new("", dna"", "", G[], DataFrame(feature = String[], locus = Locus[]))
    end
end

Chromosome(args...) = Chromosome{Gene}(args...)

struct Gene <: AbstractGene
    index::UInt
    parent::Chromosome
end

"""
    addgene!(chr::Chromosome, feature, locus; kw...)

Add gene to `chr`. `locus` can be a `Locus`, a UnitRange, or a StepRange (for
decreasing ranges, which will be annotated on the complementary strand).

# Example
```julia
addgene!(chr, "CDS", 1:756;
    locus_tag = "gene0001",
    product = "Chromosomal replication initiator protein dnaA")
```
"""
function addgene!(chr::Chromosome, feature, locus; kw...)
    locus = convert(Locus, locus)
    push!(chr.genedata, vcat([feature, locus], fill(missing, size(chr.genedata, 2) - 2)))
    index = UInt32(length(chr.genes) + 1)
    gene = Gene(index, chr)
    for (k, v) in kw
        Base.setproperty!(gene, k, v)
    end
    push!(chr.genes, gene)
    return gene
end

"""
    delete!(gene::Gene)

Delete `gene` from `gene.parent`. Warning: does not work when broadcasted! Use
`delete!(::AbstractVector{Gene}) instead.`
"""
function Base.delete!(gene::Gene)
    i = gene.index
    deleterows!(gene.parent.genedata, i)
    deleteat!(gene.parent.genes, lastindex(gene.parent.genes))
    nothing
end

"""
    delete!(genes::AbstractArray{Gene, 1})

Delete all genes in `genes` from `genes[1].parent`.

# Example
```julia
delete!(@genes(chr, length(gene) <= 60))
```
"""
function Base.delete!(genes::AbstractArray{Gene, 1})
    indices = genes.index
    DataFrames.deleterows!(genes.parent.genedata, Int.(indices))
    lastindices = length(genes.parent.genes) - length(indices) + 1 : length(genes.parent.genes)
    deleteat!(genes.parent.genes, lastindices)
    nothing
end

function Base.propertynames(gene::AbstractGene)
    names(gene.parent.genedata)
end

"""
    vectorise(A::AbstractArray{Union{Missing, T}, 1}) where T

Convert an array of type Vector{Union{Missing, T}} to Vector{Union{Missing, Vector{T}}}.
"""
function vectorise(A::AbstractArray{Union{Missing, T}, 1}) where T
    B = Union{Missing, Vector{T}}[]
    for a in A
        if ismissing(a)
            push!(B, missing)
        else
            push!(B, [a])
        end
    end
    return B
end

"""
    pushproperty!(gene::AbstractGene, qualifier::Symbol, value::T)

Add a property to `gene`, similarly to `Base.setproperty!(::gene)`, but if the
property is not missing in `gene`, it will be transformed to store a vector
instead of overwriting existing data.

```julia
julia> eltype(chr.genedata[!, :EC_number])
Union{Missing,String}

julia> chr.genes[1].EC_number = "EC:1.2.3.4"
"EC:1.2.3.4"

julia> pushproperty!(chr.genes[1], :EC_number, "EC:4.3.2.1"); chr.genes[1].EC_number
2-element Array{String,1}:
 "EC:1.2.3.4"
 "EC:4.3.2.1"

julia> eltype(chr.genedata[!, :EC_number])
Union{Missing, Array{String,1}}
```
"""
function pushproperty!(gene::AbstractGene, qualifier::Symbol, value::T; forceany = true) where T
    if hasproperty(gene.parent.genedata, qualifier)
        C = eltype(gene.parent.genedata[!, qualifier])
        if T <: C
            if ismissing(gene.parent.genedata[gene.index, qualifier])
                gene.parent.genedata[gene.index, qualifier] = value
            else
                gene.parent.genedata[!, qualifier] = vectorise(gene.parent.genedata[!, qualifier])
                push!(gene.parent.genedata[gene.index, qualifier], value)
            end
        elseif Vector{T} <: C
            if ismissing(gene.parent.genedata[gene.index, qualifier])
                gene.parent.genedata[gene.index, qualifier] = [value]
            else
                push!(gene.parent.genedata[gene.index, qualifier], value)
            end
        elseif forceany && !(C <: AbstractVector)
            if ismissing(gene.parent.genedata[gene.index, qualifier])
                gene.parent.genedata[!, qualifier] = convert(Vector{Any}, gene.parent.genedata[!, qualifier])
                gene.parent.genedata[gene.index, qualifier] = value
            else
                gene.parent.genedata[!, qualifier] = vectorise(convert(Vector{Any}, gene.parent.genedata[! ,qualifier]))
                push!(gene.parent.genedata[gene.index, qualifier], value)
            end
        elseif forceany && C <: AbstractVector
            if ismissing(gene.parent.genedata[gene.index, qualifier])
                gene.parent.genedata[!, qualifier] = convert(Vector{Any}, gene.parent.genedata[!, qualifier])
                gene.parent.genedata[gene.index, qualifier] = [value]
            else
                @error "This shouldn't happen"
            end
        else
            @error "Tried to add a '$T' to '$qualifier::$(typeof(gene.parent.genedata[!, qualifier]))'"
        end
    else
        s = size(gene.parent.genedata, 1)
        gene.parent.genedata[!, qualifier] = Vector{Union{Missing, T}}(missing, s)
        gene.parent.genedata[gene.index, qualifier] = value
    end
    return value
end

"""
    get(g::Gene, key, default)

Retrieve `key` from `g`. If `key` is missing, return `default`.
"""
function Base.get(g::Gene, key, default)
    p = getproperty(g, key)
    if ismissing(p)
        return default
    else
        return p
    end
end

"""
    sequence(gene::AbstractGene)

Return genomic sequence for `gene`.

# Example
The following code will write the translated sequences of all genes in `chr` to a file:
```julia
using FASTX
writer = FASTA.Writer(open("proteins.fasta", "w"))
for gene in @genes(chr, iscds)
    aaseq = translate(sequence(gene))
    write(writer, FASTA.record(gene.locus_tag, get(gene, :product, ""), aaseq))
end
close(writer)
```
"""
function sequence(gene::AbstractGene)
    ifelse(gene.locus.strand == '-',
        reverse_complement(gene.parent.sequence[gene.parent.genedata[gene.index, :locus].position]),
        gene.parent.sequence[gene.parent.genedata[gene.index, :locus].position])
end

function Base.length(gene::AbstractGene)
    length(gene.parent.genedata[gene.index, :locus].position)
end

"""
    iscomplement(gene::Abstract)

Return `true` if `gene.locus.compliment == '-'`, otherwise return `false`.
"""
iscomplement(gene::AbstractGene) = gene.parent.genedata[gene.index, :locus].strand == '-'

function appendstring(field, v)
    s = ""
    if v isa Bool
        s *= "\n" * join(fill(' ', 21)) * "/$field"
    elseif v isa Union{Number, Symbol}
        s *= "\n" * join(fill(' ', 21)) * "/$field=" * string(v)
    else
        v = replace(v, "\n" => "\n" * join(fill(' ', 21)))
        s *= "\n" * join(fill(' ', 21)) * "/$field=\"" * v * "\""
    end
    return s
end

function Base.show(io::IO, gene::AbstractGene)
    buf = IOBuffer()
    print(buf, "     " * rpad(string(gene.feature), 16, ' '))
    print(buf, gene.locus)
    for field in names(gene.parent.genedata)
        if field in [:feature, :locus]
            continue
        end
        v = gene.parent.genedata[gene.index, field]
        if !ismissing(v)
            if v isa AbstractVector
                for i in eachindex(v)
                    print(buf, appendstring(field, v[i]))
                end
            else
                print(buf, appendstring(field, v))
            end
        end
    end
    println(io, String(take!(buf)))
end

Base.show(io::IO, mime::MIME"text/plain", genes::Array{Gene, 1}) = Base.show(io, genes)
function Base.show(io::IO, genes::Array{Gene, 1})
    for gene in genes
        show(io, gene)
    end
end

function Base.show(io::IO, chr::Chromosome)
    s = "Chromosome '$(chr.name)' ($(length(chr.sequence)) bp) with $(length(chr.genes)) annotations\n"
    print(io, s)
end

function Base.show(io::IO, genes::Vector{AbstractGene})
    print(join(genes, "\n"))
end

function Base.show(io::IO, locus::Locus)
    s = ""
    locus.strand == '-'     && (s *= "complement(")
    !locus.complete_left    && (s *= ">")
    if length(locus.order) > 0
        if locus.join
            s *= "join(" * join([join((r.start, r.stop), "..") for r in locus.order], ",") * ")"
        else
            s *= "order(" * join([join((r.start, r.stop), "..") for r in locus.order], ",") * ")"
        end
    else
        s *= join((locus.position.start, locus.position.stop), "..")
    end
    !locus.complete_right   && (s *= "<")
    locus.strand == '-'     && (s *= ")")
    print(io, s)
end

function formatsequence(sequence, io = IOBuffer)
    p = length(string(length(sequence))) + 2
    if length(sequence) > 60
        intervals = [i:i+60 for i in range(1; step = 60, stop = length(sequence)-60)]
        for interval in intervals
            println(io, lpad(string(first(interval)), p, ' '), " ", sequence[interval[1:10]],
                " ", sequence[interval[11:20]], " ", sequence[interval[21:30]],
                " ", sequence[interval[31:40]], " ", sequence[interval[41:50]],
                " ", sequence[interval[51:60]])
        end
    else
        intervals = [1:1]
    end
    i = intervals[end].stop
    if i <= length(sequence)
        print(io, lpad(i, p, ' '), " ")
        j = 0
        while i+j <= length(sequence)
            print(io, sequence[i+j])
            (j+1) % 10 == 0 && print(io, " ")
            j += 1
        end
    end
    return io
end

"""
    printgbk([io], chr)

Print `chr` in GenBank format.
"""
function printgbk(chrs::AbstractVector{C}) where {C <: Chromosome}
    io = IOBuffer()
    printgbk(io, chrs)
end
function printgbk(io::IO, chrs::AbstractVector{C}) where {C <: Chromosome}
    for chr in chrs
        printgbk(io, chr)
    end
    return io
end
function printgbk(chr::C) where {C <: Chromosome}
    io = IOBuffer()
    printgbk(io, chr)
end
function printgbk(io::IO, chr::C) where {C <: Chromosome}
    println(io, chr.header)
    println(io, rpad("FEATURES", 21, ' '), "Location/Qualifiers")
    println(io, chr.genes)
    println(io, "ORIGIN")
    formatsequence(chr.sequence, io)
    println(io)
    println(io, "//")
    return io
end

Base.copy(chr::Chromosome) = deepcopy(chr)
Base.copy(gene::AbstractGene) = deepcopy(gene)

function Base.isless(l1::Locus, l2::Locus)
    if (l1.position.start < l2.position.start) || (l1.position.start == l2.position.start && l1.position.stop > l2.position.stop)
        return true
    end
    return false
end
Base.isless(g1::Gene, g2::Gene) = ((g1.locus == g2.locus) && (g1.feature == "gene" && g2.feature != "gene")) || (g1.locus < g2.locus)

function Base.:(==)(x::Locus, y::Locus)
    (x.position == y.position) && (x.strand == y.strand) && (x.complete_left == y.complete_left) && (x.complete_right == y.complete_right) && (x.order == y.order)
end

function Base.sort!(G::AbstractArray{Gene}; args...)
    sort!(G[1].parent.genedata; cols = (:locus, :feature), args...)
end


text:
using Crux
using POMDPModels
using Test
using CUDA
using Flux
using Random

## mdp_data
d1 = mdp_data(ContinuousSpace(3), ContinuousSpace(4), 100)
d2 = mdp_data(ContinuousSpace(3), ContinuousSpace(4), 100, [:weight, :t, :advantage, :return, :logprob])
# @test_throws ErrorException  mdp_data(ContinuousSpace(3), ContinuousSpace(4), 100, [:bad_key]) # Now this just throws a warning
d3 = mdp_data(ContinuousSpace(3), ContinuousSpace(4), 100, ArrayType = USE_CUDA ? CuArray : Array)

@test d1[:s] == zeros(Float32, 3, 100) && d2[:s] == zeros(Float32, 3, 100) && size(d3[:s]) == (3, 100)
@test d1[:sp] == zeros(Float32, 3, 100) && d2[:sp] == zeros(Float32, 3, 100) && size(d3[:sp]) == (3, 100)
@test d1[:a] == zeros(Bool, 4, 100) && size(d2[:a]) == (4, 100) && size(d3[:a]) == (4, 100)
@test d1[:r] == zeros(Float32, 1, 100) && size(d2[:r]) == (1, 100) && size(d3[:r]) == (1, 100)
@test d1[:done] == zeros(Bool, 1, 100) && size(d2[:done]) == (1, 100) && size(d3[:done]) == (1, 100)
@test d2[:return] == zeros(1, 100) && d2[:advantage] == zeros(1, 100) && d2[:logprob] == zeros(1, 100) && d2[:t] == zeros(1, 100) && d2[:weight] == ones(1, 100)
@test !haskey(d1, :return) && !haskey(d1, :advantage) && !haskey(d1, :logprob) && !haskey(d1, :t) && !haskey(d1, :weight)
USE_CUDA && @test d3[:s] isa CuArray

## circular_indices
circ_inds(start, Nsteps, l) = mod1.(start:start + Nsteps - 1, l)
@test circ_inds(4, 60, 100) == [4:63 ...]
@test circ_inds(1, 100, 100) == [1:100 ...]
@test circ_inds(1, 101, 100) == [1:100 ..., 1]
@test circ_inds(1, 120, 100) == [1:100 ..., 1:20 ...]
@test circ_inds(90, 20, 100) == [90:100 ..., 1:9 ...]

## Get last N indices test
# A not full buffer
b = ExperienceBuffer(ContinuousSpace(2), ContinuousSpace(1), 100)
d = Dict(:s => 2*ones(2,50), :a => ones(Bool, 1,50), :sp => ones(2,50), :r => ones(1,50), :done => zeros(1,50), :weight=>zeros(1,50))
push!(b, d)

@test get_last_N_indices(b, 10) == collect(41:50)
@test get_last_N_indices(b, 1) == [50]
@test get_last_N_indices(b, 50) == collect(1:50)
@test get_last_N_indices(b, 51) == collect(1:50)
@test get_last_N_indices(b, 1000) == collect(1:50)

# full buffer
push!(b, d)
push!(b, d)
@test get_last_N_indices(b, 10) == collect(41:50)
@test get_last_N_indices(b, 1) == [50]
@test get_last_N_indices(b, 50) == collect(1:50)
@test get_last_N_indices(b, 51) == [100, collect(1:50)...]
@test get_last_N_indices(b, 100) == [collect(51:100)..., collect(1:50)...]
@test get_last_N_indices(b, 1000) == [collect(51:100)..., collect(1:50)...]

## Priority Param Construction
p = PriorityParams(100, α=0.7)
@test length(p.priorities) == 100
@test p.α == 0.7f0
@test p.max_priority == 1.0f0

p2 = PriorityParams(1000, p)
@test length(p2.priorities)==1000
@test p2.α == p.α
@test p2.max_priority == p.max_priority

## Construction
b = ExperienceBuffer(ContinuousSpace(2), DiscreteSpace(4), 100,)
bpriority = ExperienceBuffer(ContinuousSpace(2), DiscreteSpace(4), 50, prioritized=true)
USE_CUDA && (b_gpu = b |> gpu)

bcopy = ExperienceBuffer(deepcopy(b.data), elements=0)
@test b.data == bcopy.data
@test b.priority_params == bcopy.priority_params
bcopy = ExperienceBuffer(deepcopy(b.data), prioritized=true)
@test isprioritized(bcopy)

@test b isa ExperienceBuffer{Array}
if USE_CUDA
    @test b_gpu isa ExperienceBuffer{CuArray}
    @test length(keys(b.data)) == length(keys(b_gpu.data))
end

@test :s in keys(b.data)
@test :a in keys(b.data)
@test :sp in keys(b.data)
@test :r in keys(b.data)
@test :done in keys(b.data)

@test !haskey(b, :weight)
@test haskey(bpriority, :weight)
@test bpriority.data[:weight] == ones(Float32, 1, 50)
@test size(b[:a]) == (4,0)

# Buffer_like function
bsmall = buffer_like(b, capacity=3, device=cpu)
@test Crux.device(bsmall) == cpu
@test keys(bsmall) == keys(b)
@test capacity(bsmall) == 3
@test length(bsmall) == 0

## Base functions
@test keys(b) == keys(b.data)

@test size(b[:s]) == (2,0)
@test length(b) == 0
@test capacity(b) == 100

USE_CUDA && @test size(b_gpu[:s]) == (2,0)
USE_CUDA && @test length(b_gpu) == 0
USE_CUDA && @test capacity(b_gpu) == 100

@test isprioritized(bpriority)
@test Crux.device(b) == cpu
@test Crux.device(bpriority) == cpu
USE_CUDA && @test Crux.device(b_gpu) == gpu

# Convert to and from GPU
USE_CUDA && @test Crux.device(gpu(b)) == gpu
USE_CUDA && @test Crux.device(cpu(b_gpu)) == cpu

## push!
#push dictionary with one element
d = Dict(:s => 2*ones(2,1), :a => ones(Bool, 4,1), :sp => ones(2,1), :r => ones(1,1), :done => zeros(1,1), :weight=>zeros(1,1))
push!(b, d)
@test length(b) == 1
@test b[:s] == 2*ones(2,1)
@test b[:a] == ones(Int, 4,1)
@test b[:sp] == ones(2,1)
@test b[:r] == ones(1,1)
@test b[:done] == zeros(1,1)

# push dictionary with more than one element
d = Dict(:s => 3*ones(2,3), :a => rand(4,3) .< 0.5, :sp => 5*ones(2,3), :r => 6*ones(1,3), :done => ones(1,3), :weight=>zeros(1,3))
push!(b, d)
@test length(b) == 4
@test b[:s][:,2:end] ==  3*ones(2,3)
@test b[:a][:,2:end] == d[:a]
@test b[:sp][:,2:end] == 5*ones(2,3)
@test b[:r][:,2:end] == 6*ones(1,3)
@test b[:done][:,2:end] == ones(1,3)

# push a buffer
push!(b, b)
@test length(b) == 8
for k in keys(b)
    @test b[k][:, 1:4] == b[k][:, 5:8]
end

## Reservoir storage
b = ExperienceBuffer(ContinuousSpace(2), DiscreteSpace(4), 10,)
d3 = Dict(:s => 3*ones(2,3), :a => rand(4,3) .< 0.5, :sp => 5*ones(2,3), :r => 6*ones(1,3), :done => ones(1,3), :weight=>zeros(1,3), :episode_end=> ones(1,3))
d1 = Dict(:s => 2*ones(2,1), :a => ones(Bool, 4,1), :sp => ones(2,1), :r => ones(1,1), :done => zeros(1,1), :weight=>zeros(1,1), :episode_end => zeros(1,1))

for i=1:10
    push_reservoir!(b, d1)
    @test length(b) == i
end

push_reservoir!(b, d3)

## clear!
bcopy = deepcopy(b)
clear!(bcopy)
@test length(bcopy) == 0
@test bcopy.next_ind == 1

## minibatch
I = [1,2,4]
d2 = minibatch(b, I)
for k in keys(d2)
    @test all(d2[k] .== b[k][:, I])
end

## Split
@test Crux.split_batches(100, [0.5, 0.5]) == [50,50]
@test Crux.split_batches(100, [1.0]) == [100]
@test Crux.split_batches(100, [1/3, 1/3, 1/3]) == [34,33,33]
@test_throws AssertionError Crux.split_batches(100, 0.4)

v = rand(3,100)
bsplit = ExperienceBuffer(Dict(:a => v))

b1, b2 = split(bsplit, [0.5, 0.5])
@test length(b1) == 50
@test b1[:a] == v[:,1:50]
@test length(b2) == 50
@test b2[:a] == v[:,51:100]

## update_priorities!
update_priorities!(bpriority, [1,2,3], [1., 2., 3.])
@test bpriority.priority_params.max_priority == 3.0
@test bpriority.priority_params.priorities[1] ≈ 1f0^bpriority.priority_params.α
@test bpriority.priority_params.priorities[2] ≈ 2f0^bpriority.priority_params.α
@test bpriority.priority_params.priorities[3] ≈ 3f0^bpriority.priority_params.α

push!(bpriority, d)
push!(bpriority, d)
@test bpriority.priority_params.max_priority == 3.0
for i=1:6
    @test bpriority.priority_params.priorities[i] ≈ 3f0^bpriority.priority_params.α
end

## sampling

# uniform sample
t = ExperienceBuffer(ContinuousSpace(2), DiscreteSpace(4), 10)
rand!(t, b)

t = ExperienceBuffer(ContinuousSpace(2), DiscreteSpace(4), 3)
Random.seed!(0)
ids = rand(1:length(b), 3)
Random.seed!(0)
rand!(t, b)

for k in keys(t)
    @test t[k] == b[k][:,ids]
end

# Test the multi-buffer sampling
t1 = ExperienceBuffer(ContinuousSpace(2), DiscreteSpace(4), 10)
d = Dict(:s => ones(2,1), :a => ones(Bool, 4,1), :sp => ones(2,1), :r => ones(1,1), :done => zeros(1,1))
push!(t1, d)

t2 = ExperienceBuffer(ContinuousSpace(2), DiscreteSpace(4), 10)
d = Dict(:s => 2*ones(2,1), :a => ones(Bool, 4,1), :sp => ones(2,1), :r => ones(1,1), :done => zeros(1,1))
push!(t2, d)

t3 = ExperienceBuffer(ContinuousSpace(2), DiscreteSpace(4), 10)
d = Dict(:s => 3*ones(2,1), :a => ones(4,1), :sp => ones(2,1), :r => ones(1,1), :done => zeros(1,1))
push!(t3, d)

t = ExperienceBuffer(ContinuousSpace(2), DiscreteSpace(4), 10)
rand!(t, t1, t2, t3)

@test all(t[:s][:, 1:4] .== 1.0)
@test all(t[:s][:, 5:7] .== 2.0)
@test all(t[:s][:, 8:10] .== 3.0)

# Priority samples
bpriority[:s] .= rand(Float32, 2, 6)
update_priorities!(bpriority, [1:6...], [1.:6. ...])
t = ExperienceBuffer(ContinuousSpace(2), DiscreteSpace(4), 1000, [:weight])
priorities = [bpriority.priority_params.priorities[i] for i=1:length(bpriority)]

rand!(t, bpriority)

# Test the the frequency of samples is proportional to their priority
freqs = [sum(t.indices .== i) for i=1:6] ./ length(t)
probs = priorities ./ sum(priorities)
relerr = abs.(freqs .- probs) ./ probs
@test all(relerr .< 0.01)

@test all(t[:s] .== bpriority[:s][:,t.indices])
@test all(t[:weight] .<= 1.)

## merge
# bmerge = merge(b, b, capacity = 300)
# @test length(bmerge) == 200
# @test capacity(bmerge) == 300
# @test bmerge[:s][:, 1:100] == bmerge[:s][:, 101:200]

## Multi-D states
b2d = ExperienceBuffer(ContinuousSpace((2,2), UInt8), DiscreteSpace(4), 100;)
@test ndims(b2d[:s]) == 3
b2d[:s]

d = Dict(:s => 3*ones(2,2,1), :a => ones(4,1), :sp => ones(2,2,1), :r => ones(1,1), :done => zeros(1,1))
push!(b2d, d)

@test all(b2d[:s] .== 3)


text:
include("Types.jl")
using Types


text:
include("entity.jl")
include("decal.jl")
include("trigger.jl")
include("tiles.jl")
include("object_tiles.jl")
include("nodes.jl")
include("entity_map.jl")
include("enums.jl")

@valueequals Base.@kwdef mutable struct Room
    fgDecals::Array{Decal, 1} = Decal[]
    bgDecals::Array{Decal, 1} = Decal[]

    fgTiles::Tiles = Tiles()
    bgTiles::Tiles = Tiles()
    objTiles::ObjectTiles = ObjectTiles()

    entities::Array{Entity, 1} = Entity[]

    triggers::Array{Trigger, 1} = Trigger[]

    name::String = "lvl_1"
    color::Integer = 0

    position::Tuple{Integer, Integer} = (0, 0)
    size::Tuple{Integer, Integer} = (320, 184)

    musicLayer1::Bool = true
    musicLayer2::Bool = true
    musicLayer3::Bool = true
    musicLayer4::Bool = true

    musicProgress::String = ""
    ambienceProgress::String = ""

    delayAltMusicFade::Bool = false

    dark::Bool = false
    space::Bool = false
    underwater::Bool = false
    whisper::Bool = false

    music::String = "music_oldsite_awake"
    altMusic::String = ""

    disableDownTransition::Bool = false

    windPattern::String = "None"

    cameraOffsetX::Number = 0.0
    cameraOffsetY::Number = 0.0
end

function tileDataCopy!(to::Matrix, from::Matrix, width::Int, height::Int)
    indices = CartesianIndices((1:height, 1:width))

    copyto!(to, indices, from, indices)
end

function getCopySize(data::Matrix, width::Int, height::Int)
    h, w = size(data)

    return min(h, height), min(w, width)
end

# Set the fg and bg tiles to the size of the room
function updateTileSize!(room::Room, center::Char='0', border::Char='0')
    tw, th = ceil.(Int, room.size ./ 8)

    fgTiles = fill(center, (th, tw))
    bgTiles = fill(center, (th, tw))
    objTiles = fill(-1, (th, tw))

    fth, ftw = getCopySize(room.fgTiles.data, tw, th)
    bth, btw = getCopySize(room.bgTiles.data, tw, th)
    oth, otw = getCopySize(room.objTiles.data, tw, th)

    tileDataCopy!(fgTiles, room.fgTiles.data, ftw, fth)
    tileDataCopy!(bgTiles, room.bgTiles.data, btw, bth)
    tileDataCopy!(objTiles, room.objTiles.data, otw, oth)

    room.fgTiles = Tiles(fgTiles)
    room.bgTiles = Tiles(bgTiles)
    room.objTiles = ObjectTiles(objTiles)
end

const blacklistedRoomAttrs = String["position", "size", "color", "fgDecals", "bgDecals", "fgTiles", "bgTiles", "objTiles", "entities", "triggers"]

function Base.Dict(r::Room)
    res = Dict{String, Any}()

    for f in fieldnames(typeof(r))
        fs = String(f)
        value = getfield(r, f)

        if !(fs in blacklistedRoomAttrs)
            res[fs] = value
        end
    end

    res["__name"] = "level"

    res["x"] = r.position[1]
    res["y"] = r.position[2]

    res["c"] = r.color

    res["width"] = r.size[1]
    res["height"] = r.size[2]

    res["__children"] = Dict{String, Any}[
        Dict{String, Any}(
            "__name" => "solids",
            "innerText" => minimizeTilesString(r.fgTiles.data, '0', "")
        ),

        Dict{String, Any}(
            "__name" => "bg",
            "innerText" => minimizeTilesString(r.bgTiles.data, '0', ""),
        ),

        Dict{String, Any}(
            "__name" => "objtiles",
            "innerText" => minimizeTilesString(r.objTiles.data, -1, ','),
        ),

        Dict{String, Any}(
            "__name" => "fgtiles",
            "tileset" => "Scenery",
        ),

        Dict{String, Any}(
            "__name" => "bgtiles",
            "tileset" => "Scenery",
        ),

        Dict{String, Any}(
            "__name" => "entities",
            "__children" => Dict.(r.entities)
        ),

        Dict{String, Any}(
            "__name" => "triggers",
            "__children" => Dict.(r.triggers)
        ),

        Dict{String, Any}(
            "__name" => "fgdecals",
            "tileset" => "Scenery",
            "__children" => Dict.(r.fgDecals)
        ),

        Dict{String, Any}(
            "__name" => "bgdecals",
            "tileset" => "Scenery",
            "__children" => Dict.(r.bgDecals)
        )
    ]

    return res
end

max_stars_repo_name: CelestialCartographers/Maple

max_issues_repo_name: CelestialCartographers/Maple

max_forks_repo_name: CelestialCartographers/Maple


text:
using CSV, DataFrames, Dates

root_dir = dirname(pwd())
unzipped_files_path = *(root_dir, "/Datasets/Dataset_Bovespa/Unzipped/")
println(unzipped_files_path)
file_prefix = "COTAHIST_A"
file_begin_counter = 1995
file_end_counter = 2019

file_counter = file_begin_counter

time_inicio = Dates.now()

while (file_counter <= file_end_counter)
    file_path_concat_txt = unzipped_files_path*file_prefix*string(file_counter)*".TXT"
    file_path_concat_csv = unzipped_files_path*file_prefix*string(file_counter)
    println(file_path_concat_txt)
    file = open(file_path_concat_txt,"r")

    local count = 0

    local dataFrame = DataFrame(
        DATPRG = [], CODBDI = [], CODNEG = [],
        TPMERC = [], NOMRES = [], ESPECI = [],
        PRAZOT = [], MODREF = [], PREABE = [],
        PREMAX = [], PREMIN = [], PREMED = [],
        PREULT = [], PREOFC = [], PREOFV = [],
        TOTNEG = [], QUATOT = [], VOLTOT = [],
        PREEXE = [], INDOPC = [], DATVEN = [],
        FATCOT = [], PTOEXE = [], CODISI = [],
        DISMES = []
    )

    csv = open((file_path_concat_csv*".csv"),"w")

    for line in eachline(file)
        count += 1
        if (SubString(line,1,2) == "00")
            println("header")
        elseif (SubString(line,1,2) == "01")

            data_pregao = SubString(line,3,10)
            codbdi = SubString(line,11,12)
            codneg = SubString(line,13,24)
            tpmerc = SubString(line,25,27)
            nomres = SubString(line,28,39)
            especi = SubString(line,40,49)
            prazot = SubString(line,50,52)
            modref = SubString(line,53,56)
            preabe = SubString(line,57,69)
            premax = SubString(line,70,82)
            premin = SubString(line,83,95)
            premed = SubString(line,96,108)
            preult = SubString(line,109,121)
            preofc = SubString(line,122,134)
            preofv = SubString(line,135,147)
            totneg = SubString(line,148,152)
            quatot = SubString(line,153,170)
            voltot = SubString(line,171,188)
            preexe = SubString(line,189,201)
            indopc = SubString(line,202,202)
            datven = SubString(line,203,210)
            fatcot = SubString(line,211,217)
            ptoexe = SubString(line,218,230)
            codisi = SubString(line,231,242)
            dismes = SubString(line,243,245)

            push!(dataFrame,
                (
                    data_pregao, codbdi, codneg, tpmerc, nomres,
                    especi, prazot, modref, preabe, premax,
                    premin, premed, preult, preofc, preofv,
                    totneg, quatot, voltot, preexe, indopc,
                    datven, fatcot, ptoexe, codisi, dismes
                )
            )
        else
            println("fim")
        end
    end

    CSV.write(csv, dataFrame)

    close(csv) # TODO: verificar se funciona, fiz sem

    global file_counter += 1
end

time_fim = Dates.now()
println("TI: "*Dates.format(time_inicio, "S:s"))
println("TF: "*Dates.format(time_fim, "S:s"))
println(time_fim-time_inicio)


text:
function bea_query(url, querydict)
    r = try
        HTTP.request("GET", url, query = querydict)
    catch e
        if isa(e, Base.IOError)
            println("Caught $e")
            HTTP.request("GET", url, query = querydict)
        else
            rethrow(e)
        end
    end
    r_body = String(r.body)
    r_json = JSON.parse(r_body)
    if haskey(r_json["BEAAPI"]["Results"], "Error")
        error(r_json["BEAAPI"]["Results"]["Error"]["APIErrorDescription"])
    end
    return r_json
end

"""
    get_bea_datasets(b::Bea)

Return a `Dict` of dataset names and descriptions.

Arguments
---------
* `b` -- a [`Bea`](@ref) connection

"""
function get_bea_datasets(b::Bea)
    url = b.url
    key = b.key
    bea_method = "GetDataSetList"

    querydict = Dict("UserID" => key,
                     "Method" => bea_method,
                     "ResultFormat" => "JSON")

    response_json = bea_query(url, querydict)
    response_dict = response_json["BEAAPI"]["Results"]["Dataset"]
    datasets = Dict(d["DatasetName"] => d["DatasetDescription"] for d in response_dict)

    return datasets
end


text:
using BinaryBuilder

name = "NNPACK"
version = v"2018.06.22"

# Collection of sources required to build Ogg
sources = [
    "https://github.com/Maratyszcza/NNPACK.git" =>
    "c039579abe21f5756e0f0e45e8e767adccc11852",

    # Buckets of deps
    "https://github.com/Maratyszcza/cpuinfo.git" =>
    "4e8f04355892c5deb64a51731a6afdb544a4294d",
    "https://github.com/Maratyszcza/FP16.git" =>
    "4b37bd31c9cc1380ef9f205f7dd031efe0e847ab",
    "https://github.com/Maratyszcza/FXdiv.git" =>
    "811b482bcd9e8d98ad80c6c78d5302bb830184b0",
    "https://pypi.python.org/packages/e8/59/8c2e293c9c8d60f206fd5d0f6c8236a2e0a97832379ac319077441552c6a/opcodes-0.3.13.tar.gz" =>
    "1859c23143fe20daa4110be87a947cbf3eefa048da71dde642290213f251590c",
    "https://github.com/Maratyszcza/psimd.git" =>
    "3d8bfe7318423462a6d9e0c6537e75efd4822c49",
    "https://github.com/Maratyszcza/PeachPy.git" =>
    "01d15157a973a4ae16b8046313ddab371ea582db",
    "https://pypi.python.org/packages/16/d8/bc6316cf98419719bd59c91742194c111b6f2e85abac88e496adefaf7afe/six-1.11.0.tar.gz" =>
    "70e8a77beed4562e7f14fe23a786b54f6296e34344c23bc42f07b15018ff98e9",
    "https://github.com/Maratyszcza/pthreadpool.git" =>
    "3fb19c58b46f3cbc78a27c7b207a6eb7946633c0",

    "https://github.com/google/googletest/archive/release-1.8.0.tar.gz" =>
    "58a6f4277ca2bc8565222b3bbd58a177609e9c488e8a72649359ba51450db7d8",
]

# Bash recipe for building across all platforms
script = raw"""
# Prepare peachpy
cd ${WORKSPACE}/srcdir/PeachPy
python3 setup.py develop

mkdir ${WORKSPACE}/srcdir/NNPACK/build
cd ${WORKSPACE}/srcdir/NNPACK/build

# If we're on Linux, we need to add in `-lrt` as a library to link against
EXTRA_FLAGS=()
if [[ "${target}" == *-linux-* ]]; then
    EXTRA_FLAGS+=("-DCMAKE_CXX_STANDARD_LIBRARIES=-lrt")
fi

# On ARM/AArch64, use `clang` instead of `gcc`.
TOOLCHAIN="/opt/${target}/${target}.cmake"
if [[ "${target}" == arm-* ]] || [[ "${target}" == aarch64-* ]]; then
    TOOLCHAIN="/opt/${target}/${target}_clang.cmake"
fi

# NNPACK wants "armv7l", not just "arm", so GIVE IT WHAT IT WANTS
# Disabled for now because it thinks we don't have NEON
#if [[ "${target}" == arm-* ]]; then
#    sed -i.bak -e "s/set(CMAKE_SYSTEM_PROCESSOR arm)/set(CMAKE_SYSTEM_PROCESSOR armv7l)/g" "${TOOLCHAIN}"
#fi

cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE="${TOOLCHAIN}" \
    -DPYTHON_EXECUTABLE=/usr/bin/python3 \
    -DPYTHON_LIBRARY=/usr/lib/libpython3.so \
    -DPYTHON_SIX_SOURCE_DIR=$(echo ${WORKSPACE}/srcdir/six*/) \
    -DPYTHON_OPCODES_SOURCE_DIR=$(echo ${WORKSPACE}/srcdir/opcodes*/) \
    -DPYTHON_PEACHPY_SOURCE_DIR=${WORKSPACE}/srcdir/PeachPy \
    -DCPUINFO_SOURCE_DIR=${WORKSPACE}/srcdir/cpuinfo \
    -DFP16_SOURCE_DIR=${WORKSPACE}/srcdir/FP16 \
    -DFXDIV_SOURCE_DIR=${WORKSPACE}/srcdir/FXdiv \
    -DPSIMD_SOURCE_DIR=${WORKSPACE}/srcdir/psimd \
    -DPTHREADPOOL_SOURCE_DIR=${WORKSPACE}/srcdir/pthreadpool \
    -DGOOGLETEST_SOURCE_DIR=$(echo ${WORKSPACE}/srcdir/googletest*/) \
    -DCMAKE_POSITION_INDEPENDENT_CODE=ON \
    -DNNPACK_LIBRARY_TYPE="shared" \
    "${EXTRA_FLAGS[@]}" \
    ..
make -j${nproc} VERBOSE=1
make install

install_license ${WORKSPACE}/srcdir/NNPACK/LICENSE
"""

# Build only Linux and MacOS
platforms = filter(p -> (p isa Linux || p isa MacOS), supported_platforms())

# Build only for AArch64, x86_64 and i686 (armv7l disabled for now until NEON support is figured out)
platforms = filter(p -> arch(p) in (:aarch64, :x86_64, :i686), platforms)

# AArch64 musl seems to have problems linking against libgcc_s, admit defeat for now
platforms = filter(p -> !(arch(p) == :aarch64 && libc(p) == :musl), platforms)

# The products that we will ensure are always built
products = Product[
    LibraryProduct("libnnpack", :libnnpack),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; allow_unsafe_flags=true)


text:
function end_search_tasks!(results::Array{RemoteChannel})
    for channel in results
        close(channel)
    end
end


text:
### A Pluto.jl notebook ###
# v0.17.3

using Markdown
using InteractiveUtils

# ╔═╡ 64b72493-6cac-4b92-9a81-6a5df12f7875
begin
    import Pkg
    Pkg.develop("PoreMatMod")
end

# ╔═╡ 5401e009-923e-4a7f-9f3a-fd534f06d8b0
# load required packages (Pluto.jl will automatically install them)
using PoreMatMod, PlutoUI

# ╔═╡ 0f99b225-db96-485c-9e2c-1b4179601e53
using PoreMatMod.ExampleHelpers

# ╔═╡ 8d523993-6e85-443a-9949-12030552b457
md"""
## example: performing a substructure search to find the linkers in a MOF
"""

# ╔═╡ 74b45651-21d5-4332-a4a5-866ea1bb02b8
input_file_message()

# ╔═╡ 3ff0cd63-8bdf-4ba5-92b6-e2a52f7573c2
xtal_folder_message()

# ╔═╡ ccc7fc86-5f48-4cc9-bd5c-349fc1d55b0c
rc[:paths][:crystals]

# ╔═╡ b66dd0d5-5bac-4b23-a33f-17305b0d4728
moiety_folder_message()

# ╔═╡ a76e79ee-c82c-4771-818d-5380a5fb4c18
rc[:paths][:moieties]

# ╔═╡ be86f07e-0669-46c2-8c79-80e65dfcc6f2
md"""
!!! example \"the task\"
    we have the IRMOF-1 crystal structure, and wish to explore appending the acetamido functional group its BDC linkers using different replacement styles.

**Parent crystal structure**: first, we read in the `.cif` file describing the parent IRMOF-1 crystal structure.
"""

# ╔═╡ b53e7c38-d8f5-4f28-a9dc-f7902a86fdb2
begin
    # read in the parent xtal
    parent = Crystal("IRMOF-1.cif")    # load .cif file
    infer_bonds!(parent, true)      # infer bonds
    view_structure(parent)          # view structure
end

# ╔═╡ b402e79a-784b-4d8b-82f1-df4fe1cedad1
fragment_construction_note()

# ╔═╡ 9b2d534a-4f78-4950-add1-9ba645669bb9
md"""
**Query fragment**: next, we construct a query fragment as a _p_-phenylene fragment to match that on the BCD linker of the IRMOF-1 parent structure.
"""

# ╔═╡ 4be03110-61ab-4cd6-b3fe-7d51ac5ee771
query = moiety("p-phenylene.xyz");

# ╔═╡ e5eaaef4-13a0-48bd-9f2b-5040b2d10ac1
view_query_or_replacement("p-phenylene.xyz")

# ╔═╡ 559ef3c3-a176-4d65-8958-810c9b0b32c5
with_terminal() do
    display_query_or_replacement_file("p-phenylene.xyz")
end

# ╔═╡ 127a2bef-3903-4801-bc75-00a6dde2bc6e
md"""
### Substructure Search

Find substructures in the `parent` that match the `query`:
"""

# ╔═╡ b7b46022-aee0-4d51-8a5e-0c8c005f341a
search = query ∈ parent

# ╔═╡ 536d3474-bc9c-4576-a120-826020f8d772
hits = isomorphic_substructures(search)

# ╔═╡ ec84ba70-b83e-4a46-9037-aed54bb07b07
view_structure(hits)

# ╔═╡ Cell order:
# ╟─64b72493-6cac-4b92-9a81-6a5df12f7875
# ╟─8d523993-6e85-443a-9949-12030552b457
# ╠═5401e009-923e-4a7f-9f3a-fd534f06d8b0
# ╠═0f99b225-db96-485c-9e2c-1b4179601e53
# ╟─74b45651-21d5-4332-a4a5-866ea1bb02b8
# ╟─3ff0cd63-8bdf-4ba5-92b6-e2a52f7573c2
# ╠═ccc7fc86-5f48-4cc9-bd5c-349fc1d55b0c
# ╟─b66dd0d5-5bac-4b23-a33f-17305b0d4728
# ╠═a76e79ee-c82c-4771-818d-5380a5fb4c18
# ╟─be86f07e-0669-46c2-8c79-80e65dfcc6f2
# ╠═b53e7c38-d8f5-4f28-a9dc-f7902a86fdb2
# ╟─b402e79a-784b-4d8b-82f1-df4fe1cedad1
# ╟─9b2d534a-4f78-4950-add1-9ba645669bb9
# ╠═4be03110-61ab-4cd6-b3fe-7d51ac5ee771
# ╟─e5eaaef4-13a0-48bd-9f2b-5040b2d10ac1
# ╟─559ef3c3-a176-4d65-8958-810c9b0b32c5
# ╟─127a2bef-3903-4801-bc75-00a6dde2bc6e
# ╠═b7b46022-aee0-4d51-8a5e-0c8c005f341a
# ╠═536d3474-bc9c-4576-a120-826020f8d772
# ╠═ec84ba70-b83e-4a46-9037-aed54bb07b07


text:
using SemanticWords
using JSON
using KernelMethods.KMap
using KernelMethods.Scores
using KernelMethods.Kernels
using KernelMethods.Supervised
using SimilaritySearch
using NearNeighborGraph
using TextModel

function parse_vector(vstring, fdelim=',')
    arr = split(vstring, fdelim)
    m = length(arr)
    vec = Vector{Float32}(m)
    @inbounds for i in 1:m
        num = parse(Float32, arr[i])
        isnan(num) && error("parsed NaN object $(arr[i]) -- $vstring")
        isinf(num) && error("parsed Inf object $(arr[i]) -- $vstring")
        vec[i] = num
    end

    DenseCosine(vec)
end

function parse_vector_list(text, vdelim=' ')
    h = DenseCosine{Float32}[]
    for vstring in split(strip(text), vdelim)
        try
            v = parse_vector(vstring)
            push!(h, v)
        catch
            info("*** ignoring:", vstring)
            continue
        end
    end

    h
end

function read_vectors(filename, key_vectors)
    f = open(filename)
    db = DenseCosine{Float32}[]
    words = String[]

    for (i, line) in enumerate(eachline(f))
        v = JSON.parse(line)
        vecs = parse_vector_list(v[key_vectors])
        append!(db, vecs)
        if i % 50000 == 1
            info("read $i vectors from $filename")
        end
    end
    close(f)

    [string(i) for i in 1:length(db)], db
end

function create_config(nlist, qlist, skiplist=[])
    config = TextConfig()
    config.del_usr = false
    config.del_punc = false
    config.del_num = false
    config.del_url = true
    config.del_diac = false
    config.nlist = collect(Int, nlist)
    config.qlist = collect(Int, qlist)
    config.skiplist = collect(Int, skiplist)
    config
end

function create_encoder(centroids, k, recall)

    function encode(vlist::String)::String
        X = Char[]
        for v in parse_vector_list(vlist)
            res = search(centroids, v, KnnResult(k))
            for p in res
                push!(X, Char(p.objID+33))
            end
            push!(X, ' ')
        end

        X |> join
    end

    encode
end

function encode_corpus(encode::Function, corpus, key_text)
    for tweet in corpus
        # tweet["original-text"] = tweet[key_text]
        tweet[key_text] = encode(tweet[key_text])  #join(encode(tweet[key_text]), ' ')
    end

    corpus
end

function vectorize_collection(model, col, key_klass, key_text, allow_zero_length)
	X = Vector{VBOW}(length(col))
	y = Vector{String}(length(col))

	i = 0
	for tweet in col
		vec = vectorize(tweet[key_text], model)
		if allow_zero_length || length(vec) > 0
			i += 1
			X[i] = vec
			y[i] = get(tweet, key_klass, "")
		end
	end

	if !allow_zero_length
		resize!(X, i)
		resize!(y, i)
	end
	X, y
end

function rocchio_model(rocchioname, model, gettrain, test; key_klass="klass", key_text="text")
    if isfile(rocchioname)
        @load rocchioname X̂ ŷ
    else
        info("vectorizing train")
        X, y = vectorize_collection(model, gettrain(), key_klass, key_text, false)

        info("computing centroids")
        ŷ = unique(y)
        X̂ = [centroid!(X[y .== label]) for label in ŷ]
        @save rocchioname X̂ ŷ
    end

    index = Sequential(X̂, angle_distance)
    info("vectorizing test")
	Xtest, ytest = vectorize_collection(model, test, key_klass, key_text, true)

    info("predicting test")
    ypred = String[]
    Y = []
    for (i, x) in enumerate(Xtest)
        res = search(index, x, KnnResult(length(ŷ)))
        label = ŷ[first(res).objID]
        push!(ypred, label)
        m = test[i]
        delete!(m, "text")
        m[key_klass] = label
        A = [(p.objID, p.dist) for p in res]
        sort!(A)
        m["decision_function"] = [a[2] for a in A]
        push!(Y, m)
    end

    _scores = scores(ytest, ypred)
    Y, _scores
end

function create_codebook()
    numcenters = parse(Int, get(ENV, "numcenters", "10000"))
    k = parse(Int, get(ENV, "k", "7"))
    # kcut = parse(Int, get(ENV, "kcut", "")
    kind = get(ENV, "kind", "random")
    maxiter = parse(Int, get(ENV, "maxiter", "10"))
    tol = parse(Float64, get(ENV, "tol", "0.001"))
    @assert (kind in ("random", "fft", "dnet")) "kind=random|fft|dnet; specifies the algorithm to compute the codebook"
    vecfile = get(ENV, "vectors", "")
    @assert (length(vecfile) > 0) "vectors=file; the embedding file is mandatory"
    modelname = get(ENV, "model", "")
    @assert (length(modelname) > 0) "model=outfile"
    recall = parse(Float64, get(ENV, "recall", "0.9"))

    info("creating codebook: numcenters:$numcenters, k:$k, kind:$kind, recall=$recall, maxiter=$maxiter, tol:$tol, vecfile=$vecfile")
    cc = ApproxKDCentroids(numcenters, k, maxiter=maxiter, tol=tol, recall=recall)

    if get(ENV, "vformat", "dense") == "dense"
        begin
            words, X = read_vectors(vecfile, "text")
            info(length(words), "--", length(X))
            centroids, codes, rlist = codebook(cc, X)
            centroids = SemanticWords.create_index(centroids, recall)
            @save modelname words centroids codes
        end
    else
        error("run-images.jl works only with dense vectors")
        # begin
        #     words, X = read_sparse(vecfile)
        #     centroids, codes, rlist = codebook(cc, X)
        #     @save modelname words centroids codes
        # end
    end
end

function main()
    action = get(ENV, "action", "")
    key_klass = get(ENV, "klass", "klass")
    key_text = get(ENV, "text", "text")
    weighting = get(ENV, "weighting", "tfidf")
    @assert (weighting in ["tfidf", "tf", "idf", "freq", "entropy"]) "unknown weighting scheme $weighting"
    smoothing = parse(Int, get(ENV, "smoothing", "3"))
    nlist = [parse(Int, i) for i in split(get(ENV, "nlist", "1"), ',') if length(i) > 0]
    qlist = [parse(Int, i) for i in split(get(ENV, "qlist", "1,3,5"), ',') if length(i) > 0]
    filter_low = parse(Int, get(ENV, "filter_low", "1"))
    filter_high = parse(Float64, get(ENV, "filter_high", "1.0"))
    recall = parse(Float64, get(ENV, "recall", "0.9"))
    k = parse(Int, get(ENV, "k", "7"))

    if action == "model"
        create_codebook()
    elseif action == "translate"
        modelname = get(ENV, "model", "")
        @assert (length(modelname) > 0) "argument model=modelfile is mandatory"
        @load modelname words centroids codes
        encode = create_encoder(centroids, k, recall)

        key_text = get(ENV, "text", "text")
        info("*** translating $ARGS with model $modelname")
        for arg in ARGS
            outname = replace(modelname, ".jld2", "") * "." * basename(arg)
            if !isfile(outname)
                f = open(outname, "w")
                for tweet in encode_corpus(encode, [JSON.parse(line) for line in readlines(arg)], key_text)
                    println(f, JSON.json(tweet))
                end
                close(f)
            end
        end
    elseif action == "evaluate"
        trainfile, testfile = ARGS[1], ARGS[2]
        train = []
        function gettrain()
            if length(train) == 0
                for line in readlines(trainfile)
                    push!(train, JSON.parse(line))
                end
            end

            train
        end
        test = [JSON.parse(line) for line in readlines(testfile)]

		config = create_config(nlist, qlist)

        _file = trainfile * ".weighting=$(weighting).filter_low=$(filter_low).filter_high=$(filter_high).nlist=$(join(map(string, nlist), ',')).qlist=$(join(map(string, qlist), ','))"
        rocchiofile = _file * ".rocchio.jld2"
        vmodelfile = _file * ".vmodel.jld2"

        if weighting == "entropy"
            X = [item[key_text] for item in gettrain()]
            y = [item[key_klass] for item in gettrain()]
            le = LabelEncoder(y)
            model = DistModel(config, X, transform.(le, y))
            model = EntModel(model, smoothing)
        else
            if isfile(vmodelfile)
                @load vmodelfile vmodel
            else
                vmodel = VectorModel(config)
                vmodel.filter_low = filter_low
                vmodel.filter_high = filter_high
                X = [item[key_text] for item in gettrain()]
                TextModel.fit!(vmodel, X)
                @save vmodelfile vmodel
            end

            info("vocabulary size: ", length(vmodel.W))
            if weighting == "tfidf"
                model = TfidfModel(vmodel)
            elseif weighting == "tf"
                model = TfModel(vmodel)
            elseif weighting == "idf"
                model = IdfModel(vmodel)
            elseif weighting == "freq"
                model = FreqModel(vmodel)
            end
        end

        Y, _scores = rocchio_model(rocchiofile, model, gettrain, test, key_klass=key_klass, key_text=key_text)
        println(JSON.json(_scores))
        f = open(testfile * ".predicted", "w")
        for y in Y
            println(f, JSON.json(y))
        end
        close(f)
    else
        exit(1)
    end
end

if !isinteractive()
	main()
end


text:
module Deldir

using DataFrames
using Deldir_jll

export
    deldir,
    voronoiarea,
    edges

include("wrapper.jl")
include("misc.jl")

end # module


text:
using OpenQuantumBase, Test, Random

coupling = ConstantCouplings(["Z"], unit=:ħ)
jfun(t₁, t₂) = 1.0
jfun(τ) = 1.0
# TODO: add test for unitary using StaticArrays
# const Sx = SMatrix{2,2}(σx)
unitary(t) = exp(-1.0im * t * σx)
tf = 5.0
u0 = PauliVec[1][1]
ρ = u0 * u0'
kernels = [(((1, 1),), coupling, OpenQuantumBase.SingleFunctionMatrix(jfun))]

L = OpenQuantumBase.quadgk((x) -> unitary(x)' * σz * unitary(x), 0, 5)[1]
ulind = OpenQuantumBase.ULELiouvillian(kernels, unitary, tf, 1e-8, 1e-6)
p = ODEParams(nothing, 5.0, (tf, t) -> t / tf)
dρ = zero(ρ)
ulind(dρ, ρ, p, 5.0)
@test dρ ≈ L * ρ * L' - 0.5 * (L' * L * ρ + ρ * L' * L) atol = 1e-6 rtol = 1e-6

# test for EᵨEnsemble
u0 = EᵨEnsemble([0.5, 0.5], [PauliVec[3][1], PauliVec[3][2]])
@test all((x)->x∈[PauliVec[3][1], PauliVec[3][2]], [sample_state_vector(u0) for i in 1:4])

Lz = Lindblad((s) -> 0.5, (s) -> σz)
Lx = Lindblad((s) -> 0.2, (s) -> σx)
Ł = OpenQuantumBase.lindblad_from_interactions(InteractionSet(Lz))[1]
p = ODEParams(nothing, 5.0, (tf, t) -> t / tf)
dρ = zero(ρ)
Ł(dρ, ρ, p, 5.0)
@test dρ ≈ [0 -0.5; -0.5 0]

cache = zeros(2, 2)
update_cache!(cache, Ł, p, 5.0)
@test cache == -0.25 * σz' * σz

Ł = OpenQuantumBase.lindblad_from_interactions(InteractionSet(Lz, Lx))[1]
p = ODEParams(nothing, 5.0, (tf, t) -> t / tf)
dρ = zero(ρ)
Ł(dρ, PauliVec[2][1] * PauliVec[2][1]', p, 5.0)
@test dρ ≈ [0 0.7im; -0.7im 0]

cache = zeros(2, 2)
update_cache!(cache, Ł, p, 5.0)
@test cache == -0.25 * σz' * σz - 0.1 * σx' * σx

Random.seed!(1234)
sample_res = [OpenQuantumBase.lind_jump(Ł, PauliVec[1][1], p, 0.5) for i in 1:4]
@test all((x)->x∈[σz, σx], sample_res)


text:
using DelimitedFiles
using Distributed
using DifferentialEquations
using ProgressMeter
using TaylorIntegration
using Base.Threads

include("PSS_Taylor.jl")
include("Yfind.jl")
include("is_it.jl")
using PyCall
pygui(:qt)
using PyPlot
pygui(true)

Energy=.3
t_end=100.
n_iter_P=101;
n_iter_Q=50;

ArrP=range(1e-8,stop=.99,length=n_iter_P)
ArrQ=range(1e-8,stop=1,length=n_iter_Q)
Q_0=zeros(0)
P_0=zeros(0)

Q_1=zeros(0)
P_1=zeros(0)

Q_2=zeros(0)
P_2=zeros(0)

Q_3=zeros(0)
P_3=zeros(0)

Threads.@threads for j=1:n_iter_P
    for k=1:n_iter_Q
    Q0,P0,Q1,P1,Q2,P2,Q3,P3=PSS_function_taylor(ArrQ[k],ArrP[j], Energy, t_end)
        if Q0!=0 && P0!=0
            push!(Q_0,Q0)
            push!(Q_1,Q1)
            push!(Q_2,Q2)
            push!(Q_3,Q3)

            push!(P_0,P0)
            push!(P_1,P1)
            push!(P_2,P2)
            push!(P_3,P3)
        end
    end
end
figure()
# axis([ -2.5, 2.5,-1.1, 1.1])
Q0=vcat(Q_0,-Q_0,Q_0,-Q_0)
P0=vcat(P_0,P_0,-P_0,-P_0)
T0=hcat(Q0,P0)
plot(P0,Q0 ,color="r",".",markersize=2, markeredgewidth=.1)

Q1=vcat(Q_1,-Q_1,Q_1,-Q_1)
P1=vcat(P_1,P_1,-P_1,-P_1)
T1=hcat(Q1,P1)
plot(P1,Q1 ,color="b",".",markersize=2, markeredgewidth=.1)

Q2=vcat(Q_2,-Q_2,Q_2,-Q_2)
P2=vcat(P_2,P_2,-P_2,-P_2)
T2=hcat(Q2,P2)
plot(P2,Q2 ,color="g",".",markersize=2, markeredgewidth=.1)

Q3=vcat(Q_3,-Q_3,Q_3,-Q_3)
P3=vcat(P_3,P_3,-P_3,-P_3)
T3=hcat(Q3,P3)
plot(P3,Q3 ,color="m",".",markersize=2, markeredgewidth=.1)

outfile_0 = "Taylor_outfile0_"*string(Energy)*"_nP_"*string(n_iter_P)*"_nQ_"*string(n_iter_Q)*".dat"
outfile_1 = "Taylor_outfile1_"*string(Energy)*"_nP_"*string(n_iter_P)*"_nQ_"*string(n_iter_Q)*".dat"
outfile_2 = "Taylor_outfile2_"*string(Energy)*"_nP_"*string(n_iter_P)*"_nQ_"*string(n_iter_Q)*".dat"
outfile_3 = "Taylor_outfile3_"*string(Energy)*"_nP_"*string(n_iter_P)*"_nQ_"*string(n_iter_Q)*".dat"

a = open(outfile_0, "w")
b = open(outfile_1, "w")
c = open(outfile_2, "w")
d = open(outfile_3, "w")

writedlm(outfile_0, T0)
writedlm(outfile_1, T1)
writedlm(outfile_2, T2)
writedlm(outfile_3, T3)

close(a)
close(b)
close(c)
close(d)


text:
@testset "trace" begin
    stpara = (12345.6, 10.123, 10.0, 5, 1.2)
    st = SolverState(1, stpara...)
    @test sprint(show, st) == """
         iter=1         f=1.23e+04  rss=10.1      xnorm=10        idx=5         step=1.2
        """
    @test sprint(show, st, true) == """
         thread  1:  iter=1         f=1.23e+04  rss=10.1      xnorm=10        idx=5         step=1.2
        """
    tr = SolverTrace(Float64)
    update!(tr, 1, stpara..., true, true, true)
    @test length(tr.states) == 1
    @test tr.states[1] == st
    update!(tr, 2, stpara..., true, false, true)
    @test length(tr.states) == 2
    @test tr[1] == st
    @test sprint(show, tr) == """
         iter=1         f=1.23e+04  rss=10.1      xnorm=10        idx=5         step=1.2
         iter=2         f=1.23e+04  rss=10.1      xnorm=10        idx=5         step=1.2
        """
    update!(tr, 3, stpara..., false, false, false)
    @test length(tr.states) == 2
    @test update!(nothing, 4, stpara..., false, false, false) === nothing
end

@testset "fit" begin
    A = [1.0 0.0; 0.0 1.0]
    b = [1.0, 1.0]
    x0 = [1.0, 0.0]
    ca = Cache(A)
    r = convexfit(A, b, x0=x0, show_trace=true, store_trace=true, cache=ca)
    @test r.x0 === x0
    @test r.sol == [0.5, 0.5]
    @test r.fit == [0.5, 0.5]
    @test r.iter == 2
    @test r.f == 0.5
    @test r.rss == 0.5
    @test r.λ == 0
    @test r.xnorm == 0.5
    @test r.dfnorm == 0
    @test r.dxnorm == 0
    @test r.f_converged && r.x_converged
    @test sprint(show, r) == "SolverResult(converged=true, iter=2, rss=0.5, λ=0.0, xnorm=0.5)"
    @test sprint(show, MIME("text/plain"), r) == """
        convexfit converged after 2 iterations:
         |f(x) - f(x')| = 0.0e+00 < 1.0e-06: true
         |x - x'| = 0.0e+00 < 1.0e-06: true
         f=0.5       rss=0.5       λ=0         xnorm=0.5     """
    @test sprint(show, r.trace) == """
         iter=1         f=1         rss=1         xnorm=1         idx=2         step=0.5
         iter=2         f=0.5       rss=0.5       xnorm=0.5       idx=1         step=0
        """

    r = convexfit(A, b, x0=[1,0], show_trace=true)
    @test r.x0 == [1.0, 0.0]
    r = convexfit(A, b, x0=view(A,:,1))
    @test r.x0 == [1.0, 0.0]

    r = convexfit(A, b, ftol=0)
    @test r.x_converged
    @test !r.f_converged

    # The case when convergence is reached by dx all being zero
    A = [1 0; 1 0]
    r = convexfit(A, b, x0=[1.0,0.0], store_trace=true)
    @test r.x_converged && r.f_converged
    @test iszero(r.dfnorm)
    @test iszero(r.dxnorm)
    @test r.iter == 0
    @test isempty(r.trace.states)

    A = [1 0 1; 0 1 1]
    b = [1, 1]
    ca = Cache(convert(Matrix{Float64}, A))
    r = convexfit(A, b, show_trace=true, cache=ca)
    @test all(r.x0.==1/3)
    @test r.sol[3] ≈ 1
    @test r.fit == b
    @test r.iter == 1
    @test r.f == r.rss ≈ 0
    @test r.dfnorm ≈ 0
    @test r.dxnorm ≈ 0
    @test r.f_converged && r.x_converged

    r = convexfit(A, b, 0.5, xtol=0, ftol=eps(Float64), show_trace=true)
    @test r.sol ≈ [0.2, 0.2, 0.6]
    @test r.f ≈ 0.3

    ca = Cache(convert(Matrix{Float64}, A))
    fill!(ca.Ax, 0.0)
    ca1 = deepcopy(ca)
    r = convexfit(A, A, cache=ca, 1)
    @test r[1].sol[1] ≈ r[2].sol[2]
    @test r[1].fit[1] ≈ r[2].fit[2]
    @test r[3].sol ≈ [0.25, 0.25, 0.5]
    @test r[3].f ≈ 0.5
    # Check that the passed Cache is used when multithreads=false
    @test ca.Ax != ca1.Ax

    if testthreads
        ca = Cache(convert(Matrix{Float64}, A))
        fill!(ca.Ax, 0.0)
        ca1 = deepcopy(ca)
        r1 = convexfit(A, A, 1, cache=ca, multithreads=true)
        @test r1 == r
        # Check that the passed Cache is not used when multithreads=true
        @test ca.Ax == ca1.Ax
    end

    r = convexfit(A, 0.5, show_trace=true)
    @test length(r) == 3
    @test r[1].sol ≈ [0.25, 0.75]
    @test r[2].sol ≈ [0.25, 0.75]
    @test r[3].sol ≈ [0.5, 0.5]

    if testthreads
        r1 = convexfit(A, 0.5, show_trace=true, multithreads=true)
        @test r1 == r
    end

    r = convexfit(A, x0=[1,0])
    @test r[1].x0 == [1.0, 0.0]
    # Check that x0 is not copied after being converted
    @test r[1].x0 === r[2].x0

    r = convexfit(A, loo=false, 0.5)
    @test length(r) == 3
    @test r[3].sol ≈ [0.2, 0.2, 0.6]

    if testthreads
        r1 = convexfit(A, loo=false, 0.5, multithreads=true)
        @test r1 == r
    end

    @test_logs (:warn, "sum of elements in x0 (3.0) is normalized to 1") convexfit(A, b, x0=[2.0,0.0,1.0])
    @test_throws ArgumentError convexfit(A, b, -1.0)
    @test_throws ArgumentError convexfit(reshape([1.0,0.0],2,1))
end


text:
using StructEquality

const SUFFIX = "_Parsed"
"""
```julia
EP.@exprparser struct MySymbol
  symbol = EP.anything = :default_parsed_value
end
```
is transformed to
```julia
Base.@kwdef struct MySymbol <: EP.ExprParserWithParsed
  symbol = anything
end
Base.@kwdef mutable struct MySymbol_Parsed{T} <: EP.ExprParsed
  symbol = :default_parsed_value
end
EP.ExprParsed(::Base.Type{MySymbol}) = MySymbol_Parsed
```

It defines the basics for an usual ExprParser, namely
- the Parser type itself like specified in the original struct. It is always immutable - if you feel the need of
  mutating a parser, try to construct a new parser instead.
- a corresponding Parsed type which will be used to hold parsed values. This is intentionally mutable as a usual
  workflow consists of adapting the parsed values to the needs of your macro, and if everything is changed, transform
  it back to an Expr using `to_expr(parsed)`.
- a mapping from the `ExprParser` to the `ExprParsed`

Additionally, the created MySymbol Parser supports the following default `parse_expr` functionality
```julia
parser = MySymbol()
parse_expr(parser, symbol = :hi)
```
which translates to
```julia
parser = MySymbol()
MySymbol_Parsed(symbol = parse_expr(parser.symbol, :hi))
```
This is generic, and works similar if you have multiple fields.

----------------

Finally, in order to finish your custom ExprParser definition, you just need to specialize the two main functions
- `parse_expr(mysymbolparser::MySymbol, expr)`
- `to_expr(mysymbolparsed::MySymbol_Parsed)`

```julia
function EP.parse_expr(mysymbolparser::MySymbol, expr)
  # do your custom parsing
  # use @passert for checking parse assertions (it will have a nice and detailed default error message)
  # construct your parsed result
  MySymbol_Parsed(symbol = ...)
end

function EP.to_expr(parsed::MySymbol_Parsed)
  # create a proper `Base.Expr` from your parsed result
  # in this case it is simple
  parsed.symbol
end
```
"""
macro exprparser(combined_struct_expr)
  # TODO performance improvement?: add typeparameters for every field (?)
  if combined_struct_expr.head == :block
    combined_struct_expr = combined_struct_expr.args[1]
  end
  @assert combined_struct_expr.head == :struct "expecting struct"
  @assert combined_struct_expr.args[1] == false "expecting immutable struct as parser"
  @assert(!isa(combined_struct_expr.args[2], Base.Expr) || combined_struct_expr.args[2].head != :(<:),
    "please ommit the inheritance notation, it will be provided by the macro")

  parser_struct_name::Base.Symbol = _get_struct_name(combined_struct_expr)
  parser_struct_expr = _parser_struct_from_combined_struct(combined_struct_expr, parser_struct_name)

  parsed_struct_name = Base.Symbol(parser_struct_name, SUFFIX)
  parsed_struct_expr = _parsed_struct_from_combined_struct(combined_struct_expr, parsed_struct_name)

  # add inheritance
  parser_struct_expr.args[2] = :($(parser_struct_expr.args[2]) <: ExprParsers.ExprParserWithParsed)
  parsed_struct_expr.args[2] = :($(parsed_struct_expr.args[2]) <: ExprParsers.ExprParsed)

  # TODO Parsed(..) calls also for typeparameters?
  esc(quote
    ExprParsers.StructEquality.@def_structequal Base.@kwdef $parser_struct_expr
    ExprParsers.StructEquality.@def_structequal Base.@kwdef $parsed_struct_expr
    ExprParsers.ExprParsed(::Base.Type{<:$parser_struct_name}) = $parsed_struct_name
  end)
end

"""
  returns name of struct definition as symbol
"""
_get_struct_name(expr::Base.Expr) = _get_struct_name(Val{expr.head}(), expr.args)
_get_struct_name(::Val{:struct}, args) = _get_struct_name(args[2])
_get_struct_name(::Union{Val{:(<:)}, Val{:(::)}, Val{:curly}}, args) = _get_struct_name(args[1])
_get_struct_name(symbol::Base.Symbol) = symbol

_rename_struct_name(newname, expr::Base.Expr) = _rename_struct_name(newname, Val{expr.head}(), expr.args)
_rename_struct_name(newname, symbol::Base.Symbol) = newname
_rename_struct_name(newname, any) = any
_rename_struct_name(newname, head::Union{Val{:(<:)}, Val{:(::)}, Val{:curly}}, args) = Base.Expr(get(head), _rename_struct_name(newname, args[1]), args[2:end]...)
_rename_struct_name(newname, ::Val{:struct}, args) = _rename_struct_name(newname, args[2])

"""
- makes it mutable
- adds SUFFIX to struct name
- deletes all default values
"""
function _parser_struct_from_combined_struct(struct_expr::Base.Expr, newname)
  @assert struct_expr.head == :struct

  # take first default value and ommit a possible second one
  _default(expr::Base.Expr) = _default(expr, Val{expr.head}(), expr.args)
  _default(any) = any
  _default(expr, head::Val{:block}, args) = Base.Expr(:block, _default.(args)...)
  _default(expr, head::Val{:(=)}, args) = _default_assignment(expr, args[1], args[2])
  _default_assignment(expr, a, b) = expr # if there is no double assignment, but hence only a single, we take the single
  _default_assignment(expr, a, b::Base.Expr) = _default_assignment(expr, a, Val(b.head), b.args)
  _default_assignment(expr, a, head::Val{:(=)}, bargs) = :($a = $(bargs[1]))  # if there are two assignments, take the first
  _default_assignment(expr, a, head, bargs) = expr  # if there is only one assignment, take it

  Base.Expr(:struct, true, _rename_struct_name(newname, struct_expr.args[2]), _default(struct_expr.args[3]))
end

"""
- makes it mutable
- adds SUFFIX to struct name
- deletes all default values
"""
function _parsed_struct_from_combined_struct(struct_expr::Base.Expr, newname)
  @assert struct_expr.head == :struct

  # take first default value and ommit a possible second one
  _default(expr::Base.Expr) = _default(expr, Val{expr.head}(), expr.args)
  _default(any) = any
  _default(expr, head::Val{:block}, args) = Base.Expr(:block, _default.(args)...)
  _default(expr, head::Val{:(=)}, args) = _default_assignment(expr, args[1], args[2])
  _default_assignment(expr, a, b) = a # if there is no double assignment, but hence only a single, we take none
  _default_assignment(expr, a, b::Base.Expr) = _default_assignment(expr, a, Val(b.head), b.args)
  _default_assignment(expr, a, head::Val{:(=)}, bargs) = :($a = $(bargs[2]))  # if there are two assignments, take the second
  _default_assignment(expr, a, head, bargs) = a  # if there is only one assignment, we take none

  Base.Expr(:struct, true, _rename_struct_name(newname, struct_expr.args[2]), _default(struct_expr.args[3]))
end


text:
import LinearAlgebra

"""
Abstract decomposition of N-mode tensor.
"""
abstract type TensorDecomposition{T<:Number, N} end

"""
Returns relative error between the re-composed and the original tensors.
"""
function rel_residue(recomposed::AbstractArray{T, N}, tensor::AbstractArray{T, N}) where {T,N}
	nans = .!isnan.(tensor)
    LinearAlgebra.norm(recomposed[nans] .- tensor[nans]) / LinearAlgebra.norm(tensor[nans])
end

"""
Returns relative error between re-composed and the original tensor.
"""
rel_residue(decomp::TensorDecomposition{T, N}, tensor::AbstractArray{T, N}) where {T,N} =
    rel_residue(compose(decomp), tensor)

"""
Returns relative tensor decomposition error, NaN if not available.
"""
rel_residue(decomp::TensorDecomposition) = get(decomp.props, :rel_residue, NaN)

_set_rel_residue(decomp::TensorDecomposition, error::AbstractFloat) = setindex!(decomp.props, error, :rel_residue)

_set_rel_residue(decomp::TensorDecomposition{T,N}, tensor::AbstractArray{T,N}) where {T,N} =
    _set_rel_residue(decomp, rel_residue(decomp, tensor))


text:
using Test

@time @testset "private method test" begin include("util_test.jl") end
@time @testset "Sound play test" begin include("play_test.jl") end


text:
@doc raw"""
    probNFE(neuralFieldEquation)

Pre-process the neural field equation in order to be computed by `solveNFE`

`probNFE` receives as input the parameters and functions wrapped in the structures `Input1D` and `Input2D`.
The arguments of `Input1D` and `Input2D` are in the following order:
# Input1D/Input2D arguments
- `α  :: AbstractFloat`
- `v  :: AbstractFloat`
- `V0 :: fV0`: Can be a number or a function
- `L  :: Number`
- `N  :: Integer`
- `T  :: AbstractFloat`
- `n  :: Integer`
- `I  :: fI`: External input function
- `K  :: fK`: Connectivity function
- `S  :: fS`: Firing rate function

Declaration of `I`, `K` and `S` functions:
- The function I must have the arguments: in 1D - `I(x,t)`, in 2D - `I(x,y,t)`
- The function K must have the arguments: in 1D - `K(x)`, in 2D - `K(x,y)`
- The function S must have the arguments: in 1D/2D - `S(V)`

# Examples
```julia-repl
julia> # 2D example
julia> I(x,y,t) = 1                      # External input
julia> K(x,y)   = exp(-x^2+y^2)          # Connectivity function
julia> S(V)     = convert(Float64,V>0.0) # Firing rate. Heavyside function H(V)

julia> α  = 1.0  # Constant decay
julia> v  = 20.0 # Finite axonal speed
julia> V0 = 0.0  # Initial condition
julia> L  = 100  # Domain length
julia> N  = 512  # Number of nodes to discretise space
julia> T  = 20.0 # Time span
julia> n  = 200  # Number of nodes to discretise time

julia> nf   = Input2D(α,v,V0,L,N,T,n,I,K,S);
julia> prob = probNFE(nf)
├─ Domain:       Ω × [0,T] = [-50.0,49.8046875]^2 × [0,20.0]
├─ Spatial step: dx   = 0.1953125
├─ Time step:    dt   = 0.1
├─ Velocity:     v    = 20.0
├─ Delay rings:  umax = 36
```

For 1D domains the procedure is the same but the inputs are wrapped using `Input1D`.
"""
function probNFE(in::Input1D)
    N  = in.N
    hN = N÷2+1 # Due the real Fourier Transform, deals with half of the dimension

    # Contruct our discretised domain
    Ω = Domain{typeof(in.L),typeof(in.T),typeof(N)}(in.L,N,in.T,in.n,in.v)

    P    = plan_rfft(zeros(N),flags=FFTW.MEASURE) # Real IFFT operator for vector of dim N
    Pinv = plan_irfft(zeros(ComplexF64,hN),N,flags=FFTW.MEASURE) # Real IFFT operator for hN vectors
    # Apply rfft:  Â = P * A (==mul!(Â,P,A)). Dim of Â is N÷2+1
    # Apply irfft: A = Pinv * Â. Do not use mul! in the irfft case

    @inbounds K = [in.kernel(i) for i in Ω.x] # Initialise K matrix
    mul!(K,Ω.dx,K)       # Multiply K by the step size due to the discretisation of the integral operator
    krings = peel(P,K,Ω) # The rfft of K computed on the several rings of delay (k_delay)

    # Initialization of past S(V) values (from t=-Tmax to t=0, where Tmax = rings*dt)
    # Each delay ring block (hN × N) stores the S(V) rfft with the corresp delayed value
    sv = Vector{ComplexF64}(undef,Ω.rings1D*hN)
    v0 = Vector{ComplexF64}(undef,hN)
    V0 = Vector{Float64}(undef,N)
    # Initialise V0 in time domain and v0 and sv in the Fourier domain
    init!(v0,sv,V0,P,in.firingRate,in.V0,Ω)

    prob = ProbOutput1D(P,Pinv,krings,V0,v0,sv,Ω,in.α,in.extInput,in.firingRate)

return prob
end

function probNFE(in::Input2D)
    N  = in.N
    hN = N÷2+1 # Due the real Fourier Transform, deals with half of the dimension

    # Contruct our discretised domain
    Ω = Domain{typeof(in.L),typeof(in.T),typeof(N)}(in.L,N,in.T,in.n,in.v)

    P    = plan_rfft(zeros(N,N),flags=FFTW.MEASURE) # Real FFT operator for NxN matrices
    Pinv = plan_irfft(zeros(ComplexF64,hN,N),N,flags=FFTW.MEASURE) # Real IFFT operator for hNxN matrices
    # Apply rfft:  Â = P * A (==mul!(Â,P,A)). Dim of Â is ((N÷2+1) × N)
    # Apply irfft: A = Pinv * Â. Do not use mul! in the irfft case

    @inbounds K = [in.kernel(i,j) for j in Ω.y, i in Ω.x] # Initialise K matrix
    mul!(K,Ω.dx*Ω.dy,K)  # Multiply K by the step size due to the discretisation of the integral operator
    krings = peel(P,K,Ω) # The rfft of K computed on the several rings of delay (k_delay)

    # Initialization of past S(V) values (from t=-Tmax to t=0, where Tmax = rings*dt)
    # Each delay ring block (rings*hN) × N stores the S(V) rfft with the corresp delayed value
    sv = Matrix{ComplexF64}(undef,Ω.rings2D*hN,N)
    v0 = Matrix{ComplexF64}(undef,hN,N)
    V0 = Matrix{Float64}(undef,N,N)
    # Initialise V0 in time domain and v0 and sv in the Fourier domain
    init!(v0,sv,V0,P,in.firingRate,in.V0,Ω)

    prob = ProbOutput2D(P,Pinv,krings,V0,v0,sv,Ω,in.α,in.extInput,in.firingRate)

    return prob
end