text:
# # ClimateMARGO.jl tutorial

# The following schematic shows the full formulation of the MARGO model, which is described in detail in the accompanying (but not yet peer-reviewed) manuscript, available for free at [EarthArXiv.org/5bgyc](https://eartharxiv.org/5bgyc). No-policy baseline emissions $q(t)$ are prescribed, leading to changing greenhouse gas (CO$_{2e}$) concentrations, radiative forcing, temperatures, and climate damages. Emissions can be decreased by **M**itigation; concentrations can be decreased by carbon dioxide **R**emoval; forcing can be decreased by solar **G**eo-engineering; and the "adapted" temperature that leads to climate damages can be reduced by **A**dapting to the changed climate.
# <img src="../docs/src/MARGO_schematic.png" alt="drawing" width="60%" style="float:left"/>

# ### Import software

# The MARGO model has been implemented and registered as the Julia package `ClimateMARGO.jl`, which can be installed with the command `Pkg.add("ClimateMARGO")`.

# If you are running this tutorial via Binder, there is no need to install the package; just import it using the command below.

using ClimateMARGO # Julia implementation of the MARGO model
using PyPlot # A basic plotting package

#

using ClimateMARGO.Models
using ClimateMARGO.Utils
using ClimateMARGO.Diagnostics

# ## Model configuration

# An instance of MARGO `ClimateModelParameters` requires specifying a `name` and the parameters for three model subcomponents:
# 1. `Domain`: determines the time period of interest, as well as the timestep and– in the policy response mode– the "present" decision-making year.
# 2. `Economics`: determines exogenous economic growth, the cliamte damage cost function, and the cost function for each of the four controls, and the discount rate.
# 3. `Physics`: determines the parameters that govern the aiborn fraction of emissions that remain in the atmosphere, the strength of climate feedbacks, and the rate of ocean heat uptake.

# An instance of the MARGO `ClimateModel` is determined by combining the prescribed set of `ClimateModelParameters` with timeseries of the four `Controls`, which will act as the control variables in the optimization framework:
# 1. **M**itigation of baseline emissions
# 2. **R**emoval of carbon dioxide from the atmopshere
# 3. **G**eongineering by solar radiation management, and
# 4. **A**daptating to the changed climate to reduce the damages of climate change impacts.

# Here, we run through the full construction of the model with its "default" parameter values.

name = "default";

# #### 1. Setting up the temporal grid
# First, we need to set up a time-frame for our experiment. Let's begin our simulation in the year 2020– present day– and consider out to 2200, with a 5-year timestep for computational efficiency.

initial_year = 2020. # [yr]
final_year = 2200. # [yr]
dt = 5. # [yr]
t_arr = t(initial_year, final_year, dt);

# While the model allows for shifting the "present" year forward or backward in time to simulate past and future policy decision-making processes, we will only consider the simplest case in which we take the perspective of a policy decision maker in the year 2020, which is also the initial year of our simulation.

present_year = initial_year
dom = Domain(dt, initial_year, initial_year, final_year);

#### 2. Configuring the carbon cycle and energy balance models

# CO$_{2e}$ concentrations are given by:
# \begin{equation}
# c_{M,R}(t) = c_{0} + \int_{t_{0}}^{t} rq(t')(1-M(t')) \text{ d}t' - q_{0} \int_{t_{0}}^{t} R(t')\text{ d}t'.
# \end{equation}
# where $q(t)$ is a specific timeseries of no-climate-policy baseline emissions (with $q(t_{0}) = q_{0}$), $c_{0}$ is the initial CO$_{2e}$ concentration, and $r$ is the fraction of emitted CO$_{2e}$ that remains in the atmosphere net of short-term uptake by the ocean and biosphere. $M(t')$ and $R(t')$ represent the effects emissions mitigation and carbon dioxide removal, respectively, on CO$_{2e}$ concentrations relative to the baseline; they will be initially set to zero and later determined by solving an optimization problem.

# The default no-policy scenario is one of rapid, fossil-fueled growth which leads to an accumulation of $c(t) = 1400$ ppm of CO$_{2e}$ in the atmosphere by 2150, when emissions are assumed to finally reach net-zero.

c0 = 460. # [ppm]
r = 0.5; # [1] fraction of emissions remaining after biosphere and ocean uptake (Solomon 2009)

q0 = 7.5
q0mult = 3.
t_peak = 2100.
t_zero = 2150.
q = ramp_emissions(t_arr, q0, q0mult, t_peak, t_zero);

#

figure(figsize=(9,3.5))

subplot(1,2,1)
plot(t_arr, ppm_to_GtCO2(q))
xlabel("year")
ylabel(L"baseline emissions $q(t)$ [GtCO2 / year]")
xlim([2020, 2200])
grid(true)

subplot(1,2,2)
q_effective = effective_emissions(r, q, 0., 0.) # No mitigation, no carbon removal
c_baseline = c(c0, q_effective, dt)
plot(t_arr, c_baseline)
xlabel("year")
ylabel(L"baseline concentrations $c(t)$ [ppm]")
xlim([2020, 2200])
tight_layout();
grid(true)
gcf()

# These CO$_{2e}$ concentrations drive an anomalous greenhouse effect, which is represented by the radiative forcing
# \begin{equation}
# F_{M,R,G} = a \ln\left(\frac{c_{M,R}(t)}{c_{0}}\right) - G(t)F_{\infty},
# \end{equation}
# where $a$ is an empirically determined coefficient, $G(t)$ represents the effects of geoengineering, and $F_{\infty} = 8.5$ W/m$^2$ is a scaling factor for the effects of geoengineering by Solar Radiation Modification (SRM).

a = (6.9/2.)/log(2.); # F4xCO2/2 / log(2) [W m^-2]
Finf = 8.5;

#

figure(figsize=(4.5,3.2))
F0 = 3.0
F_baseline = F(a, c0, Finf, c_baseline, 0.)
plot(t_arr, F_baseline .+ F0)
xlabel("year")
ylabel(L"baseline radiative forcing $F(t)$ [W/m$^2$]")
xlim([2020, 2200])
grid(true)
ylim([0,10.]);
gcf()

# Next, we configure MARGO's energy balance model, which is forced by the controlled forcing $F_{M,R,G}$. The two-layer energy balance model can be solved, approximately, as:
# \begin{equation}
#     T_{M,R,G}(t) - T_{0} = \frac{F_{M,R,G}(t)}{B + \kappa} + \frac{\kappa}{B} \int_{t_{0}}^{t} \frac{e^{\frac{t'-t}{\tau_{D}}}}{\tau_{D}} \frac{F_{M,R,G}(t')}{B+\kappa} \, \text{d}t',
# \end{equation}
# where $T_{0}$ is the initial temperature change relative to preindustrial, $B$ is the feedback parameter, $\kappa$ is , and the timescale $\tau_{D} = \dfrac{C_{D}}{B} \dfrac{B + \kappa}{\kappa}$ is a more convenient expression for the deep ocean heat capacity $C_{D}$.

# Default parameter values are taken from [Geoffroy et al. 2013](https://journals.ametsoc.org/doi/full/10.1175/JCLI-D-12-00195.1)'s physically-based calibration of the two-layer model to CMIP5.

# Two-layer EBM parameters
B = 1.13; # Feedback parameter [J yr^-1 m^-2 K^-1]
Cd = 106.; # Deep ocean heat capacity [J m^-2 K^-1]
κ = 0.73; # Heat exchange coefficient [J yr^-1 m^2 K^-1]

# Initial condition: present-day temperature, relative to pre-industrial
T0 = 1.1; # [degC] Berkeley Earth Surface Temperature (Rohde 2013)

print("τD = ", Int64(round(τd(Cd, B, κ))), " years")

# These physical parameters can be used to diagnose the climate sensitivity to a doubling of CO$_{2}$ ($ECS$).

print("ECS = ", round(ECS(a, B),digits=1) ,"ºC")

# Combined, these parameters define the physical model, which is instantiated by the calling the `Physics` constructor method:

Phys = Physics(c0, T0, a, B, Cd, κ, r);

# #### 3. Configuring the simple economic model
# Economic growth in MARGO (in terms of Gross World Product, GWP) is exogenous $E(t) = E_{0} (1 + \gamma)^{(t-t_{0})}$ and is entirely determined by the growth rate $\gamma$. By default, we set $\gamma = 2\%$.

E0 = 100. # Gross World Product at t0 [10^12$ yr^-1]
γ = 0.02 # economic growth rate

E_arr = E(t_arr, E0, γ)

figure(figsize=(4, 2.5))
plot(t_arr, E_arr)
xlabel("year")
ylabel("GWP [trillion USD]")
xlim([2020, 2200])
ylim([0, 3000]);
grid(true)
gcf()

# Economic damages $D_{M,R,G,A} = \tilde{\beta}E(t) (\Delta T_{M,R,G})^{2} (1 - A(t))$, expressed as a fraction $\tilde{\beta}(\Delta T_{M,R,G})^{2}$ of the instantaneous Gross World Product, increase quadratically with temperature change relative to preindustrial. They can be decrease by adaptation controls $A(t)$. The default value of the damage parameter $\tilde{\beta}$ corresponds to damages of 2\% of GWP at 3ºC of warming.

βtilde = 0.01; # damages [%GWP / celsius^2]

# The total cost of climate controls is simply the sum of their individual costs, each of which follows a parabolic cost function:
# \begin{equation}
#     \mathcal{C}_{M,R,G,A} = \mathcal{C}_{M}M^2 + \mathcal{C}_{R}R^2 +  \mathcal{C}_{G}G^2 + \mathcal{C}_{A}A^2
# \end{equation}

# The calculation of the reference control costs $\mathcal{C}_{M}, \mathcal{C}_{R}, \mathcal{C}_{G}, \mathcal{C}_{A}$ are somewhat more complicated; see our Methods in [the preprint](https://eartharxiv.org/5bgyc/) and `defaults.jl` for details. Here, we simply provide their default numerical values, where the costs of mitigation $\mathcal{C}_{M} = \tilde{\mathcal{C}}_{M} E(t)$ and geoengineering $\mathcal{C}_{G} = \tilde{\mathcal{C}}_{G} E(t)$ grow with the size of the global economy and are thus specified as a fraction of GWP, while adaptaiton and removal costs are in trillions of USD per year.

ti = findall(x->x==2100, t_arr)[1]
mitigate_cost_percentGWP = 0.02
mitigate_cost = mitigate_cost_percentGWP*E_arr[ti]/ppm_to_GtCO2(q[ti]); # [trillion USD / year / GtCO2]

# Costs of negative emissions technologies [US$/tCO2]
costs = Dict(
    "BECCS" => 150.,
    "DACCS" => 200.,
    "Forests" => 27.5,
    "Weathering" => 125.,
    "Biochar" => 70.,
    "Soils" => 50.
)

# Upper-bound potential for sequestration (GtCO2/year)
potentials = Dict(
    "BECCS" => 5.,
    "DACCS" => 5.,
    "Forests" => 3.6,
    "Weathering" => 4.,
    "Biochar" => 2.,
    "Soils" => 5.
)

mean_cost = sum(values(potentials) .* values(costs)) / sum(values(potentials)) # [$ tCO2^-1] weighted average
CDR_potential = sum(values(potentials))
CDR_potential_fraction = CDR_potential / ppm_to_GtCO2(q0)

# Estimate cost from Fuss 2018 (see synthesis Figure 14 and Table 2)
remove_cost = (mean_cost * CDR_potential*1.e9) / (CDR_potential_fraction^3) * 1.e-12; # [trillion USD / year]

adapt_cost = 0.115; # [%GWP / year] directly from de Bruin, Dellink, and Tol (2009)

geoeng_cost = βtilde*(Finf/B)^2; # [% GWP]

# Climate damages and control costs are discounted at the relatively low rate of $\rho = 2\%$, such that future damages and costs are reduced by a multiplicative discount factor $(1 - \rho)^{(t-t_{0})}$.

ρ = 0.02;

figure(figsize=(4, 2.5))
plot(t_arr, discount(t_arr, ρ, present_year)*100)
xlabel("year")
ylabel("discount factor (%)")
xlim([2020, 2200])
ylim([0, 100]);
grid(true)
gcf()

# These parameters, in addition to a no-policy baseline emissions time-series and present-day control values, define the economic model.

Econ = Economics(
    E0, γ, βtilde, ρ, Finf,
    mitigate_cost, remove_cost, geoeng_cost, adapt_cost,
    0., 0., 0., 0., # Initial condition on control deployments at t[1]
    q
);

# #### A MARGO model configuration is uniquely determined by the parameters defined above

params = ClimateModelParameters(
    name,
    dom,
    Econ,
    Phys
);

# ## Instanciating the MARGO climate model

# Along with economic and physical model components, the timeseries for each of the four controls must be specified. By default, we simply set these to zero.

Cont = Controls(
    zeros(size(t_arr)), # mitigate
    zeros(size(t_arr)), # remove
    zeros(size(t_arr)), # geoeng
    zeros(size(t_arr))  # adapt
);

# The above parameters determine the full configuration of the MARGO model, which is instanced using the `ClimateModel` constructor method:

m = ClimateModel(
    params,
    Cont
);

# ## Model optimization

# #### Formulating the optimization problem
# By default, the optimization problem we solve is for the most cost-effective combination of controls, as determined by minimization of the discounted net present value,

# \begin{equation}
#     \min\left\{\int_{t_{0}}^{t_{f}} \mathcal{C}_{M,R,G,A} (1 + \rho)^{-(t-t_{0})} \text{ d}t\right\}
# \end{equation}

