text:
"""
`init_subspec!(m::Model1010_1spread)`

Initializes a model subspecification by overwriting parameters from
the original model object with new parameter objects. This function is
called from within the model constructor.
"""
function init_subspec!(m::Model1010_1spread)
    if subspec(m) == "ss1"
        return
    elseif subspec(m) == "ss2"
        return ss2!(m)
    elseif subspec(m) == "ss3"
        return ss3!(m)
    elseif subspec(m) == "ss4"
        return ss4!(m)
    elseif subspec(m) == "ss5"
        return ss5!(m)
    elseif subspec(m) == "ss6"
        return ss6!(m)
    elseif subspec(m) == "ss7"
        return ss7!(m)
    elseif subspec(m) == "ss8"
        return ss8!(m)
    elseif subspec(m) == "ss9"
        return ss9!(m)
    elseif subspec(m) == "ss10"
        return ss10!(m)
    elseif subspec(m) == "ss11"
        return ss11!(m)
    elseif subspec(m) == "ss12"
        return ss12!(m)
    elseif subspec(m) == "ss13"
        return ss13!(m)
    elseif subspec(m) == "ss14"
        return ss14!(m)
    elseif subspec(m) == "ss15"
        return ss15!(m)
    elseif subspec(m) == "ss16"
        return ss16!(m)
    elseif subspec(m) == "ss17"
        return ss17!(m)
    elseif subspec(m) == "ss18"
        return ss18!(m)
    elseif subspec(m) == "ss19"
        return ss19!(m)
    elseif subspec(m) == "ss20"
        return ss20!(m)
    elseif subspec(m) == "ss21"
        return ss21!(m)
    elseif subspec(m) == "ss22"
        return ss22!(m)
    elseif subspec(m) == "ss23"
        return ss23!(m)
    else
        error("This subspec is not defined.")
    end
end

function ss4!(m::Model1010_1spread)
    # ss2, with AR(1) process for AAA spread measurement error

    m <= parameter(:ρ_AAA, 0.5, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1),
                   fixed=false,
                   description="ρ_AAA: AR(1) coefficient in the AAA spread process.",
                   tex_label="\\rho_{AAA}")

    m <= parameter(:lnb_liq, 0.47, (1e-5, 10.),   (1e-5, 10.), ModelConstructors.Exponential(), GammaAlt(0.47, 0.1),
                   fixed=false, scaling = x -> (1 + x/100)^0.25,
                   description="ln(b_liq_*): Liquidity premium (percent annualized).",
                   tex_label="ln(b_{liq})")

    m <= parameter(:lnb_safe, 0.26,  (1e-5, 10.),   (1e-5, 10.), ModelConstructors.Exponential(), GammaAlt(0.26, 0.1),
                   fixed=false, scaling = x -> (1 + x/100)^0.25,
                   description="ln(b_safe_*): Safety premium (percent annualized).",
                   tex_label="ln(b_{safe})")
end

function ss5!(m::Model1010_1spread)

    # ss4, with AR(1) process for BBB spread measurement error in addition to AAA

    m <= parameter(:ρ_AAA, 0.5, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1),
                   fixed=false,
                   description="ρ_AAA: AR(1) coefficient in the AAA spread process.",
                   tex_label="\\rho_{AAA}")

    m <= parameter(:ρ_BBB, 0.5, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1),
                   fixed=false,
                   description="ρ_BBB: AR(1) coefficient in the BBB spread process.",
                   tex_label="\\rho_{BBB}")

    m <= parameter(:σ_BBB, 0.1, (1e-8, 5.),(1e-8, 5.),ModelConstructors.Exponential(),RootInverseGamma(2., 0.10),
                   fixed=false,
                   description="σ_BBB: Standard deviation on the AR(1) process for measurement error on the BBB spread.",
                   tex_label="\\sigma_{BBB}")

    m <= parameter(:lnb_liq, 0.47, (1e-5, 10.),   (1e-5, 10.), ModelConstructors.Exponential(), GammaAlt(0.47, 0.1),
                   fixed=false, scaling = x -> (1 + x/100)^0.25,
                   description="ln(b_liq_*): Liquidity premium (percent annualized).",
                   tex_label="ln(b_{liq})")

    m <= parameter(:lnb_safe, 0.26,  (1e-5, 10.),   (1e-5, 10.), ModelConstructors.Exponential(), GammaAlt(0.26, 0.1),
                   fixed=false, scaling = x -> (1 + x/100)^0.25,
                   description="ln(b_safe_*): Safety premium (percent annualized).",
                   tex_label="ln(b_{safe})")
end

function ss6!(m::Model1010_1spread)

    # ss1, with iid measurement error on BBB spread

    m <= parameter(:σ_BBB, 0.1, (1e-8, 5.),(1e-8, 5.),ModelConstructors.Exponential(),RootInverseGamma(2., 0.10),
                   fixed=false,
                   description="σ_BBB: Standard deviation on the AR(1) process for measurement error on the BBB spread.",
                   tex_label="\\sigma_{BBB}")
end

function ss7!(m::Model1010_1spread)
    # ss3, with tight prior on σ_b_liqp and σ_b_safep

    m <= parameter(:lnb_liq, 0.47, (1e-5, 10.),   (1e-5, 10.), ModelConstructors.Exponential(), GammaAlt(0.47, 0.1),
                   fixed=false, scaling = x -> (1 + x/100)^0.25,
                   description="ln(b_liq_*): Liquidity premium (percent annualized).",
                   tex_label="ln(b_{liq})")

    m <= parameter(:lnb_safe, 0.26,  (1e-5, 10.),   (1e-5, 10.), ModelConstructors.Exponential(), GammaAlt(0.26, 0.1),
                   fixed=false, scaling = x -> (1 + x/100)^0.25,
                   description="ln(b_safe_*): Safety premium (percent annualized).",
                   tex_label="ln(b_{safe})")

    m <= parameter(:σ_BBB, 0.1, (1e-8, 5.),(1e-8, 5.),ModelConstructors.Exponential(),RootInverseGamma(2., 0.10),
                   fixed=false,
                   description="σ_BBB: Standard deviation on the AR(1) process for measurement error on the BBB spread.",
                   tex_label="\\sigma_{BBB}")

    m <= parameter(:σ_b_liqp, 0.0269, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(100., sqrt(1/400)),
                   fixed=false,
                   description="σ_b_liqp: Standard deviation of stationary component of liquid asset preference shifter process.",
                   tex_label="\\sigma_{b^p, liq}")

    m <= parameter(:σ_b_safep, 0.0269, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(100., sqrt(1/400)),
                   fixed=false,
                   description="σ_b_safep: Standard deviation of stationary component of safe asset preference shifter process.",
                   tex_label="\\sigma_{b^p, safe}")
end

function ss8!(m::Model1010_1spread)

    # Leave lnb_liq, lnb_safe fixed. AR(1) measurement error on both spreads (like ss5).
    m <= parameter(:ρ_AAA, 0.5, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1),
                   fixed=false,
                   description="ρ_AAA: AR(1) coefficient in the AAA spread process.",
                   tex_label="\\rho_{AAA}")

    m <= parameter(:ρ_BBB, 0.5, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1),
                   fixed=false,
                   description="ρ_BBB: AR(1) coefficient in the BBB spread process.",
                   tex_label="\\rho_{BBB}")

    m <= parameter(:σ_BBB, 0.1, (1e-8, 5.),(1e-8, 5.),ModelConstructors.Exponential(),RootInverseGamma(2., 0.10),
                   fixed=false,
                   description="σ_BBB: Standard deviation on the AR(1) process for measurement error on the BBB spread.",
                   tex_label="\\sigma_{BBB}")
end

function ss9!(m::Model1010_1spread)

    # Leave lnb_liq, lnb_safe fixed. AR(1) measurement error on both spreads (like ss5),
    # with prior variance centered at 0.25.

    m <= parameter(:ρ_AAA, 0.5, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1),
                   fixed=false,
                   description="ρ_AAA: AR(1) coefficient in the AAA spread process.",
                   tex_label="\\rho_{AAA}")

    m <= parameter(:σ_AAA, 0.25, (1e-8, 5.),(1e-8, 5.),ModelConstructors.Exponential(),RootInverseGamma(2., 0.25),
                   fixed=false,
                   description="σ_AAA: Standard deviation on the AR(1) process for measurement error on the AAA spread.",
                   tex_label="\\sigma_{AAA}")

    m <= parameter(:ρ_BBB, 0.5, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1),
                   fixed=false,
                   description="ρ_BBB: AR(1) coefficient in the BBB spread process.",
                   tex_label="\\rho_{BBB}")

    m <= parameter(:σ_BBB, 0.25, (1e-8, 5.),(1e-8, 5.),ModelConstructors.Exponential(),RootInverseGamma(2., 0.25),
                   fixed=false,
                   description="σ_BBB: Standard deviation on the AR(1) process for measurement error on the BBB spread.",
                   tex_label="\\sigma_{BBB}")
end

function ss10!(m::Model1010_1spread)

    # Leave lnb_liq, lnb_safe fixed. AR(1) measurement error on both spreads (like ss5),
    # with prior variance centered at 0.25.
    # basically ss9 mixed with ss7

    m <= parameter(:ρ_AAA, 0.5, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1),
                   fixed=false,
                   description="ρ_AAA: AR(1) coefficient in the AAA spread process.",
                   tex_label="\\rho_{AAA}")

    m <= parameter(:σ_AAA, 0.25, (1e-8, 5.),(1e-8, 5.),ModelConstructors.Exponential(),RootInverseGamma(2., 0.25),
                   fixed=false,
                   description="σ_AAA: Standard deviation on the AR(1) process for measurement error on the AAA spread.",
                   tex_label="\\sigma_{AAA}")

    m <= parameter(:ρ_BBB, 0.5, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1),
                   fixed=false,
                   description="ρ_BBB: AR(1) coefficient in the BBB spread process.",
                   tex_label="\\rho_{BBB}")

    m <= parameter(:σ_BBB, 0.25, (1e-8, 5.),(1e-8, 5.),ModelConstructors.Exponential(),RootInverseGamma(2., 0.25),
                   fixed=false,
                   description="σ_BBB: Standard deviation on the AR(1) process for measurement error on the BBB spread.",
                   tex_label="\\sigma_{BBB}")

    m <= parameter(:σ_b_liqp, 0.0269, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2., sqrt(1/400)),
                   fixed=false,
                   description="σ_b_liqp: Standard deviation of stationary component of liquid asset preference shifter process.",
                   tex_label="\\sigma_{b^p, liq}")

    m <= parameter(:σ_b_safep, 0.0269, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2., sqrt(1/400)),
                   fixed=false,
                   description="σ_b_safep: Standard deviation of stationary component of safe asset preference shifter process.",
                   tex_label="\\sigma_{b^p, safe}")
end

function ss11!(m::Model1010_1spread)

    # Leave lnb_liq, lnb_safe fixed. AR(1) measurement error on both spreads (like ss5),
    # with prior variance centered at 0.15.
    # basically ss10 adjusted for quarterly rather than annualized

    m <= parameter(:ρ_AAA, 0.5, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1),
                   fixed=false,
                   description="ρ_AAA: AR(1) coefficient in the AAA spread process.",
                   tex_label="\\rho_{AAA}")

    m <= parameter(:σ_AAA, 0.15, (1e-8, 5.),(1e-8, 5.),ModelConstructors.Exponential(),RootInverseGamma(2., 0.25),
                   fixed=false,
                   description="σ_AAA: Standard deviation on the AR(1) process for measurement error on the AAA spread.",
                   tex_label="\\sigma_{AAA}")

    m <= parameter(:ρ_BBB, 0.5, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1),
                   fixed=false,
                   description="ρ_BBB: AR(1) coefficient in the BBB spread process.",
                   tex_label="\\rho_{BBB}")

    m <= parameter(:σ_BBB, 0.15, (1e-8, 5.),(1e-8, 5.),ModelConstructors.Exponential(),RootInverseGamma(2., 0.25),
                   fixed=false,
                   description="σ_BBB: Standard deviation on the AR(1) process for measurement error on the BBB spread.",
                   tex_label="\\sigma_{BBB}")

    m <= parameter(:σ_b_liqp, 0.0269, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2., sqrt(1/400)/4),
                   fixed=false,
                   description="σ_b_liqp: Standard deviation of stationary component of liquid asset preference shifter process.",
                   tex_label="\\sigma_{b^p, liq}")

    m <= parameter(:σ_b_safep, 0.0269, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2., sqrt(1/400)/4),
                   fixed=false,
                   description="σ_b_safep: Standard deviation of stationary component of safe asset preference shifter process.",
                   tex_label="\\sigma_{b^p, safe}")
end

