text:
#=

This module provides code to compute Ramsey equilibria in a LQ economy with
distortionary taxation.  The program computes allocations (consumption,
leisure), tax rates, revenues, the net present value of the debt and other
related quantities.

Functions for plotting the results are also provided below.

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date: 2014-08-21

References
----------

Simple port of the file examples/lqramsey.py

http://quant-econ.net/lqramsey.html

=#
using QuantEcon
using PlotlyJS

abstract AbstractStochProcess

type ContStochProcess <: AbstractStochProcess
    A::Matrix
    C::Matrix
end

type DiscreteStochProcess <: AbstractStochProcess
    P::Matrix
    x_vals::Array
end

type Economy{SP <: AbstractStochProcess}
    bet::Real
    Sg::Matrix
    Sd::Matrix
    Sb::Matrix
    Ss::Matrix
    is_discrete::Bool
    proc::SP
end

type Path
    g
    d
    b
    s
    c
    l
    p
    tau
    rvn
    B
    R
    pi
    Pi
    xi
end

function compute_exog_sequences(econ::Economy, x)
    # Compute exogenous variable sequences
    Sg, Sd, Sb, Ss = econ.Sg, econ.Sd, econ.Sb, econ.Ss
    g, d, b, s = [squeeze(S * x, 1) for S in (Sg, Sd, Sb, Ss)]

    #= Solve for Lagrange multiplier in the govt budget constraint
    In fact we solve for nu = lambda / (1 + 2*lambda).  Here nu is the
    solution to a quadratic equation a(nu^2 - nu) + b = 0 where
    a and b are expected discounted sums of quadratic forms of the state. =#
    Sm = Sb - Sd - Ss

    return g, d, b, s, Sm
end

function compute_allocation(econ::Economy, Sm, nu, x, b)
    Sg, Sd, Sb, Ss = econ.Sg, econ.Sd, econ.Sb, econ.Ss

    # Solve for the allocation given nu and x
    Sc = 0.5 .* (Sb + Sd - Sg - nu .* Sm)
    Sl = 0.5 .* (Sb - Sd + Sg - nu .* Sm)
    c = squeeze(Sc * x, 1)
    l = squeeze(Sl * x, 1)
    p = squeeze((Sb - Sc) * x, 1)  # Price without normalization
    tau = 1 .- l ./ (b .- c)
    rvn = l .* tau

    return Sc, Sl, c, l, p, tau, rvn
end

function compute_nu(a0, b0)
    disc = a0^2 - 4a0*b0

    if disc >= 0
        nu = 0.5 *(a0 - sqrt(disc)) / a0
    else
        println("There is no Ramsey equilibrium for these parameters.")
        error("Government spending (economy.g) too low")
    end

    # Test that the Lagrange multiplier has the right sign
    if nu * (0.5 - nu) < 0
        print("Negative multiplier on the government budget constraint.")
        error("Government spending (economy.g) too low")
    end

    return nu
end

function compute_Pi(B, R, rvn, g, xi)
    pi = B[2:end] - R[1:end-1] .* B[1:end-1] - rvn[1:end-1] + g[1:end-1]
    Pi = cumsum(pi .* xi)
    return pi, Pi
end

