text:
#=
Functions for generating a single replication of a particular experiment
Dependencies: util.jl, sgd.jl, reformulation.jl, random_forests_po.jl, validation_set.jl
=#

using Parameters, DataFrames, Statistics, Random, Gurobi

## Type for strong the results of an experiment replication
@with_kw mutable struct replication_results
    SPOplus_spoloss_test::Union{Float64, Missing} = missing # SPO loss of SPO+ SGD on test set
    LS_spoloss_test::Union{Float64, Missing} = missing # SPO loss of LS SGD on test set
    SSVM_spoloss_test::Union{Float64, Missing} = missing # SPO loss of SSVM-Hamming SGD on test set
    RF_spoloss_test::Union{Float64, Missing} = missing # SPO loss of random forests on test set
    Absolute_spoloss_test::Union{Float64, Missing} = missing # SPO loss of ell_1 on test set
    Huber_spoloss_test::Union{Float64, Missing} = missing  # SPO loss of Huber on test set
    Baseline_spoloss_test::Union{Float64, Missing} = missing # SPO loss of baseline that returns mean of training data, on test set
    SPOplus_hammingloss_test::Union{Float64, Missing} = missing # Hamming loss of SPO+ SGD on test set
    LS_hammingloss_test::Union{Float64, Missing} = missing # Hamming loss of LS SGD on test set
    SSVM_hammingloss_test::Union{Float64, Missing} = missing # Hamming loss of SSVM SGD on test set
    RF_hammingloss_test::Union{Float64, Missing} = missing # Hamming loss of RF on test set
    Absolute_hammingloss_test::Union{Float64, Missing} = missing # Hamming loss of ell_1 on test set
    Huber_hammingloss_test::Union{Float64, Missing} = missing  # Hamming loss of Huber on test set
    Baseline_hammingloss_test::Union{Float64, Missing} = missing # Hamming loss of baseline on test set
    zstar_avg_test::Union{Float64, Missing} = missing # Average of z^*(c_i) on testing set
    r_squared::Union{Float64, Missing} = missing #r2 on training set from least squares
    r_squared_adj::Union{Float64, Missing} = missing #adj r2 on training set from least squares
end

"""
Make a blank data frame for storing the results of the experiment.
"""
function make_blank_complete_df()
    results_df = DataFrame(
        grid_dim = Int64[],
        n_train = Int64[],
        n_holdout = Int64[],
        n_test = Int64[],
        p_features = Int64[],
        polykernel_degree = Int64[],
        polykernel_noise_half_width = Float64[],
        SPOplus_spoloss_test = Union{Float64, Missing}[], # SPO loss of SPO+ SGD on test set
        LS_spoloss_test = Union{Float64, Missing}[],  # SPO loss of LS SGD on test set
        SSVM_spoloss_test = Union{Float64, Missing}[], # SPO loss of SSVM-Hamming SGD on test set
        RF_spoloss_test = Union{Float64, Missing}[], # SPO loss of random forests on test set
        Absolute_spoloss_test = Union{Float64, Missing}[], # SPO loss of ell_1 on test set
        Huber_spoloss_test = Union{Float64, Missing}[],  # SPO loss of Huber on test set
        Baseline_spoloss_test = Union{Float64, Missing}[], # SPO loss of baseline that returns mean of training data, on test set
        SPOplus_hammingloss_test = Union{Float64, Missing}[], # Hamming loss of SPO+ SGD on test set
        LS_hammingloss_test = Union{Float64, Missing}[], # Hamming loss of LS SGD on test set
        SSVM_hammingloss_test = Union{Float64, Missing}[], # Hamming loss of SSVM SGD on test set
        RF_hammingloss_test = Union{Float64, Missing}[], # Hamming loss of RF on test set
        Absolute_hammingloss_test = Union{Float64, Missing}[], # Hamming loss of ell_1 on test set
        Huber_hammingloss_test = Union{Float64, Missing}[],  # Hamming loss of Huber on test set
        Baseline_hammingloss_test = Union{Float64, Missing}[], # Hamming loss of baseline on test set
        zstar_avg_test = Union{Float64, Missing}[],  # Average of z^*(c_i) on testing set
        r_squared = Union{Float64, Missing}[], #r2 on training set from least squares
        r_squared_adj = Union{Float64, Missing}[]) #adj r2 on training set from least squares

    return results_df
end

"""
Make a new data frame that contains one row of results for an experiment instance that was just run.
"""
function build_complete_row(grid_dim,
    n_train, n_holdout, n_test,
    p_features, polykernel_degree, polykernel_noise_half_width, results_struct)

    df = make_blank_complete_df()

    push!(df, [
        grid_dim,
        n_train,
        n_holdout,
        n_test,
        p_features,
        polykernel_degree,
        polykernel_noise_half_width,
        results_struct.SPOplus_spoloss_test,
        results_struct.LS_spoloss_test,
        results_struct.SSVM_spoloss_test,
        results_struct.RF_spoloss_test,
        results_struct.Absolute_spoloss_test,
        results_struct.Huber_spoloss_test,
        results_struct.Baseline_spoloss_test,
        results_struct.SPOplus_hammingloss_test,
        results_struct.LS_hammingloss_test,
        results_struct.SSVM_hammingloss_test,
        results_struct.RF_hammingloss_test,
        results_struct.Absolute_hammingloss_test,
        results_struct.Huber_hammingloss_test,
        results_struct.Baseline_hammingloss_test,
        results_struct.zstar_avg_test,
        results_struct.r_squared,
        results_struct.r_squared_adj])

    return df
end

"""
Run one instance of the shortest path experiment. Uses the reformulation approach.
`num_lambda` is the number of lambdas used on the grid for each method with regularization.
Returns a replication_results struct.
"""
function shortest_path_replication(grid_dim,
    n_train, n_holdout, n_test,
    p_features, polykernel_degree, polykernel_noise_half_width;
    num_lambda = 10, lambda_max = missing, lambda_min_ratio = 0.0001, regularization = :ridge,
    gurobiEnvOracle = Gurobi.Env(), gurobiEnvReform = Gurobi.Env(), different_validation_losses = false,
    include_rf = true)

    # Get oracle
    sources, destinations = convert_grid_to_list(grid_dim, grid_dim)
    d_feasibleregion = length(sources)
    sp_oracle = sp_flow_jump_setup(sources, destinations, 1, grid_dim^2; gurobiEnv = gurobiEnvOracle)
    sp_graph = shortest_path_graph(sources = sources, destinations = destinations,
        start_node = 1, end_node = grid_dim^2, acyclic = true)

    # Generate Data
    B_true = rand(Bernoulli(0.5), d_feasibleregion, p_features)

    (X_train, c_train) = generate_poly_kernel_data_simple(B_true, n_train, polykernel_degree, polykernel_noise_half_width)
    (X_validation, c_validation) = generate_poly_kernel_data_simple(B_true, n_holdout, polykernel_degree, polykernel_noise_half_width)
    (X_test, c_test) = generate_poly_kernel_data_simple(B_true, n_test, polykernel_degree, polykernel_noise_half_width)

    # Add intercept in the first row of X
    X_train = vcat(ones(1,n_train), X_train)
    X_validation = vcat(ones(1,n_holdout), X_validation)
    X_test = vcat(ones(1,n_test), X_test)

    # Get Hamming labels
    (z_train, w_train) = oracle_dataset(c_train, sp_oracle)
    (z_validation, w_validation) = oracle_dataset(c_validation, sp_oracle)
    (z_test, w_test) = oracle_dataset(c_test, sp_oracle)

    c_ham_train = ones(d_feasibleregion, n_train) - w_train
    c_ham_validation = ones(d_feasibleregion, n_holdout) - w_validation
    c_ham_test = ones(d_feasibleregion, n_test) - w_test

    # Put train + validation together
    X_both = hcat(X_train, X_validation)
    c_both = hcat(c_train, c_validation)
    c_ham_both = hcat(c_ham_train, c_ham_validation)
    train_ind = collect(1:n_train)
    validation_ind = collect((n_train + 1):(n_train + n_holdout))

    # Set validation losses
    if different_validation_losses
        spo_plus_val_loss = :spo_loss
        ls_val_loss = :least_squares_loss
        ssvm_val_loss = :hamming_loss
        absolute_val_loss = :absolute_loss
        huber_val_loss = :huber_loss
    else
        spo_plus_val_loss = :spo_loss
        ls_val_loss = :spo_loss
        ssvm_val_loss = :spo_loss
        absolute_val_loss = :spo_loss
        huber_val_loss = :spo_loss
    end

    ### Algorithms ###

    # SPO+
    best_B_SPOplus, best_lambda_SPOplus = validation_set_alg(X_both, c_both, sp_oracle; sp_graph = sp_graph,
        train_ind = train_ind, validation_ind = validation_ind,
        val_alg_parms = val_parms(algorithm_type = :sp_spo_plus_reform, validation_loss = spo_plus_val_loss),
        path_alg_parms = reformulation_path_parms(num_lambda = num_lambda, lambda_max = lambda_max, regularization = regularization,
            gurobiEnv = gurobiEnvReform, lambda_min_ratio = lambda_min_ratio, algorithm_type = :SPO_plus))

    # Least squares
    best_B_leastSquares, best_lambda_leastSquares = validation_set_alg(X_both, c_both, sp_oracle; sp_graph = sp_graph,
        train_ind = train_ind, validation_ind = validation_ind,
        val_alg_parms = val_parms(algorithm_type = :ls_jump, validation_loss = ls_val_loss),
        path_alg_parms = reformulation_path_parms(num_lambda = num_lambda, lambda_max = lambda_max, regularization = regularization,
            gurobiEnv = gurobiEnvReform, lambda_min_ratio = lambda_min_ratio, algorithm_type = :leastSquares))

    # SSVM Hamming
    best_B_SSVM, best_lambda_SSVM = validation_set_alg(X_both, c_ham_both, sp_oracle; sp_graph = sp_graph,
        train_ind = train_ind, validation_ind = validation_ind,
        val_alg_parms = val_parms(algorithm_type = :sp_spo_plus_reform, validation_loss = ssvm_val_loss),
        path_alg_parms = reformulation_path_parms(num_lambda = num_lambda, lambda_max = lambda_max, regularization = regularization,
            gurobiEnv = gurobiEnvReform, lambda_min_ratio = lambda_min_ratio, algorithm_type = :SSVM_hamming))

    # RF
    if include_rf
        rf_mods = train_random_forests_po(X_train, c_train;
            rf_alg_parms = rf_parms())
    end

    # Absolute
    best_B_Absolute, best_lambda_Absolute = validation_set_alg(X_both, c_both, sp_oracle; sp_graph = sp_graph,
        train_ind = train_ind, validation_ind = validation_ind,
        val_alg_parms = val_parms(algorithm_type = :ls_jump, validation_loss = absolute_val_loss),
        path_alg_parms = reformulation_path_parms(num_lambda = num_lambda, lambda_max = lambda_max, regularization = regularization,
            po_loss_function = :absolute, gurobiEnv = gurobiEnvReform, lambda_min_ratio = lambda_min_ratio, algorithm_type = :Absolute))

    # Huber
    best_B_fake, best_lambda_fake = validation_set_alg(X_both, c_both, sp_oracle; sp_graph = sp_graph,
        train_ind = train_ind, validation_ind = validation_ind,
        val_alg_parms = val_parms(algorithm_type = :ls_jump, validation_loss = ls_val_loss),
        path_alg_parms = reformulation_path_parms(lambda_max = 0.0001, regularization = regularization,
            num_lambda = 2, gurobiEnv = gurobiEnvReform, algorithm_type = :Huber_LS_fake))

    fake_ls_list = abs.(vec(c_train - best_B_fake*X_train))
    delta_from_fake_ls = median(fake_ls_list)

    best_B_Huber, best_lambda_Huber = validation_set_alg(X_both, c_both, sp_oracle; sp_graph = sp_graph,
        train_ind = train_ind, validation_ind = validation_ind,
        val_alg_parms = val_parms(algorithm_type = :ls_jump, validation_loss = huber_val_loss),
        path_alg_parms = reformulation_path_parms(num_lambda = num_lambda, lambda_max = lambda_max, regularization = regularization, po_loss_function = :huber,
            huber_delta = delta_from_fake_ls, lambda_min_ratio = lambda_min_ratio, gurobiEnv = gurobiEnvReform, algorithm_type = :Huber))

    # Baseline
    c_bar_train = mean(c_train, dims=2)

    c_bar_vector = zeros(d_feasibleregion, n_train)
    for i = 1:n_train
        c_bar_vector[:, i] = c_bar_train
    end

    # Caclulate R_2 from least squares to measure fit

    (r_squared, r_squared_adj) = r_2(best_B_leastSquares, X_train, c_train, c_bar_vector, n_train, p_features)

    ### Populate final results ###
    final_results = replication_results()

    final_results.SPOplus_spoloss_test = spo_loss(best_B_SPOplus, X_test, c_test, sp_oracle)
    final_results.LS_spoloss_test = spo_loss(best_B_leastSquares, X_test, c_test, sp_oracle)
    final_results.SSVM_spoloss_test = spo_loss(best_B_SSVM, X_test, c_test, sp_oracle)

    if include_rf
        rf_preds_test = predict_random_forests_po(rf_mods, X_test)
        final_results.RF_spoloss_test = spo_loss(Matrix(1.0I, d_feasibleregion, d_feasibleregion), rf_preds_test, c_test, sp_oracle)
    else
        final_results.RF_spoloss_test = missing
    end

    final_results.Absolute_spoloss_test = spo_loss(best_B_Absolute, X_test, c_test, sp_oracle)
    final_results.Huber_spoloss_test = spo_loss(best_B_Huber, X_test, c_test, sp_oracle)

    c_bar_test_preds = zeros(d_feasibleregion, n_test)
    for i = 1:n_test
        c_bar_test_preds[:, i] = c_bar_train
    end
    final_results.Baseline_spoloss_test = spo_loss(Matrix(1.0I, d_feasibleregion, d_feasibleregion), c_bar_test_preds, c_test, sp_oracle)

    # Now Hamming
    final_results.SPOplus_hammingloss_test = spo_loss(best_B_SPOplus, X_test, c_ham_test, sp_oracle)
    final_results.LS_hammingloss_test = spo_loss(best_B_leastSquares, X_test, c_ham_test, sp_oracle)
    final_results.SSVM_hammingloss_test = spo_loss(best_B_SSVM, X_test, c_ham_test, sp_oracle)

    if include_rf
        final_results.RF_hammingloss_test = spo_loss(Matrix(1.0I, d_feasibleregion, d_feasibleregion), rf_preds_test, c_ham_test, sp_oracle)
    else
        final_results.RF_hammingloss_test = missing
    end

    final_results.Absolute_hammingloss_test = spo_loss(best_B_Absolute, X_test, c_ham_test, sp_oracle)
    #final_results.Huber_hammingloss_test = (best_B_Huber, X_test, c_ham_test, sp_oracle)
    final_results.Baseline_hammingloss_test = spo_loss(Matrix(1.0I, d_feasibleregion, d_feasibleregion), c_bar_test_preds, c_ham_test, sp_oracle)

    final_results.zstar_avg_test = mean(z_test)

    final_results.r_squared = r_squared
    final_results.r_squared_adj = r_squared_adj

    return final_results
