text:
module Leabra

######################################################################################################
# Imports
#

using Distributions
using Random
using LinearAlgebra

######################################################################################################
# Utility Functions
#

"""
Prints variables and arrays in a more structured form that's easier to read
"""
function pretty_print(variable, prefix="")
    print(prefix)
    show(IOContext(stdout, :limit => false), "text/plain", variable)
    println()
end

"""
Expands an array along a specific dimension for a specific number of times
Found here: https://stackoverflow.com/questions/62520884/expanding-array-along-a-dimension
"""
function expand(x, dim, copies)
    sz = size(x)
    rep = ntuple(d->d==dim ? copies : 1, length(sz)+1)
    new_size = ntuple(d->d<dim ? sz[d] : d == dim ? 1 : sz[d-1], length(sz)+1)
    return repeat(reshape(x, new_size), outer=rep)
end

######################################################################################################
# NXX1 Functions
#

function XX1(x::Float64)::Float64
    return x / (x + 1)
end

function XX1GainCor(x::Float64, GainCorRange::Float64, NVar::Float64, Gain::Int64, GainCor::Float64)::Float64
    gainCorFact::Float64 = (GainCorRange - (x / NVar)) / GainCorRange
    if gainCorFact < 0
        return XX1(Gain * x)
    end
    newGain::Float64 = Gain * (1.0 - GainCor * gainCorFact)
    return XX1(newGain * x)
end

struct NXX1Parameters
    Thr::Float64            # = 0.5
    Gain::Int64             # = 100
    NVar::Float64           # = 0.005
    VmActThr::Float64       # = 0.01
    SigMult::Float64        # = 0.33
    SigMultPow::Float64     # = 0.8
    SigGain::Float64        # = 3.0
    InterpRange::Float64    # = 0.01
    GainCorRange::Float64   # = 10.0
    GainCor::Float64        # = 0.1

	SigGainNVar::Float64    # = SigGain / NVar
	SigMultEff::Float64     # = SigMult * ((Gain* NVar) ^ SigMultPow)
	SigValAt0::Float64      # = 0.5 * SigMultEff
	InterpVal::Float64      # = XX1GainCor(InterpRange, GainCorRange, NVar, Gain, GainCor) - SigValAt0

    function NXX1Parameters()
        return new(0.5, 100, 0.005, 0.01, 0.33, 0.8, 3.0, 0.01, 10.0, 0.1,
                   3.0/0.005, (0.33 * ((100*0.005)^0.8)), 0.5 * (0.33 * ((100*0.005)^0.8)),
                   XX1GainCor(0.01,10.0,0.005,100,0.1) - (0.5 * (0.33 * ((100*0.005)^0.8))) )
    end
end

nxx1p = NXX1Parameters() # default set of parameters to pass in to functions

function XX1GainCor(x::Float64, xp::NXX1Parameters = nxx1p)::Float64
    gainCorFact::Float64 = (xp.GainCorRange - (x / xp.NVar)) / xp.GainCorRange
    if gainCorFact < 0
        return XX1(xp.Gain * x)
    end
    newGain::Float64 = xp.Gain * (1.0 - xp.GainCor * gainCorFact)
    return XX1(newGain * x)
end

function NoisyXX1(x::Float64, xp::NXX1Parameters = nxx1p)::Float64
	if x < 0
		return xp.SigMultEff / (1.0 + exp(-(x * xp.SigGainNVar)))
	elseif x < xp.InterpRange
		interp::Float64 = 1.0 - ((xp.InterpRange - x) / xp.InterpRange)
		return xp.SigValAt0 + interp * xp.InterpVal
	else
		return XX1GainCor(x)
    end
end

function XX1GainCorGain(x::Float64, gain::Float64, xp::NXX1Parameters = nxx1p)::Float64
	gainCorFact::Float64 = (xp.GainCorRange - (x / xp.NVar)) / xp.GainCorRange
	if gainCorFact < 0
		return XX1(gain * x)
    end
	newGain::Float64 = gain * (1.0 - xp.GainCor * gainCorFact)
	return XX1(newGain * x)
end

function NoisyXX1Gain(x::Float64, gain::Float64, xp::NXX1Parameters = nxx1p)::Float64
	if x < xp.InterpRange
		sigMultEffArg::Float64 = xp.SigMult * (gain * xp.NVar ^ xp.SigMultPow)
		sigValAt0Arg::Float64 = 0.5 * sigMultEffArg

		if x < 0 # sigmoidal for < 0
			return sigMultEffArg / (1.0 + exp(-(x * xp.SigGainNVar)))
		else # else x < interp_range
			interp::Float64 = 1.0 - ((xp.InterpRange - x) / xp.InterpRange)
			return sigValAt0Arg + interp * xp.InterpVal
        end
	else
		return XX1GainCorGain(x, gain)
    end
end

function nxx1(points, xp::NXX1Parameters = nxx1p)::Array{Float64}
    return [NoisyXX1(point, xp) for point in points]
end

######################################################################################################
# Unit Functions
#

# g = conductance
#       ge = excitatory conductance - net excitatory input to the neuron
#       gi = inhibitory conductance - net inhibitory input to the neuron
# E = driving potential
#       Ee = excitatory driving potential
# V = voltage potential
#       Vm = membrane potential
#       Vmr = reset potential = v_reset
# 0 = action potential threshold = thr = spk_thr
# I = current
# dtvm = rate constant that determines how fast membrance potential changes
#       dtvm = vm_dt = ac.Dt.VmDt = v_dt

# Summary of Neuron Equations and Normalized Parameters
# Table 2.1 From CCN 4th Ed.
# Parameter     Bio Val     Norm Val
# Time	        0.001 sec	1ms
# Current	    1x10−8 A 	10 nA
# Capacitance   1x10−12 F 	1pF
# GbarL (leak)  10 nS	    0.1
# GbarE (excite)100 nS 	    1
# ErevI (inhib)	-75mV 	    0.25
# θ (Thr)	     -50mV	    0.5
# Voltage 	    100mV	    0..2
# Conductance 	1x10−9 S	1 nS
# C (memb cap) 	281 pF	    Dt = .355
# GBarI (inhib) 100 nS	    1
# ErevL (leak) 	-70mV	    0.3
# ErevE (excite)0mV	        1
# SpikeThr	    20mV	    1.2

mutable struct Unit # rate code approximation
    act::Float64        # = 0.2           # "firing rate" of the unit
    avg_ss::Float64     # = act           # super-short time-scale activation average
    avg_s::Float64      # = act           # short time-scale activation average
    avg_m::Float64      # = act           # medium time-scale activation average
    avg_l::Float64      # = 0.1           # long time-scale average of medium-time scale (trial level) activation, used for the BCM-style floating threshold in XCAL
    g_e::Float64        # = 0.0           # net input. Asymptotically approaches g_e_raw (see cycle).
    v_m::Float64        # = 0.3           # membrane potential
    vm_eq::Float64      # = 0.3           # a version of v_m that doesn't reset with spikes
    adapt::Float64      # = 0.0           # adaptation, as in the AdEx model
    spike::Bool         # = false         # a flag that indicates spiking threshold was crossed

    # constants
    # const g_bar_e::Float64     # = 0.3
    # const g_bar_l::Float64     # = 0.3
    const g_e_dt::Float64     # = 1/1.4   # time step constant for update of 'g_e'
    const integ_dt::Float64   # = 1.0     # time step constant for integration of cycle dynamics
    const vm_dt::Float64      # = 1/3.3   # time step constant for membrane potential
    const l_dn_dt::Float64    # = 1/2.5   # time step constant for avg_l decrease   # TODO never used
    const adapt_dt::Float64   # = 1/144   # time step constant for adaptation
    const ss_dt::Float64      # = 0.5     # time step for super-short average
    const s_dt::Float64       # = 0.5     # time step for short average
    const m_dt::Float64       # = 0.1     # time step for medium-term average
    const avg_l_dt::Float64   # = 0.1     # time step for long-term average
    const avg_l_max::Float64  # = 1.5     # max value of avg_l
    const avg_l_min::Float64  # = 0.1     # min value of avg_l
    const avg_l_gain::Float64 # = 2.5

    const e_rev_e::Float64    # = 1.0     # excitatory reversal potential
    const e_rev_i::Float64    # = 0.25    # inhibitory reversal potential
    const e_rev_l::Float64    # = 0.3     # leak reversal potential
    const g_bar_l::Float64    # = 0.1     # leak conductance
    const g_bar_e::Float64    # = 1.0     # excitatory conductance
    const thr::Float64        # = 0.5     # normalized "rate threshold"
    const spk_thr::Float64    # = 1.2     # normalized spike threshold
    const vm_r::Float64       # = 0.3     # reset potential after spike
    const vm_gain::Float64    # = 0.04    # gain that voltage produces on adaptation
    const spike_gain::Float64 # = 0.00805 # effect of spikes on adaptation
    const l_up_inc::Float64   # = 0.2     # increase in avg_l if avg_m has been 'large'   # TODO never used

    function Unit()
        return new( 0.2, 0.2, 0.2, 0.2, 0.1, 0.0, 0.3, 0.3, 0.0, false,
                    1/1.4, 1.0, 1/3.3, 1/2.5, 1/144, 0.5, 0.5, 0.1, 0.1, 1.5, 0.1, 2.5, 1.0, 0.25, 0.3, 0.1, 1.0, 0.5, 1.2, 0.3, 0.04, 0.00805, 0.2)
     end
end

function rel_avg_l(u::Unit)::Float64
    return (u.avg_l - u.avg_l_min)/(u.avg_l_max - u.avg_l_min)
end

function cycle!(u::Unit, g_e_raw::Float64, g_i::Float64, xp::NXX1Parameters = nxx1p)
    # Does one Leabra cycle. Called by the layer cycle method.
    # g_e_raw = instantaneous, scaled, received input
    # g_i = fffb inhibition

    ## updating net input
    # Ge +=      DtParams.Integ * (1/ DtParams.GTau) * (GeRaw - Ge)
    u.g_e = u.g_e + u.integ_dt * u.g_e_dt * (g_e_raw - u.g_e)

    ## Finding membrane potential
    #Inet = Ge *    (Erev.E  - Vm)    + Gbar.L * (Erev.L - Vm)    +    Gi * (Erev.I - Vm) + Noise
    i_e::Float64 = u.g_e     * (u.e_rev_e - u.v_m)
    i_l::Float64 = u.g_bar_l * (u.e_rev_l - u.v_m)
    i_i::Float64 =   g_i     * (u.e_rev_i - u.v_m)
    i_net::Float64 = i_e + i_l + i_i #+ rand() # noise?

    # almost half-step method for updating v_m (adapt doesn't half step)
    v_m_half::Float64 = u.v_m + 0.5 * u.integ_dt * u.vm_dt * (i_net - u.adapt)
    i_e_h::Float64    = u.g_e     * (u.e_rev_e - v_m_half)
    i_l_h::Float64    = u.g_bar_l * (u.e_rev_l - v_m_half)
    i_i_h::Float64    =   g_i     * (u.e_rev_i - v_m_half)
    i_net_h::Float64  = i_e_h + i_l_h + i_i_h
    u.v_m    = u.v_m + u.integ_dt * u.vm_dt * (i_net_h - u.adapt)

    # new rate coded version of i_net
    i_e_r::Float64 = u.g_e     * (u.e_rev_e - u.vm_eq)
    i_l_r::Float64 = u.g_bar_l * (u.e_rev_l - u.vm_eq)
    i_i_r::Float64 =   g_i     * (u.e_rev_i - u.vm_eq)
    i_net_r::Float64 = i_e_r + i_l_r + i_i_r
    u.vm_eq = u.vm_eq + u.integ_dt * u.vm_dt * (i_net_r - u.adapt)

    # finding whether there's an action potential
    if u.v_m > u.spk_thr
        u.spike = true
        u.v_m = u.vm_r
        i_net = 0.0
    else
        u.spike = false
    end

    # finding instantaneous rate due to input
    # if Act < XX1Params.VmActThr && Vm <= X11Params.Thr:
    if u.act < xp.VmActThr && u.vm_eq <= u.thr
        # nwAct = NoisyXX1(Vm - Thr)
        nw_act::Float64 = nxx1(u.vm_eq - u.thr)[1]
    else
        ## Finding activation
        # finding threshold excitatory conductance
        # geThr = (Gi * (Erev.I -  Thr)   + Gbar.L * (Erev.L - Thr) / (Thr - Erev.E)
        g_e_thr::Float64 = (g_i * (u.e_rev_i - u.thr) + u.g_bar_l * (u.e_rev_l - u.thr) - u.adapt) / (u.thr - u.e_rev_e)
        # nwAct = NoisyXX1(Ge * Gbar.E - geThr)
        nw_act = nxx1(u.g_e * u.g_bar_e - g_e_thr)[1]
    end

    # update activity
    # y(t) = y(t − 1) + dtvm (y∗(x) − y(t − 1))
    # Act += (1 / DTParams.VmTau) * (nwAct - Act)
    # vm_dt is 1/vm_tau
    u.act = u.act + u.integ_dt * u.vm_dt * (nw_act - u.act)

    ## Updating adaptation current
    u.adapt = u.adapt + u.integ_dt * (u.adapt_dt * (u.vm_gain * (u.v_m - u.e_rev_l)  - u.adapt) + u.spike * u.spike_gain)

    ## updating averages
    u.avg_ss = u.avg_ss + u.integ_dt * u.ss_dt * (u.act    - u.avg_ss)
    u.avg_s  =  u.avg_s + u.integ_dt *  u.s_dt * (u.avg_ss - u.avg_s)
    u.avg_m  =  u.avg_m + u.integ_dt *  u.m_dt * (u.avg_s  - u.avg_m)