# which keep controlled damages below the level corresponding to a chosen temperature threshold $T^{\star}$,

# \begin{equation}
#     \tilde{\beta} E(t) (T_{M,R,G})^{2} (1 - A(t)) < \tilde{\beta} E(t) (T^{\star})^{2}.
# \end{equation}

temp_goal = 2.0;

# #### Additional policy / inertia constraints
# The optimization is also subject to the constraints described below.

# First, we set a upper bounds on the maximum plausible deployment of each control.

max_deployment = Dict("mitigate"=>1., "remove"=>1., "geoeng"=>1., "adapt"=>0.4);

# Second, we set upper limits on how quickly each control can be ramped up or down.

# (Adaptation is treated differently since we it interpret it as buying insurance against future climate damages, although the financing is spread evenly over the entire period.)

max_slope = Dict("mitigate"=>1. /40., "remove"=>1. /40., "geoeng"=>1. /40., "adapt"=> 1. /40.);

# Third, we impose restrictions on when controls can be first deployed. In particular, since carbon dioxide removal and solar radiation modification do not yet exist at scale, we delay these until 2030 and 2050, respectively, at the earliest.

delay_deployment = Dict(
    "mitigate"=>0.,
    "remove"=>0.,
    "geoeng"=>0.,
    "adapt"=>0.
);

# #### Running the optimization
# The optimization takes about ~40 seconds the first time it is run as the code compiles, but runs virtually instantly afterwards, even if model parameter values are changed.

using ClimateMARGO.Optimization

@time msolve = optimize_controls!(
    m,
    obj_option = "adaptive_temp", temp_goal = temp_goal,
    max_deployment = max_deployment, max_slope = max_slope, delay_deployment = delay_deployment
);

# ## Plotting MARGO results

# We provide some convenient functions for plotting basic model results.

using ClimateMARGO.Plotting
plot_state(m, temp_goal=temp_goal)
gcf()

# ## Customizing MARGO

# Try changing some of the parameters above and re-running the model!

# For example, in the simulation below we set a more stringent temperature goal of $T^{\star} = 1.5$ and omit solar geoengineering and adaptation completely (by setting their maximum deployment to zero).

temp_goal = 1.5
max_deployment["geoeng"] = 0.
max_deployment["adapt"] = 0.

@time optimize_controls!(
    m,
    obj_option = "adaptive_temp", temp_goal = temp_goal,
    max_deployment = max_deployment, max_slope = max_slope, delay_deployment = delay_deployment
);

plot_state(m, temp_goal=temp_goal)
gcf()

# ## Saving (and loading) MARGO simulations (or parameter configurations)

using ClimateMARGO.IO
export_path = tempname() * ".json"

#

export_parameters(export_path, params);

#

new_params = import_parameters(export_path)

# Let's say we want to see what happens if climate feedbacks $B$ are twice as strong as they are now (which means a lower feedback parameter).

new_params.name = "stronger_feedbacks"
new_params.physics.B /= 2.;

new_m = ClimateModel(new_params, Cont)

# The new equilibrium climate sensitivity (recall $ECS = F_{2x}/B$) is now:

print("ECS ≈ $(round(ECS(new_m), digits=1))°C")


text:
using Gen
using GenDirectionalStats
import Rotations
import Random
using Plots

###########################################
# plots of distributions on 3D directions #
###########################################

function show_uniform_3d_direction()
    p = plot(aspect_ratio=:equal, title="uniform_3d_dir()")
    for i in 1:1000
        x = uniform_3d_direction()
        scatter!(p, [x.v[1]], [x.v[2]], [x.v[3]], w=1, label=nothing, color="red",
            markersize=5, markerstrokealpha=0.0, alpha=0.1)
    end
    xlims!(-1.0, 1.0)
    ylims!(-1.0, 1.0)
    zlims!(-1.0, 1.0)
    return p
end

function show_vmf_3d_direction(mu, k)
    p = plot([0.0, mu.v[1]], [0.0, mu.v[2]], [0.0, mu.v[3]], w=3, label=nothing, color="blue", aspect_ratio=:equal, title="vmf_3d_dir($(mu.v), $k)")
    for i in 1:1000
        x = vmf_3d_direction(mu, k)
        scatter!(p, [x.v[1]], [x.v[2]], [x.v[3]], w=1, label=nothing, color="red",
            markersize=5, markerstrokealpha=0.0, alpha=0.1)
    end
    xlims!(-1.0, 1.0)
    ylims!(-1.0, 1.0)
    zlims!(-1.0, 1.0)
    return p
end

function show_3d_direction_dists()
    println("generating plots for 3D directions...")
    mu = UnitVector3(1.0, 0.0, 0.0)
    p1 = show_vmf_3d_direction(mu, 0.1)
    p2 = show_vmf_3d_direction(mu, 1.0)
    p3 = show_vmf_3d_direction(mu, 10.0)
    p4 = show_vmf_3d_direction(mu, 100.0)
    mu = UnitVector3(0.0, -1.0, 0.0)
    p5 = show_vmf_3d_direction(mu, 0.1)
    p6 = show_vmf_3d_direction(mu, 1.0)
    p7 = show_vmf_3d_direction(mu, 10.0)
    p8 = show_vmf_3d_direction(mu, 100.0)
    plot(
        show_uniform_3d_direction(), p1, p2, p3, p4,
        show_uniform_3d_direction(), p5, p6, p7, p8, layout=(2, 5), size=(2000, 800))
    savefig("vmf_3d_direction.png")
end

show_3d_direction_dists()

##########################################
# plots of distributions on 3D rotations #
##########################################

function scatter_rot3!(p, rot::Rot3)
    x_axis = rot[:,1]
    y_axis = rot[:,2]
    z_axis = rot[:,3]
    scatter!(p, [x_axis[1]], [x_axis[2]], [x_axis[3]], w=1, label=nothing, color="red",
        markersize=5, markerstrokealpha=0.0, alpha=0.1)
    scatter!(p, [y_axis[1]], [y_axis[2]], [y_axis[3]], w=1, label=nothing, color="green",
        markersize=5, markerstrokealpha=0.0, alpha=0.1)
    scatter!(p, [z_axis[1]], [z_axis[2]], [z_axis[3]], w=1, label=nothing, color="blue",
        markersize=5, markerstrokealpha=0.0, alpha=0.1)
end

function plot_rot3!(p, rot::Rot3; alpha=1.0, w=5)
    x_axis = rot[:,1]
    y_axis = rot[:,2]
    z_axis = rot[:,3]
    plot!(p, [0.0, x_axis[1]], [0.0, x_axis[2]], [0.0, x_axis[3]], label=nothing, color="red", w=w, alpha=alpha)
    plot!(p, [0.0, y_axis[1]], [0.0, y_axis[2]], [0.0, y_axis[3]], label=nothing, color="green", w=w, alpha=alpha)
    plot!(p, [0.0, z_axis[1]], [0.0, z_axis[2]], [0.0, z_axis[3]], label=nothing, color="blue", w=w, alpha=alpha)
end

function show_uniform_3d_rotation()
    p = plot(aspect_ratio=:equal, title="uniform_rot3()")
    for i in 1:250
        rot = uniform_rot3()
        scatter_rot3!(p, rot)
        if i == 1
            plot_rot3!(p, rot; alpha=0.5, w=2)
        end
    end
    xlims!(-1.0, 1.0)
    ylims!(-1.0, 1.0)
    zlims!(-1.0, 1.0)
    return p
end

function show_vmf_3d_rotation(mu, k)
    p = plot(aspect_ratio=:equal, title="vmf_rot3(I, $k)")
    plot_rot3!(p, mu)
    for i in 1:250
        rot = vmf_rot3(mu, k)
        scatter_rot3!(p, rot)
        if i == 1
            plot_rot3!(p, rot; alpha=0.5, w=2)
        end
    end
    xlims!(-1.0, 1.0)
    ylims!(-1.0, 1.0)
    zlims!(-1.0, 1.0)
    return p
end

function show_uniform_vmf_3d_rotation(mu, k, prob_outlier)
    p = plot(aspect_ratio=:equal, title="uniform_vmf_rot3(I, $k, $prob_outlier)")
    plot_rot3!(p, mu)
    for i in 1:250
        rot = uniform_vmf_rot3(mu, k, prob_outlier)
        scatter_rot3!(p, rot)
        if i == 1
            plot_rot3!(p, rot; alpha=0.5, w=2)
        end
    end
    xlims!(-1.0, 1.0)
    ylims!(-1.0, 1.0)
    zlims!(-1.0, 1.0)
    return p
end

function show_3d_rotation_dists()
    println("generating plots for 3D rotations...")
    mu = Rot3(1.0, 0.0, 0.0, 0.0)
    p1 = show_vmf_3d_rotation(mu, 0.1)
    p2 = show_vmf_3d_rotation(mu, 1.0)
    p3 = show_vmf_3d_rotation(mu, 10.0)
    p4 = show_vmf_3d_rotation(mu, 100.0)
    p5 = show_vmf_3d_rotation(mu, 1000.0)
    p6 = show_uniform_vmf_3d_rotation(mu, 1000.0, 0.5)
    plot(
        p1, p2, p3, p4, p5,
        layout=(1, 5), size=(2000, 400))
    savefig("vmf_3d_rotation_1.png")
    plot(
        show_uniform_3d_rotation(), p6,
        layout=(1, 2), size=(1000, 400))
    savefig("vmf_3d_rotation_2.png")
end

show_3d_rotation_dists()

#############################################
# plots of distributions on plane rotations #
#############################################

function scatter_rot2!(p, rot::Rot2)
    rotmat = Rotations.RotMatrix{2}(rot)
    x_axis = rotmat[:,1]
    y_axis = rotmat[:,2]
    scatter!(p, [x_axis[1]], [x_axis[2]], w=1, label=nothing, color="red",
        markersize=5, markerstrokealpha=0.0, alpha=0.1)
    scatter!(p, [y_axis[1]], [y_axis[2]], w=1, label=nothing, color="green",
        markersize=5, markerstrokealpha=0.0, alpha=0.1)
end

function plot_rot2!(p, rot::Rot2; alpha=1.0, w=5)
    x_axis = rot[:,1]
    y_axis = rot[:,2]
    plot!(p, [0.0, x_axis[1]], [0.0, x_axis[2]], label=nothing, color="red", w=w, alpha=alpha)
    plot!(p, [0.0, y_axis[1]], [0.0, y_axis[2]], label=nothing, color="green", w=w, alpha=alpha)
end

function show_uniform_2d_rotation()
    p = plot(aspect_ratio=:equal, title="uniform_rot2()")
    for i in 1:250
        rot = uniform_rot2()
        scatter_rot2!(p, rot)
        if i == 1
            plot_rot2!(p, rot; alpha=0.5, w=2)
        end
    end
    xlims!(-1.1, 1.1)
    ylims!(-1.1, 1.1)
    return p
end

function show_vmf_2d_rotation(mu, k)
    p = plot(aspect_ratio=:equal, title="von_mises_rot2(I, $k)")
    plot_rot2!(p, mu)
    for i in 1:250
        rot = von_mises_rot2(mu, k)
        scatter_rot2!(p, rot)
        if i == 1
            plot_rot2!(p, rot; alpha=0.5, w=2)
        end
    end
    xlims!(-1.1, 1.1)
    ylims!(-1.1, 1.1)
    return p
end

function show_2d_rotation_dists()
    println("generating plots for 2D rotations...")
    mu = Rotations.RotMatrix{2}(0.0)
    @time p1 = show_vmf_2d_rotation(mu, 0.1)
    @time p2 = show_vmf_2d_rotation(mu, 1.0)
    @time p3 = show_vmf_2d_rotation(mu, 10.0)
    @time p4 = show_vmf_2d_rotation(mu, 100.0)
    @time p5 = show_vmf_2d_rotation(mu, 1000.0)
    @time plot(
        show_uniform_2d_rotation(), p1, p2, p3, p4, p5,
        layout=(2, 3), size=(900, 600))
    @time savefig("vmf_2d_rotation.png")
end

show_2d_rotation_dists()

ext: jl


text:
using BoardGames
using Cairo

mutable struct  OthelloBoard
    v::Matrix{Int}
    turn::Int
end

struct Othello <:Game{OthelloBoard, Tuple{Int,Int}} end

function OthelloBoard()
    v = zeros(Int, (8,8))
    v[4,4] = -1
    v[4,5] =  1
    v[5,4] =  1
    v[5,5] = -1

    return OthelloBoard(v, 1)
end

function BoardGames.getmoves(board::OthelloBoard)
    ans = Tuple{Int,Int}[]

    for i in 1:8, j in 1:8
        if board.v[i,j] != 0
            continue
        end
        for direction in ((i,j) for i in -1:1 for j in -1:1
                                    if i != 0 || j != 0)
            pos = (i,j) .+ direction
            if checkmove(board.v, pos, direction, board.turn, false)
                push!(ans, (i,j))
                break
            end
        end
    end
    return ans
end

function hasmoves(board::OthelloBoard)
    for i in 1:8, j in 1:8
        if board.v[i,j] != 0
            continue
        end
        for direction in ((i,j) for i in -1:1 for j in -1:1
                                    if i != 0 || j != 0)
            pos = (i,j) .+ direction
            if checkmove(board.v, pos, direction, board.turn, false)
                return true
            end
        end
    end
    return false
end