end

"""
Run multiple replications of the shortest path experiment and return a data frame.
`rng_seed` and `num_trials` are integers
`n_train_vec`, `polykernel_degree_vec`, and `polykernel_noise_half_width_vec` are
vectors so that we try all combinations of these values.
"""
function shortest_path_multiple_replications(rng_seed, num_trials, grid_dim,
    n_train_vec, n_test,
    p_features, polykernel_degree_vec, polykernel_noise_half_width_vec;
    num_lambda = 10, lambda_max = missing, lambda_min_ratio = 0.0001, holdout_percent = 0.25, regularization = :ridge,
    different_validation_losses = false,
    include_rf = true)

    Random.seed!(rng_seed)

    big_results_df = make_blank_complete_df()

    envOracle = setup_gurobi_env(method_type = :default, use_time_limit = false)
    envReform = setup_gurobi_env(method_type = :method3, use_time_limit = false)

    for n_train in n_train_vec
        for polykernel_degree in polykernel_degree_vec
            for polykernel_noise_half_width in polykernel_noise_half_width_vec
                println("Moving on to n_train = $n_train, polykernel_degree = $polykernel_degree, polykernel_noise_half_width = $polykernel_noise_half_width")
                for trial = 1:num_trials
                    println("Current trial is $trial")
                    n_holdout = round(Int, holdout_percent*n_train)

                    current_results = shortest_path_replication(grid_dim,
                        n_train, n_holdout, n_test,
                        p_features, polykernel_degree, polykernel_noise_half_width;
                        num_lambda = num_lambda, lambda_max = lambda_max, lambda_min_ratio = lambda_min_ratio, regularization = regularization,
                        gurobiEnvOracle = envOracle, gurobiEnvReform = envReform,
                        different_validation_losses = different_validation_losses,
                        include_rf = include_rf)

                    current_df_row = build_complete_row(grid_dim,
                        n_train, n_holdout, n_test,
                        p_features, polykernel_degree, polykernel_noise_half_width, current_results)

                    big_results_df = vcat(big_results_df, current_df_row)
                end
            end
        end
    end

    return big_results_df
end

max_stars_repo_name: katerinagiannoutsou/dso-coding-project

max_issues_repo_name: katerinagiannoutsou/dso-coding-project

max_forks_repo_name: katerinagiannoutsou/dso-coding-project


text:
# 1 -> 2 <- 5
# ^    |    ^
# |    |    |
# |    v    |
# 3 <- 4 -> 6

station1 = EuclideanStation(1, EuclideanCoordinate(1, 4), "test", false)
station2 = EuclideanStation(2, EuclideanCoordinate(2, 4), "test", false)
station3 = EuclideanStation(3, EuclideanCoordinate(1, 1), "test", false)
station4 = EuclideanStation(4, EuclideanCoordinate(2, 2), "test", false)
station5 = EuclideanStation(5, EuclideanCoordinate(3, 4), "test", false)
station6 = EuclideanStation(6, EuclideanCoordinate(3, 1), "test", false)
line1 = Line("S1")
line2 = Line("S2")
edge = ProcessedEdge(station1, station2, line1, 0, 1, false)
edge3 = ProcessedEdge(station2, station4, line1, 6, 1, false)
edge4 = ProcessedEdge(station4, station3, line1, 4, 1, false)
edge5 = ProcessedEdge(station3, station1, line1, 2, 2, false)
edge6 = ProcessedEdge(station5, station2, line2, 4, 1, false)
edge7 = ProcessedEdge(station4, station6, line2, 0, 1, false)
edge8 = ProcessedEdge(station6, station5, line2, 2, 1, false)
transit_map = InputGraph([station1, station2, station3, station4, station5, station6],
    [edge, edge3, edge4, edge5, edge6, edge7, edge8],
    [line1, line2])

result = optimize(GLPKSolverMIP(msg_lev = 3), transit_map, true)

@test length(result.faces) == 2
@test Set([edge, edge3, edge4, edge5]) in result.faces
@test Set([edge6, edge7, edge8, edge3]) in result.faces

ext: jl

max_stars_repo_name: brainfo/MetroMap

max_issues_repo_name: brainfo/MetroMap

max_forks_repo_name: brainfo/MetroMap


text:
module Finometrics
using Revise

if pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

##################Dependencies

using  DataFrames, Distributions, StatsBase, GLM, CategoricalArrays,
  Dates, NLopt, ForwardDiff, Formatting, TexTables, DataStructures,
  LinearAlgebra
#=NOTE: Replace with LinearAlgebra=#

######################Methods####################
export FMLM, #Regression methods
  FM2SLS,
  project!,
  getResid!,
  pullModelMatrix,
  getCoeff!,
  getHomoskedΣ!,
  getModWhiteΣ!,
  getWhiteΣ!,
  getNeweyWest!,
  getModWhiteΣSlow,
  getHomoskedΣSlow,
  getWhiteΣSlow,
  getNeweyWestSlow,
  dropNullsFromDF,
  FMExpr,
  FMSym,
  FMQR,
  get1stStage,
  getTerm,
  getR,
  getR²,
  getNeweyWestFunc,
  getClustered!,

  texTable, #IO Mthods
  writeTables2File,
  writeNakedTable,
  array2String,
  num2Str,
  vec2String,
  vbind,

  skewnessStat, #FMstat methods
  kurtosisStat,
  JBStat,
  UnivariateTest,
  autoρ,
  ljungBox,
  test,
  pAutoρ,
  ARMA,
  ARMAFillYHat!,
  ARMALogLike,
  ARMAEstimate,
  generalizeARMA!,
  narrowARMA!,
  lagDF!,

  FMΦInv, #FMNumerical
  FMΦ,

  VAR, #FMVAR
  lagDF!,
  VAREstimate,
  VARPrediction,
  varΣ,
  propagateImpulse,

  wrdsQuery,  #FMWRDS, doesn't work
  wrdsGet,
  wrdsDrop,

  NDict,#Exported types
  NString,
  NSymbol,
  NType,

  MString,
  MBool,
  MFloat64,
  MInt,
  MSymbol

##################Custom types
###NOTE: 0.7 Compat Hacks, delte on upgrade
#const Nothing = Void
#const LinearAlgebra = Base.LinAlg #NOTE:Compat Hack
#const AbstractRange = Range #NOTE compat hack

###FMReg Types
const FMExpr = Union{Symbol,Expr,Nothing}
const FMData = Union{Float64,Int,Date,Symbol, CategoricalValue, Missing, Nothing}#,
abstract type FMModel end

###Convenience Types for Export
const NSymbol = Union{Nothing, Symbol}
const NType = Union{Nothing, Type}
const NString = Union{Nothing, String}

const MBool = Union{Bool, Missing}
const MInt = Union{Int, Missing}
const MFloat64 = Union{Float64, Missing}
const MSymbol = Union{Symbol, Missing}
const MString = Union{String, Missing}

###VAR Types
const Structure = Union{Nothing, Vector{Float64}}
const CoefficientErrors = Union{Nothing, Vector{Float64}}
const RegressionSS = Union{Nothing, Float64}
const RegressionN = Union{Nothing, Int}

################Constants

#shared cosntants
const DEFAULT_STAR_LEGEND = "*p<0.1, **p<0.05, ***p<0.01"
const DEBUG_FMMOD = false
const DEFAULT_DECIMALS = 4

#Stat constants
const MAX_TIME = 30.
const F_TOL_ABS = 10. ^ -12.
const F_TOL_REL = 10. ^ -10.
const OPT_BOUND = 2.

#numerical constants
const sqrt2Inv = 1.0/(2.0^.5)
const sqrt2PiInv = 1.0/(2.0*π)^0.5
#  MFloat64, MInt, MBool, MSymbol}

const UPPER_PAD = "%This table was programatically generated using the TexTables package and my own code.\n"
const LOWER_PAD = "\n"

###################Files############

include("FinometricsCode\\FMReg.jl")
include("FinometricsCode\\FMIO.jl")
include("FinometricsCode\\FMStat.jl")
include("FinometricsCode\\FMVAR.jl")
include("FinometricsCode\\FMNumerical.jl")
include("FinometricsCode\\FMWRDS.jl")

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
## -------------------------------------------------------------------
const _UPDATE_RECORD = Dict()
const _WATCHING_TASKS = Dict()

## -------------------------------------------------------------------
mutable struct OnDiskMonitor
    last_save_time::Real
    dt::Real
    file::String
    cache::Dict # Only for ram data
    lk::ReentrantLock

    function OnDiskMonitor(p, ps...; dt = 10, lock = ReentrantLock())
        file = joinpath(p, ps...)
        # clear && rm(file; force = true)
        last_save_time = -1.0 # ensure to save the first time
        cache = Dict()
        new(last_save_time, dt, file, cache, lock)
    end
end

_worw_lock(f::Function, lk, dolk::Bool) = dolk ? lock(f, lk) : f()

function sync_from_disk!(m::OnDiskMonitor; dolk = true)
    !isfile(m.file) && return m
    _worw_lock(m.lk, dolk) do
        m.cache = deserialize(m.file)
    end
    m
end

"""
`OnDiskMonitor(file::String; dt = 10, lock = ReentrantLock())`

Example
```
n = 100
m = OnDiskMonitor("path/to/file")
# Watch task (It can be in another process)
watch(m) do dat
    xs = dat[:xs]
    ys = dat[:ys]
    p = scatter(xs, yx;
        xlabel = "x", ylabel = "y", label = "",
        color = rand([:red, :green, :blue])
    )
    figfile = joinpath(dirname(m.file), "fig1.png")
    sfig(p, figfile)
end

# Fake job to monitor (All is tread save by default)
@threads for i in 1:n
    x = 10*pi*i/n
    y = sin(x)^2
    record!(m) do dat
        push!(get!(dat, :xs, []), x)
        push!(get!(dat, :ys, []), y)
    end
    sleep(10.0/n)
end
sync_to_disk(m)
```
"""
OnDiskMonitor

function sync_to_disk(m::OnDiskMonitor; dolk = true)
    _worw_lock(m.lk, dolk) do
        serialize(m.file, m.cache)
    end
    m
end

function clear_file(m::OnDiskMonitor; dolk = true)
    !isfile(m.file) && return m
    _worw_lock(m.lk, dolk) do
        rm(m.file; force = true)
    end
    m
end

"""
`watch(onupdate::Function, m::OnDiskMonitor; dolk = true, wt = m.dt, iters = typemax(Int))`.

Will start a task (using `@async`) that checks every `wt` secunds if the
monitored file changed.
If so, will load the data and execute `onupdate` passing that data as only argument.
`dolk` controls wheather to use the monitor lock or not.
`iters` controls how many check iters to do (mostly for dev).

Example
```
m = OnDiskMonitor("path/to/file")
watch(m) do dat
    xs = dat[:xs]
    ys = dat[:ys]
    p = plot(xs, yx;
        xlabel = "x", ylabel = "y", label = "",
        color = rand([:red, :green, :blue])
    )
    figfile = joinpath(dirname(m.file), "fig1.png")
    sfig(p, figfile)
end
```
"""
function watch(onupdate::Function, m::OnDiskMonitor;
        dolk = true, wt = m.dt, iters = typemax(Int),
        onerr = (err) -> @warn("Error on loading/update", err)
    )

    # register task
    funhash = hash(onupdate)
    _worw_lock(m.lk, dolk) do
        _WATCHING_TASKS[m.file] = funhash
    end

    for it in 1:iters
        sleep(wt)

        # Register task
        !isfile(m.file) && continue
        _worw_lock(m.lk, dolk) do
            lmtime = get!(_UPDATE_RECORD, m.file, -1.0)
            fmtime = mtime(m.file)
            if lmtime != fmtime
                try
                    dat = deserialize(m.file)
                    onupdate(dat)
                    catch err; onerr(err)
                end
            end
            _UPDATE_RECORD[m.file] = fmtime
        end

        reghash = get(_WATCHING_TASKS, m.file, nothing)
        reghash != funhash && break
    end
    return nothing