end

function clamped_cycle!(u::Unit, input::Float64)
    # This function performs one cycle of the unit when its activty
    # is clamped to an input value. The activity is set to be equal
    # to the input, and all the averages are updated accordingly.

    ## Clamping the activty to the input
    u.act = input
    ## updating averages
    u.avg_ss = u.avg_ss + u.integ_dt * u.ss_dt * (u.act    - u.avg_ss)
    u.avg_s  =  u.avg_s + u.integ_dt * u.s_dt  * (u.avg_ss - u.avg_s)
    u.avg_m  =  u.avg_m + u.integ_dt * u.m_dt  * (u.avg_s  - u.avg_m)
end

function updt_avg_l!(u::Unit)
    # This function updates the long-term average 'avg_l'

    # AvgL += (1 / Tau) * (Gain * AvgM - AvgL); AvgL = MAX(AvgL, Min)
    # Tau = 10, Gain = 2.5 (this is a key param -- best value can be lower or higher) Min = .2
    u.avg_l = u.avg_l + u.avg_l_dt * (u.avg_l_gain * u.avg_m - u.avg_l)
    u.avg_l = max(u.avg_l, 0.2) #note does the 0.2 nullify u.avg_l_min?

    # matlab code
    # if u.avg_m > 0.2
    #     u.avg_l = u.avg_l + u.avg_l_dt*(u.avg_l_max - u.avg_m)
    # else
    #     u.avg_l = u.avg_l + u.avg_l_dt*(u.avg_l_min - u.avg_m)
    # end
end

function reset!(u::Unit, random::Bool=false)
    # This function sets the activity to a random value, and sets
    # all activity time averages equal to that value.
    # Used to begin trials from a random stationary point.
    if random
        u.act = rand(Uniform(0.05,0.95))
    else
        u.act = 0.0
    end
    u.avg_ss = u.act
    u.avg_s = u.act
    u.avg_m = u.act
    u.avg_l = u.act
    u.g_e   = 0.0
    u.v_m   = 0.3
    u.vm_eq = 0.3
    u.adapt = 0.0
    u.spike = 0.0
end

######################################################################################################
#  Layer Functions
#

mutable struct Layer
    units::Array{Unit}
    pct_act_scale::Float64
    scaled_acts_array::Vector{Float64}
    acts_p_avg::Float64
    ge_avg::Float64
    ge_max::Float64
    maxvsage::Float64
    wt::Array{Float64, 2}
    ce_wt::Array{Float64, 2}
    N::Int64
    fbi::Float64

    const fb_dt::Float64        # = 1/1.4 # Integration constant for feedback inhibition (fb_tau=1.4)
    const ff::Float64           # = 1.0   # feedforward scaling of inhibition
    const ff0::Float64          # = 0.1   # threshold
    const fb::Float64           # = 0.5   # feedback scaling of inhibition
    const gi::Float64           # = 1.8   # inhibition multiplier
    const avg_act_dt::Float64   # = 0.01
end

Base.show(io::IO, lay::Layer) = print(io,
                                      "Layer - ", length(lay.units), " Units\n",
                                      "pct_act_scale: ",lay.pct_act_scale,"\n",
                                      "acts_p_avg: ",   lay.acts_p_avg,"\n",
                                      "ge_avg: ",    lay.ge_avg,"\n",
                                      "wt: ",           lay.wt,"\n",
                                      "ce_wt: ",        lay.ce_wt,"\n",
                                      "N: ",            lay.N,"\n",
                                      "fbi: ",          lay.fbi,"\n",
                                      "ff: ",           lay.ff,"\n",
                                      "ff0: ",          lay.ff0,"\n",
                                      "fb: ",           lay.fb,"\n",
                                      "fb_dt: ",        lay.fb_dt,"\n",
                                      "gi: ",           lay.gi,"\n",
                                      "avg_act_dt: ",   lay.avg_act_dt,"\n",
                                     )

function layer(dims::Tuple{Int64, Int64} = (1,1))::Layer
    N::Int64 = dims[1]*dims[2]

    units::Array{Unit} = [Unit() for i in 1:N]
    scaled_acts_array = zeros(Float64, N)

    acts_avg::Float64  = 0.2
    avg_act_n::Float64  =  1.0
    pct_act_scale::Float64  = 1.0/(avg_act_n + 2.0)
    fb::Float64  = 0.5

    lay = Layer(units,
                pct_act_scale,
                scaled_acts_array,
                acts_avg,
                0.0, # ge_avg
                0.0, # maxvsavg
                0.0, # ge_max
                zeros(Float64, (1,1)),
                zeros(Float64, (1,1)),
                N,
                fb * acts_avg,
                1/1.4, 1.0, 0.1, fb, 1.8, 0.01
    )

    return lay
end

function acts_avg(lay::Layer)::Float64
    return mean([unit.act for unit in lay.units])
end

function activities(lay::Layer)::Array{Float64}
    # returns a vector with the activities of all units
    return transpose([unit.act for unit in lay.units])
end

function scaled_acts(lay::Layer)::Array{Float64}
    # returns a vector with the scaled activities of all units
    return [(lay.pct_act_scale * unit.act) for unit in lay.units]
end

function cycle!(lay::Layer, raw_inputs::Array{Float64}, ext_inputs::Array{Float64})
    # this function performs one Leabra cycle for the layer
    #raw_inputs = A vector of size I, where I is the total number of inputs
    #             from all layers. Each input has already been scaled
    #             by the pct_act_scale of its layer of origin and by
    #             the wt_scale_rel factor.
    #ext_inputs = A vector of size N denoting inputs that don't come
    #             from another layer, where N is the number of
    #             units in this layer. An empty vector indicates
    #             that there are no external inputs.

    ## obtaining the net inputs
    # GeRaw += Sum_(recv) Prjn.GScale * Send.Act * Wt
    ge_raw::Array{Float64} = lay.ce_wt * raw_inputs  # contrast-enhanced weights are used
    if any(ext_inputs) .> 0.0
        ge_raw = ge_raw .+ ext_inputs
    end

    ## obtaining inhibition
    lay.ge_avg::Float64 = mean(ge_raw)

    # it seems for a majority of the time, this won't be done, but it's here just in case
    if lay.maxvsage > 0.0
        for value in ge_raw
            lay.ge_max = 0.0 # TODO does this go here?
            if value > lay.ge_max
                lay.ge_max = value
            end
        end
        # ffNetin = avgGe + FFFBParams.MaxVsAvg * (maxGe - avgGe)
        ffNetin::Float64 = lay.ge_avg + lay.maxvsage * (lay.ge_max - lay.ge_avg)
    else
        ffNetin = lay.ge_avg
    end

    # ffi = FFFBParams.FF * MAX(ffNetin - FFBParams.FF0, 0)
    ffi::Float64 = lay.ff * max(ffNetin - lay.ff0, 0)

    # fbi += (1 / FFFBParams.FBTau) * (FFFBParams.FB * avgAct - fbi
    lay.fbi = lay.fbi + lay.fb_dt * (lay.fb * acts_avg(lay) - lay.fbi)
    # Gi = FFFBParams.Gi * (ffi + fbi)
    g_i::Float64 = lay.gi * (ffi + lay.fbi)

    ## calling the cycle method for all units
    for i::Int64 in 1:lay.N
        cycle!(lay.units[i], ge_raw[i], g_i)
    end
end

function clamped_cycle!(lay::Layer, input::Array{Float64})
    # sets all unit activities equal to the input and updates all
    # the variables as in the cycle function.
    # input = vector specifying the activities of all units
    for i::Int64 = 1:lay.N
        clamped_cycle!(lay.units[i], input[i])
    end
    lay.fbi = lay.fb * acts_avg(lay)
end

function averages(lay::Layer)::Tuple{Array{Float64}, Array{Float64}, Array{Float64}}
    # Returns the s,m,l averages in the layer as vectors.
    avg_s::Array{Float64} = [unit.avg_s for unit in lay.units]
    avg_m::Array{Float64} = [unit.avg_m for unit in lay.units]
    avg_l::Array{Float64} = [unit.avg_l for unit in lay.units]
    return (avg_s, avg_m, avg_l)
end

function rel_avg_l(lay::Layer)::Array{Float64}
    # Returns the relative values of avg_l. These are the dependent
    # variables rel_avg_l in all units used in latest XCAL
    return [rel_avg_l(unit) for unit in lay.units]
end

function updt_avg_l!(lay::Layer)
    # updates the long-term average (avg_l) of all the units in the
    # layer. Usually done after a plus phase.
    for i in 1:lay.N
        updt_avg_l!(lay.units[i])
    end
end

function updt_long_avgs!(lay::Layer)
    # updates the acts_p_avg and pct_act_scale variables.
    # These variables update at the end of plus phases instead of
    # cycle by cycle.
    # This version assumes full connectivity when updating
    # pct_act_scale. If partial connectivity were to be used, this
    # should have the calculation in WtScaleSpec::SLayActScale, in
    # LeabraConSpec.cpp
    lay.acts_p_avg = lay.acts_p_avg + lay.avg_act_dt * (acts_avg(lay) - lay.acts_p_avg)
    # r_avg_act_n::Float64 = max(round(lay.acts_p_avg * lay.N), 1.0)
    lay.pct_act_scale = 1.0/(max(round(lay.acts_p_avg * lay.N), 1.0) + 2.0)
end

function reset!(lay::Layer)
    # This function sets the activity of all units to random values,
    # and all other dynamic variables are also set accordingly.
    # Used to begin trials from a random stationary point.
    # The activity values may also be set to zero (see unit.reset)
    for i in 1:lay.N
        reset!(lay.units[i])
    end
end

######################################################################################################
#  Network Functions
#

mutable struct Network
    layers::Array{Layer}
    connections::Array{Float64, 2}
    n_lays::Int64
    n_units::Int64
    lrate::Float64

    const avg_l_lrn_max::Float64    # = 0.01    # max amount of "BCM" learning in XCAL
    const avg_l_lrn_min::Float64    # = 0.0     # min amount of "BCM" learning in XCAL
    const lrn_m::Float64            # = 0.1     # proportion of medium to short term avgs. in XCAL
    const m_lrn::Float64            # = 1       # proportion of error-driven learning in XCAL
    const d_thr::Float64            # = 0.0001  # threshold for XCAL "check mark" function
    const d_rev::Float64            # = 0.1     # reversal value for XCAL "check mark" function
    const off::Float64              # = 1.0     # 'offset' in the SIG function for contrast enhancement
    const gain::Float64             # = 6.0     # gain in the SIG function for contrast enhancement
end

function network(dim_layers::Vector{Tuple{Int64, Int64}}, connections::Matrix{Float64}, w0::Array{Union{Array{Float64}, Float64}, 2})::Network
    # TODO get the w0 parameter down to a single type.  A union is better than an Any, but it could be improved.  However, this may require reworking
    #      the algorithm.
    # constructor to the network class.
    # dim_layers = 1D array (vector). dim_layers[i] is a vector [j,k],
    #              where j is the number of rows in layer i, and k is
    #              the number of columns in layer i.
    # connections = 2D array (matrix). If layer j sends projections to
    #               layer i, then connections[i,j] = c > 0; 0 otherwise
    # w0 = 2D array. w0[i,j] is the weight matrix with the
    #      initial weights for the connections from layer j to i

    n_lay::Int64 = length(dim_layers)  # number of layers
    (nrc::Int64, ncc::Int64) = size(connections)

    @assert nrc == ncc "Non-square connections matrix in network constructor"
    @assert nrc == n_lay "Number of layers inconsistent with connection matrix in network constructor"
    @assert sum(size(w0) .== size(connections)) >= 2 "Inconsistent dimensions between initial weights and connectivity specification in network constructor"
    @assert all(connections .>= 0) "Negative projection strengths between layers are not allowed in connections matrix"

    net = Network(Array{Layer, 1}(undef, n_lay),
                  zeros(Float64, (n_lay, n_lay)),
                  n_lay,
                  0,   # counting number of units in the network
                  0.1, # set default learning rate
                  0.01, 0.0, 0.1, 1, 0.0001, 0.1, 1.0, 6.0 #constants
    )

    ## Normalizing the rows of 'connections' so they add to 1
    row_sums::Array{Float64,2} = sum(connections, dims=2)
    for row::Int64 in 1:n_lay
        if row_sums[row] > 0
            connections[row,:] = connections[row,:] ./ row_sums[row]
        end
    end

    for i::Int64 in 1:n_lay
        net.layers[i] = layer(dim_layers[i])
        net.n_units = net.n_units + net.layers[i].N
    end

    ## Second test of argument dimensions
    for i::Int64 in 1:n_lay
        for j::Int64 in 1:n_lay
            if w0[i,j] == 0.0
                if connections[i,j] > 0.0
                    throw("Connected layers have no connection weights")
                end
            else
                if connections[i,j] == 0.0
                    throw("Non-empty weight matrix for unconnected layers")
                end

                (r,c) = size(w0[i,j])
                if net.layers[i].N != r || net.layers[j].N != c
                    throw("Initial weigths are inconsistent with layer dimensions")
                end
            end
        end
    end
    net.connections = connections # assigning layer connection matrix

    ## Setting the inital weights for each layer
    # first find how many units project to the layer in all the network
    lay_inp_size::Array{Int64} = zeros(1,n_lay)
    for i::Int64 in 1:n_lay
        for j::Int64 in 1:n_lay
            if connections[i,j] > 0  # if layer j projects to layer i
                lay_inp_size[i] = lay_inp_size[i] + net.layers[j].N
            end
        end
    end

    # add the weights for each entry in w0 as a group of columns in wt
    for i in 1:n_lay
        net.layers[i].wt = zeros(Float64, (net.layers[i].N, lay_inp_size[i]))
        index::Int64 = 1
        for j::Int64 = 1:n_lay
            if connections[i,j] > 0  # if layer j projects to layer i
                nj::Int64 = net.layers[j].N
                net.layers[i].wt[:,index:index+nj-1] = w0[i,j]
                index = index + nj
            end
        end
    end

    # set the contrast-enhanced version of the weights
    for lay::Int64 in 1:n_lay
        net.layers[lay].ce_wt = 1 ./ (1 .+ (net.off * (1 .- net.layers[lay].wt) ./ net.layers[lay].wt) .^ net.gain)
    end

    return net