function BoardGames.play(board::OthelloBoard, move::Tuple{Int,Int})
    validmove = false
    v = copy(board.v)
    if board.v[move...] != 0
        error("illegal play $play")
    end

    v[move...] = board.turn

    for direction in ((i,j) for i in -1:1 for j in -1:1
                                if i != 0 || j != 0)
        if checkmove(v, move .+ direction, direction, board.turn, true)
            validmove = true
        end
    end

    if !validmove
        error("illegal play $play")
    end

    if !hasmoves(OthelloBoard(v, -board.turn))
        return OthelloBoard(v, board.turn)
    else
        return OthelloBoard(v, -board.turn)
    end
end

function BoardGames.play!(board::OthelloBoard, move::Tuple{Int,Int})
    validmove = false

    if board.v[move...] != 0
        error("illegal play $play")
    end

    board.v[move...] = board.turn

    for direction in ((i,j) for i in -1:1 for j in -1:1
                                if i != 0 || j != 0)
        if checkmove(board.v, move .+ direction, direction, board.turn, true)
            validmove = true
        end
    end

    if !validmove
        error("illegal play $play")
    end

    board.turn = -board.turn
    if hasmoves(board)
        return board
    else
        board.turn = -board.turn
        return board
    end
end

function checkmove(b, x, direction, player, fill::Bool, count=1)
    if !(1 <= x[1] <= 8) || !(1 <= x[2] <= 8) || b[x...] == 0
        return false
    end
    if b[x...] == player
        if count == 1
            return false
        else
            return true
        end
    end
    ans = checkmove(b, x.+direction, direction, player, fill, count+1)
    if ans && fill
        b[x...] = player
    end
    return ans
end

function BoardGames.playerturn(board::OthelloBoard)
    if board.turn == 1
        return 1
    else
        return 2
    end
end

function BoardGames.winner(board::OthelloBoard)
    if hasmoves(board)
        error("call winner on ongoing game")
    end
    p1 = sum(board.v .== 1)
    p2 = sum(board.v .== -1)
    if p1 > p2
        return 1
    elseif p2 > p1
        return 2
    else
        return 0
    end
end

function Base.show(io::IO, b::OthelloBoard)
    print(io, "OtheloBoard($(b.v), $(b.turn))")
end

function Base.show(io::IO, ::MIME"text/plain", b::OthelloBoard)
    icon = Dict(-1=>"o", 0=>".", 1=>"x")

    s = "\n"

    for i in 1:8
        s *= join(map(x->icon[x], b.v[i,:])) * "\n"
    end
    print(io,s)
end

function Base.show(io::IO, m::MIME{Symbol("image/png")}, b::OthelloBoard)
    c, cr = baseboard(800., 800.)
    for i in 1:8, j in 1:8
        if b.v[i,j] != 0
            draw_circle(c, cr, (i,j), b.v[i,j])
        end
    end
    show(io, m, c)
    nothing
end

function baseboard(w::Real=800., h::Real=800.)
    c = CairoRGBSurface(w,h)
    cr = CairoContext(c)

    save(cr) #TODO: what does this do?
    set_source_rgb(cr, 0., 153/255, 51/255)
    rectangle(cr,0.0,0.0,w,h)
    fill(cr)
    restore(cr) #TODO: what does this do?

    ## original example, following here

    for x in 1:7
        set_source_rgba(cr, 0, 0, 0, 1)
        set_line_width(cr, 3.0)
        move_to(cr, 0., x/8*w); line_to(cr, w, x/8*w);
        move_to(cr, x/8*w, 0.); line_to(cr, x/8*w, h);
        stroke(cr);
    end
    (c, cr)
end

function draw_circle(c, cr, i::Tuple{Int,Int}, player::Int, w::Real=800., h::Real=800.)
    if player == 1
        set_source_rgb(cr,0.0,0.0,0.0)
    else
        set_source_rgb(cr,1.0,1.0,1.0)
    end

    x = w/16 + (i[1]-1)*w/8
    y = h/16 + (i[2]-1)*h/8
    circle(cr, x, y, min(w,h)/24)
    fill_preserve(cr)
    stroke(cr)
end


text:
@testset "knights" begin
    b = Bobby.set_board()

    nv = Bobby.gen_night_valid(Bobby.INT2UINT[58])
    @test length(nv) == 3
    nvw = [0,0,0,0,0,0,0,0,
           0,0,0,0,0,0,0,0,
           0,0,0,0,0,0,0,0,
           0,0,0,0,0,0,0,0,
           0,0,0,0,0,0,0,0,
           1,0,1,0,0,0,0,0,
           0,0,0,1,0,0,0,0,
           0,0,0,0,0,0,0,0]
    @test all(Int.(Bobby.cvt_to_bitarray(nv)) .== nvw)

    anv = Bobby.gen_all_night_valid_moves()
    @test all(Int.(Bobby.cvt_to_bitarray(anv[Bobby.INT2UINT[58]])) .== nvw)

    @test test_fen("k7/8/8/8/8/8/8/N6K w - - 0 1", 1, [5])
    @test test_fen("k7/8/8/8/8/8/1N6/7K w - - 0 1", 1, [7])
    @test test_fen("k7/8/8/8/P1P5/3P4/1N6/7K w - - 0 1", 1, [7])
    @test test_fen("k7/8/8/8/p1p5/3p4/1N6/7K w - - 0 1", 1, [7])
    @test test_fen("k7/8/8/8/3n4/8/2P1P3/7K b - - 0 1", 1, [11])
    @test test_fen("k7/8/8/3n4/1p3p2/1Pp1pP2/2P1P3/7K b - - 0 1", 1, [7])
end


text:
function forward(model::Model, q::judiVector, dobs::judiVector; options=Options(), ps=16, modelPy=nothing)
    # Python model
    isnothing(modelPy) && (modelPy = devito_model(model, options))
    # Interpolate input data to computational grid
    q_data = time_resample_data(q[1], get_dt(model))
    d_data = time_resample_data(dobs[1], get_dt(model))

    # Set up coordinates with devito dimensions
    src_coords = setup_grid(q.geometry[1], modelPy.shape)
    rec_coords = setup_grid(dobs.geometry[1], modelPy.shape)

    # QR probing vector
    ts = trunc(Int, abs(src_coords[1, end] - rec_coords[1, end])/(1.5f0*get_dt(model)))
    Q = qr_data(d_data, ps; ts=ts)

    rec, eu, _ = forward(modelPy, src_coords, rec_coords, q_data, Q, options.space_order, options.isic)
    rec = time_resample_data(rec, dobs.geometry.dt[1], get_dt(model), dobs.geometry.t[1])
    return judiVector(dobs.geometry, rec), Q, eu
end

function born(model::Model, q::judiVector, dobs::judiVector, dm; options=Options(), ps=16, modelPy=nothing)
    # Python model
    isnothing(modelPy) && (modelPy = devito_model(model, options; dm=dm))

    # Interpolate input data to computational grid
    q_data = time_resample_data(q[1], get_dt(model))
    d_data = time_resample_data(dobs[1], get_dt(model))

    # Set up coordinates with devito dimensions
    src_coords = setup_grid(q.geometry[1], modelPy.shape)
    rec_coords = setup_grid(dobs.geometry[1], modelPy.shape)

    # QR probing vector
    ts = trunc(Int, abs(src_coords[1, end] - rec_coords[1, end])/(1.5f0*get_dt(model)))
    Q = qr_data(d_data, ps; ts=ts)

    recnl, recl, eu, _ = born(modelPy, src_coords, rec_coords, q_data, Q, options.space_order, options.isic)
    recnl = time_resample_data(recnl, dobs.geometry.dt[1], get_dt(model), dobs.geometry.t[1])
    recl = time_resample_data(recl, dobs.geometry.dt[1], get_dt(model), dobs.geometry.t[1])
    return judiVector(dobs.geometry, recnl), judiVector(dobs.geometry, recl), Q, eu
end

function born_with_back(model::Model, q::judiVector, dobs::judiVector, dm; options=Options(), ps=16, modelPy=nothing)
    # Python model
    isnothing(modelPy) && (modelPy = devito_model(model, options; dm=dm))
    # Interpolate input data to computational grid
    q_data = time_resample_data(q[1], get_dt(model))
    d_data = time_resample_data(dobs[1], get_dt(model))

    # Set up coordinates with devito dimensions
    src_coords = setup_grid(q.geometry[1], modelPy.shape)
    rec_coords = setup_grid(dobs.geometry[1], modelPy.shape)

    # QR probing vector
    ts = trunc(Int, abs(src_coords[1, end] - rec_coords[1, end])/(1.5f0*get_dt(model)))
    Q = qr_data(d_data, ps; ts=ts)

    recl, eu, ev, _  = born_with_back(modelPy, src_coords, rec_coords, q_data, Q, d_data[end:-1:1, :],
                                      options.space_order, options.isic)
    recl = time_resample_data(recl, dobs.geometry.dt[1], get_dt(model), dobs.geometry.t[1])
    g = combine(eu, ev, options.isic)
    return remove_padding(g.data, modelPy.padsizes; true_adjoint=options.sum_padding), judiVector(dobs.geometry, recl)
end

function backprop(model::Model, residual::judiVector, Q::Array{Float32}, eu::PyObject;
                 options=Options(), ps=16, modelPy=nothing)
    # Python mode
    isnothing(modelPy) && (modelPy = devito_model(model, options))
    # Interpolate input data to computational grid
    d_data = time_resample_data(residual[1], get_dt(model))

    # Set up coordinates with devito dimensions
    rec_coords = setup_grid(residual.geometry[1], modelPy.shape)

    ev, _ = backprop(modelPy, d_data, rec_coords, Q, options.space_order)
    g = combine(eu, ev, options.isic)
    return remove_padding(g.data, modelPy.padsizes; true_adjoint=options.sum_padding)
end

max_stars_repo_name: slimgroup/TimeProbeSeismic

max_issues_repo_name: slimgroup/TimeProbeSeismic

max_forks_repo_name: slimgroup/TimeProbeSeismic


text:
"""
    circular_layout(g) -> x, y

Return two vectors representing the positions of the nodes of graph `g`
when placed on the unit circonference of radius one centered at the origin.
"""
function circular_layout(g::AGraphOrDiGraph)
    # Discard the extra angle since it matches 0 radians.
    θ = range(0, stop=2pi, length=nv(g)+1)[1:end-1]
    return cos.(θ), sin.(θ)
end

"""
    spring_layout(g;
                  k = 1 / √nv(g),
                  maxiter = 100,
                  inittemp = 2.0,
                  [x0, y0]) -> x, y

Return the positions of the nodes in graph `g` according to
Fruchterman and Reingold's spring/repulsion model.
The forces as function of the distance `d` beetween two nodes
are given by

    f_a(d) =  d / k # attractive force
    f_r(d) = -k^2 / d^2 # repulsive force:

`maxiter` is the number of updates of the positions.
`inittemp` controls displacement per iteration.

Initial positions can passed though the argument `x0` and `y0`.

The positions are rescaled to fit the [-1, +1]^2 box.

This function is adapted from [GraphLayout.jl](https://github.com/IainNZ/GraphLayout.jl).
"""
function spring_layout(g::AGraphOrDiGraph;
                       x0 = 2*rand(nv(g)) .- 1.0,
                       y0 = 2*rand(nv(g)) .- 1.0,
                       k = 1/√nv(g),
                       maxiters = 100,
                       inittemp = 2.0)

    x = copy(x0)
    y = copy(y0)
    nvg = nv(g)
    adj_matrix = adjacency_matrix(g)
    k² = k * k

    # Store forces and apply at end of iteration all at once
    force_x = zeros(nvg)
    force_y = zeros(nvg)

    # Iterate MAXITER times
    @inbounds for iter = 1:maxiters
        # Calculate forces
        for i = 1:nvg
            force_vec_x = 0.0
            force_vec_y = 0.0
            for j = 1:nvg
                i == j && continue
                d_x = x[j] - x[i]
                d_y = y[j] - y[i]
                dist²  = (d_x * d_x) + (d_y * d_y)
                dist = sqrt(dist²)

                if !(iszero(adj_matrix[i,j]) && iszero(adj_matrix[j,i]) )
                    # Attractive + repulsive force
                    # F_d = dist² / k - k² / dist # original FR algorithm
                    F_d = dist / k - k² / dist²
                else
                    # Just repulsive
                    # F_d = -k² / dist  # original FR algorithm
                    F_d = -k² / dist²
                end
                force_vec_x += F_d * d_x
                force_vec_y += F_d * d_y
            end
            force_x[i] = force_vec_x
            force_y[i] = force_vec_y
        end
        # Cool down
        temp = inittemp / iter
        # Now apply them, but limit to temperature
        for i = 1:nvg
            fx = force_x[i]
            fy = force_y[i]
            force_mag  = sqrt((fx * fx) + (fy * fy))
            scale      = min(force_mag, temp) / force_mag
            x[i] += force_x[i] * scale
            y[i] += force_y[i] * scale
        end
    end

    # Scale to unit square
    min_x, max_x = minimum(x), maximum(x)
    min_y, max_y = minimum(y), maximum(y)
    function scaler(z, a, b)
        2.0*((z - a) / (b - a)) - 1.0
    end
    map!(z -> scaler(z, min_x, max_x), x, x)
    map!(z -> scaler(z, min_y, max_y), y, y)

    return x, y
end

"""
    shell_layout(g, nlist) -> x, y

Position the nodes of `g` in concentric circles.

`nlist` is a vector of vectors containing the nodes
for each shell.
"""
function shell_layout(g::AGraphOrDiGraph, nlist::Vector{Vector{T}}) where T<:Integer
    if nv(g) == 1
        return [0.0], [0.0]
    end
    radius = length(nlist[1]) > 1 ? 1.0 : 0.0
    x = Float64[]
    y = Float64[]
    for nodes in nlist
        θ = range(0, stop=2pi, length=length(nodes)+1)[1:end-1]
        append!(x, radius * cos.(θ))
        append!(y, radius * sin.(θ))
        radius += 1.0
    end
    return x, y
