text:
using Distributions
using EEG
using DataFrames
using DSP
using Logging
using Docile
using Compat
using Gadfly

"""
    hotelling(::SSR; kwargs...)

Hotelling test on SSR data

Performs a hotelling test do determine the presenence and properties of a steady state response.
Results are stored as a data frame.

#### Example

```julia
s = SSR(...)
s = extract_epochs(s)
s = hotelling(s)
show(s.processing["statistics"])
```

"""
function hotelling(s::SSR; freq_of_interest::Union{Real, AbstractArray} = modulationrate(s), ID::AbstractString = "",
    spectrum_type::Function = _hotelling_spectrum, data_type::AbstractString="epochs",
    fs::Number=samplingrate(s), results_key::AbstractString="statistics", kwargs...)

    # Calculate spectrum of each epoch
    # Do calculation here once, instead of in each low level call
    spectrum, frequencies = spectrum_type(s.processing[data_type], fs, freq_of_interest)
    spectrum  = compensate_for_filter(s.processing, spectrum, fs)

    for freq in freq_of_interest

        snrDb, phase, signal, noise, statistic = hotelling(spectrum, frequencies, freq)

        actual_freq = frequencies[_find_closest_number_idx(frequencies, freq)]

        result = DataFrame( ID                  = vec(repmat([ID], length(channelnames(s)), 1)),
                            Channel             = copy(channelnames(s)),
                            ModulationRate      = copy(modulationrate(s)),
                            AnalysisType        = "Hotelling",
                            AnalysisFrequency   = actual_freq,
                            SignalAmplitude     = vec(sqrt(signal)),
                            SignalPhase         = vec(phase),
                            NoiseAmplitude      = vec(sqrt(noise)),
                            SNRdB               = vec(snrDb),
                            Statistic           = vec(statistic)  )

        result = add_dataframe_static_rows(result, kwargs)

        if haskey(s.processing, results_key)
            s.processing[results_key] = vcat(s.processing[results_key], result)
        else
            s.processing[results_key] = result
        end

    end

    return s
end

@deprecate hotelling(a::SSR, freq_of_interest::Union{Real, AbstractArray}; kwargs...) hotelling(a; freq_of_interest=freq_of_interest, kwargs...)