end

function XCAL_learn!(net::Network)
    # XCAL_learn() applies the XCAL learning equations in order to
    # modify the weights in the network. This is typically done at
    # the end of a plus phase. The equations used come from:
    # https://github.com/emer/leabra

    ## updating the long-term averages
    for lay::Int64 in 1:net.n_lays
        updt_avg_l!(net.layers[lay])
    end

    ## Extracting the averages for all layers
    avg_s = Array{Vector{Float64}, 1}(undef, net.n_lays)
    avg_m = Array{Vector{Float64}, 1}(undef, net.n_lays)
    avg_l = Array{Vector{Float64}, 1}(undef, net.n_lays)
    avg_s_eff = Array{Vector{Float64}, 1}(undef, net.n_lays)

    for lay::Int64 in 1:net.n_lays
        (avg_s[lay], avg_m[lay], avg_l[lay]) = averages(net.layers[lay]) # layer.averages() function
        # AvgSLrn = (1-LrnM) * AvgS + LrnM * AvgM
        avg_s_eff[lay] = net.lrn_m * avg_m[lay] + (1 - net.lrn_m) * avg_s[lay]
    end

    ## obtaining avg_l_lrn
    avg_l_lrn = Array{Vector{Float64}, 1}(undef, net.n_lays)
    for lay::Int64 in 1:net.n_lays
        # AvgLLrn = ((Max - Min) / (Gain - Min)) * (AvgL - Min)
        avg_l_lrn[lay] = net.avg_l_lrn_min .+ rel_avg_l(net.layers[lay]) .* (net.avg_l_lrn_max - net.avg_l_lrn_min)
    end

    ## For each connection matrix, calculate the intermediate vars.
    srs = Array{Array{Float64, 2}, 2}(undef, (net.n_lays,net.n_lays)) # srs[i,j] = matrix of short-term averages
                                                                        # where the rows correspond to the
                                                                        # units of the receiving layer, columns
                                                                        # to units of the sending layer.
    srm = Array{Array{Float64, 2}, 2}(undef, (net.n_lays,net.n_lays)) # ditto for avg_m
    for rcv::Int64 in 1:net.n_lays
        for snd::Int64 in rcv:net.n_lays
            # notice we only calculate the 'upper triangle' of the
            # arrays because of symmetry
            if net.connections[rcv,snd] > 0 || net.connections[snd,rcv] > 0
                # srs = Send.AvgSLrn * Recv.AvgSLrn
                srs[rcv,snd] = avg_s_eff[rcv] * transpose(avg_s_eff[snd])
                # srm = Send.AvgM * Recv.AvgM
                srm[rcv,snd] = avg_m[rcv] * transpose(avg_m[snd])

                if snd != rcv # using symmetry
                    srs[snd,rcv] = transpose(srs[rcv,snd])
                    srm[snd,rcv] = transpose(srm[rcv,snd])
                end
            end
        end
    end

    ## calculate the weight changes
    dwt = Array{Array{Float64, 2}, 2}(undef, (net.n_lays,net.n_lays)) # dwt[i,j] is the matrix of weight changes
                                                                      # for the weights from layer j to i
    for rcv::Int64 in 1:net.n_lays
        for snd::Int64 in 1:net.n_lays
            if net.connections[rcv,snd] > 0
                sndN::Int64 = net.layers[snd].N
                # dwt = XCAL(srs, srm) + Recv.AvgLLrn * XCAL(srs, Recv.AvgL)
                temp_dwt::Array{Float64, 2} = ( net.m_lrn .* xcal(net, srs[rcv,snd], srm[rcv,snd]) .+ ((expand(avg_l_lrn[rcv], 2, size(srs[rcv,snd])[2])) .* xcal(net, srs[rcv,snd], transpose(repeat(transpose(avg_l[rcv]), sndN)))))
                # DWt = Lrate * dwt
                dwt[rcv,snd] = net.lrate .* temp_dwt
            end
        end
    end

    ## update weights (with weight bounding)
    for rcv::Int64 in 1:net.n_lays
        # DWt = 0
        DW = zeros(Float64, (net.n_lays,net.n_lays))
        isempty_DW::Bool = true
        for snd::Int64 in 1:net.n_lays
            if net.connections[rcv,snd] > 0
                if isempty_DW
                    DW = dwt[rcv,snd]
                else
                    DW = hcat(DW, dwt[rcv,snd])
                end
                isempty_DW = false
            end
        end
        if !isempty_DW
            # Here's the weight bounding part, as in the CCN book
            # DWt *= (DWt > 0) ? Wb.Inc * (1-LWt) : Wb.Dec * LWt
            # LWt += DWt
            idxp::BitArray = net.layers[rcv].wt .> 0
            idxn::BitArray = .!idxp # maps the "!" function onto the BitArray
            net.layers[rcv].wt[idxp] = net.layers[rcv].wt[idxp] .+ (1.0 .- net.layers[rcv].wt[idxp]) .* DW[idxp]
            net.layers[rcv].wt[idxn] = net.layers[rcv].wt[idxn] .+ net.layers[rcv].wt[idxn] .* DW[idxn]
        end
    end

    ## set the contrast-enhanced version of the weights
    for lay::Int64 in 1:net.n_lays
        # Wt = SIG(LWt)
        # SIG(w) = 1 / (1 + (Off * (1-w)/w)^Gain)
        net.layers[lay].ce_wt = 1.0 ./ (1.0 .+ (net.off .* (1.0 .- net.layers[lay].wt) ./ net.layers[lay].wt) .^ net.gain)
    end
end

function updt_long_avgs!(net::Network)
    # updates the acts_p_avg and pct_act_scale variables for all layers.
    # These variables update at the end of plus phases instead of
    # cycle by cycle. The avg_l values are not updated here.
    # This version assumes full connectivity when updating
    # pct_act_scale. If partial connectivity were to be used, this
    # should have the calculation in WtScaleSpec::SLayActScale, in
    # LeabraConSpec.cpp
    for lay in 1:net.n_lays
        updt_long_avgs!(net.layers[lay])
    end
end

function m1(net::Network)::Float64
    # obtains the m1 factor: the slope of the left-hand line in the
    # "check mark" XCAL function. Notice it includes the negative
    # sign.
    return (net.d_rev - 1.0) / net.d_rev
end

function xcal(net::Network, x, th)::Array{Float64, 2}
    @assert size(x) == size(th)
    # this function implements the "check mark" function in XCAL.
    # x = an array of abscissa values.
    # th = an array of threshold values, same size as x

    # XCAL(x, th) = (x < DThr) ? 0 : (x > th * DRev) ? (x - th) : (-x * ((1-DRev)/DRev))
    # DThr = 0.0001, DRev = 0.1 defaults, and x ? y : z terminology is C syntax for: if x is true, then y, else z

    # python approach
    # if x < net.d_thr
    #     return 0
    # elseif x > (net.d_rev * th)
    #     return (x - th)
    # else
    #     return (-x * ((1 - net.d_rev)/net.d_rev))
    # end

    f::Array{Float64, 2} = zeros(size(x))
    idx1::BitArray = (x .> net.d_thr) .& (x .< (net.d_rev * th))
    idx2::BitArray = x .>= (net.d_rev * th)

    f[idx1] = m1(net) * x[idx1]
    f[idx2] = x[idx2] - th[idx2]
    return f
end

function reset!(net::Network)
    # This function sets the activity of all units to random values,
    # and all other dynamic variables are also set accordingly.
    # Used to begin trials from a random stationary point.
    for lay = 1:net.n_lays
        reset!(net.layers[lay])
    end
end

function set_weights!(net::Network, w::Array{Matrix{Float64}, 2})
    # This function receives an array w, which is like the
    # array w0 in the constructor: w[i,j] is the weight matrix with
    # the initial weights for the connections from layer j to layer
    # i. The weights are set to the values of w.
    # This whole function is a slightly modified copypasta of the
    # constructor.

    ## First we test the dimensions of w
    @assert sum(size(w) == size(net.connections)) >= 2 "Inconsistent dimensions between weights and connectivity specification in set_weights"

    for i::Int64 in 1:net.n_lays
        for j::Int64 in 1:net.n_lays
            if all(w[i,j] .== 0.0) # if isempty(w[i,j])
                if net.connections[i,j] > 0
                    throw("Connected layers have no connection weights")
                end
            else
                if net.connections[i,j] == 0
                    throw("Non-empty weight matrix for unconnected layers")
                end
                (r::Int64, c::Int64) = size(w[i,j])
                if net.layers[i].N != r || net.layers[j].N != c
                    throw("Initial weights are inconsistent with layer dimensions")
                end
            end
        end
    end

    ## Now we set the weights
    # first find how many units project to the layer in all the network
    lay_inp_size = zeros(1,net.n_lays)
    for i::Int64 in 1:net.n_lays
        for j::Int64 in 1:net.n_lays
            if net.connections[i,j] > 0  # if layer j projects to layer i
                lay_inp_size[i] = lay_inp_size[i] + net.layers[j].N
            end
        end
    end

    # add the weights for each entry in w0 as a group of columns in wt
    for i::Int64 in 1:net.n_lays
        net.layers[i].wt = zeros(net.layers[i].N,lay_inp_size[i])
        index::Int64 = 1
        for j in 1:net.n_lays
            if net.connections[i,j] > 0  # if layer j projects to layer i
                nj::Int64 = net.layers[j].N
                net.layers[i].wt[:,index:index+nj-1] = w[i,j]
                index = index + nj
            end
        end
    end

    # set the contrast-enhanced version of the weights
    for lay::Int64 in 1:net.n_lays
        net.layers[lay].ce_wt = 1 ./ (1 .+ (net.off .* (1 .- net.layers[lay].wt) ./ net.layers[lay].wt) .^ net.gain)
    end
end

function get_weights(netnet::Network)::Array{Array{Float64, 2}, 2}
    # This function returns a 2D array w.
    # w[rcv,snd] contains the weight matrix for the projections from
    # layer snd to layer rcv.
    w = Array{Array{Float64, 2}, 2}(undef, (net.n_lays,net.n_lays))
    for rcv::Int64 in 1:net.n_lays
        idx1::Int64 = 1 # column where the weights from layer 'snd' start
        for snd::Int64 in 1:net.n_lays
            if net.connections[rcv,snd] > 0
                Nsnd::Int64 = net.layers[rcv].N
                w[rcv,snd] = net.layers[rcv].wt[:,idx1:idx1+Nsnd-1]
                idx1 = idx1 + Nsnd
            end
        end
    end
    return w
end

function cycle!(net::Network, inputs::Vector{Array{Float64}}, clamp_inp::Bool)
    # this function calls the cycle method for all layers.
    # inputs[i] is  a matrix that for layer i
    #          specifies the external input to each of its units.
    #          An empty matrix denotes no input to that layer.
    # clamp_inp = a binary flag. 1 -> layers are clamped to their
    #             input value. 0 -> inputs summed to ge_raw.

    ## Testing the arguments and reshaping the input
    @assert length(inputs) == net.n_lays "Number of layers inconsistent with number of inputs in network cycle"

    for inp::Int64 in 1:net.n_lays  # reshaping inputs into column vectors
        if any(inputs[inp] .> 0.0) #if ~isempty(inputs[inp])
            inputs[inp] = reshape(inputs[inp], net.layers[inp].N, 1)
        end
    end

    ## First we set all clamped layers to their input values
    clamped_lays = zeros(Bool, 1, net.n_lays)
    if clamp_inp
        for lay::Int64 in 1:net.n_lays
            if any(inputs[lay] .> 0.0) # if ~isempty(inputs[lay])
                clamped_cycle!(net.layers[lay], inputs[lay])
                clamped_lays[lay] = true
            end
        end
    end

    ## We calculate the scaled activity for all layers
    for lay::Int64 in 1:net.n_lays
        net.layers[lay].scaled_acts_array = scaled_acts(net.layers[lay])
    end

    ## For each unclamped layer, we put all its scaled inputs in one
    #  column vector, and call its cycle function with that vector
    for recv::Int64 in 1:net.n_lays
        if !clamped_lays[recv] # if the layer is not clamped
            # for each 'recv' layer we find its input vector
            long_input = zeros(Float64, 1, net.n_units)
            n_inps::Int64 = 0  # will have the # of input units to 'recv'
            n_sends::Int64 = 0 # will have the # of layers sending to 'recv'
            conns::BitArray = (net.connections[recv, :] .> 0.0)
            wt_scale_rel::Array{Float64} = net.connections[recv, conns]
            for send::Int64 in 1:net.n_lays
                if net.connections[recv, send] > 0
                    n_sends = n_sends + 1
                    long_input[(1 + n_inps):(n_inps + net.layers[send].N)] = wt_scale_rel[n_sends] .* net.layers[send].scaled_acts_array
                    n_inps = n_inps + net.layers[send].N
                end
            end
            # now we call 'cycle'
            cycle!(net.layers[recv], long_input[1:n_inps], inputs[recv])
        end
    end