end


text:
"""
```
simulate_switching(m, scen::SwitchingScenario; verbose = :low)
```

Simulate switching in and out of a default scenario for the
SwitchingScenario scen. See `SwitchingScenario` for more
info. Returns a dictionary of results and the proportion
of times a switch actually occured.
"""
function simulate_switching(m::AbstractModel, scen::SwitchingScenario;
                            verbose::Symbol = :low)

    if VERBOSITY[verbose] >= VERBOSITY[:low]
        info("Simulating switching for " * string(scen.key) * "...")
        println("Start time: " * string(now()))
        println("Outputs will be saved in " * rawpath(m, "scenarios"))
        tic()
    end

    # Revert model alt policy to historical rule
    m <= Setting(:alternative_policy, AltPolicy(:historical, solve, eqcond), false, "apol",
                 "Alternative policy")

    # Get output file names
    original_output_files = get_scenario_output_files(m, scen.original,
                                                      [:forecastobs, :forecastpseudo])
    default_output_files = get_scenario_output_files(m, scen.default,
                                                     [:forecastobs, :forecastpseudo])

    results = Dict{Symbol, Array{Float64}}()
    switching_results = [0.0, 0.0]

    for (i, output_var) in enumerate([:forecastobs, :forecastpseudo])
        # Read in original and default draws
        original_draws = load(original_output_files[output_var], "arr")
        default_draws  = load(default_output_files[output_var], "arr")

        n_draws, n_vars, n_periods = size(original_draws)
        n_default_draws = size(default_draws, 1)

        @assert n_vars == size(default_draws, 2)
        @assert length(scen.probs_enter) == n_periods <= size(default_draws, 3)

        # Simulate switching
        results[output_var] = zeros(n_draws, n_vars, n_periods)
        proportion_switched = 0.0
        for j = 1:n_draws
            k = rand(1:n_default_draws)
            results[output_var][j, :, :], switched = switch(original_draws[j, :, :], default_draws[k, :, :],
                                                            scen.probs_enter, scen.probs_exit)
            proportion_switched += Float64(switched)
        end
        switching_results[i] = proportion_switched / n_draws
    end

    results[:proportion_switched] = [mean(switching_results)]

    # Write outputs
    output_files = get_scenario_output_files(m, scen, [:forecastobs, :forecastpseudo])
    write_scenario_forecasts(m, output_files, results, verbose = verbose)

    # Print
    if VERBOSITY[verbose] >= VERBOSITY[:low]
        switching_time = toq()
        switching_time_min = switching_time/60
        println("\nTime elapsed: " * string(switching_time_min) * " minutes")
        println("Switching complete: " * string(now()))
    end

    return results, switching_results
end

"""
```
switch(original, default, probs_enter, probs_exit)
```

Simulate entry and exit from the original vector into the default one
according to the specified probabilities.
"""
function switch(original::Matrix{Float64}, default::Matrix{Float64},
                probs_enter::Vector{Float64}, probs_exit::Vector{Float64})

    n_vars, n_periods = size(original)
    output = zeros(n_vars, n_periods)

    # Find period of switch from default to original
    period_in = choose_switch_period(probs_enter, 1)

    # Record whether original is actually switched into
    switched = period_in < n_periods

    # Find period of switch from original to default
    # Start in period_in + 1 s.t. at least 1 period is spent in the original scenario
    period_out = choose_switch_period(probs_exit, period_in + 1)

    output[:, 1:(period_in-1)]          = default[:, 1:(period_in-1)]
    output[:, period_in:(period_out-1)] = original[:, 1:(period_out-period_in)]
    output[:, period_out:n_periods]     = default[:, period_out:n_periods]

    return output, switched
end

function choose_switch_period(probs::Vector{Float64}, first_period::Int)
    n_periods = length(probs)

    for t = first_period:n_periods
        if rand() < probs[t]
            return t
        end
    end

    return n_periods + 1
end


text:
using .URIs

const absent = SubString("absent", 1, 0)

function URIs.URI(p::AbstractPath; query=absent, fragment=absent)
    if isempty(p.root)
        throw(ArgumentError("$p is not an absolute path"))
    end

    b = IOBuffer()
    print(b, "file://")

    if !isempty(p.drive)
        print(b, "/")
        print(b, p.drive)
    end

    for s in p.segments
        print(b, "/")
        print(b, URIs.escapeuri(s))
    end

    return URIs.URI(URIs.URI(String(take!(b))); query=query, fragment=fragment)
end


text:
using Functors, Test

@testset "Nested" begin
  struct Foo
    x
    y
  end

  @functor Foo

  struct Bar
    x
  end

  @functor Bar

  model = Bar(Foo(1, [1, 2, 3]))

  model′ = fmap(float, model)

  @test model.x.y == model′.x.y
  @test model′.x.y isa Vector{Float64}
end

@testset "Property list" begin
  struct Baz
    x
    y
    z
  end

  @functor Baz (y,)

  model = Baz(1, 2, 3)
  model′ = fmap(x -> 2x, model)

  @test (model′.x, model′.y, model′.z) == (1, 4, 3)
end