end

"""
`record!(add_to_cache!::Function, m::OnDiskMonitor; dolk = true)`

Execute the given function `add_to_cache!` passing the monitor internal
cache (a `Dict`) as its only argument.
After that, it uses `m.dt` to check is the data must be saved
to disk.
If so, it create/overwrite `m.file`.
`dolk` controls wheather to use the monitor lock or not.

Example
```
n = 100
m = OnDiskMonitor("path/to/file")
for i in 1:n
    x = 10*pi*i/n
    y = sin(x)^2
    record!(m) do dat
        push!(get!(dat, :xs, []), x)
        push!(get!(dat, :ys, []), y)
    end
    sleep(10.0/n)
end
```
"""
function record!(add_to_cache!::Function, m::OnDiskMonitor; dolk = true)
    _worw_lock(m.lk, dolk) do
        add_to_cache!(m.cache)

        curr_time = time()
        dosave = (curr_time - m.last_save_time) > m.dt
        if dosave
            m.last_save_time = curr_time
            sync_to_disk(m; dolk = false)
        end
    end
end

"""
`reset!(m; dolk = true)`

Reset (empty!) the internal monitor cache and also the data at `m.file`.
`dolk` controls wheather to use the monitor lock or not.
"""
function reset!(m::OnDiskMonitor; dolk = true)
    _worw_lock(m.lk, dolk) do
        empty!(m.cache)
        isfile(m.file) && sync_to_disk(m; dolk = false)
        finish_watch_task(m)
    end
    m
end

function finish_watch_task(m::OnDiskMonitor)
    _WATCHING_TASKS[m.file] = nothing
end

get_cache(m::OnDiskMonitor) = m.cache


text:
"""
`smooth(y, span=5)`

Smooth a vector using a moving average filter. Endpoints are handled by
collapsing the length of the filter as showed below.
```
yy[1] = y[1]
yy[2] = (y[1] + y[2] + y[3]) / 3
yy[3] = (y[1] + y[2] + y[3] + y[4] + y[5]) / 5
yy[4] = (y[2] + y[3] + y[4] + y[4] + y[6]) / 5
```

**Arguments**
- `y`: Array of values.

**Optional Arguments**
- `span`: Filter length.

**Returns**
- `yy`: Smoothed array of values.
"""
function smooth(y, span=5)

    halfspan = floor(span/2)
    n = length(y)
    yy = fill(NaN, n)

    for i = 1:n

        ispan = min(i-1, n-i, halfspan)
        ileft = convert(Int, i-ispan)
        iright = convert(Int, i+ispan)

        yy[i] = mean(y[ileft:iright])

    end

    return yy

end

max_stars_repo_name: vkumpost/stoosc

max_issues_repo_name: vkumpost/stoosc

max_forks_repo_name: vkumpost/stoosc


text:
## Neuron name
# name = "spindle/whole_neurons/03a_spindle2aFI"
# name = "spindle/whole_neurons/03b_spindle4aACC"
# name = "spindle/whole_neurons/03b_spindle7aACC"
# name = "spindle/whole_neurons/04b_spindle3aFI"
# name = "spindle/whole_neurons/06b_spindle8aACC"
# name = "spindle/whole_neurons/16o_spindle13aFI"
# name = "spindle/whole_neurons/16o_spindle13aFI.stl"
# name = "spindle/separated_neurons/03a_spindle2aFI_soma"
# name = "spindle/separated_neurons/03a_spindle2aFI_dendrites_1"
# name = "spindle/separated_neurons/03a_spindle2aFI_dendrites_2"
# name = "spindle/separated_neurons/03b_spindle4aACC_dendrites_1"
# name = "spindle/separated_neurons/03b_spindle4aACC_dendrites_2"
# name = "spindle/separated_neurons/03b_spindle4aACC_soma"
# name = "spindle/separated_neurons/03b_spindle6aACC_dendrites_1"
# name = "spindle/separated_neurons/03b_spindle6aACC_dendrites_2"
name = "spindle/separated_neurons/03b_spindle6aACC_soma"
# name = "spindle/separated_neurons/03b_spindle7aACC_dendrites_1"
# name = "spindle/separated_neurons/03b_spindle7aACC_dendrites_2"
# name = "spindle/separated_neurons/03b_spindle7aACC_soma"
# name = "spindle/separated_neurons/04b_spindle3aFI_dendrites_1"
# name = "spindle/separated_neurons/04b_spindle3aFI_dendrites_2"
# name = "spindle/separated_neurons/04b_spindle3aFI_soma"
# name = "spindle/separated_neurons/19o_spindle14aFI_dendrites_2.stl"
# name = "pyramidal/separated_neurons/02b_pyramidal1aACC_dendrites"
# name = "pyramidal/separated_neurons/02b_pyramidal1aACC_soma"
# name = "pyramidal/whole_neurons/02b_pyramidal1aACC"
# name = "pyramidal/whole_neurons/02a_pyramidal2aFI"
# name = "pyramidal/whole_neurons/04a_pyramidal4aACC"
# name = "pyramidal/whole_neurons/04a_pyramidal5aACC"
# name = "pyramidal/whole_neurons/04b_pyramidal6aACC"
# name = "pyramidal/whole_neurons/04b_pyramidal6aFI"
# name = "pyramidal/whole_neurons/25o_pyramidal18aFI"

# Floating point type for simulations
T = Float64

# Geometrical setup
setup = NeuronSetup{T}(;
    name,
    ecs_shape = :no_ecs,
    ecs_ratio = 0.3,
    # refinement = 0.5,
)

# Get compartimentalized coefficient vectors
coeffs = coefficients(
    setup;
    D = (; neuron = 0.002 * I(3), ecs = 0.002 * I(3)),
    T₂ = (; neuron = Inf, ecs = Inf),
    ρ = (; neuron = 1.0, ecs = 1.0),
    κ = (; neuron_ecs = 1e-4, in = 0.0, neuron = 0.0, ecs = 0.0),
    γ = 2.67513e-4,
)


text:
# This file is part of the TaylorIntegration.jl package; MIT licensed

# Load necessary components of Espresso
using Espresso: subs, simplify, ExGraph, ExH, to_expr, sanitize, genname,
    find_vars, findex, find_indices, isindexed

"""
    BookKeeping

Mutable struct that contains all the bookkeeping vectors/dictionaries used within
`_make_parsed_jetcoeffs`:
    - `d_indx`     : Dictionary mapping new variables (symbols) to old (perhaps indexed) symbols
    - `d_assign`   : Dictionary with the numeric assignments (that are substituted)
    - `d_decl`     : Dictionary declared arrays
    - `v_newvars`  : Symbols of auxiliary indexed vars
    - `v_arraydecl`: Symbols which are explicitly declared as Array or Vector
    - `v_preamb`   : Symbols or Expr used in the preamble (declarations, etc)
    - `retvar`     : *Guessed* returned variable, which defines the LHS of the ODEs

"""
mutable struct BookKeeping
    d_indx   :: Dict{Symbol, Expr}
    d_assign :: Dict{Union{Symbol,Expr}, Number}
    d_decl   :: Dict{Symbol, Expr}
    v_newvars   :: Vector{Symbol}
    v_arraydecl :: Vector{Symbol}
    v_preamb :: Vector{Union{Symbol,Expr}}
    retvar   :: Symbol

    function BookKeeping()
        return new(Dict{Symbol, Expr}(), Dict{Union{Symbol,Expr}, Number}(),
            Dict{Symbol, Expr}(), Symbol[], Symbol[], Union{Symbol,Expr}[], :nothing)
    end
end

"""
`inbookkeeping(v, bkkeep::BookKeeping)`

Checks if `v` is declared in `bkkeep`, considering the `d_indx`, `v_newvars` and
`v_arraydecl` fields.

"""
@inline inbookkeeping(v, bkkeep::BookKeeping) = (v ∈ keys(bkkeep.d_indx) ||
    v ∈ bkkeep.v_newvars || v ∈ bkkeep.v_arraydecl)

# Constants to create the structure of the new jetcoeffs! and _allocate_jetcoeffs! methods.
# The (irrelevant) `nothing` below is there to have a `:block` Expr; it is deleted later
const _HEAD_PARSEDFN_SCALAR = sanitize(:(
    function TaylorIntegration.jetcoeffs!(::Val{__fn}, __tT::Taylor1{_T}, __x::Taylor1{_S},
            __params, __tmpTaylor, __arrTaylor) where {_T<:Real, _S<:Number}
        order = __tT.order
        nothing
    end)
);

const _HEAD_PARSEDFN_VECTOR = sanitize(:(
    function TaylorIntegration.jetcoeffs!(::Val{__fn}, __tT::Taylor1{_T},
            __x::AbstractArray{Taylor1{_S}, _N}, __dx::AbstractArray{Taylor1{_S}, _N},
            __params, __tmpTaylor, __arrTaylor) where {_T<:Real, _S<:Number, _N}
        order = __tT.order
        nothing
    end)
);

const _HEAD_ALLOC_TAYLOR1_SCALAR = sanitize(:(
    function TaylorIntegration._allocate_jetcoeffs!( ::Val{__fn}, __tT::Taylor1{_T},
            __x::Taylor1{_S}, __params) where {_T<:Real, _S<:Number}
        order = __tT.order
        nothing
    end)
);

const _HEAD_ALLOC_TAYLOR1_VECTOR = sanitize(:(
    function TaylorIntegration._allocate_jetcoeffs!( ::Val{__fn}, __tT::Taylor1{_T},
            __x::AbstractArray{Taylor1{_S}, _N}, __dx::AbstractArray{Taylor1{_S}, _N},
            __params) where {_T<:Real, _S<:Number, _N}
        order = __tT.order
        nothing
    end)
);

# Constants for the initial declaration and initialization of arrays
const _DECL_ARRAY = sanitize( Expr(:block,
    :(__var1 = Array{Taylor1{_S}}(undef, __var2)),
    :(__var1 .= Taylor1( zero(_S), order )))
);

"""
`_make_parsed_jetcoeffs( ex, debug=false )`

This function constructs the expressions of two new methods, the first equivalent to the
differential equations (jetcoeffs!), which exploits the mutating functions of TaylorSeries.jl,
and the second one (_allocate_jetcoeffs) preallocates any auxiliary `Taylor1` or
`Vector{Taylor1{T}}` needed.

"""
function _make_parsed_jetcoeffs(ex::Expr, debug=false)
    # Extract the name, args and body of the function
    # `fn` name of the function having the ODEs
    # `fnargs` arguments of the function
    # `fnbody` future transformed function body
    fn, fnargs, fnbody = _extract_parts(ex)
    debug && (println("****** _extract_parts ******");
        @show(fn, fnargs, fnbody); println())

    # Set up the Expr for the new functions
    new_jetcoeffs, new_allocjetcoeffs = _newhead(fn, fnargs)

    # Expr for the for-loop block for the recursion (of the `x` variable)
    forloopblock = Expr(:for, :(ord = 1:order-1), Expr(:block, :(ordnext = ord + 1)) )

    # Transform the graph representation of the body of the functions:
    # defspreamble: inicializations used for the zeroth order (preamble)
    # defsprealloc: definitions (declarations) of auxiliary Taylor1's
    # fnbody: transformed function body, using mutating functions from TaylorSeries;
    #         used later within the recursion loop
    # bkkeep: book-keeping structure having info of the variables
    defspreamble, defsprealloc, fnbody, bkkeep = _preamble_body(fnbody, fnargs, debug)
    debug && (println("****** _preamble_body ******");
        @show(defspreamble); println(); @show(defsprealloc); println();
        @show(fnbody); println(); @show(bkkeep); println())

    # Create body of recursion loop; temporary assignements may be needed.
    # rec_preamb: recursion loop for the preamble (first order correction)
    # rec_fnbody: recursion loop for the body-function (recursion loop for higher orders)
    rec_preamb, rec_fnbody = _recursionloop(fnargs, bkkeep)
    debug && (println("****** _recursionloop ******");
        @show(rec_preamb); println(); @show(rec_fnbody); println())

    # Add rec_fnbody to forloopblock
    push!(forloopblock.args[2].args, fnbody.args[1].args..., rec_fnbody)

    # Add preamble and recursion body to `new_allocjetcoeffs`
    push!(new_jetcoeffs.args[2].args, defspreamble..., rec_preamb)

    # Push preamble and forloopblock to `new_allocjetcoeffs` and return line
    push!(new_jetcoeffs.args[2].args, forloopblock, Meta.parse("return nothing"))

    # Define the expressions of the returned vectors in `new_allocjetcoeffs`
    push!(new_allocjetcoeffs.args[2].args, defsprealloc...)
    if isempty(bkkeep.v_newvars)
        ret_defsprealloc = :(Taylor1{_S}[])
    else
        ret_defsprealloc = :([$(bkkeep.v_newvars...),])
    end
    if isempty(bkkeep.v_arraydecl)
        ret_tmparrays = :([Vector{Taylor1{_S}}(undef, 0),])
    else
        ret_tmparrays = :([$(bkkeep.v_arraydecl...),])
    end
    ret_ret = :(return $(ret_defsprealloc), $(ret_tmparrays))

    # Add return line to `new_allocjetcoeffs`
    push!(new_allocjetcoeffs.args[2].args, ret_ret)

    # Rename variables in the calling form of the new methods
    if length(fnargs) == 3
        new_jetcoeffs = subs(new_jetcoeffs, Dict(:__tT => fnargs[3],
            :__params => fnargs[2], :(__x) => fnargs[1], :(__fn) => fn ))
        new_allocjetcoeffs = subs(new_allocjetcoeffs,
            Dict(:__tT => fnargs[3], :__params => fnargs[2],
                :(__x) => fnargs[1], :(__fn) => fn ))

    elseif length(fnargs) == 4
        new_jetcoeffs = subs(new_jetcoeffs, Dict(:__tT => fnargs[4],
            :__params => fnargs[3], :(__x) => fnargs[2], :(__dx) => fnargs[1],
            :(__fn) => fn ))
        new_allocjetcoeffs = subs(new_allocjetcoeffs,
            Dict(:__tT => fnargs[4], :__params => fnargs[3],
                :(__x) => fnargs[2], :(__dx) => fnargs[1], :(__fn) => fn ))

    else
        # A priori this is not needed
        throw(ArgumentError("Wrong number of arguments in `fnargs`"))
    end

    # Bring back old variable-names to `new_jecoeffs`; includes proper renaming
    # of the temporal arrays
    dd = Dict{Symbol, Expr}()
    @inbounds for (ind, vnew) in enumerate(bkkeep.v_newvars)
        merge!(dd, Dict(vnew => Expr(:ref, :__tmpTaylor, ind)))
    end
    @inbounds for (ind, vnew) in enumerate(bkkeep.v_arraydecl)
        merge!(dd, Dict(vnew => Expr(:ref, :__arrTaylor, ind)))
    end
    new_jetcoeffs = subs(new_jetcoeffs, dd)

    return new_jetcoeffs, new_allocjetcoeffs