end

######################################################################################################
# Higher level functions for building networks
#

function build_network(dim_lays::Vector{Tuple{Int64, Int64}}, connections::Matrix{Float64})::Network
    n_lays::Int64 = length(dim_lays)
    n_units::Array{Int64} = zeros(Int64,1,n_lays)
    for i::Int64 in 1:n_lays
        n_units[i] = dim_lays[i][1] * dim_lays[i][2]
    end
    w0 = Array{Union{Array{Float64}, Float64}, 2}(undef, (n_lays,n_lays))
    for rcv::Int64 in 1:n_lays
        for snd::Int64 in 1:n_lays
            if connections[rcv,snd] > 0
                w0[rcv,snd] = 0.3 .+ 0.4*rand(Uniform(),n_units[rcv],n_units[snd])
                # w0[rcv,snd] = rand(Uniform(),n_units[rcv],n_units[snd])
            else
                w0[rcv,snd] = 0.0
            end
        end
    end
    net = Leabra.network(dim_lays, connections, w0)

    return net
end

function create_random_inputs(n_inputs::Int64)::Array{Array{Float64}, 2}
    inputs = Array{Array{Float64}, 2}(undef, (n_inputs,2))
    for i::Int64 in 1:n_inputs
        inputs[i,1] = rand(Binomial(1,0.5),(n_inputs,n_inputs))
        inputs[i,2] = inputs[i,1];
    end

    return inputs
end

function create_patterns(n_inputs::Int64)::Array{Array{Float64}, 2}
    inputs = Array{Array{Float64}, 2}(undef, (n_inputs,2))
    inputs[1,1] = repeat([0 0 1 0 0],5,1);   # vertical line
    inputs[2,1] = [1 1 1 1 1;zeros(4,5)];    # horizontal line
    inputs[3,1] = [0 0 0 0 1;0 0 0 1 0;0 0 1 0 0;0 1 0 0 0; 1 0 0 0 0]; # diagonal 1
    inputs[4,1] = reverse(inputs[3,1], dims=2)
    inputs[5,1] =  1.0 .* (!=(0.0).(inputs[3,1]) .|| !=(0.0).(inputs[4,1]))   # two diagonals
    for i::Int64 in 1:n_inputs # outputs are the same as inputs (an autoassociator)
        inputs[i,1] = 0.01 .+ 0.98 .* inputs[i,1];
        inputs[i,2] = inputs[i,1];
    end
    return inputs
end

function clamp_data(data::Matrix{Array{Float64}}, num_layers::Int64, layers_to_clamp::Vector{Int64}, selected_data::Int64)::Array{Array{Float64}}
    resulting_data = Array{Array{Float64}}(undef, num_layers)
    output_layer::Bool = false

    for layer::Int64 in 1:num_layers
        if layer == num_layers
            output_layer = true
        end

        if layer in layers_to_clamp
            if output_layer
                resulting_data[layer] = data[selected_data,2]
            else
                resulting_data[layer] = data[selected_data,1]
            end
        else
            resulting_data[layer] = []
        end
    end

    return resulting_data
end

function train_network!(net::Network, inputs::Matrix{Array{Float64}}, n_epochs::Int64)::Array{Float64}
    n_minus::Int64 = 50 # number of minus cycles per trial
    n_plus::Int64 = 25  # number of plus cycles per trial
    n_trials::Int64 = size(inputs)[1]

    lrate_sched::Vector{Float64} = collect(LinRange(0.8, 0.2, n_epochs))
    errors = zeros(Float64, n_epochs,n_trials)
    for epoch::Int64 in 1:n_epochs
        randomized_input::Vector{Int64} = randperm(n_trials)::Vector{Int64}
        net.lrate = lrate_sched[epoch]
        for trial::Int64 in 1:n_trials
            Leabra.reset!(net)
            selected_data::Int64 = randomized_input[trial] # pick a random set of data

            # Minus phase
            # Clamp input layer
            # println(typeof(net.n_lays))
            sel_data_array::Array{Array{Float64}} = clamp_data(inputs, net.n_lays, [1], selected_data)
            for _ in 1:n_minus
                Leabra.cycle!(net, sel_data_array, true)
            end
            outs::Matrix{Float64} = Leabra.activities(net.layers[net.n_lays])

            # Plus phase
            # Clamp input and output layer
            sel_data_array = clamp_data(inputs, net.n_lays, [1,net.n_lays], selected_data)
            for _ in 1:n_plus
                Leabra.cycle!(net, sel_data_array, true);
            end
            Leabra.updt_long_avgs!(net)

            # Learning
            Leabra.XCAL_learn!(net)

            # Errors
            errors[epoch, selected_data] = 1.0 - sum(outs .* transpose(inputs[selected_data, 2][:])) / ( norm(outs) * norm(transpose(inputs[selected_data, 2][:])) );
        end
    end

    return errors
end

######################################################################################################
# Hacky Export All Function
#

for n in names(@__MODULE__; all=true)
    if Base.isidentifier(n) && n ∉ (Symbol(@__MODULE__), :eval, :include)
        @eval export $n
    end
end

end # module Leabra


text:
using ForwardDiff
using Distributions

# Real

x_real = randn(5)

dists = [Normal(0, 1)]

for dist in dists

    f(x::Vector) = sum(logpdf.(Ref(dist), x))

    ForwardDiff.gradient(f, x_real)

end

# Postive

x_positive = randn(5).^2

dists = [Gamma(2, 3)]

for dist in dists

    f(x::Vector) = sum(logpdf.(Ref(dist), x))

    g = x -> ForwardDiff.gradient(f, x)

end

ext: jl


text:
using Distributed
@everywhere using Mamba

## Data
asthma = Dict{Symbol, Any}(
  :y =>
    [210 60 0 1  1
     88 641 0 4 13
     1    0 0 0  1],
  :M =>
    [272, 746, 2]
)

## Model Specification
model = Model(

  y = Stochastic(2,
    (M, q) ->
      MultivariateDistribution[
        Multinomial(M[i], vec(q[i, :]))
        for i in 1:length(M)
      ],
    false
  ),

  q = Stochastic(2,
    M ->
      MultivariateDistribution[
        Dirichlet(ones(5))
        for i in 1:length(M)
      ],
    true
  )

)