text:
module CapeOpen
  export constantstrings, constantfloats, tempreturedependents, pressuredependents
  export MaterialObject, PropertyPackage, PhysicalPropertyCalculator

  """
    Constant string properties
      1- casRegistryNumber => Chemical Abstract Service Registry Number
      2- chemicalFormula => Chemical formula
      3- iupacName => Complete IUPAC Name
      4- SMILESformula => SMILES chemical structure formula
  """
  constantstrings=ASCIIString[
    "casRegistryNumber", # Chemical Abstract Service Registry Number
    "chemicalFormula", # Chemical formula
    "iupacName", # Complete IUPAC Name
    "SMILESformula" # SMILES chemical structure formula
  ]

  """
    Constant float properties
      1-  acentricFactor => Pitzer acentric factor
      2-  associationParameter => association-parameter (Hayden-O’Connell)
      3-  bornRadius => in m
      4-  charge =>
      5-  criticalCompressibilityFactor => critical compressibility factor Z
      6-  criticalDensity => critical density in mol/m3
      7-  criticalPressure => critical pressure in Pa
      8-  criticalTemperature => critical temperature in K
      9-  criticalVolume => critical volume in m3/mol
      10- diffusionVolume => diffusion volume in m3
      11- dipoleMoment => dipole moment in Cm
      12- energyLennardJones => Lennard-Jones energy parameter in K (divided by Boltzmann constant)
      13- gyrationRadius => radius of gyration in m
      14- heatOfFusionAtNormalFreezingPoint => enthalpy change on melting at normal freezing point (101325 Pa) in J/mol
      15- heatOfVaporizationAtNormalBoilingPoint => enthalpy change on vaporization at normal boiling point (101325 Pa) in J/mol
      16- idealGasGibbsFreeEnergyOfFormationAt25C => in J/mol
      17- liquidDensityAt25C => liquid density at 25 ºC in mol/m3
      18- liquidVolumeAt25C => liquid volume at 25 ºC in m3/mol
      19- lengthLennardJones => Lennard-Jones length parameter in m
      20- molecularWeight => relative molar mass
      21- normalBoilingPoint => boiling point temperature at 101325 Pa in K
      22- normalFreezingPoint => melting point temperature at 101325 Pa in K
      23- parachor => Parachor in m3 kg0.25/(s0.5 mol)
      24- standardEntropyGas => Standard entropy of gas in J/mol
      25- standardEntropyLiquid => standard entropy of liquid in J/mol
      26- standardEntropySolid => standard entropy of solid in J/mol
      27- standardEnthalpyAqueousDilution => Standard aqueous infinite dilution enthalpy in J/mol
      28- standardFormationEnthalpyGas => standard enthalpy change on formation of gas in J/mol
      29- standardFormationEnthalpyLiquid => standard enthalpy change on formation of liquid in J/mol
      30- standardFormationEnthalpySolid => standard enthalpy change on formation of solid in J/mol
      31- standardFormationGibbsEnergyGas => standard Gibbs energy change on formation of gas in J/mol
      32- standardFormationGibbsEnergyLiquid => standard Gibbs energy change on formation of liquid in J/mol
      33- standardFormationGibbsEnergySolid => standard Gibbs energy change on formation of solid in J/mol
      34- standardGibbsAqueousDilution => Standard aqueous infinite dilution Gibbs energy in J/mol
      35- triplePointPressure => triple point pressure in Pa
      36- triplePointTemperature => triple point temperature in K
      37- vanderwaalsArea => van der Waals area in m3/mol
      38- vanderwaalsVolume => van der Waals volume in m3/mol
  """
  constantfloats=ASCIIString[
    "acentricFactor", # Pitzer acentric factor
    "associationParameter", # association-parameter (Hayden-O’Connell)
    "bornRadius", # in m
    "charge",
    "criticalCompressibilityFactor", # critical compressibility factor Z
    "criticalDensity", # critical density in mol/m3
    "criticalPressure", # critical pressure in Pa
    "criticalTemperature", # critical temperature in K
    "criticalVolume", # critical volume in m3/mol
    "diffusionVolume", # diffusion volume in m3
    "dipoleMoment", # dipole moment in Cm
    "energyLennardJones", # Lennard-Jones energy parameter in K (divided by Boltzmann constant)
    "gyrationRadius", # radius of gyration in m
    "heatOfFusionAtNormalFreezingPoint", # enthalpy change on melting at normal freezing point (101325 Pa) in J/mol
    "heatOfVaporizationAtNormalBoilingPoint", # enthalpy change on vaporization at normal boiling point (101325 Pa) in J/mol
    "idealGasGibbsFreeEnergyOfFormationAt25C", # in J/mol
    "liquidDensityAt25C", # liquid density at 25 ºC in mol/m3
    "liquidVolumeAt25C", # liquid volume at 25 ºC in m3/mol
    "lengthLennardJones", # Lennard-Jones length parameter in m
    "molecularWeight", # relative molar mass
    "normalBoilingPoint", # boiling point temperature at 101325 Pa in K
    "normalFreezingPoint", # melting point temperature at 101325 Pa in K
    "parachor", # Parachor in m3 kg0.25/(s0.5 mol)
    "standardEntropyGas", # Standard entropy of gas in J/mol
    "standardEntropyLiquid", # standard entropy of liquid in J/mol
    "standardEntropySolid", # standard entropy of solid in J/mol
    "standardEnthalpyAqueousDilution", # Standard aqueous infinite dilution enthalpy in J/mol
    "standardFormationEnthalpyGas", # standard enthalpy change on formation of gas in J/mol
    "standardFormationEnthalpyLiquid", # standard enthalpy change on formation of liquid in J/mol
    "standardFormationEnthalpySolid", # standard enthalpy change on formation of solid in J/mol
    "standardFormationGibbsEnergyGas", # standard Gibbs energy change on formation of gas in J/mol
    "standardFormationGibbsEnergyLiquid", # standard Gibbs energy change on formation of liquid in J/mol
    "standardFormationGibbsEnergySolid", # standard Gibbs energy change on formation of solid in J/mol
    "standardGibbsAqueousDilution", # Standard aqueous infinite dilution Gibbs energy in J/mol
    "triplePointPressure", # triple point pressure in Pa
    "triplePointTemperature", # triple point temperature in K
    "vanderwaalsArea", # van der Waals area in m3/mol
    "vanderwaalsVolume" # van der Waals volume in m3/mol
  ]

  """
    Temprature dependent properties
      1-  cpAqueousInfiniteDilution => Heat capacity of a solute in an infinitely dilute aqueous solution. in J/(mol K)
      2-  dielectricConstant => The ratio of the capacity of a condenser with a particular substance as dielectric to the capacity of the same condenser with a vacuum for dielectric.
      3-  expansivity => Coefficient of linear expansion for a solid at 1 atm in 1/K
      4-  fugacityCoefficientOfVapor => Fugacity coefficient of vapour on the saturation line
      5-  glassTransitionPressure => Glass transition pressure in Pa
      6-  heatCapacityOfLiquid => Heat capacity (Cp) of liquid on the saturation line in J/(mol K)
      7-  heatCapacityOfSolid => Solid heat capacity (Cp) at 1 atm in J/(mol K)
      8-  heatOfFusion => Enthalpy change on fusion fot the solid on the melting line in J/mol
      9-  heatOfSublimation => Enthalpy change on evaporation of the solid on the sublimation line in J/mol
      10- heatOfSolidSolidPhaseTransition => Enthalpy change on phase transition in J/mol
      11- heatOfVaporization => Enthalpy change on evaporation of the liquid on the saturation line in J/mol
      12- idealGasEnthalpy => Enthalpy of ideal gas in J/mol
      13- idealGasEntropy => Temperature-dependent part of entropy of ideal gas in J/(mol K)
      14- idealGasHeatCapacity => Heat capacity (Cp) of ideal gas in J/(mol K)
      15- meltingPressure => Pressure on melting line in Pa
      16- selfDiffusionCoefficientGas => Self-diffusion coefficient in gas phase at 1 atm in m2/s
      17- selfDiffusionCoefficientLiquid => self-diffusion coefficient in liquid phase on saturation line in m2/s
      18- solidSolidPhaseTransitionPressure => Pressure at phase transition in Pa
      19- sublimationPressure => Vapour pressure of solid on the sublimation line in Pa
      20- surfaceTensionSatLiquid => Surface tension of liquid on the saturation line in N/m
      21- thermalConductivityOfLiquid => Thermal conductivity of liquid on saturation line in W/(m K)
      22- thermalConductivityOfSolid => Thermal conductivity of solid at 1 atm in W/(m K)
      23- thermalConductivityOfVapor => Thermal conductivity of dilute gas in W/(m K)
      24- vaporPressure => Vapour pressure of saturated liquid in Pas
      25- virialCoefficient => Second virial coefficient of gas in m3/mol
      26- viscosityOfLiquid => Viscosity of liquid on saturation line in Pas
      27- viscosityOfVapor => Viscosity in dilute gas state in Pas
      28- volumeChangeUponMelting => Volume change for the solid on the melting line in m3/mol
      29- volumeChangeUponSolidSolidPhaseTransition => Volume change upon solid-solid phase transition in m3/mol
      30- volumeChangeUponSublimation => Volume change for the solid on the sublimation line in m3/mol
      31- volumeChangeUponVaporization => Volume change for the liquid on the saturation line in m3/mol
      32- volumeOfLiquid => Volume of liquid on saturation line in m3/mol
      33- volumeOfSolid => Volume of solid at 1 atm in m3/mol

  """
  tempreturedependents=ASCIIString[
    "cpAqueousInfiniteDilution", # Heat capacity of a solute in an infinitely dilute aqueous solution. in J/(mol K)
    "dielectricConstant", # The ratio of the capacity of a condenser with a particular substance as dielectric to the capacity of the same condenser with a vacuum for dielectric.
    "expansivity", # Coefficient of linear expansion for a solid at 1 atm in 1/K
    "fugacityCoefficientOfVapor", # Fugacity coefficient of vapour on the saturation line
    "glassTransitionPressure", # Glass transition pressure in Pa
    "heatCapacityOfLiquid", # Heat capacity (Cp) of liquid on the saturation line in J/(mol K)
    "heatCapacityOfSolid", # Solid heat capacity (Cp) at 1 atm in J/(mol K)
    "heatOfFusion", # Enthalpy change on fusion fot the solid on the melting line in J/mol
    "heatOfSublimation", # Enthalpy change on evaporation of the solid on the sublimation line in J/mol
    "heatOfSolidSolidPhaseTransition", # Enthalpy change on phase transition in J/mol
    "heatOfVaporization", # Enthalpy change on evaporation of the liquid on the saturation line in J/mol
    "idealGasEnthalpy", # Enthalpy of ideal gas in J/mol
    "idealGasEntropy", # Temperature-dependent part of entropy of ideal gas in J/(mol K)
    "idealGasHeatCapacity", # Heat capacity (Cp) of ideal gas in J/(mol K)
    "meltingPressure", # Pressure on melting line in Pa
    "selfDiffusionCoefficientGas", # Self-diffusion coefficient in gas phase at 1 atm in m2/s
    "selfDiffusionCoefficientLiquid", # self-diffusion coefficient in liquid phase on saturation line in m2/s
    "solidSolidPhaseTransitionPressure", # Pressure at phase transition in Pa
    "sublimationPressure", # Vapour pressure of solid on the sublimation line in Pa
    "surfaceTensionSatLiquid", # Surface tension of liquid on the saturation line in N/m
    "thermalConductivityOfLiquid", # Thermal conductivity of liquid on saturation line in W/(m K)
    "thermalConductivityOfSolid", # Thermal conductivity of solid at 1 atm in W/(m K)
    "thermalConductivityOfVapor", # Thermal conductivity of dilute gas in W/(m K)
    "vaporPressure", # Vapour pressure of saturated liquid in Pas
    "virialCoefficient", # Second virial coefficient of gas in m3/mol
    "viscosityOfLiquid", # Viscosity of liquid on saturation line in Pas
    "viscosityOfVapor", # Viscosity in dilute gas state in Pas
    "volumeChangeUponMelting", # Volume change for the solid on the melting line in m3/mol
    "volumeChangeUponSolidSolidPhaseTransition", # Volume change upon solid-solid phase transition in m3/mol
    "volumeChangeUponSublimation", # Volume change for the solid on the sublimation line in m3/mol
    "volumeChangeUponVaporization", # Volume change for the liquid on the saturation line in m3/mol
    "volumeOfLiquid", # Volume of liquid on saturation line in m3/mol
    "volumeOfSolid" # Volume of solid at 1 atm in m3/mol
  ]
  """
    Pressure dependent properties
      1- boilingPointTemperature => Temperature at liquid-vapour transition in K
      2- glassTransitionTemperature => Glass transition temperature in K
      3- meltingTemperature => Temperature on melting line in K
      4- solidSolidPhaseTransitionTemperature => Temperature at phase transition in K
  """
  pressuredependents=ASCIIString[
    "boilingPointTemperature", # Temperature at liquid-vapour transition in K
    "glassTransitionTemperature", # Glass transition temperature in K
    "meltingTemperature", # Temperature on melting line in K
    "solidSolidPhaseTransitionTemperature"  # Temperature at phase transition in K
  ]

  """
    Material Object refers to the software object that implements the
    representation of a Material.
    Material Object must:
    * Allow a client to set any basis-dependent Physical Property on any basis.
    * Allow a client to get any basis-dependent Physical Property using the basis with which it was stored.
    * Perform basis conversions, or delegate basis conversion as necessary. If basis conversion is not meaningful (e.g. in the case of cement), the Material Object must be able to return the quantity in its original basis and to return an error should the quantity be requested in a different basis.
    * Ensure that quantities set in one basis are consistent with quantities set in another basis, or delegate that function as necessary. Where the basis conversion on a quantity is not  feasible, the Material Object must only store the quantity in the basis with which it was set most recently
  """
  type MaterialObject
  end

  typealias PropertyMap Dict{ASCIIString,ASCIIString}
  propmap(stvec,t)=[cs=>t[ind] for (ind,cs) in enumerate(stvec)]

  """
    Property Package – a software component that is both a Physical Property Calculator and
    an Equilibrium Calculator for Materials containing a specific set of Compounds
    occurring in a specific number of physical states. A Property Package will make use of
    certain models to perform these calculations. A Property Package can be configured to
    make use of external Physical Property Calculators and/or Equilibrium Calculators.
    Alternatively, it can provide the functionality of these two components internally without
    making use of external components. Configuring the Compounds, Phases, models and
    external components used in a Property Package is outside the scope of this CAPEOPEN interface specification.
  """
  type PropertyPackage
    """
      Convert from perry units to cape-open standard units, for other properties default value of 1.0 is used
      6-  criticalDensity =>  mol/m3 mol/dm3 1E+3
      7-  criticalPressure =>  Pa MPa 1E+6
      9-  criticalVolume => m3/mol m3/Kmol 1E-3
      15- heatOfVaporizationAtNormalBoilingPoint => J/mol J/kmol 1E-3
      16- idealGasGibbsFreeEnergyOfFormationAt25C => J/mol J/kmol 1E-3
      18- liquidVolumeAt25C => m3/mol dm3/mol 1E-3
      24- standardEntropyGas => J/mol J/Kmol  1E-3
      28- standardFormationEnthalpyGas => J/mol  J/Kmol 1E-3
      31- standardFormationGibbsEnergyGas => J/mol J/Kmol 1E-3
      6-  heatCapacityOfLiquid => J/(mol K) J/(Kmol K) 1E-3
      11- heatOfVaporization => J/mol J/kmol 1E-3
      12- idealGasEnthalpy => J/mol J/Kmol 1E-3
      13- idealGasEntropy =>  J/(mol K) J/(Kmol K) 1E-3
      14- idealGasHeatCapacity =>  J/(mol K) J/(Kmol K) 1E-3
      32- volumeOfLiquid => m3/mol dm3/mol 1E-3
    """
    function PropertyPackage(constantstrings,t1,constantfloats,t2,tempreturedependents,t3,pressuredependents,t4,propertytable)
      new(
        propmap(constantstrings,t1),
        propmap(constantfloats,t2),
        propmap(tempreturedependents,t3),
        propmap(pressuredependents,t4),
        propertytable,
        Dict(
          "criticalDensity"                         =>1E+3, #mol/m3 <= mol/dm3
          "criticalPressure"                        =>1E+6, #Pa <= MPa
          "criticalVolume"                          =>1E-3, #m3/mol <= m3/Kmol
          "heatOfVaporizationAtNormalBoilingPoint"  =>1E-3, #J/mol <= J/kmol
          "idealGasGibbsFreeEnergyOfFormationAt25C" =>1E-3, #J/mol <= J/kmol
          "liquidVolumeAt25C"                       =>1E-3, #m3/mol <= dm3/mol
          "standardEntropyGas"                      =>1E-3, #J/mol <= J/Kmol /1E3
          "standardFormationEnthalpyGas"            =>1E-3, #J/mol <=  J/Kmol /1E3
          "standardFormationGibbsEnergyGas"         =>1E-3, #J/mol <= J/Kmol /1E3
          "heatCapacityOfLiquid"                    =>1E-3, #J/(mol K) <= J/(Kmol K) /1E3
          "heatOfVaporization"                      =>1E-3, #J/mol <= J/kmol /1E3
          "idealGasEnthalpy"                        =>1E-3, #J/mol <= J/Kmol /1E3
          "idealGasEntropy"                         =>1E-3, #J/(mol K) <= J/(Kmol K)
          "idealGasHeatCapacity"                    =>1E-3, #J/(mol K) <= J/(Kmol K)
          "volumeOfLiquid"                          =>1E-3 #m3/mol <= dm3/mol
        ),
        Dict{ASCIIString, Array}(
          "Criticals"=>       Any[[0,0,0], Vector{Float64}(6)],
          "LiquidsDensities"=>Any[[0,0,0], Vector{Float64}(10), 0%UInt8],
          "FormationEnergy"=> Any[[0,0,0], Vector{Float64}(5)],
          "Compounds"=>       Any[[0,0,0], Vector{UInt8}(29), Vector{UInt8}(9), Vector{UInt8}(11), 0.0, 0.0],
          "Cp"=>              Any[[0,0,0], Vector{Float64}(10), 0%UInt8],
          "LiquidsCp"=>       Any[[0,0,0], Vector{Float64}(10), 0%UInt8],
          "VaporizHeat"=>     Any[[0,0,0], Vector{Float64}(10)],
          "LiquidThermal"=>   Any[[0,0,0], Vector{Float64}(10)],
          "VaporPressure"=>   Any[[0,0,0], Vector{Float64}(10)],
          "LiquidViscos"=>    Any[[0,0,0], Vector{Float64}(10)],
          "VaporThermal"=>    Any[[0,0,0], Vector{Float64}(9)],
          "VaporViscos"=>     Any[[0,0,0], Vector{Float64}(9)]
        ),
        "perryanalytic"
      )
    end
    constantstrings::PropertyMap
    constantfloats::PropertyMap
    tempreturedependents::PropertyMap
    pressuredependents::PropertyMap
    propertytable::Vector{ASCIIString}
    convertions::Dict{ASCIIString,Float64}
    tableaddreses::Dict{ASCIIString, Array}
    datafilename::ASCIIString
  end
  perryanalytic=PropertyPackage(
    getindex(constantstrings,1:3),
    ["Compounds" for i=1:3],
    getindex(constantfloats,[ 5,6,7,8,9,15,16,17,18,
                                20,21,24,28,31]),
    ["Criticals", "Criticals", "Criticals", "Criticals", "Criticals", "VaporizHeat", "FormationEnergy", "LiquidsDensities", "LiquidsDensities",
     "Compounds", "Compounds", "FormationEnergy", "FormationEnergy", "FormationEnergy"],
    getindex(tempreturedependents,[6,11,12,13,14,21,23,24,26,27,32]),
    ["LiquidsCp", "VaporizHeat", "Cp", "Cp", "Cp", "LiquidThermal", "VaporThermal", "VaporPressure", "LiquidViscos", "VaporViscos", "LiquidsDensities"],
    ["boilingPointTemperature"],
    ["VaporPressure"],
    ASCIIString[
      "Compounds"        ,
      "VaporPressure"    ,
      "LiquidsDensities",
      "Criticals"        ,
      "VaporizHeat"      ,
      "LiquidsCp"        ,
      "Cp"                ,
      "FormationEnergy" ,
      "VaporViscos"      ,
      "LiquidViscos"     ,
      "VaporThermal"     ,
      "LiquidThermal"
    ]
  )
end


text:
@testset "test_utils.jl" begin

    estimates = (PointEstimate, DistributionEstimate)
    outputs = (SingleOutput, MultiOutput)

    @testset "$est, $out, PointPredictInput" for (est, out) in Iterators.product(estimates, outputs)
        temp = FakeTemplate{est, out}()
        test_interface(temp)
    end

    @testset "Vector inputs case" begin
        temp = FakeTemplate{PointEstimate, SingleOutput}()
        test_interface(temp; inputs=[rand(5), rand(5)], outputs=rand(1, 2))
    end

    @testset "$est, $out, PointOrDistributionPredictInput" for (est, out) in Iterators.product(estimates, outputs)
        Models.predict_input_type(m::Type{<:FakeTemplate}) = PointOrDistributionPredictInput
        Models.predict_input_type(m::Type{<:FakeModel}) = PointOrDistributionPredictInput

        temp = FakeTemplate{est, out}()
        test_interface(temp)
    end

    @testset "Vector inputs case" begin
        temp = FakeTemplate{PointEstimate, SingleOutput}()
        test_interface(
            temp;
            inputs=[rand(5), rand(5)],
            outputs=rand(1, 2),
            distribution_inputs=[MvNormal(5, m) for m in 1:2]
        )
    end