function ss12!(m::Model1010_1spread)
    # ss7 with lnb_safe and lnb_liq fixed and std dev of permanent processes centered at number adjusted for quarterly

    m <= parameter(:σ_BBB, 0.1, (1e-8, 5.),(1e-8, 5.),ModelConstructors.Exponential(),RootInverseGamma(2., 0.10),
                   fixed=false,
                   description="σ_BBB: Standard deviation on the AR(1) process for measurement error on the BBB spread.",
                   tex_label="\\sigma_{BBB}")

    m <= parameter(:σ_b_liqp, sqrt(1/400)/4, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(100., sqrt(1/400)/4),
                   fixed=false,
                   description="σ_b_liqp: Standard deviation of stationary component of liquid asset preference shifter process.",
                   tex_label="\\sigma_{b^p, liq}")

    m <= parameter(:σ_b_safep, sqrt(1/400)/4, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(100., sqrt(1/400)/4),
                   fixed=false,
                   description="σ_b_safep: Standard deviation of stationary component of safe asset preference shifter process.",
                   tex_label="\\sigma_{b^p, safe}")
end

function ss13!(m::Model1010_1spread)
    # ss12 with AR(1) on measurement error

    m <= parameter(:ρ_AAA, 0.5, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1),
                   fixed=false,
                   description="ρ_AAA: AR(1) coefficient in the AAA spread process.",
                   tex_label="\\rho_{AAA}")

    m <= parameter(:ρ_BBB, 0.5, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1),
                   fixed=false,
                   description="ρ_BBB: AR(1) coefficient in the BBB spread process.",
                   tex_label="\\rho_{BBB}")

    m <= parameter(:σ_BBB, 0.1, (1e-8, 5.),(1e-8, 5.),ModelConstructors.Exponential(),RootInverseGamma(2., 0.10),
                   fixed=false,
                   description="σ_BBB: Standard deviation on the AR(1) process for measurement error on the BBB spread.",
                   tex_label="\\sigma_{BBB}")

    m <= parameter(:σ_b_liqp, sqrt(1/400)/4, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(100., sqrt(1/400)/4),
                   fixed=false,
                   description="σ_b_liqp: Standard deviation of stationary component of liquid asset preference shifter process.",
                   tex_label="\\sigma_{b^p, liq}")

    m <= parameter(:σ_b_safep, sqrt(1/400)/4, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(100., sqrt(1/400)/4),
                   fixed=false,
                   description="σ_b_safep: Standard deviation of stationary component of safe asset preference shifter process.",
                   tex_label="\\sigma_{b^p, safe}")
end

function ss14!(m::Model1010_1spread)
    # ss13, but shutting down permanent components of b shocks

    m <= parameter(:ρ_AAA, 0.5, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1),
                   fixed=false,
                   description="ρ_AAA: AR(1) coefficient in the AAA spread process.",
                   tex_label="\\rho_{AAA}")

    m <= parameter(:ρ_BBB, 0.5, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1),
                   fixed=false,
                   description="ρ_BBB: AR(1) coefficient in the BBB spread process.",
                   tex_label="\\rho_{BBB}")

    m <= parameter(:σ_BBB, 0.1, (1e-8, 5.),(1e-8, 5.),ModelConstructors.Exponential(),RootInverseGamma(2., 0.10),
                   fixed=false,
                   description="σ_BBB: Standard deviation on the AR(1) process for measurement error on the BBB spread.",
                   tex_label="\\sigma_{BBB}")

    m <= parameter(:σ_b_liqp, 0., (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(100., sqrt(1/400)/4),
                   fixed=true,
                   description="σ_b_liqp: Standard deviation of stationary component of liquid asset preference shifter process.",
                   tex_label="\\sigma_{b^p, liq}")

    m <= parameter(:σ_b_safep, 0., (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(100., sqrt(1/400)/4),
                   fixed=true,
                   description="σ_b_safep: Standard deviation of stationary component of safe asset preference shifter process.",
                   tex_label="\\sigma_{b^p, safe}")
end

function ss15!(m::Model1010_1spread)
    # ss5, with fixed lnb_safe and lnb_liq

    m <= parameter(:ρ_AAA, 0.5, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1),
                   fixed=false,
                   description="ρ_AAA: AR(1) coefficient in the AAA spread process.",
                   tex_label="\\rho_{AAA}")

    m <= parameter(:ρ_BBB, 0.5, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1),
                   fixed=false,
                   description="ρ_BBB: AR(1) coefficient in the BBB spread process.",
                   tex_label="\\rho_{BBB}")

    m <= parameter(:σ_BBB, 0.1, (1e-8, 5.),(1e-8, 5.),ModelConstructors.Exponential(),RootInverseGamma(2., 0.10),
                   fixed=false,
                   description="σ_BBB: Standard deviation on the AR(1) process for measurement error on the BBB spread.",
                   tex_label="\\sigma_{BBB}")
end

function ss16!(m::Model1010_1spread)
    # ss13, but not adjusting for maturities in the measurement equation

    m <= parameter(:ρ_AAA, 0.5, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1),
                   fixed=false,
                   description="ρ_AAA: AR(1) coefficient in the AAA spread process.",
                   tex_label="\\rho_{AAA}")

    m <= parameter(:ρ_BBB, 0.5, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1),
                   fixed=false,
                   description="ρ_BBB: AR(1) coefficient in the BBB spread process.",
                   tex_label="\\rho_{BBB}")

    m <= parameter(:σ_BBB, 0.1, (1e-8, 5.),(1e-8, 5.),ModelConstructors.Exponential(),RootInverseGamma(2., 0.10),
                   fixed=false,
                   description="σ_BBB: Standard deviation on the AR(1) process for measurement error on the BBB spread.",
                   tex_label="\\sigma_{BBB}")

    m <= parameter(:σ_b_liqp, sqrt(1/400)/4, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(100., sqrt(1/400)/4),
                   fixed=false,
                   description="σ_b_liqp: Standard deviation of stationary component of liquid asset preference shifter process.",
                   tex_label="\\sigma_{b^p, liq}")

    m <= parameter(:σ_b_safep, sqrt(1/400)/4, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(100., sqrt(1/400)/4),
                   fixed=false,
                   description="σ_b_safep: Standard deviation of stationary component of safe asset preference shifter process.",
                   tex_label="\\sigma_{b^p, safe}")
end

function ss17!(m::Model1010_1spread)
    # ss13 with :σ_b_safep fixed at 0

    m <= parameter(:ρ_AAA, 0.5, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1),
                   fixed=false,
                   description="ρ_AAA: AR(1) coefficient in the AAA spread process.",
                   tex_label="\\rho_{AAA}")

    m <= parameter(:ρ_BBB, 0.5, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1),
                   fixed=false,
                   description="ρ_BBB: AR(1) coefficient in the BBB spread process.",
                   tex_label="\\rho_{BBB}")

    m <= parameter(:σ_BBB, 0.1, (1e-8, 5.),(1e-8, 5.),ModelConstructors.Exponential(),RootInverseGamma(2., 0.10),
                   fixed=false,
                   description="σ_BBB: Standard deviation on the AR(1) process for measurement error on the BBB spread.",
                   tex_label="\\sigma_{BBB}")

    m <= parameter(:σ_b_liqp, sqrt(1/400)/4, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(100., sqrt(1/400)/4),
                   fixed=false,
                   description="σ_b_liqp: Standard deviation of stationary component of liquid asset preference shifter process.",
                   tex_label="\\sigma_{b^p, liq}")

    m <= parameter(:σ_b_safep, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(100., sqrt(1/400)/4),
                   fixed=true,
                   description="σ_b_safep: Standard deviation of stationary component of safe asset preference shifter process.",
                   tex_label="\\sigma_{b^p, safe}")
end

function ss18!(m::Model1010_1spread)
    # ss13 with ρ_z_p fixed at 0.99

    m <= parameter(:ρ_BBB, 0.5, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1),
                   fixed=false,
                   description="ρ_BBB: AR(1) coefficient in the BBB spread process.",
                   tex_label="\\rho_{BBB}")

    m <= parameter(:σ_BBB, 0.1, (1e-8, 5.),(1e-8, 5.),ModelConstructors.Exponential(),RootInverseGamma(2., 0.10),
                   fixed=false,
                   description="σ_BBB: Standard deviation on the AR(1) process for measurement error on the BBB spread.",
                   tex_label="\\sigma_{BBB}")

    m <= parameter(:σ_b_p, sqrt(1/400)/4, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(100., sqrt(1/400)/4),
                   fixed=false,
                   description="σ_b_p: Standard deviation of stationary component of asset preference shifter process.",
                   tex_label="\\sigma_{b^p}")

    m <= parameter(:ρ_z_p,     0.99, (0.0, 1.0),      (0.0, 1.0), ModelConstructors.SquareRoot(),    BetaAlt(0.5, 0.2),
                   fixed=true, description="ρ_z_p: No description available.", tex_label="\\rho_{z^p}")
end

function ss19!(m::Model1010_1spread)
    # ss19 with a tight prior on σ_z_p

    m <= parameter(:ρ_AAA, 0.5, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1),
                   fixed=false,
                   description="ρ_AAA: AR(1) coefficient in the AAA spread process.",
                   tex_label="\\rho_{AAA}")

    m <= parameter(:ρ_BBB, 0.5, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1),
                   fixed=false,
                   description="ρ_BBB: AR(1) coefficient in the BBB spread process.",
                   tex_label="\\rho_{BBB}")

    m <= parameter(:σ_BBB, 0.1, (1e-8, 5.),(1e-8, 5.),ModelConstructors.Exponential(),RootInverseGamma(2., 0.10),
                   fixed=false,
                   description="σ_BBB: Standard deviation on the AR(1) process for measurement error on the BBB spread.",
                   tex_label="\\sigma_{BBB}")

    m <= parameter(:σ_b_liqp, sqrt(1/400)/4, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(100., sqrt(1/400)/4),
                   fixed=false,
                   description="σ_b_liqp: Standard deviation of stationary component of liquid asset preference shifter process.",
                   tex_label="\\sigma_{b^p, liq}")

    m <= parameter(:σ_b_safep, sqrt(1/400)/4, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(100., sqrt(1/400)/4),
                   fixed=false,
                   description="σ_b_safep: Standard deviation of stationary component of safe asset preference shifter process.",
                   tex_label="\\sigma_{b^p, safe}")

    m <= parameter(:ρ_z_p, 0.99, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.2),
                   fixed=true,
                   description="ρ_z_p: No description available.",
                   tex_label="\\rho_{z^p}")

    m <= parameter(:σ_z_p, sqrt(1/400)/4, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(100., sqrt(1/400)/4),
                   fixed=false,
                   description="σ_z_p: No description available.",
                   tex_label="\\sigma_{z^p}")
end

function ss20!(m::Model1010_1spread)
    # ss18 with betabar defined (correctly) with σ_c instead of σ_ω_star
    ss18!(m)
end

function ss21!(m::Model1010_1spread)
    ss13!(m)

    # standard deviations of the anticipated policy shocks
    for i = 1:n_anticipated_shocks_padding(m)
        m <= parameter(Symbol("σ_r_m$i"), .0, (1e-7, 100.), (1e-5, 0.), ModelConstructors.Exponential(), RootInverseGamma(4., .2), fixed=true,
                       description="σ_r_m$i: Standard deviation of the $i-period-ahead anticipated policy shock.",
                       tex_label=@sprintf("\\sigma_{%d,r}",i))
    end

end

function ss22!(m::Model1010_1spread)
    # ss20 but with inflation, long run inflation, etc. as pseudo-observables
    ss20!(m)
end

function ss23!(m::Model1010_1spread)
    # no m.e. on the Baa spread

    m <= parameter(:ρ_BBB, 0.5, (0.0, 1.0), (0.0, 1.0), ModelConstructors.SquareRoot(), BetaAlt(0.5, 0.1),
                   fixed=false,
                   description="ρ_BBB: AR(1) coefficient in the BBB spread process.",
                   tex_label="\\rho_{BBB}")

    m <= parameter(:σ_BBB, 0., (1e-8, 5.),(1e-8, 5.),ModelConstructors.Exponential(),RootInverseGamma(2., 0.10),
                   fixed=true,
                   description="σ_BBB: Standard deviation on the AR(1) process for measurement error on the BBB spread.",
                   tex_label="\\sigma_{BBB}")

    m <= parameter(:σ_b_p, sqrt(1/400)/4, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(100., sqrt(1/400)/4),
                   fixed=false,
                   description="σ_b_p: Standard deviation of stationary component of asset preference shifter process.",
                   tex_label="\\sigma_{b^p}")

    m <= parameter(:ρ_z_p,     0.99, (0.0, 1.0),      (0.0, 1.0), ModelConstructors.SquareRoot(),    BetaAlt(0.5, 0.2),
                   fixed=true, description="ρ_z_p: No description available.", tex_label="\\rho_{z^p}")
end


text:
using Gadfly, Compose

set_default_plot_size(4inch,9inch)

p=plot(y=[1,2,3], Geom.line)
c=Compose.context()
vstack(Union{Plot,Compose.Context}[p,c,p,c,p])

ext: jl


text:
using Poker
using Documenter

makedocs(;
    modules=[Poker],
    authors="itscoe",
    repo="https://github.com/itscoe/Poker.jl/blob/{commit}{path}#L{line}",
    sitename="Poker.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://itscoe.github.io/Poker.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/itscoe/Poker.jl.git",
)


text:
using Lithium
using Documenter

DocMeta.setdocmeta!(Lithium, :DocTestSetup, :(using Lithium); recursive=true)

