text:
using Random

mutable struct Lattice2D
	L              ::Int
	p_occupation   ::Float64
	occupied       ::BitArray{2}
	cluster_ids    ::Matrix{Int}
	equivalent_ids ::Vector{Int}
	clusters       ::Dict{Int, Set{NTuple{2, Int}}}
	is_percolating ::Bool
	rng            ::AbstractRNG

	function Lattice2D(L::Int, p_occupation::Float64, seed::Int)
		rng            = MersenneTwister(seed)
		occupied       = rand(rng, L, L) .< p_occupation
		cluster_ids    = zeros(Int, size(occupied))
		clusters       = Dict{Int, Set{NTuple{2, Int}}}()
		equivalent_ids = collect(1:L^2)
		is_percolating = false

		lattice = new(
			L,
			p_occupation,
			occupied,
			cluster_ids,
			equivalent_ids,
			clusters,
			is_percolating,
			rng
		)

		label_clusters!(lattice)
		create_cluster_sets!(lattice)
		is_lattice_percolating!(lattice)

		return lattice
	end
end

function uf_find(lattice::Lattice2D, cluster_id::Int)
	while lattice.equivalent_ids[cluster_id] ≠ cluster_id
		cluster_id = lattice.equivalent_ids[cluster_id]
	end
	return cluster_id
end

function uf_union!(lattice::Lattice2D, cluster_id_1::Int, cluster_id_2::Int)
	lattice.equivalent_ids[uf_find(lattice, cluster_id_1)] = uf_find(lattice, cluster_id_2)
	return uf_find(lattice, cluster_id_1)
end

function label_clusters!(lattice::Lattice2D)
	new_cluster_id = 1
	for j in 1:lattice.L, i in 1:lattice.L
		if lattice.occupied[i, j]
			left_cluster_id = (j == 1 ? 0 : lattice.cluster_ids[i, j - 1])
			up_cluster_id   = (i == 1 ? 0 : lattice.cluster_ids[i - 1, j])

			if left_cluster_id == up_cluster_id == 0
				lattice.cluster_ids[i, j] = new_cluster_id
				new_cluster_id += 1
			elseif left_cluster_id == 0 || up_cluster_id == 0
				lattice.cluster_ids[i, j] = maximum((left_cluster_id, up_cluster_id))
			else
				lattice.cluster_ids[i, j] = uf_union!(lattice, left_cluster_id, up_cluster_id)
			end
		end
	end

	for j in 1:lattice.L, i in 1:lattice.L
		if lattice.occupied[i, j]
			lattice.cluster_ids[i, j] = uf_find(lattice, lattice.cluster_ids[i, j])
		end
	end
end

function create_cluster_sets!(lattice::Lattice2D)
	clusters    = Dict{Int, Set{NTuple{2, Int}}}()
	for j in 1:lattice.L, i in 1:lattice.L
		if lattice.occupied[i, j]
			site = (i, j)
			if haskey(clusters, lattice.cluster_ids[i, j])
				union!(clusters[lattice.cluster_ids[i, j]], [site])
			else
				clusters[lattice.cluster_ids[i, j]] = Set([site])
			end
		end
	end

	lattice.clusters = clusters
end

function is_lattice_percolating!(lattice::Lattice2D)
	L = lattice.L
	top_boundary    = Set{NTuple{2, Int}}([(1, j) for j in 1:L])
	bottom_boundary = Set{NTuple{2, Int}}([(L, j) for j in 1:L])
	left_boundary   = Set{NTuple{2, Int}}([(i, 1) for i in 1:L])
	right_boundary  = Set{NTuple{2, Int}}([(i, L) for i in 1:L])

	for cluster in values(lattice.clusters)
		if length(cluster ∩ top_boundary) > 0 && length(cluster ∩ bottom_boundary) > 0
			lattice.is_percolating = true
			break
		elseif length(cluster ∩ left_boundary) > 0 && length(cluster ∩ right_boundary) > 0
			lattice.is_percolating = true
			break
		end
	end
end

max_stars_repo_name: cameronperot/percolation-classification

max_issues_repo_name: cameronperot/percolation-classification

max_forks_repo_name: cameronperot/percolation-classification


text:
using AbstractPlotting
using CairoMakie

array = [(t, y₁, y₂) -> (y₁ * (1.0 + y₂)),
         (t, y₁, y₂) -> (y₂ * (y₁ + 1.0))]

function phase_portrait(𝐹::Array{Function, 1}, P::Number, title::String, xlab::String, ylab::String)
    f(x, y) = Point2f0(𝐹[1](0, x, y), 𝐹[2](0, x, y))
    xs = ys = LinRange(0, P, 3 * ceil(P))
    scene = streamplot(f, xs, ys, arrow_size=0.5, linewidth=1, colormap=:magma)  # arrow_size, colormap
    axis = scene[Axis]
    axis.names.axisnames = (xlab, ylab)
    axis.names.title = title
    axis.names.font = ("Latin Modern Math", "Latin Modern Math")
    axis.ticks.font = ("Latin Modern Math", "Latin Modern Math")
    scene.resolution = (700, 450)
    return scene
end

a = phase_portrait(array, 10, "Phase portrait for mutualism", "Population 1", "Population 2")

save("phase_plot.svg", a)  # pdf, ...

ext: jl


text:
###
# using Revise
###
using DrWatson
using GenerativeAD
using StatsBase: fit!, predict, sample
using Statistics
using BSON
using Flux
using NPZ
using Random

parameters = (
    lambda_rat 	= 1,
    lr 			= 1f-4,
    batchsize 	= 1024,
    wreg 		= 0.0,
)

# loading KDDCUP data from npz file
file = datadir("kdd_cup.npz") # taken from https://github.com/mperezcarrasco/PyTorch-DAGMM
pydata = NPZ.npzread(file)

labels = pydata["kdd"][:, end];
features = copy(pydata["kdd"][:, 1:end-1]')

normal_data = features[:, labels .== 0];
anomalous_data = features[:, labels .== 1];

# split data with fixed seed
data = GenerativeAD.Datasets.train_val_test_split(normal_data, anomalous_data; seed=1)
idim = size(data[1][1], 1)

# when running in batch randomize model seed (reset after training - deterministic random batches)
model_seed = rand(1:1000)
Random.seed!(model_seed);
encoder = Chain(Dense(idim, 60, tanh), Dense(60, 30, tanh), Dense(30, 10, tanh), Dense(10, 1))
decoder = Chain(Dense(1, 10, tanh), Dense(10, 30, tanh), Dense(30, 60, tanh), Dense(60, idim))
estimator = Chain(Dense(3, 10, tanh), Dropout(0.5), Dense(10, 4))

model = GenerativeAD.Models.DAGMM(encoder, decoder, estimator)
info, fit_t, _, _, _ = @timed fit!(model, data; patience=20, check_interval=10, parameters...)
Random.seed!()

trained_model = info.model

testmode!(trained_model, true)
_, _, z, gamma = trained_model(data[1][1])
phi, mu, cov = GenerativeAD.Models.compute_params(z, gamma)
testmode!(trained_model, false)

trn_scores = predict(trained_model, data[1][1], phi, mu, cov);
val_scores = predict(trained_model, data[2][1], phi, mu, cov);
tst_scores = predict(trained_model, data[3][1], phi, mu, cov);

using EvalMetrics: ConfusionMatrix, recall, precision, f1_score, roccurve, auc_trapezoidal

# test only
threshold = quantile(vcat(trn_scores, val_scores, tst_scores), 0.8)
# threshold = quantile(vcat(trn_scores, val_scores), 0.8)

tst_ŷ = tst_scores .> threshold;
cm = ConfusionMatrix(data[3][2], tst_ŷ)
rec, prc, f1s = recall(cm), precision(cm), f1_score(cm)

roc = roccurve(data[3][2], tst_scores)
auc = auc_trapezoidal(roc)

# test + validation
# val_ŷ = val_scores .> threshold;
# cm = ConfusionMatrix(vcat(data[2][2], data[3][2]), vcat(val_ŷ, tst_ŷ))
# recall(cm), precision(cm), f1_score(cm)

results = Dict(
    :history    => info.history,
    :niter      => info.niter,
    :model      => info.model,
    :recall     => rec,
    :precision  => prc,
    :f1         => f1s,
    :auc        => auc,
    :phi        => phi,
    :mu         => mu,
    :cov        => cov,
    :threshold  => threshold,
    :trn_scores => trn_scores,
    :trn_labels => data[1][2],
    :val_scores => val_scores,
    :val_labels => data[2][2],
    :tst_scores => tst_scores,
    :tst_labels => data[3][2]
)

@info("", model_seed, threshold, rec, prc, f1s, auc)

save(datadir("dumpster/dagmm_kddcup99_$(model_seed).bson"), results)

### analysis
using DataFrames, ValueHistories
files = readdir(datadir("dumpster"), join=true)
results = [load(f) for f in files]

df = DataFrame(
    "model_seed" => [parse(Int, split(split(f, '_')[end], '.')[1]) for f in files],
    "niter" => [r[:niter] for r in results],
    "precision" => [r[:precision] for r in results],
    "recall" => [r[:recall] for r in results],
    "f1" => [r[:f1] for r in results],
    "auc" => [r[:auc] for r in results],
    "threshold" => [r[:threshold] for r in results],
)

sort!(df, :niter)
sort!(df, :auc)

###


text:
using PyPlot

function plot_data(cls)
    figure()
    #mins, maxs, (lens, steps, grid) = MPM.get_grid(vcat(data,datab))

    #close("all")

    #ga = MPM.calc_g(grid, mymh_a.db, 20)
    #ga = exp(ga.-maximum(ga))
    #imshow(reshape(ga,lens...)', extent=[mins[1],maxs[1],mins[2],maxs[2]], aspect="equal", origin="lower")
    #colorbar()
    #contour(reshape(ga,n1,n2)', extent=gext, aspect="equal", origin="lower")
    #plot(data[:,1], data[:,2], "g.", alpha=0.8)

    #allnodes = collect(r)
    #for n in allnodes
        #if length(n.vals)!=0
            #plot(n.mins[1], n.mins[2], "g.")
        #end
    #end

    #figure()
    #gavg = ga - gb
    #imshow(reshape(gavg,n1,n2)', extent=gext, aspect="equal", origin="lower")
    #colorbar()
    #contour(reshape(gavg,n1,n2)', (0.0,), extent=gext, aspect="equal", origin="lower")
    data = cls.cls1.data
    datab = cls.cls2.data
    plot(data[:,1], data[:,2], "g.", alpha=0.8)
    plot(datab[:,1], datab[:,2], "r.", alpha=0.8)
end

function plot_traces(db, names; save=false)
    if !issubtype(typeof(names), Array)
        names = [names]
    end
    for name in names
        figure()
        if issubtype(typeof(getfield(db[1],name)), Array)
            vals = [vec(getfield(x,name)) for x in db]
        else
            vals = [getfield(x,name) for x in db]
        end
        plot(vals)
        legend(map(string,1:length(vals[1])))
        ylabel(string(name))
        title(string(name))
        xlabel("samples")
        if save
            savefig("$(string(name)).png")
        end
    end
end


text:
@testset "Day9.jl" begin
    # Write your tests here.
    lines = read_infile(Day9_1, "$(@__DIR__)/assets/day9_test1.txt")
    @test solve(Day9_1, lines, 5) == 127

    @test solve(Day9_2, lines, 5) == 62
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

using Unitful

const TEMP_REF = 20.0
const PRESSURE_REF = 1.0
const TEMP_REF_UNITFUL = TEMP_REF * u"°C"
export TEMP_REF, PRESSURE_REF

const DISPFORM_NAMES = ["Schott", "Sellmeier1", "Herzberger", "Sellmeier2", "Conrady", "Sellmeier3", "HandbookOfOptics1", "HandbookOfOptics2", "Sellmeier4", "Extended1", "Sellmeier5", "Extended2", "Extended3"]
const STATUS = ["Standard", "Preferred", "Obsolete", "Special", "Melt"]

# paths for GlassCat source file builds
const GLASSCAT_DIR = @__DIR__ # contains GlassCat.jl (pre-existing)
const AGF_DIR = joinpath(GLASSCAT_DIR, "data", "agf") # contains SCHOTT.agf, SUMITA.agf, etc.
const JL_DIR = joinpath(GLASSCAT_DIR, "data", "jl") # contains AGFGlasscat.jl, SCHOTT.jl, etc.

const SOURCES_PATH = joinpath(GLASSCAT_DIR, "data", "sources.txt")
const AGFGLASSCAT_PATH = joinpath(JL_DIR, "AGFGlassCat.jl")

# NOTE: if you change JL_DIR or AGFGLASSCAT_PATH, you also need to change the include statement in GlassCat.jl


text:
using UnPack, OrderedCollections, LinearAlgebra, JLD2, LabelledArrays, SparseArrays, RiskAdjustedLinearizations, NLsolve

if !isdefined(Main, :BansalYaron2004) || !isdefined(Main, :bansal_yaron_2004)
    include("bansal_yaron_2004.jl")
end

# BansalYaron2004 with heterogeneous risk aversion
# FIRST, we implement w/out extra forward expectations of the wealth share, to reduce the complexity.
mutable struct HetRiskAversion{T <: Real}
    p::LArray{T, 1, Array{T, 1}, (:μ_y, :ρ_x, :σ_x, :ρ_σ, :σ_y, :ς, :λ₁, :β, :ψ, :γ₁, :γ₂,
                                  :δ, :τ̅₁)}
    N_approx::LArray{Int64, 1, Array{Int64, 1}, (:q1, :q2, :ω1, :ω2)}
    S::OrderedDict{Symbol, Int}
    J::OrderedDict{Symbol, Int}
    E::OrderedDict{Symbol, Int}
    SH::OrderedDict{Symbol, Int}
end

# Parameters are based off Schorfheide et al. (2016) and converted from a monthly frequency to quarterly
function HetRiskAversion(; μ_y::T = 0.0016 * 3., ρ_x::T = 0.99^3, σ_x::T = sqrt((0.74 * sqrt(1. - ρ_x^2))^2 * 3.), ρ_σ::T = 0.99^3,
                         σ_y::T = sqrt(0.0021^2 * 3.), ς::T = sqrt(0.0014^2 * 3.), λ₁::T = .5,
                         β::T = 0.999^3, ψ::T = 2., γ₁::T = 8.9, γ₂::T = (9. - λ₁ * γ₁) / (1. - λ₁), # target 9 for average risk aversion
                         δ::T = .9, τ̅₁::T = 1e-8, N_approx::LArray{Int64, 1, Array{Int64, 1}, (:q1, :q2, :ω1, :ω2)} =
                         LVector(q1 = 1, q2 = 1, ω1 = 1, ω2 = 1)
                         ) where {T <: Real}

    @assert all([N_approx[k] > 0 for k in keys(N_approx)]) "N_approx must be at least 1 for all variables."

    # Create indexing dictionaries
    S_init  = [:yy, :x, :σ²_y, :W1, :r₋₁, :b1₋₁, :s1₋₁]            # State variables
    J_init  = [:q, :v1, :v2, :ce1, :ce2, :ω1, :ω2, :r, :c1, :c2,   # Jump variables
               :b1, :b2, :s1, :s2, :Θ1, :Θ2, :logΘ1, :Θ1_t1]
#=    J_init  = [:q, :v1, :v2, :ce1, :ce2, :ω1, :ω2, :r, :c1, :c2,   # Jump variables
               :b1, :b2, :s1, :s2, :Θ1, :Θ2, :logΘ1, :logΘ2,
               :Θ1_t1, :Θ2_t1]=#
    SH_init = [:ε_y, :ε_x, :ε_σ²]                                  # Exogenous shocks
    E_init  = [:value_fnct1, :value_fnct2,                         # Equations
               :certainty_equiv1, :certainty_equiv2,
               :ez_fwd_diff1, :ez_fwd_diff2,
               :euler1, :euler2, :cap_ret1, :cap_ret2,
               :budget_constraint1, # :budget_constraint2,
               :wealth_per_agent1, :wealth_per_agent2,
               :eq_logΘ1, # :eq_logΘ2,
               :expected_wealth_per_agent1,
               # :expected_wealth_per_agent2,
               # :consumption_mc, :bond_mc, :share_mc]
               :consumption_mc, :bond_mc, :share_mc] # bond market clears by Walras' law

    for var in [:q1, :q2, :ω1, :ω2]
        inds = (var in [:q1, :q2]) ? (1:N_approx[var]) : (0:(N_approx[var] - 1))
        push!(J_init, [Symbol(:d, var, "_t$(i)") for i in inds]...)
        push!(J_init, [Symbol(:p, var, "_t$(i)") for i in 1:N_approx[var]]...)
        push!(E_init, [Symbol(:eq_d, var, "_t$(i)") for i in inds]...)
        push!(E_init, [Symbol(:eq_p, var, "_t$(i)") for i in 1:N_approx[var]]...)
    end

    S  = OrderedDict{Symbol, Int}(k => i for (i, k) in enumerate(S_init))
    J  = OrderedDict{Symbol, Int}(k => i for (i, k) in enumerate(J_init))
    E  = OrderedDict{Symbol, Int}(k => i for (i, k) in enumerate(E_init))
    SH = OrderedDict{Symbol, Int}(k => i for (i, k) in enumerate(SH_init))

    p = @LArray [μ_y, ρ_x, σ_x, ρ_σ, σ_y, ς, λ₁, β, ψ, γ₁, γ₂, δ, τ̅₁] (:μ_y, :ρ_x, :σ_x, :ρ_σ, :σ_y, :ς, :λ₁, :β, :ψ, :γ₁, :γ₂, :δ, :τ̅₁)

    return HetRiskAversion{T}(p, N_approx, S, J, E, SH)