end

"""
`_extract_parts(ex::Expr)`

Returns the function name, the function arguments, and the body of
a function passed as an `Expr`. The function may be provided as
a one-line function, or in the long form (anonymous functions
do not work).

"""
function _extract_parts(ex::Expr)
    # Capture name, args and body
    # @capture( shortdef(ex), ffn_(ffnargs__) = fnbody_ ) ||
    #     throw(ArgumentError("It must be a function call:\n $ex"))
    dd = Dict{Symbol, Any}()
    _capture_fn_args_body!(ex, dd)
    ffn, ffnargs, fnbody = dd[:fname], dd[:fnargs], dd[:fnbody]

    # Clean-up `ffn`
    fn = isa(ffn, Symbol) ? ffn : ffn.args[1]

    # Clean-up (from declarations) `ffnargs`
    fnargs = []
    for ff in ffnargs
        if isa(ff, Symbol)
            push!(fnargs, ff)
        elseif isa(ff, Expr)
            push!(fnargs, ff.args[1])
        else
            throw(ArgumentError("$ff is neither a `Symbol` nor an `Expr`"))
        end
    end

    # Standarize fnbody, same structure for one-line or long-form functions
    if length(fnbody.args) > 1
        fnbody.args[1] = Expr(:block, copy(fnbody.args)...)
        deleteat!(fnbody.args, 2:length(fnbody.args))
    end

    # Special case: the last arg of the function is a simple
    # assignement (symbol), a numerical value or simply `nothing`
    if isa(fnbody.args[1].args[end], Symbol)
        # This accouns the case of having `nothing` at the end of the function
        if fnbody.args[1].args[end] == :nothing
            pop!(fnbody.args[1].args)
        else
            fnbody.args[1].args[end] = :( identity($(fnbody.args[1].args[end])) )
        end
    elseif isa(fnbody.args[1].args[end], Number)
        fnbody.args[1].args[end] = :( $(fnbody.args[1].args[end])+zero($(fnargs[1])) )
    end

    return fn, fnargs, fnbody
end

"""
`_capture_fn_args_body!(ex, dd::Dict{Symbol, Any})`

Captures the name of a function, arguments, body and other properties,
returning them as the values of the dictionary `dd`, which is updated
in place.

"""
function _capture_fn_args_body!(ex::Expr, dd::Dict{Symbol, Any} = Dict())
    exh = ExH(ex)
    if exh.head == :(=) || exh.head == :function
        _capture_fn_args_body!.(ex.args, Ref(dd))
    elseif exh.head == :call
        push!(dd, Pair(:fname, exh.args[1]))
        push!(dd, Pair(:fnargs, exh.args[2:end]))
    elseif exh.head == :block
        push!(dd, Pair(:fnbody, ex))
    elseif exh.head == :where
        push!(dd, Pair(:where, ex.args[2:end]))
        _capture_fn_args_body!(ex.args[1], dd)
    end

    return nothing
end

"""
`_newhead(fn, fnargs)`

Creates the head of the new method of `jetcoeffs!` and `_allocate_jetcoeffs`.
`fn` is the name of the passed function and `fnargs` is a vector with its arguments
defning the function (which are either three or four).

"""
function _newhead(fn, fnargs)
    # Construct common elements of the new expression
    if length(fnargs) == 3
        newfunction = copy(_HEAD_PARSEDFN_SCALAR)
        newdeclfunc = copy(_HEAD_ALLOC_TAYLOR1_SCALAR)

    elseif length(fnargs) == 4
        newfunction = copy(_HEAD_PARSEDFN_VECTOR)
        newdeclfunc = copy(_HEAD_ALLOC_TAYLOR1_VECTOR)

    else
        throw(ArgumentError(
        "Wrong number of arguments in the definition of the function $fn"))
    end

    # Delete the irrelevant `nothing`
    pop!(newfunction.args[2].args)
    pop!(newdeclfunc.args[2].args)

    return newfunction, newdeclfunc
end

"""
`_preamble_body(fnbody, fnargs, debug=false)`

Returns expressions for the preamble, the declaration of
arrays, the body and the bookkeeping struct, which will be used to build
the new functions. `fnbody` is the expression with the body of
the original function (already adapted), `fnargs` is a vector of symbols
of the original diferential equations function.

"""
function _preamble_body(fnbody, fnargs, debug=false)
    # Inicialize BookKeeping struct
    bkkeep = BookKeeping()

    # Rename vars to have the body in non-indexed form; bkkeep has different entries
    # for bookkeeping variables/symbolds, including indexed ones
    fnbody, bkkeep.d_indx = _rename_indexedvars(fnbody)
    debug && (println("------ _rename_indexedvars ------");
        @show(fnbody); println(); @show(bkkeep.d_indx); println())

    # Create `newfnbody` which corresponds to `fnbody`, cleaned (without irrelevant comments)
    # and with all new variables in place; bkkeep.d_indx is updated
    newfnbody = _newfnbody!(fnbody, fnargs, bkkeep)
    debug && (println("------ _newfnbody ------");
        @show(newfnbody); println(); @show(bkkeep); println())

    # Parse `newfnbody` and create `prepreamble` and `prealloc`, updating `bkkeep`.
    # These objects use the mutating functions from TaylorSeries.
    preamble = Expr(:block,)
    prealloc = Expr(:block,)
    _parse_newfnbody!(newfnbody, preamble, prealloc, bkkeep, false)
    # Get rid of the initial `:block`
    preamble = preamble.args[1]
    prealloc = prealloc.args[1]

    # Substitute numeric assignments in new function's body
    newfnbody = subs(newfnbody, bkkeep.d_assign)
    preamble = subs(preamble, bkkeep.d_assign)
    prealloc = subs(prealloc, bkkeep.d_assign)
    debug && (println("------ _parse_newfnbody! ------");
        @show(newfnbody); println(); @show(preamble); println();
        @show(prealloc); println(); @show(bkkeep); println())

    # Include the assignement of indexed auxiliary variables
    defsprealloc = _defs_allocs!(prealloc, fnargs, bkkeep, false)
    preamble = subs(preamble, bkkeep.d_indx)
    defspreamble = Expr[preamble.args...]
    # Bring back substitutions
    newfnbody = subs(newfnbody, bkkeep.d_indx)

    # Define retvar; for scalar eqs is the last entry included in v_newvars
    bkkeep.retvar = length(fnargs) == 3 ? subs(bkkeep.v_newvars[end], bkkeep.d_indx) : fnargs[1]

    debug && (println("------ _defs_allocs! ------");
        @show(defsprealloc); println(); @show(preamble); println();
        @show(newfnbody); println(); @show(bkkeep); println())

    return defspreamble, defsprealloc, newfnbody, bkkeep
end

"""
`_rename_indexedvars(fnbody)`

Renames the indexed variables (using `Espresso.genname()`) that
exists in `fnbody`. Returns `fnbody` with the renamed variables
and a dictionary that links the new variables to the old
indexed ones.

"""
function _rename_indexedvars(fnbody)
    # Thanks to Andrei Zhabinski (@dfdx, see #31) for this
    # implementation, now slightly modified
    indexed_vars = findex(:(_X[_i...]), fnbody)
    st = Dict(ivar => genname() for ivar in indexed_vars)
    new_fnbody = subs(fnbody, st)
    return new_fnbody, Dict{Symbol, Expr}(v => k for (k, v) in st)
end

"""
`_newfnbody!(fnbody, fnargs, bkkeep)`

Returns a new (modified) body of the function, a priori unfolding
the expression graph (AST) as unary and binary calls, and updates
the bookkeeping structure bkkeep.

"""
function _newfnbody!(fnbody, fnargs, bkkeep::BookKeeping)
    # `fnbody` is assumed to be a `:block`-Expr
    newfnbody = Expr(:block,)

    # Local definition for possible args for `@threads` call
    local v_threads_args = (Expr(:., :Threads, QuoteNode(Symbol("@threads"))),
        Symbol("@threads"))

    # Magic happens HERE!!
    # Each line of fnbody (fnbody.args) is parsed separately
    for (i, ex) in enumerate(fnbody.args)
        if isa(ex, Expr)
            ex_head = ex.head

            # The `return` statement is treated separately
            if ex.head == :return
                @assert length(ex.args) == 1  # one sole returned value

                # Do nothing if it is `nothing` or if its in fnargs
                (ex.args[1] == :nothing || ex.args[1] in fnargs) && continue

                # Otherwise, create a new assignment, and process it
                new_aux_var = genname()
                ex = :($new_aux_var = $(ex.args[1]))
                ex_head = ex.head
            end

            # Treat `for` loops, `Threads.@threads for` and `if` blocks separately
            if ex_head == :block
                newblock = _newfnbody!(ex, fnargs, bkkeep)
                push!(newfnbody.args, newblock )

            elseif ex_head == :for
                push!(newfnbody.args, Expr(:for, ex.args[1]))
                loopbody = _newfnbody!( ex.args[2], fnargs, bkkeep )
                push!(newfnbody.args[end].args, loopbody)

            elseif ex_head == :macrocall
                # TODO: Deal with `@inbounds`
                # Deal with `Threads.@threads` and `@threads` cases
                if ex.args[1] in v_threads_args
                    push!(newfnbody.args, Expr(:macrocall, ex.args[1]))
                    # Although here `ex.args[2]` is a `LineNumberNode`,
                    # we add it to `newfnbody` because `@threads` call expressions
                    # require 3 args
                    push!(newfnbody.args[end].args, ex.args[2])
                    # Since `@threads` is called before a `for` loop, we deal
                    # with `ex.args[3]` as a `for` loop
                    push!(newfnbody.args[end].args, Expr(:for, ex.args[3].args[1]))

                    atthreadsbody = _newfnbody!( ex.args[3].args[2], fnargs, bkkeep )
                    push!(newfnbody.args[end].args[end].args, atthreadsbody)

                else
                    # If macro is not implemented, throw an `ArgumentError`
                    throw(ArgumentError("Macro $(ex.args[1]) is not yet implemented"))
                end

            elseif ex_head == :if
                # The first argument of an `if` expression is the condition, the
                # second one is the block the condition is true, and the third
                # is the `else`-block or an `ifelse`-block.
                push!(newfnbody.args, Expr(:if, ex.args[1]))

                for exx in ex.args[2:end]
                    if exx.head == :elseif
                        exxx = Expr(:block, Expr(:if, exx.args[1].args[2], exx.args[2:end]...))
                    else
                        exxx = exx
                    end

                    ifbody = _newfnbody!(exxx, fnargs, bkkeep)
                    push!(newfnbody.args[end].args, ifbody)
                end

            elseif ex_head == :(=) || ex_head == :call
                # Assignements or function calls
                ex_lhs = ex.args[1]
                ex_rhs = ex.args[2]

                # TODO: Include cases where the definition uses `similar`
                # Case of explicit declaration of Array or Vector
                if !isempty(findex(:(_Array{_TT...}), ex)) ||
                        !isempty(findex(:(Vector{_TT...}), ex))
                    inbookkeeping(ex_lhs, bkkeep) || push!(bkkeep.v_arraydecl, ex_lhs)
                    push!(newfnbody.args, ex)
                    continue
                end

                # Case of multiple assignments, e.g., ex = :((σ, ρ, β) = par)
                # Distinguishes if the RHS is the parameter or the dependent vars
                if isa(ex_lhs, Expr)
                    if ex_lhs.head == :tuple && length(fnargs) == 4 && isa(ex_rhs, Symbol)
                        # Construct new expressions
                        if ex_rhs == fnargs[3]  # params
                            for i in eachindex(ex_lhs.args)
                                ex_indx = Expr(:local, Expr(:(=), ex_lhs.args[i],
                                    Expr(:ref, ex_rhs, i)))
                                push!(newfnbody.args, ex_indx)
                            end

                        elseif ex_rhs == fnargs[2] # `x`
                            for i in eachindex(ex_lhs.args)
                                ex_indx = Expr(:(=), ex_lhs.args[i], Expr(:ref, ex_rhs, i))
                                push!(newfnbody.args, ex_indx)
                            end

                        else
                            throw(ArgumentError("Assignment no implemented in $(ex)"))
                        end
                    end
                    continue
                end

                # Unfold AST graph
                # nex = deepcopy(ex)
                # try
                #     nex = to_expr(ExGraph(simplify(ex)))
                # catch
                #     # copy `ex` as it is, if it is not "recognized"
                #     push!(newfnbody.args, ex)
                #     continue
                # end
                nex = to_expr(ExGraph(simplify(ex)))
                push!(newfnbody.args, nex.args[2:end]...)

                # Bookkeeping of indexed vars, to define assignements
                isindx_lhs = haskey(bkkeep.d_indx, ex_lhs)
                for nexargs in nex.args[2:end]
                    haskey(bkkeep.d_indx, nexargs.args[1]) && continue
                    vars_nex = find_vars(nexargs)

                    # any(haskey.(bkkeep.d_indx, Ref(vars_nex))) &&
                    #     !in(vars_nex[1], bkkeep.v_newvars) &&
                    (isindx_lhs || vars_nex[1] != ex_lhs) &&
                        !inbookkeeping(vars_nex[1], bkkeep) &&
                            push!(bkkeep.v_newvars, vars_nex[1])
                end

            elseif (ex_head == :local) || (ex_head == :continue) || (ex_head == :break)
                # If declared as `local`, or it is a `continue` or `break`, copy `ex` as it is.
                # In some cases using `local` helps performance. Very
                # useful for including (numeric) constants
                push!(newfnbody.args, ex)

            else
                # If not implemented, throw an `ArgumentError`
                throw(ArgumentError("$(ex_head) is not yet implemented; $(typeof(ex))"))
            end

        elseif isa(ex, LineNumberNode)
            continue  # Ignore `LineNumberNode`s

        else
            # In any other case, throw an `ArgumentError`
            throw(ArgumentError("$ex is not an `Expr`; $(typeof(ex))"))
        end
    end

    # Needed, if `newfnbody` consists of a single assignment (unary call)
    if newfnbody.head == :(=)
        newfnbody = Expr(:block, newfnbody)
    end

    return newfnbody