end


text:
include("roblib.jl")
function main(r::Robot)
    path = go_to_west_south_corner_and_return_path!(r; go_around_barriers = true)

    is_marker_now = iseven(get_path_length(path))

    side = East
    while !(isborder(r, North) && isborder(r, side))
        if (is_marker_now)
            putmarker!(r)
        end

        if !(isborder(r, side))
            move!(r,side)
            is_marker_now = !is_marker_now
        end

        if move_up_condition(r)
            if (is_marker_now)
                putmarker!(r)
            end
            move!(r, North)
            is_marker_now = !is_marker_now

            side = inverse_side(side)
        end
    end
    if is_marker_now
        putmarker!(r)
    end

    go_to_west_south_corner_and_return_path!(r)
    go_by_path!(r, path)

end

function move_up_condition(r::Robot)::Bool
    return (isborder(r, East) || isborder(r, West)) && !(isborder(r,North))
end

r=Robot( animate=true)

max_stars_repo_name: Savko-fokin/mirea-progs

max_issues_repo_name: Savko-fokin/mirea-progs

max_forks_repo_name: Savko-fokin/mirea-progs


text:
using Optim

struct ValueIterationLog
    header::Array{String, 1}
    keywords::Array{Symbol, 1}
    entries::Array{Any, 1}
end

function ValueIterationLog()
    header = [ "It", "ηₙ=|xₙ-xₙ₋₁|", "νₙ=|vₙ-vₙ₋₁|", "Time", "Eval steps"]
    keywords = [:it, :sa_x, :sa_v, :time, :nit]
    ValueIterationLog(header, keywords, [])
end

mutable struct ValueIterationResult
    dr::AbstractDecisionRule
    drv::AbstractDecisionRule
    iterations::Int
    complementarities::Bool
    dprocess::AbstractDiscretizedProcess
    x_converged::Bool
    x_tol::Float64
    x_err::Float64
    v_converged::Bool
    v_tol::Float64
    v_err::Float64
    log::ValueIterationLog
    trace::Union{Nothing,IterationTrace}
end

converged(r::ValueIterationResult) = r.x_converged && r.v_converged
function Base.show(io::IO, r::ValueIterationResult)
    @printf io "Results of Value Iteration Algorithm\n"
    @printf io " * Complementarities: %s\n" string(r.complementarities)
    @printf io " * Decision Rule type: %s\n" string(typeof(r.dr))
    @printf io " * Discretized Process type: %s\n" string(typeof(r.dprocess))
    @printf io " * Number of iterations: %s\n" string(r.iterations)
    @printf io " * Convergence: %s\n" converged(r)
    @printf io "   * |x - x'| < %.1e: %s\n" r.x_tol r.x_converged
    @printf io "   * |v - v'| < %.1e: %s\n" r.v_tol r.v_converged

end

function append!(log::ValueIterationLog; verbose=true, entry...)
    d = Dict(entry)
    push!(log.entries, d)
    verbose && show_entry(log, d)
end

function initialize(log::ValueIterationLog; verbose=true)
    verbose && show_start(log)
end

function finalize(log::ValueIterationLog; verbose=true)
    verbose && show_end(log)
end

function show(log::ValueIterationLog)
    show_start(log)
    for entry in log.entries
        show_entry(log,entry)
    end
    show_end(log)
end

function show_start(log::ValueIterationLog)
    println(repeat("-", 66))
    @printf "%-6s%-16s%-16s%-16s%-5s\n" "It" "|xₙ-xₙ₋₁|" "|vₙ-vₙ₋₁|" "Time" "Eval steps"
    println(repeat("-", 66))
end

function show_entry(log::ValueIterationLog, entry)
    it = entry[:it]
    sa_x = entry[:sa_x]
    sa_v = entry[:sa_v]
    time = entry[:time]
    nit = entry[:nit]
    @printf "%-6i%-16.2e%-16.2e%-16.2e%-5i\n" it sa_x sa_v time nit
end

function show_end(log::ValueIterationLog)
    println(repeat("-", 66))
end

"""
Evaluate the value function under the given decision rule.

# Arguments
* `model::NumericModel`: Model object that describes the current model environment.
* `dr`: Current guess for the decision rule.
# Returns
* `drv`: Value function.
"""
function evaluate_policy(model, dprocess::AbstractDiscretizedProcess, grid, dr;
                            verbose::Bool=true, maxit::Int=5000, tol=1e-6)

    # extract parameters
    p = SVector(model.calibration[:parameters]...)
    β = model.calibration.flat[:beta]

    # states today are the grid
    s = nodes(ListOfPoints, grid)

    # Number of endogenous nodes
    N = length(s)

    # number of smooth decision rules
    nsd = max(n_nodes(dprocess), 1)

    n_u = length(model.calibration[:rewards])

    res = [zeros(Point{n_u},N) for i=1:nsd]

    # Value function : v_t = u + β*E[V_{t+1}]
    u = deepcopy(res)
    v = deepcopy(res)
    v0 = deepcopy(res)

    # Expected utility
    E_V = deepcopy(res)

    # evaluate u at time t. This is constant because dr is constant within this
    # function
    x = [dr(i, s) for i=1:nsd]
    for i = 1:size(res, 1)
        m = node(Point, dprocess, i)
        u[i][:] = Dolo.felicity(model, m, s, x[i], p)
    end

    # Initial guess for the value function
    v0 = deepcopy(u)

    #Preparation for a loop
    err = 10.0
    it = 0

    drv = CachedDecisionRule(dprocess, grid, v0)

    verbose && @printf "%-6s%-12s\n" "It" "SA"
    verbose && println(repeat("-", 14))

    while err > tol && it < maxit
        it += 1
        # Interpolate v0 on the grid
        # Compute value function for each smooth decision rule
        for i = 1:nsd
            m = node(Point, dprocess, i)
            for (w, MM, j) in get_integration_nodes(dprocess,i)
                    M = SVector(MM...)
                     # Update the states
                     S = Dolo.transition(model, m, s, x[i], M, p)
                     # Update value function
                      E_V[i] += w*drv(i, j, S)
            end
        end

        err = 0.0
        for i in 1:length(v)
            v[i][:] = u[i] + β*E_V[i]
            err = max(err, maxabs(v[i] - v0[i]))
            copyto!(v0[i], v[i])
            E_V[i] *= 0.0
        end

        verbose && @printf "%-6i%-12.2e\n" it err

        set_values!(drv, v0)
    end

    return (drv, it)
end

function evaluate_policy(model, dr; grid=Dict(), kwargs...)
    grid = get_grid(model, options=grid)
    dprocess = discretize(model.exogenous)
    return evaluate_policy(model, dprocess, grid, dr; kwargs...)[1]

end

"""
Evaluate the right hand side of the value function at given values of states, controls, and exogenous variables.

# Arguments
* `model::NumericModel`: Model object that describes the current model environment.
* `β::Float64`: Value of the discount factor.
* `dprocess::`: Discretized exogenous process.
* `drv`: Current guess for the value function object.
* `i::Int64`: Index of node for exogeous variable(s).
* `s::ListOfPoints`: List of state variables.
* `x0::ListOfListOfPoints`: List of control variables.
* `p::Vector{Float64}`: Model parameters.
# Returns
* `E_V::`: Right hand side of the value function.
"""
function update_value(model, β::Float64, dprocess, drv, i, s::Point{d},
                      x0::Point{n_x}, p::Point{n_p}) where d where n_x where n_p
    m = node(Point,dprocess, i)
    N = length(s)
    E_V = Point{1}(0.0)
    for (w, M, j) in get_integration_nodes(Point,dprocess,i)
        # Update the states
        S = Dolo.transition(model, m, s, x0, M, p)
        E_V += w*drv(i, j, S)[1]
    end
    u = Dolo.felicity(model, m, s, x0, p)[1]
    V = u + β*E_V
    return V
end

# function update_value(model, β::Float64, dprocess, drv, i, s::Point{d},
#                       x0::Float64, p::Point{n_p}) where d where n_p
#     update_value(model, β, dprocess, drv, i, s, SVector{1,Float64}(x0), p)
# end
#

function call_optim(fobj, initial_x, lower, upper, optim_opts)
    if length(initial_x) == 1
        return optimize(fobj, lower[1], upper[1])
    end
    return optimize(fobj, lower, upper, initial_x)
end

"""
Solve for the value function and associated decision rule using value function iteration.

# Arguments
* `model::NumericModel`: Model object that describes the current model environment.
* `pdr`: Initial guess for the decision rule.
# Returns
* `dr`: Solved decision rule object.
* `drv`: Solved value function object.
"""
function value_iteration(
        model, dprocess::AbstractDiscretizedProcess, grid, init_dr;
        discount_symbol=:beta,
        maxit::Int=1000, tol_x::Float64=1e-8, tol_v::Float64=1e-8,
        optim_options=Dict(), eval_options=Dict(),
        verbose::Bool=true, trace::Bool=false
    )

    β = model.calibration.flat[discount_symbol]

    # compute the value function

    p = SVector(model.calibration[:parameters]...)

    endo_nodes = nodes(ListOfPoints,grid)

    # Number of endogenous nodes
    N = length(endo_nodes)

    # number of smooth decision rules
    n_m = nsd = max(n_nodes(dprocess), 1)

    res = [zeros(Point{1},N) for i=1:nsd]

    # States at time t+1
    S = copy(endo_nodes)

    # Controls at time t
    x = [init_dr(i, endo_nodes) for i=1:nsd]
    x0 = deepcopy(x)

    dr = CachedDecisionRule(dprocess, grid, x0)

    drv = dr # never used. Just to escape loop scope
    v0 = [zeros(Point{1},N) for i=1:nsd]
    v = deepcopy(v0)

    ti_trace = trace ? IterationTrace([x0, v0]) : nothing

    # fourth time I repeat that one...
    n_x = length(model.calibration[:controls])
    x_lb = [controls_lb(model, node(Point,dprocess,i),endo_nodes,p) for i=1:n_m]
    x_ub = [controls_ub(model, node(Point,dprocess,i),endo_nodes,p) for i=1:n_m]
    BIG = 100000
    for i=1:n_m
      for n=1:N
        x_lb[i][n] = max.(x_lb[i][n],-BIG)
        x_ub[i][n] = min.(x_ub[i][n], BIG)
      end
    end

    #Preparation for a loop
    err_v = 10.0
    err_x = 10.0
    err_eval = 10.0

    it = 0
    it_eval = 0

    maxit_eval = get(eval_options, :maxit, 1000)
    tol_eval = get(eval_options, :tol, 1e-8)

    optim_opts = Optim.Options(optim_options...)

    log = ValueIterationLog()
    initialize(log, verbose=verbose)

    mode = :improve
    done = false

    while !done

        tic()
        it += 1

        drv, it_eval = evaluate_policy(model, dprocess, grid, dr; verbose=false, eval_options...)
        mode = :improve

        # else
        for i = 1:size(res, 1)
            m = node(dprocess, i)
            for n = 1:N
                s = endo_nodes[n]
                # optimize vals
                fobj(u) = -update_value(model, β, dprocess, drv, i, s, SVector(u...), p)[1]*1000
                lower = Float64[x_lb[i][n]...]
                upper = Float64[x_ub[i][n]...]
                upper = clamp.(upper, -Inf, 1000000)
                lower = clamp.(lower, -1000000, Inf)
                initial_x = Float64[x0[i][n]...] #, :]
                test = fobj(initial_x)
                results = call_optim(fobj, initial_x, lower, upper, optim_opts)
                xn = Optim.minimizer(results)
                nv = -Optim.minimum(results)/1000.0
                x[i][n] = SVector(xn...)
                v[i][n] = SVector(nv)
            end
        end
        # compute diff in values
        err_v = maxabs(v-v0)
        err_x = maxabs(x-x0)
        for i in 1:nsd
            copyto!(v0[i], v[i])
            copyto!(x0[i], x[i])
        end

        # update values and policies
        set_values!(drv, v0)
        set_values!(dr, x0)

        # terminate only if policy didn't move
        done = (err_x<tol_x) || (err_v<tol_v) || (it>=maxit)

        trace && push!(ti_trace.trace, [x0, v0]) # this is ridiculous since only one of them is updated !

        elapsed = toq()

        append!(log; verbose=verbose, it=it, sa_x=err_x, sa_v=err_v, time=elapsed, nit=it_eval)

    end

    finalize(log, verbose=verbose)

    converged_x = err_x<tol_x
    converged_v = err_v<tol_v

    ValueIterationResult(dr.dr, drv.dr, it, true, dprocess, converged_x, tol_x, err_x, converged_v, tol_v, err_v, log, ti_trace)

end

# get stupid initial rule
function value_iteration(model, dprocess::AbstractDiscretizedProcess, init_dr; grid=Dict(), kwargs...)
    grid = get_grid(model, options=grid)
    return value_iteration(model, dprocess, grid, init_dr;  kwargs...)
end

# get stupid initial rule
function value_iteration(model, dprocess::AbstractDiscretizedProcess; grid=Dict(), kwargs...)
    grid = get_grid(model, options=grid)
    init_dr = ConstantDecisionRule(model.calibration[:controls])
    return value_iteration(model, dprocess, grid, init_dr;  kwargs...)
end

function value_iteration(model, init_dr; grid=Dict(), kwargs...)
    grid = get_grid(model, options=grid)
    dprocess = discretize( model.exogenous )
    return value_iteration(model, dprocess, grid, init_dr; kwargs...)
end