end

function het_risk_aversion(m::HetRiskAversion{T}; sparse_arrays::Bool = false,
                           sparse_jacobian::Vector{Symbol} = Symbol[],
                           m_rep = nothing, algorithm::Symbol = :relaxation) where {T <: Real}

    # Unpack parameters and indexing dictionaries
    @unpack p, N_approx, S, J, E, SH = m

    @unpack yy, x, σ²_y, W1, r₋₁, b1₋₁, s1₋₁ = S
    @unpack q, v1, v2, ce1, ce2, ω1, ω2, r, c1, c2, b1, b2, s1, s2 = J
    # @unpack Θ1, Θ2, logΘ1, logΘ2, Θ1_t1, Θ2_t1, logΘ1, logΘ2 = J
    @unpack Θ1, Θ2, logΘ1, Θ1_t1, logΘ1 = J
    @unpack ε_y, ε_x, ε_σ² = SH
    @unpack value_fnct1, value_fnct2, certainty_equiv1, certainty_equiv2 = E
    @unpack ez_fwd_diff1, ez_fwd_diff2, cap_ret1, cap_ret2 = E
    # @unpack euler1, euler2, budget_constraint1, budget_constraint2 = E
    @unpack euler1, euler2, budget_constraint1 = E
    # @unpack wealth_per_agent1, wealth_per_agent2, eq_logΘ1, eq_logΘ2 = E
    @unpack wealth_per_agent1, wealth_per_agent2, eq_logΘ1 = E
    # @unpack expected_wealth_per_agent1, expected_wealth_per_agent2 = E
    @unpack expected_wealth_per_agent1 = E
    @unpack consumption_mc, bond_mc, share_mc = E
    # @unpack consumption_mc, share_mc = E

    Nz = length(S)
    Ny = length(J)
    Nε = length(SH)

    ## Define nonlinear equations

    # Some helper functions
    m1_ξ(z, y) = log(p.β) + p.γ₁ * y[c1] - (p.ψ - p.γ₁) * y[ce1] - p.γ₁ * p.μ_y
    m2_ξ(z, y) = log(p.β) + p.γ₂ * y[c2] - (p.ψ - p.γ₂) * y[ce2] - p.γ₂ * p.μ_y
    function m1_fwd!(row, Γ₅, Γ₆)
        Γ₅[row, yy] = -p.γ₁
        Γ₆[row, c1] = -p.γ₁
        Γ₆[row, v1] = (p.ψ  - p.γ₁)
    end
    function m2_fwd!(row, Γ₅, Γ₆)
        Γ₅[row, yy] = -p.γ₂
        Γ₆[row, c2] = -p.γ₂
        Γ₆[row, v2]  = (p.ψ  - p.γ₂)
    end

    function μ(F, z, y) # note that y here refers to jump variables
        # Exogenous states
        F[yy]   = z[x]
        F[x]    = p.ρ_x * z[x]
        F[σ²_y] = (1. - p.ρ_σ) * p.σ_y^2 + p.ρ_σ * z[σ²_y]

        # Endogenous states
        F[W1]   = p.λ₁ * y[Θ1_t1]
        F[r₋₁]  = y[r]
        F[b1₋₁] = y[b1]
        F[s1₋₁] = y[s1]
    end

    function ξ(F, z, y)

        ## Preferences
        m1_ξv = m1_ξ(z, y) # evaluate SDF for type 1
        m2_ξv = m2_ξ(z, y) # evaluate SDF for type 2
        F[value_fnct1] = 1. / (1. - p.ψ) * (log(1. - p.β) + y[ω1]) - y[v1]
        F[value_fnct2] = 1. / (1. - p.ψ) * (log(1. - p.β) + y[ω2]) - y[v2]
        F[certainty_equiv1] = 1. / (1. - p.ψ) * (log(1. - p.β) - log(p.β) + log(exp(y[ω1]) - 1.)) - y[ce1]
        F[certainty_equiv2] = 1. / (1. - p.ψ) * (log(1. - p.β) - log(p.β) + log(exp(y[ω2]) - 1.)) - y[ce2]

        ## Euler equations
        F[euler1] = y[r] + m1_ξv
        F[euler2] = y[r] + m2_ξv

        ## Market-clearing equations and budget constraints
        Q = exp(y[q])
        C1 = exp(y[c1])
        C2 = exp(y[c2])
        B1 = -exp(y[b1]) # note that y[b1] = log(-B1) since B1 < 0
        B2 = exp(y[b2])
        S1 = exp(y[s1])
        S2 = exp(y[s2])
        λ₂ = 1. - p.λ₁
        F[consumption_mc] = log(p.λ₁ * C1 + λ₂ * C2)
        F[bond_mc]        = log(-p.λ₁ * B1) - log(λ₂ * B2)
        F[share_mc]       = log(p.λ₁ * S1 + λ₂ * S2)

        # W2 = 1. - z[W1]
        F[budget_constraint1] = log(z[W1] / p.λ₁ * (1. + Q) / (C1 + B1 + Q * S1))
        # F[budget_constraint1] = log(z[W1]) - log(p.λ₁) + log(1. + Q) - log(C1 + B1 + Q * S1)
        # F[budget_constraint2] = log(W2 / λ₂ * (1. + Q) / (C2 + B2 + Q * S2))
        # F[budget_constraint2] = log(W2) - log(λ₂) + log(1. + Q) - log(C2 + B2 + Q * S2)

        ## Wealth per agent (Θ) equations
        S1₋₁ = exp(z[s1₋₁])
        S2₋₁ = (1. - p.λ₁ * S1₋₁) / λ₂
        B1₋₁ = -exp(z[b1₋₁]) # note that z[b1₋₁] = log(-B1₋₁) since B1₋₁ < 0
        B2₋₁ = -p.λ₁ * B1₋₁ / λ₂
        R₋₁  = exp(z[r₋₁])
        τ₂   = -(p.τ̅₁ * p.λ₁ * (R₋₁ * B1₋₁ + (1. + Q) * S1₋₁)) / (λ₂ * (R₋₁ * B2₋₁ + (1. + Q) * S2₋₁))
        F[wealth_per_agent1] = log(1. - p.τ̅₁) + log(S1₋₁ + R₋₁ * B1₋₁ / (1. + Q)) - log(y[Θ1])
        # F[wealth_per_agent2] = log(R₋₁ * B2₋₁ / (1. + Q) + S2₋₁ + p.τ̅₁ * (R₋₁ * B1₋₁ / (1. + Q) + S1₋₁)) - log(y[Θ2])
        F[wealth_per_agent2] = log(1. - τ₂) + log(S2₋₁ + R₋₁ * B2₋₁ / (1. + Q)) - log(y[Θ2])
        F[eq_logΘ1]          = y[logΘ1] - log(y[Θ1])
        # F[eq_logΘ2]          = y[logΘ2] - log(y[Θ2])
        F[expected_wealth_per_agent1] = -log(y[Θ1_t1])
        # F[expected_wealth_per_agent2] = -log(y[Θ2_t1])

        ## Forward-difference equations separately handled b/c recursions
        F[cap_ret1]     = y[q] - log(sum([exp(y[J[Symbol("dq1_t$(i)")]]) for i in 1:N_approx[:q1]]) +
                                     exp(y[J[Symbol("pq1_t$(N_approx[:q1])")]]))
        F[cap_ret2]     = y[q] - log(sum([exp(y[J[Symbol("dq2_t$(i)")]]) for i in 1:N_approx[:q2]]) +
                                     exp(y[J[Symbol("pq2_t$(N_approx[:q2])")]]))
        F[ez_fwd_diff1] = y[ω1] - log(sum([exp(y[J[Symbol("dω1_t$(i)")]]) for i in 0:(N_approx[:ω1] - 1)]) +
                                     exp(y[J[Symbol("pω1_t$(N_approx[:ω1])")]]))
        F[ez_fwd_diff2] = y[ω2] - log(sum([exp(y[J[Symbol("dω2_t$(i)")]]) for i in 0:(N_approx[:ω2] - 1)]) +
                                     exp(y[J[Symbol("pω2_t$(N_approx[:ω2])")]]))

        # Set initial boundary conditions
        F[E[:eq_dq1_t1]] = p.μ_y - y[J[:dq1_t1]] + m1_ξv
        F[E[:eq_pq1_t1]] = p.μ_y - y[J[:pq1_t1]] + m1_ξv
        F[E[:eq_dq2_t1]] = p.μ_y - y[J[:dq2_t1]] + m2_ξv
        F[E[:eq_pq2_t1]] = p.μ_y - y[J[:pq2_t1]] + m2_ξv
        F[E[:eq_dω1_t0]] = y[J[:dω1_t0]]
        F[E[:eq_pω1_t1]] = p.μ_y - y[J[:pω1_t1]] + m1_ξv
        F[E[:eq_dω2_t0]] = y[J[:dω2_t0]]
        F[E[:eq_pω2_t1]] = p.μ_y - y[J[:pω2_t1]] + m2_ξv

        # Recursions for forward-difference equations
        for i in 2:N_approx[:q1]
            F[E[Symbol("eq_dq1_t$(i)")]] = p.μ_y - y[J[Symbol("dq1_t$(i)")]] + m1_ξv
            F[E[Symbol("eq_pq1_t$(i)")]] = p.μ_y - y[J[Symbol("pq1_t$(i)")]] + m1_ξv
        end
        for i in 2:N_approx[:q2]
            F[E[Symbol("eq_dq2_t$(i)")]] = p.μ_y - y[J[Symbol("dq2_t$(i)")]] + m2_ξv
            F[E[Symbol("eq_pq2_t$(i)")]] = p.μ_y - y[J[Symbol("pq2_t$(i)")]] + m2_ξv
        end
        for i in 2:N_approx[:ω1]
            F[E[Symbol("eq_dω1_t$(i-1)")]] = p.μ_y - y[J[Symbol("dω1_t$(i-1)")]] + m1_ξv
            F[E[Symbol("eq_pω1_t$(i)")]]   = p.μ_y - y[J[Symbol("pω1_t$(i)")]]   + m1_ξv
        end
        for i in 2:N_approx[:ω2]
            F[E[Symbol("eq_dω2_t$(i-1)")]] = p.μ_y - y[J[Symbol("dω2_t$(i-1)")]] + m2_ξv
            F[E[Symbol("eq_pω2_t$(i)")]]   = p.μ_y - y[J[Symbol("pω2_t$(i)")]]   + m2_ξv
        end
    end

    # The cache is initialized as zeros so we only need to fill non-zero elements
    function Λ(F, z)
        # Heteroskedastic risk arises b/c the wealth share is a state variable
        F[W1, Θ1] = p.λ₁
    end

    # The cache is initialized as zeros so we only need to fill non-zero elements
    function Σ(F, z)
        F[yy,   ε_y]  = sqrt(z[σ²_y])         # Take square root b/c Σ is not a variance-covariance matrix.
        F[x,    ε_x]  = sqrt(z[σ²_y]) * p.σ_x # It is instead the "volatility" loading on the martingale difference sequences.
        F[σ²_y, ε_σ²] = sqrt(z[σ²_y]) * p.ς
    end

    function ccgf(F, A, z)
        # F .= .5 * diag(α * α') # slower but this is the underlying math
        sum!(F, A.^2) # faster implementation
        F .*= .5
    end

    if sparse_arrays
        Γ₅ = spzeros(T, Ny, Nz)
        Γ₆ = spzeros(T, Ny, Ny)
    else
        Γ₅ = zeros(T, Ny, Nz)
        Γ₆ = zeros(T, Ny, Ny)
    end

    # Equations other than forward difference equations
    m1_fwd!(euler1, Γ₅, Γ₆)
    m2_fwd!(euler2, Γ₅, Γ₆)
    Γ₆[expected_wealth_per_agent1, logΘ1] = one(T)

    # Forward difference equations: boundary conditions
    m1_fwd!(E[:eq_dq1_t1], Γ₅, Γ₆)
    Γ₅[E[:eq_dq1_t1], yy] = one(T)

    m1_fwd!(E[:eq_pq1_t1], Γ₅, Γ₆)
    Γ₅[E[:eq_pq1_t1], yy] = one(T)
    Γ₆[E[:eq_pq1_t1], q]  = one(T)

    m1_fwd!(E[:eq_pω1_t1], Γ₅, Γ₆)
    Γ₅[E[:eq_pω1_t1], yy] = one(T)
    Γ₆[E[:eq_pω1_t1], c1] = one(T)
    Γ₆[E[:eq_pω1_t1], ω1] = one(T)

    m2_fwd!(E[:eq_dq2_t1], Γ₅, Γ₆)
    Γ₅[E[:eq_dq2_t1], yy] = one(T)

    m2_fwd!(E[:eq_pq2_t1], Γ₅, Γ₆)
    Γ₅[E[:eq_pq2_t1], yy] = one(T)
    Γ₆[E[:eq_pq2_t1], q]  = one(T)

    m2_fwd!(E[:eq_pω2_t1], Γ₅, Γ₆)
    Γ₅[E[:eq_pω2_t1], yy] = one(T)
    Γ₆[E[:eq_pω2_t1], c2] = one(T)
    Γ₆[E[:eq_pω2_t1], ω2] = one(T)

    # Forward difference equations: recursions
    for i in 2:N_approx[:q1]
        m1_fwd!(E[Symbol("eq_dq1_t$(i)")], Γ₅, Γ₆)
        Γ₅[E[Symbol("eq_dq1_t$(i)")], yy] = one(T)
        Γ₆[E[Symbol("eq_dq1_t$(i)")], J[Symbol("dq1_t$(i-1)")]] = one(T)

        m1_fwd!(E[Symbol("eq_pq1_t$(i)")], Γ₅, Γ₆)
        Γ₅[E[Symbol("eq_pq1_t$(i)")], yy] = one(T)
        Γ₆[E[Symbol("eq_pq1_t$(i)")], J[Symbol("pq1_t$(i-1)")]] = one(T)
    end
    for i in 2:N_approx[:ω1]
        m1_fwd!(E[Symbol("eq_dω1_t$(i-1)")], Γ₅, Γ₆)
        Γ₅[E[Symbol("eq_dω1_t$(i-1)")], yy] = one(T)
        Γ₆[E[Symbol("eq_dω1_t$(i-1)")], c1] = one(T)
        Γ₆[E[Symbol("eq_dω1_t$(i-1)")], J[Symbol("dω1_t$(i-2)")]] = one(T)

        m1_fwd!(E[Symbol("eq_pω1_t$(i)")], Γ₅, Γ₆)
        Γ₅[E[Symbol("eq_pω1_t$(i)")], yy] = one(T)
        Γ₆[E[Symbol("eq_pω1_t$(i)")], c1] = one(T)
        Γ₆[E[Symbol("eq_pω1_t$(i)")], J[Symbol("pω1_t$(i-1)")]] = one(T)
    end
    for i in 2:N_approx[:q2]
        m2_fwd!(E[Symbol("eq_dq2_t$(i)")], Γ₅, Γ₆)
        Γ₅[E[Symbol("eq_dq2_t$(i)")], yy] = one(T)
        Γ₆[E[Symbol("eq_dq2_t$(i)")], J[Symbol("dq2_t$(i-1)")]] = one(T)

        m2_fwd!(E[Symbol("eq_pq2_t$(i)")], Γ₅, Γ₆)
        Γ₅[E[Symbol("eq_pq2_t$(i)")], yy] = one(T)
        Γ₆[E[Symbol("eq_pq2_t$(i)")], J[Symbol("pq2_t$(i-1)")]] = one(T)
    end
    for i in 2:N_approx[:ω2]
        m2_fwd!(E[Symbol("eq_dω2_t$(i-1)")], Γ₅, Γ₆)
        Γ₅[E[Symbol("eq_dω2_t$(i-1)")], yy] = one(T)
        Γ₆[E[Symbol("eq_dω2_t$(i-1)")], c2] = one(T)
        Γ₆[E[Symbol("eq_dω2_t$(i-1)")], J[Symbol("dω2_t$(i-2)")]] = one(T)

        m2_fwd!(E[Symbol("eq_pω2_t$(i)")], Γ₅, Γ₆)
        Γ₅[E[Symbol("eq_pω2_t$(i)")], yy] = one(T)
        Γ₆[E[Symbol("eq_pω2_t$(i)")], c2] = one(T)
        Γ₆[E[Symbol("eq_pω2_t$(i)")], J[Symbol("pω2_t$(i-1)")]] = one(T)
    end

    z, y, Ψ = create_guess(m, (isnothing(m_rep) ?
                               BansalYaron2004(μ_y = p.μ_y, ρ_x = p.ρ_x, σ_x = p.σ_x,
                                               ρ_σ = p.ρ_σ, σ_y = p.σ_y, ς = p.ς,
                                               β = p.β, ψ = p.ψ, γ = p.λ₁ * p.γ₁ + (1. - p.λ₁) * p.γ₂,
                                               N_approx = LVector(q = N_approx[:q1], ω = N_approx[:ω1])) : m_rep); algorithm = algorithm)

    if sparse_arrays
        return RiskAdjustedLinearization(μ, Λ, Σ, ξ, Γ₅, Γ₆, ccgf, vec(z), vec(y), Ψ, Nε; sparse_jacobian = sparse_jacobian,
                                         Λ_cache_init = dims -> spzeros(dims...),
                                         Σ_cache_init = dims -> Diagonal(Vector{eltype(z)}(undef, dims[1])))
    else
        return RiskAdjustedLinearization(μ, Λ, Σ, ξ, Γ₅, Γ₆, ccgf, vec(z), vec(y), Ψ, Nε; sparse_jacobian = sparse_jacobian)
    end