## Initial Values
inits = [
  Dict{Symbol, Any}(
    :y => asthma[:y],
    :q => vcat([rand(Dirichlet(ones(5)))' for i in 1:3]...)
  )
  for i in 1:3
]

## Sampling Scheme
scheme = [SliceSimplex(:q)]
setsamplers!(model, scheme)

## MCMC Simulations
sim = mcmc(model, asthma, inits, 10000, burnin=2500, thin=2, chains=3)
describe(sim)


text:
"""
Loads dependencies and bootstraps a Genie app.
"""
module Genie

push!(LOAD_PATH, joinpath(Pkg.dir("Genie"), "src"))

include(joinpath(Pkg.dir("Genie"), "src", "configuration.jl"))

include(joinpath(Pkg.dir("Genie"), "src", "constants.jl"))
if haskey(ENV, "GENIE_ENV") && isfile(abspath(joinpath(ENV_PATH, ENV["GENIE_ENV"] * ".jl")))
  isdefined(:config) || include(abspath(joinpath(ENV_PATH, ENV["GENIE_ENV"] * ".jl")))
  isfile(abspath(joinpath(CONFIG_PATH, "app.jl"))) && include(abspath(joinpath(CONFIG_PATH, "app.jl")))

  const IS_IN_APP = true
else
  const IS_IN_APP = false
  const config = Configuration.Settings(app_env = Configuration.DEV)
end

const SEARCHLIGHT_ON = isdir(Pkg.dir("SearchLight")) && IS_IN_APP ? true : false

export IS_IN_APP, SEARCHLIGHT_ON

isfile(abspath(joinpath(CONFIG_PATH, "plugins.jl"))) && include(abspath(joinpath(CONFIG_PATH, "plugins.jl")))

push!(LOAD_PATH,  joinpath(Pkg.dir("Genie"), "src", "cache_adapters"),
                  joinpath(Pkg.dir("Genie"), "src", "session_adapters"),
                  RESOURCES_PATH, HELPERS_PATH)

include(joinpath(Pkg.dir("Genie"), "src", "genie_types.jl"))
include(joinpath(Pkg.dir("Genie"), "src", "file_templates.jl"))
include(joinpath(Pkg.dir("Genie"), "src", "generator.jl"))
include(joinpath(Pkg.dir("Genie"), "src", "REPL.jl"))

using Macros, Logger, App, Commands, AppServer, Millboard, Renderer

SEARCHLIGHT_ON && eval(:(using SearchLight))

if IS_IN_APP
  @eval parse("@dependencies")
  include(joinpath(Pkg.dir("Genie"), "src", "deprecations.jl"))
end

"""
    run() :: Void

Runs the Genie app by parsing the command line args and invoking the corresponding actions.
"""
function run() :: Void
  Commands.execute(Genie.config)

  nothing
end

end


text:
using Images
using ArtBooster

img = load(Pkg.dir("ArtBooster") * "/test/julia-logo.jpg")

param = Dict("max_depth" => 2,
             "eta" => 1.,
             "objective" => "reg:linear")

boostimage(img, param, iterations = 1, res_x = 200, res_y = 200)


text:
const POSITIVES = 1:∞

struct Skip{Synthesizer}
    synthesizer::Synthesizer
    time::TIME
end

function make_series(skip::Skip, sample_rate)
    make_series(skip.synthesizer, sample_rate)[
        (round(Int, skip.time * sample_rate)+1):end
    ]
end

"""
    Map(a_function, synthesizers...)

Map `a_function` over `synthesizers`. Supports [`make_series`](@ref).

```jldoctest
julia> using AudioSchedules

julia> using Unitful: Hz

julia> first(make_series(Map(sin, Cycles(440Hz)), 44100Hz))
0.06264832417874369
```
"""
struct Map{AFunction,Synthesizers}
    a_function::AFunction
    synthesizers::Synthesizers
    function Map(a_function::AFunction, synthesizers...) where {AFunction}
        new{AFunction,typeof(synthesizers)}(a_function, synthesizers)
    end
end
export Map

"""
    make_series(synthesizer, sample_rate)

Return an iterator that will the play the `synthesizer` at `sample_rate` (with frequency units, like `Hz`).
The iterator should yield `Float64`s between -1 and 1.
Assumes that iterators will never end while they are scheduled.
"""
function make_series(a_map::Map, sample_rate)
    broadcast(
        a_map.a_function,
        map(let sample_rate = sample_rate
            function (synthesizer)
                make_series(synthesizer, sample_rate)
            end
        end, a_map.synthesizers)...,
    )
end
export make_series

"""
    Line(start_level, slope)

A line from `start_level` (unitless) with `slope` (with units per time like `1/s`). Supports
[`make_series`](@ref) and [`segments`](@ref).

```jldoctest
julia> using AudioSchedules

julia> using Unitful: Hz, s

julia> first(make_series(Line(0, 1 / s), 44100Hz))
2.2675736961451248e-5
```
"""
struct Line
    start_level::Float64
    slope::RATE
end
export Line

function make_series(line::Line, sample_rate)
    line.start_level .+ (line.slope / sample_rate) .* POSITIVES
end

"""
    segments(shape, start_level, duration, end_level)

Called for each envelope segment passed to [`push!`](@ref). Return a tuple of pairs in the form `(segment, duration)`,
where duration has units of time (like `s`), for a segment of shape `shape`.

```jldoctest
julia> using AudioSchedules

julia> using Unitful: s

julia> segments(Grow, 1, 1s, ℯ)
((Grow(1.0, 1.0 s⁻¹), 1 s),)
```
"""
function segments(::Type{Line}, start_level, duration, end_level)
    ((Line(start_level, (end_level - start_level) / duration), duration),)
end
export segments

const τ = 2 * π

# TODO: No need to cycle here; sin is periodic?
"""
    Cycles(frequency)

Cycles from 0 to 2π to repeat at a `frequency` (with frequency units, like `Hz`).
Supports [`make_series`](@ref).

```jldoctest
julia> using AudioSchedules

julia> using Unitful: Hz

julia> first(make_series(Cycles(440Hz), 44100Hz))
0.06268937721449021
```
"""
struct Cycles
    frequency::FREQUENCY
end
export Cycles

function make_series(cycles::Cycles, sample_rate)
    ((cycles.frequency / sample_rate * τ) .* POSITIVES) .% τ
end

"""
    Grow(start_level, rate)

Exponentially grow or decay from `start_level` (unitless), at a continuous `rate` (with units per time like `1/s`).
Supports [`make_series`](@ref) and [`segments`](@ref).

```jldoctest
julia> using AudioSchedules

julia> using Unitful: Hz, s

julia> first(make_series(Grow(1, 1 / s), 44100Hz))
1.0000226759940578
```
"""
struct Grow
    start_level::Float64
    rate::RATE
end
export Grow

function make_series(grow::Grow, sample_rate)
    grow.start_level .* exp(grow.rate / sample_rate) .^ POSITIVES
end

function segments(::Type{Grow}, start_level, duration, end_level)
    ((Grow(start_level, log(end_level / start_level) / duration), duration),)
end

max_stars_repo_name: bramtayl/AudioSchedulers

max_issues_repo_name: bramtayl/AudioSchedulers

max_forks_repo_name: bramtayl/AudioSchedulers


text:
import Base.Random: rand, srand
import RNG: AbstractRNG, gen_seed, split_uint

"""
```julia
AbstractXoroshiro128 <: AbstractRNG{UInt64}
```

The base abstract type for `Xoroshiro128`, `Xoroshiro128Star` and `Xoroshiro128Plus`.
"""
abstract AbstractXoroshiro128 <: AbstractRNG{UInt64}

@inline xorshift_rotl(x::UInt64, k) = (x << k) | (x >> (64 - k))

for (star, plus) in (
        (false, false),
        (false, true),
        (true, false),
    )
    rng_name = Symbol(string("Xoroshiro128", star ? "Star" : plus ? "Plus" :""))
    @eval begin
        type $rng_name <: AbstractXoroshiro128
            x::UInt64
            y::UInt64
            function $rng_name(seed::NTuple{2, UInt64}=gen_seed(UInt64, 2))
                r = new(0, 0)
                srand(r, seed)
                r
            end
        end

        $rng_name(seed::Integer) = $rng_name(split_uint(seed % UInt128))

        @inline function xorshift_next(r::$rng_name)
            $(plus ? :(p = r.x + r.y) : nothing)
            s1 = r.y $ r.x
            r.x = xorshift_rotl(r.x, 55) $ s1 $ (s1 << 14)
            r.y = xorshift_rotl(s1, 36)
            $(star ? :(r.y * 2685821657736338717) :
              plus ? :(p) : :(r.y))
        end
    end
end

srand(r::AbstractXoroshiro128, seed::Integer) = srand(r, split_uint(seed % UInt128))
function srand(r::AbstractXoroshiro128, seed::NTuple{2, UInt64}=gen_seed(UInt64, 2))
    r.x = seed[1]
    r.y = seed[2]
    xorshift_next(r)
    xorshift_next(r)
    r
end

@inline rand(r::AbstractXoroshiro128, ::Type{UInt64}) = xorshift_next(r)


text:
#using IteratorInterfaceExtensions, TableTraits, TableTraitsUtils

"""
ResultSet is the primary object that represents data returned from
reading a SAS data file.  ResultSet implements the Base.Iteration
interface as well as the IterableTables.jl interface.

*Fields*
- `columns`: a vector of columns, each being a vector itself
- `names`: a vector of column symbols
- `size`: a tuple (nrows, ncols)

*Accessors*
- `columns(::ResultSet)`
- `names(::ResultSet)`
- `size(::ResultSet)`
- `size(::ResultSet, dim::Integer)`

*Single Row/Column Indexing*
- `rs[i]` returns a tuple for row `i`
- `rs[:c]` returns a vector for column `c`

*Multiple Row/Column Indexing*
- `rs[i:j]` returns a view of ResultSet with rows between `i` and `j`
- `rs[:c...]` returns a view of ResultSet with columns specified e.g. `rs[:A, :B]`

*Cell Indexing*
- `rs[i,j]` returns a single value for row `i` column `j`
- `rs[i,:c]` returns a single value for row `i` column `c`
- Specific cell can be assigned using the above indexing methods
"""
struct ResultSet
    columns::AbstractVector{AbstractVector}
    names::AbstractVector{Symbol}
    size::NTuple{2, Int}
    ResultSet() = new([], [], (0,0))  # special case
    ResultSet(c,n,s) = new(c, n, s)
end

# exports
export columns

# accessors
columns(rs::ResultSet) = getfield(rs, :columns)
Base.names(rs::ResultSet) = getfield(rs, :names)

Base.size(rs::ResultSet) = getfield(rs, :size)
Base.size(rs::ResultSet, i::Integer) = getfield(rs, :size)[i]
Base.length(rs::ResultSet) = getfield(rs, :size)[1]

# Size displayed as a string
sizestr(rs::ResultSet) = string(size(rs, 1)) * " rows x " * string(size(rs, 2)) * " columns"

# find index for the column symbol
function symindex(rs::ResultSet, s::Symbol)
    n = findfirst(x -> x == s, names(rs))
    n == 0 && error("column symbol not found: $s")
    n
end

# Direct cell access
Base.getindex(rs::ResultSet, i::Integer, j::Integer) = columns(rs)[j][i]
Base.getindex(rs::ResultSet, i::Integer, s::Symbol) = columns(rs)[symindex(rs, s)][i]
Base.setindex!(rs::ResultSet, val, i::Integer, j::Integer) = columns(rs)[j][i] = val
Base.setindex!(rs::ResultSet, val, i::Integer, s::Symbol) = columns(rs)[symindex(rs, s)][i] = val

# Return a single row as a named tuple
Base.getindex(rs::ResultSet, i::Integer) =
    NamedTuple{Tuple(names(rs))}([c[i] for c in columns(rs)])

# Return a single column
Base.getindex(rs::ResultSet, c::Symbol) = columns(rs)[symindex(rs, c)]

# index by row range => returns ResultSet object
function Base.getindex(rs::ResultSet, r::UnitRange{Int})
    ResultSet(map(x -> view(x, r), columns(rs)), names(rs), (length(r), size(rs, 2)))
end

# index by columns => returns ResultSet object
function Base.getindex(rs::ResultSet, ss::Symbol...)
    v = Int[]
    for (idx, nam) in enumerate(names(rs))
        nam in ss && push!(v, idx)
    end
    ResultSet(columns(rs)[v], names(rs)[v], (size(rs, 1), length(v)))
end

# Each property must represent a column to satisfy Tables.jl Columns interface
Base.propertynames(rs::ResultSet) = names(rs)

function Base.getproperty(rs::ResultSet, s::Symbol)
    s ∈ names(rs) && return rs[s]
    error("Column $s not found")
end

# Iterators
Base.iterate(rs::ResultSet, i=1) = i > size(rs,1) ? nothing : (rs[i], i+1)

# Display ResultSet object
function Base.show(io::IO, rs::ResultSet)
    println(io, "SASLib.ResultSet (", sizestr(rs), ")")
    max_rows = 5
    max_cols = 10
    n = min(size(rs, 1), max_rows)
    m = min(size(rs, 2), max_cols)
    (n < 1 || m < 1) && return
    print(io, "Columns ")
    for i in 1:m
        i > 1 && print(io, ", ")
        print(io, i, ":", names(rs)[i])
    end
    m < length(names(rs)) && print(io, " …")
    println(io)
    for i in 1:n
        print(io, i, ": ")
        for j in 1:m
            j > 1 && print(io, ", ")
            print(io, columns(rs)[j][i])
        end
        println(io)
    end
    n < size(rs, 1) && println(io, "⋮")
end


text:
using Base.Threads: @spawn

function threads(xs)
    n = length(xs)
    chunks = view(xs, 1:n÷2), view(xs, n÷2+1:n)
    p = bufferfrom([0.0, 0.0])
    @sync begin
        for i = 1:2
            @spawn begin
                s = zero(eltype(chunks[i]))
                for j = 1:length(chunks[i])
                    s += chunks[i][j]
                end
                p[i] = s
            end
        end
    end
    return p[1] + p[2]
end

@testset "threads" begin
    @test gradient(threads, [1, 2, 3, 4]) == ([1, 1, 1, 1],)
end


text:
# main.jl

include("FFGmain.jl")
include("lattices.jl")

# define a model
dims = (dl[i],) # dimension of the lattice
lattice = Lattices.get_SquareLattice(dims) # make a square lattice
# observable = FFG.observe_Swap()
filling = 1 # setting filling ≢ 1 means there are holes.

model = Dict(
    "dims" => dims,
    "lattice" => lattice,
    "filling" => filling,
    "hamiltonian" => Lattices.get_Tightbinding_Wavefunctions,
    "fermi_energy" => 0,
    )

# Policy
policy = FFG.SwapNeighborsPolicy()

# MC specifications
mc_spec = Dict("mc_steps" => 1E4*prod(dims),
                "mc_warmup_steps" => 100*prod(dims),
                "sample_interval" => 4*prod(dims))

# create a markov chain object representing our gutzwiller projected
# free fermion model
dg = FFG.DoubleGutzwillerChain()

# initialize the chain: sets up the initial state
FFG.init_Chain!(dg,
    model=model,observable=FFG.MeasureSwap,
    policy=policy,mc_spec=mc_spec)

FFG.runMC!(dg)

# print(FFG.get_Data(gutzwiller_chain))
# print("\n")
# print("Acceptance Ratio: ", FFG.get_Diagnostics(gutzwiller_chain)["acceptance_ratio"])
# print("\n")
print(FFG.get_Data(dg), "\n")


text:
"""
Returns list of similarly named colors, from the color_list
"""
function similarly_named_colors(
		name::AbstractString,
		color_list::AbstractDict{S} = ALL_COLORS) where {S<:AbstractString}

	levsort(name, collect(keys(color_list)))
end

struct UnknownColorError{S<:AbstractString} <:Exception
    name::S
    msg::String
end

function UnknownColorError(
		name::AbstractString,
		color_list::AbstractDict{S}) where {S<:AbstractString}

	similars = similarly_named_colors(name, color_list)
    msg=if length(similars)==0
        unescape_string("\"$name\" is not a known color name.")
    else
        suggestions = join(repr.(similars), ", ", " or ")
        unescape_string("\"$name\" is not a known color name. Perhaps you meant: $suggestions")
    end
    UnknownColorError(name, msg)
end

Base.showerror(io::IO, ex::UnknownColorError) = print(io, "UnknownColorError: $(ex.msg)")

"""
Returns a color with the given name.
If not found,a list of suggestions is provided.
"""
function named_color(name::AbstractString, color_list::AbstractDict{S} = ALL_COLORS) where {S<:AbstractString}

    if haskey(color_list, name)
        color_list[name]
    else
        throw(UnknownColorError(name, color_list))
    end
end

"""
Attempts to find a color to match any name you give it.
First by using the methods from Color.jl's @colorant_str macro,
to handle things like `colorant"#cd32cd"`
Then if that fails, does a lookup in the big list of named colors.
If that fails, then it provides suggestions for what color name you may have meant.
As an user, you basically should just use this, without worrying if a color is defined.
There are almost 3500 defined colors, and if you miss one,
then the suggestions will help you out.
(this is easier than trying to workout the name yourself)
"""
macro colorant_str(name::AbstractString)
    local col
    try
        col = parse(Colorant, name)
    catch ex
        col = named_color(name)
        :($col)
    end
    :($col)
end


text:
srand(1)

# TODO: Contribs gets updated properly when fitness(...) is called.

context("genotype.jl") do
  landscapes = [
    NK.NKLandscape(4, 1),
    NK.NKqLandscape(4, 1, 2),
    NK.NKpLandscape(4, 1, 0.90),
  ]

  for ls = landscapes
    # The ith locus is linked to the locus i places to its right, assuming a
    # periodic boundary.
    ls.links = [
      0b1000,
      0b0100,
      0b0010,
      0b0001
    ]

    context("$(ls)") do
      context("Genotype") do
        g = NK.Genotype(0b0000, ls)
        @fact g.alleles --> 0b0000
        @fact g.landscape --> ls
      end

      context("fitness(...)") do
        g = NK.Genotype(0b0000, ls)
        @fact NK.fitness(g) --> less_than_or_equal(1.0)
        @fact NK.fitness(g) --> greater_than_or_equal(0.0)
        @fact NK.fitness(g) --> roughly(NK.contribs(g) |> mean)
      end

      context("rand(::Type{Genotype}, ...)") do
        g = rand(NK.Genotype, ls)
        @fact g.alleles --> less_than_or_equal(0b1111)
        @fact g.alleles --> greater_than_or_equal(0b0000)
      end

      context("zeros(::Type{Genotype}, ...)") do
        g = zeros(NK.Genotype, ls)
        @fact g.alleles --> 0
        @fact g.landscape --> ls
      end
    end
  end
end


text:
#C:\Users\totr\.julia\dev\MyExample\src
#C:/Users/totr/.julia/dev/MyExample/src

using MyExample
using Test

@testset "MyExample.jl" begin
    # 2x + 3y
    @test my_f(2,1) == 7
    @test my_f(2,3) == 13
end


text:
using SpecialFunctions
using DataStructures
using StaticArrays
using ComplexOperations
using Memoize
using LoopVectorization

import WignerD
import ChainRulesCore
import Zygote

export γ_mn
export γ_mn_dash
export wignerdjmn_ELZOUKA
export ∂wignerdjmn_by_∂θ
export single_index_from_m_n
export get_max_single_index_from_n_max

#############################################################################################
# import "WignerD.jl" code to calculate wigner-d (https://github.com/jishnub/WignerD.jl)
# [TODO] I have already installed "WignerD" package using the following commands: # include("WignerD.jl/src/WignerD.jl") ;  import Pkg; Pkg.add("WignerD")
# include("WignerD.jl/src/WignerD.jl")
# using .WignerD

#############################################################################################
"""
    Dirac delta function
returns 1 if the two inputs are equal. Otherwise, return Zero
"""
δ(x,y) = ==(x, y)

"""
    Fuse the two indices (m,n) into a single index.
"""
function single_index_from_m_n(m, n)
    return n * (n + 1) + m
end

function get_max_single_index_from_n_max(n_max)
    return single_index_from_m_n(n_max, n_max)
end

#############################################################
"""
    defining derivative for Bessel functions.
Inside "SpecialFunctions.jl" package, the definition of derivative of Bessel functions exists. However, it returns an error when called by Zygote.
The error because the derivative of Bessel function with respect to the integer rank is not defined, and the code return error "ChainRulesCore.@thunk(error("not implemented"))"
Since adjoint calculation will never require calculation gradient with respect to the integer rank, I am defining this integer as zero or `ChainRulesCore.ZeroTangent()`
"""

#
ChainRulesCore.@scalar_rule(
    besselj(ν, x),
    (
        ChainRulesCore.ZeroTangent(), # ChainRulesCore.@thunk(error("not implemented")),
        (besselj(ν - 1, x) - besselj(ν + 1, x)) / 2
    ),
)
ChainRulesCore.@scalar_rule(
    besseli(ν, x),
    (
        ChainRulesCore.ZeroTangent(), # ChainRulesCore.@thunk(error("not implemented")),
        (besseli(ν - 1, x) + besseli(ν + 1, x)) / 2,
    ),
)
ChainRulesCore.@scalar_rule(
    bessely(ν, x),
    (
        ChainRulesCore.ZeroTangent(), # ChainRulesCore.@thunk(error("not implemented")),
        (bessely(ν - 1, x) - bessely(ν + 1, x)) / 2,
    ),
)

# Wigner-d
"""
    Wigner-d function calculated from eq. B.1
This is calculated directly, not with recurrence relation.
The Wigner-d calculated by recurrence is more numerically stable.
Use this function only as a validation for the recurrence relation, or when automatic differentiation is needed.

TODO: for large s,m,n, this function is not stable due to large integers resulted from factorial calclulation and their multiplication. Here, I am using big() to avoid the overflow error. I need to implement when I should apply big(), while not affecting performance and autodifferentiation

"""
function wignerdjmn_ELZOUKA(s::I, m::I, n::I, θ::R) where {R <: Real, I <: Integer}
    # println("s=$s, m=$m, n=$n, θ=$θ")
    if θ == zero(θ) # TODO: make the zero the same type of θ. e.g., zero(θ)
        d = δ(m, n);
    elseif θ == π
        d = (-1)^(s - n) * δ(-n, m);
    else
        d = zero(θ)

        k_min = max(0, m - n)
        k_max = min(s + m, s - n)

        if k_max >= k_min
            # TODO: find a better way to detect when do we need to use big integers
                if (
                    (max(k_max, s + m - k_max, s - n - k_max, n - m + k_max, s + abs(m), s + abs(n)) >= 21) ||
                    ((factorial(s + m) * factorial(s - m) * factorial(s + n) * factorial(s - n)) < 0) ||
                    ((factorial(k_max) * factorial(s + m - k_max) * factorial(s - n - k_max) * factorial(n - m + k_max)) < 0)
                )
                s = big(s); m = big(m); n = big(n)
                k_max = min(s + m, s - n)
            end
            for k in k_min:k_max
                d += (-1)^k *
                        (cos(θ / 2)^(2s - 2k + m - n) * sin(θ / 2)^(2k - m + n)) /
                        (factorial(k) * factorial(s + m - k) * factorial(s - n - k) * factorial(n - m + k))
            end
            d *= sqrt(factorial(s + m) * factorial(s - m) * factorial(s + n) * factorial(s - n))
            return convert(typeof(θ), d)
        else # wigner-d is zero if there is any negative factorial
            return zero(θ)
        end
    end

    return convert(typeof(θ), d)
end

"""
putting non-differentiable arguments as kwargs
"""
function wignerdjmn_ELZOUKA(θ::R; s::I=0, m::I=0, n::I=0) where {R <: Real, I <: Integer}
    return wignerdjmn_ELZOUKA(s,m,n,θ)
end

"""
    Calculate Wigner-d function, using recurrence and memoize
"""
@memoize function wignerdjmn_recurrence_memoize(s::I, m::I, n::I, θ::R) where {R <: Real, I <: Integer}

    s_min = max(abs(m), abs(n))

    if typeof(s) == Int && 2s_min >= 21
        s = big(s); m = big(m); n = big(n)
        s, s_min = promote(s, s_min)
    end

    x = cos(θ)

    if m > n # I need to memoize only m < n, because it is relevant to calculating angular functions
        return (-1)^(m - n) * wignerdjmn_recurrence_memoize(s, n, m, θ) # eq. B.5

    elseif θ < 0 # I need to memoize only θ > 0
        return (-1)^(m - n) * wignerdjmn_recurrence_memoize(s, m, n, -θ)

    elseif n < 0 # I need to memoize only n > 0
        return (-1)^(m - n) * wignerdjmn_recurrence_memoize(s, -m, -n, θ)

    elseif s < s_min
        return zero(θ)

    elseif s == s_min
        if typeof(s) == Int && (factorial((abs(m - n))) * factorial((abs(m + n)))) < 0 # check if the denominator may cause overflow
            s = big(s); m = big(m); n = big(n)
            s, s_min = promote(s, s_min)
        end
        # calculate ξ from eq. B.16
        if n >= m
            ξ_mn = one(s)
        else
            ξ_mn = (-1)^(m - n)
        end

        # calculate d^s_min__m_n(θ) from eq. B.24
        return convert(typeof(θ),
            ξ_mn * 2.0^(-s_min) * sqrt(
                factorial((2s_min)) / (factorial((abs(m - n))) * factorial((abs(m + n))))
            ) *
            (1 - x)^(abs(m - n) / 2) *
            (1 + x)^(abs(m + n) / 2)
        )

    else
        d_s_here_plus_1 = zero(θ)
        for s_here = s_min:s - 1
            d_s_here_plus_1 = 1 / (s_here * sqrt((s_here + 1)^2 - m^2) * sqrt((s_here + 1)^2 - n^2)) * (
                (2s_here + 1) * (s_here * (s_here + 1) * x - m * n) * wignerdjmn_recurrence_memoize(s_here, m, n, θ)
                - 1 * (s_here + 1) * sqrt(s_here^2 - m^2) * sqrt(s_here^2 - n^2) * wignerdjmn_recurrence_memoize(s_here - 1, m, n, θ)
            ) # eq. B.22
        end
        return convert(typeof(θ), d_s_here_plus_1)
    end
end

"""
    Calculate Wigner-d function, using recurrence
"""
function wignerdjmn_recurrence(s::I, m::I, n::I, θ::R) where {R <: Real, I <: Integer}

    s_min = max(abs(m), abs(n))

    if typeof(s) == Int && 2s_min >= 21
        s = big(s); m = big(m); n = big(n)
        s, s_min = promote(s, s_min)
    end

    x = cos(θ)

    if m > n # I need to memoize only m < n, because it is relevant to calculating angular functions
        return (-1)^(m - n) * wignerdjmn_recurrence(s, n, m, θ) # eq. B.5

    elseif θ < 0 # I need to memoize only θ > 0
        return (-1)^(m - n) * wignerdjmn_recurrence(s, m, n, -θ)

    elseif n < 0 # I need to memoize only n > 0
        return (-1)^(m - n) * wignerdjmn_recurrence(s, -m, -n, θ)

    elseif s < s_min
        return zero(θ)

    elseif s == s_min
        if typeof(s) == Int && (factorial((abs(m - n))) * factorial((abs(m + n)))) < 0 # check if the denominator may cause overflow
            s = big(s); m = big(m); n = big(n)
            s, s_min = promote(s, s_min)
        end
        # calculate ξ from eq. B.16
        if n >= m
            ξ_mn = one(s)
        else
            ξ_mn = (-1)^(m - n)
        end

        # calculate d^s_min__m_n(θ) from eq. B.24
        return convert(typeof(θ),
            ξ_mn * 2.0^(-s_min) * sqrt(
                factorial((2s_min)) / (factorial((abs(m - n))) * factorial((abs(m + n))))
            ) *
            (1 - x)^(abs(m - n) / 2) *
            (1 + x)^(abs(m + n) / 2)
        )

    else
        d_s_here_plus_1 = zero(θ)
        s_here = s_min
        d_s_here_minus_1 = wignerdjmn_recurrence(s_here - 1, m, n, θ)
        d_s_here = wignerdjmn_recurrence(s_here, m, n, θ)
        for s_here = s_min:s - 1
            d_s_here_plus_1 = 1 / (s_here * sqrt((s_here + 1)^2 - m^2) * sqrt((s_here + 1)^2 - n^2)) * (
                (2s_here + 1) * (s_here * (s_here + 1) * x - m * n) * d_s_here
                - 1 * (s_here + 1) * sqrt(s_here^2 - m^2) * sqrt(s_here^2 - n^2) * d_s_here_minus_1
            ) # eq. B.22
            d_s_here_minus_1 = d_s_here
            d_s_here = d_s_here_plus_1
        end
        return convert(typeof(θ), d_s_here_plus_1)
    end
end

wignerdjmn = wignerdjmn_ELZOUKA # I did it to make it work with auto-diff, although "wignerdjmn_ELZOUKA" is not efficient.
# wignerdjmn = wignerdjmn_recurrence_memoize
# I may need to define "ChainRulesCore.@scalar_rule" for "WignerD.wignerdjmn"
# wignerdjmn = WignerD.wignerdjmn

# derivative of wigner-D
# TODO: add the special case of θ=0
function ∂wignerdjmn_by_∂θ(s::I, m::I, n::I, θ::R; numerical_derivative=false, verysmallnumber=1e-30) where {R <: Real, I <: Integer}
    """
    derivative of wigner-d with resepect to θ. Adopted from eq. B.25 from Mishchenko, M.I., Travis, L.D., and Lacis, A.A. (2002). Scattering, absorption, and emission of light by small particles (Cambridge University Press).
    """
    if numerical_derivative
        return (wignerdjmn(s, m, n, θ + verysmallnumber) - wignerdjmn(s, m, n, θ - verysmallnumber)) / (verysmallnumber * 2)
    else
        #try
        # I don't know why this is not giving the same results as the second one.
        #    return      (m - n * cos(θ)) / sin(θ) * wignerdjmn(s, m, n, θ) + sqrt((s + n) * (s - n + 1)) * wignerdjmn(s, m, n - 1, θ) # first line of equation B.25
        #catch
        #println("I am ∂wignerdjmn_by_∂θ, s=$s; m=$m; n=$n")
        if s >= n
            return -1 * (Zygote.dropgrad(m - n) * cos(θ)) / sin(θ) * wignerdjmn(s, m, n, θ) - Zygote.dropgrad(sqrt((s - n) * (s + n + 1))) * wignerdjmn(s, m, n + 1, θ) # second line of equation B.25
        else # the second part is zero. TODO: Is there a more elegent way to do it? using short circuit for example?
            return -1 * (Zygote.dropgrad(m - n) * cos(θ)) / sin(θ) * wignerdjmn(s, m, n, θ) # second line of equation B.25
        end
        #end
    end
end

# adding this rule expedites Zygote

ChainRulesCore.@scalar_rule(
    wignerdjmn(s::Integer, m::Integer, n::Integer, θ::Real),
    (
        ChainRulesCore.ZeroTangent(),
        ChainRulesCore.ZeroTangent(),
        ChainRulesCore.ZeroTangent(),
        ∂wignerdjmn_by_∂θ(s, m, n, θ)
    )
)

#ChainRulesCore.@scalar_rule(wignerdjmn(θ::Real; s::Integer=0, m::Integer=0, n::Integer=0), (ZeroTangent(), ZeroTangent(), ZeroTangent(), ∂wignerdjmn_by_∂θ(s, m, n, θ)))

#############################################################################################
# calculate π(θ) and τ(θ)
function πₘₙ(m::I, n::I, θ::R) where {R <: Real, I <: Integer}
    return Zygote.dropgrad(m) / sin(θ) * wignerdjmn(n, 0, m, θ)
end

function τₘₙ(m::I, n::I, θ::R) where {R <: Real, I <: Integer}
    return ∂wignerdjmn_by_∂θ(n, 0, m, θ)
end

# TODO: we may create a lookup table
function γ_mn(m::I, n::I) where {I <: Integer}
    if !((n + abs(m)) <= 19 || typeof(m) == BigInt )
        # this is to solve overflow problems
        # I hope this doesn't break Zygote autodiff
        m = big(m); n = big(n)
    end

    return Zygote.dropgrad(
            sqrt(
            ((2n + 1) * factorial(n - m)) / (4π * n * (n + 1) * factorial(n + m))
        )
    ) # equation C.22
end

# TODO: we may create a lookup table
function γ_mn_dash(m::I, n::I) where {I <: Integer}
    if !((n + abs(m)) <= 19 || typeof(m) == BigInt )
        # this is to solve overflow problems
        # I hope this doesn't break Zygote autodiff
        m = big(m); n = big(n)
    end

    return Zygote.dropgrad(sqrt(
            ((2n + 1) * factorial(n - m)) / (4π * factorial(n + m))
        )
    ) # equation C.25
end

"""return the radial function that is appropriate with the type of VSWF
"""
function get_radial_function_and_special_derivative_given_kind(kind)
    if kind in ["regular", "incoming", 1]
        radial_function = spherical_Bessel_j_n
        radial_function_special_derivative = one_over_x_by_∂_x_j_n_by_∂x
    elseif kind in ["irregular", "outgoing", 2]
        radial_function = spherical_Hankel_h1_n
        radial_function_special_derivative = one_over_x_by_∂_x_h_n_by_∂x
    else
        throw(DomainError(""" 'kind' has to be one of the following: ["regular", "incoming", "irregular", "outgoing"] """))
    end

    return radial_function, radial_function_special_derivative
end

function sqrt_factorial_n_plus_m_over_factorial_n_minus_m(m::I, n::I) where {I <: Integer}
    if n + abs(m) > 20 # Is there a better way to do it? I did it because factorial fails for input of type Int that is larger than 20
        m = big(m); n = big(n)
    end
    return Zygote.dropgrad(
        sqrt(factorial(n + m) / factorial(n - m))
    )
end


text:
@testset "Simulation Build Tests" begin
    models = create_simulation_build_test_problems(get_template_basic_uc_simulation())
    sequence = SimulationSequence(
        models=models,
        feedforwards=Dict(
            "ED" => [
                SemiContinuousFeedforward(
                    component_type=ThermalStandard,
                    source=OnVariable,
                    affected_values=[ActivePowerVariable],
                ),
            ],
        ),
        ini_cond_chronology=InterProblemChronology(),
    )
    sim = Simulation(
        name="test",
        steps=1,
        models=models,
        sequence=sequence,
        simulation_folder=mktempdir(cleanup=true),
    )

    build_out = build!(sim)
    @test build_out == PSI.BuildStatus.BUILT

    for field in fieldnames(SimulationSequence)
        if fieldtype(SimulationSequence, field) == Union{Dates.DateTime, Nothing}
            @test getfield(sim.sequence, field) !== nothing
        end
    end

    @test length(findall(x -> x == 2, sequence.execution_order)) == 24
    @test length(findall(x -> x == 1, sequence.execution_order)) == 1

    state = PSI.get_simulation_state(sim)

    uc_vars = [OnVariable, StartVariable, StopVariable]
    ed_vars = [ActivePowerVariable]
    for (key, data) in state.decision_states.variables
        if PSI.get_entry_type(key) ∈ uc_vars
            count, _ = size(data.values)
            @test count == 24
        elseif PSI.get_entry_type(key) ∈ ed_vars
            count, _ = size(data.values)
            @test count == 288
        end
    end
end

@testset "Simulation with provided initial time" begin
    models = create_simulation_build_test_problems(get_template_basic_uc_simulation())
    sequence = SimulationSequence(
        models=models,
        feedforwards=Dict(
            "ED" => [
                SemiContinuousFeedforward(
                    component_type=ThermalStandard,
                    source=OnVariable,
                    affected_values=[ActivePowerVariable],
                ),
            ],
        ),
        ini_cond_chronology=InterProblemChronology(),
    )
    second_day = DateTime("1/1/2024  23:00:00", "d/m/y  H:M:S") + Hour(1)
    sim = Simulation(
        name="test",
        steps=1,
        models=models,
        sequence=sequence,
        simulation_folder=mktempdir(cleanup=true),
        initial_time=second_day,
    )
    build_out = build!(sim)
    @test build_out == PSI.BuildStatus.BUILT

    for model in PSI.get_decision_models(PSI.get_models(sim))
        @test PSI.get_initial_time(model) == second_day
    end

    for field in fieldnames(SimulationSequence)
        if fieldtype(SimulationSequence, field) == Union{Dates.DateTime, Nothing}
            @test getfield(sim.sequence, field) !== nothing
        end
    end

    @test length(findall(x -> x == 2, sequence.execution_order)) == 24
    @test length(findall(x -> x == 1, sequence.execution_order)) == 1
end

@testset "Negative Tests (Bad Parametrization)" begin
    models = create_simulation_build_test_problems(get_template_basic_uc_simulation())
    sequence = SimulationSequence(
        models=models,
        feedforwards=Dict(
            "ED" => [
                SemiContinuousFeedforward(
                    component_type=ThermalStandard,
                    source=OnVariable,
                    affected_values=[ActivePowerVariable],
                ),
            ],
        ),
        ini_cond_chronology=InterProblemChronology(),
    )

    @test_throws UndefKeywordError sim = Simulation(name="test", steps=1)

    sim = Simulation(
        name="test",
        steps=1,
        models=models,
        sequence=sequence,
        simulation_folder=mktempdir(cleanup=true),
        initial_time=Dates.now(),
    )

    @test_throws IS.ConflictingInputsError build!(sim, console_level=Logging.AboveMaxLevel)

    sim = Simulation(
        name="fake_path",
        steps=1,
        models=models,
        sequence=sequence,
        simulation_folder="fake_path",
    )

    @test_throws IS.ConflictingInputsError PSI._check_folder(sim)
end

@testset "Multi-Stage Hydro Simulation Build" begin
    sys_md = PSB.build_system(SIIPExampleSystems, "5_bus_hydro_wk_sys")

    sys_uc = PSB.build_system(SIIPExampleSystems, "5_bus_hydro_uc_sys")
    transform_single_time_series!(sys_uc, 48, Hour(24))

    sys_ed = PSB.build_system(SIIPExampleSystems, "5_bus_hydro_ed_sys")

    template = ProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, ThermalStandard, ThermalBasicUnitCommitment)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    set_device_model!(template, HydroEnergyReservoir, HydroDispatchReservoirBudget)

    template_uc = ProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template_uc, ThermalStandard, ThermalBasicUnitCommitment)
    set_device_model!(template_uc, PowerLoad, StaticPowerLoad)
    set_device_model!(template_uc, HydroEnergyReservoir, HydroDispatchRunOfRiver)

    template_ed = ProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template_ed, ThermalStandard, ThermalBasicUnitCommitment)
    set_device_model!(template_ed, PowerLoad, StaticPowerLoad)
    set_device_model!(template_ed, HydroEnergyReservoir, HydroDispatchRunOfRiver)

    models = SimulationModels([
        DecisionModel(
            template,
            name="MD",
            sys_md,
            initialize_model=false,
            system_to_file=false,
            optimizer=HiGHS_optimizer,
        ),
        DecisionModel(
            template_uc,
            name="UC",
            sys_uc,
            initialize_model=false,
            system_to_file=false,
            optimizer=HiGHS_optimizer,
        ),
        DecisionModel(
            template_ed,
            name="ED",
            sys_ed,
            initialize_model=false,
            system_to_file=false,
            optimizer=HiGHS_optimizer,
        ),
    ])

    feedforwards = Dict(
        "UC" => [
            EnergyLimitFeedforward(
                source=ActivePowerVariable,
                affected_values=[ActivePowerVariable],
                component_type=HydroEnergyReservoir,
                number_of_periods=24,
            ),
        ],
        "ED" => [
            EnergyLimitFeedforward(
                source=ActivePowerVariable,
                affected_values=[ActivePowerVariable],
                component_type=HydroEnergyReservoir,
                number_of_periods=12,
            ),
        ],
    )

    test_sequence = SimulationSequence(
        models=models,
        ini_cond_chronology=InterProblemChronology(),
        feedforwards=feedforwards,
    )

    sim = Simulation(
        name="test_md",
        steps=2,
        models=models,
        sequence=test_sequence,
        simulation_folder=mktempdir(cleanup=true),
    )
    @test build!(sim, serialize=false) == PSI.BuildStatus.BUILT