function value_iteration(model; grid=Dict(), kwargs...)
    grid = get_grid(model; options=grid)
    dprocess = discretize( model.exogenous )
    init_dr = ConstantDecisionRule(model.calibration[:controls])
    return value_iteration(model, dprocess, grid, init_dr; kwargs...)
end

# compatibility
const solve_policy = value_iteration


text:
"""
Module providing extensions of linalg to BigFloat and Complex{BigFloat}
"""
module LinAlgBigFloat

using LinearAlgebra

include("util.jl")
#include("mpfrmodify.jl")
include("deflationcrit.jl")
include("refineprecision.jl")
include("separate.jl")
include("transformhess.jl")
#include("hessenbergbf.jl")
include("schurbf.jl")
include("sylvester.jl")
include("spaces.jl")
include("jordan.jl")
include("matrices.jl")
include("randommatrices.jl")

export is_hessenberg

end # module

max_stars_repo_name: KlausC/LinAlgBigFloat

max_issues_repo_name: KlausC/LinAlgBigFloat

max_forks_repo_name: KlausC/LinAlgBigFloat


text:
#-----------------------------------------------------------------------------------------------
# 1st run - packages installation

import Pkg;
Pkg.add("CSV")
Pkg.add("DataFrames")
Pkg.add("StatsBase")

file_path = "../../../../../../externals/Core.Math.Data/data/HolisticWare.Core.Sample.Data/CSV/Pejcic_318.csv"

using DataFrames, CSV, Statistics, StatsBase

data = CSV.read(file_path, copycols = true)

# @show data.ATV
# println()

for col in eachcol(data,true);
    if  (
        eltype(col[2]) <: Real)
        println(mean(col[2])
        );
    end;
end;


text:
  # Some utilities to modify nested QuoteNode's

  "takes :A,:B,:C and return :(A.B.C)"
  qn(a,b) = Expr(:(.),a,QuoteNode(b))
  qn(a,b,rest...) = qn(qn(a,b), rest...)

  "parse :(A.B.C) and returns an array of symbols with :A,:B,:C"
  qs(ex::Expr) = qs(ex,Symbol[])
  qs(s::Symbol) = [s]
  function qs(ex::Expr,out::Vector{Symbol})

      if ex.head == :(.)
          qs(ex.args[1],out)
          qs(ex.args[2],out)
      end
      out
  end
  qs(s::Symbol,out::Vector{Symbol}) = push!(out,s)
  qs(s::QuoteNode,out::Vector{Symbol}) = qs(s.value,out)

  #qn(:A, qs(:(B.C))...) == :(A.B.C)
  #qn(:A,:B,:C) == :(A.B.C)

  """
      eval_symbol(s, eval_in::Module)

  eval s in module eval_in, used for data hint.

  """
  function eval_symbol(ex::Union{Expr,Symbol},eval_in::Module)

      #this prepends eval_in so we can evaluate in Main directly :(eval_in.ex)
      ex = qn(eval_in, qs(ex)...)

      evalout = try eval(Main,ex) catch err "" end
  end

  function get_doc(s::Symbol,eval_in::Module)
      Base.Docs.doc(
          Base.Docs.Binding(eval_in,s)
      )
  end
  function get_doc(ex::Expr,eval_in::Module)
      try eval(eval_in,:( @doc $ex )) catch "" end
  end


text:
using SatisfiabilityInterface
using Symbolics
using Symbolics: Sym

function all_different(v)
    vv = vec(v)
    n = length(vv)

    return [vv[i] ≠ vv[j] for i in 1:n for j in i+1:n]
end

function make_matrix(name, m, n)
    return [Num(Sym{Real}(name, i, j)) for i in 1:m, j in 1:n]
end

n = 9
M = make_matrix(:M, n, n)

## initial condition from https://www.juliaopt.org/notebooks/JuMP-Sudoku.html

initial = [
    5 3 -1 -1 7 -1 -1 -1 -1
    6 -1 -1 1 9 5 -1 -1 -1
    -1 9 8 -1 -1 -1 -1 6 -1
    8 -1 -1 -1 6 -1 -1 -1 3
    4 -1 -1 8 -1 3 -1 -1 1
    7 -1 -1 -1 2 -1 -1 -1 6
    -1 6 -1 -1 -1 -1 2 8 -1
    -1 -1 -1 4 1 9 -1 -1 5
    -1 -1 -1 -1 8 -1 -1 7 9
    ]

constraints = [
    [M[i, j] ∈ 1:9 for i in 1:n for j in 1:n]

    [all_different(M[:, i]) for i in 1:n]
    [all_different(M[i, :]) for i in 1:n]

    [all_different(M[ 3i+1 : 3i+3, 3j+1 : 3j+3 ]) for i in 0:2 for j in 0:2]

    [M[i, j] == initial[i, j] for i in 1:9 for j in 1:9 if initial[i, j] > 0]

] |> Iterators.flatten |> collect

prob = DiscreteCSP(constraints)

@time status, results = solve(prob);

prob.varmap

results

output = reshape([last(x) for x in results], 9, 9)'

output == [
    5  3  4  6  7  8  9  1  2
    6  7  2  1  9  5  3  4  8
    1  9  8  3  4  2  5  6  7
    8  5  9  7  6  1  4  2  3
    4  2  6  8  5  3  7  9  1
    7  1  3  9  2  4  8  5  6
    9  6  1  5  3  7  2  8  4
    2  8  7  4  1  9  6  3  5
    3  4  5  2  8  6  1  7  9
]

all_solutions(prob)


text:
using BinaryBuilder

name = "pigz"
version = v"2.7"