end

function create_guess(m::HetRiskAversion{T}, m_rep::BansalYaron2004; algorithm::Symbol = :relaxation) where {T <: Real}

    ## Set up

    # Get guesses from BansalYaron2004 representative agent model
    ral_rep = bansal_yaron_2004(m_rep)
    solve!(ral_rep, algorithm = algorithm, verbose = :none)
    yrep = ral_rep.y
    zrep = ral_rep.z
    Ψrep = ral_rep.Ψ
    Jr   = m_rep.J
    Sr   = m_rep.S

    # Unpack parameters and indexing dictionaries
    @unpack p, N_approx, S, J, E, SH = m
    @unpack yy, x, σ²_y, W1, r₋₁, b1₋₁, s1₋₁ = S
    @unpack q, v1, v2, ce1, ce2, ω1, ω2, r = J
    @unpack c1, c2, b1, b2, s1, s2 = J
    # @unpack Θ1, Θ2, logΘ1, logΘ2, Θ1_t1, Θ2_t1, logΘ1, logΘ2 = J
    @unpack Θ1, Θ2, logΘ1, Θ1_t1, logΘ1 = J

    # Initialize deterministic steady state guesses
    Nz = length(S)
    Ny = length(J)
    z = Vector{T}(undef, Nz)
    y = Vector{T}(undef, Ny)
    Ψ = zeros(T, Ny, Nz)

    ## Compute guesses for z, y

    # Steady state values of state variables known ex-ante
    z[yy]   = 0.
    z[x]    = 0.
    z[σ²_y] = p.σ_y^2

    # Guesses based on representative agent model's solution
    y[ω1]  = yrep[Jr[:ω]]
    y[ω2]  = yrep[Jr[:ω]]
    y[q]   = yrep[Jr[:q]]
    y[v1]  = yrep[Jr[:v]]
    y[ce1] = yrep[Jr[:ce]]
    y[v2]  = yrep[Jr[:v]]
    y[ce2] = yrep[Jr[:ce]]

    # Guesses for consumption and portfolio choice
    S1 = (1. / p.γ₁) / (1. / p.γ₁ + 1. / p.γ₂) / p.λ₁
    S2 = (1. - p.λ₁ * S1) / (1. - p.λ₁)
    C1 = .99 * S1 # .99 just for now b/c heterogeneity is small
    C2 = (1. - p.λ₁ * C1) / (1. - p.λ₁)
    M0 = m_rep.p.β * (yrep[Jr[:v]] / yrep[Jr[:ce]])^(m_rep.p.ψ - m_rep.p.γ) *
        exp(-m_rep.p.γ * m_rep.p.μ_y) # make a guess for SDF
    R = 1. / M0 # just guess it's 1 / SDF
    Q = exp(yrep[Jr[:q]])
    B1 = (S1 * ((1. - p.τ̅₁) - p.τ̅₁ * Q) - C1) / (1. - (1. - p.τ̅₁) * R)
    B2 = -p.λ₁ * B1 / (1. - p.λ₁)
    τ₂ = -(p.τ̅₁ * (R * B1 + (1 + Q) * S1)) / (R * B2 + (1 + Q) * S2)
    Θ10 = (1. - p.τ̅₁) * (S1 + R * B1 / (1. + Q))
    Θ20 = (1. - τ₂) * (S2 + R * B2 / (1. + Q))

    y[s1] = log(S1)
    y[s2] = log(S2)
    y[c1] = log(C1)
    y[c2] = log(C2)
    y[b1] = log(-B1) # agent 1 leveraged -> B10 < 0
    y[b2] = log(B2)
    y[Θ1] = Θ10
    y[Θ2] = Θ20
    y[logΘ1] = log(Θ10)
    # y[logΘ2] = log(Θ20)
    y[Θ1_t1] = Θ10 # in a steady state, the one-period ahead expectation agrees with the current day value.
    # y[Θ2_t1] = Θ20
    y[r] = log(R)
    z[W1] = p.λ₁ * Θ10
    z[r₋₁] = y[r]
    z[b1₋₁] = y[b1]
    z[s1₋₁] = y[s1]

    M10 = m_rep.p.β * (yrep[Jr[:v]] / yrep[Jr[:ce]])^(m_rep.p.ψ - p.γ₁) * # use agent 1's risk aversion here
        exp(-m_rep.p.γ * m_rep.p.μ_y) # make a guess for SDF
    M20 = m_rep.p.β * (yrep[Jr[:v]] / yrep[Jr[:ce]])^(m_rep.p.ψ - p.γ₂) * # use agent 1's risk aversion here
        exp(-m_rep.p.γ * m_rep.p.μ_y) # make a guess for SDF
    Y0 = 1. # steady state growth rate in endowment
    y[J[:dq1_t1]] = convert(T, log(exp(p.μ_y) * M10 * Y0))
    y[J[:pq1_t1]] = convert(T, log(exp(p.μ_y) * M10 * Y0 * Q))
    y[J[:dω1_t0]] = zero(T)
    y[J[:pω1_t1]] = convert(T, log(exp(p.μ_y) * M10 * Y0 * exp(y[ω1])))
    y[J[:dq2_t1]] = convert(T, log(exp(p.μ_y) * M20 * Y0))
    y[J[:pq2_t1]] = convert(T, log(exp(p.μ_y) * M20 * Y0 * Q))
    y[J[:dω2_t0]] = zero(T)
    y[J[:pω2_t1]] = convert(T, log(exp(p.μ_y) * M20 * Y0 * exp(y[ω2])))

    for i in 2:N_approx[:q1]
        y[J[Symbol("dq1_t$(i)")]] = convert(T, log(M10) + μ_y + log(Y0) + y[J[Symbol("dq1_t$(i-1)")]])
        y[J[Symbol("pq1_t$(i)")]] = convert(T, log(M10) + μ_y + log(Y0) + y[J[Symbol("pq1_t$(i-1)")]])
    end
    for i in 2:N_approx[:q2]
        y[J[Symbol("dq2_t$(i)")]] = convert(T, log(M20) + μ_y + log(Y0) + y[J[Symbol("dq2_t$(i-1)")]])
        y[J[Symbol("pq2_t$(i)")]] = convert(T, log(M20) + μ_y + log(Y0) + y[J[Symbol("pq2_t$(i-1)")]])
    end
    for i in 2:N_approx[:ω1]
        y[J[Symbol("dω1_t$(i-1)")]] = convert(T, μ_y + log(M10) + log(Y0) + y[J[Symbol("dω1_t$(i-2)")]])
        y[J[Symbol("pω1_t$(i)")]] = convert(T, μ_y + log(M10) + log(Y0) + y[J[Symbol("pω1_t$(i-1)")]])
    end
    for i in 2:N_approx[:ω2]
        y[J[Symbol("dω2_t$(i-1)")]] = convert(T, μ_y + log(M20) + log(Y0) + y[J[Symbol("dω2_t$(i-2)")]])
        y[J[Symbol("pω2_t$(i)")]] = convert(T, μ_y + log(M20) + log(Y0) + y[J[Symbol("pω2_t$(i-1)")]])
    end

    ## Populate Ψ

    # Parts inferred from representative agent model
    Ψ[q, 1:3] = Ψrep[Jr[:q], 1:3] # indexes 1-3 are the common states for yy, x, and σ²_y
    Ψ[v1, 1:3] = Ψrep[Jr[:v], 1:3]
    Ψ[v2, 1:3] = Ψrep[Jr[:v], 1:3]
    Ψ[ce1, 1:3] = Ψrep[Jr[:ce], 1:3]
    Ψ[ce2, 1:3] = Ψrep[Jr[:ce], 1:3]
    Ψ[ω1, 1:3] = Ψrep[Jr[:ω], 1:3]
    Ψ[ω2, 1:3] = Ψrep[Jr[:ω], 1:3]

    for i in 2:N_approx[:q1]
        Ψ[J[Symbol("dq1_t$(i)")], 1:3] = Ψrep[Jr[Symbol("dq$(i)")], 1:3]
        Ψ[J[Symbol("pq1_t$(i)")], 1:3] = Ψrep[Jr[Symbol("pq$(i)")], 1:3]
    end
    for i in 2:N_approx[:q2]
        Ψ[J[Symbol("dq2_t$(i)")], 1:3] = Ψrep[Jr[Symbol("dq$(i)")], 1:3]
        Ψ[J[Symbol("pq2_t$(i)")], 1:3] = Ψrep[Jr[Symbol("pq$(i)")], 1:3]
    end
    for i in 2:N_approx[:ω1]
        Ψ[J[Symbol("dω1_t$(i-1)")], 1:3] = Ψrep[Jr[Symbol("dω$(i-1)")], 1:3]
        Ψ[J[Symbol("pω1_t$(i)")], 1:3] = Ψrep[Jr[Symbol("pω$(i)")], 1:3]
    end
    for i in 2:N_approx[:ω2]
        Ψ[J[Symbol("dω2_t$(i-1)")], 1:3] = Ψrep[Jr[Symbol("dω$(i-1)")], 1:3]
        Ψ[J[Symbol("pω2_t$(i)")], 1:3] = Ψrep[Jr[Symbol("pω$(i)")], 1:3]
    end

    # Guess dependence on the wealth share

    return z, y, Ψ
end


text:
const __m128i = NTuple{2, VecElement{UInt64}}
@assert ccall(
    "llvm.x86.aesni.aeskeygenassist", llvmcall, __m128i, (__m128i, UInt8),
    __m128i((0x0123456789123450, 0x9876543210987654)), 0x1
) ≡ __m128i((0x857c266f7c266e85, 0x2346382146382023))


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Testing

#=
test subjects:
- surface pressure load in curved surface
- verification of elements wedge6 and wedge15

from Code Aster:
N38      -8.85861895037377E-01 -3.46944695195361E-18 -3.46944695195361E-18
coords of N38 = (1.0, 0.0, 0.0)

Analytical solution
uᵣ(r) = b³p/(2Er²(a³-b³)) * (a³(ν+1) + r³(-4ν+2)), where
a = inner surface radial distance, b = outer surface ...

if a=0.9, b=1.0, ν=1/3, E = 24580 and p = 7317 equation yields
-9/10 for radial displacement
=#