end

"""
_parse_newfnbody!(ex::Expr, preex::Expr, prealloc::Expr, bkkeep::BookKeeping, inloop::Bool)

Parses `ex` (the new body of the function) replacing the expressions
to use the mutating functions of TaylorSeries, and building the preamble
`preex` and `prealloc` expressions. This is done by traversing recursively (again)
the args of `ex`, updating the bookkeeping struct `bkkeep`, in particular
the fieldnames `v_newvars` and `d_assign`.

"""
function _parse_newfnbody!(ex::Expr, preex::Expr, prealloc::Expr, bkkeep::BookKeeping,
        inloop::Bool)
    # Numeric assignements to be deleted
    indx_rm = Int[]

    # Definition for possible args for `@threads` call
    local v_threads_args = (Expr(:., :Threads, QuoteNode(Symbol("@threads"))), Symbol("@threads"))

    # Magic happens HERE
    # Each line of `ex` (i.e., `ex.args`) is parsed separately
    for (i, aa) in enumerate(ex.args)
        # Treat for-loops, @threads macro, blocks and if-blocks separately
        # (through `inloop=true`)
        if (aa.head == :for)
            push!(preex.args, Expr(:for, aa.args[1]))
            push!(prealloc.args, Expr(:for, aa.args[1]))

            _parse_newfnbody!(aa, preex.args[end], prealloc.args[end], bkkeep, true)

        elseif (aa.head == :macrocall && aa.args[1] in v_threads_args)
            push!(preex.args, Expr(:macrocall, aa.args[1]))
            push!(preex.args[end].args, aa.args[2])
            push!(preex.args[end].args, Expr(:for, aa.args[3].args[1]))
            push!(prealloc.args, Expr(:macrocall, aa.args[1]))
            push!(prealloc.args[end].args, aa.args[2])
            push!(prealloc.args[end].args, Expr(:for, aa.args[3].args[1]))

            _parse_newfnbody!(aa.args[3], preex.args[end].args[end], prealloc.args[end].args[end],
                bkkeep, true)

        elseif (aa.head == :block)
            push!(preex.args, Expr(:block))
            push!(prealloc.args, Expr(:block))

            # `inloop` needs to be `false` here
            _parse_newfnbody!(aa, preex.args[end], prealloc.args[end], bkkeep, false)

        elseif (aa.head == :if)
            push!(preex.args, Expr(:if, aa.args[1]))
            push!(prealloc.args, Expr(:if, aa.args[1]))

            for exx in aa.args[2:end]
                push!(preex.args[end].args, Expr(:block))
                push!(prealloc.args[end].args, Expr(:block))
                _parse_newfnbody!(exx, preex.args[end].args[end],
                    prealloc.args[end].args[end], bkkeep, inloop)
            end

        elseif (aa.head == :(=))
            # Skip parsing the first for-loop expression (`inloop=true`), which is of
            # the form `:(i = ...)`
            i == 1 && inloop && continue

            # Main case
            aa_lhs = aa.args[1]
            aa_rhs = aa.args[2]

            # Replace expressions when needed, and bookkeeping
            if isa(aa_rhs, Expr)
                (aa_rhs.args[1] == :eachindex || aa_rhs.head == :(:)) && continue

                # Replace expression
                _replace_expr!(ex, preex, prealloc, i, aa_lhs, aa_rhs, bkkeep)

                # Remove new aa_lhs declaration from `ex`; it is already declared
                in(aa_lhs, bkkeep.v_arraydecl) && push!(indx_rm, i)

            elseif isa(aa_rhs, Symbol) # occurs when there is a simple assignment
                bb = subs(aa, Dict(aa_rhs => :(identity($aa_rhs))))
                bb_lhs = bb.args[1]
                bb_rhs = bb.args[2]

                # Replace expression
                _replace_expr!(ex, preex, prealloc, i, bb_lhs, bb_rhs, bkkeep)

            elseif isa(aa_rhs, Number)  # case of numeric values
                push!(bkkeep.d_assign, aa_lhs => aa_rhs)
                push!(indx_rm, i)
                if aa_lhs ∈ bkkeep.v_newvars
                    iis = findall(x->x==aa_lhs, bkkeep.v_newvars)
                    deleteat!(bkkeep.v_newvars, iis)
                end

            else # needed?
                error("Either $aa or $typeof(aa_rhs[2]) are not an `Expr`, `Symbol` or `Number`")
            end

        elseif aa.head == :local
            # If declared as `local`, copy `ex` as it is, and delete it from the recursion body
            push!(preex.args, aa)
            push!(prealloc.args, aa)
            push!(indx_rm, i)   # delete associated expr in body function

        else
            # Unrecognized head; pass the expression as it is
            push!(preex.args, aa)
            push!(prealloc.args, aa)
        end
    end

    # Delete trivial numeric assignement statements
    isempty(indx_rm) || deleteat!(ex.args, indx_rm)

    return nothing
end

"""
`_replace_expr!(ex::Expr, preex::Expr, , prealloc::Expr, i::Int, aalhs, aarhs,
    bkkeep::BookKeeping)`

Replaces the calls in `ex.args[i]`, and updates `preex` and `prealloc` with the
appropriate expressions, based on the the LHS (`aalhs`) and RHS (`aarhs`) of the
base assignment. The bookkeeping struct is updated (`v_newvars`) within `_replacecalls!`.
`d_indx` is used to bring back the indexed variables.

"""
function _replace_expr!(ex::Expr, preex::Expr, prealloc::Expr, i::Int,
        aalhs, aarhs, bkkeep::BookKeeping)
    # Replace calls
    fnexpr, def_fnexpr, auxfnexpr, def_alloc, aux_alloc = _replacecalls!(bkkeep, aarhs, aalhs)

    # Bring back indexed variables
    fnexpr     = subs(fnexpr, bkkeep.d_indx)
    def_fnexpr = subs(def_fnexpr, bkkeep.d_indx)
    auxfnexpr  = subs(auxfnexpr, bkkeep.d_indx)
    def_alloc  = subs(def_alloc, bkkeep.d_indx)
    aux_alloc  = subs(aux_alloc, bkkeep.d_indx)

    # Update `ex`, `preex` and `prealloc`
    if (fnexpr != def_fnexpr)
        ex.args[i] = fnexpr
    end

    if def_fnexpr.head != :nothing
        if def_fnexpr.head == :block
            append!(preex.args, def_fnexpr.args)
        else
            push!(preex.args, def_fnexpr)
        end
    end

    if def_alloc.head != :nothing
        if def_alloc.head == :block
            append!(prealloc.args, def_alloc.args)
        else
            push!(prealloc.args, def_alloc)
        end
    end

    # Same for the auxiliary expressions
    if auxfnexpr.head != :nothing
        if auxfnexpr.head == :block
            append!(preex.args, auxfnexpr.args)
        else
            push!(preex.args, auxfnexpr)
        end
        aux_alloc.head == :nothing || push!(prealloc.args, aux_alloc)
    end

    return nothing
end

"""
`_replacecalls!(bkkeep, fnold, newvar)`

Replaces the symbols of unary and binary calls
of the expression `fnold`, which defines `newvar`,
by the mutating functions in TaylorSeries.jl.
The vector `bkkeep.v_vars` is updated if new auxiliary
variables are introduced (bookkeeping).

"""
function _replacecalls!(bkkeep::BookKeeping, fnold::Expr, newvar::Symbol)
    ll = length(fnold.args)
    dcall = fnold.args[1]
    newarg1 = fnold.args[2]

    # If call is not in the mutating functions dictionaries, copy original one
    # to def_fnexpr and def_alloc, except if it is an Array/Vector definition
    if !( in(dcall, keys(TaylorSeries._dict_unary_calls)) ||
          in(dcall, keys(TaylorSeries._dict_binary_calls)) )

        fnexpr = :($newvar = $fnold)
        if isa(dcall, Expr) && (dcall.args[1] == :Array || dcall.args[1] == :Vector)
            def_fnexpr = Expr(:nothing)
            inbookkeeping(newvar, bkkeep) || push!(bkkeep.v_arraydecl, newvar)
        else
            def_fnexpr = fnexpr
            inbookkeeping(newvar, bkkeep) || push!(bkkeep.v_newvars, newvar)
        end
        def_alloc = fnexpr

        return fnexpr, def_fnexpr, Expr(:nothing), def_alloc, Expr(:nothing)
    end

    # Bookkeeping
    inbookkeeping(newvar, bkkeep) || push!(bkkeep.v_newvars, newvar)

    # Initializing
    def_alloc = Expr(:nothing)
    aux_alloc = Expr(:nothing)

    if ll == 2
        # Unary call
        # Replacements
        fnexpr, def_fnexpr, aux_fnexpr = TaylorSeries._dict_unary_calls[dcall]
        fnexpr = subs(fnexpr, Dict(:_res => newvar, :_arg1 => newarg1, :_k => :ord))

        def_alloc = :( _res = Taylor1($(def_fnexpr.args[2]), order) )
        def_alloc = subs(def_alloc,
            Dict(:_res => newvar, :_arg1 => :(constant_term($(newarg1))), :_k => :ord))

        def_fnexpr = Expr(:block,
            :(_res.coeffs[1] = $(def_fnexpr.args[2])),
            :(_res.coeffs[2:order+1] .= zero(_res.coeffs[1])) )
        def_fnexpr = subs(def_fnexpr,
            Dict(:_res => newvar, :_arg1 => :(constant_term($(newarg1))), :_k => :ord))
        # def_fnexpr = Expr(:block,
        #     :(_res[0] = $(def_fnexpr.args[2])),
        #     :(_res[1:order] .= zero(_res[0])) )
        # def_fnexpr = subs(def_fnexpr,
        #     Dict(:_res => newvar, :_arg1 => :(constant_term($(newarg1))), :_k => :ord))

        # Auxiliary expression
        if aux_fnexpr.head != :nothing
            newaux = genname()

            aux_alloc = :( _res = Taylor1($(aux_fnexpr.args[2]), order) )
            aux_alloc = subs(aux_alloc,
                Dict(:_res => newaux, :_arg1 => :(constant_term($(newarg1))), :_aux => newaux))

            aux_fnexpr = Expr(:block,
                :(_res.coeffs[1] = $(aux_fnexpr.args[2])),
                :(_res.coeffs[2:order+1] .= zero(_res.coeffs[1])) )
            aux_fnexpr = subs(aux_fnexpr,
                Dict(:_res => newaux, :_arg1 => :(constant_term($(newarg1))), :_aux => newaux))

            fnexpr = subs(fnexpr, Dict(:_aux => newaux))
            if newvar ∈ bkkeep.v_arraydecl
                push!(bkkeep.v_arraydecl, newaux)
            else
                push!(bkkeep.v_newvars, newaux)
            end
        end

    elseif ll == 3
        # Binary call; no auxiliary expressions needed
        newarg2 = fnold.args[3]

        # Replacements
        fnexpr, def_fnexpr, aux_fnexpr = TaylorSeries._dict_binary_calls[dcall]
        fnexpr = subs(fnexpr,
            Dict(:_res => newvar, :_arg1 => newarg1, :_arg2 => newarg2, :_k => :ord))

        def_alloc = :(_res = Taylor1($(def_fnexpr.args[2]), order) )
        def_alloc = subs(def_alloc, Dict(:_res => newvar, :_arg1 => :(constant_term($(newarg1))),
            :_arg2 => :(constant_term($(newarg2))), :_k => :ord) )

        def_fnexpr = Expr(:block,
            :(_res.coeffs[1] = $(def_fnexpr.args[2])),
            :(_res.coeffs[2:order+1] .= zero(_res.coeffs[1])) )
        def_fnexpr = subs(def_fnexpr,
            Dict(:_res => newvar, :_arg1 => :(constant_term($(newarg1))),
                :_arg2 => :(constant_term($(newarg2))), :_k => :ord))
        # def_fnexpr = Expr(:block,
        #     :(_res[0] = $(def_fnexpr.args[2])),
        #     :(_res[1:order] .= zero(_res[0])) )
        # def_fnexpr = subs(def_fnexpr,
        #     Dict(:_res => newvar, :_arg1 => :(constant_term($(newarg1))),
        #         :_arg2 => :(constant_term($(newarg2))), :_k => :ord))

    else
        # Recognized call, but not a unary or binary call; copy expression
        fnexpr = :($newvar = $fnold)
        def_fnexpr = fnexpr
        aux_fnexpr = Expr(:nothing)
    end

    return fnexpr, def_fnexpr, aux_fnexpr, def_alloc, aux_alloc