makedocs(;
    modules=[Lithium],
    authors="Abraham Nunes",
    repo="https://github.com/cpsylab/Lithium.jl/blob/{commit}{path}#{line}",
    sitename="Lithium.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://cpsylab.github.io/Lithium.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/cpsylab/Lithium.jl",
)


text:
# Autogenerated wrapper script for acl_jll for armv6l-linux-gnueabihf
export chacl, getfacl, libacl, setfacl

using Attr_jll
JLLWrappers.@generate_wrapper_header("acl")
JLLWrappers.@declare_executable_product(chacl)
JLLWrappers.@declare_executable_product(getfacl)
JLLWrappers.@declare_library_product(libacl, "libacl.so.1")
JLLWrappers.@declare_executable_product(setfacl)
function __init__()
    JLLWrappers.@generate_init_header(Attr_jll)
    JLLWrappers.@init_executable_product(
        chacl,
        "bin/chacl",
    )

    JLLWrappers.@init_executable_product(
        getfacl,
        "bin/getfacl",
    )

    JLLWrappers.@init_library_product(
        libacl,
        "lib/libacl.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        setfacl,
        "bin/setfacl",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: applicationcostprofiler
using AWS.Compat
using AWS.UUIDs

"""
    delete_report_definition(report_id)
    delete_report_definition(report_id, params::Dict{String,<:Any})

Deletes the specified report definition in AWS Application Cost Profiler. This stops the
report from being generated.

# Arguments
- `report_id`: Required. ID of the report to delete.

"""
function delete_report_definition(
    reportId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return applicationcostprofiler(
        "DELETE",
        "/reportDefinition/$(reportId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_report_definition(
    reportId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return applicationcostprofiler(
        "DELETE",
        "/reportDefinition/$(reportId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_report_definition(report_id)
    get_report_definition(report_id, params::Dict{String,<:Any})

Retrieves the definition of a report already configured in AWS Application Cost Profiler.

# Arguments
- `report_id`: ID of the report to retrieve.

"""
function get_report_definition(reportId; aws_config::AbstractAWSConfig=global_aws_config())
    return applicationcostprofiler(
        "GET",
        "/reportDefinition/$(reportId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_report_definition(
    reportId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return applicationcostprofiler(
        "GET",
        "/reportDefinition/$(reportId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    import_application_usage(source_s3_location)
    import_application_usage(source_s3_location, params::Dict{String,<:Any})

Ingests application usage data from Amazon Simple Storage Service (Amazon S3). The data
must already exist in the S3 location. As part of the action, AWS Application Cost Profiler
copies the object from your S3 bucket to an S3 bucket owned by Amazon for processing
asynchronously.

# Arguments
- `source_s3_location`: Amazon S3 location to import application usage data from.

"""
function import_application_usage(
    sourceS3Location; aws_config::AbstractAWSConfig=global_aws_config()
)
    return applicationcostprofiler(
        "POST",
        "/importApplicationUsage",
        Dict{String,Any}("sourceS3Location" => sourceS3Location);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function import_application_usage(
    sourceS3Location,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return applicationcostprofiler(
        "POST",
        "/importApplicationUsage",
        Dict{String,Any}(
            mergewith(
                _merge, Dict{String,Any}("sourceS3Location" => sourceS3Location), params
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_report_definitions()
    list_report_definitions(params::Dict{String,<:Any})

Retrieves a list of all reports and their configurations for your AWS account. The maximum
number of reports is one.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to return.
- `"nextToken"`: The token value from a previous call to access the next page of results.
"""
function list_report_definitions(; aws_config::AbstractAWSConfig=global_aws_config())
    return applicationcostprofiler(
        "GET", "/reportDefinition"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function list_report_definitions(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return applicationcostprofiler(
        "GET",
        "/reportDefinition",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    put_report_definition(destination_s3_location, format, report_description, report_frequency, report_id)
    put_report_definition(destination_s3_location, format, report_description, report_frequency, report_id, params::Dict{String,<:Any})

Creates the report definition for a report in Application Cost Profiler.

# Arguments
- `destination_s3_location`: Required. Amazon Simple Storage Service (Amazon S3) location
  where Application Cost Profiler uploads the report.
- `format`: Required. The format to use for the generated report.
- `report_description`: Required. Description of the report.
- `report_frequency`: Required. The cadence to generate the report.
- `report_id`: Required. ID of the report. You can choose any valid string matching the
  pattern for the ID.

"""
function put_report_definition(
    destinationS3Location,
    format,
    reportDescription,
    reportFrequency,
    reportId;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return applicationcostprofiler(
        "POST",
        "/reportDefinition",
        Dict{String,Any}(
            "destinationS3Location" => destinationS3Location,
            "format" => format,
            "reportDescription" => reportDescription,
            "reportFrequency" => reportFrequency,
            "reportId" => reportId,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function put_report_definition(
    destinationS3Location,
    format,
    reportDescription,
    reportFrequency,
    reportId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return applicationcostprofiler(
        "POST",
        "/reportDefinition",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "destinationS3Location" => destinationS3Location,
                    "format" => format,
                    "reportDescription" => reportDescription,
                    "reportFrequency" => reportFrequency,
                    "reportId" => reportId,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    update_report_definition(destination_s3_location, format, report_description, report_frequency, report_id)
    update_report_definition(destination_s3_location, format, report_description, report_frequency, report_id, params::Dict{String,<:Any})

Updates existing report in AWS Application Cost Profiler.

# Arguments
- `destination_s3_location`: Required. Amazon Simple Storage Service (Amazon S3) location
  where Application Cost Profiler uploads the report.
- `format`: Required. The format to use for the generated report.
- `report_description`: Required. Description of the report.
- `report_frequency`: Required. The cadence to generate the report.
- `report_id`: Required. ID of the report to update.

"""
function update_report_definition(
    destinationS3Location,
    format,
    reportDescription,
    reportFrequency,
    reportId;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return applicationcostprofiler(
        "PUT",
        "/reportDefinition/$(reportId)",
        Dict{String,Any}(
            "destinationS3Location" => destinationS3Location,
            "format" => format,
            "reportDescription" => reportDescription,
            "reportFrequency" => reportFrequency,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function update_report_definition(
    destinationS3Location,
    format,
    reportDescription,
    reportFrequency,
    reportId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return applicationcostprofiler(
        "PUT",
        "/reportDefinition/$(reportId)",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "destinationS3Location" => destinationS3Location,
                    "format" => format,
                    "reportDescription" => reportDescription,
                    "reportFrequency" => reportFrequency,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
using RasterDataSources, URIs, Test, Dates
using RasterDataSources: rastername, rasterpath, zipurl, zipname, zippath

@testset "AWAP" begin
    using RasterDataSources: rastername, zipurl, zipname, zippath

    raster_file = joinpath(ENV["RASTERDATASOURCES_PATH"], "AWAP", "vprp", "vprph09", "20010101.grid")
    @test rasterpath(AWAP, :vprpress09; date=Date(2001, 1)) == raster_file
    @test rastername(AWAP, :vprpress09; date=Date(2001, 1)) == "20010101.grid"

    @test zipurl(AWAP, :vprpress09; date=Date(2001, 1)) ==
        URI(scheme="http", host="www.bom.gov.au", path="/web03/ncc/www/awap/vprp/vprph09/daily/grid/0.05/history/nat/2001010120010101.grid.Z")
    @test zippath(AWAP, :vprpress09; date=Date(2001, 1)) ==
        joinpath(ENV["RASTERDATASOURCES_PATH"], "AWAP", "vprp", "vprph09", "20010101.grid.Z")
    @test zipname(AWAP, :vprpress09; date=Date(2001, 1)) == "20010101.grid.Z"

    if Sys.islinux()
        @test getraster(AWAP, :vprpress09; date=DateTime(2001, 01, 01)) == raster_file
        @test getraster(AWAP, (:vprpress09,); date=DateTime(2001, 01, 01)) == (vprpress09=raster_file,)
        @test getraster(AWAP, [:vprpress09]; date=DateTime(2001, 01, 01)) == (vprpress09=raster_file,)
        @test isfile(raster_file)
    end
end

max_forks_repo_name: cesaraustralia/RasterDataSources


text:
"""
Based on:
https://medium.com/vooban-ai/hyperopt-tutorial-for-optimizing-neural-networks-hyperparameters-e3102814b919
"""
module TestFindMin

using Test
using TreeParzen

function objective_min(space)
    x = space[:x]
    y = space[:y]

    return x^2 + y^2
end

space = Dict(
    :x => HP.Uniform(:x, -5.0, 5.0),
    :y => HP.Uniform(:y, -5.0, 5.0)
)

best = fmin(objective_min, space, 100)
@test isapprox(best[:x], 0.013181950926553512, rtol = 1e2)
@test isapprox(best[:y], 0.0364742933684085, rtol = 1e2)

end
true


text:
@testset "Is integral" begin
  Qx, x = FlintQQ["x"]
  f = x^2 + 1
  K, a = number_field(f, "a")
  Ky, y = K["y"]

  g = y^3 + 3
  L, b = number_field(g, "b")

  @test Hecke.isintegral(b) == true
  @test Hecke.isintegral(fmpq(1, 2)*b) == false

  h = y^4 + 3
  M, c = number_field([g, h], "c")

  @test Hecke.isintegral(c[1]) == true
  @test Hecke.isintegral(fmpq(1, 2)*c[1]) == false
end


text:
# test the time restriction methods
t = 0:.25:10
sig = Signal(t, [t])
b1 = before(sig, 5.0s)
@test b1.time == 0s:0.25s:5s
@test all(b1[1] .== 0:0.25:5)
a1 = after(sig, 5.0s)
@test a1.time == 5s:.25s:10s
@test all(a1[1] .== 5:.25:10)
w1 = within(sig, 2s, 3s)
@test w1.time == 2s:.25s:3s
@test all(w1[1] .== 2:.25:3)

# Matrix Signal
t = 0:.25:10
sig = Signal(t, [t 2t])
b1 = before(sig, 5.0s)
@test b1.time == 0s:0.25s:5s
@test all(b1[1] .== 0:0.25:5)
@test all(b1[2] .== 2*(0:0.25:5))
a1 = after(sig, 5.0s)
@test a1.time == 5s:.25s:10s
@test all(a1[1] .== 5:.25:10)
@test all(a1[2] .== 2*(5:.25:10))
w1 = within(sig, 2s, 3s)
@test w1.time == 2s:.25s:3s
@test all(w1[1] .== 2:.25:3)
@test all(w1[2] .== 2*(2:.25:3))


text:
##  Copyright 2021 Myrl Marmarelis
##
##  Licensed under the Apache License, Version 2.0 (the "License");
##  you may not use this file except in compliance with the License.
##  You may obtain a copy of the License at
##
##    http://www.apache.org/licenses/LICENSE-2.0
##
##  Unless required by applicable law or agreed to in writing, software
##  distributed under the License is distributed on an "AS IS" BASIS,
##  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
##  See the License for the specific language governing permissions and
##  limitations under the License.

# the following is a timeline of my original research..

# see "The Soft Multivariate Truncated Normal Distribution", circa 2018.
# a Gibbs sampler for the linearly-constrained/truncated multivariate normal,
# without the approximated cutoff, could probably be figured out as well.

# for implementational simplicity, and probably because our dimensionality is
# not high enough to warrant the above, we go with the procedure described in
# "Efficient Sampling Methods for Truncated Multivariate Normal and Student-t
# Distributions Subject to Linear Inequality Constraints" (2015), wherein the
# straightforward Gibbs approach is validated theoretically and fleshed out.

# "Exact Hamiltonian Monte Carlo for Truncated Multivariate Gaussians" (2014)
# with a Gibbs step in between Hamiltonian iterations to possibly switch
# between mixture components. All components are truncated the same way, and
# the partition function takes a value over the entire mixture---not per
# component; hence some may have much lower total mass than others. Can I
# therefore form a categorical likelihood over normalized component masses,
# and sample the resultant Dirichlet-multinomial (n=1) posterior? I think so.
# See "Gibbs sampling for fitting finite and infinite Gaussian mixture models"
# from 2013. Though it is not directly useful.

using StatsFuns

# p(membership|utility) = p(membership) * p(utility|membership) / (sum over possible memberships)
# used for sampling the mixture, not any parameter posterior.
# this sampler could have trouble mixing if cohort spreads are not wide enough (via large precision_scale and dof?)
function sample_realization_membership(utility::Utility{N,T},
    mixture::VoterMixture{N,M,T})::Int where {N,M,T}
  log_prior = log.(mixture.shares)
  log_likelihood = SVector{M}( logpdf(cohort, utility)
    for cohort in mixture.cohorts )
  energy = log_prior + log_likelihood
  log_probability = energy .- logsumexp(energy) # for numerical stability?
  exp.(log_probability) |> Categorical |> rand
end

# for walls, g is positive but basically zero (for strict inequality constraints,)
# and F is populated with {-1, 0, 1}. I might want a tiny `threshold` somewhere in here.
# further, we must add Sum(f_j^i mu_i) to g_j to account for the centered version. paper was unclear here.
# `walls` could possibly have fewer columns `R`, but we simply pad with zeros
# `sample` is an allocated staging area for the velocity vector.
# fully deterministic..
function simulate_particle_trajectory!(sample::AbstractVector{T},
    utility::Utility{N,T}, mean_utility::SVector{N,T}, inv_cov::SMatrix{N,N,T,L}, # `inv_cov` for this mixture component
    cohort::VoterCohort{N,T,L}, walls::SMatrix{N,N,Int,L}, remaining_time::T,
    collision_limit::Int, small_gap::T=T(1e-5))::Tuple{Utility{N,T}, T} where {N,T,L}
  rand!(cohort, sample) # ^ needs to be like 1e-2 for Float32, although said type appears to be slower for some reason... (casting at runtime?)
  velocity = SVector{N,T}(sample) .- mean_utility
  position = utility
  leftover_time = T(Inf)
  for t in 1:collision_limit
    momentum::SVector = inv_cov * velocity
    alpha = velocity
    beta = position - mean_utility
    #no_collision = iszero.(walls * alpha) .& iszero.(walls * beta) # Eq. 2.24 in the paper, simply check for any nonzero per j
    condition = sqrt.( (walls * alpha).^2 + (walls * beta).^2 )
    shift = walls * mean_utility
    collision_time = T(Inf)
    collision_wall = @SVector zeros(T, N)
    for wall in 1:N
      if condition[wall] < abs(shift[wall]) #no_collision[wall]
        continue
      end
      coefficients = walls[wall, :] # should be an SVector
      tangent = -(coefficients' * alpha) / (coefficients' * beta)
      phase = atan(tangent)
      # cos(t + phase) = 0 for t = pi/2 - phase + n*pi. find smallest n for t>0
      # cos(t + phase) = -shift/condition for t = acos(-shift/condition) - phase + n*pi
      target = -shift[wall] / condition[wall]
      time = acos(target) - phase - 2pi # MAKE SURE we don't start too large, by subtracting an extra 2pi
      while time <= 0 # todo a better closed-form solution?
        time += pi
      end
      if time < collision_time
        collision_time = time
        collision_wall = coefficients
      end
    end
    @assert collision_time > 0
    collision_time *= 1 - small_gap # leave an air gap so that we are for sure not overstepping the boundary
    eval_time = min(remaining_time, collision_time)
    position = mean_utility + alpha*sin(eval_time) + beta*cos(eval_time)
    next_velocity = alpha*cos(eval_time) - beta*sin(eval_time)
    leftover_time = remaining_time - collision_time
    if leftover_time <= 0
      break
    end
    reflection = (
      (collision_wall' * next_velocity) / sum(abs2, collision_wall) )
    reflected_velocity = next_velocity - 2reflection * collision_wall
    # Julia cares not for tail-call optimization, so I unwound into a loop.
    # alas, the recursion was much more elegant.
    velocity = reflected_velocity
    remaining_time = leftover_time
  end
  position, max(leftover_time, 0)
end

# actually, we have an `n_candidates-1` maximum number of active walls per issue,
# so they will never occupy the entire `N` columns. but the impact should be marginal.
function build_walls(::Val{N}, voter_index::Int, vote::MultiIssueVote,
    indifference::Bool)::SMatrix{N,N,Int} where N
  construction = @MMatrix zeros(Int, N, N) # transposed for memory efficiency
  seen = @MVector fill(false, N)
  wall = 1
  offset = 0
  for issue in vote
    ranking = issue.choices[voter_index]
    n_ranked = length(ranking)
    @assert ranking[1] > 0
    last_ranking = ranking[1]
    seen[last_ranking + offset] = true # in case the below iterand skips entirely
    for constraint in 2:n_ranked
      if ranking[constraint] == 0
        continue
      end
      greater = last_ranking
      lesser = ranking[constraint]
      construction[greater+offset, wall] = 1
      construction[lesser+offset, wall] = -1
      seen[greater+offset] = seen[lesser+offset] = true
      last_ranking = lesser
      wall += 1
    end
    if !indifference
      for candidate in 1:issue.n_candidates
        if seen[candidate+offset]
          continue
        end
        construction[last_ranking+offset, wall] = 1
        construction[candidate+offset, wall] = -1
        wall += 1
      end
    end
    offset += issue.n_candidates
  end
  convert(SVector, construction')
end

# use `voters` as MCMC initializations for the Hamiltonian part of the overall,
# hierarchical, Gibbs sampler for the posterior.
# initialize either with a call to `sample_utilities!`, or by zeroing utilities
# across the board. the latter works because we aren't doing strict inequalities. (yet)
# walls do not couple different issues on an absolute scale, so it is important to
# restrain the means by having a large `mean_scale`.
function hamiltonian_sample_utilities!(voters::Vector{VoterRealization{N,M,T}},
    vote::MultiIssueVote, mixture::VoterMixture{N,M,T},
    n_trajectories::Int, collision_limit::Int, indif::Bool=false) where {N,M,T}
  trajectory_length = pi/2 |> T
  n_voters = length(voters)
  n_issues = length(vote)
  @assert isapprox(sum(mixture.shares), 1)
  check_issues(vote, n_voters, N)
  # `inv(cohort.Σ)` does preserve PDMat structure
  inv_covs = [ invcov(cohort) |> SMatrix{N,N,T} for cohort in mixture.cohorts ] # all this comprehension *could* be an @SVector
  n_threads = nthreads()
  sample = zeros(T, N, n_threads) # intermediary for initial velocity
  incompletenesses = zeros(T, n_threads)
  @threads for voter_index in 1:n_voters
    walls = build_walls(Val(N), voter_index, vote, indif)
    thread_id = threadid()
    this_sample = @view sample[:, thread_id]
    voter = voters[voter_index]
    utility = voter.utility
    membership = voter.membership
    for trajectory in 1:n_trajectories
      cohort = mixture.cohorts[membership]
      cohort_mean = mean(cohort) |> SVector{N}
      utility, leftover_time = simulate_particle_trajectory!(this_sample,
        utility, cohort_mean, inv_covs[membership], cohort, walls,
        trajectory_length, collision_limit) # as prescribed by the paper
      membership = sample_realization_membership(utility, mixture)
      incompleteness = leftover_time / trajectory_length # aggregate and average these..
      incompletenesses[thread_id] += incompleteness
    end
    realization = VoterRealization{N,M,T}(membership, utility)
    voters[voter_index] = realization
  end
  incompleteness = sum(incompletenesses) / (n_trajectories * n_voters)
  incompleteness
end

using Random: shuffle!

# is this guaranteed to converge? it should, methinks..
function shoehorn_ranking(utility::Utility{N,T},
    voter_index::Int, vote::MultiIssueVote, indif::Bool)::Utility{N,T} where {N,T}
  malleable = convert(MVector, utility)
  walls = build_walls(Val(N), voter_index, vote, indif)
  search_order = 1:N |> collect
  satisfied = false
  while !satisfied
    satisfied = true
    shuffle!(search_order)
    for wall_index in search_order
      wall = walls[wall_index, :]
      if wall' * malleable < 0 # pair is wrong
        pair = findall( wall .!= 0 )
        @assert length(pair) == 2
        temp = malleable[pair[1]]
        malleable[pair[1]] = malleable[pair[2]]
        malleable[pair[2]] = temp
        satisfied = false
      end
    end
  end
  convert(SVector, malleable)
end


text:
@test SBA.calc_epsilon([1 2 3; 3 1 2], [4 3 2; 5 0 1]) == [-3 -1 1; -2 1 1]

A = Array{Float64}(undef, 2, 4, 4)

A[:, :, 1] = [ 1  3  0 -1;
               2 -1 -1  4]
A[:, :, 2] = [ 3  0  4  1;
               2  1 -1 -3]
A[:, :, 3] = [-1  4  2  0;
               3  1  3  3]
A[:, :, 4] = [ 1  0  3  2;
               0  2  1  1]

AtA = Array{Float64}(undef, 4, 4, 4)
for index in 1:size(A, 3)
    AtA[:, :, index] = A[:, :, index]' * A[:, :, index]
end

B = Array{Float64}(undef, 2, 3, 4)

B[:, :, 1] = [ 2  3  1;
               1  0  2]
B[:, :, 2] = [ 1  0  2;
               1  1  3]
B[:, :, 3] = [ 1  2  3;
               1  0  1]
B[:, :, 4] = [ 0  0  5;
               4  0  3]

BtB = Array{Float64}(undef, 3, 3, 4)
for index in 1:size(B, 3)
    BtB[:, :, index] = B[:, :, index]' * B[:, :, index]
end

# 2 points, 3 viewpoints
# i = 1:2, j = 1:3
# N = 4 (number of visible points)
# X    = [x_11 x_12 x_22 x_23]
# mask = [   1    1    0;
#            0    1    1]

viewpoint_indices = [1, 2, 2, 3]
point_indices = [1, 1, 2, 2]

indices = SBA.Indices(viewpoint_indices, point_indices)

U = SBA.calc_U(indices, A)
@test size(U) == (4, 4, 3)  # (n_pose_paramas, n_pose_paramas, n_viewpoints)
@test U[:, :, 1] == AtA[:, :, 1]
@test U[:, :, 2] == AtA[:, :, 2] + AtA[:, :, 3]
@test U[:, :, 3] == AtA[:, :, 4]

V_inv = SBA.calc_V_inv(indices, B)
@test size(V_inv) == (3, 3, 2)  # (n_point_params, n_point_params, n_points)
@test V_inv[:, :, 1] == inv(BtB[:, :, 1] + BtB[:, :, 2])
@test V_inv[:, :, 2] == inv(BtB[:, :, 3] + BtB[:, :, 4])

W = SBA.calc_W(indices, A, B)
@test size(W) == (4, 3, 4)  # (n_pose_params, n_point_params, N)
@test W[:, :, 1] == A[:, :, 1]' * B[:, :, 1]
@test W[:, :, 2] == A[:, :, 2]' * B[:, :, 2]
@test W[:, :, 3] == A[:, :, 3]' * B[:, :, 3]
@test W[:, :, 4] == A[:, :, 4]' * B[:, :, 4]

epsilon = [ 1  2 -3  1;
            3 -1  3  0]

epsilon_a = SBA.calc_epsilon_a(indices, A, epsilon)
@test size(epsilon_a) == (4, 3)  # (n_pose_params, n_viewpoints)
@test epsilon_a[:, 1] == A[:, :, 1]' * epsilon[:, 1]
@test epsilon_a[:, 2] == A[:, :, 2]' * epsilon[:, 2] +
                         A[:, :, 3]' * epsilon[:, 3]
@test epsilon_a[:, 3] == A[:, :, 4]' * epsilon[:, 4]

epsilon_b = SBA.calc_epsilon_b(indices, B, epsilon)
@test size(epsilon_b) == (3, 2)  # (n_point_params, n_points)
@test epsilon_b[:, 1] == B[:, :, 1]' * epsilon[:, 1] +
                         B[:, :, 2]' * epsilon[:, 2]
@test epsilon_b[:, 2] == B[:, :, 3]' * epsilon[:, 3] +
                         B[:, :, 4]' * epsilon[:, 4]

Y = SBA.calc_Y(indices, W, V_inv)
@test size(Y) == (4, 3, 4)  # (n_pose_params, n_point_params, N)
@test Y[:, :, 1] == W[:, :, 1] * V_inv[:, :, 1]  # (i, j) = (1, 1)
@test Y[:, :, 2] == W[:, :, 2] * V_inv[:, :, 1]  # (i, j) = (1, 2)
@test Y[:, :, 3] == W[:, :, 3] * V_inv[:, :, 2]  # (i, j) = (2, 2)
@test Y[:, :, 4] == W[:, :, 4] * V_inv[:, :, 2]  # (i, j) = (2, 3)

e = SBA.calc_e(indices, Y, epsilon_a, epsilon_b)
@test size(e) == (4, 3)  # (n_pose_params, n_viewpoints)
@test e[:, 1] ==  - Y[:, :, 1] * epsilon_b[:, 1] + epsilon_a[:, 1]
@test e[:, 2] == (- Y[:, :, 2] * epsilon_b[:, 1] - Y[:, :, 3] * epsilon_b[:, 2]
                  + epsilon_a[:, 2])
@test e[:, 3] ==  - Y[:, :, 4] * epsilon_b[:, 2] + epsilon_a[:, 3]

S = SBA.calc_S(indices, U, Y, W)
# (n_pose_params * n_viewpoints, n_pose_params * n_viewpoints)
@test size(S) == (4 * 3, 4 * 3)
# (j, k) == (1, 1)
@test S[1: 4, 1: 4] ==  - Y[:, :, 1] * W[:, :, 1]' + U[:, :, 1]
# (j, k) == (1, 2)
@test S[1: 4, 5: 8] ==  - Y[:, :, 1] * W[:, :, 2]'
# (j, k) == (1, 3)
@test S[1: 4, 9:12] == zeros(Float64, 4, 4)
# (j, k) == (2, 1)
@test S[5: 8, 1: 4] ==  - Y[:, :, 2] * W[:, :, 1]'
# (j, k) == (2, 2)
@test S[5: 8, 5: 8] == (- Y[:, :, 2] * W[:, :, 2]' - Y[:, :, 3] * W[:, :, 3]'
                        + U[:, :, 2])
# (j, k) == (2, 3)
@test S[5: 8, 9:12] ==  - Y[:, :, 3] * W[:, :, 4]'
# (j, k) == (3, 1)
@test S[9:12, 1: 4] == zeros(Float64, 4, 4)
# (j, k) == (3, 2)
@test S[9:12, 5: 8] ==  - Y[:, :, 4] * W[:, :, 3]'
# (j, k) == (3, 3)
@test S[9:12, 9:12] ==  - Y[:, :, 4] * W[:, :, 4]' + U[:, :, 3]

delta_a = SBA.calc_delta_a(S, e)
@test size(delta_a) == (4, 3)  # (n_pose_params, n_viewpoints)
@test isapprox(S * vec(delta_a), vec(e))

delta_b = SBA.calc_delta_b(indices, V_inv, W, epsilon_b, delta_a)

@test size(delta_b) == (3, 2)  # (n_point_params, n_points)
@test isapprox(delta_b[:, 1],
               V_inv[:, :, 1] * (epsilon_b[:, 1]
                              - (W[:, :, 1]' * delta_a[:, 1] +
                                 W[:, :, 2]' * delta_a[:, 2])))
@test isapprox(delta_b[:, 2],
               V_inv[:, :, 2] * (epsilon_b[:, 2]
                              - (W[:, :, 3]' * delta_a[:, 2] +
                                 W[:, :, 4]' * delta_a[:, 3])))


text:
module SparseFactorizations

include("ldlt_api.jl")
include("selinv_api.jl")

export
    ldlt_structure,
    level_of_fill,
    simplicial_ldlt,
    ildlt,
    ildlt_lof,
    ildlt_tol,
    selinv

end # module

max_stars_repo_name: ettersi/SparseFactorizations

max_issues_repo_name: ettersi/SparseFactorizations

max_forks_repo_name: ettersi/SparseFactorizations


text:
#=
    ^, pow, root

    exp, expm1, log, log1p, log2, log10, logbase,
    sin, sinpi, cos, cospi, tan, tanpi, cot, cotpi,
    sinh, cosh, tanh, coth,
    asin, acos, atan, asinh, acosh, atanh,
    sincos, sincospi, sinhcosh,
    sinc,
    gamma, lgamma, zeta
=#

for (op,cfunc) in ((:exp,:arb_exp), (:expm1, :arb_expm1),
    (:log,:arb_log), (:log1p, :arb_log1p),
    (:sin, :arb_sin), (:sinpi, :arb_sinpi), (:cos, :arb_cos), (:cospi, :arb_cospi),
    (:tan, :arb_tan), (:cot, :arb_cot),
    (:sinh, :arb_sinh), (:cosh, :arb_cosh), (:tanh, :arb_tanh), (:coth, :arb_coth),
    (:asin, :arb_asin), (:acos, :arb_asin), (:atan, :arb_atan),
    (:asinh, :arb_asinh), (:acosh, :arb_asinh), (:atanh, :arb_atanh),
    (:sinc, :arb_sinc),
    (:gamma, :arb_gamma), (:lgamma, :arb_lgamma), (:zeta, :arb_zeta),
    (:digamma, :arb_digamma), (:rgamma, :arb_rgamma)
    )
  @eval begin
    function ($op){P}(x::ArbFloat{P})
      z = ArbFloat{P}()
      ccall(@libarb($cfunc), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}, Int), &z, &x, P)
      z
    end
  end
end

function logbase{P}(x::ArbFloat{P}, base::Int)
    b = UInt(abs(base))
    z = ArbFloat{P}()
    ccall(@libarb(arb_log_base_ui), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}, UInt, Int), &z, &x, b, P)
    z
end

log2{P}(x::ArbFloat{P}) = logbase(x, 2)
log10{P}(x::ArbFloat{P}) = logbase(x, 10)

for (op,cfunc) in ((:sincos, :arb_sin_cos), (:sincospi, :arb_sin_cos_pi), (:sinhcosh, :arb_sinh_cosh))
  @eval begin
    function ($op){P}(x::ArbFloat{P})
        sz = ArbFloat{P}()
        cz = ArbFloat{P}()
        ccall(@libarb($cfunc), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}, Ptr{ArbFloat}, Int), &sz, &cz, &x, P)
        sz, cz
    end
  end
end

function atan2{P}(a::ArbFloat{P}, b::ArbFloat{P})
    z = ArbFloat{P}()
    ccall(@libarb(arb_atan2), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}, Ptr{ArbFloat}, Int), &z, &a, &b, P)
    z
end

for (op,cfunc) in ((:root, :arb_root_ui),)
  @eval begin
    function ($op){P}(x::ArbFloat{P}, y::UInt)
      z = ArbFloat{P}()
      ccall(@libarb($cfunc), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}, UInt, Int), &z, &x, &y, P)
      z
    end
  end
end

for (op,cfunc) in ((:^,:arb_pow), (:pow,:arb_pow))
  @eval begin
    function ($op){P,I<:Integer}(x::I, y::ArbFloat{P})
      xx = ArbFloat{P}(x)
      z = ArbFloat{P}()
      ccall(@libarb($cfunc), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}, Ptr{ArbFloat}, Int), &z, &xx, &y, P)
      return z
    end
    function ($op){P,I<:Integer}(x::ArbFloat{P}, y::I)
      sy,ay = signbit(y), abs(y)
      yy = ArbFloat{P}(ay)
      z = ArbFloat{P}()
      ccall(@libarb($cfunc), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}, Ptr{ArbFloat}, Int), &z, &x, &yy, P)
      return sy ? 1/z : z
    end
    function ($op){P,R<:Rational}(x::R, y::ArbFloat{P})
      xx = ArbFloat{P}(x)
      z = ArbFloat{P}()
      ccall(@libarb($cfunc), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}, Ptr{ArbFloat}, Int), &z, &xx, &y, P)
      return z
    end
    function ($op){P,R<:Rational}(x::ArbFloat{P}, y::R)
      yy = ArbFloat{P}(y)
      z = ArbFloat{P}()
      ccall(@libarb($cfunc), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}, Ptr{ArbFloat}, Int), &z, &x, &yy, P)
      return z
    end
    function ($op){P,R<:Real}(x::R, y::ArbFloat{P})
      xx = ArbFloat{P}(x)
      z = ArbFloat{P}()
      ccall(@libarb($cfunc), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}, Ptr{ArbFloat}, Int), &z, &xx, &y, P)
      return z
    end
    function ($op){P,R<:Real}(x::ArbFloat{P}, y::R)
      yy = ArbFloat{P}(y)
      z = ArbFloat{P}()
      ccall(@libarb($cfunc), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}, Ptr{ArbFloat}, Int), &z, &x, &yy, P)
      return z
    end
    function ($op){P}(x::ArbFloat{P}, y::ArbFloat{P})
      z = ArbFloat{P}()
      ccall(@libarb($cfunc), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}, Ptr{ArbFloat}, Int), &z, &x, &y, P)
      z
    end
    ($op){P,Q}(x::ArbFloat{P}, y::ArbFloat{Q}) = ($op)(promote(x,y)...)
  end
end

root{P}(x::ArbFloat{P}, y::ArbFloat{P}) = pow(x, inv(y))
root{P}(x::Integer, y::ArbFloat{P}) = pow(ArbFloat{P}(x), inv(y))
function root{P}(x::ArbFloat{P}, y::Integer)
   return (
     if y>=0
       yy = UInt64(y)
       z = ArbFloat{P}()
       ccall(@libarb(arb_root_ui), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}, Ptr{UInt64}, Int), &z, &x, &yy, P)
       z
    else
      pow(ArbFloat{P}(x), inv(y))
    end )
end


text:
@testset "Atomic measures" begin
    Mod.@polyvar x y
    δ1 = WeightedDiracMeasure([1, 0],     2.0)
    δ2 = WeightedDiracMeasure([1/2, 1/2], 3.0)
    η = AtomicMeasure([x, y], [δ1, δ2])
    @test η ≈ AtomicMeasure([x, y], [δ2, δ1])
    @test !(η ≈ AtomicMeasure([x, y], [WeightedDiracMeasure([1, 0],     3.0),
                                       WeightedDiracMeasure([1/2, 1/2], 2.0)]))
    @test string(η) == "Atomic measure on the variables x, y with 2 atoms:\n at [1.0, 0.0] with weight 2.0\n at [0.5, 0.5] with weight 3.0"
    p = x^2 + x + y + x*y^2 + 1
    @test dot(η, p) == 2 * p(x=>1, y=>0) + 3 * p(x => 1/2, y => 1/2)
    @test dot(η, p) == dot(p, η)
    μ = measure(η, monomials(p))
    @test dot(μ, p) == 2 * p(x=>1, y=>0) + 3 * p(x => 1/2, y => 1/2)
end


text:
export QRscaling, updateQRscaling,
                  fqrnegative_log_likelihood!,
                  fqrnegative_log_likelihood,
                  gqrnegative_log_likelihood!,
                  gqrnegative_log_likelihood,
                  qrnegative_log_likelihood!,
                  qrnegative_log_likelihood,
                  qrprecond!
                  # qrprecond

"""
$(TYPEDEF)

An immutable structure to hold a preconditioner and its Cholesky factorization.

## Fields
$(TYPEDFIELDS)

"""
struct QRscaling
    R::UpperTriangular{Float64,Array{Float64,2}}
    Rinv::UpperTriangular{Float64,Array{Float64,2}}
    D::Diagonal{Float64, Array{Float64,1}}
    Dinv::Diagonal{Float64, Array{Float64,1}}
    U::UpperTriangular{Float64,Array{Float64,2}}
    Uinv::UpperTriangular{Float64,Array{Float64,2}}
    L2Uinv::Array{Float64,2}
end

function QRscaling(R, ψnorm)
    D = Diagonal(ψnorm)
    Dinv = inv(D)
    R = UpperTriangular(R)
    Rinv = inv(R)
    U = UpperTriangular(R*D)
    Uinv = Dinv*Rinv
    L2Uinv = Uinv'*Uinv
    return QRscaling(R, Rinv, D, Dinv, U, Uinv, L2Uinv)
end

function QRscaling(S::Storage)
    D = Diagonal(S.ψnorm)
    Dinv = inv(D)
    R = UpperTriangular(qr(S.ψoffψd * Dinv).R)
    Rinv = inv(R)
    U = UpperTriangular(R*D)
    Uinv = Dinv*Rinv
    L2Uinv = Uinv'*Uinv
    return QRscaling(R, Rinv, D, Dinv, U, Uinv, L2Uinv)
end

"""
$(TYPEDSIGNATURES)

Updates the `QRscaling` variable `F` with the appended column of the `Storage` variable `S`.
"""
function updateQRscaling(F::QRscaling, S::Storage)
    Nψ = S.Nψ
    # Update F with the last column of QR scaling
    D = Diagonal(vcat(F.D.diag, S.ψnorm[Nψ]))
    Dinv = Diagonal(vcat(F.Dinv.diag, 1 ./ S.ψnorm[Nψ]))
    mul!(S.ψoffψd, S.ψoffψd, Dinv)
    R = UpperTriangular(qraddcol(view(S.ψoffψd,:,1:Nψ-1), F.R, view(S.ψoffψd,:, Nψ)))
    mul!(S.ψoffψd, S.ψoffψd, D)
    # Use Schur complement to efficiently compute the inverse of the UpperDiagonal Matrix
    # https://en.wikipedia.org/wiki/Schur_complement
    Rinv = UpperTriangular(hcat(vcat(F.Rinv.data, zeros(1, Nψ-1)), vcat(-1.0/R[Nψ,Nψ]*F.Rinv*view(R,1:Nψ-1,Nψ), 1.0/R[Nψ,Nψ])))

    U = UpperTriangular(hcat(vcat(F.U.data, zeros(1, Nψ-1)), S.ψnorm[Nψ]*view(R,:,Nψ)))
    # Use Schur complement to efficiently compute the inverse of the UpperDiagonal Matrix
    Uinv = UpperTriangular(hcat(vcat(F.Uinv.data, zeros(1, Nψ-1)), vcat(-1.0/U[Nψ,Nψ]*F.Uinv*view(U,1:Nψ-1,Nψ), 1.0/U[Nψ,Nψ])))

    L2Uinvlower = F.Uinv'*view(Uinv, 1:Nψ-1, Nψ)
    L2Uinv = hcat(vcat(F.L2Uinv, L2Uinvlower'), vcat(L2Uinvlower, 0.0))
    L2Uinv[end,end] = norm(view(Uinv,:,Nψ))^2
    F = QRscaling(R, Rinv, D, Dinv, U, Uinv, L2Uinv)
    return F
end

"""
$(TYPEDSIGNATURES)

Computes in-place the negative log likelihood and its gradient for `c̃oeff`
(rescaled coefficients) for the `HermiteMapComponent` `C` and for the ensemble matrix `X`.
"""
function qrnegative_log_likelihood!(J̃, dJ̃, c̃oeff, F::QRscaling, S::Storage, C::HermiteMapComponent, X)
    # In this version, c̃oeff is expressed in the rescaled space
    NxX, Ne = size(X)
    m = C.m
    Nx = C.Nx
    Nψ = C.Nψ
    # @assert NxX == Nx "Wrong dimension of the sample X"
    # @assert size(S.ψoff, 1) == Ne
    # @assert size(S.ψoff, 2) == Nψ

    # Output objective, gradient
    xlast = view(X,NxX,:)
    # Get back to the coeff in the unscaled space
    mul!(c̃oeff, F.Uinv, c̃oeff)

    fill!(S.cache_integral, 0.0)

    # Integrate at the same time for the objective, gradient
    function integrand!(v::Vector{Float64}, t::Float64)
        repeated_grad_xk_basis!(S.cache_dcψxdt, S.cache_gradxd, C.I.f, t*xlast)

         # This computing is also reused in the computation of the gradient, no interest to skip it
        @avx @. S.cache_dcψxdt *= S.ψoff
        # mul!(S.cache_dcψxdt, S.cache_dcψxdt, F.Uinv)
        mul!(S.cache_dψxd, S.cache_dcψxdt, c̃oeff)

        # Integration for J̃
        vJ = view(v,1:Ne)
        evaluate!(vJ, C.I.g, S.cache_dψxd)

        # Integration for dcJ̃
        grad_x!(S.cache_dψxd, C.I.g, S.cache_dψxd)

        v[Ne+1:Ne+Ne*Nψ] = reshape(S.cache_dψxd .* (S.cache_dcψxdt), (Ne*Nψ))
    end

    quadgk!(integrand!, S.cache_integral, 0.0, 1.0)#; order = 9, rtol = 1e-10)

    # Multiply integral by xlast (change of variable in the integration)
    @inbounds for j=1:Nψ+1
        @. S.cache_integral[(j-1)*Ne+1:j*Ne] *= xlast
    end

    # Add f(x_{1:d-1},0) i.e. (S.ψoffψd0 .* S.ψd0)*coeff to S.cache_integral
    mul!(view(S.cache_integral,1:Ne),S.ψoffψd0, c̃oeff, 1.0, 1.0)

    # Store g(∂_{xk}f(x_{1:k})) in S.cache_g
    mul!(S.cache_g, S.ψoffdψxd, c̃oeff)

    # Return back to the scaled space for the coefficients
    mul!(c̃oeff, F.U, c̃oeff)

    # Formatting to use with Optim.jl
    if dJ̃ != nothing
        reshape_cacheintegral = reshape(S.cache_integral[Ne+1:end], (Ne, Nψ))
        fill!(dJ̃, 0.0)
        @inbounds for i=1:Ne
            for j=1:Nψ
            dJ̃[j] += gradlog_pdf(S.cache_integral[i])*(reshape_cacheintegral[i,j] + S.ψoffψd0[i,j]) + # dsoftplus(S.cache_g[i])*S.ψoff[i,j]*S.dψxd[i,j]*(1/softplus(S.cache_g[i]))
                     grad_x(C.I.g, S.cache_g[i])*S.ψoffdψxd[i,j]/C.I.g(S.cache_g[i])
            end
        end
        # Use chain rule ∂J/∂c̃ = ∂J/∂c ∂c/∂c̃ = ∂J/∂c U^{-1} = U^{-T} ∂J/∂c
        lmul!(F.Uinv', dJ̃)
        # Add derivative of the L2 penalty term ∂_c α ||U^{-1}c||^2 = 2 α U^{-1}c
        mul!(dJ̃, Symmetric(F.L2Uinv), c̃oeff, 2*C.α, -1/Ne)
    end

    if J̃ != nothing
        J̃ = 0.0
        @avx for i=1:Ne
            J̃ += log_pdf(S.cache_integral[i]) + log(C.I.g(S.cache_g[i]))
        end
        J̃ *=(-1/Ne)
        J̃ += C.α*norm(F.Uinv*c̃oeff)^2
        # J̃ = 0.0
        return J̃
    end
end

"""
$(TYPEDSIGNATURES)

Computes the negative log likelihood and its gradient for `c̃oeff`
(in the rescaled coefficient space) for the `HermiteMapComponent` `C` and for the ensemble matrix `X`.
"""
qrnegative_log_likelihood(F::QRscaling, S::Storage, C::HermiteMapComponent, X) = (J̃, dJ̃, c̃oeff) -> qrnegative_log_likelihood!(J̃, dJ̃, c̃oeff, F, S, C, X)

# F.Uinv'*precond(c)*F.Uinv

# qrprecond! expects c̃oeff to be expressed in the QR space
"""
$(TYPEDSIGNATURES)

Computes in-place  Hessian preconditioner based on the Gauss-Newton approximation (outer-product of gradient of the cost function) in the rescaled space.
"""
function qrprecond!(P, c̃oeff, F::QRscaling, S::Storage, C::HermiteMapComponent, X)
    Nψ = C.Nψ
    NxX, Ne = size(X)
    @assert NxX == C.Nx "Wrong dimension of the sample X"
    @assert size(S.ψoff, 1) == Ne
    @assert size(S.ψoff, 2) == Nψ

    # Output objective, gradient
    xlast = view(X,NxX,:)

    # Get back to the coeff in the unscaled space
    mul!(c̃oeff, F.Uinv, c̃oeff)

    fill!(S.cache_integral, 0)

    # Integrate at the same time for the objective, gradient
    function integrand!(v::Vector{Float64}, t::Float64)
        repeated_grad_xk_basis!(S.cache_dcψxdt, S.cache_gradxd, C.I.f, t*xlast)

        # @avx @. S.cache_dψxd = (S.cache_dcψxdt .* S.ψoff) *ˡ coeff
        @avx @. S.cache_dcψxdt *= S.ψoff
        # mul!(S.cache_dcψxdt, S.cache_dcψxdt, F.Uinv)
        mul!(S.cache_dψxd, S.cache_dcψxdt, c̃oeff)

        # Integration for J
        vJ = view(v,1:Ne)
        evaluate!(vJ, C.I.g, S.cache_dψxd)

        # Integration for dcJ

        grad_x!(S.cache_dψxd, C.I.g, S.cache_dψxd)

        # @avx for j=2:Nψ+1
        #     for i=1:Ne
        #         v[(j-1)*Ne+i] = S.cache_dψxd[i]*S.cache_dcψxdt[i,j-1]
        #     end
        # end

        v[Ne+1:Ne+Ne*Nψ] .= reshape(S.cache_dψxd .* S.cache_dcψxdt , (Ne*Nψ))
        # v[Ne+1:Ne+Ne*Nψ] .= reshape(grad_x(C.I.g, S.cache_dψxd) .* S.cache_dcψxdt , (Ne*Nψ))
    end

    quadgk!(integrand!, S.cache_integral, 0.0, 1.0)#; order = 9, rtol = 1e-10)

    # Multiply integral by xk (change of variable in the integration)
    @inbounds for j=1:Nψ+1
        @. S.cache_integral[(j-1)*Ne+1:j*Ne] *= xlast
    end
    # @avx for j=1:Nψ+1
    #     for i=1:Ne
    #         S.cache_integral[(j-1)*Ne+i] *= xlast[i]
    #     end
    # end

    # Add f(x_{1:d-1},0) i.e. (S.ψoff .* S.ψd0)*coeff to S.cache_integral
    @avx for i=1:Ne
        f0i = zero(Float64)
        for j=1:Nψ
            f0i += S.ψoffψd0[i,j]*c̃oeff[j]
        end
        S.cache_integral[i] += f0i
    end

    # Store g(∂_{xk}f(x_{1:k})) in S.cache_g
    @avx for i=1:Ne
        prelogJi = zero(Float64)
        for j=1:Nψ
            prelogJi += S.ψoffdψxd[i,j]*c̃oeff[j]
        end
        S.cache_g[i] = prelogJi
    end

    reshape_cacheintegral = reshape(S.cache_integral[Ne+1:Ne+Ne*Nψ], (Ne, Nψ))

    # Return back to the scaled space for the coefficients
    mul!(c̃oeff, F.U, c̃oeff)

    fill!(P, 0.0)
    @inbounds for l=1:Ne
        # Exploit symmetry of the Hessian
        for i=1:Nψ
            for j=i:Nψ
            # P[i,j] +=  reshape2_cacheintegral[l,i,j]*S.cache_integral[l]
            P[i,j] +=  (reshape_cacheintegral[l,i] + S.ψoffψd0[l,i]) * (reshape_cacheintegral[l,j] + S.ψoffψd0[l,j])
            P[i,j] -=  ( (S.ψoffdψxd[l,i]) * (S.ψoffdψxd[l,j])*(
                            hess_x(C.I.g, S.cache_g[l]) * C.I.g(S.cache_g[l]) -
                            grad_x(C.I.g, S.cache_g[l])^2))/C.I.g(S.cache_g[l])^2

            P[j,i] = P[i,j]
            end
        end
    end
    rmul!(P, 1/Ne)
    # Add derivative of the L2 penalty term ∂^2_c̃ α ||Uinv c̃||^2 = ∂^2_c̃ (α c̃' Uinv' Uinv c̃) = 2*α Uinv'*Uinv
    lmul!(F.Uinv', P)
    rmul!(P, F.Uinv)

    P .+= 2*C.α*F.L2Uinv
end

"""
$(TYPEDSIGNATURES)

Computes in-place  Hessian preconditioner based on the Gauss-Newton approximation (outer-product of gradient of the cost function) in the rescaled space.
"""
qrprecond!(S::Storage, C::HermiteMapComponent, X) = (P, c̃oeff) -> qrprecond!(P, c̃oeff, S, C, X)
#

#
# function fqrnegative_log_likelihood!(c̃oeff, F::QRscaling, S::Storage, C::HermiteMapComponent, X)
#     # In this version, c̃oeff is expressed in the rescaled space
#     NxX, Ne = size(X)
#     m = C.m
#     Nx = C.Nx
#     Nψ = C.Nψ
#
#     # Output objective, gradient
#     xlast = view(X,NxX,:)
#
#     fill!(S.cache_integral, 0.0)
#
#     # Integrate at the same time for the objective, gradient
#     function integrand!(v::Vector{Float64}, t::Float64)
#         repeated_grad_xk_basis!(S.cache_dcψxdt, S.cache_gradxd, C.I.f, t*xlast)
#
#          # This computing is also reused in the computation of the gradient, no interest to skip it
#         @avx @. S.cache_dcψxdt *= S.ψoff
#         mul!(S.cache_dcψxdt, S.cache_dcψxdt, F.Uinv)
#         mul!(S.cache_dψxd, S.cache_dcψxdt, c̃oeff)
#
#         # Integration for J̃
#         vJ = view(v,1:Ne)
#         evaluate!(vJ, C.I.g, S.cache_dψxd)
#     end
#
#     quadgk!(integrand!, S.cache_integral[1:Ne], 0.0, 1.0; rtol = 1e-3)#; order = 9, rtol = 1e-10)
#
#     # Multiply integral by xlast (change of variable in the integration)
#     @avx for i=1:Ne
#         S.cache_integral[i] *= xlast[i]
#     end
#
#     # Add f(x_{1:d-1},0) i.e. (S.ψoffψd0 .* S.ψd0)*coeff to S.cache_integral
#     mul!(view(S.cache_integral,1:Ne),S.ψoffψd0, c̃oeff, 1.0, 1.0)
#
#     # Store g(∂_{xk}f(x_{1:k})) in S.cache_g
#     mul!(S.cache_g, S.ψoffdψxd, c̃oeff)
#
#     J̃ = 0.0
#     @avx for i=1:Ne
#         J̃ += log_pdf(S.cache_integral[i]) + log(C.I.g(S.cache_g[i]))
#     end
#     J̃ *=(-1/Ne)
#     J̃ += C.α*norm(F.Uinv*c̃oeff)^2
#     return J̃
# end
#
# fqrnegative_log_likelihood(F::QRscaling, S::Storage, C::HermiteMapComponent, X) = (c̃oeff) -> fqrnegative_log_likelihood!(c̃oeff, F, S, C, X)
#
# function gqrnegative_log_likelihood!(dJ̃, c̃oeff, F::QRscaling, S::Storage, C::HermiteMapComponent, X)
#     # In this version, c̃oeff is expressed in the rescaled space
#     NxX, Ne = size(X)
#     m = C.m
#     Nx = C.Nx
#     Nψ = C.Nψ
#     # @assert NxX == Nx "Wrong dimension of the sample X"
#     # @assert size(S.ψoff, 1) == Ne
#     # @assert size(S.ψoff, 2) == Nψ
#
#     # Output objective, gradient
#     xlast = view(X,NxX,:)
#
#     fill!(S.cache_integral, 0.0)
#
#     # Integrate at the same time for the objective, gradient
#     function integrand!(v::Vector{Float64}, t::Float64)
#         repeated_grad_xk_basis!(S.cache_dcψxdt, S.cache_gradxd, C.I.f, t*xlast)
#
#          # This computing is also reused in the computation of the gradient, no interest to skip it
#         @avx @. S.cache_dcψxdt *= S.ψoff
#         mul!(S.cache_dcψxdt, S.cache_dcψxdt, F.Uinv)
#         mul!(S.cache_dψxd, S.cache_dcψxdt, c̃oeff)
#
#         # Integration for J̃
#         vJ = view(v,1:Ne)
#         evaluate!(vJ, C.I.g, S.cache_dψxd)
#
#         # Integration for dcJ̃
#         grad_x!(S.cache_dψxd, C.I.g, S.cache_dψxd)
#
#         v[Ne+1:Ne+Ne*Nψ] .= reshape(S.cache_dψxd .* S.cache_dcψxdt , (Ne*Nψ))
#     end
#
#     quadgk!(integrand!, S.cache_integral, 0.0, 1.0; rtol = 1e-3, order = 3)#; order = 9, rtol = 1e-10)
#
#     # Multiply integral by xlast (change of variable in the integration)
#     @inbounds for j=1:Nψ+1
#         @. S.cache_integral[(j-1)*Ne+1:j*Ne] *= xlast
#     end
#
#     # Add f(x_{1:d-1},0) i.e. (S.ψoffψd0 .* S.ψd0)*coeff to S.cache_integral
#     mul!(view(S.cache_integral,1:Ne),S.ψoffψd0, c̃oeff, 1.0, 1.0)
#
#     # Store g(∂_{xk}f(x_{1:k})) in S.cache_g
#     mul!(S.cache_g, S.ψoffdψxd, c̃oeff)
#
#     # Formatting to use with Optim.jl
#     reshape_cacheintegral = reshape(S.cache_integral[Ne+1:end], (Ne, Nψ))
#     fill!(dJ̃, 0.0)
#     @inbounds for i=1:Ne
#         for j=1:Nψ
#         dJ̃[j] += gradlog_pdf(S.cache_integral[i])*(reshape_cacheintegral[i,j] + S.ψoffψd0[i,j]) + # dsoftplus(S.cache_g[i])*S.ψoff[i,j]*S.dψxd[i,j]*(1/softplus(S.cache_g[i]))
#                  grad_x(C.I.g, S.cache_g[i])*S.ψoffdψxd[i,j]/C.I.g(S.cache_g[i])
#         end
#     end
#     # Add derivative of the L2 penalty term ∂_c α ||U^{-1}c||^2 = 2 α U^{-1}c
#     mul!(dJ̃, F.L2Uinv, c̃oeff, 2*C.α, -1/Ne)
#     # dJ̃
#     nothing
# end
#
# gqrnegative_log_likelihood(F::QRscaling, S::Storage, C::HermiteMapComponent, X) = (dJ̃, c̃oeff) -> gqrnegative_log_likelihood!(dJ̃, c̃oeff, F, S, C, X)
#


text:
# Imports
using EpistemicNetworkAnalysis
using DataFrames
using CSV
using Plots
using GLM
using Statistics
using MixedModels
using HypothesisTests

# Data
RSData = ena_dataset("RS.data")

# Specification
conversations = [:Condition, :GameHalf, :GroupName]
units = [:Condition, :GameHalf, :UserName]
codes = [
    :Data,
    :Technical_Constraints,
    :Performance_Parameters,
    :Client_and_Consultant_Requests,
    :Design_Reasoning,
    :Collaboration
]

# Rotations
means = MeansRotation(:Condition, "FirstGame", "SecondGame", flipsvd=true)
means2 = MeansRotation(:GameHalf, "First", "Second")
interaction = EpistemicNetworkAnalysis.Formula2Rotation(
    LinearModel, 2, @formula(y ~ 1 + Condition + GameHalf + Condition&GameHalf), Dict(:Condition => EffectsCoding(), :GameHalf => EffectsCoding()),
    LinearModel, 2, @formula(y ~ 1 + GameHalf), nothing #Dict(:CondGroupName => EffectsCoding())
)

# Artist
tv = EpistemicNetworkAnalysis.TVRemoteArtist(
    :Condition, "FirstGame", "SecondGame",
    :GameHalf, "First", "Second"
)

tv2 = EpistemicNetworkAnalysis.TVRemoteArtist(
    :GameHalf, "First", "Second",
    :Condition, "FirstGame", "SecondGame"
)

# Model 1
ena = ENAModel(
    RSData, codes, conversations, units,
    rotateBy=means, windowSize=5
)

display(ena)
p = plot(
    ena,
    artist=tv,
    xaxisname="First/Second Game",
    yaxisname="SVD",
    displayFilter=RSRow -> RSRow[:Condition] in ["FirstGame"]
)

xlims!(p, -1, +1)
ylims!(p, -1, +1)
title!(p, "Means Rotation, First Game")
savefig(p, "first-game.svg")
savefig(p, "first-game.png")
run(`inkscape first-game.svg -E first-game.eps --export-ignore-filters --export-ps-level=3`)
display(p)

p = plot(
    ena,
    artist=tv,
    xaxisname="First/Second Game",
    yaxisname="SVD",
    displayFilter=RSRow -> RSRow[:Condition] in ["SecondGame"]
)

xlims!(p, -1, +1)
ylims!(p, -1, +1)
title!(p, "Means Rotation, Second Game")
savefig(p, "second-game.svg")
savefig(p, "second-game.png")
run(`inkscape second-game.svg -E second-game.eps --export-ignore-filters --export-ps-level=3`)
display(p)

p = plot(
    ena,
    artist=tv,
    xaxisname="First/Second Game",
    yaxisname="SVD"
)

xlims!(p, -1, +1)
ylims!(p, -1, +1)
title!(p, "Means Rotation, Subtraction")
savefig(p, "subtraction.svg")
savefig(p, "subtraction.png")
run(`inkscape subtraction.svg -E subtraction.eps --export-ignore-filters --export-ps-level=3`)
display(p)

# Model 2
ena = ENAModel(
    RSData, codes, conversations, units,
    rotateBy=interaction, windowSize=5,
)

display(ena)
p = plot(
    ena,
    artist=tv,
    xaxisname="First/Second Game",
    yaxisname="First/Second Half",
    displayFilter=RSRow -> RSRow[:Condition] in ["FirstGame"]
)

xlims!(p, -1, +1)
ylims!(p, -1, +1)
title!(p, "Custom Rotation, First Game")
savefig(p, "interaction-first.svg")
savefig(p, "interaction-first.png")
run(`inkscape interaction-first.svg -E interaction-first.eps --export-ignore-filters --export-ps-level=3`)
display(p)

p = plot(
    ena,
    artist=tv,
    xaxisname="First/Second Game",
    yaxisname="First/Second Half",
    displayFilter=RSRow -> RSRow[:Condition] in ["SecondGame"]
)

xlims!(p, -1, +1)
ylims!(p, -1, +1)
title!(p, "Custom Rotation, Second Game")
savefig(p, "interaction-second.svg")
savefig(p, "interaction-second.png")
run(`inkscape interaction-second.svg -E interaction-second.eps --export-ignore-filters --export-ps-level=3`)
display(p)

p = plot(
    ena,
    artist=tv,
    xaxisname="First/Second Game",
    yaxisname="First/Second Half"
)

xlims!(p, -1, +1)
ylims!(p, -1, +1)
title!(p, "Custom Rotation, Subtraction")
savefig(p, "interaction.svg")
savefig(p, "interaction.png")
run(`inkscape interaction.svg -E interaction.eps --export-ignore-filters --export-ps-level=3`)
display(p)

# Model 3
ena = ENAModel(
    RSData, codes, conversations, units,
    rotateBy=means2, windowSize=5
)

display(ena)
p = plot(
    ena,
    artist=tv2,
    xaxisname="First/Second Half",
    yaxisname="SVD",
    displayFilter=RSRow -> RSRow[:GameHalf] in ["First"]
)

xlims!(p, -1, +1)
ylims!(p, -1, +1)
title!(p, "Means Rotation, First Half")
savefig(p, "first-half.svg")
savefig(p, "first-half.png")
run(`inkscape first-half.svg -E first-half.eps --export-ignore-filters --export-ps-level=3`)
display(p)

p = plot(
    ena,
    artist=tv2,
    xaxisname="First/Second Half",
    yaxisname="SVD",
    displayFilter=RSRow -> RSRow[:GameHalf] in ["Second"]
)

xlims!(p, -1, +1)
ylims!(p, -1, +1)
title!(p, "Means Rotation, Second Half")
savefig(p, "second-half.svg")
savefig(p, "second-half.png")
run(`inkscape second-half.svg -E second-half.eps --export-ignore-filters --export-ps-level=3`)
display(p)

p = plot(
    ena,
    artist=tv2,
    xaxisname="First/Second Half",
    yaxisname="SVD"
)

xlims!(p, -1, +1)
ylims!(p, -1, +1)
title!(p, "Means Rotation, Subtraction, Halves")
savefig(p, "subtraction2.svg")
savefig(p, "subtraction2.png")
run(`inkscape subtraction2.svg -E subtraction2.eps --export-ignore-filters --export-ps-level=3`)
display(p)


text:
module BertWordPieceTokenizer

# import dependencies
import HTTP
import IsURL

# export functions and variables
export init, tokenize, enable_truncation, encode, decode, token_to_id, id_to_token

# include local julia scripts
include("utils.jl")

# define global variables
# initialized vocabulary
global const _token2id = Ref{Dict{String, Int}}(Dict{String, Int}())
global const _id2token = Ref{Dict{Int, String}}(Dict{Int, String}())
global const _max_word_length = Ref{Int}(200)
global const _do_lowercase = Ref{Bool}(false)
# max tokenized length
global const _max_length = Ref{Union{Int, Nothing}}(nothing)
# truncation direction, post or pre
global const _truncation = Ref{String}("post")
# special tokens
global const _PAD = "[PAD]"
global const _UNK = "[UNK]"
global const _MASK = "[MASK]"
global const _CLS = "[CLS]"
global const _SEP = "[SEP]"
global const _SPECIAL_TOKENS = Vector{String}([_PAD, _UNK, _MASK, _CLS, _SEP])

function init(vocab_path::String; cache_path::Union{String, Nothing} = nothing, do_lowercase::Bool = false, max_word_length::Union{Int, Nothing} = nothing)
    #= Initialize tokenizer
    Args:
      vocab_path: String, local path or http url to vocabulary
      cache_path: Union{String, Nothing}, specify the path to cache remote file in local
      do_lowercase: do lowercase
    =#
    global _token2id[] = load_vocab(vocab_path, cache_path=cache_path)
    global _id2token[] = Dict{Int, String}((v, k) for (k, v) in _token2id[])
    global _do_lowercase[] = do_lowercase
    if max_word_length != nothing
        global _max_word_length[] = max_word_length
    end
end

function enable_truncation(max_length::Int; truncation::String = "post")
    #= Enable truncation
    Args:
      max_length: Int, specify max token length
      truncation: String, specify truncation from ["post", "pre"]
    =#
    @assert truncation in ["pre", "post"] "please specify truncation from [`pre`, `post`]"
    global _max_length[] = max_length
    global _truncation[] = truncation
end

function decode(token_ids::Vector{Int})::String
    #= Recover text from token indices
    Args:
      token_ids: Vector{Int}
    =#
    tokens = [id_to_token(id) for id in token_ids]
    tokens = [token for token in tokens if token ∉ _SPECIAL_TOKENS]
    text = ""
    for (i, token) in enumerate(tokens)
        if startswith(token, "##")
            text *= token[3:end]
        elseif length(token) == 1 && is_cjk_character(first(token))
            text *= token
        elseif length(token) == 1 && ispunct(first(token))
            text *= token * " "
        elseif i > 1 && is_cjk_character(first(text[end]))
            text *= token
        else
            text *= " " * token
        end
    end

    return strip(text)
end

function encode(first_text::String; second_text::Union{String, Nothing} = nothing)::Tuple{Vector{Int}, Vector{Int}}
    #= Encode text to token indices and segment indices
    Args:
      first_text: String, first text
      second_text: Union{String, Nothing} = nothing, second text
    =#
    tokens = tokenize(first_text, second_text=second_text)
    token_ids = []
    segment_ids = []
    segment_value = 0
    for token in tokens
        push!(token_ids, token_to_id(token))
        push!(segment_ids, segment_value)
        if token == _SEP
            segment_value = 1
        end
    end
    return token_ids, segment_ids
end

function tokenize(first_text::String; second_text::Union{String, Nothing} = nothing)::Vector{String}
    #= Tokenize text
    Args:
      first_text: String, first text
      second_text: Union{String, Nothing} = nothing, second text
    =#

    first_tokens = _tokenize(first_text)
    pushfirst!(first_tokens, _CLS)

    if second_text != nothing
        second_tokens = _tokenize(second_text)
        pushfirst!(second_tokens, _SEP)
        texts = [first_tokens, second_tokens]
    else
        texts = [first_tokens]
    end

    if _max_length[] != nothing
        max_length = _max_length[] - 1
        if _truncation[] == "post"
            offset = 0
        elseif _truncation[] == "pre"
            offset = 1
        end
        truncate_text!(texts, max_length, offset, _truncation[])
    end

    if second_text != nothing
        append!(first_tokens, second_tokens)
    end
    push!(first_tokens, _SEP)
    return first_tokens
end

function _tokenize(text::String)::Vector{String}
    #= Tokenize text
    Args:
      text: String, input text
    =#
    if _do_lowercase[]
        text = Base.lowercase(text)
    end

    tmp = ""
    for ch in text
        ch = first(ch)  # to Char
        if ispunct(ch) || is_cjk_character(ch)
            tmp *= " " * ch * " "
        elseif isspace(ch)
            tmp *= " "
        elseif Int(ch) == 0 || Int(ch) == 65533 || iscntrl(ch)
            continue
        else
            tmp *= ch
        end
    end

    tokens = Vector{String}()
    for word in split(strip(tmp))
        append!(tokens, _wordpiece_tokenize(String(word)))
    end

    return tokens
end

function _wordpiece_tokenize(word::String)::Vector{String}
    #= Word piece tokenize
    Args:
      word: String, a word
    =#
    if length(word) > _max_word_length[]
        return [word]
    end

    # to list of string
    word = split(word, "")

    tokens = Vector{String}()
    i = j = 0
    while i < length(word)
        j = length(word)
        sub = nothing
        while j > i
            sub = join(word[i+1:j])
            if i > 0
                sub = "##" * sub
            end
            if haskey(_token2id[], sub)
                break
            end
            j -= 1
        end
        if i == j
            j += 1
        end
        if sub != nothing
            push!(tokens, sub)
        end
        i = j
    end

    return tokens
end

function token_to_id(token::String)::Int
    #= Convert token to the corresponding index
    Args:
      token: String
    =#
    if haskey(_token2id[], token)
        return _token2id[][token]
    end
    return _token2id[][_UNK]
end

function id_to_token(id::Int)::String
    #= Convert index to the corresponding token
    Args:
      id: Int
    =#
    if haskey(_id2token[], id)
        return _id2token[][id]
    end
    return _UNK
end

function load_vocab(vocab_path::String; cache_path::Union{String, Nothing} = nothing)::Dict{String, Int}
    #= Load vocabulary from local file or remote http url
    Args:
      vocab_path: String, path to vocabulary
      cache_path: Union{String, nothing} = nothing, specify the path to cache remote file in local
    =#
    vocab = Dict{String, Int}()

    if IsURL.isurl(vocab_path)
        println("downloading vocabulary from $vocab_path")
        r = HTTP.get(vocab_path)
        @assert r.status == 200
        lines = strip(String(r.body))
        if cache_path != nothing
            open(cache_path, "w") do writer
                write(writer, lines)
            end
        end
        for line in split(lines, "\n")
            token = split(line)
            if isempty(token)
                token = strip(line)
            else
                token = token[1]
            end
            vocab[token] = length(vocab)
        end
    else
        open(vocab_path) do reader
            for line in eachline(reader)
                token = split(line)
                if isempty(token)
                    token = strip(line)
                else
                    token = token[1]
                end
                vocab[token] = length(vocab)
            end
        end
    end
    return vocab
end

end # module


text:
using StaticArrays
using ConstrainedDynamics

# Parameters
ex = [1.;0.;0.]
ey = [0.;1.;0.]
ez = [0.;0.;1.]

l1 = 2.0
x, y = .1, .1
r = l1 / 4
h = r / 10
box = Box(x, l1, y, l1, color = RGBA(1., 1., 0.))
cyl = Cylinder(r, h, 2 * r, color = RGBA(1., 0., 0.))

p0 = [0;l1 / 2;0]
p1 = -p0
p2 = [0;r;0]
p3 = -p2

# Initial orientation
q1 = UnitQuaternion(RotY(pi / 2))

# Links
origin = Origin{Float64}()
link1 = Body(box)
link2 = Body(cyl)

# Constraints
joint0to12 = EqualityConstraint(CylindricalFree(origin, link1, ey; p2=p0), Revolute(origin, link2, ez; p1=p3 + p1))
joint1to2 = EqualityConstraint(Cylindrical(link1, link2, ez; p1=p1, p2=p3))

links = [link1; link2]
constraints = [joint0to12;joint1to2]
shapes = [box,cyl]

function slider_crank_control!(mechanism, k)
    τ = SA[0;0.3]
    setForce!(mechanism, joint1to2, τ * 0)
    return
end

mech = Mechanism(origin, links, constraints, shapes = shapes)
setPosition!(link1,x = -p0)
setPosition!(link2,x = p3 + p1)

# a = [0;0;0;0;0.3]

# setForce!(mech,joint0to12,a)
# setForce!(link2,τ=[0;0;0.3])


text:
using IssueReporter
using Base.Test
using URIParser, GitHub

delete!(ENV, "GITHUB_ACCESS_TOKEN")
secrets_path = joinpath(Pkg.dir("IssueReporter"), "secrets.jl")
isfile(secrets_path) && mv(secrets_path, secrets_path * "_")

@testset "Basic features" begin
  @testset "Looking up an existing package returns a proper repo URI" begin
    @test IssueReporter.package_uri("DataFrames") |> URIParser.isvalid
  end

  @testset "Looking up a non-existent package should throw an error" begin
    @test_throws ErrorException IssueReporter.package_uri("DataFramess")
  end
end

@testset "GitHub integration" begin
  @testset "An undefined token should return false" begin
    @test ! IssueReporter.token_is_defined()
  end

  @testset "Attempting to access a token that is not set will error" begin
    @test_throws ErrorException IssueReporter.token()
  end

  # setup a mock token
  ENV["GITHUB_ACCESS_TOKEN"] = "1234"

  @testset "Token is defined" begin
    @test IssueReporter.token_is_defined()
  end

  @testset "A valid token is a non empty string and has the set value" begin
    token = IssueReporter.token()
    @test isa(token, String) && ! isempty(token)
    @test token == "1234"
  end

  delete!(ENV, "GITHUB_ACCESS_TOKEN")
  open(secrets_path, "w") do io
    write(io, """ENV["GITHUB_ACCESS_TOKEN"] = "1234" """)
  end

  @testset "Token is loaded from the `secrets.jl` file" begin
    @test IssueReporter.token_is_defined()
    @test IssueReporter.token() == "1234"
  end

  rm(secrets_path)

  @testset "The configure method writes the token to the secrets file" begin
    IssueReporter.configure("1234")
    @test readline(secrets_path) == """ENV["GITHUB_ACCESS_TOKEN"] = "1234" """
    rm(secrets_path)
  end
end

isfile(secrets_path * "_") && mv(secrets_path * "_", secrets_path)

delete!(ENV, "GITHUB_ACCESS_TOKEN")
@testset "Adding GitHub issues" begin
  @testset "Successful authentication should return a GitHub.OAuth2 instance" begin
    @test_skip isa(IssueReporter.github_auth(), GitHub.OAuth2)
  end
  @testset "Converting package name to GitHub id" begin
    @test IssueReporter.repo_id("Trie") == "JuliaArchive/Trie.jl"
  end
  @testset "Submitting an issue should result in a GitHub.Issue object" begin
    @test_skip isa(IssueReporter.report("Trie", "I found a bug", "Here is how you can reproduce the problem: ..."), GitHub.Issue)
  end
end


text:
using DSP, DelimitedFiles, LsqFit

function DT_model(t, β::Array{Float64,1}, ρ::Float64, ndet::Float64, nmed::Float64)
	#optimzed semi-infinite model

    n::Float64 = nmed/ndet
    μa::Float64 = β[1]
    μsp::Float64 = β[2]
    D::Float64 = 1/3μsp
	ν::Float64 = 29.9792345/nmed

    Rt = Array{Float64}(undef, length(t))

	if n > 1.0
		A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
		 5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
	elseif n < 1.0
		A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
	else
		A = 1.0
	end

	zs::Float64 = 1/μsp
	ze::Float64 = 2A*D

	z3m::Float64 = - zs
    z4m::Float64 = 2ze +zs

	Threads.@threads for n in eachindex(t)

   		Rt[n] = -exp(-(ρ^2/(4D*ν*t[n])) - μa*ν*t[n])
   		Rt[n] = Rt[n]/(2*(4π*D*ν)^(3/2)*t[n]^(5/2))

   		Rt[n] = Rt[n]*(z3m*exp(-(z3m^2/(4D*ν*t[n]))) - z4m*exp(-(z4m^2/(4D*ν*t[n]))))

	end

   replace!(Rt, NaN => 0)
   return Rt./maximum(Rt)

end

function get_fit_window(counts, risefactor, tailfactor)

    maxval = maximum(counts)
    ind1::Int = findfirst(x -> x > maxval*risefactor, counts)
    ind2::Int =  findlast(x -> x >maxval*tailfactor, counts)

    return ind1, ind2
end

function pads(u::Array{Float64,1})
    outsize = 2 .* length(u) .- 1
    out = zeros(eltype(u), nextfastfft(outsize))
    copyto!(out, u)

end

function conv_DT(t, β::Array{Float64,1}, data::input_data, out, fftIRF, pl, pli;
    RtDT = Array{Float64}(undef, length(data.IRF))
    )

        RtDT = DT_model(data.t, β, data.ρ, data.nmed, data.ndet)
        copyto!(out, RtDT)

        fftRtDT = pl * out
        fftRtDT .*= fftIRF
        out = pli * fftRtDT

        out = out./maximum(out)

	    tidx = findfirst(x -> x == t[1], data.t)

	    out = out[tidx:tidx+length(t)-1]
	return log.(out)
end

function getfit(input_data, model_params)
    ## get FFT padding and plans
    out = pads(input_data.IRF)

    pl = plan_rfft(out)
    fftIRF = pl * out
    pli = plan_irfft(fftIRF, length(out))

    #get fit windows
    ind1, ind2 = get_fit_window(input_data.DTOF, model_params.risefactor, model_params.tailfactor)

    curve_fit((t, β) -> conv_DT(t, β, input_data, out, fftIRF, pl, pli), input_data.t[ind1:ind2], log.(input_data.DTOF[ind1:ind2]),
    model_params.initparams, lower=model_params.lb, upper=model_params.ub)

end

#=
@benchmark getfit(data,modelp)
BenchmarkTools.Trial:
  memory estimate:  21.43 MiB
  allocs estimate:  5717
  --------------
  minimum time:     8.924 ms (0.00% GC)
  median time:      9.845 ms (0.00% GC)
  mean time:        12.204 ms (4.42% GC)
  maximum time:     54.299 ms (3.45% GC)
  --------------
  samples:          410
  evals/sample:     1

this is if you use only 1024 points
    data1 = input_data(tDTOF[1:1024], RtDTOF[500:1524], RtIRF[500:1524], rand(1024), 2.56, 1.5, 1.51, 700, false)

julia> @benchmark getfit(data1,modelp)
BenchmarkTools.Trial:
  memory estimate:  6.13 MiB
  allocs estimate:  5501
  --------------
  minimum time:     3.154 ms (0.00% GC)
  median time:      3.261 ms (0.00% GC)
  mean time:        5.035 ms (3.86% GC)
  maximum time:     88.410 ms (0.00% GC)
  --------------
  samples:          993
  evals/sample:     1

using DelimitedFiles, LsqFit, DSP, FFTW
filename1 = "/home/heltonmc/Desktop/DTOF.asc"
filename = "/home/heltonmc/Desktop/IRF.asc"

  data_array = readdlm(filename,skipstart=10)
data_array1 = readdlm(filename1,skipstart=10)

RtIRF = data_array[1:end-1,2]
RtIRF = RtIRF./maximum(RtIRF)
tIRF = data_array[1:end-1,1]

RtDTOF = data_array1[1:end-1,2]
RtDTOF = RtDTOF./maximum(RtDTOF)
tDTOF = data_array1[1:end-1,1]

struct input_data
    t::Array{Float64,1}
	DTOF::Array{Float64,1}
	IRF::Array{Float64,1}
	DTOFerrors::Array{Float64,1}
	ρ::Float64
	nmed::Float64
	ndet::Float64
	lambda::Float64
	normpeaks::Bool
end

struct model_params
	#yerrors::Array{Float64,1}          #uncertainties in dependent variable
    model::Function                    #model to be fitted against
	initparams::Array{Float64,1}
	lb::Array{Float64,1}
	ub::Array{Float64,1}
	risefactor::Float64  			#percent of peak to fit on rising edge
	tailfactor::Float64				#perfecnt of peak to fit on failing tail
	dof::Int64
end

struct fitresult
    #Fit inputs
    xdata::Array{Float64,1}     #independent variable
    ydata::Array{Float64,1}     #dependent variable
    yerrors::Array{Float64,1}   #uncertainties in dependent variable
	model::Function             #model to be fitted against

    #Fit results
    residuals::Array{Float64,1}  	 #weighted residuals
    params::Array{Float64,1}    		#best fit parameters for model
	marginerror::Array{Float64,1}        #product of standard error and the critical value of each parameter at a certain significance level (default is 5%) from t-distribution.
    perrors::Array{Float64,1}            #uncertainties in best fit parameters
    dof::Int64                  		#number of degrees of freedom in fit
    chisq::Float64             		 #chi^2 value of fit
    probability::Float64      		  #probability of chi^2 value
end

data = input_data(tDTOF, RtDTOF, RtIRF, rand(length(RtDTOF)), 2.56, 1.5, 1.51, 700, false)

modelp = model_params(DT_model, [0.2,10], [0.001,1], [1,20], 0.9, 0.1, 2)

max_stars_repo_name: sammzer/LightPropagation

max_issues_repo_name: sammzer/LightPropagation

max_forks_repo_name: sammzer/LightPropagation