@testset """1/8 hollow sphere with surface load""" begin
    mesh_file = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
    mesh = aster_read_mesh(mesh_file, "HOLLOWSPHERE8_WEDGE6")
    body = Problem(Elasticity, "hollow sphere 1/8 model", 3)
    body.elements = create_elements(mesh, "HOLLOWSPHERE8")
    update!(body, "youngs modulus", 24580.0)
    update!(body, "poissons ratio", 1/3)
    bc = Problem(Dirichlet, "symmetry bc", 3, "displacement")
    el1 = create_elements(mesh, "FACE1")
    update!(el1, "displacement 3", 0.0)
    el2 = create_elements(mesh, "FACE2")
    update!(el2, "displacement 2", 0.0)
    el3 = create_elements(mesh, "FACE3")
    update!(el3, "displacement 1", 0.0)
    bc.elements = [el1; el2; el3]
    lo = Problem(Elasticity, "pressure load", 3)
    lo.elements = create_elements(mesh, "OUTER")
    update!(lo, "surface pressure", 7317.0)
    solver = LinearSolver(body, bc, lo)
    solver()

    X = lo("geometry")
    u = lo("displacement")
    nids = sort(collect(keys(X)))
    umag = Float64[norm(u[id]) for id in nids]
    um = mean(umag)
    info("mean umag = $um")
    info("std umag = ", std(umag))
    rtol = norm(um - 0.9) / max(norm(um), 0.9) * 100.0
    info("rtol = $rtol")
    @test rtol < 1.5 # percents
    u_CA = [-8.85861895037377E-01, -3.46944695195361E-18, -3.46944695195361E-18]
    rtol = norm(u[38] - u_CA) / max(norm(u[38]), norm(u_CA)) * 100.0
    info("rel diff to CA = $rtol %")
    @test isapprox(u[38], u_CA)
end


text:
# This file contains essentials functions to define
# tables to geometric, inner and outer products.

"""
```
bbgeometric(a::BasisBlade, b::BasisBlade) (Liga function)
```
This function is used to determine the geometric product
of basis blade. Essentially, this function is used when layout
function is trigged, but it can be used separated too.
## Example
```julia-repl
julia> bbgeometric(e1,e12)
```
returns e2.
"""
function bbgeometric(a::BasisBlade, b::BasisBlade)
		len = GIdLayout[1]+GIdLayout[2] #p+q
		Ea = copy(LogicalBasisBlade[a.index])
		Eb = copy(LogicalBasisBlade[b.index])
		abscalar = 1
		abindex = 1
		ablogical = Vector{Bool}(undef,len)
		cont = 0
		for i=1:len-1
			if Eb[i] == true
				for j=i+1:len
					if Ea[j] == true
						cont = cont + 1
					end
				end
			end
		end
		for i=1:len
			if Ea[i] != Eb[i]
				ablogical[i] = true
			else
					ablogical[i] = false
					if Ea[i] == true && i>GIdLayout[1]
						cont =cont +1
					end
			end
		end
		abscalar = ((-1)^cont)
		ind=1
		while ablogical != LogicalBasisBlade[ind]
			ind = ind+1
		end
		abindex = ind
		ab=BasisBlade(abscalar*abindex)
		return ab
end

"""
```
bbinner(a::BasisBlade,b::BasisBlade) (Liga function)
```
This function is used to determine the inner product
of basis blades. Essentially, this function is used when layout
function is trigged, but it can be used separated too.
## Example
```julia-repl
julia> bbinner(e12,e2)
```
returns the Basis Blade e1.
"""
function bbinner(a::BasisBlade,b::BasisBlade)
		k=grade(a)
		l=grade(b)
		if k != 0 && l !=0
			return gradeprojection(bbgeometric(a,b),abs(k-l))
		else
			return BasisBlade(0)
		end
end

"""
```
bbouter(a::BasisBlade,b::BasisBlade) (Liga function)
```
This function is used to determine the outer product
of basis blades. Essentially, this function is used when layout
function is trigged, but it can be used separated too.
## Example
```julia-repl
julia> bbouter(e12,e3)
```
returns the Basis Blade e123.
"""
function bbouter(a::BasisBlade,b::BasisBlade)
	k=grade(a)
	l=grade(b)
	return gradeprojection(bbgeometric(a,b),k+l)
end

"""
```
OperationTable() (Liga function)
```
OperationTable is function trigged by layout function and
it is used to define the operational table related to geometric,
inner and outer product. As result, constants array are created
to optimize these products to multivectors.
## Example
```julia-repl
julia> A,B,C,D=OperationTable()
```
returns A a table for geometric product, B a table for inner product,
C a table for outer product and D a tensor.
"""
function OperationTable()
	len=2^(GIdLayout[1]+GIdLayout[2])
	OT=Array{Any,2}(undef,len,len)
	IOT=Array{Any,2}(undef,len,len)
	OOT=Array{Any,2}(undef,len,len)
	Tensor=Vector{SparseMatrixCSC}(undef,len)
	for i=1:len
        for j=1:len
			OT[i,j]=bbgeometric(BasisBlade(i),BasisBlade(j))
			IOT[i,j]=bbinner(BasisBlade(i),BasisBlade(j))
			OOT[i,j]=bbouter(BasisBlade(i),BasisBlade(j))
		end
	end
	for k=1:len
		Tensor[k]=spzeros(len,len)
		for i=1:len
			for j=1:len
				if OT[i,j]==BasisBlade(k)
					Tensor[k][i,j]=1.0
				end
				if OT[i,j]==-BasisBlade(k)
				   Tensor[k][i,j]=-1.0
			    end
			end
		end
	end
	#display(Tensor)
	return OT,IOT,OOT,Tensor
end


text:
function build_globalmin{T <: Union(ASCIIString, UTF8String)}(list::Array{T})

  returns  = returns(list)
  params  = modelparams(returns)
  weights = rowsums(Array(inv(params[5]))) / sum(inv(params[5]))

  port = Portfolio(params[1],
                   params[2],
                   params[3],
                   params[4],
                   params[5],
                   weights)
  port
end

function build_efficient(x)
  #code here
end

function build_tangency(x)
  #code here
end


text:
# logical operators
n = 3
println(n % 2 == 0 || n % 3 == 0)
println(n % 2 == 0 && n % 3 == 0)
println(!(n % 2 == 0) && n % 3 == 0)
println(!true)
println("This is valid: 0 < n < 10 = $(0 < n < 10)")

# Conditional operators
if n > 0
    println("n is greater than 0")
end

a, b = 3, 4
if a % 2 == 0
    println("a is even")
elseif b == 4
    println("b is 4, die!")
elseif a < b
    println("a is smaller than b")
else
    println("a is greater than b")
end

# Nested conditional
if a < 10
    if a % 2 != 0
        println("a is odd and less than 10")
    end
end

if a < 10 && b > 3
    println("a is less than 10 and b is greater than 3")
end


text:
module Ogg
using Ogg_jll
using FileIO
import Base: show, convert
export load, save, OggDecoder, OggEncoder

include("types.jl")
include("decoder.jl")
include("encoder.jl")

end # module


text:
"""
    transform(se::SummarizedExperiment, assay_name::String, method::Function, kwargs...)

Applies a transformation to a given assay and returns it as a `Matrix{<:Real}` object.

# Arguments
- `se::SummarizedExperiment`: the experiment object of interest.
- `assay_name::String`: the name of the assay to base the evaluation on.
- `method::Function`: the transformation to apply. It can entail the following values:
  `log10`, `pa`, `ztransform`, `clr`.
- `kwargs...`: additional parameters specific to some methods.
"""
function transform(se::SummarizedExperiment, assay_name::String, method::Function, kwargs...)

    method(assay(se, assay_name), kwargs...)

end

"""
    transform(se::SummarizedExperiment, assay_name::String, method::Function, kwargs...; output_name::String = "transformed_" * assay_name)

Applies a transformation to a given assay and stores it into `se` as a new assay.

# Arguments
- `se::SummarizedExperiment`: the experiment object of interest.
- `assay_name::String`: the name of the assay to base the evaluation on.
- `method::Function`: the transformation to apply. It can entail the following values:
  `log10`, `pa`, `ztransform`, `clr`.
- `kwargs...`: additional parameters specific to some methods.
- `output_name::String`: custom name for the newly created assay.
"""
function transform!(se::SummarizedExperiment, assay_name::String, method::Function, kwargs...; output_name::String = "transformed_" * assay_name)

    se.assays[output_name] = method(assay(se, assay_name), kwargs...)

    return se

end

function log10(assay::Matrix, pseudocount::Real = 0)

    assay .+= pseudocount

    return log10.(assay)

end

pa(assay::Matrix, threshold::Real) = map(x -> ifelse(x > threshold, 1, 0), assay)

function relabund(assay::Matrix)

    mat = zeros(size(assay))

    for (idx, col) in enumerate(eachcol(assay))

        tot = sum(col)
        mat[:, idx] = map(x -> x / tot, col)

    end

    return mat

end

function ztransform(assay::Matrix)

    mat = zeros(size(assay))

    for (idx, col) in enumerate(eachcol(assay))

        avg = mean(col)
        stdev = std(col)

        mat[:, idx] = map(x -> (x - avg) / stdev, col)

    end

    return mat

end

function clr(assay::Matrix, pseudocount::Real = 0)

    assay .+= pseudocount

    clog = log.(assay)
    mat = zeros(size(assay))

    for (idx, col) in enumerate(eachcol(clog))

        log_mean = mean(col)

        mat[:, idx] = map(x -> x - log_mean, col)

    end

    return mat

end


text:
using ApproxFun, SpecialFunctions, LinearAlgebra

include("../ODEs/Airy equation.jl")
include("../ODEs/Blasius and Falkner-Skan.jl")
include("../ODEs/Lee & Greengard equation.jl")
include("../ODEs/Multivariate Newton iteration.jl")
include("../ODEs/Newton iteration for nonlinear BVPs.jl")
include("../ODEs/Picard iteration.jl")
include("../ODEs/Quantum states.jl")
include("../ODEs/Tenth order ODE.jl")

include("../PDEs/Linear KdV.jl")

include("../Sampling/2D Cauchy Distribution.jl")
include("../Sampling/GUE Sampling.jl")

include("../Extras/Delay Equation.jl")
include("../Extras/Dual Numbers.jl")
include("../Extras/Lanczos.jl")
include("../Extras/Maual 3D interpolation.jl")


text:
module DataValueInterfaces

"""
    nondatavaluetype(T)

For a type `T`, return the corresponding non-`DataValue` type, translating between
`Union{T, Missing}` and `DataValue{T}`.

For example, `nondatavaluetype(Int64)` returns `Int64`, while
`nondatavaluetype(DataValue{Int64})` returns `Union{Int64, Missing}`.

This generic function is owned by DataValueInterfaces.jl itself, which is the sole provider of the
default definition.
"""
function nondatavaluetype end

nondatavaluetype(::Type{T}) where {T} = T

"""
    datavaluetype(T)

For a type `T`, return the corresponding `DataValue` type, translating between
`Union{T, Missing}` and `DataValue{T}`.

For example, `datavaluetype(Int64)` returns `Int64`, while
`datavaluetype(Union{Int64, Missing})` returns `DataValue{Int64}`.

This generic function is owned by DataValueInterfaces.jl itself, which is the sole provider of the
default definition.
"""
function datavaluetype end

datavaluetype(::Type{T}) where {T} = T

end # module


text:
# Build up functions
"""
    function filter_kwargs(; kwargs...)

Gets and filters keyword arguments to build the example system for case 4
"""
function filter_kwargs(; kwargs...)
    system_kwargs = filter(x -> in(first(x), PSY.SYSTEM_KWARGS), kwargs)
    return (system_kwargs)
end

# Buses of the System
"""
    function nodes4_tcrd()

Returns the data with the information of the buses of the example system case 4
"""
function nodes4_tcrd()
    nodes4 = [
        Bus(
            1,
            "Bus 1",
            "REF",
            0.0,
            1.0,
            (min = 0.9, max = 1.1),
            220,
            nothing,
            nothing),
        Bus(
            2,
            "Bus 2",
            "PV",
            0.0,
            1.0,
            (min = 0.9, max = 1.1),
            220,
            nothing,
            nothing,
        ),
        Bus(
            3,
            "Bus 3",
            "PV",
            0.0,
            1.0,
            (min = 0.9, max = 1.1),
            220,
            nothing,
            nothing,
        ),
        Bus(
            4,
            "Bus 4",
            "PV",
            0.0,
            1.0,
            (min = 0.9, max = 1.1),
            220,
            nothing,
            nothing,
        ),
    ]
    return nodes4
end

# Lines of the System
"""
    function branches4_tcrd(nodes4_tcrd)

Returns the data with the information of the Lines of the system for the example case 4

# Arguments
- `nodes4_tcrd`:                            Array with the information of the buses of
                                            the example system case 4
"""
function branches4_tcrd(nodes4_tcrd)
    branches4 = [
        Line(
            "Line1",
            true,
            0.0,
            0.0,
            Arc(from = nodes4_tcrd[1], to = nodes4_tcrd[2]),
            0.0,
            0.1,
            (from = 0.0, to = 0.0),
            0.20,
            1.0,
        ),
        Line(
            "Line2",
            true,
            0.0,
            0.0,
            Arc(from = nodes4_tcrd[1], to = nodes4_tcrd[3]),
            0.0,
            0.1,
            (from = 0.0, to = 0.0),
            2.0,
            1.0,
        ),
        Line(
            "Line3",
            true,
            0.0,
            0.0,
            Arc(from = nodes4_tcrd[2], to = nodes4_tcrd[3]),
            0.0,
            0.1,
            (from = 0.0, to = 0.0),
            2.0,
            1.0,
        ),
        Line(
            "Line4",
            true,
            0.0,
            0.0,
            Arc(from = nodes4_tcrd[2], to = nodes4_tcrd[4]),
            0.0,
            0.1,
            (from = 0.0, to = 0.0),
            2.0,
            1.0,
        ),
        Line(
            "Line5",
            true,
            0.0,
            0.0,
            Arc(from = nodes4_tcrd[3], to = nodes4_tcrd[4]),
            0.0,
            0.1,
            (from = 0.0, to = 0.0),
            2.0,
            1.0,
        ),
    ]
    return branches4
end