end

@testset "Test SemiContinuous Feedforward with Active and Reactive Power variables" begin
    template_uc = get_template_basic_uc_simulation()
    set_network_model!(template_uc, NetworkModel(DCPPowerModel, use_slacks=true))
    # network slacks added because of data issues
    template_ed =
        get_template_nomin_ed_simulation(NetworkModel(ACPPowerModel, use_slacks=true))
    c_sys5_hy_uc = PSB.build_system(PSITestSystems, "c_sys5_hy_uc")
    c_sys5_hy_ed = PSB.build_system(PSITestSystems, "c_sys5_hy_ed")
    models = SimulationModels(
        decision_models=[
            DecisionModel(
                template_uc,
                c_sys5_hy_uc;
                name="UC",
                optimizer=HiGHS_optimizer,
                initialize_model=false,
            ),
            DecisionModel(
                template_ed,
                c_sys5_hy_ed;
                name="ED",
                optimizer=ipopt_optimizer,
                initialize_model=false,
            ),
        ],
    )

    sequence = SimulationSequence(
        models=models,
        feedforwards=Dict(
            "ED" => [
                SemiContinuousFeedforward(
                    component_type=ThermalStandard,
                    source=OnVariable,
                    affected_values=[ActivePowerVariable, ReactivePowerVariable],
                ),
            ],
        ),
        ini_cond_chronology=InterProblemChronology(),
    )

    sim = Simulation(
        name="reactive_feedforward",
        steps=2,
        models=models,
        sequence=sequence,
        simulation_folder=mktempdir(cleanup=true),
    )
    build_out = build!(sim)
    @test build_out == PSI.BuildStatus.BUILT
    ac_power_model = PSI.get_simulation_model(PSI.get_models(sim), :ED)
    c = PSI.get_constraint(
        PSI.get_optimization_container(ac_power_model),
        FeedforwardSemiContinousConstraint(),
        ThermalStandard,
        "ActivePowerVariable_ub",
    )
    @test !isempty(c)
    c = PSI.get_constraint(
        PSI.get_optimization_container(ac_power_model),
        FeedforwardSemiContinousConstraint(),
        ThermalStandard,
        "ActivePowerVariable_lb",
    )
    @test !isempty(c)
    c = PSI.get_constraint(
        PSI.get_optimization_container(ac_power_model),
        FeedforwardSemiContinousConstraint(),
        ThermalStandard,
        "ReactivePowerVariable_ub",
    )
    @test !isempty(c)
    c = PSI.get_constraint(
        PSI.get_optimization_container(ac_power_model),
        FeedforwardSemiContinousConstraint(),
        ThermalStandard,
        "ReactivePowerVariable_lb",
    )
    @test !isempty(c)