function hotelling{T <: AbstractFloat}(spectrum::Array{Complex{T},3}, frequencies::AbstractArray, freq_of_interest::Real)

    Logging.info("Calculating hotelling statistic on $(size(spectrum)[end]) channels at $freq_of_interest Hz with $(size(spectrum)[2]) epochs")

    idx  = _find_closest_number_idx(frequencies, freq_of_interest)
    for_stats = [real(spectrum[idx, :, :]) ; imag(spectrum[idx, :, :])]

    signal_amplitude = Array(T, size(spectrum, 3))
    signal_power     = Array(T, size(spectrum, 3))
    signal_phase     = Array(T, size(spectrum, 3))
    noise_power      = Array(T, size(spectrum, 3))
    snrDb            = Array(T, size(spectrum, 3))
    statistic        = Array(T, size(spectrum, 3))

    for c in 1:size(spectrum,3)

        # Signal metrics
        signal_amplitude[c] = abs(mean(spectrum[idx, :, c], 2))[1]
        signal_power[c]     = signal_amplitude[c]^2
        signal_phase[c]     = angle(mean(spectrum[idx, :, c], 2))[1]

        # Calculate noise
        n = [squeeze(spectrum[idx, :, c], 1); mean(squeeze(spectrum[idx, :, c], 1))]
        n = std(n)                      # Noise amplitude per epoch
        n = n / sqrt(size(spectrum, 2)) # Noise amplitude of recording

        # Save metrics
        noise_power[c] = real(n .^2)    # Noise power of recording
        snrDb[c]       = 10 * log10( signal_power[c] / noise_power[c] )
        statistic[c]   = _hotelling_T2_1sample(for_stats[:, :, c]')
    end

    return snrDb, signal_phase, signal_power, noise_power, statistic
end

function _hotelling_T2_1sample(data; corr::Number=1)

    n, p = size(data)

    m = mean(data, 1)

    if sum(cov(data)) == 0

        statistic = 1.0

    else

        T2 = (corr * n) * (m / cov(data)) * m'
        continuous_distribution = Chisq(p)
        statistic = ccdf(continuous_distribution, T2)[1]

    end

    return statistic
end

function _apes_spectrum{T <: AbstractFloat}(sweep::Array{T,3}, fs::Number, freq_of_interest)

    spectrum = apes(sweep, freq_of_interest / (fs / 2))

    return spectrum, freq_of_interest
end

# Calculate spectrum and return the associated frequencies
# Ignore the freq_of_interest variable as this is used for other spectrum estimation types
function _hotelling_spectrum{T <: AbstractFloat}(sweep::Array{T,3}, fs::Number, freq_of_interest)

    sweepLen = size(sweep)[1]

    p = plan_rfft(sweep, 1, flags = FFTW.MEASURE)

    spectrum = (2 / sweepLen) * (p * sweep)

    frequencies = linspace(0, 1, size(spectrum, 1)) * fs / 2

    return spectrum, frequencies
end

function plot_hotelling{T <: AbstractFloat}(spectrum::Array{Complex{T},3}, frequencies::AbstractArray, freq_of_interest::Real; c::Int=1, fig_name="hot.pdf")

    idx  = _find_closest_number_idx(frequencies, freq_of_interest)
    for_plots = [vec(real(spectrum[idx, :, :])) ; vec(imag(spectrum[idx, :, :]))]

    # Signal metrics
    signal_amplitude = abs(mean(spectrum[idx, :, c], 2))[1]
    signal_power     = signal_amplitude[c]^2
    signal_phase     = angle(mean(spectrum[idx, :, c], 2))[1]

    # Calculate noise
    n = [squeeze(spectrum[idx, :, c], 1), mean(squeeze(spectrum[idx, :, c], 1))]
    n = std(n)                      # Noise amplitude per epoch
    n = n / sqrt(size(spectrum, 2)) # Noise amplitude of recording

    # Save metrics
    noise_power = real(n .^2)    # Noise power of recording
    snrDb       = 10 * log10( signal_power / noise_power )

    d = layer(x=for_plots, Geom.histogram)
    s = layer(xintercept=[signal_amplitude], Geom.vline(color="black"))
    n = layer(xintercept=[-n, n] .+ signal_amplitude, Geom.vline(color="red"))
    p = plot(n, s, d,
        Guide.xlabel("Real and Imag FFT Bin (uV)"),
        Guide.title("SNR = $(round(snrDb[1], 3)) (dB)"))

    draw(PDF(fig_name, 26cm, 17cm), p)

    return p
end


text:
# Try to approximate squared-exponential kernel:
using PyPlot
using Optim

function matern5_2_approx_2()

ntime = 1000
time = collect(linspace(0,10,ntime))
matern52=(1.0+time+time.^2./3.).*exp(-time)
nexp = 3

function compute_chi(x)
model = zeros(ntime)
for k=1:nexp
  model += x[1+(k-1)*2].*exp(-x[2+(k-1)*2].*time)
end
return sum((matern52-model).^2)
end

# If we want to approximate Gaussian, then we just get a cosine... that doesn't work well!

# How about sum of two?  This I solved numerically with Mathematically, and unfortunately I find
# a growing solution in each case... darn.

npar = 2*nexp
abest=zeros(npar)
abest[1:6]=[-9.,1.,10.,0.9,0.1,0.1]
atrial = zeros(npar)
for k=1:npar
  atrial[k] = abest[k]
end
model = zeros(ntime)
for k=1:nexp
  model += atrial[1+(k-1)*2].*exp(-atrial[2+(k-1)*2].*time)
end
clf()
plot(time,matern52)
plot(time,model)
#chibest=sum((matern52-model).^2)
chibest=compute_chi(atrial)

println("Initial chi-square: ",chibest)
read(STDIN,Char)

result = optimize(compute_chi, atrial, BFGS(), OptimizationOptions(autodiff = true))

println(Optim.minimizer(result),Optim.minimum(result))
atrial = Optim.minimizer(result)
clf()
plot(time,(1+time+time.^2/3).*exp(-time))
plot(time,model)
model = zeros(ntime)
for k=1:nexp
  model += atrial[1+(k-1)*2].*exp(-atrial[2+(k-1)*2].*time)
end
chi=sum((matern52-model).^2)
plot(time,(matern52-model))
println("New minimum: ",chi," ",atrial," ",std(matern52-model))
read(STDIN,Char)

atrial = zeros(npar)
for i=1:1000
  for j=1:npar-2
    atrial[j] = abest[j] + 0.01*randn()
    atrial[7]= 1.0-atrial[1]-atrial[3]-atrial[5]
    atrial[8]=-(atrial[1]*atrial[2]+atrial[3]*atrial[4]+atrial[5]*atrial[6])/atrial[7]
#    atrial[5]= 1.0-atrial[1]-atrial[3]
#    atrial[6]=-(atrial[1]*atrial[2]+atrial[3]*atrial[4])/atrial[5]
    model = zeros(ntime)
    for k=1:nexp
      model += atrial[1+(k-1)*2].*exp(-atrial[2+(k-1)*2].*time)
    end
    chi=sum(((1.0+time+time.^2./3.0).*exp(-time)-model).^2)
    if chi < chibest
      chibest = chi
      for k=1:npar
        abest[k] = atrial[k]
      end
      println("New minimum: ",i," ",chibest," ",abest)
      clf()
      plot(time,(1+time+time.^2/3).*exp(-time))
      plot(time,model)
      plot(time,(1+time+time.^2/3).*exp(-time)-model)
    end
  end
end

end

ext: jl

max_stars_repo_name: dfm/ess

max_issues_repo_name: dfm/ess

max_forks_repo_name: dfm/ess


text:
# # Replication of the Quantitative Analysis in Afrouzi and Yang (2020)

# This notebook replicates the qunatitative analysis in [Afrouzi and Yang (2020)](http://www.afrouzi.com/dynamic_inattention.pdf) using the methods and the [solver](https://github.com/afrouzi/DRIPs.jl) from [Afrouzi and Yang (2020)](http://www.afrouzi.com/dynamic_inattention.pdf).

# [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/afrouzi/DRIPs.jl/binder?filepath=examples) to run and modify the following code (no software is needed on the local machine).

# ## Contents
# * [Setup](@ref setup)
# * [A Three-Equation RI Model](@ref three_equation)
# * [Matrix Representation](@ref matrix_representation)
# * [Initialization](@ref param)
# * [Functions](@ref Functions)
# * [Solution](@ref solve_model)
#     * [Post-Volcker Calibration](@ref model_post_volcker)
#     * [Pre-Volcker Calibration](@ref model_pre_volcker)
# * [Simulation](@ref simulation)
# * [Impulse Response Functions](@ref irfs)

# ## [Setup](@id setup)

# ### Households

# Households are fully rational and maximize their life-time utilities:
# ```math
# \begin{aligned}
# \max & \ \ \mathbb{E}_{t}^{f}\left[\sum_{t=0}^{\infty}\beta^{t}\left(\frac{C_{t}^{1-\sigma}}{1-\sigma}-\frac{\int_{0}^{1}L_{i,t}^{1+\psi}di}{1+\psi}\right)\right]\\
# \text{s.t.} & \ \ \int P_{i,t}C_{i,t} di+B_{t}\leq R_{t-1}B_{t-1}+\int_{0}^{1}W_{i,t}L_{i,t}di+\Pi_{t}, \quad \text{for all } t
# \end{aligned}
# ```
# where
# ```math
# \begin{aligned}
# C_{t} = \left(\int C_{i,t} ^{\frac{\theta-1}{\theta}}di\right)^{\frac{\theta}{\theta-1}}.
# \end{aligned}
# ```
# Here $\mathbb{E}_{t}^{f}\left[\cdot\right]$ is the full information rational expectation operator at time $t$. Since the main purpose of this paper is to study the effects of nominal rigidity and rational inattention among firms, I assume that the household is fully informed about all prices and wages. $B_{t}$ is the demand for nominal bond and $R_{t-1}$ is the nominal interest rate. $L_{i,t}$ is firm-specific labor supply of the household, $W_{i,t}$ is the firm-specific nominal wage, and $\Pi_{t}$ is the aggregate profit from the firms. $C_{t}$ is the aggregator over the consumption for goods produced by firms. $\theta$ is the constant elasticity of substitution across different firms.
#
#
# ### Firms
# There is a measure one of firms, indexed by $i$, that operate in monopolistically competitive markets. Firms take wages and demands for their goods as given, and choose their prices $P_{i,t}$ based on their information set, $S_{i}^{t}$, at that time.  After setting their prices, firms hire labor from a competitive labor market and produce the realized level of demand that their prices induce with a production function,
# ```math
# \begin{aligned}
#     Y_{i,t} = A_{t} L_{i,t},
# \end{aligned}
# ```
# where $L_{i,t}$ is firm $i$'s demand for labor. I assume that shocks to $A_{t}$ are independently and identically distributed and the log of the productivity shock, $a_{i,t}\equiv\log(A_{t})$, follows a AR(1) process:
# ```math
# \begin{aligned}
#     a_{t} = \rho_a a_{t-1} + \varepsilon_{a,t}, \ \   \varepsilon_{a,t} \sim N(0,\sigma_{a}^2).
# \end{aligned}
# ```
# Then, firm $i$'s nominal profit from sales of all goods at prices $\{P_{i,j,t}\}_{j=1}^N$ is given by
# ```math
# \begin{aligned}
# \Pi_{i,t} (P_{i,t},A_{t},W_{i,t},P_t,Y_t) =\left(P_{i,t}-W_{i,t}A_{t}\right) \left(\frac{P_{i,t}}{P_{t}}\right)^{-\theta}Y_{t},
# \end{aligned}
# ```
# where $Y_t$ is the nominal aggregate demand.
#
# At each period, firms optimally decide their prices and signals subject to costs of processing information. Firms are rationally inattentive in a sense that they choose their optimal information set by taking into account the cost of obtaining and processing information. At the beginning of period $t$, firm $i$ wakes up with its initial information set, $S_{i}^{t-1}$. Then it chooses optimal signals, $s_{i,t}$, from a set of available signals, $\mathcal{S}_{i,t}$, subject to the cost of information which is linear in Shannon's mutual information function. Denote $ \omega $ as the marginal cost of information processing. Firm $i$ forms a new information set, $S_i^t = S_i^{t-1} \cup s_{i,t}$, and sets its new prices, $P_{i,t}$, based on that.
#
# The firm $i$ chooses a set of signals to observe over time $(s_{i,t} \in  \mathcal{S}_{i,t})_{t=0}^\infty$ and a pricing strategy that maps the set of its prices at $t-1$ and its information set at $t$ to its optimal price at any given period, $P_{i,t}:(S_i^t)\rightarrow \mathbb{R}$ where $S_i^t=S_i^{t-1}\cup s_{i,t} = S_i^{-1} \cup \{s_{i,\tau}\}_{\tau=0}^t$ is the firm's information set at time $t$.  Then, the firm $i$'s problem is to maximize the net present value of its life time profits given an initial information set:
# ```math
# \begin{aligned}
# \max_{\{s_{i,t}\in \mathcal{S}_{i,t},P_{i,t} (S_i^t) \}_{t\geq0}} &  \mathbb{E}\left[\sum_{t=0}^{\infty}\beta^t \Lambda_t \left\{ \Pi_{i,t} (P_{i,t},A_{t},W_{i,t},P_t,Y_t)  - \omega \mathbb{I} ( S_{i}^t ;  (A_{\tau},W_{i,\tau},P_\tau,Y_\tau)_{\tau \leq t} | S_i^{t-1} ) \right\} \Bigg| S_{i}^{-1}\right]  \\
# \text{s.t.} \quad & {S}_{i}^{t}={S}_{i}^{t-1}\cup s_{i,t}
# \end{aligned}
# ```
# where $\Lambda_t$ is the stochastic discount factor and $\mathbb{I} ( S_{i}^t ;  (A_{\tau},W_{i,\tau},P_\tau,Y_\tau)_{\tau \leq t} | S_i^{t-1} )$ is the Shannon's mutual information function.
#
#
# ### Monetary Policy
# Monetary policy is specified as a standard Talor rule:
# ```math
# \begin{aligned}
#     R_t = (R_{t-1})^{\rho} \left( \Pi_t^{\phi_\pi} \left( \frac{Y_t}{Y_t^n} \right)^{\phi_x} \left(\frac{Y_t}{Y_{t-1}} \right)^{\phi_{\Delta y}} \right)^{1-\rho} \exp(u_t)
# \end{aligned}
# ```
# where $u_t \sim N(0,\sigma_{u}^2)$ is the monetary policy shock.
#
# ## [A Three-Equation GE Rational Inattention Model](@id three_equation)
#
# Our general equlibrium model is characterized by the following three equations with two stochastic processes of technology ($a_{t}$) and monetary policy shocks ($u_t$):
# ```math
# \begin{aligned}
# x_{t} & =\mathbb{E}_{t}^{f}\left[x_{t+1}-\frac{1}{\sigma}\left(i_{t}-\pi_{t+1}\right)\right]+\mathbb{E}_{t}^{f}\left[y_{t+1}^{n}\right]-y_{t}^{n} \\
# p_{i,t} & =\mathbb{E}_{i,t}\left[p_{t} + \alpha x_{t}\right] \\
# i_{t} & =\rho i_{t-1} + \left( 1 - \rho \right) \left(\phi_{\pi} \pi_{t}+\phi_{x} x_{t} - \phi_{\Delta y} \Delta y_t \right) + u_{t}
# \end{aligned}
# ```
# where $\mathbb{E}_{i,t}[\cdot]$ is the firm $i$'s expectation operator conditional on her time $t$ information set, $x_t = y_t - y_t^n$ is the output gap, $y_{t}^{n}=\frac{1+\psi}{\sigma+\psi}a_{t}$ is the natural level of output, $i_t$ is the nominal interest rate, and $\alpha = \frac{\sigma+\psi}{1+\psi\theta}$ is the degree of strategic complementarity.
#
#
# ## [Matrix Representation](@id matrix_representation)
#
# Firms wants to keep track of their ideal price, $p_{i,t}^* = p_t + \alpha x_t$. Notice that the state space representation for $p_{i,t}^*$ is no longer exogenous and is determined in the equilibrium. However, we know that this is a Guassian process and by Wold's theorem we can decompose it to its $MA(\infty)$ representation:
# ```math
# \begin{aligned}
#     p_{i,t}^*=\Phi_a(L)\varepsilon_{a,t} + \Phi_u(L)\varepsilon_{u,t}
# \end{aligned}
# ```
# where $\Phi_a(.)$ and $\Phi_u(.)$ are lag polynomials. Here, we have basically guessed that the process for $p_{i,t}^*$ is determined uniquely by the history of monetary shocks which requires that rational inattention errors of firms are orthogonal (See [Afrouzi (2020)](www.afrouzi.com/strategic_inattetion.pdf)).
#
# Since we cannot put $MA(\infty)$ processes in the computer and have to truncate them. However, we know that for stationary processes we can arbitrarily get close to the true process by truncating $MA(\infty)$ processes. Our problem here is that $p_{i,t}^*$ has a unit root and is not stationary. We can bypass this issue by re-writing the state space in the following way:
# ```math
# \begin{aligned}
#     p_{i,t}^*=\Phi_a(L)\varepsilon_{a,t} + \phi_u(L)\tilde{\varepsilon}_{u,t},\quad \tilde{\varepsilon}_{u,t}=(1-L)^{-1}\varepsilon_{u,t} =\sum_{j=0}^\infty \varepsilon_{u,t-j}
# \end{aligned}
# ```
# here $\tilde{\varepsilon}_{u,t}$ is the unit root of the process and basically we have differenced out the unit root from the lag polynomial, and $\phi_u(L)=(1-L)\Phi_u(L)$. Notice that since the original process was difference stationary, differencing out the unit root means that $\phi_u(L)$ is now in $\ell_2$, and the process can now be approximated arbitrarily precisely with truncation.
#
# For ease of notation, let $z_t = (\varepsilon_{a,t}, \varepsilon_{u,t})$ and $\tilde{z}_t = (\varepsilon_{a,t}, \tilde{\varepsilon}_{u,t})$. For a length of truncation $L$, let $\vec{x}_t' \equiv ({z}_t,{z}_{t-1},\dots,{z}_{t-(L+1)})\in\mathbb{R}^{2L}$ and $\vec{\mathbf{x}}_t' \equiv (\tilde{z}_t,\tilde{z}_{t-1},\dots,\tilde{z}_{t-(L+1)})\in\mathbb{R}^{2L}$. Notice that
# ```math
# \begin{aligned}
#     \vec{x}_t & = (\mathbf{I} - \mathbf{\Lambda} \mathbf{M}')\vec{\mathbf{x}}_t \\
#     \vec{\mathbf{x}}_t & = (\mathbf{I} - \mathbf{\Lambda} \mathbf{M}')^{-1}\vec{x}_t
# \end{aligned}
# ```
# where $\mathbf{I}$ is a $2 \times 2$ identity matrix, $\mathbf{\Lambda}$ is a diagonal matrix where $\mathbf{\Lambda}_{(2i,2i)} = 1$ and $\mathbf{\Lambda}_{(2i-1,2i-1)}=0$ for all $i={1,2,\cdots,L}$, and $\mathbf{M}$ is a shift matrix:
# ```math
# \begin{aligned}
#    \mathbf{M} = \left[\begin{array}{ccccc}
#         0 & 0 & \cdots & 0 & 0 & 0 & 0\\
#         0 & 0 & \cdots & 0 & 0 & 0 & 0\\
#         1 & 0 & \cdots & 0 & 0 & 0 & 0\\
#         0 & 1 & \cdots & 0 & 0 & 0 & 0\\
#         \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots\\
#         0 & 0 & \cdots & 1 & 0 & 0 & 0\\
#         0 & 0 & \cdots & 0 & 1 & 0 & 0
#         \end{array}\right]
# \end{aligned}
# ```
# Then, note that $p_{i,t}^*\approx \mathbf{H} '\vec{\mathbf{x}}_{t}$ where $\mathbf{H}\in \mathbb{R}^{2L}$ is the truncated matrix analog of the lag polynominal, and is endogenous to the problem. Our objective is to find the general equilibrium $\mathbf{H}$ along with the optimal information structure that it implies.
#
# Moreover, note that
# ```math
# \begin{aligned}
# a_t & = \mathbf{H}_a'\vec{x}_t,\quad \mathbf{H}_a'=(1,0,\rho_a,0,\rho_a^2,0,\dots,\rho_a^{L-1},0) \\
# u_t & = \mathbf{H}_u'\vec{x}_t,\quad \mathbf{H}_u'=(0,1,0,0,0,0,\dots,0,0)
# \end{aligned}
# ```
# We will solve for $\mathbf{H}$ by iterating over the problem. In particular, in iteration $n\geq 1$, given the guess $\mathbf{H}_{(n-1)}$, we have the following state space representation for the firm's problem
# ```math
# \begin{aligned}
#    \vec{\mathbf{x}}_{t}&
#     = \underset{\mathbf{A}}{\underbrace{\left[\begin{array}{ccccc}
#         0 & 0 & 0 & \cdots & 0 & 0 & 0 & 0\\
#         0 & 1 & 0 & \cdots & 0 & 0 & 0 & 0\\
#         1 & 0 & 0 & \cdots & 0 & 0 & 0 & 0\\
#         0 & 1 & 0 & \cdots & 0 & 0 & 0 & 0\\
#         0 & 0 & 1 & \cdots & 0 & 0 & 0 & 0\\
#         \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots\\
#         0 & 0 & 0 & \cdots & 1 & 0 & 0 & 0\\
#         0 & 0 & 0 & \cdots & 0 & 1 & 0 & 0
#         \end{array}\right]}}\, \vec{\mathbf{x}}_{t-1}
#     + \underset{\mathbf{Q}}{\underbrace{\left[\begin{array}{cc}
#         1 & 0\\
#         0 & 1\\
#         0 & 0\\
#         \vdots & \vdots \\
#         0 & 0
#         \end{array}\right]}}\, z_t, \\
#         p_{i,t}^*&=\mathbf{H}_{(n-1)}'\vec{\mathbf{x}}_{t}
# \end{aligned}
# ```
# Now, note that
# ```math
# \begin{aligned}
#     p_{t} &= \int_{0}^{1} p_{i,t}di = \mathbf{H}_{(n-1)}' \int_{0}^{1}\mathbb{E}_{i,t}[\vec{\mathbf{x}}_{t}]di \\
#                   &= \mathbf{H}_{(n-1)}' \sum_{j=0}^{\infty}[(\mathbf{I}-\mathbf{K}_{(n)}\mathbf{Y}'_{(n)})\mathbf{A}]^{j}\mathbf{K}_{(n)}\mathbf{Y}'_{(n)}\vec{\mathbf{x}}_{t-j} \\
#                   &\approx\underset{\equiv\mathbf{X}_{(n)}}{\mathbf{H}_{(n-1)}' \underbrace{\left[\sum_{j=0}^{\infty}[(\mathbf{I}-\mathbf{K}_{(n)}\mathbf{Y}'_{(n)})\mathbf{A}]^{j}\mathbf{K}_{(n)}\mathbf{Y}'_{(n)}\mathbf{M}'^{j}\right]}}\vec{\mathbf{x}}_{t} \\
#                   &= \mathbf{H}_{(n-1)}' \mathbf{X}_{(n)}\vec{\mathbf{x}}_{t} = \mathbf{H}_{p}' \vec{\mathbf{x}}_{t}
# \end{aligned}
# ```
#
#
#
# Let $x_t = \mathbf{H}_{x}'\vec{x}_{t} $, $i_t = \mathbf{H}_{i}'\vec{x}_{t} $, and $\pi_t = \mathbf{H}_{\pi}'\vec{x}_{t} = \mathbf{H}_{p}'(\mathbf{I} - \mathbf{\Lambda} \mathbf{M}')^{-1} (\mathbf{I} - \mathbf{M}') \vec{x}_t $. Then from the households Euler equation, we have:
# ```math
# \begin{aligned}
#     & x_{t} =\mathbb{E}_{t}^{f}\left[x_{t+1}-\frac{1}{\sigma}\left(i_{t}-\pi_{t+1}\right)\right]+\mathbb{E}_{t}^{f}[y_{t+1}^n] - y_t^n \\
#     \Longrightarrow \ \ & \mathbf{H}_{i} =\sigma\left(\mathbf{M}'-\mathbf{I}\right)\mathbf{H}_{x} + \frac{\sigma(1+\psi)}{\sigma+\psi}\left(\mathbf{M}'-\mathbf{I}\right) \mathbf{H}_{a}+\mathbf{M}'\mathbf{H}_{\pi}
# \end{aligned}
# ```
# Also, the Talyor rule gives:
# ```math
# \begin{aligned}
# & i_{t} =\rho_{1}i_{t-1}+\rho_{2}i_{t-2}+\left(1-\rho_{1}-\rho_{2}\right)\left(\phi_{\pi}\pi_{t}+\phi_{x}x_{t}+\phi_{\Delta y}\left(y_{t}-y_{t-1}\right)\right)+u_{t} \\
# \Longrightarrow \ \ & \left(\mathbf{I}-\rho_{1}\mathbf{M}-\rho_{2}\mathbf{M}^{2}\right)\mathbf{H}_{i} =\left(1-\rho_{1}-\rho_{2}\right)\phi_{\pi}\mathbf{H}_{\pi}+\left(1-\rho_{1}-\rho_{2}\right)\phi_{x}\mathbf{H}_{x}\\
#  & +\left(1-\rho_{1}-\rho_{2}\right)\phi_{\Delta y}\left(\mathbf{I}-\mathbf{M}\right)\left(\mathbf{H}_{x} + \frac{1+\psi}{\sigma+\psi} \mathbf{H}_{a} \right) + \mathbf{H}_{u}
# \end{aligned}
# ```
#
# These give us $\mathbf{H}_{x}$ and $\mathbf{H}_{i}$ and we update new $\mathbf{H}_{(n)}$ using:
# ```math
# \begin{aligned}
#     & \mathbf{H}_{(n)} = \mathbf{H}_{p} + \alpha (\mathbf{I} - \mathbf{M} \mathbf{\Lambda}') \mathbf{H}_{x}
# \end{aligned}
# ```
# We iterate until convergence of $\mathbf{H}_{(n)}$.
#
# ## [Initialization](@id param)
# Load DRIPs solver and other packages:

using DRIPs;
using BenchmarkTools, LinearAlgebra, GLM, Statistics, Suppressor, Printf;
using PyPlot; rc("text", usetex="True") ;
              rc("font",family="serif",serif=:"Palatino") ;
using Plots, LaTeXStrings; pyplot() ;

# Set parameters:
struct param
    β; σ; ψ; θ; α; # Deep parameters
    ϕ_π; ϕ_x; ϕ_dy; ρ; # Monetary policy parameters
    ρa; ρu; σa; σu; # Shock paramters
end

# ###Assign model parameters
σ       = 2.5       ; #Risk aversion
β       = 0.99      ; #Time discount
ψ       = 2.5       ; #Inverse of Frisch elasticity of labor supply
θ       = 10        ; #Elasticity of substitution across firms
α       = (σ+ψ)/(1+ψ*θ) ; #Strategic complementarity (1-α)

# ###Assign monetary policy parameters: post-Volcker
ϕ_π     = 2.028     ; #Taylor rule response to inflation
ϕ_x     = 0.673/4   ; #Taylor rule response to output
ϕ_dy    = 3.122     ; #Taylor rule response to growth
ρ       = 0.9457    ; #interest rate smoothing

# ###Assign monetary policy parameters: pre-Volcker
ϕ_π_pre = 1.589     ; #Taylor rule response to inflation
ϕ_x_pre = 1.167/4   ; #Taylor rule response to output
ϕ_dy_pre= 1.028     ; #Taylor rule response to output growth
ρ_pre   = 0.9181    ; #interest rate smoothing

# ###Assign parameters governing shock processes
ρu      = 0.0       ; #persistence of MP shock (post-Volcker)
σu      = 0.279     ; #S.D. of MP shock (post-Volcker)
ρu_pre  = 0.0       ; #persistence of MP shock (pre-Volcker)
σu_pre  = 0.535     ; #S.D. of MP shock (pre-Volcker)

# ###Calibrated parameters
ρa      = 0.85      ; #persistence of technology shock
σa      = 1.56      ; #S.D. of technology shock
ω       = 0.773     ; #marginal cost of information

# ###Parameters for simulation/irfs
simT    = 50000     ;
nburn   = 500       ;
T       = 20        ;

# Primitives of Drip

numshock= 2         ; #number of shocks
L       = 160       ; #length of trunction
M       = [zeros(1,L-1) 0; Matrix(I,L-1,L-1) zeros(L-1,1)];
M       = M^2       ;
J       = zeros(L,L); J[2,2] = 1 ;
Lambda  = zeros(L,L);
for i = 1:Int64(L/numshock); Lambda[i*numshock,i*numshock] = 1; end
A       = M + J     ;
eye     = Matrix(I,L,L);

# ## [Functions](@id Functions)
# We start with a function that solves the GE problem and returns the solution in a `Drip` structure:

function agg_drip(p::param,ω,A,Q,M,Lambda,
                  Ha,               #state space rep. of a
                  Hu;               #state space rep. of u
                  H0    = Hu+Ha,    #optional: initial guess for H0
                  Sigma = A*A'+Q*Q',#optional: initial guess for Σ_0
                  Omega = H0*H0',   #optional: initial guess for Ω
                  maxit = 10000,    #optional: max. iterations for GE code
                  maxit_in = 100,   #optional: max. iterations for solving DRIP
                  tol   = 1e-4,     #optional: tolerance for iterations
                  w     = 1)        #optional: update weight for RI
    err   = 1;
    iter  = 1;
    L     = length(H0);
    eye   = Matrix(I,L,L);

    temp0 = (eye-M)*inv(eye - M*Lambda') ;
    Htemp1 = (p.σ*(eye-p.ρ*M)*(M'-eye) - (1-p.ρ)*(p.ϕ_x*eye + p.ϕ_dy*(eye-M))) ;
    Htemp2 = (1-p.ρ)*p.ϕ_π*eye - (eye-p.ρ*M)*M' ;
    Htemp3 = (1+p.ψ)/(p.σ+p.ψ)*((1-p.ρ)*p.ϕ_dy*(eye-M) - p.σ*(eye-p.ρ*M)*(M'-eye)) ;

    while (err > tol) & (iter < maxit)

        if iter == 1
            global ge  = Drip(ω,p.β,A,Q,H0;
                              Ω0=Omega, Σ0=Sigma, w=w, maxit=maxit_in);
        else
            global ge  = Drip(ω,p.β,A,Q,H0;
                              Ω0=ge.ss.Ω, Σ0=ge.ss.Σ_1, w=w, maxit=maxit_in);
        end

        XFUN(jj) = ((eye-ge.ss.K*ge.ss.Y')*ge.A)^jj * (ge.ss.K*ge.ss.Y') * (M')^jj
        X = DRIPs.infinitesum(XFUN; maxit=200, start = 0);  #E[x⃗]=X×x⃗

        global Hp = X'*H0 ;
        global Hπ = temp0*Hp ;

        Hπ[L-20:end,:] .= 0

        global Hx = (Htemp1)\(Htemp2*Hπ + Htemp3*Ha + Hu) ;

        global H1 = Hp + p.α*(eye - M*Lambda')*Hx ;

        err= 0.5*norm(H1-H0,2)/norm(H0)+0.5*err;

        H0 = H1;

        cap = DRIPs.capacity(ge, unit = "bit")

        if iter == maxit
            print("***GE loop hit maxit -- no convergence\n")
        elseif mod(iter,50) == 0
            println("  Iteration $iter. Error: $err. Capacity: $cap.")
        end

        iter += 1;
    end

    return(ge,H1,Hx,Hp,Hπ)
end;

# ## [Model Solution](@id solve_model)
# ### [Post-Volcker Calibration](@id model_post_volcker)
# Start with post-Volcker calibration:

Q       = zeros(L,2); Q[1,1]=σa; Q[2,2]=σu;

Ha_post = ρa.^(0:1:L/2-1) ;
Hu_post = ρu.^(0:1:L/2-1) ;
Ha_post = kron(Ha_post,[1,0])[:,:] ;
Hu_post = kron(Hu_post,[0,1])[:,:] ;

p       = param(β,σ,ψ,θ,α,ϕ_π,ϕ_x,ϕ_dy,ρ,ρa,ρu,σa,σu) ;

# Get initial guess
print("\nGet initial guess for post-Volcker solutions\n")
@time begin
    @suppress agg_drip(p,ω,A,Q,M,Lambda,Ha_post,Hu_post;
                       w=0.95, maxit=300, maxit_in=5);
end;

# Solve the model: Post-Volcker
print("\nSolve for the post-Volcker model:\n");
@time (ge_post,H1_post,Hx_post,Hp_post,Hπ_post) =
            agg_drip(p,ω,A,Q,M,Lambda,Ha_post,Hu_post;
                     H0    = H1,
                     Sigma = ge.ss.Σ_1,
                     Omega = ge.ss.Ω,
                     w=0.95, maxit=5000, maxit_in=500) ;

Hy_post = Hx_post + (1+ψ)/(σ+ψ)*Ha_post     ;
Hi_post = M'*Hπ_post + σ*(M'-eye)*Hy_post   ;
Hr_post = Hi_post - (M')*Hπ_post            ;

# ### [Pre-Volcker Calibration](@id model_pre_volcker)
# Then, pre-Volcker calibration:

Q_pre    = zeros(L,2); Q_pre[1,1]=σa; Q_pre[2,2]=σu_pre;

Ha_pre   = ρa.^(0:1:L/2-1)          ;
Hu_pre   = ρu_pre.^(0:1:L/2-1)      ;
Ha_pre   = kron(Ha_pre,[1,0])[:,:]  ;
Hu_pre   = kron(Hu_pre,[0,1])[:,:]  ;

p_pre    = param(β,σ,ψ,θ,α,ϕ_π_pre,ϕ_x_pre,ϕ_dy_pre,
                 ρ_pre,ρa,ρu_pre,σa,σu_pre) ;

# Get initial guess: Pre-Volcker
print("\nGet initial guess for pre-Volcker solutions\n")
@time begin
    @suppress agg_drip(p_pre,ω,A,Q_pre,M,Lambda,Ha_pre,Hu_pre;
                       H0    = H1_post,
                       Sigma = ge_post.ss.Σ_1,
                       Omega = ge_post.ss.Ω,
                       w=0.95, maxit=300, maxit_in=5) ;
end;

# Solve the model: Pre-Volcker
print("\nSolve for the pre-Volcker model:\n");
@time begin (ge_pre,H1_pre,Hx_pre,Hp_pre,Hπ_pre) =
            agg_drip(p_pre,ω,A,Q_pre,M,Lambda,Ha_pre,Hu_pre;
                     H0    = H1,
                     Sigma = ge.ss.Σ_1,
                     Omega = ge.ss.Ω,
                     w=0.95, maxit=5000, maxit_in=500) ;
end

Hy_pre    = Hx_pre + (1+ψ)/(σ+ψ)*Ha_pre     ;
Hi_pre    = M'*Hπ_pre + σ*(M'-eye)*Hy_pre   ;
Hr_pre    = Hi_pre - (M')*Hπ_pre            ;

# ## [Model Simulation](@id simulation)
# Start with the simulation for post-Volcker:

print("\nSimulate the models:\n");

@time begin
    sim_post    = simulate(ge_post; T=simT, burn=nburn, seed=1) ;

    x_shock     = sim_post.x        ;
    xhat_avg    = sim_post.x_hat    ;

    sim_π       = (Hπ_post'*(eye-Lambda*M')*x_shock)' ;
    sim_y       = (Hy_post'*(eye-Lambda*M')*x_shock)' ;
    sim_x       = (Hx_post'*(eye-Lambda*M')*x_shock)' ;

    mat_sim     = [sim_π sim_y sim_x]   ;
    cor_sim     = cor(mat_sim)          ;

    stat_post   = vec([std(sim_π/100) std(sim_y/100) cor_sim[2,1]]) ;

    s = @sprintf("==> Post-Volcker: std(π)=%5.3f, std(y)=%5.3f, corr(π,y)=%5.3f",
                 stat_post[1], stat_post[2], stat_post[3])  ;
    println(s) ;
end

# Now simulate for pre-Volcker:

@time begin
    sim_pre     = simulate(ge_pre; T=simT, burn=nburn, seed=1) ;

    x_shock_pre = sim_pre.x     ;
    xhat_avg_pre= sim_pre.x_hat ;

    sim_π_pre   = (Hπ_pre'*(eye-Lambda*M')*x_shock_pre)';
    sim_y_pre   = (Hy_pre'*(eye-Lambda*M')*x_shock_pre)';
    sim_x_pre   = (Hx_pre'*(eye-Lambda*M')*x_shock_pre)';

    mat_sim_pre = [sim_π_pre sim_y_pre sim_x_pre]   ;
    cor_sim_pre = cor(mat_sim_pre)                  ;

    stat_pre    = vec([std(sim_π_pre/100) std(sim_y_pre/100) cor_sim_pre[2,1]]);

    s = @sprintf("==> Pre-Volcker : std(π)=%5.3f, std(y)=%5.3f, corr(π,y)=%5.3f",
                 stat_pre[1],stat_pre[2],stat_pre[3])  ;
    println(s) ;
end

# ## [Impulse Response Functions](@id irfs)
# Reshape variables for IRFs

pi      = reshape(Hπ_post,numshock,Int64(L/numshock))' ;
x       = reshape(Hx_post,numshock,Int64(L/numshock))' ;
y       = reshape(Hy_post,numshock,Int64(L/numshock))' ;
i       = reshape(Hi_post,numshock,Int64(L/numshock))' ;
r       = reshape(Hr_post,numshock,Int64(L/numshock))' ;
a       = reshape(Ha_post,numshock,Int64(L/numshock))' ;
u       = reshape(Hu_post,numshock,Int64(L/numshock))' ;

pi_pre  = reshape(Hπ_pre,numshock,Int64(L/numshock))' ;
x_pre   = reshape(Hx_pre,numshock,Int64(L/numshock))' ;
y_pre   = reshape(Hy_pre,numshock,Int64(L/numshock))' ;
i_pre   = reshape(Hi_pre,numshock,Int64(L/numshock))' ;
r_pre   = reshape(Hr_pre,numshock,Int64(L/numshock))' ;
a_pre   = reshape(Ha_pre,numshock,Int64(L/numshock))' ;
u_pre   = reshape(Hu_pre,numshock,Int64(L/numshock))' ;

# Plot IRFs:

# #Draw IRFs
title1 = Plots.plot(ylabel     = "IRFs to Technoloygy Shock",
             guidefont  = font(10),
             grid       = false, showaxis = false,
             bottom_margin = 20Plots.px)
title2 = Plots.plot(ylabel     = "IRFs to Monetary Shock",
             guidefont  = font(10),
             grid       = false, showaxis = false,
             bottom_margin = -50Plots.px,
             top_margin = 30Plots.px)

p1 = Plots.plot(0:T,[σa*pi[1:T+1,1],σa*pi_pre[1:T+1,1]],
        title           = L"Inflation ($\pi_t$)",
        ylabel          = "Percent",
        guidefont       = font(9),
        yticks          = -1:0.2:0.3,
        label           = ["Post-Volcker" "Pre-Volcker"],
        legend          = :bottomright,
        legendfont      = font(8),
        color           = [:black :darkgray],
        linestyle       = [:solid :dashdot],
        )
p2 = Plots.plot(0:T,[σa*y[1:T+1,1],σa*y_pre[1:T+1,1]],
        title           = L"Output ($y_t$)",
        legend          = false,
        color           = [:black :darkgray],
        linestyle       = [:solid :dashdot],
        )
p3 = Plots.plot(0:T,[σa*i[1:T+1,1],σa*i_pre[1:T+1,1]],
        title           = L"Nominal Rate ($i_t$)",
        legend          = false,
        yticks          = -0.15:0.03:0.0,
        color           = [:black :darkgray],
        linestyle       = [:solid :dashdot],
        )
p4 = Plots.plot(0:T,[σa*r[1:T+1,1],σa*r_pre[1:T+1,1]],
        title           = L"Real Rate ($r_t$)",
        legend          = false,
        color           = [:black :darkgray],
        linestyle       = [:solid :dashdot],
        )
p5 = Plots.plot(0:T,[-σu*pi[1:T+1,2],-σu_pre*pi_pre[1:T+1,2]],
        title           = L"Inflation ($\pi_t$)",
        ylabel          = "Percent",
        xlabel          = "Time",
        guidefont       = font(9),
        label           = ["Post-Volcker" "Pre-Volcker"],
        legend          = :topright,
        legendfont      = font(8),
        color           = [:black :darkgray],
        linestyle       = [:solid :dashdot],
        )
p6 = Plots.plot(0:T,[-σu*y[1:T+1,2],-σu_pre*y_pre[1:T+1,2]],
        title           = L"Output ($y_t$)",
        xlabel          = "Time",
        guidefont       = font(9),
        legend          = false,
        color           = [:black :darkgray],
        linestyle       = [:solid :dashdot],
        )
p7 = Plots.plot(0:T,[-σu*i[1:T+1,2],-σu_pre*i_pre[1:T+1,2]],
        title           = L"Nominal Rate ($i_t$)",
        xlabel          = "Time",
        guidefont       = font(9),
        legend          = false,
        color           = [:black :darkgray],
        linestyle       = [:solid :dashdot],
        )
p8 = Plots.plot(0:T,[-σu*r[1:T+1,2],-σu_pre*r_pre[1:T+1,2]],
        title           = L"Real Rate ($r_t$)",
        xlabel          = "Time",
        guidefont       = font(9),
        legend          = false,
        color           = [:black :darkgray],
        linestyle       = [:solid :dashdot],
        )

l = @layout [
            a{0.001w} Plots.grid(1,4)
            a{0.001w} Plots.grid(1,4)
            ]

Plots.plot(title1,p1,p2,p3,p4,title2,p5,p6,p7,p8,
        layout      = l,
        gridstyle   = :dot,
        gridalpha   = 0.2,
        lw          = [2.5 2],
        titlefont   = font(10),
        xticks      = (0:4:T),
        xlim        = (0,T),
        tickfont    = font(8),
        size        = (900,550),
        framestyle  = :box)


text:
module test_jive_skip_exprs

using Jive # Skipped @skip
using Test # @test

empty!(Jive.Skipped.expressions)

@skip if 3 > 2
end

@skip quote
end

@skip begin
end

@skip macro m()
end

@skip using Jive, Test

@skip return

@skip f() do x
end

@skip struct A
end

@skip mutable struct B
end

@skip :hello
@skip true
@skip pi
@skip 42
@skip nothing
@skip a = 2
@skip [1,2,3]

@skip for i in 1:2
end

@test Jive.Skipped.expressions == [:if, :quote, :block, :macro=>:m, :using, :return, :do, :struct=>:A, :struct=>:B, :hello, true, :pi, 42, nothing, :(=), :vect, :for]
@test !isdefined(@__MODULE__, :a)

empty!(Jive.Skipped.expressions)
ENV["JIVE_SKIP"] = "0"

@skip if 10 > 2
end

@skip :hello
@skip true
@skip pi
@skip 42
@skip nothing
@skip a = 2
@skip [1,2,3]

@test isempty(Jive.Skipped.expressions)
@test isdefined(@__MODULE__, :a)

ENV["JIVE_SKIP"] = "1"

end # module test_jive_skip_exprs


text:
using HackMoji
using Test
using Suppressor

@testset "HackMoji.jl" begin
    # Write your tests here.
    out = @capture_out runme("++++++++[>++++[>++>+++>+++>+<<<<-]>+>+>->>+[<]<-]>>.>---.+++++++..+++.>>.<-.<.+++.------.--------.>>+.>++.")
    @test out == "Hello World!\n"

    out = @capture_out runme("++++++++++++++++++++++++++++++++++++++++++++++++.")
    @test out == "0"

    twoplusfive = """
    ++       Cell c0 = 2
    > +++++  Cell c1 = 5
    [        Start your loops with your cell pointer on the loop counter (c1 in our case)
    < +      Add 1 to c0
    > -      Subtract 1 from c1
    ]        End your loops with the cell pointer on the loop counter
    At this point our program has added 5 to 2 leaving 7 in c0 and 0 in c1
    but we cannot output this value to the terminal since it is not ASCII encoded!
    To display the ASCII character "7" we must add 48 to the value 7
    48 = 6 * 8 so let's use another loop to help us!
    ++++ ++++  c1 = 8 and this will be our loop counter again
    [
    < +++ +++  Add 6 to c0
    > -        Subtract 1 from c1
    ]
    < .        Print out c0 which has the value 55 which translates to "7"!
    """
    out = @capture_out runme(twoplusfive)
    @test out == "7"

end


text:
## basic modifier expressors

"""
	A(v::Float64, z)

Accelirate `z`::Expi; multiplies all atom velocities by `v`.
"""
A(v::Float64, z) = 0 < v ? accel(v,atomOrExp(z)) : throw(DomainError())

"""
	D(v::IR, z)

Dilate `z`; multiplies all atom durations and offsets by `v`.
"""
D(v::IR, z) = 0 < v ? dilate(Rat(v),atomOrExp(z)) : throw(DomainError())

"""
	F(v::IR, z)

Shift `z`: adds `v` to all atom offsets
"""
F(v::IR, z) = sshift(Rat(v),atomOrExp(z))

"""
	T(v::Int, z)

Translate `z`; adds `v` to all note interval values.
"""
T(v::Int, z) = transl(v,atomOrExp(z))

"""
	V(v::Int, z)

Octave `z`; sets all note octave values to `v`.
"""
V(v::Int, z) = octave(v,atomOrExp(z))

"""
	Cha(v::Int, z)

Channel `z`; sets all channel values to `v`.
"""
Cha(v::Int, z) = channel(v,atomOrExp(z))

"""
	Sca(v::Function, z)

Scale `z`; sets all note scales to `v`.  See SirenSeq.Scales
"""
Sca(v::Scales.Scale, z) = sscale(v,atomOrExp(z))

## used by constructor N
defaultScale = Scales.cMaj
defaultOctave = 3

"""
	setDefaultScale(sca::Function)

Sets the default scale for `N` constructor (defualt value is cMaj).
"""
setDefaultScale(sca::Scales.Scale) = ( global defaultScale = sca ; )

"""
	setDefaultOctave(ocv::Function)

Sets the default octave for `N` constructor (defualt value is 3).
"""
setDefaultOctave(ocv::Function) = ( global defaultOctave = ocv ; )

## basic constructor expressors

"""
	B(dur::IR=1//1)

An empty `Exp` with specified duration
"""
B(dur::IR=1//1) = ( @assert dur >= 0 ; Exp(dur,Atom[]) )

"""
	N(itv::Int)

A whole-note `Atom`; `itv` is the interval value on its scale. \\
Duration is `1//1`, offset is `0//1`, channel is `1`, velocity is `1.0`. \\
Uses default octave and scale (see `setDefaultScale` and `setDefaultOctave`).
"""
N(itv::Int) = Note(0//1,1//1,1,itv,defaultOctave,1.0,defaultScale)

## integers are converted to the default note
Base.convert(::Type{Atom}, v::Int) = N(v)

"""
	Bank(v::Int)

Bank select `Atom`; bank number is `v`.
"""
Bank(v::Int) = BankSel(0//1,1,v)

"""
	Prog(v::Int)

Program select `Atom`; program number is `v`.
"""
Prog(v::Int) = ProgSel(0//1,1,v)

"""
	Vol(v::Float64)

Volume set `Atom`; volume is `v`.
"""
Vol(v::Float64) = VolSet(0//1,1,v)

"""
	Wh()
	Wh(v::Float64)

Picth wheel bend `Atom`; bend value is `v` [-1,1], use Wh() for no bend (i.e. reset).
"""
Wh() = PitchWheel(0//1,1,2^13)
Wh(v::Float64) = ( v = clamp(v,-1,1) ; w = v==0 ? 2^13 : clamp(convert(Int,round(2^13*(1.0+v))),0,2^14-1) ; PitchWheel(0//1,1,w) )

"""
	Aft(v::Float64)

Channel aftertouch `Atom`; `v` [0,1] is the value.
"""
Aft(v::Float64) = ( w = clamp(convert(Int,round(v*127)),0,127) ; ChanAfter(0//1,1,w) )

"""
	Cn7(num::Int, v::Float64)
	Cn7(num::Int, v::Int)

7bit channel control `Atom`; `num` is the contol number, `v` is value ([0,1] if Float, {0..127} if Int).
"""
Cn7(num::Int, v::Float64) = ( w = clamp(convert(Int,round(v*127)),0,127) ; Control7(0//1,1,num,w) )
Cn7(num::Int, v::Int) = Control7(0//1,1,num,clamp(v,0,127))

"""
	Cn14(num::Int, v::Float64)
	Cn14(num::Int, v::Int)

14bit channel control `Atom`; `num` is the contol number, `v` is value ([0,1] if Float, {0..2^14-1} if Int).
"""
Cn14(num::Int, v::Float64) = ( w = clamp(convert(Int,round(v*(2^14-1))),0,2^14-1) ; Control14(0//1,1,num,w) )
Cn14(num::Int, v::Int) = Control14(0//1,1,num,clamp(v,0,2^14-1))

"""
	Pan()
	Pan(v::Float64)
	Pan(v::Int)

Channel pan `Atom`; use `Pan()` for center, `v` is the value ([-1,1] if Float, {0..2^14-1} if Int).
"""
Pan() = Ch14(10,2^13)
Pan(v::Float64) = ( v = clamp(v,-1,1) ; w = v==0 ? 2^13 : clamp(convert(Int,round(2^13*(1.0+v))),0,2^14-1) ;  Cn14(10,w) )
Pan(v::Int) = Cn14(10,v)

"channel sustain off `Atom`"
Sus0() = Cn7(64,0)

"channel sustain on `Atom`"
Sus1() = Cn7(64,127)

"all sounds off system message `Atom`"
SoundsOff() = Cn7(120,0)

"all contols reset system message `Atom`"
CntrlOff() = Cn7(121,0)

"all notes off system message `Atom`"
NotesOff() = Cn7(123,0)

"""
	Mark(v::AbstractString)

Track marker `Atom`; `v` is the marker name.
"""
Mark(v::AbstractString) = Marker(0//1,v)

"""
	Sig(num::Int=4, den::Int=3, clocks::Int=24, sub::Int=8)

Track time signature set `Atom`; \\
`num` is the numerator (defaults to 4), \\
`den` is the denumerator as a logarithm of 2 (defuaults to 3), \\
`clocks` is the divisions per quarter note (defualts to 24), \\
`sub` us the subdivisions of `clocks` (defualts to 8).
"""
Sig(num::Int=4, den::Int=3, clocks::Int=24, sub::Int=8) = TimeSignature(0//1,num,den,clocks,sub)

"""
	Tpo(bpm::Int)

Track tempo `Atom`; `bpm` is the beats per minute.
"""
Tpo(bpm::Int) = Tempo(0//1,bpm)

"""
	Inst(pg::Tuple{Int,Int}, [v::Float64])

Sets instrument to bank `pg[1]` program `pg[2]`. \\
If `v` is specified it also sets the volume to `v` which shoud be in the range `[0,1]`. \\
Named (bank,program) tuple constants can be found in the `SirenSeq.GenFluidSF144` and `SirenSeq.GenMuseScoreSF144` modules.
For example, `Instr(SirenSeq.GenMuseScoreSF144.iTubularBells,0.8)` selects the tubular-bells instrument from the general musescore 144 soundfont file and sets its volume to `0.8`.
"""
Inst(pg::Tuple{Int,Int}) = S(Bank(pg[1]),Prog(pg[2]))
Inst(pg::Tuple{Int,Int}, v::Float64) = S(Bank(pg[1]),Prog(pg[2]),Vol(v))

"""
	Da(v1::IR, v2::Float64, z) = D(v1,A(v2,z))

Dilate `z` by `v1` and accelirate it by `v2`.
"""
Da(v1::IR, v2::Float64, z) = D(v1,A(v2,z))

"""
	Con(z, b1::IR, b2::IR)
	Con(z, b2::IR) = Con(0,b2)
	Con(z) = Con(0,z.dur)

Contain `z`; this function makes sure that no events happen outside the time interval [`b1`,`b2`].
Events that cross this interval may be removed or shortened.
"""
Con(z) = ( x = atomOrExp(z); Con(x,0//1,x.dur) )
Con(z, b2::IR) = ( @assert b2>0 ; Con(z,0//1,b2) )
Con(z, b1::IR, b2::IR) = ( @assert b2>b1 ; contain(convert(Exp,atomOrExp(z)),Rat(b1),Rat(b2)) )

contain(x::Exp, b1::Rational{Int}, b2::Rational{Int}) = Exp(x.dur,atomsOp(a->contain(a,b1,b2),x))

function contain(a::Atom, b1::Rational{Int}, b2::Rational{Int})
	dur = a.dur ; ofs = a.ofs
	t1 = ofs ; t2 = ofs+dur
	if t1 < b1 ; dt = b1-t1 ; ofs += dt ; dur -= dt ; end
	if t2 > b2 ; dt = t2-b2 ; dur -= dt ; end
	if dur < 0 ; return Exp(0//1,Atom[]) ; end
	if dur == 0 && isa(a,Duratom) ; return Exp(0//1,Atom[]) ; end
	y = deepcopy(a) ; y.ofs=ofs ; y.dur=dur
	y
end

"""
	Kl(v::Rational{Int}, z)

Klip `z` by `v`.  The duration of all notes in `z` is multiplied by `v`.
A rest is inserted after each note where its length is `(1-v)` times the duration of the original note.
`0 < v < 1` must hold.
"""
Kl(v::Rational{Int}, z) = ( @assert 0<v<1 ; klip(v,atomOrExp(z)) )

klip(v::Rational{Int}, a::Note) = ( y = deepcopy(a) ; y.dur*=v ; Exp(a.dur,Atom[y]) )
klip(v::Rational{Int}, a::Atom) = a
klip(v::Rational{Int}, x::Exp) = Exp(x.dur,atomsOp(z->klip(v,z),x))

"""
	Spl(v::Int, z)

Split `z` `v` times.  Replaces each note in `z` by the same note repeated `v` times in the same time interval.
"""
Spl(v::Int, z) = ( @assert 1<v ; splitt(v,atomOrExp(x)) )

splitt(v::Int, a::Note) = ( y = deepcopy(a) ; y.dur/=v ; y.ofs=0//1 ; F(x.ofs,R(v,y)) )
splitt(v::Int, a::Atom) = a
splitt(v::Int, x::Exp) = atomsOp(z->splitt(v,z),x)


text:
"""
Return a decoded JSON file.
"""
function read_json(filename::AbstractString)
    open(filename, "r") do io
        JSON.parse(io)
    end
end

"""
Return a DataFrame from a CSV file.
"""
function read_dataframe(filename::AbstractString)
    open(filename, "r") do io
        DataFrames.DataFrame(CSV.File(io))
    end
end

"""
Return the SHA 256 hash of a file.
"""
function compute_sha256(filename::AbstractString)
    return open(filename) do io
        return bytes2hex(SHA.sha256(io))
    end
end

"""
Return the key for the given value
"""
function find_key_with_value(d, value)
    for (k, v) in d
        v == value && return k
    end
    error("dict does not have value == $value")
end

function compute_file_hash(path::String, files::Vector{String})
    data = Dict("files" => [])
    for file in files
        file_path = joinpath(path, file)
        file_info = Dict("filename" => file_path, "hash" => compute_sha256(file_path))
        push!(data["files"], file_info)
    end

    open(joinpath(path, HASH_FILENAME), "w") do io
        write(io, JSON.json(data))
    end
end

function compute_file_hash(path::String, file::String)
    return compute_file_hash(path, [file])
end

function read_file_hashes(path)
    data = open(joinpath(path, HASH_FILENAME), "r") do io
        JSON.parse(io)
    end

    return data["files"]
end

# writing a dictionary of dataframes to files

function write_data(vars_results::Dict, save_path::String; kwargs...)
    if :duals in keys(kwargs)
        name = "dual_"
    elseif :params in keys(kwargs)
        name = "parameter_"
    else
        name = ""
    end
    for (k, v) in vars_results
        file_path = joinpath(save_path, "$name$k.csv")
        if isempty(vars_results[k])
            @debug "$name$k is empty, not writing $file_path"
        else
            CSV.write(file_path, vars_results[k])
        end
    end
end

# writing a dictionary of dataframes to files and appending the time

function write_data(
    vars_results::Dict,
    time::DataFrames.DataFrame,
    save_path::AbstractString,
)
    for (k, v) in vars_results
        var = DataFrames.DataFrame()
        if size(time, 1) == size(v, 1)
            var = hcat(time, v)
        else
            var = v
        end
        file_path = joinpath(save_path, "$(k).csv")
        CSV.write(file_path, var)
    end
end

function write_data(
    data::DataFrames.DataFrame,
    save_path::AbstractString,
    file_name::String,
)
    if isfile(save_path)
        save_path = dirname(save_path)
    end
    file_path = joinpath(save_path, "$(file_name).csv")
    CSV.write(file_path, data)
    return
end

#Given the changes in syntax in ParameterJuMP and the new format to create anonymous parameters
function add_parameter(model::JuMP.Model, val::Number)
    param = JuMP.@variable(model, variable_type = PJ.Param())
    PJ.set_value(param, val)
    return param
end

function write_data(base_power::Float64, save_path::String)
    JSON.write(joinpath(save_path, "base_power.json"), JSON.json(base_power))
end

function _jump_value(input::JuMP.VariableRef)
    return JuMP.value(input)
end

function _jump_value(input::PJ.ParameterRef)
    return PJ.value(input)
end

function _jump_value(input::JuMP.ConstraintRef)
    return JuMP.dual(input)
end

"""
Creates a DataFrame from a JuMP DenseAxisArray or SparseAxisArray.

# Arguments
- `input_array`: JuMP DenseAxisArray or SparseAxisArray to convert
- `columns::Vector{Symbol}`: Required when there is only one axis which is data. Ignored if
  `input_array` includes an axis for device names.
"""
function axis_array_to_dataframe(
    input_array::JuMP.Containers.DenseAxisArray{Float64},
    columns = nothing,
)
    if length(axes(input_array)) == 1
        result = Vector{Float64}(undef, length(first(input_array.axes)))

        for t in input_array.axes[1]
            result[t] = input_array[t]
        end

        @assert columns !== nothing
        return DataFrames.DataFrame(columns[1] => result)

    elseif length(axes(input_array)) == 2
        result = Array{Float64, length(input_array.axes)}(
            undef,
            length(input_array.axes[2]),
            length(input_array.axes[1]),
        )
        names = Array{Symbol, 1}(undef, length(input_array.axes[1]))

        for t in input_array.axes[2], (ix, name) in enumerate(input_array.axes[1])
            result[t, ix] = input_array[name, t]
            names[ix] = Symbol(name)
        end

        return DataFrames.DataFrame(result, names)

    elseif length(axes(input_array)) == 3
        extra_dims = sum(length(axes(input_array)[2:(end - 1)]))
        extra_vars = [Symbol("S$(s)") for s in 1:extra_dims]
        result_df = DataFrames.DataFrame()
        names = vcat(extra_vars, Symbol.(axes(input_array)[1]))

        for i in input_array.axes[2]
            third_dim = collect(fill(i, size(input_array)[end]))
            result = Array{Float64, 2}(
                undef,
                length(last(input_array.axes)),
                length(first(input_array.axes)),
            )
            for t in last(input_array.axes),
                (ix, name) in enumerate(first(input_array.axes))

                result[t, ix] = input_array[name, i, t]
            end
            res = DataFrames.DataFrame(hcat(third_dim, result), :auto)
            result_df = vcat(result_df, res)
        end

        return DataFrames.rename!(result_df, names)

    else
        error("Dimension Number $(length(axes(input_array))) not Supported")
    end
end

function axis_array_to_dataframe(
    input_array::JuMP.Containers.DenseAxisArray{},
    columns = nothing,
)
    if length(axes(input_array)) == 1
        result = Vector{Float64}(undef, length(first(input_array.axes)))
        for t in input_array.axes[1]
            result[t] = _jump_value(input_array[t])
        end

        @assert columns !== nothing
        return DataFrames.DataFrame(columns[1] => result)
    elseif length(axes(input_array)) == 2
        result = Array{Float64, length(input_array.axes)}(
            undef,
            length(input_array.axes[2]),
            length(input_array.axes[1]),
        )
        names = Array{Symbol, 1}(undef, length(input_array.axes[1]))

        for t in input_array.axes[2], (ix, name) in enumerate(input_array.axes[1])
            result[t, ix] = _jump_value(input_array[name, t])
            names[ix] = Symbol(name)
        end

        return DataFrames.DataFrame(result, names)

    elseif length(axes(input_array)) == 3
        extra_dims = sum(length(axes(input_array)[2:(end - 1)]))
        extra_vars = [Symbol("S$(s)") for s in 1:extra_dims]
        result_df = DataFrames.DataFrame()
        names = vcat(extra_vars, Symbol.(axes(input_array)[1]))

        for i in input_array.axes[2]
            third_dim = collect(fill(i, size(input_array)[end]))
            result = Array{Float64, 2}(
                undef,
                length(last(input_array.axes)),
                length(first(input_array.axes)),
            )
            for t in last(input_array.axes),
                (ix, name) in enumerate(first(input_array.axes))

                result[t, ix] = _jump_value(input_array[name, i, t])
            end
            res = DataFrames.DataFrame(hcat(third_dim, result), :auto)
            result_df = vcat(result_df, res)
        end
        return DataFrames.rename!(result_df, names)
    else
        @warn(
            "Dimension Number $(length(axes(input_array))) not Supported, returning empty DataFrame"
        )
        return DataFrames.DataFrame()
    end
end

function axis_array_to_dataframe(
    input_array::JuMP.Containers.SparseAxisArray,
    columns = nothing,
)
    column_names = unique([(k[1], k[3]) for k in keys(input_array.data)])
    array_values = Vector{Vector{Float64}}(undef, length(column_names))
    for (ix, col) in enumerate(column_names)
        res = values(filter(v -> first(v)[[1, 3]] == col, input_array.data))
        array_values[ix] = PSI._jump_value.(res)
    end
    return DataFrames.DataFrame(array_values, Symbol.(column_names))
end

function axis_array_to_dataframe(input_array::Matrix{Float64}, columns)
    return DataFrames.DataFrame(input_array, columns)
end

function to_array(array::JuMP.Containers.DenseAxisArray)
    ax = axes(array)
    len_axes = length(ax)
    if len_axes == 1
        data = _jump_value.((array[x] for x in ax[1]))
    elseif len_axes == 2
        data = Array{Float64, 2}(undef, length(ax[2]), length(ax[1]))
        for t in ax[2], (ix, name) in enumerate(ax[1])
            data[t, ix] = _jump_value(array[name, t])
        end
        # TODO: this needs a better plan
        #elseif len_axes == 3
        #    extra_dims = sum(length(axes(array)[2:(end - 1)]))
        #    arrays = Vector{Array{Float64, 2}}()

        #    for i in ax[2]
        #        third_dim = collect(fill(i, size(array)[end]))
        #        data = Array{Float64, 2}(undef, length(last(ax)), length(first(ax)))
        #        for t in last(ax), (ix, name) in enumerate(first(ax))
        #            data[t, ix] = _jump_value(array[name, i, t])
        #        end
        #        push!(arrays, data)
        #    end
        #    data = vcat(arrays)
    else
        error("array axes not supported: $(axes(array))")
    end

    return data
end

function to_array(array::JuMP.Containers.DenseAxisArray{<:Number})
    length(axes(array)) > 2 && error("array axes not supported: $(axes(array))")
    return permutedims(array.data)
end

function to_array(array::JuMP.Containers.SparseAxisArray)
    columns = unique([(k[1], k[3]) for k in keys(array.data)])
    # PERF: can we determine the 2-d array size?
    tmp_data = Dict{Any, Vector{Float64}}()
    for (ix, col) in enumerate(columns)
        res = values(filter(v -> first(v)[[1, 3]] == col, array.data))
        tmp_data[col] = PSI._jump_value.(res)
    end

    data = Array{Float64, 2}(undef, length(first(values(tmp_data))), length(columns))
    for (i, column) in enumerate(columns)
        data[:, i] = tmp_data[column]
    end

    return data
end

to_array(array::Array) = array

# this ensures that the time_stamp is not double shortened
function find_var_length(es::Dict, e_list::Array)
    return size(es[Symbol(splitext(e_list[1])[1])], 1)
end

""" Returns the correct container spec for the selected type of JuMP Model"""
function container_spec(::Type{T}, axs...) where {T <: Any}
    return JuMP.Containers.DenseAxisArray{T}(undef, axs...)
end

""" Returns the correct container spec for the selected type of JuMP Model"""
function sparse_container_spec(::Type{T}, axs...) where {T <: Any}
    indexes = Base.Iterators.product(axs...)
    contents = Dict{eltype(indexes), Any}(indexes .=> 0)
    return JuMP.Containers.SparseAxisArray(contents)
end

function middle_rename(original::Symbol, split_char::String, addition::String)
    parts = split(String(original), split_char)
    parts[1] = parts[1] * "_" * addition
    return Symbol(join(parts, split_char))
end

"Replaces the string in `char` with the string`replacement`"
function replace_chars(s::String, char::String, replacement::String)
    return replace(s, Regex("[$char]") => replacement)
end

convert_for_path(x::Dates.DateTime) = replace(string(x), ":" => "-")

"Removes the string `char` from the original string"
function remove_chars(s::String, char::String)
    return replace_chars(s::String, char::String, "")
end

function get_available_components(::Type{T}, sys::PSY.System) where {T <: PSY.Component}
    return PSY.get_components(T, sys, x -> PSY.get_available(x))
end

function get_available_components(
    ::Type{PSY.RegulationDevice{T}},
    sys::PSY.System,
) where {T <: PSY.Component}
    return PSY.get_components(
        PSY.RegulationDevice{T},
        sys,
        x -> (PSY.get_available(x) && PSY.has_service(x, PSY.AGC)),
    )
end

"""
    check_file_integrity(path::String)

Checks the hash value for each file made with the file is written with the new hash_value to verify the file hasn't been tampered with since written

# Arguments
- `path::String`: this is the folder path that contains the results and the check.sha256 file
"""
function check_file_integrity(path::String)
    matched = true
    for file_info in read_file_hashes(path)
        filename = file_info["filename"]
        @info "checking integrity of $filename"
        expected_hash = file_info["hash"]
        actual_hash = compute_sha256(filename)
        if expected_hash != actual_hash
            @error "hash mismatch for file" filename expected_hash actual_hash
            matched = false
        end
    end

    if !matched
        throw(
            IS.HashMismatchError(
                "The hash value in the written files does not match the read files, results may have been tampered.",
            ),
        )
    end
end

to_namedtuple(val) = (; (x => getfield(val, x) for x in fieldnames(typeof(val)))...)

make_system_filename(sys::PSY.System) = "system-$(IS.get_uuid(sys)).json"
make_system_filename(sys_uuid::Base.UUID) = "system-$(sys_uuid).json"

function encode_symbol(::Type{T}, name1::AbstractString, name2::AbstractString) where {T}
    return Symbol(join((name1, name2, IS.strip_module_name(T)), PSI_NAME_DELIMITER))
end

function encode_symbol(
    ::Type{T},
    name1::AbstractString,
    name2::AbstractString,
) where {T <: PSY.Reserve}
    T_ = replace(IS.strip_module_name(T), "{" => "_")
    T_ = replace(T_, "}" => "")
    return Symbol(join((name1, name2, T_), PSI_NAME_DELIMITER))
end

function encode_symbol(::Type{T}, name1::Symbol, name2::Symbol) where {T}
    return encode_symbol(IS.strip_module_name(T), string(name1), string(name2))
end

function encode_symbol(::Type{T}, name::AbstractString) where {T}
    return Symbol(join((name, IS.strip_module_name(T)), PSI_NAME_DELIMITER))
end

function encode_symbol(::Type{T}, name::AbstractString) where {T <: PSY.Reserve}
    T_ = replace(IS.strip_module_name(T), "{" => "_")
    T_ = replace(T_, "}" => "")
    return Symbol(join((name, T_), PSI_NAME_DELIMITER))
end

function encode_symbol(::Type{T}, name::Symbol) where {T}
    return encode_symbol(T, string(name))
end

function encode_symbol(name::AbstractString)
    return Symbol(name)
end

function encode_symbol(name1::AbstractString, name2::AbstractString)
    return Symbol(join((name1, name2), PSI_NAME_DELIMITER))
end

function encode_symbol(name::Symbol)
    return name
end

function decode_symbol(name::Symbol)
    return split(String(name), PSI_NAME_DELIMITER)
end


text:
using JLLWrappers
using Pkg
using Test

# We use preferences only in Julia v1.6+
@static if VERSION >= v"1.6.0-DEV"
    using Preferences
end

module TestJLL end

@testset "JLLWrappers.jl" begin
    mktempdir() do dir
        Pkg.activate(dir)

        # actually use the development version of JLLWrappers
        Pkg.develop(PackageSpec(path=joinpath(@__DIR__, "..")))

        # Prepare some overrides for various products
        @static if VERSION >= v"1.6.0-DEV"
            set_preferences!(joinpath(dir, "LocalPreferences.toml"), "Vulkan_Headers_jll", "vulkan_hpp_path" => "foo")
            set_preferences!(joinpath(dir, "LocalPreferences.toml"), "HelloWorldC_jll", "goodbye_world_path" => "goodbye")
            set_preferences!(joinpath(dir, "LocalPreferences.toml"), "OpenLibm_jll", "libnonexisting_path" => "libreallynonexisting")
        end

        # Package with a FileProduct
        Pkg.develop(PackageSpec(path=joinpath(@__DIR__, "Vulkan_Headers_jll")))
        @test_nowarn @eval TestJLL using Vulkan_Headers_jll
        @test @eval TestJLL Vulkan_Headers_jll.is_available()
        @test isfile(@eval TestJLL vk_xml)
        @test isfile(@eval TestJLL Vulkan_Headers_jll.vk_xml_path)
        @test isfile(@eval TestJLL Vulkan_Headers_jll.get_vk_xml_path())
        @static if VERSION >= v"1.6.0-DEV"
            @test !isfile(@eval TestJLL vulkan_hpp)
            @test !isfile(@eval TestJLL Vulkan_Headers_jll.vulkan_hpp_path)
            @test !isfile(@eval TestJLL Vulkan_Headers_jll.get_vulkan_hpp_path())
        else
            @test isfile(@eval TestJLL vulkan_hpp)
            @test isfile(@eval TestJLL Vulkan_Headers_jll.vulkan_hpp_path)
            @test isfile(@eval TestJLL Vulkan_Headers_jll.get_vulkan_hpp_path())
        end
        @test isdir(@eval TestJLL Vulkan_Headers_jll.artifact_dir)
        @test isempty(@eval TestJLL Vulkan_Headers_jll.PATH[])
        @test occursin(Sys.BINDIR, @eval TestJLL Vulkan_Headers_jll.LIBPATH[])

        # Package with an ExecutableProduct
        Pkg.develop(PackageSpec(path=joinpath(@__DIR__, "HelloWorldC_jll")))
        @test_nowarn @eval TestJLL using HelloWorldC_jll
        if Sys.isfreebsd()
            @test @eval TestJLL !HelloWorldC_jll.is_available()
        else
            @test @eval TestJLL HelloWorldC_jll.is_available()
            @test "Hello, World!" == @eval TestJLL hello_world(h->readchomp(`$h`))
            @test isfile(@eval TestJLL HelloWorldC_jll.hello_world_path)
            @test isfile(@eval TestJLL HelloWorldC_jll.get_hello_world_path())
            @test isdir(@eval TestJLL HelloWorldC_jll.artifact_dir)
            @test !isempty(@eval TestJLL HelloWorldC_jll.PATH[])
            @test occursin(Sys.BINDIR, @eval TestJLL HelloWorldC_jll.LIBPATH[])
            @test !isfile(@eval TestJLL HelloWorldC_jll.goodbye_world_path)
            @test !isfile(@eval TestJLL HelloWorldC_jll.get_goodbye_world_path())
            @static if VERSION >= v"1.6.0-DEV"
                @test basename(@eval TestJLL HelloWorldC_jll.get_goodbye_world_path()) == "goodbye"
            else
                @test basename(@eval TestJLL HelloWorldC_jll.get_goodbye_world_path()) == "goodbye_world"
            end
        end

        # Package with a LibraryProduct
        Pkg.develop(PackageSpec(path=joinpath(@__DIR__, "OpenLibm_jll")))
        @test_nowarn @eval TestJLL using OpenLibm_jll
        @test @eval TestJLL OpenLibm_jll.is_available()
        @test exp(3.14) ≈ @eval TestJLL ccall((:exp, libopenlibm), Cdouble, (Cdouble,), 3.14)
        @test isfile(@eval TestJLL OpenLibm_jll.libopenlibm_path)
        @test isfile(@eval TestJLL OpenLibm_jll.get_libopenlibm_path())
        @test isdir(@eval TestJLL OpenLibm_jll.artifact_dir)
        @test isempty(@eval TestJLL OpenLibm_jll.PATH[])
        @test occursin(Sys.BINDIR, @eval TestJLL OpenLibm_jll.LIBPATH[])
        @test C_NULL == @eval TestJLL OpenLibm_jll.libnonexisting_handle

        @static if VERSION >= v"1.6.0-DEV"
            @test @eval TestJLL OpenLibm_jll.libnonexisting_path == "libreallynonexisting"
        else
            @test startswith(basename(@eval TestJLL OpenLibm_jll.libnonexisting_path), "libnonexisting")
        end

        # Issue #20
        if Sys.iswindows()
            @test Sys.BINDIR ∈ JLLWrappers.get_julia_libpaths()
        end

        # Package with a augment_platform!
        @static if VERSION >= v"1.6.0-DEV"
            Pkg.develop(PackageSpec(path=joinpath(@__DIR__, "LAMMPS_jll")))
            if VERSION < v"1.7"
                # On 1.6 we lazyily download the missing artifact, causing output to stderr
                @eval TestJLL using LAMMPS_jll
            else
                @test_nowarn @eval TestJLL using LAMMPS_jll
            end
            if @eval TestJLL LAMMPS_jll.is_available()
                @test @eval TestJLL LAMMPS_jll.abi == (Sys.iswindows() ? :MicrosoftMPI : :MPICH)
                @test isfile(@eval TestJLL LAMMPS_jll.liblammps_path)
                @test isfile(@eval TestJLL LAMMPS_jll.get_liblammps_path())
                @test isdir(@eval TestJLL LAMMPS_jll.artifact_dir)
                @test occursin(Sys.BINDIR, @eval TestJLL LAMMPS_jll.LIBPATH[])

                artifact_dir = @eval TestJLL LAMMPS_jll.artifact_dir

                if !Sys.iswindows()
                    @eval TestJLL LAMMPS_jll.set_abi(:MPItrampoline)
                    script = """
                    using Test
                    using LAMMPS_jll
                    @test LAMMPS_jll.abi == "MPItrampoline"
                    if !LAMMPS_jll.is_available()
                        @warn "LAMMPS_jll is not available, skipped test" LAMMPS_jll.host_platform
                        exit()
                    end
                    @test LAMMPS_jll.is_available()
                    @test isdir(LAMMPS_jll.artifact_dir)
                    @test LAMMPS_jll.artifact_dir !== $(repr(artifact_dir))
                    """
                    @test success(pipeline(`$(Base.julia_cmd()) --project=$dir -e $script`, stderr=stderr, stdout=stdout))
                end
            else
                platform = @eval TestJLL LAMMPS_jll.host_platform
                @warn "LAMMPS_jll is not available, skipped test" platform
            end
        end
    end
end

max_stars_repo_name: giordano/JLLWrappers

max_issues_repo_name: giordano/JLLWrappers

max_forks_repo_name: giordano/JLLWrappers


text:
## constructors of null versions of types for passing to C
for T in (:SuperMatrix,:NCformat,:NRformat,:SCformat,:SCPformat,:NCPformat,
          :DNformat,NRformat_loc,:e_node,:ExpHeader,:LU_stack_t,
          :mem_usage_t,:GlobalLU_t)
    @eval begin
        nv{}(::Type{$T}) = apply($T,map(zero,$T.types))  # null type T value
    end
end

function nv(::Type{superlu_options_t})
    v = apply(superlu_options_t,map(zero,superlu_options_t.types))
    ccall((:set_default_options,:libsuperlu),Void,(Ptr{superlu_options_t},),&v)
    v
end

function nv(::Type{SuperLUStat_t})
    v = apply(SuperLUStat_t,map(zero,SuperLUStat_t.types))
    ccall((:StatInit,:libsuperlu),Void,(Ptr{SuperLUStat_t},),&v)
    v
end

size(A::SuperLUMat,d::Integer) = d == 1 ? A.sm.nrow : (d == 2 ? A.sm.ncol : 1)
size(A::SuperLUMat) = (A.sm.nrow,A.sm.ncol)

Stype(A::SuperMatrix) = A.Stype # storage type (compressed column, dense, etc. - see Stype_t)
Dtype(A::SuperMatrix) = A.Dtype # data type (Cdouble, complex, etc. - see Dtype_t)
Mtype(A::SuperMatrix) = A.Mtype # matrix type (general, tril, triu, sym, etc. - see Mtype_t)

for T in (:NCformat,:NRformat,:SCformat,:SCPformat,:NCPformat,:DNformat,:NRformat_loc)
    @eval begin
        countnz(A::$T) = A.nnz
    end
end

function matstore(A::SuperMatrix)
    st = Stype(A)
    st == SLU_NC && return unsafe_load(convert(Ptr{NCformat},A.Store))
    st == SLU_NCP && return unsafe_load(convert(Ptr{NCPformat},A.Store))
    st == SLU_NR && return unsafe_load(convert(Ptr{NRformat},A.Store))
    st == SLU_SC && return unsafe_load(convert(Ptr{SCformat},A.Store))
    st == SLU_SCP && return unsafe_load(convert(Ptr{SCPformat},A.Store))
    st == SLU_SR && return unsafe_load(convert(Ptr{SRformat},A.Store))
    st == SLU_DN && return unsafe_load(convert(Ptr{DNformat},A.Store))
    st == SLU_NR_loc && return unsafe_load(convert(Ptr{NRformat_loc},A.Store))
    error("Unknown Stype = $st")
end

## could cheat here because nnz is always the first member of the type
countnz(A::SuperLUMat) = A.sm.Store == zero(Ptr{Void}) ? zero(Cint) : countnz(matstore(A))

function NCMat(A::SparseMatrixCSC{Cdouble})
    m,n = map(int32,size(A));
    ## turns out that you can't pass a matrix stored in symmetric form to dgssv
    ## symA = issym(A); AA = symA ? tril(A) : A
    ## mtyp = symA ? SLU_SYL : (istril(A) ? SLU_TRL : (istriu(A) ? SLU_TRU : SLU_GE))
    nz = int32(nfilled(A)); nzval = copy(A.nzval); rv = A.rowval; cp = A.colptr
    rowind = Array(int_t,nz); colptr = Array(int_t,n+i1)
                                        # use zero-based indices
    for i in 1:length(rowind); rowind[i] = int32(rv[i]) - i1; end
    for i in 1:length(colptr); colptr[i] = int32(cp[i]) - i1; end
    cres = nv(SuperMatrix)
    ccall((:dCreate_CompCol_Matrix,:libsuperlu),Void,
          (Ptr{SuperMatrix},Cint,Cint,Cint,Ptr{Cdouble},Ptr{Cint},Ptr{Cint},Stype_t,Dtype_t,Mtype_t),
          &cres,m,n,nz,nzval,rowind,colptr,SLU_NC,SLU_D,SLU_GE)
#    ccall((:dPrint_CompCol_Matrix,:libsuperlu),Void,(Ptr{Uint8},Ptr{SuperMatrix}),"A",&cres)
    NCMat(cres,colptr,rowind,nzval)
end

function DMat(B::StridedVecOrMat{Cdouble})
    mtyp = isa(B, Vector) ? SLU_GE : (istril(B) ? SLU_TRL : (istriu(B) ? SLU_TRU : SLU_GE))
    nzval = copy(B[:])
    cres = nv(SuperMatrix)
    ccall((:dCreate_Dense_Matrix,:libsuperlu),Void,
          (Ptr{SuperMatrix},Cint,Cint,Ptr{Cdouble},Cint,Cint,Cint,Cint),
          &cres,size(B,1),size(B,2),pointer(nzval),stride(B,2),SLU_DN,SLU_D,mtyp)
#    ccall((:dPrint_Dense_Matrix,:libsuperlu),Void,(Ptr{Uint8},Ptr{SuperMatrix}),"B",&cres)
    DMat(cres,nzval)
end

## Extract compiled-in sizes from libsuperlu
function sp_ienv(i::Integer)
    0 < i < 8 || error("i must be in 1,2,...,7")
    ccall((:sp_ienv,:libsuperlu),Cint,(Cint,),i)
end

function (\)(A::NCMat,B::StridedVecOrMat)
    A.sm.Dtype == SLU_D && eltype(B) == Cdouble || error("Both A and B must be Float64 arrays")
    m,n = size(A); m == n == size(B,1) || error("Dimension mismatch")
    BB = DMat(B)
    opts = nv(superlu_options_t)
    perm_r = zeros(Cint,n)
    perm_c = zeros(Cint,n)
    sTat = nv(SuperLUStat_t)
    L = nv(SuperMatrix)
    U = nv(SuperMatrix)
    inFo = zeros(Cint,1)
    ccall((:dgssv,:libsuperlu),Void,
          (Ptr{superlu_options_t},Ptr{SuperMatrix},Ptr{Cint},Ptr{Cint},
           Ptr{SuperMatrix},Ptr{SuperMatrix},Ptr{SuperMatrix},
           Ptr{SuperLUStat_t},Ptr{Cint}),
          &opts,&A.sm,perm_c,perm_r,&L,&U,&BB.sm,&sTat,inFo)
    inFo[1] == 0 || error("dgssv returned error code $(inFo[1])")
#    ccall((:StatPrint,:libsuperlu),Void,(Ptr{SuperLUStat_t},),&sTat)
    ccall((:StatFree,:libsuperlu),Void,(Ptr{SuperLUStat_t},),&sTat)
    U.Stype == SLU_NC || error("Matrix U should be of type NCformat (Stype = 0)")
    ccall((:Destroy_CompCol_Matrix,:libsuperlu),Void,(Ptr{SuperMatrix},),&U)
    L.Stype == SLU_SC || error("Matrix L should be of type SCformat (Stype = 3)")
    ccall((:Destroy_SuperNode_Matrix,:libsuperlu),Void,(Ptr{SuperMatrix},),&L)
    reshape(BB.nzval,size(B))
end

function lufact(A::NCMat,opts::superlu_options_t)
    m,n = size(A)
    perm_c = zeros(Cint,n)
                                        # determine fill-reducing permutation
    ccall((:get_perm_c,:libsuperlu),Void,(Cint,Ptr{SuperMatrix},Ptr{Cint}),
          opts.ColPerm,&(A.sm),perm_c)
    AC = nv(SuperMatrix)
    eTree = zeros(Cint,n)
                                        # apply perm_c to create AC and form eTree
    ccall((:sp_preorder,:libsuperlu),Void,
          (Ptr{superlu_options_t},Ptr{SuperMatrix},Ptr{Cint},Ptr{Cint},Ptr{SuperMatrix}),
          &opts,&(A.sm),perm_c,eTree,&AC)
    L = nv(SuperMatrix)
    U = nv(SuperMatrix)
    perm_r = zeros(Cint,m)
    inFo = zeros(Cint,1)
    sTat = nv(SuperLUStat_t)
    ccall((:dgstrf,:libsuperlu),Void,
          (Ptr{superlu_options_t},Ptr{SuperMatrix},Cint,Cint,Ptr{Cint},
           Ptr{Void},Cint,Ptr{Cint},Ptr{Cint},Ptr{SuperMatrix},Ptr{SuperMatrix},
           Ptr{SuperLUStat_t},Ptr{Cint}),
          &opts,&(A.sm),sp_ienv(2),sp_ienv(1),eTree,
          zero(Ptr{Void}),0,perm_c,perm_r,&L,&U,&sTat,inFo)
    println("info = ",inFo[1])
    inFo[1] == 0 || error("dgstrf returned error code $(inFo[1])")
    ccall((:StatPrint,:libsuperlu),Void,(Ptr{SuperLUStat_t},),&sTat)
    ccall((:dPrint_CompCol_Matrix,:libsuperlu),Void,
          (Ptr{Uint8},Ptr{SuperMatrix}), "U", &U)
    ccall((:dPrint_SuperNode_Matrix,:libsuperlu),Void,
          (Ptr{Uint8},Ptr{SuperMatrix}), "L", &L)
    showall(perm_c)
    showall(perm_r)
end
lufact(A::NCMat) = lufact(A,nv(superlu_options_t))


text:
# developed with Julia 1.1.1
#
# functions to simulate an EMS

function simulate_sites(controller::AbstractController,
                        path_to_save_folder::String,
                        path_to_price_csv_file::String,
                        path_to_metadata_csv_file::String,
                        path_to_test_data_folder::String,
                        path_to_train_data_folder::Union{String, Nothing}=nothing)

    mkpath(path_to_save_folder)
    prices = load_prices(path_to_price_csv_file)
    sites = load_sites(path_to_metadata_csv_file, path_to_test_data_folder,
        path_to_train_data_folder, path_to_save_folder)

    elapsed = 0.0

    for site in sites

        elapsed += @elapsed simulate_site(controller, site, prices)

    end

    group_all_simulations(sites)
    println("Terminating model simulation in $(elapsed) seconds")

    return nothing

end

function simulate_sites_parallel(controller::EMSx.AbstractController,
                                 path_to_save_folder::String,
                                 path_to_price_folder::String,
                                 path_to_metadata_csv_file::String,
                                 path_to_test_data_folder::String,
                                 path_to_train_data_folder::Union{String, Nothing}=nothing)

    mkpath(path_to_save_folder)
    prices = load_prices(path_to_price_folder)
    sites = load_sites(path_to_metadata_csv_file, path_to_test_data_folder,
        path_to_train_data_folder, path_to_save_folder)

    to_do = length(sites)

    @sync begin
        for p in workers()
            @async begin
                while true
                    idx = to_do
                    to_do -= 1
                    if idx <= 0
                        break
                    end
                    println("processing a new job - jobs left in queue : $(idx-1) / $(length(sites))")
                    _ = remotecall_fetch(simulate_site, p, controller, sites[idx], prices)
                end
            end
        end
    end

    group_all_simulations(sites)

end

function simulate_site(controller::AbstractController,
                       site::Site,
                       prices::Prices)

    test_data, site = load_site_data(site)
    controller = initialize_site_controller(controller, site, prices)
    periods = unique(test_data[!, :period_id])
    simulations = Simulation[]

    @showprogress for period_id in periods

        test_data_period = test_data[test_data.period_id .== period_id, :]
        period = Period(string(period_id), test_data_period, site)
        simulation = simulate_period(controller, period, prices)
        push!(simulations, simulation)

    end

    save_simulations(site, simulations)

    return nothing

end

function simulate_period(controller::AbstractController,
                           period::Period,
                           prices::Prices)

    horizon = size(period.data, 1) - 96 # test data: 24h of history lag + period data
    id = Id(period.site.id, period.id, prices.name, string(typeof(controller)))
    state_of_charge = 0.
    result = Result(horizon)
    timer = zeros(horizon)

    for t in 1:horizon

        information = Information(t, prices, period, state_of_charge)
        timing = @elapsed control = compute_control(controller, information)

        stage_cost, state_of_charge = apply_control(t, horizon, prices, period,
            state_of_charge, control)

        result.cost[t] = stage_cost
        result.soc[t] = state_of_charge
        timer[t] = timing

    end

    return Simulation(result, timer, id)

end

function apply_control(t::Int64, horizon::Int64,
                       prices::Prices, period::Period,
                       soc::Float64, control::Float64)
    """
    note on the load and pv values:
    at the end of the period, values at t+1 cannot be accessed: replaced by values at t=horizon
    with a minor impact since the optimal control is to empty the battery anyway
    """

    load = period.data[! ,:actual_consumption][min(t+96+1, horizon+96)]
    pv = period.data[!, :actual_pv][min(t+96+1, horizon+96)]
    net_energy_demand = load-pv

    stage_cost = compute_stage_cost(period.site.battery, prices, t, control, net_energy_demand)
    new_state_of_charge = compute_stage_dynamics(period.site.battery, soc, control)

    return stage_cost, new_state_of_charge

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1ServiceAccount <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    automountServiceAccountToken::Any # spec type: Union{ Nothing, Bool } # spec name: automountServiceAccountToken
    imagePullSecrets::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1LocalObjectReference} } # spec name: imagePullSecrets
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    secrets::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1ObjectReference} } # spec name: secrets

    function IoK8sApiCoreV1ServiceAccount(;apiVersion=nothing, automountServiceAccountToken=nothing, imagePullSecrets=nothing, kind=nothing, metadata=nothing, secrets=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ServiceAccount, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiCoreV1ServiceAccount, Symbol("automountServiceAccountToken"), automountServiceAccountToken)
        setfield!(o, Symbol("automountServiceAccountToken"), automountServiceAccountToken)
        validate_property(IoK8sApiCoreV1ServiceAccount, Symbol("imagePullSecrets"), imagePullSecrets)
        setfield!(o, Symbol("imagePullSecrets"), imagePullSecrets)
        validate_property(IoK8sApiCoreV1ServiceAccount, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiCoreV1ServiceAccount, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiCoreV1ServiceAccount, Symbol("secrets"), secrets)
        setfield!(o, Symbol("secrets"), secrets)
        o
    end
end # type IoK8sApiCoreV1ServiceAccount

const _property_map_IoK8sApiCoreV1ServiceAccount = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("automountServiceAccountToken")=>Symbol("automountServiceAccountToken"), Symbol("imagePullSecrets")=>Symbol("imagePullSecrets"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("secrets")=>Symbol("secrets"))
const _property_types_IoK8sApiCoreV1ServiceAccount = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("automountServiceAccountToken")=>"Bool", Symbol("imagePullSecrets")=>"Vector{IoK8sApiCoreV1LocalObjectReference}", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("secrets")=>"Vector{IoK8sApiCoreV1ObjectReference}")
Base.propertynames(::Type{ IoK8sApiCoreV1ServiceAccount }) = collect(keys(_property_map_IoK8sApiCoreV1ServiceAccount))
Swagger.property_type(::Type{ IoK8sApiCoreV1ServiceAccount }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1ServiceAccount[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ServiceAccount }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ServiceAccount[property_name]

function check_required(o::IoK8sApiCoreV1ServiceAccount)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ServiceAccount }, name::Symbol, val)
end


text:
using PackageCompiler
create_sysimage(:DelftBlue,
  sysimage_path=joinpath(@__DIR__,"..","DelftBlue.so"),
  precompile_execution_file=joinpath(@__DIR__,"warmup.jl"))


text:
module Bertweet

using PyCall

function __init__()
	py"""
	import emoji
	import transformers
	import torch
	"""
	global transformers = py"transformers"
	global model = transformers.AutoModel.from_pretrained("vinai/bertweet-base")
	global tokenizer = transformers.AutoTokenizer.from_pretrained("vinai/bertweet-base", normalized=true)
	global torch = py"torch"
	global demojize = py"emoji.demojize"
end

function setmodel!(modelname::String)
	global model = transformers.AutoModel.from_pretrained(modelname)
end

function settokenizer!(modelname::String)
	global tokenizer = transformers.AutoTokenizer.from_pretrained(modelname, normalized=true)
end

function encode(tweet::String)
	rawtweet = demojize(tweet)
	while occursin("::", rawtweet)
		rawtweet = replace(rawtweet, "::" => ":")
	end
	tokenizer.encode(rawtweet)
end

function pad!(tokenvec::Vector{Int}; maxsize::Int)
	while length(tokenvec) < maxsize
		push!(tokenvec, 1)
	end
end

function batchpad!(tokenbatch::Vector{Vector{Int}})
	pad_to = maximum(length.(tokenbatch))
	foreach(tokenbatch) do tokenvec
		pad!(tokenvec; maxsize=pad_to)
	end
end

function embed(tokenbatch::Vector{Vector{Int}})
	batchpad!(tokenbatch)
	tensorbatch = torch.tensor(tokenbatch)
	output = model(tensorbatch)
	pydecref(tensorbatch)
	cls_vecs = output["last_hidden_state"] .|> first
	embeddings = map(vec -> vec.tolist(), cls_vecs)
	output |> values .|> pydecref
	GC.gc(true)
	torch.cuda.empty_cache()
	return embeddings
end

end


text:
# setup clean and empty directories

user_dir = abspath(dirname(@__FILE__), "..", "myMatrixDepot")
data_dir = abspath(dirname(@__FILE__),"..","data")

function save_target(path::AbstractString)
    base = basename(path)
    dir = dirname(path)
    target = abspath(dir, string(base, ".saved"))
    while isdir(target)
        target = abspath(target, string(base, ".saved"))
    end
    if isdir(path)
        mv(path, target)
    end
    if basename(path) == "data"
        mkpath(path)
        ufm = "uf_matrices.html"
        cp(abspath(target, ufm), abspath(path, ufm))
        mmm = "mm_matrices.html"
        cp(abspath(target, mmm), abspath(path, mmm))
    end
    target
end

function revert_target(saved::AbstractString, orig::AbstractString)
    savetest = abspath(dirname(orig), string(basename(orig), ".test"))
    if isdir(orig)
        if isdir(savetest)
            mv(savetest, joinpath(orig, basename(savetest)))
        end
        mv(orig, savetest)
    end
    if isdir(saved)
        mv(saved, orig)
    end
end


text:
#File treating all the prediction functions

"""
Compute the mean of the predicted latent distribution of `f` on `X_test` for the variational GP `model`

Return also the variance if `covf=true` and the full covariance if `fullcov=true`
"""
function predict_f(model::VGP,X_test::AbstractMatrix{T};covf::Bool=true,fullcov::Bool=false) where T
    k_star = kernelmatrix.([X_test],[model.X],model.kernel)
    μf = k_star.*model.invKnn.*model.μ
    if !covf
        return model.nLatent == 1 ? μf[1] : μf
    end
    A = model.invKnn.*([I].-model.Σ.*model.invKnn)
    σ²f = []
    if fullcov
        k_starstar = kernelmatrix.([X_test],model.kernel)
        σ²f = Symmetric.(k_starstar .- k_star.*A.*transpose.(k_star) .+ convert(T,Jittering()).*[I])
    else
        k_starstar = kerneldiagmatrix.([X_test],model.kernel)
        σ²f = k_starstar .- opt_diag.(k_star.*A,k_star)
    end
    return model.nLatent == 1 ? (μf[1],σ²f[1]) : (μf,σ²f)
end

"""
Compute the mean of the predicted latent distribution of f on `X_test` for a sparse GP `model`
Return also the variance if `covf=true` and the full covariance if `fullcov=true`
"""
function predict_f(model::SVGP,X_test::AbstractMatrix{T};covf::Bool=true,fullcov::Bool=false) where T
    k_star = kernelmatrix.([X_test],model.Z,model.kernel)
    μf = k_star.*model.invKmm.*model.μ
    if !covf
        return model.nLatent == 1 ? μf[1] : μf
    end
    A = model.invKmm.*([I].-model.Σ.*model.invKmm)
    if fullcov
        k_starstar = kernelmatrix.([X_test],model.kernel)
        σ²f = Symmetric.(k_starstar .- k_star.*A.*transpose.(k_star))
    else
        k_starstar = kerneldiagmatrix.([X_test],model.kernel)
        σ²f = k_starstar .- opt_diag.(k_star.*A,k_star)
    end
    return model.nLatent == 1 ? (μf[1],σ²f[1]) : (μf,σ²f)
end

function predict_f(model::VGP{<:Likelihood,<:GibbsSampling},X_test::AbstractMatrix{T};covf::Bool=true,fullcov::Bool=false) where T
    k_star = kernelmatrix.([X_test],[model.X],model.kernel)
    f = [[k_star[min(k,model.nPrior)]*model.invKnn[min(k,model.nPrior)]].*model.inference.sample_store[k] for k in 1:model.nLatent]
    μf =  [vec(mean(hcat(f[k]...),dims=2)) for k in 1:model.nLatent]
    if !covf
        return model.nLatent == 1 ? μf[1] : μf
    end
    σ²f = []
    if fullcov
        k_starstar = kernelmatrix.([X_test],model.kernel)
        σ²f = Symmetric.(k_starstar .- k_star.*model.invKnn.*transpose.(k_star) .+  cov.(f))
    else
        k_starstar = kerneldiagmatrix.([X_test],model.kernel)
        σ²f = k_starstar .- opt_diag.(k_star.*model.invKnn,k_star) .+  diag.(cov.(f))
    end
    return model.nLatent == 1 ? (μf[1],σ²f[1]) : (μf,σ²f)
end

function predict_f(model::AbstractGP,X_test::AbstractVector{T};covf::Bool=false,fullcov::Bool=false) where T
    predict_f(model,reshape(X_test,length(X_test),1),covf=covf,fullcov=fullcov)
end

function predict_y(model::AbstractGP,X_test::AbstractVector)
    return predict_y(model,reshape(X_test,length(X_test),1))
end

"""
`predict_y(model::AbstractGP{<:RegressionLikelihood},X_test::AbstractMatrix)`

Return the predictive mean of `X_test`
"""
function predict_y(model::AbstractGP{<:RegressionLikelihood},X_test::AbstractMatrix)
    return predict_f(model,X_test,covf=false)
end

"""
`predict_y(model::AbstractGP{<:ClassificationLikelihood},X_test::AbstractMatrix)`

Return the predicted most probable sign of `X_test`
"""
function predict_y(model::AbstractGP{<:ClassificationLikelihood},X_test::AbstractMatrix)
    return [sign.(f) for f in predict_f(model,X_test,covf=false)]
end

"""
`predict_y(model::AbstractGP{<:MultiClassLikelihood},X_test::AbstractMatrix)`

Return the predicted most probable class of `X_test`
"""
function predict_y(model::AbstractGP{<:MultiClassLikelihood},X_test::AbstractMatrix)
    n = size(X_test,1)
    μ_f = predict_f(model,X_test,covf=false)
    return [model.likelihood.class_mapping[argmax([μ[i] for μ in μ_f])] for i in 1:n]
end

"""
`predict_y(model::AbstractGP{<:EventLikelihood},X_test::AbstractMatrix)`

Return the expected number of events for the locations `X_test`
"""
function predict_y(model::AbstractGP{<:EventLikelihood},X_test::AbstractMatrix)
    n = size(X_test,1)
    μ_f = predict_f(model,X_test,covf=false)
    return model.likelihood.λ.*((x->logistic.(x)).(μ_f))
end

function proba_y(model::AbstractGP,X_test::AbstractVector{T}) where {T<:Real}
    return proba_y(model,reshape(X_test,length(X_test),1))
end

"""
`proba_y(model::AbstractGP,X_test::AbstractMatrix)`

Return the probability distribution p(y_test|model,X_test) :

    - Tuple of vectors of mean and variance for regression
    - Vector of probabilities of y_test = 1 for binary classification
    - Dataframe with columns and probability per class for multi-class classification
"""
function proba_y(model::AbstractGP,X_test::AbstractMatrix{T}) where {T<:Real}
    μ_f,Σ_f = predict_f(model,X_test,covf=true)
    compute_proba(model.likelihood,μ_f,Σ_f)
end

function proba_y(model::VGP{<:MultiClassLikelihood,<:GibbsSampling},X_test::AbstractMatrix{T};nSamples::Int=200) where {T<:Real}
    k_star = kernelmatrix.([X_test],[model.X],model.kernel)
    f = [[k_star[min(k,model.nPrior)]*model.invKnn[min(k,model.nPrior)]].*model.inference.sample_store[k] for k in 1:model.nLatent]
    k_starstar = kerneldiagmatrix.([X_test],model.kernel)
    K̃ = k_starstar .- opt_diag.(k_star.*model.invKnn,k_star) .+ [zeros(size(X_test,1)) for i in 1:model.nLatent]
    nf = length(model.inference.sample_store[1])
    proba = zeros(size(X_test,1),model.nLatent)
    labels = Array{Symbol}(undef,model.nLatent)
    for i in 1:nf
        res = compute_proba(model.likelihood,getindex.(f,[i]),K̃,nSamples)
        if i ==  1
            labels = names(res)
        end
        proba .+= Matrix(res)
    end
    return DataFrame(proba/nf,labels)
end

function proba_y(model::VGP{<:ClassificationLikelihood,<:GibbsSampling},X_test::AbstractMatrix{T};nSamples::Int=200) where {T<:Real}
    k_star = kernelmatrix.([X_test],[model.X],model.kernel)
    f = [[k_star[min(k,model.nPrior)]*model.invKnn[min(k,model.nPrior)]].*model.inference.sample_store[k] for k in 1:model.nLatent]
    k_starstar = kerneldiagmatrix.([X_test],model.kernel)
    K̃ = k_starstar .- opt_diag.(k_star.*model.invKnn,k_star) .+ [zeros(size(X_test,1)) for i in 1:model.nLatent]
    nf = length(model.inference.sample_store[1])
    proba = [zeros(size(X_test,1)) for i in 1:model.nLatent]
    for i in 1:nf
        proba .+= compute_proba(model.likelihood,getindex.(f,[i]),K̃)
    end
    if model.nLatent == 1
        return proba[1]/nf
    else
        return proba./nf
    end
end

function compute_proba(l::Likelihood{T},μ::AbstractVector{<:AbstractVector{T}},σ²::AbstractVector{<:AbstractVector{T}}) where {T<:Real}
    compute_proba.([l],μ,σ²)
end

### TODO I Have to think about this solution

function compute_proba(l::Likelihood{T},μ::AbstractVector{T},σ²::AbstractVector{}) where {T<:Real}
    @error "Non implemented for the likelihood $l"
end

# "Return the mean of likelihood p(y*=1|X,x*) via the probit link with a linear model"
# function probitpredictproba(model::LinearModel,X_test::AbstractArray{T}) where {T<:Real}
#     if model.Intercept
#       X_test = [ones(T,size(X_test,1)) X_test]
#     end
#     n = size(X_test,1)
#     pred = zeros(n)
#     for i in 1:nTest
#       pred[i] = cdf(Normal(),(dot(X_test[i,:],model.μ))/(dot(X_test[i,:],model.Σ*X_test[i,:])+1))
#     end
#     return pred
# end

"""Return the mean of likelihood p(y*=1|X,x*) via the probit link with a GP model"""
function probitpredictproba(model::AbstractGP,X_test::AbstractArray{<:Real})
    m_f,cov_f = predict_f(model,X_test,covf=true)
    return broadcast((m,c)->cdf(Normal(),m/(c+1)),m_f,cov_f)
end

"Return the modified softmax likelihood given the latent functions"
function sigma_max(f::Vector{T},index::Integer) where {T<:Real}
    return logit(f[index])/sum(logit.(f))
end

"Return the modified softmax likelihood given the array of 'σ' and their sum (can be given via sumsig)"
function mod_soft_max(σ::Vector{T},sumsig::T=zero(T)) where {T<:Real}
    return sumsig == 0 ? σ./(sum(σ)) : σ./sumsig
end

"Return the gradient of the modified softmax likelihood given 'σ' and their sum (can be given via sumsig)"
function grad_mod_soft_max(σ::Array{T,1},sumsig::T=zero(T)) where {T<:Real}
    sumsig = sumsig == 0 ? sum(σ) : sumsig
    shortened_sum = sumsig.-σ
    sum_square = sumsig^2
    base_grad = (σ-(σ.^2))./sum_square
    n = size(σ,1)
    grad = zeros(n,n)
    for i in 1:n
        for j in 1:n
            if i==j
                grad[i,i] = shortened_sum[i]*base_grad[i]
            else
                grad[i,j] = -σ[i]*base_grad[j]
            end
        end
    end
    return grad
end

"Return the hessian of the modified softmax likelihood given 'σ' and their sum (can be given via sumsig)"
function hessian_mod_soft_max(σ::AbstractVector{T},sumsig::T=zero(T)) where {T<:Real}
    sumsig = sumsig == 0 ? sum(σ) : sumsig
    shortened_sum = sumsig.-σ
    sum_square = sumsig^2
    sum_cube = sumsig^3
    base_grad = (σ-σ.^2).*((1.0.-2.0.*σ).*sum_square-2.0.*(σ-σ.^2)*sumsig)./sum_cube
    n = size(σ,1)
    grad = zeros(n,n)
    for i in 1:n
        for j in 1:n
            if i==j
                grad[i,i] = shortened_sum[i]*base_grad[i]
            else
                grad[i,j] = -σ[i]*base_grad[j]
            end
        end
    end
    return grad
end


text:
module PATH2

    using ForwardDiff
    using FunctionWrappers: FunctionWrapper

    if isfile(joinpath(dirname(dirname(@__FILE__)), "deps", "deps.jl"))
      include(joinpath(dirname(dirname(@__FILE__)), "deps", "deps.jl"))
    else
      error("PATH2 not properly installed. Please run Pkg.build(\"PATH2\")")
    end

    export solveMCP, solveLCP, options

    # Global function pointers for the user-supplied function and jacobian evaluators.
    const user_f = Ref(FunctionWrapper{Vector{Cdouble}, Tuple{Vector{Cdouble}}}(identity))
    # The annotated SparseMatrixCSC return type will automatically convert the
    # jacobian into the correct sparse form for PATH
    const user_j = Ref(FunctionWrapper{SparseMatrixCSC{Cdouble, Cint}, Tuple{Vector{Cdouble}}}(identity))

    const cached_J = [convert(SparseMatrixCSC{Cdouble, Cint}, zeros(0, 0))]
    const cached_J_filled = Ref(false)

    const status =
       [  :Solved,                          # 1 - solved
          :StationaryPointFound,            # 2 - stationary point found
          :MajorIterationLimit,             # 3 - major iteration limit
          :CumulativeMinorIterationLimit,   # 4 - cumulative minor iteration limit
          :TimeLimit,                       # 5 - time limit
          :UserInterrupt,                   # 6 - user interrupt
          :BoundError,                      # 7 - bound error (lb is not less than ub)
          :DomainError,                     # 8 - domain error (could not find a starting point)
          :InternalError                    # 9 - internal error
      ]

    count_nonzeros(M::AbstractSparseMatrix) = nnz(M)
    count_nonzeros(M::AbstractMatrix) = count(x -> x != 0, M) # fallback for dense matrices

    function solveMCP(f_eval::Function, lb::Vector, ub::Vector, guess::Vector, var_name=C_NULL, con_name=C_NULL; nnz = length(lb)^2)
      j_eval = x -> ForwardDiff.jacobian(f_eval, x)
      return solveMCP(f_eval, j_eval, lb, ub, guess, var_name, con_name, nnz = nnz)
    end

    function solveMCP(f_eval::Function, j_eval::Function, lb::Vector, ub::Vector,guess::Vector, var_name=C_NULL, con_name=C_NULL; nnz = length(lb)^2)
      user_f[] = f_eval
      user_j[] = j_eval
      f_user_cb = cfunction(f_user_wrap, Cint, (Cint, Ptr{Cdouble}, Ptr{Cdouble}))
      j_user_cb = cfunction(j_user_wrap, Cint, (Cint, Cint, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}))

      n = length(lb)
      z = copy(guess)
      f = zeros(n)

      J0 = j_eval(z)
      if nnz == -1
            nnz = count_nonzeros(J0)
      end
      t = ccall( (:path_main, "libpath47julia"), Cint,
              (Cint, Cint,
               Ptr{Cdouble}, Ptr{Cdouble},
               Ptr{Cdouble}, Ptr{Cdouble},
               Ptr{Ptr{Cchar}}, Ptr{Ptr{Cchar}},
               Ptr{Void}, Ptr{Void}),
               n, nnz, z, f, lb, ub, var_name, con_name, f_user_cb, j_user_cb)

      remove_option_file()
      return status[t], z, f
    end

    function solveLCP(f_eval::Function, lb::AbstractVector, ub::AbstractVector,
                      var_name=C_NULL, con_name=C_NULL)
      J = ForwardDiff.jacobian(f_eval, lb)
      solveLCP(f_eval, J, lb, ub, var_name, con_name)
    end

    function solveLCP(f_eval::Function, M::AbstractMatrix,
                      lb::AbstractVector, ub::AbstractVector,
                      var_name=C_NULL, con_name=C_NULL)
      user_f[] = f_eval
      cached_J[] = M
      cached_J_filled[] = false
      f_user_cb = cfunction(f_user_wrap, Cint, (Cint, Ptr{Cdouble}, Ptr{Cdouble}))
      j_user_cb = cfunction(cached_j_user_wrap, Cint, (Cint, Cint, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}))

      n = length(lb)
      z = copy(lb)
      f = zeros(n)

      nnz = count_nonzeros(M)
      t = ccall( (:path_main, "libpath47julia"), Cint,
              (Cint, Cint,
               Ptr{Cdouble}, Ptr{Cdouble},
               Ptr{Cdouble}, Ptr{Cdouble},
               Ptr{Ptr{Cchar}}, Ptr{Ptr{Cchar}},
               Ptr{Void}, Ptr{Void}),
               n, nnz, z, f, lb, ub, var_name, con_name, f_user_cb, j_user_cb)
      remove_option_file()
      return status[t], z, f
    end

    function remove_option_file()
      if isfile("path.opt")
        rm("path.opt")
      end
    end

    function options(;kwargs...)
      opt_file = open("path.opt", "w")
      println(opt_file, "* Generated by PATHS2.jl. Do not edit.")
      for (key, value) in kwargs
        println(opt_file, key, " ", value)
      end
      close(opt_file)
    end

    ###############################################################################
    # wrappers for callback functions
    ###############################################################################
    # static int (*f_eval)(int n, double *z, double *f);
    # static int (*j_eval)(int n, int nnz, double *z, int *col_start, int *col_len,
          # int *row, double *data);
    function f_user_wrap(n::Cint, z_ptr::Ptr{Cdouble}, f_ptr::Ptr{Cdouble})
      z = unsafe_wrap(Array{Cdouble}, z_ptr, Int(n), false)
      f = unsafe_wrap(Array{Cdouble}, f_ptr, Int(n), false)
      f .= user_f[](z)
      return Cint(0)
    end

    function j_user_wrap(n::Cint, expected_nnz::Cint, z_ptr::Ptr{Cdouble},
                         col_start_ptr::Ptr{Cint}, col_len_ptr::Ptr{Cint},
                         row_ptr::Ptr{Cint}, data_ptr::Ptr{Cdouble})

      z = unsafe_wrap(Array{Cdouble}, z_ptr, Int(n), false)
      J::SparseMatrixCSC{Cdouble, Cint} = user_j[](z)
      if nnz(J) > expected_nnz
        error("Evaluated jacobian has more nonzero entries than were initially provided in solveMCP()")
      end
      load_sparse_matrix(J, n, expected_nnz, col_start_ptr, col_len_ptr, row_ptr, data_ptr)
      return Cint(0)
    end

    function cached_j_user_wrap(n::Cint, expected_nnz::Cint, z_ptr::Ptr{Cdouble},
                         col_start_ptr::Ptr{Cint}, col_len_ptr::Ptr{Cint},
                         row_ptr::Ptr{Cint}, data_ptr::Ptr{Cdouble})
      if !(cached_J_filled[])
        load_sparse_matrix(cached_J[], n, expected_nnz, col_start_ptr, col_len_ptr, row_ptr, data_ptr)
        cached_J_filled[] = true
      end
      return Cint(0)
    end

    function load_sparse_matrix(J::SparseMatrixCSC, n::Cint, expected_nnz::Cint,
                                col_start_ptr::Ptr{Cint}, col_len_ptr::Ptr{Cint},
                                row_ptr::Ptr{Cint}, data_ptr::Ptr{Cdouble})
      # Transfer data from the computed jacobian into the sparse format that PATH
      # expects. Fortunately, PATH uses a compressed-sparse-column storage which
      # is compatible with Julia's default SparseMatrixCSC format.

      # col_start in PATH corresponds to J.colptr[1:end-1]
      col_start = unsafe_wrap(Array{Cint}, col_start_ptr, Int(n), false)
      # col_len in PATH corresponds to diff(J.colptr)
      col_len = unsafe_wrap(Array{Cint}, col_len_ptr, Int(n), false)
      # row in PATH corresponds to rowvals(J)
      row = unsafe_wrap(Array{Cint}, row_ptr, Int(expected_nnz), false)
      # data in PATH corresponds to nonzeros(J)
      data = unsafe_wrap(Array{Cdouble}, data_ptr, Int(expected_nnz), false)

      @inbounds for i in 1:n
        col_start[i] = J.colptr[i]
        col_len[i] = J.colptr[i + 1] - J.colptr[i]
      end

      rv = rowvals(J)
      nz = nonzeros(J)
      num_nonzeros = nnz(J)
      @inbounds for i in 1:num_nonzeros
        row[i] = rv[i]
        data[i] = nz[i]
      end
    end

end # Module