# Thermal Generators of the System
"""
    function thermal_generators4_tcrd(nodes4_tcrd)

Returns the data with the information of the thermal generators of the system for the
example case 4

# Arguments
- `nodes4_tcrd`:                            Array with the information of the buses of
                                            the example system case 4
"""
function thermal_generators4_tcrd(nodes4_tcrd)
    thermal4 = [
        ThermalStandard(
            name = "GBus1",
            available = true,
            status = true,
            bus = nodes4_tcrd[1],
            active_power = 0.50,
            reactive_power = 0.10,
            rating = 0.500,
            prime_mover = PrimeMovers.ST,
            fuel = ThermalFuels.COAL,
            active_power_limits = (min = 0.0, max = 1.0),
            reactive_power_limits = (min = -0.25, max = 0.25),
            time_limits = nothing,
            ramp_limits = nothing,
            operation_cost = ThreePartCost((0.175000, 10.0000), 0.0, 0.0, 0.0),
            base_power = 100.0,
        ),
        ThermalStandard(
            name = "GBus2",
            available = true,
            status = true,
            bus = nodes4_tcrd[2],
            active_power = 0.50,
            reactive_power = 0.10,
            rating = 0.500,
            prime_mover = PrimeMovers.ST,
            fuel = ThermalFuels.COAL,
            active_power_limits = (min = 0.0, max = 0.50),
            reactive_power_limits = (min = -0.25, max = 0.25),
            time_limits = nothing,
            ramp_limits = nothing,
            operation_cost = ThreePartCost((0.497000, 10.0000), 0.0, 0.0, 0.0),
            base_power = 100.0,
        ),
        ThermalStandard(
            name = "GBus3",
            available = true,
            status = true,
            bus = nodes4_tcrd[3],
            active_power = 0.50,
            reactive_power = 0.10,
            rating = 2.00,
            prime_mover = PrimeMovers.ST,
            fuel = ThermalFuels.COAL,
            active_power_limits = (min = 0.0, max = 2.00),
            reactive_power_limits = (min = -1.00, max = 1.00),
            time_limits = nothing,
            ramp_limits = nothing,
            operation_cost = ThreePartCost((0.260000, 20.000), 0.0, 0.0, 0.0),
            base_power = 100.0,
        ),
        ThermalStandard(
            name = "GBus4",
            available = true,
            status = true,
            bus = nodes4_tcrd[4],
            active_power = 0.50,
            reactive_power = 0.10,
            rating = 2.00,
            prime_mover = PrimeMovers.ST,
            fuel = ThermalFuels.COAL,
            active_power_limits = (min = 0.0, max = 2.00),
            reactive_power_limits = (min = -1.00, max = 1.00),
            time_limits = nothing,
            ramp_limits = nothing,
            operation_cost = ThreePartCost((0.325000, 20.000), 0.0, 0.0, 0.0),
            base_power = 100.0,
        ),
    ]
    return thermal4
end

# Loads of the System
"""
    function loads4_tcrd(nodes4_tcrd)

Returns the data with the information of the Loads of the system for the example case 4

# Arguments
- `nodes4_tcrd`:                            Array with the information of the buses of
                                            the example system case 4
"""
function loads4_tcrd(nodes4_tcrd)
    loads4 = [
        PowerLoad(
            "Bus3",
            true,
            nodes4_tcrd[3],
            LoadModels.ConstantPower,
            1.0,
            0.0,
            100.0,
            1.0,
            0.0,
        ),
        PowerLoad(
            "Bus4",
            true,
            nodes4_tcrd[4],
            LoadModels.ConstantPower,
            1.0,
            0.0,
            100.0,
            1.0,
            0.0,
        ),
    ]
    return loads4
end

max_stars_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_issues_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_forks_repo_name: invenia/MaximumGeneratorProfitTCRDD


text:
function test_blob_reshape(backend::Backend, T, eps)
  println("-- Testing blob reshape on $(typeof(backend)){$T}...")

  data = rand(T, 2,4,5,6)
  blob = make_blob(backend, data)
  blob2 = reshape_blob(backend, blob, length(data),1,1,1)
  data2 = rand(T, size(blob2))
  copy!(blob2, data2)
  copy!(data, blob)
  @test all(abs.(data - reshape(data2,size(data))) .< eps)
end

function test_blob_reshape(backend::Backend)
  test_blob_reshape(backend, Float32, 1e-10)
  test_blob_reshape(backend, Float64, 1e-10)
end

if test_cpu
  test_blob_reshape(backend_cpu)
end
if test_gpu
  test_blob_reshape(backend_gpu)
end


text:
function grad_y(x_hat_l,grad_x_hat_l,G, dG)
    m = size(G,1)
    return -kron(x_hat_l', I(m))*dG'-G*grad_x_hat_l'
end

function grad_x_hat(x_hat_l, dA,A,grad_x_hat_l, y, grad_y, K,dK)
    n = size(A,1)
    return kron(x_hat_l',I(n))*dA' + A*grad_x_hat_l' + kron(y',I(n))*dK' + K*grad_y'
end

function diff_K(A,P,G,S_l,dA,dG,dP,dS_l)
    m = size(G,1)
    n = size(A,1)
    P_inv = inv(P)
    Kmn = commutation_matrix(m,n)
    return kron(P_inv*G*S_l,I(n))*dA' + kron(P_inv*G,A)*dS_l'+kron(P_inv,A*S_l)*Kmn*dG' - kron(P_inv,A*S_l*G'*P_inv)*dP'
end