# Collection of sources required to build pigz
sources = [
    ArchiveSource("https://zlib.net/pigz/pigz-$(version.major).$(version.minor).tar.gz",
                  "b4c9e60344a08d5db37ca7ad00a5b2c76ccb9556354b722d56d55ca7e8b1c707"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/pigz-*

export CPPFLAGS="-I${includedir}"
if [[ "${target}" == *-freebsd* ]]; then
    # Without this, `__XSI_VISIBLE`, `S_IFMT`, `S_IFREG`, `S_IFMT`,
    # `S_IFREG`, `S_IFIFO`, and `S_IFDIR` are all undefined
    CPPFLAGS="${CPPFLAGS} -D_XOPEN_SOURCE=700"
fi
make -j${nproc} CC=${CC}

# Install
mkdir -p ${bindir}
for bin in pigz unpigz; do
    cp "${bin}" "${bindir}/${bin}${exeext}"
done

# License is embedded at the end of the README
install_license README
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    ExecutableProduct("pigz", :pigz),
    ExecutableProduct("unpigz", :unpigz),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Zlib_jll"),
]

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

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
module JITrench

include("utils/error.jl")

include("core/variable.jl")
include("core/functions_utils.jl")
include("core/function.jl")
include("core/operators.jl")
include("core/propagation.jl")
include("core/math_functions.jl")

include("arr/reshape.jl")
include("arr/flatten.jl")
include("arr/transpose.jl")
include("arr/sum.jl")
include("arr/sum_to.jl")
include("arr/broadcast_to.jl")
include("arr/broadcast.jl")
include("arr/matmul.jl")
include("arr/getindex.jl")

include("utils/out.jl")
include("utils/plot.jl")
include("utils/utils.jl")

include("nn/utils.jl")
include("nn/layers/layer.jl")
include("nn/layers/linear.jl")
include("nn/funcitons/activation/sigmoid.jl")
include("nn/funcitons/loss/mean_squared_error.jl")

include("nn/model.jl")

include("nn/optimizer/optimizer.jl")
include("nn/optimizer/sgd.jl")

export Variable, DiffableFunction, backward!, parameters, cleargrad!, flatten, matmul, sum_to, sigmoid, mean_squared_error, linear, Model, Layer, SGD, @diff!

end


text:
module TestDeprecated

using Test, DataFrames
import DataFrames: identifier

# old sort(df; cols=...) syntax
df = DataFrame(a=[1, 3, 2], b=[6, 5, 4])
@test sort(df; cols=[:a, :b]) == DataFrame(a=[1, 2, 3], b=[6, 4, 5])
sort!(df; cols=[:a, :b])
@test df == DataFrame(a=[1, 2, 3], b=[6, 4, 5])

df = DataFrame(a=[1, 3, 2], b=[6, 5, 4])
@test sort(df; cols=[:b, :a]) == DataFrame(a=[2, 3, 1], b=[4, 5, 6])
sort!(df; cols=[:b, :a])
@test df == DataFrame(a=[2, 3, 1], b=[4, 5, 6])

@test first.(collect(pairs(DataFrameRow(df, 1, :)))) == [:a, :b]
@test last.(collect(pairs(DataFrameRow(df, 1, :)))) == [df[1, 1], df[1, 2]]

@testset "identifier" begin
    @test identifier("%_B*_\tC*") == :_B_C_
    @test identifier("2a") == :x2a
    @test identifier("!") == :x!
    @test identifier("\t_*") == :_
    @test identifier("begin") == :_begin
    @test identifier("end") == :_end
end

# Check that reserved words are up to date

@testset "reserved words" begin
    f = "$(Sys.BINDIR)/../../src/julia-parser.scm"
    if isfile(f)
        r1 = r"define initial-reserved-words '\(([^)]+)"
        r2 = r"define \(parse-block s(?: \([^)]+\))?\)\s+\(parse-Nary s (?:parse-eq '\([^(]*|down '\([^)]+\) '[^']+ ')\(([^)]+)"
        body = read(f, String)
        m1, m2 = match(r1, body), match(r2, body)
        if m1 == nothing || m2 == nothing
            error("Unable to extract keywords from 'julia-parser.scm'.")
        else
            s = replace(string(m1.captures[1]," ",m2.captures[1]), r";;.*?\n" => "")
            rw = Set(split(s, r"\W+"))
            @test rw == DataFrames.RESERVED_WORDS
        end
    else
        @warn("Unable to validate reserved words against parser. ",
              "Expected if Julia was not built from source.")
    end
end

# deprecated combine

df = DataFrame(a=[1, 1, 2, 2, 2], b=1:5)
gd = groupby(df, :a)
@test combine(gd) == combine(identity, gd)

end # module


text:
include("PosEncodeTest.jl")
include("../src/MultiHead.jl")

Layers = 8
@assert Int(EmbDims / Layers) - (EmbDims / Layers) == 0

MHA = MultiHeadAttention(EmbDims, Layers)
@time output = MHA(output)


text:
# Returns the header containing all instructions
# to use alphalib
alphalib_funcs = [
        "alphaprintint",
        "alphaprintstr",
        "alphascan",
    ]

function alphalib_head()
    s = ""
    for f in alphalib_funcs
        s *= "extern $f\n"
    end

    return s
end


text:
function mark_corners!(r)
    side1 = Nord
    side2 = West
    go_to_corner_with_walls!(r, side1, side2)
    for i in [Sud, Ost, Nord, West]
        putmarker_to_corner(i)
    end
    while length(sides) > 0
        move!(r, inverse_side(pop!(sides)))
    end
end

max_stars_repo_name: dmitrii-a-piskunov/mirea-progs

max_issues_repo_name: dmitrii-a-piskunov/mirea-progs

max_forks_repo_name: dmitrii-a-piskunov/mirea-progs


text:
using AXMethods
using Test

@testset "AXLeastSquares.jl" begin
    # Generate example data
    n_i = 1000
    X   = 2 * randn((n_i, 2))
    y   = X * [1, -2] + randn((n_i, 1))

    # Estimate the least square regression
    ls_fit = AXLeastSquares(y, X)

    # Check the methods
    β_hat = coef(ls_fit)
    y_hat = predict(ls_fit)
    se_hat = infer(ls_fit, print_df = false)
    se_robust_hat = infer(ls_fit, heteroskedastic = true, print_df = false)

    # Let's check that everything is of correct type.
    @test typeof(ls_fit) == AXLeastSquares
    @test typeof(β_hat) == Array{Float64,2}
    @test typeof(y_hat) == Array{Float64,2}
    @test typeof(se_hat) == NamedTuple{(:β, :se, :t, :p),
                                       Tuple{Array{Float64,2},
                                             Array{Float64,1},
                                             Array{Float64,2},
                                             Array{Float64,2}}}
    @test typeof(se_robust_hat) == NamedTuple{(:β, :se, :t, :p),
                                          Tuple{Array{Float64,2},
                                                Array{Float64,1},
                                                Array{Float64,2},
                                                Array{Float64,2}}}
end

@testset "AX2StageLeastSquares.jl" begin
    # Generate example data
    n_i = 1000
    ν = randn((n_i, 1))
    instrument = 2 * randn((n_i, 2))
    D = instrument * [1, -2] + ν
    y = D + 0.3 * ν + randn((n_i, 1))

    # Estimate the least square regression
    tsls_fit = AX2StageLeastSquares(y, D, instrument)

    # Check the methods
    β_hat = coef(tsls_fit)
    y_hat = predict(tsls_fit)
    se_hat = infer(tsls_fit, print_df = false)
    se_hc_hat = infer(tsls_fit, heteroskedastic = true, print_df = false)

    # Let's check that everything is of correct type.
    @test typeof(tsls_fit) == AX2StageLeastSquares
    @test typeof(β_hat) == Array{Float64,2}
    @test typeof(y_hat) == Array{Float64,2}
    @test typeof(se_hat) == NamedTuple{(:β, :se, :t, :p),
                                       Tuple{Array{Float64,2},
                                             Array{Float64,1},
                                             Array{Float64,2},
                                             Array{Float64,2}}}
    @test typeof(se_hc_hat) == NamedTuple{(:β, :se, :t, :p),
                                          Tuple{Array{Float64,2},
                                                Array{Float64,1},
                                                Array{Float64,2},
                                                Array{Float64,2}}}
end


text:
"""
    phase_diff(field, dir)

Compute point-to-point difference of phase on a grid along a direction

Returns an array of size `(size(field)[1], size(field)[2], size(field)[3])`
containing differences in direction `dir`.
"""
function phase_diff(field, dir)
    out = similar(field, Float64)

    out .= angle.(circshift(field, dir) ./  field)

    out
end

"""
    phase_diff(field::PencilArray, dir)

Compute point-to-point difference of phase on a grid along a direction

Returns an array of size `(size(field)[1], size(field)[2], size(field)[3])`
containing differences in direction `dir`.
"""
function phase_diff(field::PencilArray, dir)
    out = similar(field, Float64)

    # TODO diff at boundaries
    out .= 0
    out[1:end-dir[1], 1:end-dir[2], 1:end-dir[3]] .= angle.(circshift(field[1:end-dir[1], 1:end-dir[2], 1:end-dir[3]], dir) ./ field[1+dir[1]:end, 1+dir[2]:end, 1+dir[3]:end])

    out
end

"""
    normed_max_phase_grad(grids)

Compute maximum phase gradient of a grid

Normalised to ignore large gradients in regions with low density.  These tend
to be anomalous.
"""
function max_normed_phase_diff(psi, rho, density_threshold)
    tmp = similar(psi, Float64)
    n = ndims(tmp)

    max_grads = zeros(n)

    shift = zeros(Int, n)
    shift[1] = 1

    for i in 1:n
        dir = circshift(shift, i)
        tmp .= abs.(phase_diff(psi, dir))
        tmp[rho / maximum(rho) .< density_threshold] .= NaN
        max_grads[i] = maximum(tmp)
    end

    maximum(max_grads)

end

"""
    bad_phase_diff(grids, config)

Check if the phase of the ψ field is in a trustable regime.

Returns
"""
function good_phase_diff(grids, config)
    density_threshold = config.density_threshold
    if max_normed_phase_diff(grids.ψx, grids.ρx, density_threshold) > config.phase_grad_limit
        false
    else
        true
    end
end


text:
"""
Define your model to be a subtype of `AbstractModel`. Your model has to have the following fields, but can also have other fields of your choice.

e.g.

```
mutable struct MyModel <: AbstractModel
  scheduler::Function
  space
  agents::Array{Integer}  # a list of agents ids
end
```

`scheduler` can be one of the default functions (`random_activation`), or your own function.
"""
abstract type AbstractModel end

"""
  nagents(model::AbstractModel)

Returns the number of agents.
"""
nagents(model::AbstractModel) = length(model.agents)

function dummystep(a::AbstractAgent, b::AbstractModel)
end

"""
    kill_agent!(agent::AbstractAgent, model::AbstractModel)

Removes an agent from the list of agents and from the space.
"""
function kill_agent!(agent::AbstractAgent, model::AbstractModel)
  if typeof(agent.pos) <: Tuple
    agentnode = coord_to_vertex(agent.pos, model)
  else
    agentnode = agent.pos
  end
  splice!(model.space.agent_positions[agentnode], findfirst(a->a==agent.id, model.space.agent_positions[agentnode]))  # remove from the grid
  splice!(model.agents, findfirst(a->a==agent, model.agents))  # remove from the model.agents
end

"""
    step!(agent_step::Function, model::AbstractModel)

Updates agents one step. Agents will be updated as specified by the `model.scheduler`.
"""
function step!(agent_step, model::AbstractModel)
  activation_order = return_activation_order(model)
  for index in activation_order
    agent_step(model.agents[index], model)
  end
end

"""
    step!(agent_step::Function, model::AbstractModel, nsteps::Integer)

Repeats the `step` function `nsteps` times without collecting data.
"""
function step!(agent_step, model::AbstractModel, nsteps::Integer)
  for i in 1:nsteps
    step!(agent_step, model)
  end
end

"""
    step!(agent_step::Function, model::AbstractModel, nsteps::Integer, agent_properties::Array{Symbol}, steps_to_collect_data::Array{Int64})

Repeats the `step` function `nsteps` times, and collects all agent fields in `agent_properties` at steps `steps_to_collect_data`.

"""
function step!(agent_step::Function, model::AbstractModel, nsteps::Integer, agent_properties::Array{Symbol}, steps_to_collect_data::Array{Int64})

  # Run the first step of the model to fill in the dataframe
  # step!(agent_step, model)
  df = data_collector(agent_properties, steps_to_collect_data, model, 1)

  for ss in 2:nsteps
    step!(agent_step, model)
    # collect data
    if ss in steps_to_collect_data
      df = data_collector(agent_properties, steps_to_collect_data, model, ss, df)
    end
  end
  # if 1 is not in `steps_to_collect_data`, remove the first row.
  if !in(1, steps_to_collect_data)
    df = df[2:end, :]
  end
  return df
end

"""
    step!(agent_step::Function, model::AbstractModel, nsteps::Integer, agent_properties::Array{Symbol}, aggregators::Array, steps_to_collect_data::Array{Integer})

Repeats the `step` function `nsteps` times, and applies functions in `aggregators` to values of agent fields in `agent_properties` at steps `steps_to_collect_data`.
"""
function step!(agent_step, model::AbstractModel, nsteps::Integer, agent_properties::Array{Symbol}, aggregators::Array, steps_to_collect_data::Array{Int64})

  # Run the first step of the model to fill in the dataframe
  # step!(agent_step, model)
  df = data_collector(agent_properties, aggregators, steps_to_collect_data, model, 1)

  for ss in 2:nsteps
    step!(agent_step, model)
    # collect data
    if ss in steps_to_collect_data
      df = data_collector(agent_properties, aggregators, steps_to_collect_data, model, ss, df)
    end
  end
  # if 1 is not in `steps_to_collect_data`, remove the first row.
  if !in(1, steps_to_collect_data)
    df = df[2:end, :]
  end
  return df
end

"""
    step!(agent_step::Function, model::AbstractModel, nsteps::Integer, propagg::Dict, steps_to_collect_data::Array{Integer})

Repeats the `step` function `nsteps` times, and applies functions in values of the `propagg` dict to its keys at steps `steps_to_collect_data`.
"""
function step!(agent_step, model::AbstractModel, nsteps::Integer, propagg::Dict, steps_to_collect_data::Array{Int64})

  # Run the first step of the model to fill in the dataframe
  # step!(agent_step, model)
  df = data_collector(propagg, steps_to_collect_data, model, 1)

  for ss in 2:nsteps
    step!(agent_step, model)
    # collect data
    if ss in steps_to_collect_data
      df = data_collector(propagg, steps_to_collect_data, model, ss, df)
    end
  end
  # if 1 is not in `steps_to_collect_data`, remove the first row.
  if !in(1, steps_to_collect_data)
    df = df[2:end, :]
  end
  return df
end

"""
    step!(agent_step::Function, model_step::Function, model::AbstractModel)

Updates agents one step without collecting data. This function accepts two functions, one for update agents and one for updating the whole model one after all the agents have been updated.
"""
function step!(agent_step, model_step, model::AbstractModel)
  activation_order = return_activation_order(model)
  for index in activation_order
    agent_step(model.agents[index], model)
  end
  model_step(model)
end

"""
    step!(agent_step::Function, model_step::Function, model::AbstractModel, nsteps::Integer)

Repeats the `step` function `nsteps` times without collecting data.
"""
function step!(agent_step, model_step, model::AbstractModel, nsteps::Integer)
  for ss in 1:nsteps
    step!(agent_step, model_step, model)
  end
end

"""
    step!(agent_step::Function, model_step::Function, model::AbstractModel, nsteps::Integer, agent_properties::Array{Symbol}, steps_to_collect_data::Array{Integer})

Repeats the `step` function `nsteps` times, and collects all agent fields in `agent_properties` at steps `steps_to_collect_data`.
"""
function step!(agent_step, model_step, model::AbstractModel, nsteps::Integer, agent_properties::Array{Symbol}, steps_to_collect_data::Array{Int64})

  # Run the first step of the model to fill in the dataframe
  # step!(agent_step, model_step, model)
  df = data_collector(agent_properties, steps_to_collect_data, model, 1)

  for ss in 2:nsteps
    step!(agent_step, model_step, model)
    # collect data
    if ss in steps_to_collect_data
      df = data_collector(agent_properties, steps_to_collect_data, model, ss, df)
    end
  end
  # if 1 is not in `steps_to_collect_data`, remove the first columns. TODO: remove ids that were only present in the first step
  if !in(1, steps_to_collect_data)
    first_col = length(agent_properties)+2 # 1 for id and 1 for passing these agent properties
    end_col = size(df, 2)
    df = df[:, vcat([1], collect(first_col:end_col))]
  end
  return df
end

"""
    step!(agent_step::Function, model_step::Function, model::AbstractModel, nsteps::Integer, agent_properties::Array{Symbol}, aggregators::Array, steps_to_collect_data::Array{Integer})

Repeats the `step` function `nsteps` times, and applies functions in `aggregators` to values of agent fields in `agent_properties` at steps `steps_to_collect_data`.
"""
function step!(agent_step, model_step, model::AbstractModel, nsteps::Integer, agent_properties::Array{Symbol}, aggregators::Array, steps_to_collect_data::Array{Int64})

  # Run the first step of the model to fill in the dataframe
  # step!(agent_step, model_step, model)
  df = data_collector(agent_properties, aggregators, steps_to_collect_data, model, 1)

  for ss in 2:nsteps
    step!(agent_step, model_step, model)
    # collect data
    if ss in steps_to_collect_data
      df = data_collector(agent_properties, aggregators, steps_to_collect_data, model, ss, df)
    end
  end
  # if 1 is not in `steps_to_collect_data`, remove the first row.
  if !in(1, steps_to_collect_data)
    df = df[2:end, :]
  end
  return df
end

"""
    step!(agent_step::Function, model_step::Function, model::AbstractModel, nsteps::Integer, propagg::Dict, steps_to_collect_data::Array{Integer})

Repeats the `step` function `nsteps` times, and applies functions in values of the `propagg` dict to its keys at steps `steps_to_collect_data`.
"""
function step!(agent_step, model_step, model::AbstractModel, nsteps::Integer, propagg::Dict, steps_to_collect_data::Array{Int64})

  # Run the first step of the model to fill in the dataframe
  # step!(agent_step, model_step, model)
  df = data_collector(propagg, steps_to_collect_data, model, 1)

  for ss in 2:nsteps
    step!(agent_step, model_step, model)
    # collect data
    if ss in steps_to_collect_data
      df = data_collector(propagg, steps_to_collect_data, model, ss, df)
    end
  end
  # if 1 is not in `steps_to_collect_data`, remove the first row.
  if !in(1, steps_to_collect_data)
    df = df[2:end, :]
  end
  return df
end


text:
######### Stan batch program example  ###########

using Stan, Mamba

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

const eightschools ="
data {
  int<lower=0> J; // number of schools
  real y[J]; // estimated treatment effects
  real<lower=0> sigma[J]; // s.e. of effect estimates
}
parameters {
  real mu;
  real<lower=0> tau;
  real eta[J];
}
transformed parameters {
  real theta[J];
  for (j in 1:J)
    theta[j] <- mu + tau * eta[j];
}
model {
  eta ~ normal(0, 1);
  y ~ normal(theta, sigma);
}
"

const schools8data = [
  Dict("J" => 8,
    "y" => [28,  8, -3,  7, -1,  1, 18, 12],
    "sigma" => [15, 10, 16, 11,  9, 11, 10, 18],
    "tau" => 25
  )
]

stanmodel = Stanmodel(name="schools8", model=eightschools);
sim1 = stan(stanmodel, schools8data, ProjDir, CmdStanDir=CMDSTAN_HOME)

nodesubset = ["lp__", "accept_stat__", "mu", "tau", "theta.1", "theta.2", "theta.3", "theta.4", "theta.5", "theta.6", "theta.7", "theta.8"]

## Subset Sampler Output
sim = sim1[1:size(sim1, 1), nodesubset, 1:size(sim1, 3)]

describe(sim)
println()

## Brooks, Gelman and Rubin Convergence Diagnostic
try
  gelmandiag(sim1, mpsrf=true, transform=true) |> display
catch e
  #println(e)
  gelmandiag(sim, mpsrf=false, transform=true) |> display
end

## Geweke Convergence Diagnostic
gewekediag(sim) |> display

## Highest Posterior Density Intervals
hpd(sim) |> display

## Cross-Correlations
cor(sim) |> display

## Lag-Autocorrelations
autocor(sim) |> display

## Deviance Information Criterion
#dic(sim) |> display

## Plotting

p = plot(sim, [:trace, :mean, :density, :autocor], legend=true);
draw(p, nrow=4, ncol=4, filename="$(stanmodel.name)-summaryplot", fmt=:svg)
draw(p, nrow=4, ncol=4, filename="$(stanmodel.name)-summaryplot", fmt=:pdf)

# Below will only work on OSX, please adjust for your environment.
# JULIA_SVG_BROWSER is set from environment variable JULIA_SVG_BROWSER
@static is_apple() ? if isdefined(Main, :JULIA_SVG_BROWSER) && length(JULIA_SVG_BROWSER) > 0
        for i in 1:4
          isfile("$(stanmodel.name)-summaryplot-$(i).svg") &&
            run(`open -a $(JULIA_SVG_BROWSER) "$(stanmodel.name)-summaryplot-$(i).svg"`)
        end
      end : println()

end # cd