end

"""
`_defs_allocs!(preamble, fnargs, bkkeep, [inloop=false, ex_aux::Expr(:block,)])`

Returns a vector with expressions defining the auxiliary variables
in the preamble, and the declaration of the arrays. This function
may modify `bkkeep.d_indx` if new variables are introduced.
`bkkeep.v_preamb` is for bookkeeping the introduced variables.

"""
function _defs_allocs!(preamble::Expr, fnargs, bkkeep::BookKeeping,
    inloop::Bool, ex_aux::Expr = Expr(:block,))

    # Initializations
    defspreamble = Expr[]
    defs_alloc = Expr[]   # declaration of arrays

    # Local definition for possible args for `@threads` call
    local v_threads_args = (Expr(:., :Threads, QuoteNode(Symbol("@threads"))), Symbol("@threads"))

    for ex in preamble.args
        isa(ex, Expr) || throw(ArgumentError("$ex is not an `Expr`; $(typeof(ex))"))

        # Treat block, for loops, @threads for loops, if separately
        if (ex.head == :block)
            newdefspr = _defs_allocs!(ex, fnargs, bkkeep, inloop, ex_aux)

            append!(defspreamble, newdefspr)

        elseif (ex.head == :for)
            push!(ex_aux.args, ex.args[1])

            newdefspr = _defs_allocs!(ex.args[2], fnargs, bkkeep, true, ex_aux)

            append!(defspreamble, newdefspr)
            pop!(ex_aux.args)

        elseif (ex.head == :macrocall && ex.args[1] in v_threads_args)
            push!(ex_aux.args, ex.args[3].args[1])

            newdefspr = _defs_allocs!(ex.args[3].args[2], fnargs, bkkeep, true, ex_aux)

            append!(defspreamble, newdefspr)
            pop!(ex_aux.args)

        elseif (ex.head == :if)
            for exx in ex.args[2:end]
                newdefspr = _defs_allocs!(exx, fnargs, bkkeep, inloop, ex_aux)

                append!(defspreamble, newdefspr)
            end
            continue

        elseif (ex.head == :local)
            push!(defspreamble, subs(ex, bkkeep.d_indx))
            continue

        elseif (ex.head == :(=))
            # `ex.head` is a :(=) of some kind
            alhs = ex.args[1]
            arhs = ex.args[2]
            arhs = subs(ex.args[2], bkkeep.d_indx) # substitute updated vars in rhs

            # Outside of a loop
            if !inloop
                # `alhs` is declared as an array
                if in(alhs, bkkeep.v_arraydecl) || !in(alhs, bkkeep.v_preamb)
                    push!(defspreamble, subs(ex, bkkeep.d_decl))
                    push!(bkkeep.v_preamb, alhs)
                end

                continue
            end

            # Inside a loop
            if isindexed(alhs)
                # `var1` may be a vector or a matrix; declaring it is subtle
                var1 = alhs.args[1]
                while isa(var1, Expr)
                    var1 = var1.args[1]
                end
                (in(var1, fnargs) || in(var1, bkkeep.v_arraydecl) ||
                    in(alhs, bkkeep.v_preamb)) && continue

                # indices of var1
                d_subs = Dict(veexx.args[1] => veexx.args[2] for veexx in ex_aux.args)
                indx1 = alhs.args[2:end]
                ex_tuple = :( [$(indx1...)] )
                ex_tuple = subs(ex_tuple, d_subs)
                ex_tuple = :( size( $(Expr(:tuple, ex_tuple.args...)) ))
                exx = subs(copy(_DECL_ARRAY),
                    Dict( :__var1 => :($var1), :__var2 => :($ex_tuple)) )

                # Bookkeeping
                push!(defspreamble, exx.args...)
                push!(bkkeep.v_preamb, var1)
                continue

            elseif isindexed(arhs)
                # `alhs` is an aux indexed var, so something in `arhs` is indexed.
                in(alhs, bkkeep.v_preamb) && continue

                vars_indexed = findex(:(_X[_i...]), arhs)

                # NOTE: Uses the size of the var with more indices
                # to define the declaration of the new array.
                ii_indx = argmax( [length(find_indices(aa)[1]) for aa in vars_indexed] )
                var1 = vars_indexed[ii_indx].args[1]
                indx1 = vars_indexed[ii_indx].args[2:end]

                exx_indx = ones(Int, length(indx1))
                exx = subs(copy(_DECL_ARRAY),
                    Dict(:__var1 => :($alhs), :__var2 => :(size($var1))) )

                # Bookkeeping
                inbookkeeping(alhs, bkkeep) && alhs ∈ bkkeep.v_newvars &&
                        deleteat!(bkkeep.v_newvars, findall(x->x==alhs, bkkeep.v_newvars))

                push!(bkkeep.d_indx, alhs => :($alhs[$(indx1...)]) )
                push!(bkkeep.d_decl, alhs => :($alhs[$(exx_indx...)]))
                push!(bkkeep.v_arraydecl, alhs)
                push!(bkkeep.v_preamb, alhs)

                append!(defs_alloc, exx.args)
                continue

            else
                # `alhs` is not indexed nor an aux indexed var
                in(alhs, bkkeep.v_preamb) && continue
                vars_indexed = findex(:(_X[_i...]), arhs)
                if !isempty(vars_indexed)
                    ex = subs(ex, Dict(vv => :(one(S)) for vv in vars_indexed))
                end
                ex = subs(ex, bkkeep.d_decl)
                push!(defspreamble, ex)
                push!(bkkeep.v_preamb, alhs)
                continue
            end

        # # CHECK: Is the following needed?
        # else
        #     # If this is reached, other `ex.head` was encountered (e.g., `.=`);
        #     # we include this in the `defs_tmparrays`
        #     inloop || (push!(defs_tmparrays, subs(ex, bkkeep.d_decl)); continue)
        end

        exx = subs(ex, bkkeep.d_indx)
        inloop || push!(defspreamble, exx)
    end

    # Include allocations first in `defspreamble`
    pushfirst!(defspreamble, defs_alloc...)

    return defspreamble
end

"""
`_recursionloop(fnargs, bkkeep)`

Build the expression for the recursion-loop.

"""
function _recursionloop(fnargs, bkkeep::BookKeeping)
    ll = length(fnargs)

    if ll == 3
        rec_preamb = sanitize( :( $(fnargs[1]).coeffs[2] = $(bkkeep.retvar).coeffs[1] ) )
        rec_fnbody = sanitize( :( $(fnargs[1]).coeffs[ordnext+1] = $(bkkeep.retvar).coeffs[ordnext]/ordnext ) )

    elseif ll == 4
        bkkeep.retvar = fnargs[1]
        rec_preamb = sanitize(:(
            for __idx in eachindex($(fnargs[2]))
                $(fnargs[2])[__idx].coeffs[2] = $(bkkeep.retvar)[__idx].coeffs[1]
            end))
        rec_fnbody = sanitize(:(
            for __idx in eachindex($(fnargs[2]))
                $(fnargs[2])[__idx].coeffs[ordnext+1] =
                    $(bkkeep.retvar)[__idx].coeffs[ordnext]/ordnext
            end))

    else
        throw(ArgumentError(
        "Wrong number of arguments in the definition of the function $fn"))
    end

    return rec_preamb, rec_fnbody
end

"""
`@taylorize expr`

This macro `eval`s the function given by `expr` and defines a new
method of [`jetcoeffs!`](@ref) which is specialized on that
function. Integrating via [`taylorinteg`](@ref) of
[`lyap_taylorinteg`](@ref) after using the macro yields better performance.

See the [documentation](@ref taylorize) for more details and limitations.

!!! warning
    This macro is on an experimental stage; check the integration
    results carefully.

"""
macro taylorize( ex )
    nex1, nex2 = _make_parsed_jetcoeffs(ex)
    esc(quote
        $(ex)   # evals to calling scope the passed function
        $(nex1) # evals the new method of `jetcoeffs!`
        $(nex2) # evals the new method of `_allocate_jetcoeffs`
    end)
end


text:
using RCall, JLD
include("common_functions_for_tests.jl")

grid_n = [collect(10:500:5000); 5000] |> unique
grid_σ = [collect(0.1:0.3:0.9); 0.9]
grid_β = [0.5, 1., 2.]
param_grid = expand_grid(grid_n, grid_σ, grid_β) |>
    df -> rename(df, [:n, :σ, :β])

Nit = size(param_grid,1)

# Vnk_numerical_accuracy = GibbsTypePriors.Vnk_NGG.(param_grid[!,:n], 1, param_grid[!,:β], param_grid[!,:σ]) |> x -> accuracy_bits.(x)

Vnk_numerical_accuracy = load("test/graphical_tests/saves_for_graphical_tests/accuracy_Vn1.jld", "Vnk_numerical_accuracy")
R"library(tidyverse)"
# R"library(ggpubr)"
# R"p = tibble(Vnk = $Vnk_numerical_accuracy, n = $(param_grid[!,:n]), s = $(param_grid[!,:σ]), b = $(param_grid[!,:β])) %>%
#     ggplot(aes(x = n, y = Vnk, colour = s, group = s)) +
#     theme_bw() +
#     facet_wrap(~b) +
#     geom_point() +
#     # geom_line() +
#   stat_smooth(aes(fill = s, color = s), method = 'lm', formula = formula) +
#   stat_regline_equation(
#     aes(label =  paste(..eq.label.., ..adj.rr.label.., sep = '~~~~')),
#     formula = formula
#   ) +
#     geom_hline(yintercept = 64, linetype = 'dotted', colour = 'red') +
#     viridis::scale_colour_viridis()"
#
# R"ggpar(p, palette = 'jco')"
R"my.formula <- y ~ x"
R"p = tibble(Vnk = $Vnk_numerical_accuracy, n = $(param_grid[!,:n]), s = $(param_grid[!,:σ]), b = $(param_grid[!,:β])) %>%
ggplot(aes(x=n, y = Vnk, group = s, colour = factor(s))) +
  geom_smooth(method='lm', formula= y~x) +
  facet_wrap(~b) +
  theme_classic() +
  geom_point() +
  viridis::scale_colour_viridis(discrete = T, name = expression(sigma)) +
  ggpmisc::stat_poly_eq(formula = my.formula,
               aes(label = paste(..eq.label.., ..rr.label.., sep = '~~~')),
               parse = TRUE) +
  geom_hline(yintercept = 64, linetype = 'dotted', colour = 'red')
"

R"pdf('test/graphical_tests/figures_graphical_tests/accuracy_Vn1.pdf', height = 5, width = 5)
    plot(p)
dev.off()"

max_stars_repo_name: konkam/GibbsTypePriors