function diff_S(A,P,K,S_l,dA,dP,dK,dS_l,dQ)
    m = size(P,1)
    n = size(A,1)
    Knn = commutation_matrix(n,n)
    Kmn = commutation_matrix(m,n)
    b1 = (kron(A*S_l,I(n)) + kron(I(n),S_l*A')*Knn)*dA'
    b2 = (kron(K*P,I(n)) + kron(I(n),K*P)*Kmn)*dK'
    return b1 + kron(A,A)*dS_l' -b2 - kron(K,K)*dP' +dQ'
end

function diff_P(G,S_l,dG,dS_l,dR)
    m = size(G,1)
    n = size(G,2)
    Kmn = commutation_matrix(m,n)
    b1 = (kron(G*S_l,I(m)) + kron(I(m),G*S_l)*Kmn)*dG'
    return b1 + kron(G,G)*dS_l' + dR'
end

function diff_ll(P,y,grad_y,dP)
    P_inv = inv(P)
    return -1/2*det(P)*vec(P_inv)'*dP' - y'*P_inv*grad_y' + 1/2*y'*P_inv*dP'*P_inv*y
end

function diff_S0(A,S0,dA,dQ)
    n = size(A,1)
    Knn = commutation_matrix(n,n)
    b1 = (kron(A*S0,I(n)) + kron(I(n),A*S0)*Knn)*dA'
    ba = I(n^2) - kron(A,A)
    bb = b1 + dQ'
    return inv(ba)*bb
end

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
####
#   IO

function subscriptify(n::Int)
    subscript_0 = Int(0x2080) # Char(0x2080) -> subscript 0
    return join((Char(subscript_0 + i) for i in reverse(digits(n))))
end

function superscriptify(n::Int)
    superscripts = Dict(
        0 => "⁰",
        1 => "¹",
        2 => "²",
        3 => "³",
        4 => "⁴",
        5 => "⁵",
        6 => "⁶",
        7 => "⁷",
        8 => "⁸",
        9 => "⁹",
    )
    return join(superscripts[d] for d in reverse(digits(n)))
end

function Base.show(io::IO, α::Cyclotomic{T}) where {T}
    α = reduced_embedding(α)
    ζ = "ζ" * subscriptify(conductor(α))
    if iszero(α)
        print(io, zero(T))
    else
        for (i, exp) in enumerate(exponents(α))
            coeff = α[exp]
            if iszero(exp)
                print(io, coeff)
                continue
            end
            if isone(coeff) && T<:Integer
                sign_str = i == 1 ? " " : " + "
                coeff_str = ""
            elseif isone(-coeff) && T<:Integer
                sign_str = "-"
                coeff_str = ""
            elseif coeff > zero(coeff)
                sign_str = i == 1 ? " " : " + "
                coeff_str = "$coeff*"
            else
                sign_str = i == 1 ? "" : " "
                coeff_str = "$coeff*"
            end

            exp_str = isone(exp) ? "" : "$(superscriptify(exp))"
            print(io, sign_str, coeff_str, ζ, exp_str)
        end
    end
end

function Base.print(io::IO, α::Cyclotomic)
    if iszero(α)
        print(io, zero(valtype(α)))
    else
        E_str = "E($(conductor(α)))"
        for (i, exp) in enumerate(exponents(α))
            coeff = α[exp]
            if coeff > zero(coeff)
                sign_str = i == 1 ? " " : " + "
            else
                sign_str = i == 1 ? "" : " "
            end

            print(io, sign_str, coeff, "*", E_str, "^", exp)
        end
    end
end


text:
module TestPackage1
export hello, domath

hello(who::String) = "Hello, $who"
domath(x::Number) = x + 5

end


text:
using CSV, GLM, LinearAlgebra, StatsPlots, Plots, Measures; pyplot()

df = CSV.read("../data/weightHeight.csv")
n = size(df)[1]

model = lm(@formula(Height ~ Weight), df)
MSE = deviance(model)/dof_residual(model)
pred(x) = coef(model)'*[1, x]

A = [ones(n) df.Weight]
H = A*pinv(A)
residuals = (I-H)*df.Height
studentizedResiduals = residuals ./ (sqrt.(MSE*(1 .- diag(H))))

tau = rank(H)
println("tau = ", tau, " or ", sum(diag(H)), " or ", tr(H))

cookDistances = (studentizedResiduals.^2/tau) .* diag(H) ./ (1 .- diag(H))
maxCook, indexMaxCook = findmax(cookDistances)
println("Maximal Cook's distance = ", maxCook)

p1 = plot([minimum(df.Weight),maximum(df.Weight)],[0,0],
		lw = 2, c=:black, label=:none)
     scatter!([df.Weight[indexMaxCook]],[studentizedResiduals[indexMaxCook]],
                c=:red, ms = 10, msw = 0, shape =:cross,
                label = "Point with maximal Cook's distance",legend=:topleft)

    scatter!(df.Weight, studentizedResiduals, xlabel = "Weight (kg)",
             ylabel = "Studentized Residual (cm)",c=:blue,msw=0,label=:none)

p2 = qqnorm(studentizedResiduals, msw=0, lw=2, c =[:red :blue],legend = false,
        xlabel="Normal Theoretical Quantiles",
        ylabel="Studentized Residual Quantiles")

plot(p1,p2,size=(1000,500),margin = 5mm)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
# ============================================
# Flowpipes not contiguous in time
# ============================================

"""
    MixedFlowpipe{N, D, FT<:AbstractFlowpipe, VOA<:VectorOfArray{N, D, Vector{FT}}} <: AbstractFlowpipe

Type that wraps a vector of flowpipes of the same time, such that they are
not necessarily contiguous in time.

### Fields

- `Fk`  -- vector of flowpipes
- `ext` -- (optional, default: empty) dictionary for extensions

### Notes

This type does not assume that the flowpipes are contiguous in time.
"""
struct MixedFlowpipe{N, RT<:AbstractReachSet{N}, FT<:AbstractFlowpipe} <: AbstractFlowpipe
    Fk::VectorOfArray{RT, 2, Vector{FT}}
    ext::Dict{Symbol, Any}
end

function MixedFlowpipe(Fk::Vector{FT}) where {N, RT<:AbstractReachSet{N}, FT<:Flowpipe{N, RT}}
    voa = VectorOfArray{RT, 2, Vector{FT}}(Fk)
    ext = Dict{Symbol, Any}()
    return MixedFlowpipe{N, RT, FT}(voa, ext)
end

function MixedFlowpipe(Fk::Vector{FT}, ext::Dict{Symbol, Any}) where {N, RT<:AbstractReachSet{N}, FT<:Flowpipe{N, RT}}
    voa = VectorOfArray{RT, 2, Vector{FT}}(Fk)
    return MixedFlowpipe{N, RT, FT}(voa, ext)
end

# interface functions
array(fp::MixedFlowpipe) = fp.Fk
flowpipe(fp::MixedFlowpipe) = fp
setrep(::Type{MixedFlowpipe{N, RT, FT}}) where {N, RT, FT} = RT
numrsets(fp::MixedFlowpipe) = mapreduce(length, +, fp)

# indexing: fp[j, i] returning the j-th reach-set of the i-th flowpipe
Base.getindex(fp::MixedFlowpipe, I::Int...) = getindex(fp.Fk, I...)

@inline tspan(fp::MixedFlowpipe) = error("for mixed flowpipes you should specify the index, as in `tspan(fp, i)`")
@inline tspan(fp::MixedFlowpipe, i::Int) = tspan(fp[i])

function Base.similar(fp::MixedFlowpipe{N, RT, FT}) where {N, RT, FT}
    return MixedFlowpipe(Vector{FT}())
end

function (fp::MixedFlowpipe)(t::Number)
    error("not implemented yet")
end

function (fp::MixedFlowpipe)(dt::TimeInterval)
    error("not implemented yet")
end

function overapproximate(fp::MixedFlowpipe, args...)
    return MixedFlowpipe([overapproximate(Fi, args...) for Fi in fp], fp.ext)
end

function project(fp::MixedFlowpipe, args...)
    return [project(F, args...) for F in fp]
end

# LazySets interface
function LazySets.ρ(d::AbstractVector, fp::MixedFlowpipe)
    return maximum(ρ(d, F) for F in array(fp))
end

function LazySets.σ(d::AbstractVector, fp::MixedFlowpipe)
    error("not implemented")
end

function LazySets.linear_map(M, fp::MixedFlowpipe)
    out = [linear_map(M, F) for F in fp]
    return MixedFlowpipe(out, fp.ext)
end

function LazySets.affine_map(M, b, fp::MixedFlowpipe)
    out = [affine_map(M, b, F) for F in fp]
    return MixedFlowpipe(out, fp.ext)
end


text:
using Test
using GraphOrdering: symrcm

# https://people.sc.fsu.edu/~jburkardt/m_src/rcm/rcm.html
#
#    5--7--6
#    |  | /
# 4--8--2
# |  |  |
# 9--1--3

G = Dict(
    1 => [9, 8, 3],
    2 => [6, 7, 8, 3],
    3 => [2, 1],
    4 => [8, 9],
    5 => [7, 8],
    6 => [7, 2],
    7 => [5, 6, 2],
    8 => [1, 2, 4, 5],
    9 => [4, 1])

results = symrcm(G, 9)
@test results.perm == [6, 7, 2, 5, 3, 8, 1, 4, 9]


text:
module TestAvgFilter1D

using TapirBenchmarks
using Test

@testset for n in [3, 5, 100, 2^8]
    ys0, xs = avgfilter1d_setup(n)
    ys1 = avgfilter1d_seq!(copy(ys0), xs)
    @test avgfilter1d_threads!(copy(ys0), xs) == ys1
    @test avgfilter1d_tapir_dac!(copy(ys0), xs) == ys1
    @test avgfilter1d_tapir_seq!(copy(ys0), xs) == ys1
end

end


text:
@doc raw"""
    gausschebyshev(n::Integer) -> Tuple{Vector{Float64},Vector{Float64}}
    gausschebyshev(n::Integer, 1) -> Tuple{Vector{Float64},Vector{Float64}}

Return nodes and weights of [Gauss-Chebyshev quadrature](https://en.wikipedia.org/wiki/Chebyshev%E2%80%93Gauss_quadrature) of the 1st kind.

```math
\int_{-1}^{1} \frac{f(x)}{\sqrt{1-x^2}} dx \approx \sum_{i=1}^{n} w_i f(x_i)
```

# Examples
```jldoctest
julia> x, w = gausschebyshev(3);

julia> f(x) = x^4;

julia> I = dot(w, f.(x));

julia> I ≈ 3π/8
true
```

---

    gausschebyshev(n::Integer, 2) -> Tuple{Vector{Float64},Vector{Float64}}

Return nodes and weights of [Gauss-Chebyshev quadrature](https://en.wikipedia.org/wiki/Chebyshev%E2%80%93Gauss_quadrature) of the 2nd kind.

```math
\int_{-1}^{1} f(x)\sqrt{1-x^2} dx \approx \sum_{i=1}^{n} w_i f(x_i)
```

# Examples
```jldoctest
julia> x, w = gausschebyshev(3, 2);

julia> f(x) = x^4;

julia> I = dot(w, f.(x));

julia> I ≈ π/16
true
```

---

    gausschebyshev(n::Integer, 3) -> Tuple{Vector{Float64},Vector{Float64}}

Return nodes and weights of [Gauss-Chebyshev quadrature](https://en.wikipedia.org/wiki/Chebyshev%E2%80%93Gauss_quadrature) of the 3rd kind.

```math
\int_{-1}^{1} f(x)\sqrt{1-x^2} dx \approx \sum_{i=1}^{n} w_i f(x_i)
```

# Examples
```jldoctest
julia> x, w = gausschebyshev(3, 3);

julia> f(x) = x^4;

julia> I = dot(w, f.(x));

julia> I ≈ 3π/8
true
```

---

    gausschebyshev(n::Integer, 4) -> Tuple{Vector{Float64},Vector{Float64}}

Return nodes and weights of [Gauss-Chebyshev quadrature](https://en.wikipedia.org/wiki/Chebyshev%E2%80%93Gauss_quadrature) of the 4th kind.

```math
\int_{-1}^{1} f(x)\sqrt{1-x^2} dx \approx \sum_{i=1}^{n} w_i f(x_i)
```

# Examples
```jldoctest
julia> x, w = gausschebyshev(3, 4);

julia> f(x) = x^4;

julia> I = dot(w, f.(x));

julia> I ≈ 3π/8
true
```
"""
function gausschebyshev(n::Integer, kind::Integer=1)
    # GAUSS-CHEBYSHEV NODES AND WEIGHTS.

    if n < 0
        throw(DomainError(n, "Input n must be a non-negative integer"))
    end

    # Use known explicit formulas. Complexity O(n).
    if kind == 1
        # Gauss-ChebyshevT quadrature, i.e., w(x) = 1/sqrt(1-x^2)
        return ([cos((2 * k - 1) * π / (2 * n)) for k = n:-1:1], fill(π / n, n))
    elseif kind == 2
        # Gauss-ChebyshevU quadrature, i.e., w(x) = sqrt(1-x^2)
        return ([cos(k * π / (n + 1)) for k = n:-1:1],
                [π/(n + 1) * sin(k / (n + 1) * π)^2 for k = n:-1:1])
    elseif kind == 3
        # Gauss-ChebyshevV quadrature, i.e., w(x) = sqrt((1+x)/(1-x))
        return ([cos((k - .5) * π / (n + .5)) for k = n:-1:1],
                [2π / (n + .5) * cos((k - .5) * π / (2 * (n + .5)))^2 for k = n:-1:1])
    elseif kind == 4
        # Gauss-ChebyshevW quadrature, i.e., w(x) = sqrt((1-x)/(1+x))
        return ([cos(k * π / (n + .5)) for k = n:-1:1],
                [2π / (n + .5) * sin(k * π / (2 * (n + .5)))^2 for k = n:-1:1])
    else
        throw(ArgumentError("Chebyshev kind should be 1, 2, 3, or 4"))
    end
end


text:
"""
    @enter ex

Step into `ex`.
"""
macro enter(ex)
  Main.Atom.enter(__module__, ex)
end

"""
    @run ex

Run `ex` with the Interpreter and break when a breakpoint is encountered.
"""
macro run(ex)
  Main.Atom.enter(__module__, ex; initial_continue = true)
end

"""
    @trace ex

Analyse `ex` for type instabilities using Traceur.jl.
"""
macro trace(ex, args...)
  Main.Atom.trace(ex, args...)
end

function connect(args...; kws...)
  activate()
  eval(:(Main.Atom.connect($args...; $kws...)))
  return
end


text:
module LoggingFacilities

using Logging
using LoggingExtras
using Dates
using JSON

using Logging: Info
import Logging: shouldlog, min_enabled_level, catch_exceptions, handle_message

export TimestampLoggingFormat, OneLineLoggingFormat, JSONLoggingFormat
export logger
export SimplestLogger
export InjectByPrependingToMessage, InjectByAddingToKwargs

include("formats.jl")
include("simple.jl")

end # module


text:
using Dates, Statistics

##
function agg_mean(df::DataFrame, d::Vector{Symbol}, y::String)
	if y=="day"
		df[:,:dateTime] = floor.(df[:,:dateTime], Day(1))
	elseif y=="week"
		df[:,:dateTime] = floor.(df[:,:dateTime], Week(1))
	elseif y=="month"
		df[:,:dateTime] = floor.(df[:,:dateTime], Month(1))
	elseif y=="year"
		df[:,:dateTime] = floor.(df[:,:dateTime], Year(1))
	else
		println("Pick year,month or day")
	end
combine(groupby(df, [:dateTime]), d .=> mean)
end
###

##
function agg_sum(df::DataFrame, d::Vector{Symbol}, y::String)
	if y=="day"
		df[:,:dateTime] = floor.(df[:,:dateTime], Day(1))
	elseif y=="week"
		df[:,:dateTime] = floor.(df[:,:dateTime], Week(1))
	elseif y=="month"
		df[:,:dateTime] = floor.(df[:,:dateTime], Month(1))
	elseif y=="year"
		df[:,:dateTime] = floor.(df[:,:dateTime], Year(1))
	else
		println("Pick year,month or day")
	end
combine(groupby(df, [:dateTime]), d .=> sum)
end
###

##
function agg_max(df::DataFrame, d::Vector{Symbol}, y::String)
	if y=="day"
		df[:,:dateTime] = floor.(df[:,:dateTime], Day(1))
	elseif y=="week"
		df[:,:dateTime] = floor.(df[:,:dateTime], Week(1))
	elseif y=="month"
		df[:,:dateTime] = floor.(df[:,:dateTime], Month(1))
	elseif y=="year"
		df[:,:dateTime] = floor.(df[:,:dateTime], Year(1))
	else
		println("Pick year,month or day")
	end
combine(groupby(df, [:dateTime]), d .=> maximum)
end
############################################################
function agg_min(df::DataFrame, d::Vector{Symbol}, y::String)
	if y=="day"
		df[:,:dateTime] = floor.(df[:,:dateTime], Day(1))
	elseif y=="week"
		df[:,:dateTime] = floor.(df[:,:dateTime], Week(1))
	elseif y=="month"
		df[:,:dateTime] = floor.(df[:,:dateTime], Month(1))
	elseif y=="year"
		df[:,:dateTime] = floor.(df[:,:dateTime], Year(1))
	else
		println("Pick year,month or day")
	end
combine(groupby(df, [:dateTime]), d .=> minimum)
end

#############################################################

max_stars_repo_name: atott/Fitbit-Julia

max_issues_repo_name: atott/Fitbit-Julia

max_forks_repo_name: atott/Fitbit-Julia


text:
"""
    Identity(:anything)

Identity is a simple wrapper, which works as a single-element container.

It can be used as the trivial Monad, and as such can be helpful in monadic
abstractions. For those who don't know about Monads, just think of it like
container-abstractions.
"""
struct Identity{T}
  value::T
  Identity(value) = new{typeof(value)}(value)
end
# == controversy https://github.com/JuliaLang/julia/issues/4648
Base.:(==)(a::Identity, b::Identity) = a.value == b.value

function Base.show(io::IO, x::Identity)
  print(io, "Identity($(repr(x.value)))")
end

"""
    isidentity(Identity(3)) -> true
    isidentity("anythingelse") -> false

returns true only if given an instance of  [Identity](@ref)
"""
isidentity(::Identity) = true
isidentity(other) = false

Base.length(::Identity) = 1

Base.get(a::Identity) = a.value
Base.getindex(a::Identity) = a.value

Base.eltype(::Type{<:Identity{T}}) where T = T
Base.eltype(::Type{<:Identity}) = Any

Base.iterate(a::Identity) = a.value, nothing
Base.iterate(a::Identity, state::Nothing) = state
Base.foreach(f, a::Identity) = begin f(a.value); nothing; end
Base.map(f, a::Identity) = Identity(f(a.value))
# for convenience, Identity does not use convert, whatever monad is returned is valid, providing maximum flexibility.
Base.Iterators.flatten(a::Identity) = a.value

Base.convert(::Type{Identity{T}}, x::Identity{T}) where {T} = x
Base.convert(::Type{Identity{T}}, x::Identity) where {T} = Identity(convert(T, x.value))

# Identity is covariate
# promote_rule only works on concrete types, more general checks Type{<:Const} may overwrite
# unintentionally more specific promote_rule types
Base.promote_rule(::Type{Identity{T}}, ::Type{Identity{S}}) where {T, S<:T} = Identity{T}
Base.promote_rule(::Type{Identity{T}}, ::Type{Identity}) where T = Identity
Base.promote_rule(::Type{Identity}, ::Type{Identity}) = Identity

# we need this for safety, if someone overwrites typejoin for Unions with Identiy
Base.promote_typejoin(::Type{Identity{T}}, ::Type{Identity{T}}) where T = Identity{T}
Base.promote_typejoin(::Type{<:Identity}, ::Type{<:Identity}) = Identity


text:
### LLVM

# sample code
julia> function fib(n::Int64)
           if n == 1
               return n
           elseif n == 0
               return 0
           else
               return fib(n-1) + fib(n-2)
           end
       end
fib (generic function with 1 method)

julia> fib(10)
55

# Inspect in Julia
julia> @code_llvm(fib(10))

define i64 @julia_fib_62805.2(i64) #0 !dbg !5 {
top:
  %1 = icmp eq i64 %0, 1
  br i1 %1, label %if, label %L4

if:                                               ; preds = %L4, %top
  %merge = phi i64 [ 1, %top ], [ 0, %L4 ]
  ret i64 %merge

L4:                                               ; preds = %top
  %2 = icmp eq i64 %0, 0
  br i1 %2, label %if, label %L9

L9:                                               ; preds = %L4
  %3 = add i64 %0, -1
  %4 = call i64 @julia_fib_62805(i64 %3)
  %5 = add i64 %0, -2
  %6 = call i64 @julia_fib_62805(i64 %5)
  %7 = add i64 %6, %4
  ret i64 %7
}

### Parallel Computing
# current processes
julia> nprocs()
2

# add 2 more processes
julia> addprocs(2)
2-element Array{Int64,1}:
 3
 4

# check the updated processes again
julia> nprocs()
4

# Create a remotecall
julia> remote_call = remotecall(+, 2, 2, 3)
Future(2,1,5,Nullable{Any}())

julia> interim_res = @spawnat 2 1 + fetch(remote_call)
Future(2,1,6,Nullable{Any}())

julia> fetch(interim_res)
6

# running on a single process, i.e when workers are just 1.

julia> workers()
1-element Array{Int64,1}:
 1

# simple for loop, works perfectly
julia> for i in 1:5
           println(i)
       end
1
2
3
4
5

# parallel for loop running on sigle worker doesn't give a
# synchronised output
julia> @parallel for i in 1:5
           println(i)
       end
1
1-element Array{Future,1}2:
3
4Future(1,1,1,#NULL)
5

julia> addprocs(2)
2-element Array{Int64,1}:
 2
 3

julia> nprocs()
3

julia> workers()
2-element Array{Int64,1}:
 2
 3

julia> result = @parallel for i in 1:5
           println(i)
       end
2-element Array{Future,1}:
 Future(3,1,5,#NULL)
 Future(2,1,6,#NULL)

julia>  From worker 3:  1
        From worker 3:  2
        From worker 3:  3
        From worker 2:  4
        From worker 2:  5
julia>

julia> pmap(sum, [1 2 3 4 5])
5-element Array{Any,1}:
 1
 2
 3
 4
 5

 ### TCP Sockets & Servers
 # define the PORT
PORT = 7575

# create the server
server = listen(PORT)

# create the connection
connection = accept(server)

# read the input from the user
aline = readline(connection)

# write the following to the receiver's screen
write(connection, "Hey, how are you?")

# close the connection immediatly
close(connection)

# Sockets
julia> u = UDPSocket()
UDPSocket(init)

julia> bind(u, IPv4(127,0,0,1), 5001)
true

julia> v = UDPSocket()
UDPSocket(init)

julia> bind(v, IPv4(127,0,0,1), 5002)
true

julia> send(v,ip"127.0.0.1",5001,"Hey there!")

julia> msg = recv(u)
10-element Array{UInt8,1}:
 0x48
 0x65
 0x79
 0x20
 0x74
 0x68
 0x65
 0x72
 0x65
 0x21

julia> String(msg)
"Hey there!"

julia> close(u)

julia> close(v)

max_stars_repo_name: PacktPublishing/Learning-Julia

max_issues_repo_name: PacktPublishing/Learning-Julia

max_forks_repo_name: PacktPublishing/Learning-Julia


text:
module NSBEntropy
import GSL
import Base.Test
import StatsBase
using Distributions,StatsBase

s_0(B,K) = polygamma(0,K*B+1) - polygamma(0,B+1) # <S(n=0)>
sigma_0(B,K) = (B+1)/(B*K+1)*polygamma(1,B+1) - polygamma(1,B*K+1)
d_s_0(B,K) = K*polygamma(1,K*B+1) - polygamma(1,B+1)
#s_1(N,K,n) = sum((n+1)./(N+K).*(map(x-> polygamma(1,x+2),n)-polygamma(1,N+K+1))) #<S(n)>/<S(n=0)>
#S_1(N,K,n,B) = s_0(B,K)*s_1(N,K,n) # <S(n)>
#rhoB(B,n) = gamma(BigFloat(length(n)*B))/gamma(BigFloat(sum(n)+length(n)*B)).*prod(map(x->gamma(BigFloat(x+B))/gamma(B),n))
#rho(B,n) = gamma(length(n)*B)/gamma(sum(n)+length(n)*B).*prod(map(x->gamma(x+B)/gamma(B),n))
#I(B,K,N,n) = rho(B,n)*S1(N,K,n,B)*d_s_0(B,K)

I{T<:Integer}(B::Real,n::Array{T,1}) = I(B,n,length(n))

function series{T<:Real}(sz::Integer, a::Array{T,1},x::Real,n::Integer)
	results = 0.0
	for i in sz:-1:1
		results += a[i]*x^(i+n)
	end
	return results
end

function I{T<:Integer}(B::Real,n::Array{T,1},K::Integer)
    N = sum(n)
    r = rho(B,n)
    s = s1(n,B)
    d = d_s_0(B,K) #
    return r*s*d
end

function simplefunc{T<:Real}(n::Array{T,1})
    #gives the entropy, s1 and the variance of the estimate, s2
    #these are valid up to zeroth order in 1/K and 1/N, i.e. for both N and K large
    N = sum(n)
	K1 = sum(n.>=1)
    Delta = N-K1 #number of coincidences
    s1 = (eulergamma - log(2)) + 2*log(N) - polygamma(0,Delta) #Nemenman(2002, Eq. 30)
    s2 = polygamma(1,Delta) #Nemenman(2002, Eq. 31)
    return s1/log(2), s2/log(2)
end

function saddlepoint2{T}(n::Array{T,1}, K::Integer)
    N = sum(n)
    K1 = sum(n.>=1) #number of coincidences
    K2 = sum(n.>=2)
    d = (N - K1)/N #relative number of coincidences
	#println(d)
	order = 10
	N2 = N*N
	N3 = N2*N
	N4 = N3*N
	N5 = N4*N
	N6 = N5*N
	N7 = N6*N
	N8 = N7*N
	N9 = N8*N
	N10 = N9*N
	N11 = N10*N
	N12 = N11*N

	Nm = N-1.0
	overN = 1/N
	Nm2 = Nm^2
	Nm3 = Nm^3
	Nm4 = Nm^4
	Nm5 = Nm^5
	Nm6 = Nm^6
	Nm7 = Nm^7
	Nm8 = Nm^8
	Nm9 = Nm^9
	Nm10 = Nm^10

    #b_1 =  (N-1.0)/(2*N)
    #b0 = (-2.0+1.0/N)/3
    #b1 = (N^2.0 - N - 2.0)/(9.0*(N^2 - N))
	b = [b_1 =  (N-1.0)/(2*N),
        (-2.0+1.0/N)/3,
        (N^2.0 - N - 2.0)/(9.0*(N^2 - N)),
	    (2.0*(2.0 - 3.0*N - 3.0*N2 + 2.0*N3))/(135.0*Nm2*N),
	    (4.0*(22.0 + 13.0*N - 12*N2 - 2.0*N3 + N4))/(405.0*Nm3*N),
	    (4.0*(-40.0 + 58.0*N + 65.0*N2 - 40.0*N3 - 5.0*N4 +
	      2.0*N5))/(1701.0*Nm4*N),
	    (4.0*(-9496.0 - 6912.0*N + 5772.0*N2 + 2251.0*N3 -
	      1053.0*N4 - 87.0*N5 + 29.0*N6))/(42525.0*Nm5*N),
	    (16.0*(764.0 - 1030.0*N - 1434.0*N2 + 757.0*N3 + 295.0*N4
	       - 111.0*N5 - 7.0*N6 + 2.0*N7))/(18225.0*Nm6*N),
	    (16.0*(167000.0 + 142516.0*N - 108124.0*N2 - 66284.0*N3
	       + 26921.0*N4 + 7384.0*N5 - 2326.0*N6 - 116.0*N7
	       + 29.0*N8))/(382725*Nm7*N),
	    (16.0*(-17886224.0 + 22513608.0*N + 37376676.0*N2
	       - 17041380.0*N3 - 11384883.0*N4 + 3698262.0*N5
	       + 846930.0*N6 - 229464.0*N7 - 9387.0*N8 + 2086*N9))/
			(37889775.0*Nm8*N),
	    (16.0*(-4166651072.0 - 3997913072.0*N + 2783482560.0*N2 +
	       2290151964.0*N3 - 803439834.0*N4 - 395614251.0*N5 +
	       108055443.0*N6 + 20215218.0*N7 - 4805712.0*N8 - 165395.0*N9
	       + 33079.0*N10))/(795685275*Nm9*N),
	    (32.0*(52543486208.0 - 62328059360.0*N - 118489458160.0*N2 +
	       47185442088.0*N3 + 44875379190.0*N4 - 12359832987.0*N5 -
	       5400540075.0*N6 + 1272974916.0*N7 + 200644800.0*N8 -
	       42495955.0*N9 - 1255067.0*N10 +
	       228194.0*N11))/(14105329875.0*Nm10*N)]

	B0 = N*series(order+2, b, (N-K1)/N, -1)
	@assert B0 > 0
	b_1 = b[1]
	b0 = b[2]
	b1 = b[3]
end

function mlog_evidence(n::Array{Int64,1},Β::Real)
    N = sum(n)
    n1 = n[n.>=1] #coincidences
	K1 = length(n1)
	f = 0.0
	if K1 > 0.0
		f += lgamma(n1)
	end
end

function saddlepoint{T}(n::Array{T,1}, K::Integer)
    N = sum(n)
    K1 = sum(n.>=1) #number of coincidences
    K2 = sum(n.>=2)
	b_1 =  (N-1.0)/(2*N)
    b0 =   (-2.0+1.0/N)/3
    b1 =   (N^2.0 - N - 2.0)/(9.0*(N^2 - N))
    d = (N - K1)/N #relative number of coincidences
    k0 = N*(b_1/d + b0 + b1*d)
    idx = find(n.>1)
    k1 = 0.0
    for i=1:length(idx)
        k1 += (polygamma(0,n[idx[i]]) - polygamma(0,1))/(K1/k0^2 - polygamma(1,k0) + polygamma(1,k0 + N))
    end
    k2 = (K1/k0^3 + (polygamma(2,k0) - polygamma(2,k0+N))/2)*k1^2
    for i=1:length(idx)
        k2 += k0*(polygamma(1,n[idx[i]]) - polygamma(1,1))
    end
    k2 /= (K1/k0^2 - polygamma(1,k0) + polygamma(1,k0 + N))
    ks = k0 + k1/K1 + k2/K2
    println("k₀ = $k0, k₁ = $k1, k₂ = $k2, d = $d, b_1 = $b_1, b0 = $b0, b1 = $b1")
    return ks/K #return Beta
end

rho{T<:Integer}(B::Real, n::Array{T,1}) = rho(B,n,length(n))
rho{T<:Integer}(B::Real, n::Array{T,1}) = exp(logrho(B,n,length(n)))

logrho{T<:Integer}(B::Real, n::Array{T,1}) = logrho(B, n,length(n))

function logrho{T<:Integer}(B::Real, n::Array{T,1},K::Integer)
    N = sum(n)
    C1 = GSL.sf_lngamma(B*K) - GSL.sf_lngamma(N+B*K)
    C2 = 0.0
    for i=1:K
        C2 += GSL.sf_lngamma(n[i] + B) - GSL.sf_lngamma(B)
    end
    return C1+C2
end

s1{T<:Integer}(n::Array{T,1},B::Real) = s1(n,B,length(n))

function s1{T<:Integer}(n::Array{T,1},B::Real, K::Integer)
    N = sum(n)
    C = 0.0
    for i=1:K
        C += (n[i]+B)*polygamma(0,n[i]+B+1)
    end
    C -= (N+B*K)*polygamma(0,N+K*B+1)
    return -C/(N+B*K)
end

function s1{T<:Integer}(n::Array{T,1})
    N = sum(n)
    K = length(n)
    C = 0.0
    for i=1:K
        C += (n[i]+1)*polygamma(0,n[i]+2)
    end
    C -= (K+N)*polygamma(0,N+K+1)
    return -C/(N+K)
end

function s2{T<:Integer}(n::Array{T,1})
    N = sum(n)
    K = length(n)
    C = 0.0
    deltapg(k,x,y) = polygamma(k, x) - polygamma(k,y)
    rr = (N+K)*(N+K+1)
    for i=1:K
        C += (n[i]+1)*(n[i] + 2)*(deltapg(0,n[i]+3, N+K+2)^2 + deltapg(1,n[i]+3, N+K+2))/rr
        for j=1:K
            if i != j
                C += (n[i] + 1)*(n[j]+1)*(deltapg(0,n[i] + 2, N+K +2)*deltapg(0,n[j]+1, N+K+2) - polygamma(1,N+K+2))/rr
            end
        end
    end
    return C
end

function s2{T<:Integer}(n::Array{T,1},B::Real)
    N = sum(n)
    K = length(n)
    C = 0.0
    deltapg(k,x,y) = polygamma(k, x) - polygamma(k,y)
    rr = (N+B*K)*(N+B*K+1)
    for i=1:K
        C += (n[i]+B)*(n[i] + B + 1)*(deltapg(0,n[i]+B+2, N+B*K+2)^2 + deltapg(1,n[i]+B+2, N+B*K+2))/rr
        for j=1:K
            if i != j
                C += (n[i]+ B)*(n[j]+B)*(deltapg(0,n[i] +B + 1, N+B*K +2)*deltapg(0,n[j]+B, N+B*K+2) - polygamma(1,N+B*K+2))/rr
            end
        end
    end
    return C
end

function S1{T<:Integer}(B::Real, n::Array{T,1})
    C1 = s1(n,B)
    K = length(n)
    C2 = 1#s_0(B,K)
    return C1*C2
end

entropy{T<:Integer}(n::Array{T,1}) = getEntropy(n,1e-7,3.0,length(n))
entropy{T<:Integer}(n::Array{T,1},bim::Real, bmax::Real) = getEntropy(n,bmin,bmax,length(n))

function getEntropy{T<:Integer}(n::Array{T,1},bmin::Real,bmax::Real,K::Integer)
    #compute the normalization constant
    Z = quadgk(x -> rho(x,n)*d_s_0(x,K),bmin,bmax)
    #compute the integral over S
    Q = quadgk(x -> I(x,n),bmin,bmax)
    E = Q[1]/Z[1]
    #convert from natural logarithm to bits
    E /= log(2)
    dE = sqrt(Q[2]^2 + Z[2]^2)
    return E,dE
end

function test1()
    n = rand(1:5,100)
    q = counts(n)
    E = S_1(length(n),length(q),q,1.0)
    return E
end

function testEntropy()
    n = rand(1:10,100)
    q = counts(n)
    E = getEntropy(q)
    return E
end

function test_s_0()
    B = 0.1
    K = 100
    ss = s_0(B,K)
    Base.Test.@test_approx_eq ss 2.775507529477798
    println("Test passed. S₀ = $ss")
end

function test_sigma_0()
    B = 0.1
    K = 100
    ss = sigma_0(B,K)
    Base.Test.@test_approx_eq ss 0.04816357939759014
    println("Test passed. ơ₀ = $ss")
end

function test_s1()
    s = s1([0,2])
    Base.Test.@test_approx_eq s 0.45833333333333326
    println("Test passed. s₁ = $s")
end

function test_simplefunc()
    srand(1234)
    n = rand(1:1000,100)
    q = StatsBase.counts(n,1:1000)
    s,ds = simplefunc(q)
    st = log2(1000)
    Base.Test.@test_approx_eq s 12.510509693214221
    Base.Test.@test_approx_eq ds 0.930443179942287
    println("Test passed. s = $s, ds² = $ds. N/K = 0.1. True entropy: $st")
end

end


text:
"""
Package: Forecast

    loess(xv, yv;
          d = 0,
          q = 3*sum((!ismissing).(yv))÷4,
          rho = fill(1.0,length(yv)),
          exact = [], extra = [])

    loess(yv;
          d = 0,
          q = 3*sum((!ismissing).(yv))÷4,
          rho = fill(1.0,length(yv)),
          exact = [], extra = [])

Return a funcion to smooth a vector of observations using locally weighted regressions.

Although loess can be used to smooth observations for any given number of independent variables, this implementation is univariate. The speed of loess can be greatly increased by using fast aproximations for the linear fitting calculations, however this implementation calculates allows as well for exact results.

The loess functionality and nomenclature follows the descriptions in:

"STL: A Seasonal, Trend Decomposition Procedure Based on Loess"
Robert B. Cleveland, William S. Cleveland, Jean E. McRae, and Irma Terpenning.
Journal of Official Statistics Vol. 6. No. 1, 1990, pp. 3-73 (c) Statistics Sweden.

# Arguments
- `xv`: Observations' support, if not provided 1:length(yv) is used instead.
- `yv`: Observation values.
- `d`: Degree of the linear fit, it accepts values 0, 1 or 2, if 0 an estimation of `d` is calculated.
- `q`: As q increases loess becomes smoother, when q tends to infinity loess tends to an ordinary least square poynomial fit of degree `d`. It defaults to the rounding of 3/4 of xv's non-missing values length.
- `rho`: Weights expressing the reliability of the observations (e.g. if yi had variances σ^2*ki where ki where known, the rhoi could be 1/ki). It defaults to 1.
- `exact`: Vector containing support values where exact loess will be calculated, the least values the faster will be loess. It defaults to an emtpy array. For series of 10 values or less exact loess is guaranteed.
- `extra`: Vector containing support values where, besideds those values chosen by the heuristic, exact loess will be calculated.

# Returns
A function returning the exact loess for the values contained in `exact` and either exact loess or fast approximations for the rest.

# Examples
```julia-repl
n = 1000
x = sort(rand(n)*2*pi);
y = sin.(x) + randn(n);
f = Smoothers.loess(x,y)
isapprox(f(pi),0;atol=0.1) #true
[...]
```
"""
loess

function loess(yv::AbstractVector{<:Union{Missing,T}};
               d::Integer=0, # 0 -> auto
               q::Integer=3*length(yv)÷4,
               rho::AbstractVector{<:Union{Missing,T}}=fill(T(1),length(yv)),
               exact::AbstractVector{T}=T[],
               extra::AbstractVector{T}=T[]) where {T<:Real,S<:Real}

    loess(T(1):length(yv),yv;d,q,rho,exact,extra)
end

function loess(xv::AbstractVector{R},
               yv::AbstractVector{<:Union{Missing,T}};
               d::Integer=0,
               q::Integer=3*length(yv)÷4,
               rho::AbstractVector{<:Union{Missing,T}}=fill(T(1.0),length(xv)),
               exact::AbstractVector{R}=R[],
               extra::AbstractVector{R}=R[]) where {R<:Real,T<:Real}

    @assert issorted(xv) "xv values must be sorted"
    @assert d in 0:2 "Linear Regression must be of degree 1 or 2 or 0 (auto)"
    @assert length(xv) == length(yv) "xv and yv must have the same length"

    # Removing missing values
    myi = findall(x -> !ismissing(x),yv)
    @assert length(myi) > 0 "at least one yv value must be non missing"
    q = Int(floor(q*length(myi)/length(yv)))
    xv = xv[myi]
    yv = Vector{T}(yv[myi])
    rho = rho[myi]

    length(xv) == 1 && return x -> repeat([yv[1]],length(x))
    length(xv) == 2 && begin d=(yv[2]-yv[1])/(xv[2]-xv[1]); return x -> @. d*x+(yv[1]-d*xv[1]) end

    # Promote to same Float type
    xv,yv,rho,exact,extra,_ = Base.promote(xv,yv,rho,collect(exact),collect(extra),
                                           [Base.promote_op(/,T,T)(1.0)])

    ## Ax = b prediction
    P = eltype(xv)
    A = hcat(xv,fill(P(1.0),length(xv)))
    b = yv

    ## Order Estimation
    d = (d == 0) ? min(khat(yv),2) : d  # loess accepts 1 or 2
    d == 2 && (A = hcat(xv .^ 2, A))

    # Select collection
    predictx = exact
    if isempty(exact)
        N = length(xv)
        n = 2*sr(N)*N÷q
        m,M = extrema(xv)
        gap = P((M-m)/n)
        #exact = vcat(exact, xv[Int64.(round.(LinRange(1,length(xv),20)))])
        predict = sort(unique(vcat(collect(m:gap:M),M,extra)))
        predictx = vcat(m-P(10.0)*gap,m-gap,
                        predict,
                        M+gap,M+P(10.0)*gap)
    else
        @assert length(exact) > d "length(exact) should be greater than d"
    end

    # Predict collection
    res = similar(predictx)
    for (i,xi) in enumerate(predictx)
        res[i] = ghat(xi;A,b,d,q,rho)
    end

    Spline1D(predictx, res; k=d, bc="extrapolate",s=0)

end

# loess estimation
function ghat(x::T;
              A::AbstractMatrix{T},
              b::AbstractVector{T},
              d::Integer,
              q::Integer,
              rho::AbstractVector{T}) where T<:Real

    xv = A[:,d]
    yv = b

    ## λ_q
    n = length(xv)
    xvx = @. abs(xv-x)
    qidx = sortperm(xvx)[1:min(q,n)]
    qdist = abs(xv[last(qidx)]-x)*max(T(1),q/n)

    ## upsilon
    w = zeros(n)
    for wi in qidx
        aq = abs(xv[wi]-x)/qdist
        w[wi] = max((1-aq^3)^3,T(0))
    end

    A = @. A*(w*rho)
    b = @. b*(w*rho)

    lsq_x = A\b

    d == 1 ? sum([x,T(1)].*lsq_x) : sum([x^2,x,T(1)].*lsq_x)

end

#Sturges' Rule
sr(n) = Int(round(1+3.322*log(n)))

# absolute standard deviation
function astd(x::AbstractVector{<:Union{Missing,Real}})
    astd(collect(skipmissing(x)))
end

function astd(x::AbstractVector{<:Real})
    n = length(x)
    sum(abs.(x.-sum(x)/n))/n
end

# k order estimation for Splines
function khat(y::AbstractVector{<:Union{Missing,Real}})
    khat(collect(skipmissing(y)))
end

function khat(y::AbstractVector{<:Real})
    dy = diff(y); ady = astd(dy)
    ady < 1e-10 && return 1
    mi = 1
    for i in 1:3
        dy = diff(dy)
        nady = astd(dy); nady < 1e-10 && break
        mi += (ady > nady) ? 1 : break
        ady = nady
    end
    mi+1
end


text:
using JuMP, Gurobi, Combinatorics, Cbc
import XLSX

#### read data ####
xf = XLSX.readxlsx("../data/data2_simplified2.xlsx")
l = 8     # work day duration
N = xf["General!B3"] + 1  # number of clients + 1
M = xf["General!B2"]  # number of providers
# + 1 is for including the headquarter (HS14)

# Clients
s = xf["Clients!A2:A$N"]   # client service start times
d = vcat([0], xf["Clients!B2:B$N"])   # durations for clients, set d[1] = 0 for HS14

# Providers
w = xf["Providers!B2:B$(M+1)"]  # start time of providers
f = xf["Providers!A2:A$(M+1)"]  # hiring cost of providers

t = xf["Distances"][:]  # 2-d array of travel times

L = 10000 # large number
E = collect(powerset(collect(2:N), 2, N-1)) # for subtour 'E'limination constraints
# E = collect(powerset(collect(1:N), 2, N-1)) # for subtour 'E'limination constraints
s_test = ones(M, N)

for k in 1:M
    global s_test[k,:] = vcat([w[k]], s)
end
s = s_test

model = Model();
# Gurobi.setparam!(backend(model).optimizer.model.inner, "TimeLimit", 100.0)
set_optimizer(model, Gurobi.Optimizer);
set_time_limit_sec(model, 100.0);

#(13)
@variable(model, x[1:N, 1:N, 1:M], Bin)
@variable(model, 1 <= u[1:N, 1:M] <= N-1)
@variable(model, tl[1:N, 1:M] >= 0, Int)#Tardiness hours
@variable(model, te[1:N, 1:M] >= 0, Int)# early hours
@variable(model, td[1:N, 1:N, 1:M] >= 0, Int)# actual depature time
P = 30 # later arrival penalty cost = 30/hour
######################## routing constraints ######################

# #
# each client served once by one provider
@constraint(model, served[j in 2:N], sum(x[i,j,k] for  k in 1:M, i in 1:N if j != i) == 1)

# route continuity
# @constraint(model, continuity[p in 2:N, k in 1:M], sum(x[i,p,k] for i in 1:N) == sum(x[p,j,k] for j in 1:N) )
@constraint(model, continuity[j in 1:N, k in 1:M],(sum(x[i,j,k] for i in 1:N)-sum(x[j,i,k] for i in 1:N))==0 )

#  all providers only provide service once/ each vehicle is used at most once
@constraint(model, leave_depot[k in 1:M], sum(x[1,j,k] for j in 2:N) <= 1)

# subtour elimination for each provider
for k in 1:M
    for i in 2:N
        for j in 2:N
            @constraint(model, u[i,k] - u[j,k] + (N-1) * x[i,j,k] <= (N-1) - 1 )
        end
    end
end

#Either arrive early, arrive on time or late
@constraint(model, oneinthree[i in 1:N, j in 1:N, k in 1:M], tl[j,k]*te[j,k] == 0 )
# service provider depoart from depot if they work, and they depart as soon as they are available.
@constraint(model, work_when_avai[j in 1:N, k in 1:M], td[1,j,k] == x[1,j,k]*s[k,1])
#define the actual depature time
@constraint(model, actual_dep[i in 1:N, j in 1:N, k in 1:M], x[i,j,k]*(s[k,i]+tl[i,k]+d[i])==td[i,j,k])
# come back to 14HS by hour l
@constraint(model, back_depot[i in 1:N, k in 1:M], x[i,1,k]*(td[i,1,k]+t[i,1]) <= l)
# define the next actual arrival time
@constraint(model, late[i in 1:N, j in 2:N, k in 1:M], x[i,j,k]*(td[i,j,k]+t[i,k])==x[i,j,k]*(s[k,i]+tl[j,k]-te[j,k]))
######################## end time constraints ######################

# @objective(model, Min,sum(f[k]*x[1,j,k] for j in 1:N for k in 1:M))
# @objective(model, Min, sum(f[k]*x[i,1,k]*(td[i,1,k]+t[i,1]-s[k,1]) for i in 1:N for k in 1:M)+sum(tl[j,k]*P*sum(x[i,j,k] for i in N) for j in 1:N for k in 1:M))
@objective(model, Min,sum(x[1,j,k]*f[k]*(sum(x[p,q,k]*(t[p,q] + d[p] + max(0, s[k,q] - s[k,p] - t[p,q] - d[p])) for p in 1:N for q in 1:N)) for j in 2:N for k in 1:M) +sum(tl[j,k]*P*sum(x[i,j,k] for i in N) for j in N for k in M))
optimize!(model)

# @show (value.(x[:,:,1]))

for k in 1:M
    length = 0
    if sum(value.(x[1,:,k])) >= 1
        println("Provider $k provides service ")
        # come back hr
        for i in 1:N
            if value.(x[i,1,k]) >= 1
                println("Provider $k comes back to 14HS at ",
                value(sum(x[i,j,k]*( t[i,j] + d[j]) for i in 1:N, j in 1:N)),"hr.")
                # break
            end
        end
        for i in 1:N
            for j in 1:N
                if value.(x[i,j,k])>=1
                    println(i,"-",j)
                    length = length + t[i,j]
                end
            end
        end
    end
    println("travel length:", length)
    println()
end

for k in 1:M
    for j in 1:N
        if value.(tl[j,k])>0
            println("There is penalty cost for vehicle $k at customer $j with tardiness hour ", value.(tl[j,k]))
        end
    end
end


text:
# Checks if probabilities are strictly positive and if not adjusts slightly.
# What adjustments can be justified have to be determined in each instance.

function checkTree(tree0::Tree, eps)

# For every parent node going backwards
tmp = reverse(unique(tree0.parent)); pop!(tmp)
for i in tmp
    #Note that the probs are the conditional ones
    prob = tree0.probability[tree0.children[i+1]]
    if any(x->x>0,prob)
      #do nothing
    else
      prob = prob .+ eps
      prob = prob./ sum(prob)
      println()
      println("updated tree")
    end
  end

end


text:
csg_round_lin(x::T; sigdigits::Int = 9) where {T} = round(x, sigdigits = sigdigits) # min. Δx = 1 nm
csg_round_lin(x::Float32; sigdigits::Int = 5) = round(x, sigdigits = sigdigits) # min. Δx = 10 μm
csg_round_rad(φ::T; digits::Int = 4) where {T} = round(φ, digits = digits) # min. Δφ ≈ 0.011 °

csg_default_tol(::Type{T}) where {T} = T(1e-12)

csg_isapprox(x::T, y::T; csgtol::T = T(1e-9)) where {T} = abs(x - y) < csgtol


text:
using Test
using PrePostCall

struct NNegativeException <: Exception end
struct NSmallerException <: Exception end
struct NLargerException <: Exception end
struct NZeroException <: Exception end
struct BeforeException <: Exception end
struct FixLengthException <: Exception end

@pre function positive(x::Int)
    x<0 && throw(NNegativeException())
end

@pre function larger(x,y)
    x<y && throw(NSmallerException())
end

@post function not0(res::Int)
    res==0 && throw(NZeroException())
end

@post function smaller(x,y)
    x>y && throw(NLargerException())
end

@post function fixLength(v::Vector{Float64})
    length(v)>3 && throw(FixLengthException())
end

@post function multret(x,y)
    x>y && throw(NSmallerException())
end

module Level1
using PrePostCall
using Test
struct InModException <: Exception end
@pre function inmod(x::Int)
    throw(InModException())
end
@inmod function foo(x)
    x
end
function bar(x)
    foo(x)
end
function baz(x)
    @test_throws InModException foo(x)
end
end

@testset "PrePostCall" begin
    @testset "Pre-Single" begin
        @positive function foo(x::Int,y::T) where T<:Number
            return x^y
        end
        @positive a function foo2(a::Int,y::T) where T<:Number
            return a^y
        end
        @test foo(10,10) == 10^10
        @test_throws NNegativeException foo(-1,10)
        @test foo2(10,10) == 10^10
        @test_throws NNegativeException foo2(-1,10)
    end

    @testset "Pre-Two" begin
        @larger function foo(x::Int,y::T) where T<:Number
            return x^y
        end
        @larger a,y function foo2(a::Int,y::T) where T<:Number
            return a^y
        end
        @test foo(10,10) == 10^10
        @test_throws NSmallerException foo(-1,10)
        @test foo2(10,10) == 10^10
        @test_throws NSmallerException foo2(-1,10)
    end

    @testset "Pre-Multiple" begin
        @positive @larger function foo(x::Int,y::T) where T<:Number
            return x^y
        end
        @larger x,b @positive function foo2(x::Int,b::T) where T<:Number
            return x^b
        end
        @test foo(10,10) == 10^10
        @test foo2(10,10) == 10^10
        @test_throws NNegativeException foo(-1,10)
        @test_throws NSmallerException foo2(-1,10)
        @test_throws NNegativeException foo2(-1,-2)
        @test_throws NSmallerException foo2(1,20)
    end

    @testset "Post-Single" begin
        @not0 function foo(x::Int,y::T) where T<:Number
            x^y
        end
        @not0 a function foo2(x::Int,y::T) where T<:Number
            a = x^y
            return 2*a
        end
        @not0 a function foo3(x::Int,y::T) where T<:Number
            a = x^y
            throw(BeforeException())
            return a
        end
        @test foo(10,10) == 10^10
        @test_throws NZeroException foo(0,10)
        @test foo2(10,10) == 2*10^10
        @test_throws NZeroException foo2(0,10)
        @test_throws BeforeException foo3(0,10)
    end

    @testset "Post-Two" begin
        @smaller x,a function foo(x::Int,a::T) where T<:Number
            return x^a
        end
        @smaller res,con function foo2(x::Int,y::T) where T<:Number
            con = 100
            res = x^y
            return res
        end
        @test foo(10,10) == 10^10
        @test_throws NLargerException foo(2,0)
        @test foo2(5,2) == 5^2
        @test_throws NLargerException foo2(10,10)
    end
    @testset "Pre-Post" begin
        @not0 @larger y,x function foo(x::Int,y::T) where T<:Number
            return x^y
        end
        @test foo(2,3) == 2^3
        @test_throws NSmallerException foo(2,1)
        @test_throws NZeroException foo(0,2)
    end
    @testset "Module" begin
        @test_throws Main.Level1.InModException Level1.foo(1)
        @test_throws Main.Level1.InModException Level1.bar(1)
        Main.Level1.baz(1)
    end
    @testset "No return" begin
        v = [1.0,2.0,3.0]
        @fixLength v function foo!(v)
            push!(v,1.0)
        end
        @fixLength v function baz!(v)
            v[1] = 2.0
        end
        @test baz!(v) == nothing
        @test_throws FixLengthException foo!(v)
    end

    @testset "Multiple Return" begin
        @multret function foo(x,y)
            return x,y
        end
        @test_throws NSmallerException foo(2,1)
        @test foo(1,2) == (1,2)
    end

    @testset "Short Syntax explicit Pre" begin
        @positive x @positive y z function foo(x,y,z)
            return (x,y,z)
        end
        @test foo(1,2,3) == (1,2,3)
        @test_throws NNegativeException foo(-1,2,3)
        @test_throws NNegativeException foo(1,-2,3)
        @test_throws NNegativeException foo(1,2,-3)
    end

    @testset "Short Syntax explicit Post" begin
        @not0 x @not0 y z function foo(x,y,z)
            return (x,y,z)
        end
        @test foo(1,2,3) == (1,2,3)
        @test_throws NZeroException foo(0,2,3)
        @test_throws NZeroException foo(1,0,3)
        @test_throws NZeroException foo(1,2,0)
    end

end


text:
struct ContactRegion{T}
    transform::AffineMap{RotMatrix3{T}, SVector{3, T}}
    μ::T
    μrot::T
    A::Matrix{T}
    b::Vector{T}
end


text:
n, m = parse.(Int, split(readline()))
A = zeros(Int, n, m)
for i in 1:n
    A[i, :] = parse.(Int, split(readline()))
end

v = zeros(Int, m)
for i in 1:m
    v[i] = parse(Int, readline())
end

w = A * v

for i in w
    println(i)
end


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

"""
    UnstructuredMesh(points, connec)

An unstructured mesh with `points` and connectivities
`connec`. The i-th element of the mesh is lazily built
based on the connectivity list `connec[i]`.
"""
struct UnstructuredMesh{Dim,T,Connectivity} <: Mesh{Dim,T}
  points::Vector{Point{Dim,T}}
  connec::Vector{Connectivity}
end

==(m1::UnstructuredMesh, m2::UnstructuredMesh) =
  m1.points == m2.points && m1.connec == m2.connec

vertices(m::UnstructuredMesh) = m.points

function faces(m::UnstructuredMesh{Dim}, r) where {Dim}
  @assert 0 ≤ r ≤ Dim "invalid rank for mesh"
  ps, cs = m.points, m.connec
  r == 0 && return ps
  rank(c) = paramdim(polytopetype(c))
  (materialize(c, ps) for c in cs if rank(c) == r)
end

function Base.show(io::IO, m::UnstructuredMesh{Dim,T}) where {Dim,T}
  nvert = length(m.points)
  nface = length(m.connec)
  print(io, "UnstructuredMesh($nvert vertices, $nface faces)")
end

function Base.show(io::IO, ::MIME"text/plain", m::UnstructuredMesh{Dim,T}) where {Dim,T}
  nvert = length(m.points)
  nface = length(m.connec)
  println(io, "UnstructuredMesh{$Dim,$T}")
  println(io, "  $nvert vertices")
  lines = ["    └─$p" for p in m.points]
  println(io, join(lines, "\n"))
  println(io, "  $nface faces")
  lines = ["    └─$f" for f in m.connec]
  print(  io, join(lines, "\n"))
end