end


text:
using JULES, Oceananigans
using Plots
using Printf

Nx = Ny = 1
Nz = 32
L = 10e3

grid = RegularCartesianGrid(size=(Nx, Ny, Nz), halo=(2, 2, 2),
                            x=(0, L), y=(0, L), z=(0, L))
#####
##### Dry adiabatic atmosphere
#####

gas = IdealGas()
Rᵈ, cₚ, cᵥ = gas.Rᵈ, gas.cₚ, gas.cᵥ

g  = 9.80665
pₛ = 100000
θₛ = 300
πₛ = 1

θ(x, y, z) = θₛ
π(x, y, z) = πₛ - g*z / (cₚ*θₛ)
p(x, y, z) = pₛ * π(x, y, z)^(cₚ/Rᵈ)
ρ(x, y, z) = pₛ / (Rᵈ * θ(x, y, z)) * π(x, y, z)^(cᵥ/Rᵈ)
Θ(x, y, z) = ρ(x, y, z) * θ(x, y, z)

#####
##### Create model
#####

model = CompressibleModel(grid=grid, buoyancy=gas, reference_pressure=pₛ,
                          thermodynamic_variable=ModifiedPotentialTemperature(),
                          tracers=(:Θᵐ,))

#####
##### Set initial conditions
#####

set!(model.density, ρ)
set!(model.tracers.Θᵐ, θ)

# Initial profiles
Θ₀_prof = model.tracers.Θᵐ[1, 1, 1:Nz]
ρ₀_prof = model.density[1, 1, 1:Nz]

# Arrays we will use to store a time series of ρw(z = L/2).
times = [model.clock.time]
ρw_ts = [model.momenta.ρw[1, 1, Int(Nz/2)]]

#####
##### Time step and keep plotting vertical profiles of ρθ′, ρw, and ρ′.
#####

# @animate for i=1:100
while model.clock.time < 500
    time_step!(model; Δt=0.5, Nt=10)

    @show model.clock.time
    Θ_prof = model.tracers.Θᵐ[1, 1, 1:Nz] .- Θ₀_prof
    W_prof = model.momenta.ρw[1, 1, 1:Nz+1]
    ρ_prof = model.density[1, 1, 1:Nz] .- ρ₀_prof

    Θ_plot = plot(Θ_prof, grid.zC, xlim=(-2.5, 2.5), xlabel="Theta_prime", ylabel="z (m)", label="")
    W_plot = plot(W_prof, grid.zF, xlim=(-1.2, 1.2), xlabel="rho*w", label="")
    ρ_plot = plot(ρ_prof, grid.zC, xlim=(-0.01, 0.01), xlabel="rho", label="")

    t_str = @sprintf("t = %d s", model.clock.time)
    display(plot(Θ_plot, W_plot, ρ_plot, title=["" "$t_str" ""], layout=(1, 3), show=true))

    push!(times, model.clock.time)
    push!(ρw_ts, model.momenta.ρw[1, 1, Int(Nz/2)])
end

#####
##### Plot the initial profile and the hydrostatically balanced profile.
#####