function compute_paths(econ::Economy{DiscreteStochProcess}, T)
    # simplify notation
    bet, Sg, Sd, Sb, Ss = econ.bet, econ.Sg, econ.Sd, econ.Sb, econ.Ss
    P, x_vals = econ.proc.P, econ.proc.x_vals

    state = mc_sample_path(P, 1, T)
    x = x_vals[:, state]

    # Compute exogenous sequence
    g, d, b, s, Sm = compute_exog_sequences(econ, x)

    # compute a0, b0
    ns = size(P, 1)
    F = eye(ns) - bet.*P
    a0 = (F \ ((Sm * x_vals)'.^2))[1] ./ 2
    H = ((Sb - Sd + Sg) * x_vals) .* ((Sg - Ss)*x_vals)
    b0 = (F \ H')[1] ./ 2

    # compute lagrange multiplier
    nu = compute_nu(a0, b0)

    # Solve for the allocation given nu and x
    Sc, Sl, c, l, p, tau, rvn = compute_allocation(econ, Sm, nu, x, b)

    # compute remaining variables
    H = ((Sb - Sc)*x_vals) .* ((Sl - Sg)*x_vals) - (Sl*x_vals).^2
    temp = squeeze(F*H', 2)
    B = temp[state] ./ p
    H = squeeze(P[state, :] * ((Sb - Sc)*x_vals)', 2)
    R = p ./ (bet .* H)
    temp = squeeze(P[state, :] *((Sb - Sc) * x_vals)', 2)
    xi = p[2:end] ./ temp[1:end-1]

    # compute pi
    pi, Pi = compute_Pi(B, R, rvn, g, xi)

    Path(g, d, b, s, c, l, p, tau, rvn, B, R, pi, Pi, xi)
end

function compute_paths(econ::Economy{ContStochProcess}, T)
    # simplify notation
    bet, Sg, Sd, Sb, Ss = econ.bet, econ.Sg, econ.Sd, econ.Sb, econ.Ss
    A, C = econ.proc.A, econ.proc.C

    # Generate an initial condition x0 satisfying x0 = A x0
    nx, nx = size(A)
    x0 = nullspace((eye(nx) - A))
    x0 = x0[end] < 0 ? -x0 : x0
    x0 = x0 ./ x0[end]
    x0 = squeeze(x0, 2)

    # Generate a time series x of length T starting from x0
    nx, nw = size(C)
    x = zeros(nx, T)
    w = randn(nw, T)
    x[:, 1] = x0
    for t=2:T
        x[:, t] = A *x[:, t-1] + C * w[:, t]
    end

    # Compute exogenous sequence
    g, d, b, s, Sm = compute_exog_sequences(econ, x)

    # compute a0 and b0
    H = Sm'Sm
    a0 = 0.5 * var_quadratic_sum(A, C, H, bet, x0)
    H = (Sb - Sd + Sg)'*(Sg + Ss)
    b0 = 0.5 * var_quadratic_sum(A, C, H, bet, x0)

    # compute lagrange multiplier
    nu = compute_nu(a0, b0)

    # Solve for the allocation given nu and x
    Sc, Sl, c, l, p, tau, rvn = compute_allocation(econ, Sm, nu, x, b)

    # compute remaining variables
    H = Sl'Sl - (Sb - Sc)' *(Sl - Sg)
    L = Array(Float64, T)
    for t=1:T
        L[t] = var_quadratic_sum(A, C, H, bet, x[:, t])
    end
    B = L ./ p
    Rinv = squeeze(bet .* (Sb- Sc)*A*x, 1) ./ p
    R = 1 ./ Rinv
    AF1 = (Sb - Sc) * x[:, 2:end]
    AF2 = (Sb - Sc) * A * x[:, 1:end-1]
    xi =  AF1 ./ AF2
    xi = squeeze(xi, 1)

    # compute pi
    pi, Pi = compute_Pi(B, R, rvn, g, xi)

    Path(g, d, b, s, c, l, p, tau, rvn, B, R, pi, Pi, xi)
end

function gen_fig_1(path::Path)
    T = length(path.c)

    tr1, tr2, tr4 = GenericTrace[], GenericTrace[], GenericTrace[]

    # Plot consumption, govt expenditure and revenue
    push!(tr1, scatter(; y = path.rvn, legendgroup = "rev", marker_color = "blue", name = L"$\tau_t \ell_t$"))
    push!(tr1, scatter(; y = path.g, legendgroup = "gov", marker_color = "red", name = L"$g_t$"))
    push!(tr1, scatter(; y = path.c, marker_color = "green",name = L"$c_t$"))

    # Plot govt expenditure and debt
    push!(tr2, scatter(; x = 1:T, y = path.rvn, legendgroup = "rev", marker_color = "blue", showlegend = false, name = L"$\tau_t \ell_t$"))
    push!(tr2, scatter(; x = 1:T, y = path.g, legendgroup = "gov", marker_color = "red", showlegend = false, name = L"$g_t$"))
    push!(tr2, scatter(; x = 1:T-1, y = path.B[2:end], marker_color = "orange", name = L"$B_{t+1}$"))

    # Plot risk free return
    tr3 = scatter(; x =1:T, y = path.R - 1, marker_color = "pink", name = L"$R_{t - 1}$")

    # Plot revenue, expenditure and risk free rate
    push!(tr4, scatter(; x = 1:T, y = path.rvn, legendgroup = "rev", marker_color = "blue",  showlegend = false, name = L"$\tau_t \ell_t$"))
    push!(tr4, scatter(; x = 1:T, y = path.g, legendgroup = "gov", marker_color = "red", showlegend = false, name = L"$g_t$"))
    push!(tr4, scatter(; x = 1:T-1, y = path.pi, marker_color = "violet", name = L"$\pi_{t+1}$"))

    p1 = plot(tr1)
    p2 = plot(tr2)
    p3 = plot(tr3, Layout(; xaxis_title = "Time"))
    p4 = plot(tr4, Layout(; xaxis_title = "Time"))
    p = [p1 p2; p3 p4]
    relayout!(p, height = 900)

    p
end

function gen_fig_2(path::Path)
    T = length(path.c)

    # Plot adjustment factor
    p1 = plot(scatter(; x = 2:T, y = path.xi, name = L"$\xi_t$"))

    # Plot adjusted cumulative return
    p2 = plot(scatter(; x = 2:T, y = path.Pi, name = L"$\Pi_t$"), Layout(; xaxis_title = "Time"))

    p = [p1; p2]
    relayout!(p, height = 600)

    p
end

max_stars_repo_name: mwhchen/quantecon

max_issues_repo_name: mwhchen/quantecon

max_forks_repo_name: mwhchen/quantecon


text:
module TestPartition

using Test
using ThreadsX.Implementations: adhoc_partition

@testset for input_length in 1:100, partition_size in 1:100
    @test collect(adhoc_partition(1:input_length, partition_size)) ==
        collect(Iterators.partition(1:input_length, partition_size))
end

end  # module

ext: jl


text:
# Metrics
classify(x) = argmax.(eachcol(x))

function accuracy(model, data; batches = length(data), no_gpu::Bool = false, kwargs...)
    total_correct = 0
    total = 0
    for (i, (x_, y)) in enumerate(collect(data))
        i > batches && break
        x = no_gpu ? x_ : (x_ |> gpu |> track)
        target_class = classify(cpu(y))
        predicted_class = classify(cpu(model(x; kwargs...)[1] |> untrack))
        total_correct += sum(target_class .== predicted_class)
        total += length(target_class)
        x_ = nothing
        GC.gc(true)
    end
    return total_correct * 100 / total
end

function loglikelihood(model, data; batches = length(data))
    total_loglikelihood = 0.0f0
    total = 0
    for (i, x_) in enumerate(data)
        i > batches && break
        x = x_ |> gpu
        res = model(x |> track)
        total_loglikelihood += sum(res[1])
        total += size(x, 2)
        x_ = nothing
        GC.gc(true)
    end
    return total_loglikelihood / total
end


text:
"""
Stores the calculated metrics
"""
struct SmartReducer <: MetricStore
    state::Dict{Int, Real}
    momentum::Real
    SmartReducer(momentum=0.9) = new(Dict(), momentum)
end

function update!(r::SmartReducer, step::Int, value::Real)
    if haskey(r.state, step)
        r.state[step] = r.momentum * r.state[step] + (1-r.momentum) * value
    else
        r.state[step] = value
    end
end

"""
Get the history of a metric. The provided metric has the fully qualified symbol
and the returned value is a tuple of steps and values.

```julia
h = history(workout, :val_loss)
# returns e.g ([1000, 2000, 3000, 4000], [0.81, 0.73, 0.64, 0.61])
```

"""
function history(workout::Workout, symbol::Symbol)::Tuple
      h = workout.history[symbol].state
      steps = sort(collect(keys(h)))
      return (steps, [h[step] for step in steps])
end

"""
Binary accuracy calculation
"""
struct BinaryAccuracy
    name::Symbol
    threshold

    BinaryAccuracy(;threshold=0.5, name=:acc) = new(name,threshold)
end

function (a::BinaryAccuracy)(y_pred::Tensor, y_true::Tensor)
    y_pred = y_pred .> a.threshold
    return mean(y_true .== y_pred)
end

"""
Binary accuracy for a onehot classification.
"""
struct OneHotBinaryAccuracy
    name::Symbol

    OneHotBinaryAccuracy(name=:acc) = new(name)
end

function (a::OneHotBinaryAccuracy)(y_pred::Tensor, y_true::Tensor)
    y_pred, y_true = Knet.mat(y_pred), Knet.mat(y_true)

    y_pred = mapslices(argmax,Knet.value(Array(y_pred)),dims=1)
    y_true = Array(y_true)

    mean(y_pred .== y_true)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ a62a0de6-04b8-11eb-342f-bfea16900b85
begin
	using Pkg
	cd(joinpath(dirname(@__FILE__),".."))
    Pkg.activate(pwd())
	using PlutoUI
	using JuMP, GLPK, Tulip, Optim
	# Non-linear problem
	using Ipopt, Plots
	#using JuMP, Tulip, GLPK, LinearAlgebra
	#using Distributions, Plots, StatsPlots, LaTeXStrings, Measures
end

# ╔═╡ 86702ce8-04b7-11eb-22d1-b7f1437cf740
md"""
# Linear programming II
## Maximize flow in a network

We try to maximize the flow in a network using Linear Programming.

Let $N = (V, E)$ be a directed graph, where $V$ denotes the set of vertices and $E$ is the set of edges. Let $s ∈ V$ and $t ∈ V$ be the source and the sink of $N$, respectively. The capacity of an edge is a mapping $c : E \mapsto \mathbb{R}^+$, denoted by $c_{u,v}$ or $c(u, v)$. It represents the maximum amount of flow that can pass through an edge.

A flow is a mapping $f : E \mapsto \mathbb{R}^+$ , denoted by $f_{uv}$ or  $f(u, v)$, subject to the following two constraints:

* Capacity Constraint:

```math
\forall e \in E: f_{uv} \le c_{uv}
```

* Conservation of Flows:

```math
\forall v \in V\setminus\{s,t\} : \sum_{u:(u,v)\in E}f_{uv} = \sum_{w:(v,w)\in E} f_{vw}
```

We want to maximize the flow in the network, i.e.
```math
\max |f| = \max \sum_{v:(s,v)\in E}f_{sv} = \max \sum_{v:(v,t)\in E}f_{vt}
```

#### Setting:
Consider the following network:

$(PlutoUI.LocalResource("./Exercises/img/network.png"))

We want to:
1. Determine the maximal flow in the network
2. Be able to get a troughput of 35 from the source node to the sink node, whilst keeping the costs limited. Each link has a possible increase, with an associated cost (cf. table)

$(PlutoUI.LocalResource("./Exercises/img/networkcost.png"))
"""

# ╔═╡ caf5c6d8-04b8-11eb-1c04-0966207fbf28
# given set-up
begin
	# Topology and maximum flow matrix
	C = [0 13 6 10 0 0 0;
		 0 0  0  9 5 7 0;
		 0 0  0  8 0 0 0;
		 0 0  0  0 3 0 12;
		 0 0  0  0 0 4 6;
		 0 0  0  0 0 0 9;
		 0 0  0  0 0 0 0;
	];
	# extra capacity
	xcap = [ 0 6  4  3 0 0 0;
			 0 0  0  4 5 3 0;
			 0 0  0  5 0 0 0;
			 0 0  0  0 2 0 5;
			 0 0  0  0 0 2 4;
			 0 0  0  0 0 0 5;
			 0 0  0  0 0 0 0;
	];
	# cost per increased capacity
	xcost= [ 0 2.8  2.5  2.8 0   0   0;
			 0 0    0    2.5 3.1 1.6 0;
			 0 0    0    3.9 0   0   0;
			 0 0    0    0   2.8 0   1.6;
			 0 0    0    0   0   4.6 2.9;
			 0 0    0    0   0   0   1.8;
			 0 0    0    0   0   0   0;
	];
end

# ╔═╡ 5d0d7415-38fe-4c31-9c3d-59b2cbb54a8e
let
	model = Model(Ipopt.Optimizer)
	@variable(model, F[1:7,1:7])
	@constraint(model, F .<= C) # limite capacité
	for i = 2:6
		@constraint(model, sum(F[i,:]) == sum(F[:,i]))
	end
	@constraint(model, F .>= 0)
	@objective(model, Max, sum(F[1,:]))
	optimize!(model)
	abs.(round.(value.(F),digits=6))
end

# ╔═╡ 0f118a6f-9609-44b3-9e08-2e809f516377
let
	model = Model(Ipopt.Optimizer)
	@variable(model,F[1:size(C,1), 1:size(C,2)])
	# targt
	@objective(model, Max, sum(F[1,:])) #@objective(model, Max, sum(F[:,end]))
	# simplifier un peu
	for ind in findall(iszero,C)
		JuMP.fix(F[ind], 0., force=true)
	end
	# constraints
	for i = 2:6
		@constraint(model, sum(F[i,:]) == sum(F[:,i]))
	end
	@constraint(model, F .<= C)
	@constraint(model, F .>= 0)
	#
	optimize!(model)
	value.(F)
	# actuellement débit = 27
end

# ╔═╡ e6ba22db-4fc3-4a23-bb56-fa500ef616d2
let
	model = Model(Ipopt.Optimizer)
	@variable(model,F[1:size(C,1), 1:size(C,2)])
	# targt
	@objective(model, Min, sum(F .* xcost)) # faux raisonnement :-(
	# simplifier un peu
	for ind in findall(iszero,C)
		JuMP.fix(F[ind], 0., force=true)
	end
	# constraints
	for i = 2:6
		@constraint(model, sum(F[i,:]) == sum(F[:,i]))
	end
	@constraint(model, F .<= C .+ xcap)
	@constraint(model, F .>= 0)
	@constraint(model, sum(F[1,:]) == 35.)
	#
	optimize!(model)
	value.(F)
	# actuellement débit = 27
	Δ_cap = value.(F) - C
	ind_pos = Δ_cap .> 0
	sum((Δ_cap .* xcost)[ind_pos])
	Δ_cap
end

# ╔═╡ c38feffb-245b-4812-bcec-bdecd43d6c6c
let
	model = Model(Ipopt.Optimizer)
	@variable(model,F[1:size(C,1), 1:size(C,2)]) # capacité existante
	@variable(model,G[1:size(C,1), 1:size(C,2)]) # capacité supplémentaire
	# target
	@objective(model, Min, sum(G .* xcost)) # bon raisonnement :-)
	# simplifier un peu
	for ind in findall(iszero,C)
		JuMP.fix(F[ind], 0., force=true)
		JuMP.fix(G[ind], 0., force=true)
	end
	# constraints
	for i = 2:6
		@constraint(model, sum(F[i,:] + G[i,:]) == sum(F[:,i] + G[:,i]))
	end
	@constraint(model, F .<= C)
	@constraint(model, F .>= 0)
	@constraint(model, G .<= xcap)
	@constraint(model, G .>= 0)
	@constraint(model, sum(F[1,:]+ G[1,:]) == 35.)
	#
	optimize!(model)
	value.(F)
	# actuellement débit = 27
	sum(value.(G) .* xcost)
	value.(G)
end

# ╔═╡ d82dbb55-1d31-445a-b5e4-ec11ec4fcf65

# ╔═╡ 84804f16-8dfb-4bad-afe7-fe8ccd854e39

# ╔═╡ d4f240a2-73c5-4828-b728-d922db6148b8

# ╔═╡ 27080e2e-a7df-4acb-93a9-98f2352fd8b5

# ╔═╡ 508ae591-46f9-4dd9-97b8-08a5fd90e6e1

# ╔═╡ d3dcd770-058d-11eb-2829-fbb18e85f401
begin
	model = Model(Ipopt.Optimizer)
	@variable(model, F[1:size(C,1),1:size(C,2)])
	# gekende waarde op nul zetten
	for ind in findall(iszero, C)
		JuMP.fix(F[ind],0.)
	end
	@objective(model, Max, sum(F[1,:])) # of ook @objective(model, Max, sum(F[:,end]))
	# capacity constraint
	@constraint(model, F .<= C)
	# maintain flow
	for i = 2:6
		@constraint(model, sum(F[:,i]) == sum(F[i,:]))
	end

	optimize!(model)
end

# ╔═╡ c78e64e1-df73-4c6d-b705-2e1766d8ef3a
let
	# capaciteitsverhoging
	model = Model(Ipopt.Optimizer)
	@variable(model, F[1:size(C,1),1:size(C,2)])
	# gekende waarde op nul zetten
	for ind in findall(iszero, C)
		JuMP.fix(F[ind],0.)
	end
	# minimise total cost
	@objective(model, Min, sum(F .* xcost) ) # hier zit een redeneringsprobleem (!)
	# better?
	#helpfun(F) = F .- C .> 0
	#@objective(model, Min, filter(x-> x>0, (F .- C) .* xcost))

	# capacity constraint
	@constraint(model, F .<= C .+ xcap)
	@constraint(model, F .>= 0)
	# maintain flow
	for i = 2:6
		@constraint(model, sum(F[:,i]) == sum(F[i,:]))
	end
	# demand specifc throughput
	@constraint(model, sum(F[1,:]) == 35)

	optimize!(model)
	value.(F)
	# value of extra cost?
	Δ_cap = value.(F) .- C
	indpos = Δ_cap .> 0.
	sum(Δ_cap[indpos] .* xcost[indpos])
	#helpfun(F)
end

# ╔═╡ 74a20456-a91a-44de-9e4b-3e6d86440392
let
	# capaciteitsverhoging
	model = Model(Ipopt.Optimizer)
	@variable(model, F[1:size(C,1),1:size(C,2)]) # oorspronkelijk
	@variable(model, G[1:size(C,1),1:size(C,2)]) # extra
	# gekende waarde op nul zetten
	for ind in findall(iszero, C)
		JuMP.fix(F[ind],0.)
		JuMP.fix(G[ind],0.)
	end
	# minimise total cost of additional capacity
	@objective(model, Min, sum(G .* xcost) )

	# capacity constraint
	@constraint(model, F .<= C)
	@constraint(model, G .<= xcap)
	@constraint(model, F .>= 0)
	@constraint(model, G .>= 0)
	# maintain flow
	for i = 2:6
		@constraint(model, sum(F[:,i] .+ G[:,i]) == sum(F[i,:] .+ G[i,:]))
	end
	# demand specifc throughput
	@constraint(model, sum(F[1,:] .+ G[1,:]) == 35)

	optimize!(model)
	value.(F) .+ value.(G)
	# total cost:
	sum(value.(G) .* xcost)
end

# ╔═╡ 6ea03142-b487-4459-8391-cd50dae56340

# ╔═╡ 42371ec0-55fa-459b-bcb2-a6ea4363d876
value.(F)

# ╔═╡ 0a95f12e-c5eb-490a-bc7c-0cd7a70edddb
sum(value.(F)[1,:])

# ╔═╡ ef581235-7b51-4d2d-8097-7f6637c40b24
let
	model = Model(Ipopt.Optimizer)
	@variable(model, F[1:size(C,1),1:size(C,2)] >=0)
	# gekende waarde op nul zetten
	for ind in findall(iszero, C)
		JuMP.fix(F[ind],0, force=true)
	end
	@objective(model, Max, sum(F[1,:]))
	# capacity constraint
	@constraint(model, F .<= C .+ xcap)
	# maintain flow
	for i = 2:6
		@constraint(model, sum(F[:,i]) == sum(F[i,:]))
	end

	optimize!(model)
	value.(F)
end

# ╔═╡ 1283cdae-1a6b-4a4e-be2b-a9072bfc3972

# ╔═╡ 55ac8bb2-eb86-4536-9a89-245cf3883763

# ╔═╡ 7e852fdd-4363-4e15-ab24-e76448fbe1cf

# ╔═╡ a510792e-04c0-11eb-0f6e-a7e40dfe602c
md"""
## Optimizing an investment portfolio

In 1952 [Harry Max Markowitz](https://en.wikipedia.org/wiki/Harry_Markowitz) proposed a new approach for the optimization of an investment portfolio. This ultimately led to him winning the Nobel Prize in Economic Sciences in 1990. The idea is relatively simple:

Given a portfolio with $n$ stock proportions $S_1,S_2,\dots, S_n$, we want to maximize the return (=profit) and minimize the risk. The goal is to find the values $S_i$ that lead to either a minimum risk attribution with a minimal return or that lead to a maximum return attribution with a maximal risk.

Remembering that $\sigma^{2}_{\sum_{i=1}^{n}X_i}= \sum_{i=1}^{n}\sigma^2_{X_i} + \sum_{i \ne j}\text{Cov}(X_i,X_j) $, the risk can be expressed in terms of the covariance matrix $\Sigma$:

$$S^\mathsf{T} \Sigma S $$

The return can be expressed as:
$$\mu^\mathsf{T}S$$

Consider the following portfolio problem:
You are given the covariance matrix and expected returns and you want study several approaches. For each case you should formulate a proper Linear/Quadratic Programming form.
1. Ignore the risk and go for optimal investment (i.e. maximal return)
2. Same as (1), but a single stock can be at most 40% of the portfolio
3. Minimize the risk, with a lower bound on the return e.g. with at least 35% expected return
4. Make a graph for:
    * the minimal risk in fuction of the expected return.
    * the distribution of the portfolio with the minimal risk in function of the expected return
    * the final portfolio value in function of the expected return
"""

# ╔═╡ c4eca510-04c0-11eb-0e29-a5dcc1386bd5
# data for problem
begin
	P = [60; 127; 4; 50; 150; 20] # stock prices
	μ = [0.2; 0.42; 1.; 0.5; 0.46; 0.3] # expected returns
	Σ = [0.032 0.005 0.03 -0.031 -0.027 0.01;
		 0.005 0.1 0.085 -0.07 -0.05 0.02;
		 0.03 0.085 1/3 -0.11 -0.02 0.042;
		 -0.031 -0.07 -0.11 0.125 0.05 -0.06;
		 -0.027 -0.05 -0.02 0.05 0.065 -0.02;
		 0.01 0.02 0.042 -0.06 -0.02 0.08]; # covariance matrix
end

# ╔═╡ 0589b400-06a9-4f91-99aa-97e52f191dd0
let
	# Q1
	model = Model(Ipopt.Optimizer)
	# variables
	@variable(model,S[1:length(μ)]) # % van ieder aandeel in mijn portefeuille
	# objective
	@objective(model, Max, μ'*S) # max profit, risico bestaat niet
	# constraints
	@constraint(model, sum(S) == 1.) # total 100 %
	@constraint(model, S .>= 0.) # geen negatieve bijdrages
	# los het op
	optimize!(model)
	value.(S)

	# Q2
	model = Model(Ipopt.Optimizer)
	# variables
	@variable(model,S[1:length(μ)]) # % van ieder aandeel in mijn portefeuille
	# objective
	@objective(model, Max, μ'*S) # max profit, risico bestaat niet
	# constraints
	@constraint(model, sum(S) == 1.) # total 100 %
	@constraint(model, S .>= 0.) # geen negatieve bijdrages
	@constraint(model, S .<= 0.4) # diversificatie!
	# los het op
	optimize!(model)
	value.(S)

	# Q3
	model = Model(Ipopt.Optimizer)
	# variables
	@variable(model,S[1:length(μ)]) # % van ieder aandeel in mijn portefeuille
	# objective
	@objective(model, Min, S'*Σ*S) # min risk
	# constraints
	@constraint(model, μ'*S >= 0.35) # min return van 0.35
	@constraint(model, sum(S) == 1.) # total 100 %
	@constraint(model, S .>= 0.) # geen negatieve bijdrages
	@constraint(model, S .<= 0.4) # diversificatie!
	# los het op
	optimize!(model)
	value.(S), μ'*value.(S)
	termination_status(model)
end

# ╔═╡ e15f2647-04e0-4882-a193-5a63813852d4
let
	# ignorer le risque & maximiser le gain
	model = Model(Ipopt.Optimizer)
	@variable(model, S[1:length(P)])
	@constraint(model, S .<= 1.) # pas plus de 100% par action
	@constraint(model, S .>= 0.) # pas de proportion négative
	@constraint(model, sum(S) == 1.) # Σproportion = 1
	@objective(model, Max, μ' * S)
	optimize!(model)
	value.(S)
	# ignorer le risque & maximiser le gain & <=40% par action
	model = Model(Ipopt.Optimizer)
	@variable(model, S[1:length(P)])
	@constraint(model, S .<= 0.4) # pas plus de 100% par action
	@constraint(model, S .>= 0.) # pas de proportion négative
	@constraint(model, sum(S) == 1.) # Σproportion = 1
	@objective(model, Max, μ' * S)
	optimize!(model)
	value.(S)
	# minimiser le risque & gain minimal de 0.35
	model = Model(Ipopt.Optimizer)
	@variable(model, S[1:length(P)])
	@constraint(model, S .<= 0.4) # pas plus de 100% par action
	@constraint(model, S .>= 0.) # pas de proportion négative
	@constraint(model, sum(S) == 1.) # Σproportion = 1
	@constraint(model, μ' * S >= 0.35) # gain minimal à assurer
	@objective(model, Min, S' * Σ * S)

	optimize!(model)
	value.(S)

end

# ╔═╡ c6cb56d8-058d-11eb-2f71-73a5e5ca3dd0

# ╔═╡ ac5e426e-55f6-4ed0-98e4-19fb00ff1aec
let
	model = Model(Tulip.Optimizer)
	N = 20
	@variable(model, T[1:N] >= 0) # theorie part
	@variable(model, A[1:N] >= 0) # application part
	@constraint(model, T .<= 1) # logic
	@constraint(model, A .<= 1) # logic
	# teaching model
	a = 2; b = 3;
	@constraint(model, sum(A[1:N]) <= sum( T[b+1:N]))
	for i = 1:N
		if i <= b
			JuMP.fix(T[i], 1., force=true)
			JuMP.fix(A[i], 0., force=true)
		else
			@constraint(model, T[i] + A[i] == 1.)
		end
	end
	# emotional states
	s₀ = 0.
	α = -0.1; θ = 0.05; β = 1.4;
	S = [(1-θ)*s₀ + θ*(α*T[1] + β*A[1])]
	for i = 2:N
		push!(S, (1-θ)*S[i-1] + θ*(α*T[i-1] + β*A[i-1]))
	end
	@objective(model, Max, S[end])
	optimize!(model)
	termination_status(model)
	plot(value.(T))
	plot!(value.(A))
end

# ╔═╡ cb6898d6-0565-11eb-2a3c-1514a0fa4f50
md"""
## Optimal course planning
Suppose a professor teaches a course with $N=20$ lectures. We must decide how to split each lecture between theory and applications. Let $T_i$ and $A_i$ denote the fraction of the i$^{\text{th}}$ lecture devoted to theory and applications, for $i=1,\dots,N$. We can already determine the following:

```math
\forall i: T_i \ge 0, A_i \ge 0, T_i+A_i =1.
```

As you may know from experience, you need to cover a certain amount of theory before you can start doing applications. For this application consider the following model:

$$\sum_{i=1}^{N} A_i \le \phi \left( \sum_{i=1}^{N} T_i \right)$$

We interpret $\phi(u)$ as the cumulative amount of applications that can be covered, when the cumulative amount of theory covered is $u$. We will use the simple form $\phi(u) = a(u − b)$, with $a=2, b=3$, which means that no applications can be covered until $b$ lectures of the theory are covered; after that, each lecture of theory covered opens the possibility of covering a lecture on applications.

Psychological studies have shown that the theory-applications split affects the emotional state of students differently. Let $s_i$ denote the emotional state of a student after lecture $i$, with $s_i = 0$ meaning neutral, $s_i > 0$ meaning happy, and $s_i < 0$ meaning unhappy. Careful studies have shown that $s_i$ evolves via a linear recursion dynamic:

$$s_i =(1−\theta)s_{i−1} +\theta(\alpha T_i +\beta A_i)\text{ with }\theta \in[0,1]$$

with $s_0=0$. In order to make sure that the student leave with a good feeling at the end of the course, we try to maximize $s_N$, i.e. the emotional state after the last lecture.

Questions:
1. Determine the optimal split that leads to the most positive emotional state (for $\theta = 0.05, \alpha = -0.1, \beta = 1.4$);
2. Show the course repartition graphically
3. Determine values for $\alpha$ and $\beta$ that lead to a neutral result at the end of the course. Can you give an interpretation to these values?
"""

# ╔═╡ b76bae36-058d-11eb-216c-dfb8abd85772

# ╔═╡ Cell order:
# ╠═a62a0de6-04b8-11eb-342f-bfea16900b85
# ╠═5d0d7415-38fe-4c31-9c3d-59b2cbb54a8e
# ╟─86702ce8-04b7-11eb-22d1-b7f1437cf740
# ╠═0f118a6f-9609-44b3-9e08-2e809f516377
# ╠═e6ba22db-4fc3-4a23-bb56-fa500ef616d2
# ╠═c38feffb-245b-4812-bcec-bdecd43d6c6c
# ╠═caf5c6d8-04b8-11eb-1c04-0966207fbf28
# ╠═d82dbb55-1d31-445a-b5e4-ec11ec4fcf65
# ╠═84804f16-8dfb-4bad-afe7-fe8ccd854e39
# ╠═d4f240a2-73c5-4828-b728-d922db6148b8
# ╠═27080e2e-a7df-4acb-93a9-98f2352fd8b5
# ╠═508ae591-46f9-4dd9-97b8-08a5fd90e6e1
# ╠═d3dcd770-058d-11eb-2829-fbb18e85f401
# ╠═c78e64e1-df73-4c6d-b705-2e1766d8ef3a
# ╠═74a20456-a91a-44de-9e4b-3e6d86440392
# ╠═6ea03142-b487-4459-8391-cd50dae56340
# ╠═42371ec0-55fa-459b-bcb2-a6ea4363d876
# ╠═0a95f12e-c5eb-490a-bc7c-0cd7a70edddb
# ╠═ef581235-7b51-4d2d-8097-7f6637c40b24
# ╠═1283cdae-1a6b-4a4e-be2b-a9072bfc3972
# ╠═55ac8bb2-eb86-4536-9a89-245cf3883763
# ╠═7e852fdd-4363-4e15-ab24-e76448fbe1cf
# ╟─a510792e-04c0-11eb-0f6e-a7e40dfe602c
# ╠═0589b400-06a9-4f91-99aa-97e52f191dd0
# ╠═c4eca510-04c0-11eb-0e29-a5dcc1386bd5
# ╠═e15f2647-04e0-4882-a193-5a63813852d4
# ╠═c6cb56d8-058d-11eb-2f71-73a5e5ca3dd0
# ╠═ac5e426e-55f6-4ed0-98e4-19fb00ff1aec
# ╟─cb6898d6-0565-11eb-2a3c-1514a0fa4f50
# ╠═b76bae36-058d-11eb-216c-dfb8abd85772

max_stars_repo_path: Exercises/PS05 - Linear programming II.jl


text:
Random.seed!(7103)
rng = MersenneTwister(123)

@test_throws DomainError CompleteRandomisation([-0.5, 0.5])

CR1 = CompleteRandomisation([1, 1])
@test narms(CR1) == 2
@test target(CR1) == [0.5, 0.5]
@test sequence(CR1) == zeros(0)
@test dist(CR1) == zeros(narms(CR1))
@test imbalance(CR1) == 0.0
@test predictability(CR1) == 0.0
@test prob(CR1) == target(CR1)

# Check that randomise does not alter state
res = randomise(rng, CR1)
@test narms(CR1) == 2
@test target(CR1) == [0.5, 0.5]
@test sequence(CR1) == zeros(0)
@test dist(CR1) == zeros(narms(CR1))
@test imbalance(CR1) == 0.0
@test predictability(CR1) == 0.0

# Check that randomise! does alter state
res1 = randomise!(rng, CR1)
@test narms(CR1) == 2
@test target(CR1) == [0.5, 0.5]
@test sequence(CR1) == [res1]
@test imbalance(CR1) == euclidean(target(CR1), [0, 1])
@test dist(CR1) == Int.([x == res1 for x in 1:2])

res2 = randomise!(rng, CR1, 10)
@test length(res2) == 10


text:
#
# This module creates a platform independent window with a canvas inside.
# It uses OpenGL and GR libraries for that. Has an ability to draw a color pixels
# on a canvas. Uses RGB color format for drawing. It implements lazy update
# logic. It means, that canvas will be updated only after calling update()
# method. So, you may draw many dots, but user will not see them. Only after
# calling update() all dots will appear.
#
# Issues:
#
# Usage:
#   import OpenGlWindow
#   ...
#   win = OpenGlWindow.create(300, 300, 3)
#   CanvasWindow.dot(win, 20, 20, UInt16(11197883)) # R=AA,G=DD,B=BB
#   OpenGlWindow.dot(win, 30, 30, UInt16(11197883)) # R=AA,G=DD,B=BB
#   OpenGlWindow.update(win)                        # Two dots will be shown
#   ...
#   OpenGlWindow.destroy(win)
#
# @author DeadbraiN
# TODO: check if this mosule platform independent
# TODO: describe 1256 colors limit
#
module OpenGlWindow
  import GR
  import Dots
  import Helper

  export Window
  export create
  export setColor
  export dot
  export update
  export title
  export destroy
  #
  # Height of footer in pixels. Footer is an area for text, because window
  # title is unaccessible for us :( It should be at the bottom, because it
  # simplier to draw pixels without adding _FOOTER_HEIGHT all the time.
  #
  const _FOOTER_HEIGHT = 50.0
  #
  # Data type, which contains window related data: Window, Canvas, Context,...
  # Is used in all public methods as a canvas for drawing.
  #
  type Window
    scale::Int
    xBuf::Array{Int, 1}
    yBuf::Array{Int, 1}
    width::Int
    height::Int
    ratio::Float64
    running::Bool
  end
  #
  # Creates window and returns window related data type, which is used
  # in all public methods of current module. Sets default background
  # color according to global configuration.
  # @param width Window width in pixels
  # @param height Window height in pixels
  # @param scale Scale of canvas 1:1 by default
  # @param title Window title
  # @return Window object
  #
  function create(width::Int, height::Int, scale::Int, title::String = "")
    local emptyColor::Int = Dots.INDEX_EMPTY
    local win::Window = Window(scale, Int[0], Int[0], width, height, 1.0, true)
    local wndWidth::Int = width * scale
    local wndHeight::Int = height * scale + _FOOTER_HEIGHT
    local ratio::Float64 = _createWindow(win, wndWidth, wndHeight)

    GR.setwindow(1, wndWidth, 1, wndHeight)
    GR.setviewport(0, 1, 0, ratio)
    #
    # Fill background with "empty" color
    #
    GR.setcolorrep(emptyColor, Dots.COLOR_EMPTY...)
    GR.setlinewidth(1)
    GR.setlinetype(GR.LINETYPE_SOLID)
    GR.setlinecolorind(emptyColor)
    GR.settextcolorind(Int(Dots.INDEX_TEXT))
    #
    # Settings for fillrect() function
    #
    GR.setfillcolorind(emptyColor)
    GR.setfillintstyle(GR.INTSTYLE_SOLID)
    #
    # Settings for dots drawing with polymarker() function
    #
    GR.setmarkertype(GR.MARKERTYPE_DOT)
    GR.setmarkersize(1)
    GR.clearws()
    GR.fillrect(1, wndWidth, 1, wndHeight)
    #
    # TODO: This is a fix for an issue in GR library. fillrect() doesn't
    # TODO: draw right and bottom lines correctly. So, we have to add these
    # TODO: lines "by hands".
    #
    GR.drawrect(1, wndWidth, 1, wndHeight)
    #
    # Sets text height in title, otherwise it will be sized depending
    # on window size. So now it doesn't depend on window size.
    #
    GR.setcharheight(Float64(wndHeight / 1920) * 0.00001)
    OpenGlWindow.title(win, title)
    update(win)

    win
  end
  #
  # Redefines color index with RGB components.
  # @param color Color index we want to change
  # @param r Red component
  # @param g Green component
  # @param b Blue component
  #
  function setColor(color::Int, r::Float64, g::Float64, b::Float64)
    GR.setcolorrep(color, r, g, b)
  end
  #
  # Draws one dot (point) on the canvas with specified color
  # @param win Windows type
  # @param x X coordinate of the point
  # @param y Y coordinate of the point
  # @param color Color of the dot. We use only last three bytes (24bits) of four.
  #
  function dot(win::Window, x::UInt16, y::UInt16, color::UInt16)
    #
    # We use GR library. It supports only Dots.INDEX_MAX_COLOR colors. We have to check it
    #
    if color > Dots.INDEX_MAX_COLOR
      Helper.warn(string("Unsupported color index ", color))
      color = Dots.INDEX_MAX_COLOR
    end
    if win.scale > 1
      x = (x - UInt16(1)) * UInt16(win.scale + 1)
      y = win.height * win.scale - (y - UInt16(1)) * UInt16(win.scale + 1)
      GR.setlinecolorind(Int(color))
      GR.setfillcolorind(Int(color))
      GR.fillrect(Int(x), Int(x) + win.scale, Int(y), Int(y) + win.scale)
    else
      GR.setmarkercolorind(Int(color))
      win.xBuf[1] = Int(x)
      win.yBuf[1] = Int(y)
      GR.polymarker(win.xBuf, win.yBuf)
    end
  end
  #
  # Sets title to specified string. We can't change real window's title,
  # so we have to add a header to the window, where all text messages
  # will be posted. This header is above the world as a part of a canvas.
  # @param win Window data type
  # @param title String title
  #
  function title(win::Window, title::String)
    local ymax::Float64 = Float64(win.height * win.scale + _FOOTER_HEIGHT)
    local ycoef::Float64 = win.ratio / ymax

    GR.setlinecolorind(Int(Dots.INDEX_EMPTY))
    GR.setfillcolorind(Int(Dots.INDEX_EMPTY))
    GR.fillrect(1, win.width * win.scale + win.scale, ymax - _FOOTER_HEIGHT + win.scale, ymax)
    # TODO: fix this hardcode (15.0)
    GR.text(0.01, win.ratio - (ycoef * _FOOTER_HEIGHT) + ycoef * 19.0, title)
  end
  #
  # Updates the canvas. It's not nessesary to update it after
  # every drawing. It's better to update it after several dots
  # are drown.
  # @param win Current window
  #
  function update(win::Window) GR.updatews() end
  #
  # Destroys OpenGL windows by quitting from main loop
  # @param win Window object to destroy
  #
  function destroy(win::Window)
    win.running = false
  end
  #
  # Creates window with specified width and height in pixels.
  # @param wnd Window data object
  # @param width Window width
  # @param height Window height
  # @return {Float64} ratio
  #
  function _createWindow(win::Window, width::Int, height::Int)
    local mWidth::Float64, mHeight::Float64, scrWidth::Int32, scrHeight::Int32
    local ratio::Float64
    local mSize::Float64

    mWidth, mHeight, scrWidth, scrHeight = GR.inqdspsize()
    if width >= height
    	ratio = float(height) / width
    	mSize = mWidth * width / scrWidth
    	GR.setwsviewport(0, mSize, 0, mSize * ratio)
    	GR.setwswindow(0, 1, 0, ratio)
    else
    	ratio = float(width) / height
    	mSize = mHeight * height / scrHeight
    	GR.setwsviewport(0, mSize * ratio, 0, mSize)
    	GR.setwswindow(0, ratio, 0, 1)
    end

    (win.ratio = ratio)
  end
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Neighborhood

A neighborhood is a geometry that is not attached to any specific
point in space, and is free to slide over a domain of interest.
"""
abstract type Neighborhood end

"""
    MetricBall

A metric ball is a neighborhood that can be expressed in terms
of a metric and a set of radii.
"""
abstract type MetricBall <: Neighborhood end

"""
    metric(ball)

Return the metric of the metric `ball`.
"""
function metric end

"""
    radii(ball)

Return the radii of the metric `ball`.
"""
function radii end

# ----------------
# IMPLEMENTATIONS
# ----------------

include("neighborhoods/isotropicball.jl")
include("neighborhoods/anisotropicball.jl")


text:
#################################################################################################
# Functions
#################################################################################################
function get_ϕs_eff_rng(eff_type, config)

    # Pick max values
    if eff_type == "mml"
        ϕ1_max = [-1.0,-10.0,2.0]
        ϕ2_max = [+1.0,+10.0,2.0]
    elseif eff_type == "nme"
        ϕ1_max = [0.0,0.0,0.0]
        ϕ2_max = [0.0,0.0,5.0]
    elseif eff_type == "cmd"
        ϕ1_max = [+1.0,+10.0,2.0]
        ϕ2_max = [+1.0,+10.0,2.0]
    else
        CpelNano.print_log("Effect type $(eff_type) type not recognized.)")
    end

    # Loop over percentages
    i = 1
    eff_perc = 0.0
    eff_size_arr = Vector{Float64}()
    ϕ1_arr = Vector{Vector{Float64}}()
    ϕ2_arr = Vector{Vector{Float64}}()
    next_eff_size = minimum(EF_PERC_RNG)
    eff_size_check = fill(false, length(EF_PERC_RNG))
    while ! all(eff_size_check)

        # Scale vector
        if eff_type == "mml"
            ϕ1 = [eff_perc * ϕ1_max[1],eff_perc * ϕ1_max[2],ϕ1_max[3]]
            ϕ2 = [eff_perc * ϕ2_max[1],eff_perc * ϕ2_max[2],ϕ2_max[3]]
        elseif eff_type == "nme"
            ϕ1 = eff_perc * ϕ1_max
            ϕ2 = eff_perc * ϕ2_max
        else
            ϕ1 = ϕ1_max * ( 1.0 .+ eff_perc)
            ϕ2 = ϕ2_max / ( 1.0 .+ eff_perc)
        end

        # Baseline struct
        rs_g1 = create_baseline_struct(config)
        rs_g2 = create_baseline_struct(config)

        # Set vector
        rs_g1.ϕhat = ϕ1
        rs_g2.ϕhat = ϕ2

        # Compute stats
        proc_rep_rs!(rs_g1, config)
        proc_rep_rs!(rs_g2, config)

        # Non-empty analysis regions
        kp_in = [int != 0:0 for int in rs_g1.nls_rgs.cpg_ind]

        # Effect size
        if eff_type == "mml"
            eff_size = mean(abs.(rs_g1.mml[kp_in] - rs_g2.mml[kp_in]))
        elseif eff_type == "nme"
            eff_size = mean(abs.(rs_g1.nme[kp_in] - rs_g2.nme[kp_in]))
        elseif eff_type == "cmd"
            eff_size = mean(CpelNano.comp_cmd(rs_g1, rs_g2)[kp_in])
        end

        println("eff_size=$(eff_size)")

        # Record vector
        if eff_size >= next_eff_size

            # Record vector
            push!(ϕ1_arr, ϕ1)
            push!(ϕ2_arr, ϕ2)
            push!(eff_size_arr, eff_size)

            # Break
            i == length(EF_PERC_RNG) && break

            # Next effect size
            next_eff_size = EF_PERC_RNG[i + 1]

            # Increase counter
            i += 1

        end

        # Increase effect percentage
        eff_perc += 0.0005

    end

    # Return vector
    return eff_size_arr, ϕ1_arr, ϕ2_arr

end
@everywhere function create_baseline_struct(config::CpelNano.CpelNanoConfig)

    # Initialize
    rs = CpelNano.RegStruct()

    # Read in file
    rs_data = readdlm("$(aim_dir)/estimation_region_chr_chr1_42900254_42903250.txt")

    # Get values
    rs.chrst = 42900254
    rs.chrend = 42903250
    rs.cpg_pos = Int64.(rs_data[:,1])
    rs.ρn = Float64.(rs_data[:,2])
    rs.dn = Float64.(rs_data[:,3])

    # Get rest of values
    rs.N = length(rs.cpg_pos)
    rs.cpg_grps = CpelNano.get_grps_from_cgs(rs.cpg_pos, config.min_grp_dist)
    rs.Nl = map(grp -> length(grp.cpg_ind), rs.cpg_grps)
    rs.ρl = CpelNano.get_ρl(rs.ρn, rs.cpg_grps)
    rs.dl = CpelNano.get_dl(rs.cpg_pos, rs.cpg_grps)
    rs.L = length(rs.cpg_grps)

    # Get analysis region info
    CpelNano.get_nls_reg_info!(rs, config)

    # Return struct
    return rs

end
@everywhere function cp_rs_flds!(rs::CpelNano.RegStruct, rs_aux::CpelNano.RegStruct)

    # Fill fields
    rs_aux.L = rs.L
    rs_aux.Nl = rs.Nl
    rs_aux.ρl = rs.ρl
    rs_aux.dl = rs.dl
    rs_aux.N = rs.N
    rs_aux.ρn = rs.ρn
    rs_aux.dn = rs.dn
    rs_aux.cpg_pos = rs.cpg_pos
    rs_aux.nls_rgs = rs.nls_rgs

    # Return nothing
    return nothing

end
@everywhere function gen_grp_data(n_grp_reps::Int64, ϕ1::Vector{Float64}, ϕ2::Vector{Float64}, config::CpelNano.CpelNanoConfig)

    ## Set baseline struct
    rs = create_baseline_struct(config)
    αs1, βs1 = CpelNano.get_αβ_from_ϕ(ϕ1, rs)
    αs2, βs2 = CpelNano.get_αβ_from_ϕ(ϕ2, rs)

    ## Generate group data
    ms_g1 = Vector{CpelNano.RegStruct}()
    ms_g2 = Vector{CpelNano.RegStruct}()
    for rep = 1:n_grp_reps

        ## Group 1

        # Generate calls
        rs_aux = CpelNano.cpel_samp_ont(M, αs1, βs1, POBS, LEVEL_MEAN_M, LEVEL_STD_M, LEVEL_MEAN_U, LEVEL_STD_U)

        # Fill fields
        cp_rs_flds!(rs, rs_aux)

        # Push struct
        push!(ms_g1, rs_aux)

        ## Group 2

        # Generate data
        rs_aux = CpelNano.cpel_samp_ont(M, αs2, βs2, POBS, LEVEL_MEAN_M, LEVEL_STD_M, LEVEL_MEAN_U, LEVEL_STD_U)

        # Fill fields
        cp_rs_flds!(rs, rs_aux)

        # Push struct
        push!(ms_g2, rs_aux)

    end

    # Returns group data
    return ms_g1, ms_g2

end
@everywhere function proc_rep_rs!(rs::CpelNano.RegStruct, config::CpelNano.CpelNanoConfig)

    # Estimate parameter vector ϕ
    if length(rs.ϕhat) == 0
        CpelNano.get_ϕhat!(rs, config)
        length(rs.ϕhat) > 0 || return nothing
    end

    # Re-scale to per cpg site resolution
    CpelNano.rscle_grp_mod!(rs)

    # Get matrices
    CpelNano.get_rs_lgtrck_mats!(rs)

    # Partition function
    CpelNano.get_rs_logZ!(rs)

    # Get them all
    CpelNano.get_rs_log_gs!(rs)

    # Compute E[log g(X)]
    CpelNano.get_rs_exp_log_g1!(rs)
    CpelNano.get_rs_exp_log_g2!(rs)

    # Compute E[X] and E[XX]
    CpelNano.get_rs_exps!(rs)

    # Compute μ(X)
    CpelNano.comp_mml!(rs)

    # Compute h(X)
    CpelNano.comp_nme!(rs)

    # Set processed
    rs.proc = true

    # Return nothing
    return nothing

end
function get_pvals(test_out::CpelNano.RegStatTestStruct, eff_type::String)

    # Get
    if eff_type == "mml"
        pvals = [x[2] for x in test_out.tests.tmml_test]
    elseif eff_type == "nme"
        pvals = [x[2] for x in test_out.tests.tnme_test]
    elseif eff_type == "cmd"
        pvals = [x[2] for x in test_out.tests.tcmd_test]
    else
        CpelNano.print_log("Wrong effect type $(eff_type)")
    end

    # Return p-values
    return pvals

end
@everywhere function pmap_run_roc_sim(n_grp_reps::Int64, ϕ1::Vector{Float64}, ϕ2::Vector{Float64}, config::CpelNano.CpelNanoConfig)

    # Choose whether effect or not (50% chance)
    eff_bool_j = rand() > 0.5

    # Generate data
    if eff_bool_j
        # Introduce effect
        ms_g1, ms_g2 = gen_grp_data(n_grp_reps, ϕ1, ϕ2, config)
    else
        # Do not introduce effect
        ms_g1, ms_g2 = gen_grp_data(n_grp_reps, ϕ1, ϕ1, config)
    end

    # Estimate parameters
    for rep = 1:n_grp_reps
        proc_rep_rs!(ms_g1[rep], config)
        proc_rep_rs!(ms_g2[rep], config)
    end

    # Find valid reps
    n_val = 0
    for rep = 1:n_grp_reps
        rs_1 = ms_g1[rep]
        rs_2 = ms_g2[rep]
        n_val += (rs_1.proc && rs_2.proc) ? 1 : 0
    end
    n_val == n_grp_reps || return nothing

    # Perform test
    test_out = config.matched ? CpelNano.mat_est_reg_test(ms_g1, ms_g2) : CpelNano.unmat_est_reg_test(ms_g1, ms_g2)

    # Return test out and effect bool
    return eff_bool_j, test_out

end
@everywhere function pmap_run_null_sim(n_grp_reps::Int64, ϕ::Vector{Float64}, config::CpelNano.CpelNanoConfig)

    CpelNano.print_log("Producing set of p-values...")

    # Generate data
    ms_g1, ms_g2 = gen_grp_data(n_grp_reps, ϕ, ϕ, config)

    # Estimate parameters
    for rep = 1:n_grp_reps
        proc_rep_rs!(ms_g1[rep], config)
        proc_rep_rs!(ms_g2[rep], config)
    end

    # Find valid reps
    n_val = 0
    for rep = 1:n_grp_reps
        rs_1 = ms_g1[rep]
        rs_2 = ms_g2[rep]
        n_val += (rs_1.proc && rs_2.proc) ? 1 : 0
    end
    n_val == n_grp_reps || return nothing

    # Perform test
    test_out = config.matched ? CpelNano.mat_est_reg_test(ms_g1, ms_g2) : CpelNano.unmat_est_reg_test(ms_g1, ms_g2)

    # Return test out
    return test_out

end
function run_roc_sim(n_sim_reps::Int64, n_grp_reps::Int64, eff_type::String, matched::Bool)

    # Init
    pval_thrs_rng = 0.0:0.001:1.0
    tpr = Array{Float64}(undef, length(EF_PERC_RNG), length(pval_thrs_rng))
    fpr = Array{Float64}(undef, length(EF_PERC_RNG), length(pval_thrs_rng))

    # CpelNano config
    config = CpelNano.CpelNanoConfig()
    config.max_size_subreg = 350
    config.matched = matched
    config.max_em_iters = 25
    config.max_em_init = 5
    config.verbose = false

    # Get vectors
    eff_size_arr, ϕ1_arr, ϕ2_arr = get_ϕs_eff_rng(eff_type, config)

    # Loop over percent of effect size
    for i = 1:length(ϕ1_arr)

        # Init
        ϕ1 = ϕ1_arr[i]
        ϕ2 = ϕ2_arr[i]
        CpelNano.print_log("Effect percentage $(eff_size_arr[i])")

        # Generate data for multiple reps
        pmap_out = pmap(i -> pmap_run_roc_sim(n_grp_reps, ϕ1, ϕ2, config), 1:n_sim_reps)

        # Get pvals
        pvals_neg = []
        pvals_pos = []
        for x in pmap_out
            isnothing(x) && continue
            if x[1]
                push!(pvals_pos, get_pvals(x[2], eff_type))
            else
                push!(pvals_neg, get_pvals(x[2], eff_type))
            end
        end

        # Flatten vectors
        pvals_pos = vcat(pvals_pos...)
        pvals_neg = vcat(pvals_neg...)

        # Remove NaN
        is_nan = isnan.(pvals_pos)
        pvals_pos = pvals_pos[.!is_nan]
        is_nan = isnan.(pvals_neg)
        pvals_neg = pvals_neg[.!is_nan]

        # Store
        for (j, p_thrs) in enumerate(pval_thrs_rng)

            # Count true positives and true negatives
            if length(pvals_pos) > 0
                tp = sum(pvals_pos .<= p_thrs)
                fn = sum(pvals_pos .> p_thrs)
            else
                tp = NaN
                fn = NaN
            end
            if length(pvals_neg) > 0
                tn = sum(pvals_neg .> p_thrs)
                fp = sum(pvals_neg .<= p_thrs)
            else
                tn = NaN
                fp = NaN
            end

            # Store TPR & FPR
            tpr[i,j] = tp / (tp + fn + 1)
            fpr[i,j] = fp / (fp + tn + 1)

        end

    end

    return eff_size_arr, tpr, fpr

end
function run_null_sim(n_sim_reps::Int64, n_grp_reps::Int64, matched::Bool)

    # Init
    ϕ = fill(0.0, 3)

    # CpelNano config
    config = CpelNano.CpelNanoConfig()
    config.max_size_subreg = 350
    config.matched = matched
    config.max_em_iters = 25
    config.max_em_init = 5
    config.verbose = false

    # Generate data for multiple reps
    pmap_out = pmap(i -> pmap_run_null_sim(n_grp_reps, ϕ, config), 1:n_sim_reps)

    # Get pvals
    pvals_tmml = []
    pvals_tnme = []
    pvals_tcmd = []
    for x in pmap_out
        isnothing(x) && continue
        push!(pvals_tmml, get_pvals(x, "mml"))
        push!(pvals_tnme, get_pvals(x, "nme"))
        push!(pvals_tcmd, get_pvals(x, "cmd"))
    end

    # Flatten vectors
    pvals_tmml = vcat(pvals_tmml...)
    pvals_tnme = vcat(pvals_tnme...)
    pvals_tcmd = vcat(pvals_tcmd...)

    # Remove NaN
    is_nan = isnan.(pvals_tmml)
    pvals_tmml = pvals_tmml[.!is_nan]
    is_nan = isnan.(pvals_tnme)
    pvals_tnme = pvals_tnme[.!is_nan]
    is_nan = isnan.(pvals_tcmd)
    pvals_tcmd = pvals_tcmd[.!is_nan]

    # Return pvals
    return pvals_tmml, pvals_tnme, pvals_tcmd

end
function plt_roc(eff_size_arr, tpr, fpr, ttl)

    # Add each effect size
    plt = plot(xlabel="FPR", ylabel="TPR", title=ttl, xlim=(0, 1), ylim=(0, 1), size=(700, 700), legend=:bottomright)
    for (i, eff_size) in enumerate(eff_size_arr)
        plot!(plt, fpr[i,:], tpr[i,:],  seriestype=:line, label="$(round(eff_size;digits=2))", color=blind_friend_col[i])
    end

    # Return plot
    return plt

end
function plt_null_pvals_hist(pvals_tmml, pvals_tnme, pvals_tcmd, matched)

    # Vars
    y_lim = 0.5
    n_bins = 20
    y_label = "percentage (%)"

    # Tmml
    h = fit(Histogram, pvals_tmml, 0:1 / n_bins:1.0)
    midpts = collect(h.edges[1]) .- 2 / n_bins / n_bins
    percs = h.weights ./ sum(h.weights)
    plt_tmml = plot(midpts, percs,  seriestype=:bar, label="", xlabel="",
        ylabel=y_label, title="Tmml", xlim=(0, 1), ylim=y_lim, size=(600, 500))

    # Tnme
    h = fit(Histogram, pvals_tnme, 0:1 / n_bins:1.0)
    midpts = collect(h.edges[1]) .- 2 / n_bins
    percs = h.weights ./ sum(h.weights)
    if matched
        plt_tnme = plot(midpts, percs,  seriestype=:bar, label="", xlabel="p-value",
            ylabel=y_label, title="Tnme", xlim=(0, 1), ylim=y_lim, size=(600, 500))
    else
        plt_tnme = plot(midpts, percs,  seriestype=:bar, label="", xlabel="",
            ylabel=y_label, title="Tnme", xlim=(0, 1), ylim=y_lim, size=(600, 500))
    end

    # Tcmd
    h = fit(Histogram, pvals_tcmd, 0:1 / n_bins:1.0)
    midpts = collect(h.edges[1]) .- 2 / n_bins
    percs = h.weights ./ sum(h.weights)
    plt_tcmd = plot(midpts, percs,  seriestype=:bar, label="", xlabel="p-value",
        ylabel=y_label, title="Tcmd", xlim=(0, 1), ylim=y_lim, size=(600, 500))

    # Collage
    if matched
        plt = plot(plt_tmml, plt_tnme, layout=(2, 1), size=(600, 800))
    else
        plt = plot(plt_tmml, plt_tnme, plt_tcmd, layout=(3, 1), size=(600, 1200))
    end

    # Return plot
    return plt

end
function plt_null_pvals_ecdf(pvals_tmml, pvals_tnme, pvals_tcmd, matched)

    # Vars
    y_lim = (0, 1.0)
    x_lim = (0, 1.0)
    y_label = "ecdf"
    pval_rng = 0.01:0.01:1.0

    # Tmml
    F = ecdf(pvals_tmml)
    plt_tmml = plot(xlabel=" ", ylabel=y_label, title="Tmml", xlim=x_lim, ylim=y_lim, size=(600, 500))
    plot!(plt_tmml, pval_rng, pval_rng, seriestype=:line, label="Ideal")
    plot!(pval_rng, F.(pval_rng), seriestype=:line, label="Empirical")

    # Tnme
    F = ecdf(pvals_tnme)
    if matched
        plt_tnme = plot(xlabel="p-value", ylabel=y_label, title="Tnme", xlim=x_lim, ylim=y_lim, size=(600, 500))
        plot!(plt_tnme, pval_rng, pval_rng, seriestype=:line, label="Ideal")
        plot!(pval_rng, F.(pval_rng), seriestype=:line, label="Empirical")
    else
        plt_tnme = plot(xlabel=" ", ylabel=y_label, title="Tnme", xlim=x_lim, ylim=y_lim, size=(600, 500))
        plot!(plt_tnme, pval_rng, pval_rng, seriestype=:line, label="Ideal")
        plot!(pval_rng, F.(pval_rng), seriestype=:line, label="Empirical")
    end

    # Tcmd
    F = ecdf(pvals_tcmd)
    plt_tcmd = plot(xlabel="p-value", ylabel=y_label, title="Tcmd", xlim=x_lim, ylim=y_lim, size=(600, 500))
    plot!(plt_tcmd, pval_rng, pval_rng, seriestype=:line, label="Ideal")
    plot!(pval_rng, F.(pval_rng), seriestype=:line, label="Empirical")

    # Collage
    if matched
        plt = plot(plt_tmml, plt_tnme, layout=(2, 1), size=(600, 800))
    else
        plt = plot(plt_tmml, plt_tnme, plt_tcmd, layout=(3, 1), size=(600, 1200))
    end

    # Return plot
    return plt

end


text:
using Minerva

env = MountainCar()
s, done = reset!(env)
@info s, done
for step in 1:env.maxsteps
    a = rand(env.action_space)
    s′, r, done, info = step!(env, a)
    @info step, s′, r, done, info
    done && break
end


text:
module BibTeX

using DataStructures
using BibInternal

export parse_string

mutable struct Accumulator
    from::Int
    to::Int
end

struct Content
    # comments::Dict{Int, String}
    entries::OrderedDict{String,BibInternal.Entry}
    # free::Dict{Int, String}
    # preambles::Dict{Int, String}
    strings::Dict{String,String}
end

function Content()
    # comments = Dict{Int, String}()
    entries = OrderedDict{String,BibInternal.Entry}()
    # free = Dict{Int, String}()
    # preambles = Dict{Int, String}()
    strings = Dict{String,String}([
        "jan" => "January",
        "feb" => "February",
        "mar" => "March",
        "apr" => "April",
        "may" => "May",
        "jun" => "June",
        "jul" => "July",
        "aug" => "August",
        "sep" => "September",
        "oct" => "October",
        "nov" => "November",
        "dec" => "December",
    ])
    return Content(entries, strings)
end

mutable struct Field
    braces::Int
    name::String
    quotes::Bool
    value::String

    Field() = new(0, "", false, "")
end

mutable struct Position
    row::Int
    col::Int
end

Position() = Position(0, 0)

mutable struct Storage
    delim::Union{Char,Nothing}
    fields::Vector{Field}
    key::String
    kind::String
end

Storage() = Storage(nothing, Vector{Field}(), "", "")

function make_entry(storage)
    # @info "making entry" storage
    d = Dict("_type" => storage.kind)
    foreach(field -> push!(d, field.name => field.value), storage.fields)
    return d
end

mutable struct Parser
    acc::Accumulator
    content::Content
    field::Field
    input::Vector{Char}
    pos_start::Position
    pos_end::Position
    storage::Storage
    task::Symbol

    function Parser(input;
        acc=Accumulator(1, 0),
        content=Content(),
        field=Field(),
        pos_start=Position(1, 1),
        pos_end=Position(1, 0),
        storage=Storage(),
        task=:free,
    )
        new(acc, content, field, collect(input), pos_start, pos_end, storage, task)
    end
end

rev(char) = char == '(' ? ')' : '}'

get_entries(parser) = parser.content.entries
function get_acc(parser; from = 1, to = 0)
    a = from + parser.acc.from - 1
    b = parser.acc.to - to
    return prod(parser.input[a:b])
end

set_delim!(parser, char) = parser.storage.delim = char
set_entry_kind!(parser, kind) = parser.storage.kind = kind

inc_col!(t) = t.pos_end.col += 1
function inc_row!(t)
    t.pos_end.row += 1
    t.pos_end.col = 1
end

function inc!(t, char, dumped)
    char == '\n' ? inc_row!(t) : inc_col!(t)
    t.acc.to += 1
    if dumped
        t.pos_start = t.pos_end
        t.acc.from = t.acc.to
    end
end

is_dumped(parser, char, ::Val{:free}) = char == '@'
dump!(parser, char, ::Val) = char == '@' && (parser.task = :entry)

is_dumped(parser, char, ::Val{:entry}) = occursin(r"[@{\(\n]", char)
function dump!(parser, char, ::Val{:entry})
    if char == '\n'
        parser.task = :free
    elseif char ∈ ['{', '(']
        set_delim!(parser, char)
        acc = split(lowercase(get_acc(parser; from = 2)), r"[\t ]+")
        if length(acc) ≤ 2 && !isempty(acc[1])
            set_entry_kind!(parser, acc[1])
            # parser.task = acc[1] ∈ ["comment", "preamble", "string"] ? Symbol(acc[1]) : :key
            parser.task = acc[1] ∈ ["string"] ? Symbol(acc[1]) : :key
        else
            parser.task = :free
        end
    end
end

is_dumped(parser, char, ::Val{:key}) = char ∈ ['@', ',']
function dump!(parser, char, ::Val{:key})
    if char == '@'
        parser.task = :entry
    elseif char == ','
        acc = split(get_acc(parser; from = 2), r"[\t ]+")
        if length(acc) ≤ 2 && !isempty(acc[1])
            parser.storage.key = acc[1]
            parser.task = :field_name
        else
            parser.task = :free
        end
    end
end

is_dumped(parser, char, ::Val{:field_name}) = char ∈ ['=', '@']
function dump!(parser, char, ::Val{:field_name})
    if char == '@'
        parser.task = :entry
    elseif char == '='
        acc = split(get_acc(parser; from = 2), r"[\t\r\n ]+"; keepempty=false)
        if length(acc) == 1
            parser.field.name = acc[1]
            parser.task = :field_in
        else
            parser.task = :free
        end
    end
end

is_dumped(parser, char, ::Val{:field_in}) = occursin(r"[0-9@a-zA-Z\"{]", char)
function dump!(parser, char, ::Val{:field_in})
    if char == '@'
        parser.task = :entry
    elseif char == '"'
        parser.pos_start.col += 1
        parser.acc.from += 1
        parser.task = :field_inquote
    elseif char == '{'
        parser.task = :field_inbrace
    elseif occursin(r"[0-9]", char)
        parser.task = :field_number
    elseif occursin(r"[a-zA-Z]", char)
        parser.task = :field_var
    end
end

is_dumped(parser, char, ::Val{:field_inbrace}) = char ∈ ['@', '{', '}']
function dump!(parser, char, ::Val{:field_inbrace})
    if char == '@'
        parser.task = :entry
    elseif char == '}' && parser.field.braces == 0
        parser.field.value *= get_acc(parser; from = 2)
        push!(parser.storage.fields, deepcopy(parser.field))
        parser.field.value = ""
        parser.task = :field_out
    elseif char == '{'
        parser.field.value *= get_acc(parser; from = 2) # * parser.field.value
        parser.field.braces += 1
    elseif char == '}'
        parser.field.value *= get_acc(parser; from = 1, to = -1) # * parser.field.value
        parser.field.braces -= 1
    end
end

is_dumped(parser, char, ::Val{:field_inquote}) = char ∈ ['{', '"', '}']
function dump!(parser, char, ::Val{:field_inquote})
    if char == '"' && parser.field.braces == 0
        parser.field.value *= get_acc(parser; from = 2)
        parser.task = :field_outquote
    elseif char == '{'
        parser.field.value *= get_acc(parser; from = 2) # * parser.field.value
        parser.field.braces += 1
    elseif char == '}'
        parser.field.value *= get_acc(parser; from = 1, to = -1) # * parser.field.value
        parser.field.braces -= 1
    end
end

function is_dumped(parser, char, ::Val{:field_outquote})
    return char ∈ ['@', ',', '#', rev(parser.storage.delim)]
end
function dump!(parser, char, ::Val{:field_outquote})
    if char == '@'
        parser.task = :entry
    elseif char == '#'
        parser.task = :field_concat
    else
        push!(parser.storage.fields, deepcopy(parser.field))
        parser.field.value = ""
        if char == ','
            parser.task = :field_next
        elseif char == rev(parser.storage.delim)
            entry = make_entry(parser.storage)
            push!(parser.content.entries,
            parser.storage.key => BibInternal.make_bibtex_entry(parser.storage.key, entry)
            )
            parser.task = :free
        end
    end
end

is_dumped(parser, char, ::Val{:field_concat}) = occursin(r"[a-zA-Z\"@]", char)
function dump!(parser, char, ::Val{:field_concat})
    if char == '@'
        parser.task = :entry
    elseif char == '"'
        parser.pos_start.col += 1
        parser.acc.from += 1
        parser.task = :field_inquote
    elseif occursin(r"[a-zA-Z]", char)
        parser.task = :field_var
    end
end

function is_dumped(parser, char, ::Val{:field_var})
    return char ∈ ['@', ',', '#', rev(parser.storage.delim)]
end
function dump!(parser, char, ::Val{:field_var})
    if char == '@'
        parser.task = :entry
    else
        acc = split(get_acc(parser), r"[\t\r\n ]+"; keepempty=false)
        if length(acc) == 1
            # @info "Printing" parser.content.strings[acc[1]] parser.field.value
            parser.field.value *= parser.content.strings[acc[1]]
            if char == '#'
                parser.task = :field_concat
            else
                push!(parser.storage.fields, deepcopy(parser.field))
                parser.field.value = ""
                if char == ','
                    parser.task = :field_next
                elseif char == rev(parser.storage.delim)
                    entry = make_entry(parser.storage)
                    push!(parser.content.entries, parser.storage.key =>
                        BibInternal.make_bibtex_entry(parser.storage.key, entry)
                    )
                    parser.task = :free
                end
            end
        else
            parser.task = :free
        end
    end
end

is_dumped(parser, char, ::Val{:field_number}) = char ∈ ['@', ',',rev(parser.storage.delim)]
function dump!(parser, char, ::Val{:field_number})
    if char == '@'
        parser.task = :entry
    else
        acc = split(get_acc(parser), r"[\t\r\n ]+"; keepempty=false)
        # @show acc
        if length(acc) == 1
            parser.field.value = acc[1]
            push!(parser.storage.fields, deepcopy(parser.field))
            parser.field.value = ""
            if char == ','
                parser.task = :field_next
            elseif char == rev(parser.storage.delim)
                entry = make_entry(parser.storage)
                push!(parser.content.entries,
                parser.storage.key => BibInternal.make_bibtex_entry(parser.storage.key, entry)
                )
                parser.task = :free
            end
        else
            parser.task = :free
        end
    end
end

is_dumped(parser, char, ::Val{:field_out}) = char ∈ ['@', ',', rev(parser.storage.delim)]
function dump!(parser, char, ::Val{:field_out})
    if char == '@'
        parser.task = :entry
    elseif char == ','
        parser.task = :field_next
    elseif char == rev(parser.storage.delim)
        entry = make_entry(parser.storage)
        # @show entry
        push!(parser.content.entries,
            parser.storage.key => BibInternal.make_bibtex_entry(parser.storage.key, entry)
        )
        parser.task = :free
    end
end

is_dumped(parser, char, ::Val{:field_next}) = char ∈ ['=', '@', rev(parser.storage.delim)]
function dump!(parser, char, ::Val{:field_next})
    if char == '@'
        parser.task = :entry
    elseif char == '='
        acc = split(get_acc(parser; from = 2), r"[\t\r\n ]+"; keepempty=false)
        # @show acc
        if length(acc) == 1
            parser.field.name = acc[1]
            parser.task = :field_in
        else
            parser.task = :free
        end
    elseif char == rev(parser.storage.delim)
        entry = make_entry(parser.storage)
        push!(parser.content.entries,
            parser.storage.key => BibInternal.make_bibtex_entry(parser.storage.key, entry)
        )
        parser.task = :free
    end
end

is_dumped(parser, char, ::Val{:string}) = char ∈ ['=', '@']
function dump!(parser, char, ::Val{:string})
    if char == '@'
        parser.task = :entry
    elseif char == '='
        acc = split(get_acc(parser; from = 2), r"[\t\r\n ]+"; keepempty=false)
        if length(acc) == 1
            parser.field.name = acc[1]
            parser.task = :string_inquote
        else
            parser.task = :free
        end
    else
        parser.task = :free
    end
end

is_dumped(parser, char, ::Val{:string_inquote}) = char ∈ ['"', '@']
function dump!(parser, char, ::Val{:string_inquote})
    if char == '@'
        parser.task = :entry
    elseif char == '"'
        parser.task = :string_value
    end
end

is_dumped(parser, char, ::Val{:string_value}) = char ∈ ['"']
function dump!(parser, char, ::Val{:string_value})
    if char == '"'
        parser.field.value = get_acc(parser; from = 2)
        parser.task = :string_outquote
    end
end

is_dumped(parser, char, ::Val{:string_outquote}) = char ∈ [rev(parser.storage.delim), '@']
function dump!(parser, char, ::Val{:string_outquote})
    if char == '@'
        parser.task = :entry
    elseif char == rev(parser.storage.delim)
        parser.content.strings[parser.field.name] = parser.field.value
        parser.task = :free
    end
end

is_dumped(parser, char) = is_dumped(parser, char, Val(parser.task))

function dump!(parser, char=' ')
    dump!(parser, char, Val(parser.task))
    parser.pos_start = parser.pos_end
    parser.acc.from = parser.acc.to
end

function parse!(parser, char)
    dumped = is_dumped(parser, char)
    dumped && dump!(parser, char)
    inc!(parser, char, dumped)
    # fieldout = parser.task == :field_out
    # dumped && @info parser parser char dumped
    # dumped && @warn parser.acc parser.content parser.field parser.input parser.pos_start parser.pos_end parser.storage parser.task
    # @info parser parser char dumped
    # @warn parser.acc parser.content parser.field parser.input parser.pos_start parser.pos_end parser.storage parser.task
end

function parse_string(str, ::Val{:bibtex})
    parser = Parser(str)
    foreach(char -> parse!(parser, char), parser.input)
    # @info "Dev: " parser
    return get_entries(parser)
end

function parse_file(path::String)
    entries = parse_string(read(path, String), Val(:bibtex))
    return entries
end

end # module


text:
# module PushInterface.DisplayInterface

# The display uses 16 bit per pixels in a 5:6:5 format
#     MSB                             LSB
#       b b b b|b g g g|g g g r|r r r r

# reference: https://github.com/Ableton/push2-display-with-juce/blob/master/Source/push2/Push2-Bitmap.h#L75
function rgb_to_pixel(x, alpha::UInt8, r::UInt8, g::UInt8, b::UInt8)::Tuple{UInt8,UInt8}
    if alpha == 0x00
        r, g, b = UInt8.((0, 0, 0))
    else
        r, g, b = (x -> UInt8(UInt16(x) * 0xff / UInt16(alpha))).((r, g, b))
    end
    pixel = (UInt16(b) & 0xf8) >> 3
    pixel <<= 6
    pixel += (UInt16(g) & 0xFC) >> 2
    pixel <<= 5;
    pixel += (UInt16(r) & 0xF8) >> 3
    pixel = xor(iseven(x) ? 0xffe7 : 0xf3e7, pixel)
    h = UInt8(pixel & 0xff)
    l = UInt8(pixel >> 8)
    (h, l)
end

const FillerBytes = fill(0x00, 128)

function framebuffer(canvas::Canvas)
    data = ccall((:cairo_image_surface_get_data, _jl_libcairo), Ptr{Cuchar}, (Ptr{Cvoid},), canvas.surface.ptr)
    stride = ccall((:cairo_image_surface_get_stride, _jl_libcairo), Cint, (Ptr{Cvoid},), canvas.surface.ptr) # 3840
    componentsize = Int(stride/width(canvas.surface)) # a r g b
    pixeldata = reshape(unsafe_wrap(Array, data, stride*DisplayPixelHeight), (componentsize, DisplayPixelWidth, DisplayPixelHeight))
    buf = UInt8[]
    for y in 1:DisplayPixelHeight
        for x in 1:DisplayPixelWidth
            argb = reverse(pixeldata[:, x, y])
            (h, l) = rgb_to_pixel(x, argb...)
            push!(buf, h, l)
        end
        push!(buf, FillerBytes...)
    end
    buf
end

#import Cairo: write_to_png
#write_to_png(canvas.surface, "sample_line.png")

# module PushInterface.DisplayInterface


text:
module ABCDeGraphGenerator

using Random
using StatsBase
using Base.Threads
using SentinelArrays

include("pl_sampler.jl")
include("graph_sampler.jl")

end # module


text:
function [x,err,it] = perform_conjugate_gradient(A,y)

    """
     perform_conjugate_gradient - perform (bi)-conjugate gradient

         [x,err,k] = perform_conjugate_gradient(A,y,options);

       Solves for A*x=y.
       Works both for vector x,y and matrix x,y (parralel soving of many
       linear equations).

       Important: the algorithm assumes that the matrix is symmetric definite
       positive. If it is not the case, then you must set options.is_sdp=0,
       and the algorithm will use a bi-conjugate gradient descent (x2 slower).

       A can be a matrix or a callback function y=A(x,options).
       In this case, you have to set options.ncols as the number of columns of
       A (if it is a callback).
       In the non-symmetric case, then options.transpose=1 for the computation
       of A*x and options.transpose=-1 for the computation of A'*x.

       err monitors the decreasing of |A*x-y|.
       k is the total number of iterations.

       You can set:
           options.x is an initial guess
           options.epsilon is maximum error
           options.niter_max is maximum number of iterations

       Copyright (c) 2007 Gabriel Peyre
    """

options.null = 0;
niter = 100
epsilon = getoptions(options, 'epsilon', 1e-5);
is_sdp = getoptions(options, 'is_sdp', 1);
x = getoptions(options, 'x', []);

use_callback = 0;
if not(isnumeric(A))
    use_callback = 1;
end
if not(isempty(x))
    x = options.x;
else
    if use_callback==0
        x = zeros(size(A,2),1);
    else
        if isfield(options, 'ncols')
            x = zeros(options.ncols, 1);
        else
            error('You have to specify options.ncols');
        end
    end
end

if norm(y(:), 'fro') == 0
    normb = epsilon;
else
    normb = epsilon; % * sum(y(:).^2);
end

if use_callback==0
    r  = y - A*x;
else
    options.transpose = 1;
    r  = y - feval(A,x,options);
end
p = r;
if is_sdp==0
    rr = r;
    pp = p;
end
r0 = sum(r(:).^2);

err = [sum(r0)];
for it=1:niter

    % auxiliary vector
    if use_callback==0
        w  = A*p;
    else
        options.transpose = 1;
        % options.transpose = -1;
        w  = feval(A,p,options);
    end

    if is_sdp==1
        d = sum(p(:) .* w(:));
    else
        d = sum(pp(:) .* w(:));
    end
    I = find(abs(d)<eps); d(I) = 1;
    alpha = repmat( r0 ./ d, [size(x,1) 1] );              % found optimal alpha in line search
    x = x + alpha.*p;                       % new guess
    r = r - alpha.*w;                       % the residual is in fact r=b-A*x
    if is_sdp==0
        if use_callback==0
            w  = A'*pp;
        else
            options.transpose = -1;
            w  = feval(A,pp,options);
        end
        rr = rr-alpha.*w;
    end

    rprev = r0;                             % save norm of the old residual
    if is_sdp==0
        r0 = sum(r.*rr);
    else
        r0 = sum(r.^2);                         % compute norm of new residual
    end

    err(end+1) = sqrt( sum(r.^2) );
    if err(end)<normb
        return;
    end

    % search direction
    beta = r0./rprev;
    p = r + repmat(beta, [size(x,1) 1]).*p;
    if is_sdp==0
        pp = rr + repmat(beta, [size(x,1) 1]).*pp;
    end

end


text:
function LnsCoefficients(arg0::jint, arg1::jint, arg2::Vector{Vector{jdouble}}, arg3::TreeMap, arg4::jdouble)
    return LnsCoefficients((jint, jint, Vector{Vector{jdouble}}, TreeMap, jdouble), arg0, arg1, arg2, arg3, arg4)
end

function get_lns(obj::LnsCoefficients, arg0::jint, arg1::jint)
    return jcall(obj, "getLns", jdouble, (jint, jint), arg0, arg1)
end

function getd_lnsd_gamma(obj::LnsCoefficients, arg0::jint, arg1::jint)
    return jcall(obj, "getdLnsdGamma", jdouble, (jint, jint), arg0, arg1)
end


text:
using Boltzmann
using Mocha
using HDF5

const DBN_PATH = "/tmp/dbn_params.hdf5"

# prepare training data
X = rand(100, 2000)
y = rand(2000)

# construct and train DBN
dbn_layers = [("vis", GRBM(100, 50)),
              ("hid1", BernoulliRBM(50, 25)),
              ("hid2", BernoulliRBM(25, 20))]
dbn = DBN(dbn_layers)
fit(dbn, X)

# save learn parameters (weights and biases)
if isfile(DBN_PATH)
    rm(DBN_PATH)
end
save_params(DBN_PATH, dbn)

# prepare Mocha network with similar architecture, but with additional layers
# for data and softmax classification
# names of layers should be the same as in DBN
backend = CPUBackend()
data = MemoryDataLayer(tops=[:data, :label], batch_size=500, data=Array[X, y])
vis = InnerProductLayer(name="vis", output_dim=50, tops=[:vis], bottoms=[:data])
hid1 = InnerProductLayer(name="hid1", output_dim=25, tops=[:hid1], bottoms=[:vis])
hid2 = InnerProductLayer(name="hid2", output_dim=20, tops=[:hid2], bottoms=[:hid1])
loss = SoftmaxLossLayer(name="loss",bottoms=[:hid2, :label])

net = Net("TEST", backend, [data, vis, hid1, hid2])

# finally, load pretrained parameters into corresponding layers of Mocha network
h5open(DBN_PATH) do h5
    load_network(h5, net)
    # if architecture of Mocha network doesn't match DBN exactly, one can
    # additonally pass `die_if_not_found=false` to ignore missing layers and
    # use default initializers instead
    #  load_network(h5, net, die_if_not_found=false)
end

rm(DBN_PATH)


text:
using LinearAlgebra
function vech(x)
    k = size(x,1)
    a = zeros(Int((k^2-k)/2 + k))
    m = 1
    for i = 1:k
        for j = 1:i
            a[m] = x[i,j]
            m += 1
        end
    end
    a
end

# for testing vech()
function vech()
    a = eye(3)
    b = vech(a)
    c = [1.0; 0.0; 1.0; 0.0; 0.0; 1.0]
    b == c
end

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


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

mutable struct VirtualMachineScaleSetPublicIPAddressConfigurationProperties <: SwaggerModel
    idleTimeoutInMinutes::Any # spec type: Union{ Nothing, Int32 } # spec name: idleTimeoutInMinutes
    dnsSettings::Any # spec type: Union{ Nothing, VirtualMachineScaleSetPublicIPAddressConfigurationDnsSettings } # spec name: dnsSettings

    function VirtualMachineScaleSetPublicIPAddressConfigurationProperties(;idleTimeoutInMinutes=nothing, dnsSettings=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetPublicIPAddressConfigurationProperties, Symbol("idleTimeoutInMinutes"), idleTimeoutInMinutes)
        setfield!(o, Symbol("idleTimeoutInMinutes"), idleTimeoutInMinutes)
        validate_property(VirtualMachineScaleSetPublicIPAddressConfigurationProperties, Symbol("dnsSettings"), dnsSettings)
        setfield!(o, Symbol("dnsSettings"), dnsSettings)
        o
    end
end # type VirtualMachineScaleSetPublicIPAddressConfigurationProperties

const _property_map_VirtualMachineScaleSetPublicIPAddressConfigurationProperties = Dict{Symbol,Symbol}(Symbol("idleTimeoutInMinutes")=>Symbol("idleTimeoutInMinutes"), Symbol("dnsSettings")=>Symbol("dnsSettings"))
const _property_types_VirtualMachineScaleSetPublicIPAddressConfigurationProperties = Dict{Symbol,String}(Symbol("idleTimeoutInMinutes")=>"Int32", Symbol("dnsSettings")=>"VirtualMachineScaleSetPublicIPAddressConfigurationDnsSettings")
Base.propertynames(::Type{ VirtualMachineScaleSetPublicIPAddressConfigurationProperties }) = collect(keys(_property_map_VirtualMachineScaleSetPublicIPAddressConfigurationProperties))
Swagger.property_type(::Type{ VirtualMachineScaleSetPublicIPAddressConfigurationProperties }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualMachineScaleSetPublicIPAddressConfigurationProperties[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetPublicIPAddressConfigurationProperties }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetPublicIPAddressConfigurationProperties[property_name]

function check_required(o::VirtualMachineScaleSetPublicIPAddressConfigurationProperties)
    true
end

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


text:
# Chemfiles.jl, a modern library for chemistry file reading and writing
# Copyright (C) Guillaume Fraux and contributors -- BSD license

export bonds_count, angles_count, dihedrals_count, impropers_count
export bonds, angles, dihedrals, impropers
export add_bond!, remove_bond!, bond_order, bond_orders
export add_residue!, residues, are_linked, count_residues

__ptr(topology::Topology) = __ptr(topology.__handle)
__const_ptr(topology::Topology) = __const_ptr(topology.__handle)

"""
Possible bond orders in Chemfiles:
    - ``Chemfiles.UnknownBond``: when the bond order is not specified
    - ``Chemfiles.SingleBond``: for single bonds
    - ``Chemfiles.DoubleBond``: for double bonds
    - ``Chemfiles.TripleBond``: for triple bonds
    - ``Chemfiles.QuadrupleBond``: for quadruple bonds (present in some metals)
    - ``Chemfiles.QintupletBond``: for qintuplet bonds (present in some metals)
    - ``Chemfiles.AmideBond``: for amide bonds
    - ``Chemfiles.AromaticBond``: for aromatic bonds
"""
@enum BondOrder begin
    UnknownBond = lib.CHFL_BOND_UNKNOWN
    SingleBond = lib.CHFL_BOND_SINGLE
    DoubleBond = lib.CHFL_BOND_DOUBLE
    TripleBond = lib.CHFL_BOND_TRIPLE
    QuadrupleBond = lib.CHFL_BOND_QUADRUPLE
    QintupletBond = lib.CHFL_BOND_QINTUPLET
    AmideBond = lib.CHFL_BOND_AMIDE
    AromaticBond = lib.CHFL_BOND_AROMATIC
end

"""
Create an empty ``Topology``.
"""
function Topology()
    ptr = @__check_ptr(lib.chfl_topology())
    return Topology(CxxPointer(ptr, is_const=false))
end

"""
Get a copy of the ``Topology`` of the given ``frame``.
"""
function Topology(frame::Frame)
    ptr = @__check_ptr(lib.chfl_topology_from_frame(__const_ptr(frame)))
    topology = Topology(CxxPointer(ptr, is_const=true))
    copy = deepcopy(topology)
    finalize(topology)
    return copy
end

"""
Get the ``Topology`` size, i.e. the current number of atoms.
"""
function Base.size(topology::Topology)
    count = Ref{UInt64}(0)
    __check(lib.chfl_topology_atoms_count(
        __const_ptr(topology), count
    ))
    return count[]
end

"""
Add an ``atom`` at the end of a ``topology``.
"""
function add_atom!(topology::Topology, atom::Atom)
    __check(lib.chfl_topology_add_atom(__ptr(topology), __const_ptr(atom)))
    return nothing
end

"""
Remove the atom at the given ``index`` from a ``topology``.
"""
function remove_atom!(topology::Topology, index::Integer)
    __check(lib.chfl_topology_remove(__ptr(topology), UInt64(index)))
    return nothing
end

"""
Get the number of bonds in the ``topology``.
"""
function bonds_count(topology::Topology)
    count = Ref{UInt64}(0)
    __check(lib.chfl_topology_bonds_count(__const_ptr(topology), count))
    return count[]
end

"""
Get the number of angles in the ``topology``.
"""
function angles_count(topology::Topology)
    count = Ref{UInt64}(0)
    __check(lib.chfl_topology_angles_count(__const_ptr(topology), count))
    return count[]
end

"""
Get the number of dihedral angles in the ``topology``.
"""
function dihedrals_count(topology::Topology)
    count = Ref{UInt64}(0)
    __check(lib.chfl_topology_dihedrals_count(__const_ptr(topology), count))
    return count[]
end

"""
Get the number of improper angles in the ``topology``.
"""
function impropers_count(topology::Topology)
    count = Ref{UInt64}(0)
    __check(lib.chfl_topology_impropers_count(__const_ptr(topology), count))
    return count[]
end

"""
Get the bonds in the ``topology``, in a ``2 x bonds_count(topology)`` array.
"""
function bonds(topology::Topology)
    count = bonds_count(topology)
    result = Array{UInt64}(undef, 2, count)
    __check(lib.chfl_topology_bonds(__const_ptr(topology), pointer(result), count))
    return result
end

"""
Get the angles in the ``topology``, in a ``3 x angles_count(topology)`` array.
"""
function angles(topology::Topology)
    count = angles_count(topology)
    result = Array{UInt64}(undef, 3, count)
    __check(lib.chfl_topology_angles(__const_ptr(topology), pointer(result), count))
    return result
end

"""
Get the dihedral angles in the ``topology``, in a ``4 x dihedrals_count(topology)``
array.
"""
function dihedrals(topology::Topology)
    count = dihedrals_count(topology)
    result = Array{UInt64}(undef, 4, count)
    __check(lib.chfl_topology_dihedrals(__const_ptr(topology), pointer(result), count))
    return result
end

"""
Get the improper angles in the ``topology``, in a ``4 x impropers_count(topology)``
array.
"""
function impropers(topology::Topology)
    count = impropers_count(topology)
    result = Array{UInt64}(undef, 4, count)
    __check(lib.chfl_topology_impropers(__const_ptr(topology), pointer(result), count))
    return result
end

"""
Add a bond between the atoms ``i`` and ``j`` in the ``topology``, optionaly
setting the bond ``order``.
"""
function add_bond!(topology::Topology, i::Integer, j::Integer, order=nothing)
    if order == nothing
        __check(lib.chfl_topology_add_bond(__ptr(topology), UInt64(i), UInt64(j)))
    else
        # Check that the order is a valid BondOrder
        order = BondOrder(Integer(order))
        __check(lib.chfl_topology_bond_with_order(
            __ptr(topology), UInt64(i), UInt64(j), lib.chfl_bond_order(order)
        ))
    end
    return nothing
end

"""
Remove any existing bond between the atoms ``i`` and ``j`` in the ``topology``.
"""
function remove_bond!(topology::Topology, i::Integer, j::Integer)
    __check(lib.chfl_topology_remove_bond(
        __ptr(topology), UInt64(i), UInt64(j)
    ))
    return nothing
end

"""
Get the ``BondOrder`` for the bond between atoms ``i`` and ``j`` in the
``topology``.
"""
function bond_order(topology::Topology, i::Integer, j::Integer)
    order = Ref{lib.chfl_bond_order}(lib.CHFL_BOND_UNKNOWN)
    __check(lib.chfl_topology_bond_order(
        __const_ptr(topology), UInt64(i), UInt64(j), order
    ))
    return BondOrder(order[])
end

"""
Get the ``BondOrder`` for all the bonds in the ``topology``.
"""
function bond_orders(topology::Topology)
    count = bonds_count(topology)
    result = Array{lib.chfl_bond_order}(undef, count)
    __check(lib.chfl_topology_bond_orders(
        __const_ptr(topology), pointer(result), count
    ))
    return map(BondOrder, result)
end

"""
Add a copy of ``residue`` to this ``topology``.

The residue id must not already be in the topology, and the residue must
contain only atoms that are not already in another residue.
"""
function add_residue!(topology::Topology, residue::Residue)
    __check(lib.chfl_topology_add_residue(
        __ptr(topology), __const_ptr(residue)
    ))
    return nothing
end

"""
Get the number of residues in the ``topology``.
"""
function count_residues(topology::Topology)
    count = Ref{UInt64}(0)
    __check(lib.chfl_topology_residues_count(
        __const_ptr(topology), count
    ))
    return count[]
end

"""
Check if the two residues ``first`` and ``second`` from the ``topology`` are
linked together, *i.e.* if there is a bond between one atom in the first
residue and one atom in the second one.
"""
function are_linked(topology::Topology, first::Residue, second::Residue)
    result = Ref{UInt8}(0)
    __check(lib.chfl_topology_residues_linked(
        __const_ptr(topology), __const_ptr(first), __const_ptr(second), result
    ))
    return convert(Bool, result[])
end

"""
Resize the ``topology`` to hold ``natoms`` atoms. If the new number of atoms is
bigger than the current number, new atoms will be created with an empty name
and type.

If it is lower than the current number of atoms, the last atoms will be removed,
together with the associated bonds, angles and dihedrals.
"""
function Base.resize!(topology::Topology, size::Integer)
    __check(lib.chfl_topology_resize(
        __ptr(topology), UInt64(size)
    ))
end

"""
Make a deep copy of a ``topology``.
"""
function Base.deepcopy(topology::Topology)
    ptr = lib.chfl_topology_copy(__const_ptr(topology))
    return Topology(CxxPointer(ptr, is_const=false))
end


text:
# matrices of Gensys
G0 = zeros(24,24)
G1 = zeros(24,24)
Psi = zeros(24,8)
Pi = zeros(24,6)

# Equation 1 (Euler)

G0[1,1] = -1
G0[1,2] = -(1-h)/sig
G0[1,3] = (1-h)/sig
G0[1,7] = 1+h
G0[1,8] = (1-h)/sig

G1[1,3] = (1+h)/sig
G1[1,7] = h

Pi[1,1] = -1
Pi[1,2] = -(1-h)/sig
Pi[1,3] = (1-h)/sig

# Equation 2 (goods market equilibrium)

G0[2,9]  = 1
G0[2,7]  = -(1-alfa)
G0[2,10] = -alfa*eta*(2-alfa)
G0[2,11] = -alfa*eta
G0[2,17] = -alfa

# Equation 3 (tot)

G0[3,10] = 1
G0[3,13] = -1
G0[3,12] = 1

G1[3,10] = 1

# Equation 4 (relationship q and tot)

G0[4,10] = -(1-alfa)
G0[4,14] = 1
G0[4,11] = -1

# Equation 5 (relationship q, s, pi)

G0[5,14] = 1
G0[5,15] = -1
G0[5,18] = -1

G1[5,2]  = -1
G1[5,14] = 1
G1[5,15] = -1

# Equation 6 (phillips curve - domestic)

G0[6,5]  = -bet
G0[6,12] = 1+bet*delta_D
G0[6,24] = -bet*(1-theta_D)*(1-theta_D*bet)/theta_D

G1[6,12] = delta_D

Pi[6,4] = -bet

# Equation 7 (marginal cost)

G0[7,24] = 1
G0[7,10] = -alfa
G0[7,9]  = -phi
G0[7,7]  = -sig/(1-h)
G0[7,20] = 1+phi

G1[7,7] = -(sig*h)/(1-h)

# Equation 8  (phillips curve - importing)

G0[8,6]  = -bet
G0[8,13] = 1+bet*delta_I
G0[8,11] = -(1-theta_I)*(1-theta_I*bet)/theta_I
G0[8,22] = -1

G1[8,13] = delta_I

Pi[8,5]  = -bet

# Equation 9 (relationship pi and tot)

G0[9,12] = -1
G0[9,10] = -alfa

G1[9,10] = -alfa
G1[9,2]  = -1

# Equation 10 (bugdet constraint)

G0[10,21] = 1
G0[10,10] = alfa
G0[10,11] = alfa
G0[10,9]  = -1
G0[10,7]  = 1

G1[10,21] = 1/bet

# Equation 11 (UIP)

G0[11,4]  = -1
G0[11,15] = 1
G0[11,8]  = 1
G0[11,19] = -1
G0[11,21] = chi
G0[11,23] = -1

Pi[11,6] = -1

# Equation 12 (Taylor RMI)

G0[12,8]  = 1
G0[12,16] = -rho_1
G0[12,9]  = -(1-rho_1-rho_2)*lambda_y
G0[12,15] = -(1-rho_1-rho_2)*lambda_s

G1[12,16] = rho_2
G1[12,2]  = (1-rho_1-rho_2)*lambda_pi
G1[12,15] = -(1-rho_1-rho_2)*lambda_s

Psi[12,5] = sigma_i

###################
## Foreign Block ##
###################

# System (Equations 13, 14 and 15)

G0[13,17] = 1
G0[14,17] = a_0piy
G0[14,18] = 1
G0[15,17] = a_0iy
G0[15,18] = a_0ipi
G0[15,19] = 1

G1[13,17] = a_1yy
G1[13,18] = a_1ypi
G1[13,19] = a_1yi
G1[14,17] = a_1piy
G1[14,18] = a_1pipi
G1[14,19] = a_1pii
G1[15,17] = a_1iy
G1[15,18] = a_1ipi
G1[15,19] = a_1ii

Psi[13,6] = sigma_y_star
Psi[14,7] = sigma_pi_star
Psi[15,8] = sigma_i_star

#################
## Shocks block #
#################

# Equation 16 (a)

G0[16,20] = 1

G1[16,20] = rho_a

Psi[16,1] = sigma_a

# Equation 17 (gamma)

G0[17,3] = 1

G1[17,3] = rho_gamma

Psi[17,2] = sigma_gamma

# Equation 18 (epsilon_cp)

G0[18,22] = 1

G1[18,22] = rho_cp

Psi[18,3] = sigma_cp

# Equation 19 (phi)

G0[19,23] = 1

G1[19,23] = rho_phi

Psi[19,4] = sigma_phi

####################
## Identity block ##
####################

# Equation 20 (c)

G0[20,7] = 1
G1[20,1] = 1

# Equation 21 (s)

G0[21,15] = 1
G1[21,4]  = 1

# Equation 22 (pi_D)

G0[22,12] = 1
G1[22,5]  = 1

# Equation 23 (pi_I)

G0[23,13] = 1
G1[23,6]  = 1

# Equation 24 (i)

G0[24,16] = 1
G1[24,8]  = 1

# G0[findall(G0 .== nothing)] .= 0
# G1[findall(G1 .== nothing)] .= 0
# Psi[findall(Psi .== nothing)] .= 0
# Pi[findall(Pi .== nothing)] .= 0

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
# Based on this guide: https://ase.tufts.edu/chemistry/lin/images/FortranMD_TeachersGuide.pdf
# Uses DFTK in place of LJ for the production stage as a proof-of-concept ab initio MD simulation
# Note that the choice of parameters is for demonstration purposes, and the results are non-physical

using Atomistic
using AtomsBase
using DFTK
using NBodySimulator
using Unitful
using UnitfulAtomic

setup_threading(n_blas = 4)

N = 8
element = :Ar
box_size = 1.5u"nm" # this number was chosen arbitrarily
reference_temp = 94.4u"K"
thermostat_prob = 0.1 # this number was chosen arbitrarily
Δt = 1e-2u"ps"

pspkey = list_psp(:Ar, functional = "lda")[1].identifier
initial_bodies = generate_bodies_in_cell_nodes(N, element, box_size, reference_temp)
for body ∈ initial_bodies
    body.data[:pseudopotential] = pspkey
end
initial_system = FlexibleSystem(initial_bodies, CubicPeriodicBoundaryConditions(austrip(box_size)))

eq_steps = 20000
eq_thermostat = AndersenThermostat(austrip(reference_temp), thermostat_prob / austrip(Δt))
eq_simulator = NBSimulator(Δt, eq_steps, thermostat = eq_thermostat)
potential = LennardJonesParameters(1.657e-21u"J", 0.34u"nm", 0.765u"nm")

eq_result = @time simulate(initial_system, eq_simulator, potential)

display(@time plot_temperature(eq_result, eq_simulator.steps ÷ 200))
display(@time plot_energy(eq_result, eq_simulator.steps ÷ 200))
display(@time plot_rdf(eq_result, potential.σ, 0.5))

ab_initio_steps = 200
ab_initio_simulator = NBSimulator(Δt, ab_initio_steps, t₀ = get_time_range(eq_result)[end])
dftk_potential = DFTKPotential(5u"hartree", [1, 1, 1]; damping = 0.7) # very non-physical but fast for demonstration purposes

ab_initio_result = @time simulate(get_system(eq_result), ab_initio_simulator, dftk_potential)

# Plotting on separate plots because the timespan is so much smaller than in the first phase
display(plot_temperature(ab_initio_result, 1))
display(plot_energy(ab_initio_result, 1))
display(@time plot_rdf(ab_initio_result, potential.σ))

write_nbs_animation(ab_initio_result, "artifacts/argon_ab_initio.gif")
write_ase_trajectory(ab_initio_result, "artifacts/argon_ab_initio.traj")

;


text:
function constraint_jacobian_residual!(prob::GameProblem, pdtraj::PrimalDualTraj)
    N = prob.probsize.N
    n = prob.probsize.n
    p = prob.probsize.p
    pu = prob.probsize.pu
    game_con = prob.game_con
    jac_sub = prob.core.jac_sub

    stamp = Stamp()
    # State constraints
    for i = 1:p
        for conval in game_con.state_conval[i]
            TrajectoryOptimization.evaluate!(conval, pdtraj.pr)
            TrajectoryOptimization.jacobian!(conval, pdtraj.pr)
            TrajectoryOptimization.cost_expansion!(conval)
            for (j,k) in enumerate(conval.inds)
                stampify!(stamp,  :opt, i, :x, 1, k, :x, 1, k)
                valid(stamp,N,p)  ? add2sub(jac_sub[stamp],  conval.hess[j]) : nothing
            end
        end
    end
    # Control constraints
    for conval in game_con.control_conval
        TrajectoryOptimization.evaluate!(conval, pdtraj.pr)
        TrajectoryOptimization.jacobian!(conval, pdtraj.pr)
        TrajectoryOptimization.cost_expansion!(conval)
        for (j,k) in enumerate(conval.inds)
            for i = 1:p
                stampify!(stamp,  :opt, i, :u, i, k, :u, i, k)
                valid(stamp,N,p)  ? add2sub(jac_sub[stamp],  conval.hess[j][pu[i],pu[i]]) : nothing # Be careful this assumes that the control constraints are decoupled per player
                # To introduce control constraints coupled per player we would need to have a separate set of control constraint for each player.
            end
        end
    end
    return nothing
end

function constraint_residual!(prob::GameProblem, pdtraj::PrimalDualTraj)
    N = prob.probsize.N
    n = prob.probsize.n
    p = prob.probsize.p
    pu = prob.probsize.pu
    game_con = prob.game_con
    res_sub = prob.core.res_sub

    vstamp = VStamp()
    # State constraints
    for i = 1:p
        for conval in game_con.state_conval[i]
            TrajectoryOptimization.evaluate!(conval, pdtraj.pr)
            TrajectoryOptimization.jacobian!(conval, pdtraj.pr)
            TrajectoryOptimization.cost_expansion!(conval)
            for (j,k) in enumerate(conval.inds)
                stampify!(vstamp, :opt, i, :x, 1, k)
                valid(vstamp,N,p) ? add2sub(res_sub[vstamp], conval.grad[j]) : nothing
            end
        end
    end
    # Control constraints
    for conval in game_con.control_conval
        TrajectoryOptimization.evaluate!(conval, pdtraj.pr)
        TrajectoryOptimization.jacobian!(conval, pdtraj.pr)
        TrajectoryOptimization.cost_expansion!(conval)
        for (j,k) in enumerate(conval.inds)
            for i = 1:p
                stampify!(vstamp, :opt, i, :u, i, k)
                valid(vstamp,N,p) ? add2sub(res_sub[vstamp], conval.grad[j][pu[i]]) : nothing # Be careful this assumes that the control constraints are decoupled per player
                # To introduce control constraints coupled per player we would need to have a separate set of control constraint for each player.
            end
        end
    end
    return nothing
end