max_issues_repo_name: konkam/GibbsTypePriors

max_forks_repo_name: konkam/GibbsTypePriors


text:
Statistics.mean(vartype::InferenceVariable, args...; kwargs...) = mean(getManifold(vartype), args...; kwargs...)
Statistics.std(vartype::InferenceVariable, args...; kwargs...) = std(getManifold(vartype), args...; kwargs...)
Statistics.var(vartype::InferenceVariable, args...; kwargs...) = var(getManifold(vartype), args...; kwargs...)

function Statistics.cov(vartype::InferenceVariable, ptsArr::AbstractVector; basis::Manifolds.AbstractBasis = Manifolds.DefaultOrthogonalBasis(), kwargs...)
  return cov(getManifold(vartype), ptsArr; basis, kwargs... )
end

# To replace calcCovarianceBasic
function calcStdBasicSpread(vartype::InferenceVariable, ptsArr::Vector{P}) where P
  σ = std(vartype, ptsArr)

  #if no std yet, set to 1
  msst = 1e-10 < σ ? σ : 1.0
  return msst
end

#TODO consolidate
function calcMeanCovar(vari::DFGVariable, solvekey=:default)

  pts = getSolverData(vari, solvekey).val
  μ = mean(getManifold(vari), pts)
  Σ = cov(getVariableType(vari), pts)
  return μ, Σ
end


text:
module Execute

export execute, ExecutionArgs
using ..Types

function execute(args::ExecutionArgs)
  return Dict("TODO" => "TODO")
end

end


text:
module DiscreteModelPortionsTests

using Gridap.Arrays
using Gridap.ReferenceFEs
using Gridap.Geometry

domain = (0,2,0,2,0,2)
partition = (10,10,10)
oldmodel = CartesianDiscreteModel(domain,partition)

const R = 0.7

function is_in(coords)
  n = length(coords)
  x = (1/n)*sum(coords)
  d = x[1]^2 + x[2]^2 - R^2
  d < 0
end

oldgrid = get_grid(oldmodel)
oldcell_to_coods = get_cell_coordinates(oldgrid)
oldcell_to_mask = lazy_map(is_in,oldcell_to_coods)

model = DiscreteModelPortion(oldmodel,oldcell_to_mask)
test_discrete_model(model)

#using Gridap.Visualization
#
#writevtk(model,"model")

end # module


text:
# functions related to F distribution

import .RFunctions:
    fdistpdf,
    fdistlogpdf,
    fdistcdf,
    fdistccdf,
    fdistlogcdf,
    fdistlogccdf,
    fdistinvcdf,
    fdistinvccdf,
    fdistinvlogcdf,
    fdistinvlogccdf

# pdf for numbers with generic types
fdistpdf(ν1::Real, ν2::Real, x::Number) = sqrt((ν1 * x)^ν1 * ν2^ν2 / (ν1 * x + ν2)^(ν1 + ν2)) / (x * beta(ν1 / 2, ν2 / 2))

# logpdf for numbers with generic types
fdistlogpdf(ν1::Real, ν2::Real, x::Number) = (ν1 * log(ν1 * x) + ν2 * log(ν2) - (ν1 + ν2) * log(ν1 * x + ν2)) / 2 - log(x) - lbeta(ν1 / 2, ν2 / 2)


text:
# Sytem tests

# Test constructors
@test size(LinearSystem(ones(2,2)).A) == (2,2)
@test size(LinearSystem(1.0).A) == (1,1)
@test size(LinearSystem(1.0, 1.0).Q) == (1,1)
@test size(LinearSystem(eye(3), eye(3)).A) == (3,3)
@test_throws DimensionMismatch LinearSystem(eye(2), eye(3))
@test_throws DimensionMismatch LinearSystem(eye(2,3), eye(2))
@test_throws DimensionMismatch predict(DiscreteState(ones(2)),
                                       LinearSystem(1.0), 1)

# Test state transtion matrix methods
lin_sys = LinearSystem(ones(2,2), eye(2))
state = DiscreteState(2*ones(2))
@test state_transition_matrix(lin_sys, state, 10) == ones(2,2)^10
lin_sys = LinearSystem(2.0, 1.0)
state = ContinuousState([1.0])
@test state_transition_matrix(lin_sys, state, 2.0) == exp(4.0)*ones(1,1)

# Test consistency between predict and state_transition_matrix
predict(DiscreteState(2.0), lin_sys, 2).x ==
    state_transition_matrix(lin_sys, DiscreteState(2.0), 2) * 2.0
predict(UncertainDiscreteState(2.0, 1.0), lin_sys, 2).x ==
    state_transition_matrix(lin_sys, UncertainDiscreteState(2.0, 1.0), 2) * 2.0

# Test discrete prediction methods
lin_sys = LinearSystem(ones(2,2), eye(2))
@test predict(DiscreteState(ones(2)), lin_sys, 2) == DiscreteState(4*ones(2), 2)
@test predict(UncertainDiscreteState(ones(2), eye(2)), lin_sys, 2) ==
    UncertainDiscreteState(4*ones(2), eye(2)+10*ones(2,2), 2)

# Test continuous prediction methods
lin_sys = LinearSystem([[0.0, 1.0] [-1.0, 0.0]], 2*eye(2))
@test predict(ContinuousState(ones(2)), lin_sys, 2.0) ==
    ContinuousState(expm(lin_sys.A*2)*ones(2), 2.0)
result = predict(UncertainContinuousState(ones(2), eye(2)), lin_sys, 2.0)
@test result.x == expm(lin_sys.A*2)*ones(2)
@test isapprox(result.P, [[5.0, 0.0] [0.0, 5.0]], rtol=1e-6)
@test predict(ContinuousState(ones(2)), lin_sys, 0.0) == ContinuousState(ones(2))
@test predict(UncertainContinuousState(ones(2), eye(2)), lin_sys, 0.0) ==
    UncertainContinuousState(ones(2), eye(2))

# Test in-place discrete prediction methods
lin_sys = LinearSystem(ones(2,2), eye(2))
state = DiscreteState(ones(2))
predict!(state, lin_sys, 2)
@test state == DiscreteState(4*ones(2), 2)
state = UncertainDiscreteState(ones(2), eye(2))
predict!(state, lin_sys, 2)
@test state == UncertainDiscreteState(4*ones(2), eye(2)+10*ones(2,2), 2)

# Test in-place continuous prediction methods
lin_sys = LinearSystem([[0.0, 1.0] [-1.0, 1.0]], 2*eye(2))
state = ContinuousState(ones(2))
predict!(state, lin_sys, 0.1)
@test state == ContinuousState(expm(lin_sys.A*0.1)*ones(2), 0.1)
state = UncertainContinuousState(ones(2), eye(2))
predict!(state, lin_sys, 0.1)
@test state.x == expm(lin_sys.A*0.1)*ones(2)
@test isapprox(state.P, [[1.20075, -0.0117417] [-0.0117417, 1.44201]],rtol=1e-5)

# Test simulation methods
lin_sys = LinearSystem([[0.0, 1.0] [-1.0, 1.0]], 2*eye(2))
continuous_state = ContinuousState(ones(2))
uncertain_discrete_state = UncertainDiscreteState(ones(2), eye(2))
@test typeof(simulate(continuous_state, lin_sys, 1.0)) <: ContinuousState
@test typeof(simulate(uncertain_discrete_state, lin_sys, 1)) <: DiscreteState

# Test nonlinear constructors
discrete_nl_fcn(t, x::Vector) = x.^2;
discrete_nl_jac(t, x::Vector) = diagm(2*x)
NonlinearSystem(discrete_nl_fcn, discrete_nl_jac, 1.0)
NonlinearSystem(discrete_nl_fcn, discrete_nl_jac, eye(2))
NonlinearSystem(LinearSystem(ones(2,2), eye(2)))

# Test state transtion matrix methods
continuous_nl_fcn(t, x::Vector) = -x.^2;
continuous_nl_jac(t, x::Vector) = -diagm(2*x)
continuous_state = ContinuousState(ones(2))
discrete_state = DiscreteState(ones(2))
nonlin_sys = NonlinearSystem(continuous_nl_fcn, continuous_nl_jac, eye(2))
@test size(state_transition_matrix(nonlin_sys, continuous_state, 1.0)) == (2,2)
@test size(state_transition_matrix(nonlin_sys, discrete_state, 1)) == (2,2)

# Test nonlinear discrete prediction methods
nonlin_sys = NonlinearSystem(discrete_nl_fcn, discrete_nl_jac, eye(3))
@test predict(UncertainDiscreteState([0.0, 1.0, 2.0], eye(3)), nonlin_sys, 1) ==
    UncertainDiscreteState([0.0, 1.0, 4.0], diagm([1.0, 5.0, 65.0]), 1)
@test predict(DiscreteState([0.0, 1.0, 2.0]), nonlin_sys, 1) ==
    DiscreteState([0.0, 1.0, 4.0], 1)

# Test nonlinear continuous prediction methods
continuous_nl_fcn(t, x::Vector) = [x[2], 0.0]
continuous_nl_jac(t, x::Vector) = [0.0 1.0; 0.0 0.0]
nonlin_sys = NonlinearSystem(continuous_nl_fcn, continuous_nl_jac, eye(2))
@test predict(ContinuousState(1.0), nonlin_sys, 0.0) == ContinuousState(1.0)
@test predict(UncertainContinuousState(1.0, 0.0), nonlin_sys, 0.0) ==
    UncertainContinuousState(1.0, 0.0)
@test isapprox(predict(ContinuousState([0.0, 1.0]), nonlin_sys,  2).x,
    ContinuousState([2.0, 1.0], 2.0).x)

# Test consistency between linear and nonlinear prediction
lin_sys = LinearSystem([0.0 1.0; 0.0 0.0], eye(2))
@test isapprox(
    predict(UncertainContinuousState([0.0, 1.0], eye(2)), nonlin_sys, 2.0).P,
    predict(UncertainContinuousState([0.0, 1.0], eye(2)), lin_sys, 2.0).P)

# Test simulation methods
continuous_nl_fcn(t, x::Vector) = [x[2], 0.0]
continuous_nl_jac(t, x::Vector) = [0.0 1.0; 0.0 0.0]
nonlin_sys = NonlinearSystem(continuous_nl_fcn, continuous_nl_jac, eye(2))
@test typeof(simulate(ContinuousState([0.0, 1.0]), nonlin_sys, 1.0)) <:
    ContinuousState


text:
# Test internal utilities functions