ρ∞_prof = model.density[1, 1, 1:Nz]

θ₀_prof = Θ₀_prof ./ ρ₀_prof
θ∞_prof = model.tracers.Θᵐ[1, 1, 1:Nz] ./ ρ∞_prof

θ_plot = plot(θ₀_prof, grid.zC, xlabel="theta (K)", ylabel="z (m)", label="initial", legend=:topleft)
plot!(θ_plot, θ∞_prof, grid.zC, label="balanced")

ρ_plot = plot(ρ₀_prof, grid.zC, xlabel="rho (kg/m³)", ylabel="z (m)", label="initial")
plot!(ρ_plot, ρ∞_prof, grid.zC, label="balanced")

display(plot(θ_plot, ρ_plot, show=true))

#####
##### Plot timeseries of maximum ρw
#####

plot(times, ρw_ts, linewidth=2, xlabel="time (s)", ylabel="rho*w(z=$(Int(L/2)) m)", label="", show=true)


text:
struct TreeSearchHelper
    depth::Int
    visited::Bool
end

function do_nothing(node, depth, visited)
end

function show_node(node, depth, visited)
    show(node)
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
import PWS
using PWS: SMCEstimate, name, propagate, DirectMCEstimate, marginal_configuration, ConditionalEnsemble, MarginalEnsemble, energy_difference, generate_configuration, log_marginal, logpdf, simulate
using Test
using Statistics

system = PWS.chemotaxis_system(dtimes=0:0.025:0.5)

@test system.dist.aggregator.update_map == [0, 0, 0, 0, 1, 2]

conf = generate_configuration(system)
mconf = marginal_configuration(conf)
cens = ConditionalEnsemble(system)
mens = MarginalEnsemble(system)

@test energy_difference(conf, cens) == energy_difference(mconf, mens)

algs = [SMCEstimate(16), DirectMCEstimate(16)]

result = Dict(map(algs) do alg
    mi = map(1:20) do i
        cresult = simulate(alg, conf, cens)
        mresult = simulate(alg, mconf, mens)
        log_marginal(cresult) - log_marginal(mresult)
    end
    name(alg) => hcat(mi...)
end)

# Test that all results lie within 3 standard deviations of each other.
# This should find coarse implementation errors for specific estimates.
for (name, mi) in result
    lower_bound = mean(mi, dims=2) - 3*std(mi, dims=2)
    for (name2, mi2) in result
        upper_bound = mean(mi2, dims=2) + 3*std(mi2, dims=2)
        @test all(lower_bound .<= upper_bound)
    end
end

# # Some plotting routines are commented out below to see what goes wrong.
# # Use them to see what's wrong when the tests above fail.
#--
# using Plots
# plot()
# colormap = Dict("SMC" => 1, "Direct MC" => 2)
# for (name, mi) in result
#     plot!(system.dtimes, mean(mi, dims=2), ribbon=3*std(mi, dims=2), color=colormap[name], label=name)
#     plot!(system.dtimes, mi, color=colormap[name], linewidth=0.3, label="")
# end
# plot!()
# --


text:
module Dynamics

# using TrajectoryOptimization: Model, UncertainModel, Trajectory, Constraint, Equality, Inequality, Problem
using TrajectoryOptimization: Equality, Inequality
using TrajectoryOptimization
# using RigidBodyDynamics
using LinearAlgebra
using DocStringExtensions
using StaticArrays
using Parameters
import TrajectoryOptimization: dynamics, AbstractModel, RigidBody, state_diff, state_diff_jacobian,
    state_diff_size, orientation

export
    pendulum,
    doublependulum,
    cartpole,
    cartpole_urdf,
    ballonbeam,
    acrobot_model,
    quadrotor,
    quadrotor_euler,
    kuka,
    doubleintegrator,
    car

include("../dynamics/testrigidbody.jl")
include("../dynamics/pendulum.jl")
# include("../dynamics/doublependulum.jl")
include("../dynamics/acrobot.jl")
# include("../dynamics/ballonbeam.jl")
include("../dynamics/cartpole.jl")
include("../dynamics/quadrotor.jl")
# include("../dynamics/quadrotor_euler.jl")
# include("../dynamics/kuka.jl")
include("../dynamics/double_integrator.jl")
include("../dynamics/car.jl")
include("../dynamics/rigidbody.jl")
include("../dynamics/satellite.jl")
include("../dynamics/yak_plane.jl")
include("../dynamics/noisy_rigidbody.jl")
include("../dynamics/freebody.jl")
include("../dynamics/copymodel.jl")

end


text:
module PortfolioBedtest

using Dates
using CSV
using HTTP

export TimedEvent, glue, yahoo, momentum_score, mutate, Signal, Broker, current, assets, update, issignal, sim, Recorder, AbstractOrder, order, execute
export filter_eom_indices

# From https://github.com/JuliaQuant/MarketData.jl/blob/master/src/downloads.jl
"""
    struct YahooOpt <: AbstractQueryOpt
      period1  # the start time
      period2  # the end time
      interval # "1d", "1wk" or "1mo"
      events   # currently only `:history` supported
    end
The Yahoo Finance HTTP API query object.

Examples
========
```julia
julia> t = Dates.now()
2020-08-09T01:38:04.735
julia> YahooOpt(period1 = t - Year(2), period2 = t)
YahooOpt{DateTime} with 4 entries:
  :period1  => 1533778685
  :period2  => 1596937085
  :interval => "1d"
  :events   => :history
```
"""
struct YahooOpt
    period1::DateTime
    period2::DateTime
    interval::String
    events::String
end
function YahooOpt(; period1::DateTime = DateTime(1971, 2, 8),
                    period2::DateTime = Dates.now(),
                    interval::String  = "1d",
                    events::Symbol    = "history")
    YahooOpt(period1, period2, interval, events)
end

function asquery(yo::YahooOpt)
    return ("period1" => round(Int, datetime2unix(yo.period1)),
            "period2" => round(Int, datetime2unix(yo.period2)),
            "interval" => yo.interval,
            "events" => yo.events)
end

function yahoo(sym; period1 = DateTime(1971, 2, 8), period2 = Dates.now(), interval = "1d", events = "history")
    yo = YahooOpt(period1, period2, interval, events)
    yahoo(sym, yo)
end

function yahoo(sym, yo::YahooOpt)
    host = rand(("query1", "query2"))
    url  = "https://$host.finance.yahoo.com/v7/finance/download/$sym"
    res  = HTTP.get(url, query = asquery(yo))
    res.status == 200 || throw("Error requesting yahoo $res")
    csv = CSV.File(res.body, missingstrings = ["null"])
end

########################################
# TimedEvents (idea from Timestamps.jl)
########################################
struct TimedEvent{D, T}
    ts::D
    event::T
end

function Base.show(io::IO, te::TimedEvent)
    print(io, te.ts, ": ", te.event)
end

Base.isless(v1::TimedEvent{D}, v2::TimedEvent{D}) where D = v1.ts < v2.ts

Base.:(-)(v1::TimedEvent{D}, v2::TimedEvent{D}) where D = TimedEvent(v1.ts, v1.event .- v2.event)
Base.:(+)(v1::TimedEvent{D}, v2::TimedEvent{D}) where D = TimedEvent(v1.ts, v1.event .+ v2.event)
Base.:(/)(v1::TimedEvent{D}, v2::TimedEvent{D}) where D = TimedEvent(v1.ts, v1.event ./ v2.event)
Base.:(/)(v1::TimedEvent, v2::Number) = TimedEvent(v1.ts, v1.event ./ v2)
Base.:(-)(v1::TimedEvent, v2::Number) = TimedEvent(v1.ts, v1.event .- v2)
Base.:(*)(v1::TimedEvent, v2::Number) = TimedEvent(v1.ts, v1.event .* v2)
Base.:(*)(v2::Number, v1::TimedEvent) = TimedEvent(v1.ts, v1.event .* v2)

function glue(f::Function, v1::Vector{TimedEvent{D, T1}}, v2::Vector{TimedEvent{D, T2}}) where {D, T1, T2}
    sort!(v1)
    sort!(v2)

    T = typeof(f(v1[1].event, v2[1].event, 0)) # ugly hack...
    res = TimedEvent{D, T}[]

    # double index
    i1, i2 = 1, 1
    # We should do beggars push here for performance
    @inbounds while (i1 <= length(v1)) & (i2 <= length(v2))
        if v1[i1] < v2[i2]
            while (v1[i1] < v2[i2]) & (i1 <= length(v1))
                push!(res, TimedEvent(v1[i1].ts, f(v1[i1].event, v2[i2].event, 1)))
                i1 += 1
            end
        elseif v2[i2] < v1[i1]
            while (v2[i2] < v1[i1]) & (i2 <= length(v2))
                push!(res, TimedEvent(v2[i2].ts, f(v1[i1].event, v2[i2].event, 2)))
                i2 += 1
            end
        end
        v1[i1].ts != v2[i2].ts && continue
        # there are no duplicates, so we can simplify calculations
        push!(res, TimedEvent(v1[i1].ts, f(v1[i1].event, v2[i2].event, 0)))
        i1 += 1
        i2 += 1
    end

    return res
end

function mfun(v1, v2, k)
    k == 0 && return merge(v1, v2)
    k == 1 && return merge(v1, map(x -> NaN, v2))
    return merge(map(x -> NaN, v1), v2)
end

glue(v1::Vector{TimedEvent{D, T1}}, v2::Vector{TimedEvent{D, T2}}) where {D, T1, T2} = glue(mfun, v1, v2)
glue(v1::Vector, v2::Vector...) = foldl((x, y) -> glue(x, y), v2; init = v1)
glue(f::Function, v1::Vector, v2::Vector...) = foldl((x, y) -> glue(f, x, y), v2; init = v1)

function filter_eom(v)
    res = similar(v, 0)
    for (i, c) in pairs(v)
        i == length(v) && continue
        month(v[i].ts) != month(v[i+1].ts) && push!(res, c)
    end
    return res
end

function filter_eom_indices(v)
    res = Int[]
    for (i, c) in pairs(v)
        i == length(v) && continue
        month(v[i].ts) != month(v[i+1].ts) && push!(res, i)
    end
    return res
end

Base.map(f::Function, te::TimedEvent) = TimedEvent(te.ts, map(f, te.event))
mutate(f::Function, te::TimedEvent) = TimedEvent(te.ts, f(te.event))

function momentum_score(v)
    v = filter_eom(v)
    res = similar(v)
    for i in eachindex(v)
        if i < 13
            res[i] = map(x -> NaN, v[i])
            continue
        end
        p0 = values(v[i].event)
        p1 = values(v[i - 1].event)
        p3 = values(v[i - 3].event)
        p6 = values(v[i - 6].event)
        p12 = values(v[i - 12].event)
        momentum = @. 12 * (p0 / p1 - 1) + 4 * (p0 / p3 - 1) + 2 * (p0 / p6 - 1) + (p0 / p12 - 1)
        res[i] = TimedEvent(v[i].ts, NamedTuple{keys(v[i].event)}(momentum))
    end
    res
end

########################################
# Simulation part
########################################

abstract type AbstractOrder end
function order end

struct Signal{D, T}
    idx::Int
    signal::Vector{TimedEvent{D, T}}
end
Signal(sig) = Signal(1, sig)

function update(sig::Signal, ts)
    sig.signal[sig.idx].ts >= ts && return sig
    sig.idx == length(sig.signal) && return sig
    idx = sig.idx
    while idx <= length(sig.signal)
        sig.signal[idx].ts >= ts && return Signal(idx, sig.signal)
        idx += 1
    end
    return Signal(length(signal), sig.signal)
end
issignal(sig::Signal, ts) = @inbounds sig.signal[sig.idx].ts == ts
current(sig::Signal) = sig.signal[sig.idx]

struct Recorder{D, T}
    data::Vector{TimedEvent{D, T}}
end

function record(r::Recorder, broker, prices, ts)
    event = (; shares = broker.shares, total = assets(broker, prices))
    te = TimedEvent(ts, event)
    push!(r.data, te)
    nothing
end

function Base.empty!(r::Recorder)
    empty!(r.data)
end

abstract type AbstractBroker end

struct Broker{T} <: AbstractBroker
    cash::Float64
    shares::T
end

assets(broker::Broker, prices) = sum(values(broker.shares) .* values(prices)) + broker.cash

function execute end

function sim(broker::AbstractBroker, strategy, prices, signal, recorder; recrestart = true)
    if recrestart
        empty!(recorder)
    end

    for i in eachindex(prices)
        ev = prices[i]
        ts = ev.ts
        signal = update(signal, ts)
        if issignal(signal, ts)
            ord = order(strategy, broker, current(signal).event, ev.event)
            broker = execute(broker, ord, ev.event)
        end
        record(recorder, broker, ev.event, ts)
    end

    return assets(broker, prices[end].event)
end

end # module


text:
using Distributions, BayesianNonparametricStatistics, PyPlot

sde = SDEWithConstantVariance(x->sinpi(2*x), 1.0, 0.0, 10000.0, 0.01)
X = rand(sde)
distribution = MvNormal([k^(-1.0) for k in 1:10])
Π = GaussianProcess([fourier(k) for k in 1:10], distribution)
M = SDEModel(1.0, 0.0)
postΠ = calculateposterior(Π, X, M)
clf()
x=0.0:0.01:1.0
for k in 1:10
    f = rand(postΠ)
    plot(x, f.(x))
end