make_pdmat(n::Int) = (g = randn(n, n); NMF.adddiag!(g'g, 0.1))

## adddiag!

a0 = rand(3, 3)
a = copy(a0)

NMF.adddiag!(a, 0.)
@test a == a0

NMF.adddiag!(a, 2.5)
@test a == a0 + 2.5I

## normalize1!

a = rand(5)
NMF.normalize1!(a)
@test sum(a) ≈ 1.0

## normalize1_cols!

a = rand(5, 6)
NMF.normalize1_cols!(a)
@test vec(sum(a, dims=1)) ≈ ones(6)

## projectnn!

a0 = randn(5, 5)
a = copy(a0)
NMF.projectnn!(a)
@test a == max.(a0, 0.0)

## posneg!

a = randn(5, 5)
ac = copy(a)
ap = zeros(size(a))
an = zeros(size(a))

NMF.posneg!(a, ap, an)
@test a == ac
@test ap == max.(a, 0.0)
@test an == max.(-a, 0.0)

## pdsolve!

A = make_pdmat(5)
X = rand(5, 3)
Y = A * X
NMF.pdsolve!(A, Y)
@test X ≈ Y

## pdrsolve!

B = make_pdmat(5)
X = rand(4, 5)
Y = X * B
Xr = zeros(4, 5)
NMF.pdrsolve!(Y, B, Xr)
@test Xr ≈ X


text:
module SplineFEEC

using CompactBases # Used for basis creation
using DomainSets
import DomainSets: (..), ×

using LinearAlgebra
using StaticArrays
using FastGaussQuadrature
using Kronecker
import Kronecker: ⊗
using FFTW

using Plots
using Cubature

export (..), ×, ⊗, ⊙, adjoint,
        PeriodicBSpline, DirichletBSpline, PeriodicBSplineTensorProductBasis, DirichletBSplineTensorProductBasis,
        mass_matrix, stiffness_matrix, L2_prod, col_mat, L2_proj, weighted_mass_matrix, interior_product_matrix,
        grad, curl, div, ∇, disc_diff, disc_grad, disc_curl, disc_div,
        Basis, TensorProductBasis, BSplineTensorProductBasis,
        DirichletMultipatchBSpline, mass_matrix_d, L2_prod_d, L2_proj_d, col_mat_d

include("bases.jl")
include("circulant.jl")
include("util.jl")
include("operator.jl")
include("multipatch.jl")
include("examples.jl")

end


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

include("../../src/Sparso.jl")
include("../../src/simple-show.jl")
using Sparso

set_options(SA_ENABLE, SA_VERBOSE, SA_USE_SPMP, SA_CONTEXT)

include("utils.jl")

function spmatmul_witheps{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti}, eps;
                         sortindices::Symbol = :sortcols)
    mA, nA = size(A)
    mB, nB = size(B)
    nA==mB || throw(DimensionMismatch())

    colptrA = A.colptr; rowvalA = A.rowval; nzvalA = A.nzval
    colptrB = B.colptr; rowvalB = B.rowval; nzvalB = B.nzval
    # TODO: Need better estimation of result space
    nnzC = min(mA*nB, length(nzvalA) + length(nzvalB))
    colptrC = Array(Ti, nB+1)
    rowvalC = Array(Ti, nnzC)
    nzvalC = Array(Tv, nnzC)

    @inbounds begin
        ip = 1
        xb = zeros(Ti, mA)
        x  = zeros(Tv, mA)
        for i in 1:nB
            if ip + mA - 1 > nnzC
                resize!(rowvalC, nnzC + max(nnzC,mA))
                resize!(nzvalC, nnzC + max(nnzC,mA))
                nnzC = length(nzvalC)
            end
            colptrC[i] = ip
            for jp in colptrB[i]:(colptrB[i+1] - 1)
                nzB = nzvalB[jp]
                j = rowvalB[jp]
                for kp in colptrA[j]:(colptrA[j+1] - 1)
                    nzC = nzvalA[kp] * nzB
                    k = rowvalA[kp]
                    if xb[k] != i
                        rowvalC[ip] = k
                        ip += 1
                        xb[k] = i
                        x[k] = nzC
                    else
                        x[k] += nzC
                    end
                end
            end
            idx = colptrC[i]
            for vp in colptrC[i]:(ip - 1)
                col = rowvalC[vp]
                if col == i || x[col] > eps
                  nzvalC[idx] = x[col]
                  rowvalC[idx] = col
                  idx += 1
                end
            end
            ip = idx
        end
        colptrC[nB+1] = ip
    end

    deleteat!(rowvalC, colptrC[end]:length(rowvalC))
    deleteat!(nzvalC, colptrC[end]:length(nzvalC))

    # The Gustavson algorithm does not guarantee the product to have sorted row indices.
    Cunsorted = SparseMatrixCSC(mA, nB, colptrC, rowvalC, nzvalC)
    C = Base.SparseMatrix.sortSparseMatrixCSC!(Cunsorted, sortindices=sortindices)
    return C
end

function foo()
    eps = 1e-5

    m = 10
    X = generate_symmetric_nonzerodiagonal_sparse_matrix(m)
    Y = generate_symmetric_nonzerodiagonal_sparse_matrix(m)
    # X = SparseMatrixCSC{Float64, Int32}((eMax*speye(m) - X)/(eMax - eMin))

    for i = 1:2
        A = spmatmul_witheps(X, X', eps)
        B = spmatmul_witheps(X', X, eps)
        C = spmatmul_witheps(X, Y', eps)
        D = spmatmul_witheps(X, Y, eps)
    end
end

@acc foo()

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antddrawer

"""
    'feffery'_antddrawer(;kwargs...)
    'feffery'_antddrawer(children::Any;kwargs...)
    'feffery'_antddrawer(children_maker::Function;kwargs...)

An AntdDrawer component.

Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional)
- `id` (String; optional)
- `className` (String; optional)
- `closable` (Bool; optional)
- `containerId` (String; optional)
- `destroyOnClose` (Bool; optional)
- `forceRender` (Bool; optional)
- `height` (Real; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `mask` (Bool; optional)
- `maskClosable` (Bool; optional)
- `placement` (a value equal to: 'left', 'right', 'top', 'bottom'; optional)
- `style` (Dict; optional)
- `title` (String; optional)
- `visible` (Bool; optional)
- `width` (Real; optional)
- `zIndex` (Real; optional)
"""
function 'feffery'_antddrawer(; kwargs...)
        available_props = Symbol[:children, :id, :className, :closable, :containerId, :destroyOnClose, :forceRender, :height, :loading_state, :mask, :maskClosable, :placement, :style, :title, :visible, :width, :zIndex]
        wild_props = Symbol[]
        return Component("'feffery'_antddrawer", "AntdDrawer", "feffery_antd_components", available_props, wild_props; kwargs...)
end

'feffery'_antddrawer(children::Any; kwargs...) = 'feffery'_antddrawer(;kwargs..., children = children)
'feffery'_antddrawer(children_maker::Function; kwargs...) = 'feffery'_antddrawer(children_maker(); kwargs...)

max_stars_repo_name: limmyforget/feffery-antd-components

max_issues_repo_name: limmyforget/feffery-antd-components

max_forks_repo_name: limmyforget/feffery-antd-components


text:
using MLFS
using PrettyPrint
using Setfield
using MLFS.HM
using MLStyle
import JSON
import FileIO

export smlfsCompile

macro _check_prune(a)
    quote
        let ret = pruneWithVarCheck(collect_unsolved, $a)
            !isempty(unsolved_vars) &&
                throw(MLError(ln, UnsolvedTypeVariables(unsolved_vars)))
            ret
        end
    end |> esc
end
function mkPostInfer(g::GlobalTC)
    pruneWithVarCheck = g.tcstate.pruneWithVarCheck
    unsolved_vars = Set{Var}()
    function collect_unsolved(v::Var)
        push!(unsolved_vars, v)
    end

    function postInfer(ln::LineNumberNode, root::IR.Decl)
        @match root begin
            IR.Perform(_) => IR.gTransCtx(postInfer, ln, root)
            IR.Assign(a, t, e) =>
                IR.Assign(a, @_check_prune(t), postInfer(ln, e))
        end
    end

    function postInfer(ln::LineNumberNode, root::IR.Expr)
        @match root begin
            IR.Expr(ln, ty, expr) =>
                let expr = postInfer(ln, expr),
                    ty = ty === nothing ? nothing : @_check_prune(ty)

                    IR.Expr(ln, ty, expr)
                end
        end
    end

    function postInfer(ln::LineNumberNode, root::IR.ExprImpl)
        @match root begin
            IR.EIm(expr, t, insts) =>
                let expr = postInfer(ln, expr),
                    t = @_check_prune(t)

                    IR.EApp(expr, instanceResolve(g, insts, t, ln))
                end
            IR.ETypeVal(t) => IR.ETypeVal(@_check_prune(t))
            _ => IR.gTransCtx(postInfer, ln, root)
        end
    end
    postInfer
end

const preDefinedTypeEnv = Store{Symbol, HMT}()[[
    a => T(Nom(a)) for a in Symbol[
        :i8, :i16, :i32, :i64,
        :i16, :i32, :i64,
        :f16, :f32, :f64,
        :char, :str, :bool,
        :Module,
        :field,
        :namespace,
        :Type
    ]
]]

const Path = String

moduleGensym(moduleName::Symbol) = Symbol(:|, moduleName, :|)
moduleExportName(moduleName::Symbol) = Symbol(moduleName)

PostInfer = Base._return_type(mkPostInfer, (GlobalTC, ))

const Signature = Pair{Vector{Symbol}, Vector{Pair{Nom, Vector{InstRec}}}}

"""
smlfs-c A=a.mlfs B=b.mlfs C=c.mlfs -sigdir -o pkgA
"""
function smlfsCompile(
    srcFiles::Vector{Path},
    sigFiles::Vector{Path},
    outName::String,
    outDir::String="./")

    loadedModules = Set{Symbol}()
    g = empty(GlobalTC)
    loadSigs(sigFiles, loadedModules, g)

    decls, signature = loadSrcs(srcFiles, loadedModules, g)
    show_hints(g)

    declsJSON = toVec(Vector{IR.Decl}, decls)
    open(joinpath(outDir, "$outName.mlfso"), "w") do f
        write(f, JSON.json(declsJSON))
    end

    sigJSON = toVec(Signature, signature)
    open(joinpath(outDir, "$outName.mlfsa"), "w") do f
        write(f, JSON.json(sigJSON))
    end
end

const FieldClass = Nom(:field)
const NamespaceClass = Nom(:namespace)

function loadSigs(sigFiles::Vector{Path}, loadedModules::Set{Symbol}, g::GlobalTC)
    for each in sigFiles
        sig = open(each) do f
            srcCode = read(f, String)
            fromVec(
                Signature,
                JSON.Parser.parse(srcCode))
        end
        sig === nothing && error("unrecognised signature file at $each.")
        newModuleNames, implicits = sig.value
        reloadedModules = intersect(loadedModules,  newModuleNames)
        isempty(reloadedModules) ||
            error("module name conflicts: duplicate modules with the same names: $reloadedModules")
        for (tConsHead, insts) in implicits
            append!(g.globalImplicits[tConsHead], insts)
        end
    end
end

function loadSrcs(srcFiles::Vector{Path}, loadedModules::Set{Symbol}, g::GlobalTC)
    globalImplicitDeltas = g.globalImplicitDeltas
    decls = IR.Decl[]
    moduleNames = Symbol[]
    for path in srcFiles
        mR = open(path) do f
            srcCode = read(f, String)
            mR = fromVec(Surf.ModuleRecord, JSON.Parser.parse(srcCode))
            if mR === nothing
                error("Invalid surface file at $path")
            end
            mR = mR.value
            moduleName = mR.name
            moduleName in loadedModules &&
                error("module name conflicts: duplicate modules with the same names: $moduleName")
            push!(loadedModules, moduleName)
            push!(moduleNames, moduleName)
            append!(
                decls,
                loadModule(moduleName, mR.imports, mR.decls, path, g))
        end
    end

    prune = g.tcstate.prune
    pruneRel(x::InstRec) = begin
        x.isPruned || (x.t = prune(x.t))
        x
    end
    signature = Pair{Nom, Vector{InstRec}}[
        nom =>
            pruneRel.(g.globalImplicits[nom][end-n+1:end])
        for (nom, n) in globalImplicitDeltas
    ]
    decls, (moduleNames=>signature)
end

function loadModule(
    moduleName::Symbol,
    imports::Vector{Pair{Symbol, Symbol}},
    decls::Vector{Surf.Decl},
    path::Path,
    g::GlobalTC)

    tcstate = g.tcstate
    prune = tcstate.prune

    g = @set g.moduleName = moduleName
    l = empty(LocalTC)
    l = @set l.typeEnv = preDefinedTypeEnv
    for (importModule, alias) in imports
        moduleType = App(Nom(:Module), Nom(importModule))
        l = @set l.typeEnv = l.typeEnv[alias => moduleType]
        l = @set l.symmap = l.symmap[alias => moduleGensym(importModule)]
    end

    results, localTC = inferModule(g, l, decls)
    ln = LineNumberNode(1, path)
    globalImplicits = get!(g.globalImplicits, FieldClass) do
        InstRec[]
    end

    globalImplicitDeltas = g.globalImplicitDeltas
    get!(globalImplicitDeltas, FieldClass, 0)

    moduleType = App(Nom(:Module), Nom(moduleName))

    function mkField(usersym::Symbol, genexp::ExprImpl, pruned_type::HMT)
        local instanceTy
        instanceTy = App(FieldClass, Tup(HMT[
            moduleType, Nom(usersym), pruned_type]))
        push!(globalImplicits, InstRec(instanceTy, ln, genexp, true))
    end

    fields = Dict{Symbol, Tuple{Gensym, HMT}}()
    for (usersym, type) in localTC.typeEnv.unbox
        if haskey(fields, usersym) || type === get(preDefinedTypeEnv, usersym, nothing)
            # shadowing
            continue
        end

        gensym = get(localTC.symmap, usersym, nothing)
        method_gensym = symgen(g)
        type = prune(type)
        genexp = @match type begin
            T(type) => IR.ETypeVal(type)

            # TODO: what does this mean?
            if gensym === nothing end => IR.EInt(0, 8)

            _ => IR.EVar(gensym::Symbol)
        end
        fields[usersym] = (gensym, type)
        genexp = IR.Expr(ln, nothing, IR.EFun(:_, IR.Expr(ln, nothing, genexp)))
        push!(results, IR.Assign(method_gensym, Arrow(HM.⊤, type), genexp))
        mkField(usersym, IR.EVar(method_gensym), prune(type))
    end
    globalImplicitDeltas[FieldClass] += length(fields)

    globalImplicits = get!(g.globalImplicits, Nom(:Module)) do
        InstRec[]
    end
    push!(globalImplicits, InstRec(moduleType, ln, IR.ETypeVal(moduleType), true))

    get!(globalImplicitDeltas, Nom(:Module), 0)
    globalImplicitDeltas[Nom(:Module)] += 1

    globalImplicits = get!(g.globalImplicits, NamespaceClass) do
        InstRec[]
    end

    get!(globalImplicitDeltas, NamespaceClass, 0)
    instanceTy = App(
        NamespaceClass,
        Tup(
            HMT[
                moduleType,
                Tup(HMT[Tup([Nom(k), Nom(gensym), t]) for (k, (gensym, t)) in fields])]))

    push!(globalImplicits, InstRec(instanceTy, ln, IR.ETypeVal(instanceTy), true))
    globalImplicitDeltas[NamespaceClass] += 1

    # type class, etc
    postInfer = mkPostInfer(g)
    for i in eachindex(results)
        results[i] = postInfer(ln, results[i])
    end

    push!(results,
        IR.Assign(moduleGensym(moduleName), moduleType,
            IR.Expr(ln, nothing, IR.ETypeVal(moduleType))),
        IR.Assign(moduleExportName(moduleName), instanceTy,
            IR.Expr(ln, nothing, IR.ETypeVal(instanceTy))),
    )
    results
end