text:
"""
    log_histogram(logger, name, (bins,weights); step=step(logger))

Logs a histogram under the tag `name` on given `step`. The histogram must be
passed as a tuple holding the `N+1` bin edges and the height of the `N` bins.

You can also pass the raw data, and a binning algorithm from `StatsBase.jl` will
be used to bin the data.
"""
function log_histogram(logger::TBLogger, name::AbstractString, (bins,weights)::Tuple{AbstractVector, AbstractArray};
                       step=nothing)
    weights = collect(vec(weights))
    summ = SummaryCollection(histogram_summary(name, collect(bins), weights))
    write_event(logger.file, make_event(logger, summ, step=step))
end

"""
    log_histogram(logger, name, data::Vector; step=step(logger))

Bins the values found in `data` and logs them as an histogram under the tag
`name`.
"""
function log_histogram(logger::TBLogger, name::AbstractString, data::AbstractArray;
                       step=nothing)
    data = vec(data)
    hvals = fit(Histogram, data)
    summ = SummaryCollection(histogram_summary(name, collect(hvals.edges[1]), hvals.weights))
    write_event(logger.file, make_event(logger, summ, step=step))
end

"""
    log_vector(logger, name, data::Vector; step=step(logger))

Logs the vector found in `data` as an histogram under the name `name`.
"""
function log_vector(logger::TBLogger, name::AbstractString, data::AbstractVector; step=nothing)
    summ = SummaryCollection(histogram_summary(name, collect(0:length(data)),data))
    write_event(logger.file, make_event(logger, summ, step=step))
end

function histogram_summary(name::AbstractString, edges::AbstractVector{<:Number}, hist_vals::AbstractVector{<:Number})
    @assert length(edges) == length(hist_vals)+1
    return histogram_summary(name, Histogram(edges, hist_vals))
end

function histogram_summary(name::AbstractString, hist::StatsBase.Histogram{T,1}) where T
    edges     = first(hist.edges)
    hist_vals = hist.weights

    hp = HistogramProto(min=first(edges), max=last(edges),
                        bucket_limit=collect(edges[2:end]),
                        bucket=hist_vals)
    return Summary_Value(tag=name, histo=hp)
end

# Writes to an Histogram summary the flattened version of the array.
# Also stores the shape of the array as a field in a plugin, which allows to
# reconstruct the original shape when read back into Julia
function histogram_arr_summary(name::AbstractString, tensor::AbstractArray)

    smpd = SummaryMetadata_PluginData(plugin_name=TB_PLUGIN_JLARRAY_NAME,
                    content=reinterpret(UInt8, collect(size(tensor))))
    sm = SummaryMetadata(plugin_data=smpd)

    edges = collect(0:length(tensor))
    hp = HistogramProto(min=minimum(edges), max=maximum(edges),
                        bucket_limit=edges[2:end],
                        bucket=vec(tensor))

    return Summary_Value(tag=name, histo=hp, metadata=sm)
end


text:
##=================####==============##

using NetABM
using LightGraphs

using StatsBase
using Statistics
using Bootstrap

using Plots
using LaTeXStrings

using Distributions

##=================####==============##
# SYSYTEM PARAMETERS

params = Params()

params.repetition = 1      # Repetition ID for ensemble statistics
params.num_agents = 100000   # Number of agents (32768)
params.attack_rate = 0.5   # Virus attack rate
params.p_infected_t0 = 0.001# Initial fraction of infected agents

params.σ_cop_agents = 0.1  # Std deviation of cooperation probability
params.μ_cop_agents = 0.25  # Mean of cooperation probability
params.μ_recovT_agents = 15 # Mean recovery time

##=================####==============##
# INITIALIZE MEETS, COOPERATION AND RECOVERY TIMES DISTRIBUTION

pc_μ = params.μ_cop_agents    # MEAN COOPERATION PROBABILITY
pc_σ = params.σ_cop_agents    # STD_DEV COOPERATION PROBABILITY
rt_μ = params.μ_recovT_agents # MEAN RECOVERY TIME

# meets_dist = Exponential()
meets_dist = Geometric()
recovt_dist = truncated(Poisson(rt_μ), 7, 35)

coop_dist = truncated(Normal(pc_μ, pc_σ), 0.0, 1.0)

init_I = params.p_infected_t0
init_no_I = 1.0 - params.p_infected_t0

##=================####==============##

#  g = erdos_renyi(num_agents,0.4)
g = barabasi_albert(params.num_agents, 9)
agents = [Agent(i) for i = 1:(params.num_agents)];

NetABM.mod_init_demographics!(
    agents,
    coop_dist,
    meets_dist,
    recovt_dist;
    states = ["S", "I"],
    initial = [init_no_I, init_I],
    ϵ_thr = [0.25, 0.75]
)
map(x -> assign_contacts!(g, x), agents);

# TIME AT WHICH LOCKDOWN IS ANNOUNCED
#  time_lockdown = 10

tau = []
num_S = []
num_I = []
num_R = []
t_p_cops = Array{Float64}[]

# INITIALIZE TIME AND POPULATIONS
params.now_t = 0
populations = get_populations(agents, params)
global this_I = populations[2]

#  lockdown = false
#  time_lockdown = 0

push!(tau, params.now_t)
push!(num_S, populations[1])
push!(num_I, populations[2])
push!(num_R, populations[3])
push!(t_p_cops, [ag.p_cop for ag in agents])

println(params.now_t, "|", populations[2])
println("SIM STARTS")

while populations[2] > 0.0 # WHILE THERE ARE INFECTED AGENTS

    # RECORD DATA
    push!(tau, params.now_t)
    push!(num_S, populations[1])
    push!(num_I, populations[2])
    push!(num_R, populations[3])
    push!(t_p_cops, [ag.p_cop for ag in agents])

    # FINDS NEXT STATE AND UPDATES ALL AGENTS
    update_all_agents!(agents, params)

    # OPINION DYNAMICS
    for ag_id = 1:(params.num_agents)
        NetABM.HK_od_model!(agents, ag_id)
    end
    NetABM.set_athome_agents!(agents)

    #  EACH AGENTS DECIDES TO STAY AT HOME/ NOT TO STAY AT HOME
    #  if params.now_t >= time_lockdown
    #  #  if lockdown
    #      NetABM.set_athome_agents!(agents)
    #      # OPINION DYNAMICS
    #      for ag_id = 1:(params.num_agents)
    #          NetABM.HK_od_model!(agents, ag_id)
    #      end
    #  end

    # GET CURRENT POPULATIONS AND UPDATE TIME
    populations = get_populations(agents, params)
    params.now_t += 1

    prev_I = this_I
    this_I = populations[2]

    #  if !lockdown
    #      if this_I / prev_I >= 2.0
    #          lockdown = true
    #          time_lockdown = params.now_t
    #      end
    #  end

    #  println(params.now_t, "|", populations[2], "|", this_I / prev_I)
    println(params.now_t, populations)
end

println("NO MORE INFECTED!")

##=================####==============##

plt = plot()
plot!(plt, tau, num_I ./ params.num_agents, lab = "I")
#  plot!(plt, tau, num_S ./ params.num_agents, lab = "S")
#  plot!(plt, tau, num_R ./ params.num_agents, lab = "R")
#  vline!(plt, [time_lockdown], color = :red, ls = :dash, lab = "")

##=================####==============##

od_t = transpose(hcat(t_p_cops...));
plot(od_t, alpha = 0.45, lab = "")

##=================####==============##

ext: jl


text:
# IntervalSets.jl currently has deprecation warnings on Julia 0.7

export Interval, ..

abstract type AbstractInterval{T}; end

"""
    Interval(start, stop)

Construct a closed interval, which is a collection which contains (via `in`) elements
between `start` and `stop` (inclusive) according to `isless`. The collection is abstract
in nature and doesn't support iteration, indexing, etc.

Can be constructed via the `..` function, e.g. `1..3 === Interval(1, 3)`.

# Examples

```julia
julia> 2 in Interval(1, 3)
true

julia> 3 in Interval(1, 3)
true

julia> 4 in Interval(1, 3)
false
"""
struct Interval{T} <: AbstractInterval{T}
    start::T
    stop::T
end

"""
    ..(start, stop)
    start..stop

Constructs an `Interval(start, stop)`, which represents the closed interval between `start`
and `stop`. `Interval`s are abstract collections which support `in` but not iteration,
indexing, etc.

# Examples

```julia
julia> 2 in 1..3
true

julia> 3 in 1..3
true

julia> 4 in 1..3
false
"""
..(start::T, stop::T) where {T} = Interval{T}(start, stop)
..(start, stop) = ..(promote(start, stop)...)

function in(x::T, interval::Interval{T}) where {T}
    return !isless(x, interval.start) && !isless(interval.stop, x)
end

function in(x, interval::Interval{T}) where {T}
    return !isless(promote(x, interval.start)...) && !isless(promote(interval.stop, x)...)
end

function Base.:(==)(i1::Interval, i2::Interval)
    (isequal(i1.start, i2.start) && isequal(i1.stop, i2.stop)) || isless(i1.stop, i2.start)
end

function Base.isequal(i1::Interval, i2::Interval)
    isequal(i1.start, i2.start) && isequal(i1.stop, i2.stop)
end

function Base.isless(i1::Interval, i2::Interval)
    isless(i1.start, i2.start) || (isequal(i1.start, i2.start) && isless(i1.start, i2.start))
end

Base.hash(interval::Interval, h) = hash(hash(interval.start, h), interval.stop)

function show(io::IO, interval::Interval)
    print(io, interval.start)
    print(io, "..")
    print(io, interval.stop)
end


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
Unit tests for the One type.
"""

# --- Imports

# Standard library
import InteractiveUtils.subtypes
using Test

# GeometricAlgebra.jl
using GeometricAlgebra

# --- Constructor tests

@testset "One: constructor tests" begin
    # --- Inner constructor

    for precision_type in subtypes(AbstractFloat)
        B = One{precision_type}()
        @test B isa One{precision_type}
        @test B === One{precision_type}()
    end

    # --- Outer constructor

    B = One()
    @test B isa One{Float64}
end

@testset "One: one()" begin

    # --- one(M::AbstractMultivector)

    for precision_type in subtypes(AbstractFloat)
        # --- M::AbstractScalar

        # M::Zero
        B = one(Zero{precision_type}())
        @test B isa One{precision_type}

        # M::One
        B = one(One{precision_type}())
        @test B isa One{precision_type}

        # M::Scalar
        B = one(Scalar{precision_type}(5))
        @test B isa One{precision_type}

        # --- M::AbstractBlade

        # M::Blade
        B = one(Blade{precision_type}([1 2 3]))
        @test B isa One{precision_type}

        # M::Pseudoscalar
        B = one(Pseudoscalar{precision_type}(10, 5))
        @test B isa One{precision_type}

        # --- M::AbstractMultivector

        # M::Multivector
        B = one(Multivector{precision_type}([Scalar(3)]))
        @test B isa One{precision_type}
    end

    # --- one(M::Type{<:AbstractMultivector})

    B = one(AbstractScalar)
    @test B isa One{Float64}

    B = one(Zero)
    @test B isa One{Float64}

    B = one(One)
    @test B isa One{Float64}

    B = one(Scalar)
    @test B isa One{Float64}

    B = one(AbstractBlade)
    @test B isa One{Float64}

    B = one(Blade)
    @test B isa One{Float64}

    B = one(Pseudoscalar)
    @test B isa One{Float64}

    B = one(AbstractMultivector)
    @test B isa One{Float64}

    B = one(Multivector)
    @test B isa One{Float64}

    # --- one(::Type{<:AbstractMultivector{T}})

    for precision_type in subtypes(AbstractFloat)
        # --- Type{<:AbstractScalar{T}}

        B = one(AbstractScalar{precision_type})
        @test B isa One{precision_type}

        B = one(Zero{precision_type})
        @test B isa One{precision_type}

        B = one(One{precision_type})
        @test B isa One{precision_type}

        B = one(Scalar{precision_type})
        @test B isa One{precision_type}

        # --- Type{<:AbstractBlade{T}}

        B = one(AbstractBlade{precision_type})
        @test B isa One{precision_type}

        B = one(Blade{precision_type})
        @test B isa One{precision_type}

        B = one(Pseudoscalar{precision_type})
        @test B isa One{precision_type}

        # --- Type{<:AbstractMultivector{T}}

        B = one(AbstractMultivector{precision_type})
        @test B isa One{precision_type}

        B = one(Multivector{precision_type})
        @test B isa One{precision_type}
    end
end

# --- Test comparison operations

@testset "One: isone(B)" begin
    # Basic functions
    for precision_type in subtypes(AbstractFloat)
        # Blade
        vectors = [3 3; 4 4; 0 1; 0 1]
        B = Blade{precision_type}(vectors)
        @test !isone(B)

        # Pseudoscalar
        B = Pseudoscalar{precision_type}(10, 3)
        @test !isone(B)

        # Scalar
        B = Scalar{precision_type}(1)
        @test isone(B)

        B = Scalar{precision_type}(3)
        @test !isone(B)

        # Zero
        B = Zero{precision_type}()
        @test !isone(B)

        # One
        B = One{precision_type}()
        @test isone(B)
    end
end

# --- Test attribute methods

@testset "One: AbstractMultivector attribute functions" begin
    # Basic functions
    for precision_type in subtypes(AbstractFloat)
        B = One{precision_type}()
        @test dim(B) == 0
        @test grades(B) == [0]
        @test blades(B) == [B]
        @test norm(B) isa precision_type
        @test norm(B) == 1

        @test B[0] == [B]
        @test B[1] == []
    end
end

@testset "One: AbstractBlade attribute functions" begin
    # Basic functions
    for precision_type in subtypes(AbstractFloat)
        B = One{precision_type}()
        @test grade(B) == 0
        @test basis(B) == 1
        @test volume(B) isa precision_type
        @test volume(B) == 1
        @test sign(B) == 1
    end
end

@testset "One: AbstractScalar attribute functions" begin
    for precision_type in subtypes(AbstractFloat)
        B = One{precision_type}()
        @test value(B) isa precision_type
        @test value(B) == 1
    end
end

# --- Tests for AbstractMultivector interface functions

@testset "One: -(B)" begin
    for precision_type in subtypes(AbstractFloat)
        B = One{precision_type}()

        negative_B = -B
        @test negative_B isa Scalar{precision_type}
        @test negative_B == Scalar{precision_type}(-1)
    end
end

@testset "One: reverse(B)" begin
    for precision_type in subtypes(AbstractFloat)
        B = One{precision_type}()
        reverse_B = reverse(B)
        @test reverse_B === B
        @test B * reverse_B == norm(B)^2
    end
end

@testset "One: dual(B)" begin
    for precision_type in subtypes(AbstractFloat)
        B = One{precision_type}()

        for test_dim in 5:8
            dual_B = dual(B, test_dim)

            expected_dual = mod(test_dim, 4) < 2 ?
                Pseudoscalar{precision_type}(test_dim, 1) :
                Pseudoscalar{precision_type}(test_dim, -1)
            @test dual_B isa Pseudoscalar{precision_type}
            @test dual_B == expected_dual
        end

        expected_message = "The dual of a scalar is not well-defined if " *
                           "`dim` is not specified"
        @test_throws ErrorException(expected_message) dual(B)
    end
end

@testset "One: convert(B)" begin
    for precision_type_converted in subtypes(AbstractFloat)
        for precision_type_src in subtypes(AbstractFloat)
            B = One{precision_type_src}()
            B_converted = convert(precision_type_converted, B)
            @test B_converted isa One{precision_type_converted}
        end
    end
end

# --- Tests for AbstractBlade interface functions

@testset "One: inv(B)" begin
    for precision_type in subtypes(AbstractFloat)
        B = One{precision_type}()

        inverse_B = inv(B)
        @test inverse_B === B
    end
end


text:
using ChromNet
using JLD
using Base.Test
srand(1)

# create a sample network dataset to use
jldopen("data/small.ChromNet.jld", "w") do f
    data = zeros(Int32, 4, 1000000)
    data[:,:] = rand(0:100, 4, 1000000)
    ids = ["A01", "B01", "C01", "D01"]
    f["P"] = 4
    f["N"] = 1000000
    f["ids"] = ids
    f["data"] = data
    f["means"] = mean(data, 2)
    f["norms"] = Float64[norm(data[i,:] - mean(data[i,:])) for i in 1:4]
    f["metadata"] = Dict(
        "A01" => Dict(
            "name" => "A",
            "id" => "A01",
            "description" => "A desc",
            "cellType" => "A cellType",
            "lab" => "A cellType",
            "organism" => "Homo sapiens",
            "lifeStage" => "None",
            "treatments" => "None",
            "antibody" => "Unknown"
        ),
        "B01" => Dict(
            "name" => "B",
            "id" => "B01",
            "description" => "B desc",
            "cellType" => "B cellType",
            "lab" => "B cellType",
            "organism" => "Homo sapiens",
            "lifeStage" => "None",
            "treatments" => "None",
            "antibody" => "Unknown"
        ),
        "C01" => Dict(
            "name" => "C",
            "id" => "C01",
            "description" => "C desc",
            "cellType" => "C cellType",
            "lab" => "C cellType",
            "organism" => "Homo sapiens",
            "lifeStage" => "None",
            "treatments" => "None",
            "antibody" => "Unknown"
        ),
        "D01" => Dict(
            "name" => "D",
            "id" => "D01",
            "description" => "D desc",
            "cellType" => "D cellType",
            "lab" => "D cellType",
            "organism" => "Homo sapiens",
            "lifeStage" => "None",
            "treatments" => "None",
            "antibody" => "Unknown"
        )
    )
end

include("util.jl")
include("groupgm.jl")
include("scripts.jl")


text:
using BinaryBuilder, Pkg

julia_version = v"1.5.3"

name = "HelFEM"
version = v"0.1.5"
sources = [
    GitSource("https://github.com/mortenpi/HelFEM.git", "a4d3b2e6f16f7f7953afd5f69a44257a65c5b131"),
]

script = raw"""
cp -v ${WORKSPACE}/srcdir/HelFEM/julia/CMake.system ${WORKSPACE}/srcdir/HelFEM/CMake.system

# Set up some platform specific CMake configuration. This is more or less borrowed from
# the Armadillo build_tarballs.jl script:
#   https://github.com/JuliaPackaging/Yggdrasil/blob/48d7a89b4aa46b1a8c91269bb138a660f4ee4ece/A/armadillo/build_tarballs.jl#L23-L52
#
# We need to manually set up OpenBLAS because FindOpenBLAS.cmake does not work with BB:
if [[ "${nbits}" == 64 ]] && [[ "${target}" != aarch64* ]]; then
    OPENBLAS="${libdir}/libopenblas64_.${dlext}"
else
    OPENBLAS="${libdir}/libopenblas.${dlext}"
fi

# Compile libhelfem as a static library
cd ${WORKSPACE}/srcdir/HelFEM/
cmake \
    -DCMAKE_BUILD_TYPE=Release \
    -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_POSITION_INDEPENDENT_CODE=ON \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DHELFEM_BINARIES=OFF \
    -DHELFEM_FIND_DEPS=ON \
    -DUSE_OPENMP=OFF \
    -B build/ -S .
make -C build/ -j${nproc} helfem
make -C build/ install
# Copy the HelFEM license
install_license LICENSE

# Compile the CxxWrap wrapper as a shared library
cd ${WORKSPACE}/srcdir/HelFEM/julia
cmake \
    -DBLAS_LIBRARIES=${OPENBLAS} \
    -DCMAKE_BUILD_TYPE=Release \
    -DCMAKE_FIND_ROOT_PATH=${prefix} \
    -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DJulia_PREFIX=${prefix} \
    -B build/ -S .
make -C build/ -j${nproc}
make -C build/ install
"""

# These are the platforms the libcxxwrap_julia_jll is built on.
include("../../L/libjulia/common.jl")
platforms = libjulia_platforms(julia_version)
platforms = expand_cxxstring_abis(platforms)

products = [
    LibraryProduct("libhelfem-cxxwrap", :libhelfem),
]

dependencies = [
    BuildDependency(PackageSpec(name="libjulia_jll", version=julia_version)),
    Dependency("libcxxwrap_julia_jll"),
    Dependency(PackageSpec(name = "armadillo_jll", version = "9.850.1")),
    Dependency(PackageSpec(name = "GSL_jll", version = "2.6.0")),
    Dependency(PackageSpec(name = "OpenBLAS_jll", version = "0.3.9")),
]

# preferred_gcc_version = v"8" is a requirement from libcxxwrap_julia_jll
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version=v"8",
               julia_compat = "$(julia_version.major).$(julia_version.minor)")

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
using Rethinking
using Test
using Turing
using DataFrames
using PrettyTables

@testset "univariate quap" begin
    @model function test(y)
        μ ~ Normal(10, 1)
        y .~ Normal(μ, 1)
    end
    y = randn(1000) .+ 10
    res = quap(test(y))
    μ = res.coef.μ
    @test round(μ) == 10
end

@testset "multivariate quap" begin
    @model function test(y)
        μ ~ Normal(10, 1)
        σ ~ Exponential(5)
        y .~ Normal(μ, σ)
    end
    y = randn(1000)*5 .+ 10
    res = quap(test(y))
    μ, σ = res.coef.μ, res.coef.σ
    @test round(μ) == 10
    @test round(σ) == 5
end

@testset "precise" begin
    df = DataFrame((a=rand(1000)))
    @test isnothing(precise(df))
    """
    ┌───────┬───────────────────────────────────────────┐
    │ param │ mean   std  5.5%   50%  94.5%   histogram │
    ├───────┼───────────────────────────────────────────┤
    │     a │ 0.51  0.28  0.06  0.52   0.95  ▆▇▇▇▇▇█▇▆▇ │
    └───────┴───────────────────────────────────────────┘
    """
end


text:
## Función de inflación subyacente MAI (muestra ampliada implícitamente)

## Grilla de variaciones intermensuales
const V = range(-200, 200, step=0.01f0) # -100:0.01:100

## Algoritmos de cómputo MAI

abstract type AbstractMaiMethod end

# Algoritmo de cómputo de MAI-G con n segmentos marcados en las posiciones p
"""
    MaiG{P} <: AbstractMaiMethod

    MaiG(n::Int)
    MaiG(p::AbstractArray)

Tipo para englobar la metodología de cómputo de inflación subyacente MAI-G, en
la cual se transforma la distribución ponderada de variaciones intermensuales
utilizando la distribución histórica de variaciones intermensuales ponderadas.

Se proporciona el número de segmentos de normalización, o bien, las posiciones
de los cuantiles utilizados para llevar a cabo la transformación. Si se
proporcionan las posiciones, la primera y la última posición deben ser el
cuantil 0 y 1, respectivamente.

## Ejemplos

1. Utilizar los quintiles como puntos de referencia para normalización
```jldoctest
julia> method = MaiG(5)
MaiG{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}(5, 0.0:0.2:1.0)
```

2. Utilizar tres segmentos de normalización, en el primer y tercer cuartil:
```jldoctest
julia> method = MaiG([0, 0.25, 0.75, 1])
MaiG{Vector{Float64}}(3, [0.0, 0.25, 0.75, 1.0])
```
"""
struct MaiG{P} <: AbstractMaiMethod
    n::Int
    p::P

    function MaiG(n, p)
        _checkmethod(n, p)
        new{typeof(p)}(n, p)
    end
end

# Algoritmo de cómputo de MAI-F con n segmentos marcados en las posiciones p
"""
    MaiF{P} <: AbstractMaiMethod

    MaiF(n::Int)
    MaiF(p::AbstractArray)

Tipo para englobar la metodología de cómputo de inflación subyacente MAI-F, en
la cual se transforma la distribución **de ocurrencias** (o equiponderada) de
variaciones intermensuales utilizando la distribución histórica de variaciones
intermensuales ponderadas.

Se proporciona el número de segmentos de normalización, o bien, las posiciones
de los cuantiles utilizados para llevar a cabo la transformación. Si se
proporcionan las posiciones, la primera y la última posición deben ser el
cuantil 0 y 1, respectivamente.

## Ejemplos

1. Utilizar los quintiles como puntos de referencia para normalización
```jldoctest
julia> method = MaiF(5)
MaiF{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}(5, 0.0:0.2:1.0)
```

2. Utilizar tres segmentos de normalización, en el primer y tercer cuartil:
```jldoctest
julia> method = MaiF([0, 0.25, 0.75, 1])
MaiF{Vector{Float64}}(3, [0.0, 0.25, 0.75, 1.0])
```
"""
struct MaiF{P} <: AbstractMaiMethod
    n::Int
    p::P

    function MaiF(n, p)
        _checkmethod(n, p)
        new{typeof(p)}(n, p)
    end
end

"""
    MaiFP{P} <: AbstractMaiMethod

    MaiFP(n::Int)
    MaiFP(p::AbstractArray)

Tipo para englobar la metodología de cómputo de inflación subyacente MAI-FP, en
la cual se transforma la distribución **de ocurrencias** de variaciones
intermensuales utilizando la distribución histórica de variaciones
intermensuales equiponderadas. Esta es la versión equivalente a la MAI-G, que
utiliza **todas las distribuciones de ocurrencias**.

Se proporciona el número de segmentos de normalización, o bien, las posiciones
de los cuantiles utilizados para llevar a cabo la transformación. Si se
proporcionan las posiciones, la primera y la última posición deben ser el
cuantil 0 y 1, respectivamente.

## Ejemplos

1. Utilizar los quintiles como puntos de referencia para normalización
```jldoctest
julia> method = MaiFP(5)
MaiFP{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}(5, 0.0:0.2:1.0)
```

2. Utilizar tres segmentos de normalización, en el primer y tercer cuartil:
```jldoctest
julia> method = MaiFP([0, 0.25, 0.75, 1])
MaiFP{Vector{Float64}}(3, [0.0, 0.25, 0.75, 1.0])
```
"""
struct MaiFP{P} <: AbstractMaiMethod
    n::Int
    p::P

    function MaiFP(n, p)
        _checkmethod(n, p)
        new{typeof(p)}(n, p)
    end
end

# Revisión de condiciones para aplicación de métodos MAI
function _checkmethod(n, p)
    length(p) == n+1 || error("Distribución de percentiles de tamaño incorrecto")
    issorted(p) || error("Distribución de percentiles debe estar ordenada")
    all(0 .<= p .<= 1) || error("Cuantiles deben estar entre cero y uno")
    (first(p) == 0 && last(p) == 1) || error("Primer y último cuantil deben ser 0, 1")
end

# Constructores para división equitativa de posiciones p
MaiF(n::Int) = MaiF(n, (0:n)/n)
MaiG(n::Int) = MaiG(n, (0:n)/n)
MaiFP(n::Int) = MaiFP(n, (0:n)/n)
MaiF(p::AbstractArray) = MaiF(length(p)-1, p)
MaiG(p::AbstractArray) = MaiG(length(p)-1, p)
MaiFP(p::AbstractArray) = MaiFP(length(p)-1, p)

function Base.string(method::AbstractMaiMethod)
    algorithm = method isa MaiG ? "G" : (method isa MaiF ? "F" : "FP")
    if method.p isa StepRangeLen
        return "($algorithm," * string(method.n) * ")"
    else
        return "($algorithm," * string(method.n) * "," * string(round.(method.p[2:end-1], digits=2)) * ")"
    end
end

## Definición de la función de inflación
"""
    InflationCoreMai{T <: AbstractFloat, B, M <:AbstractMaiMethod}
        <: InflationFunction

    InflationCoreMai(vspace::StepRangeLen{T, B, B}, method::M)

Función de inflación de la metodología de muestra ampliada implicítamente (MAI).
Se parametriza en el tipo flotante `T` de los datos, `B` el tipo para
representar la precisión de la grilla `vspace` en el algoritmo de cómputo y el
método `M` de cómputo.

Los métodos de cómputo disponibles son:
- Metodología MAI-G: se transforma la distribución **ponderada** de variaciones
  intermensuales utilizando la distribución histórica de variaciones
  intermensuales ponderadas. Se debe dar como argumento el método
  [`MaiG`](@ref).
- Metodología MAI-F: se transforma la distribución **de ocurrencias** de
  variaciones intermensuales utilizando la distribución histórica de variaciones
  intermensuales ponderadas. Se debe dar como argumento el método
  [`MaiF`](@ref).
- Metodología MAI-FP: se transforma la distribución **de ocurrencias** de
  variaciones intermensuales utilizando la distribución histórica de variaciones
  intermensuales **equiponderadas**. Es un equivalente directo de la metodología
  MAI-G, reemplazando todas las distribuciones por las versiones de ocurrencias.
  Se debe dar como argumento el método [`MaiFP`](@ref).
"""
Base.@kwdef struct InflationCoreMai{T <: AbstractFloat, B, M <:AbstractMaiMethod} <: InflationFunction
    vspace::StepRangeLen{T, B, B} = V
    method::M = MaiG(4)
end

# Constructor de conveniencia para especificar V por defecto
InflationCoreMai(method::AbstractMaiMethod) = InflationCoreMai(V, method)

# Nombre de la medida
measure_name(inflfn::InflationCoreMai) = "MAI " * string(inflfn.method)

# Parámetros
CPIDataBase.params(inflfn::InflationCoreMai) = (inflfn.method, )

## Métodos de cómputo intermensual

# Operación sobre CountryStructure para obtener variaciones intermensuales de la
# estructura de país. Esta función llama a los métodos que reciben el método de
# cómputo MAI
function (inflfn::InflationCoreMai)(cs::CountryStructure, ::CPIVarInterm)
    inflfn(cs, CPIVarInterm(), inflfn.method)
end

# Función de resumen intermensual para metodología MAI-G
function (inflfn::InflationCoreMai{T})(cs::CountryStructure, ::CPIVarInterm, method::MaiG) where T
    # Computar flp y glp, tomando en cuenta observaciones de años completos en
    # la última base del CountryStructure
    V_star = _get_vstar(cs)
    W_star = _get_wstar(cs)
    glp = WeightsDistr(V_star, W_star, inflfn.vspace)

    # Obtener distribuciones acumuladas y sus percentiles
    GLP = cumsum(glp)
    q_glp::Vector{T} = quantile(GLP, method.p)

    # Llamar al método de cómputo de inflación intermensual
    vm_fn = base -> inflfn(base, inflfn.method, glp, GLP, q_glp)
    vm = mapfoldl(vm_fn, vcat, cs.base)
    vm
end

# Función de resumen intermensual para metodología MAI-F
function (inflfn::InflationCoreMai{T})(cs::CountryStructure, ::CPIVarInterm, method::MaiF) where T
    # Intuitivamente, las distribuciones de largo plazo podrían computarse más
    # sencillamente de esta forma. Sin embargo, parece que hay problemas de
    # precisión en los vectores dispersos al agregar las distribuciones de cada
    # base de esta manera.

    # flp_bases = ObservationsDistr.(cs.base, Ref(vspace))
    # glp_bases = WeightsDistr.(cs.base, Ref(vspace))
    # flp = sum(flp_bases)
    # glp = sum(glp_bases)

    # Computar flp y glp, tomando en cuenta observaciones de años completos en
    # la última base del CountryStructure
    V_star = _get_vstar(cs)
    W_star = _get_wstar(cs)
    flp = ObservationsDistr(V_star, inflfn.vspace)
    glp = WeightsDistr(V_star, W_star, inflfn.vspace)

    # Obtener distribuciones acumuladas y sus percentiles
    FLP = cumsum(flp)
    GLP = cumsum(glp)
    # q_glp::Vector{T} = quantile(GLP, method.p)
    q_flp::Vector{T} = quantile(FLP, method.p)

    # Llamar al método de cómputo de inflación intermensual
    vm_fn = base -> inflfn(base, inflfn.method, glp, GLP, q_flp)
    vm = mapfoldl(vm_fn, vcat, cs.base)
    vm
end

# Función de apoyo para obtener V_star, la ventana histórica con todas las
# variaciones intermensuales. Utilizada para cómputos de distribuciones de largo
# plazo
function _get_vstar(cs::CountryStructure)
    # Revisar paquete CatViews para ahorrar un poco más de memoria, to-do...
    lastbase = cs.base[end]
    T_lp = periods(lastbase) ÷ 12
    v_last = view(lastbase.v[1 : 12*T_lp, :], :)
    v_first = map(base -> view(base.v, :), cs.base[1:end-1])
    V_star = vcat(v_first..., v_last)
    V_star
end

# Función de apoyo para obtener W_star, vector de ponderaciones asociado a las
# variaciones intermensuales históricas
function _get_wstar(cs::CountryStructure)
    # Ponderaciones de toda la base
    lastbase = cs.base[end]
    T_lp = periods(lastbase) ÷ 12
    w_first = map(base -> view(repeat(base.w', periods(base)), :), cs.base[1:end-1])
    w_last = view(repeat(lastbase.w', 12*T_lp), :)
    W_star = vcat(w_first..., w_last)
    W_star
end

## Métodos de cómputo MAI sobre VarCPIBase

# Se utiliza información de largo plazo provista por el método que opera sobre
# CountryStructure y CPIVarInterm

# Variaciones intermensuales resumen con método de MAI-G
function (inflfn::InflationCoreMai)(base::VarCPIBase{T}, method::MaiG, glp, GLP, q_glp) where T

    mai_m = Vector{T}(undef, periods(base))
    q_g_list = [zeros(T, method.n+1) for _ in 1:Threads.nthreads()]

    # Utilizar la glp y la GLP para computar el resumen intermensual por
    # metodología de inflación subyacente MAI-G
    Threads.@threads for t in 1:periods(base)

        # Obtener lista de percentiles para el hilo
        j = Threads.threadid()
        q_g = q_g_list[j]

        # Computar distribución g y acumularla
        g = WeightsDistr((@view base.v[t, :]), base.w, inflfn.vspace)
        g_acum = cumsum!(g)

        # Computar percentiles de distribución g
        quantile!(q_g, g_acum, method.p)

        # Computar resumen intermensual basado en glpₜ
        mai_m[t] = renorm_g_glp_perf(g_acum, GLP, glp, q_g, q_glp, method.n)
    end

    mai_m
end

# Variaciones intermensuales resumen con método de MAI-F
function (inflfn::InflationCoreMai)(base::VarCPIBase{T}, method::MaiF, glp, GLP, q_flp) where T

    mai_m = Vector{T}(undef, periods(base))
    q_f_list = [zeros(T, method.n+1) for _ in 1:Threads.nthreads()]

    # Utilizar la glp y (FLP, GLP) para computar el resumen intermensual por
    # metodología de inflación subyacente MAI-F
    Threads.@threads for t in 1:periods(base)

        # Obtener lista de percentiles para el hilo
        j = Threads.threadid()
        q_f = q_f_list[j]

        # Computar distribución f y acumularla
        f = ObservationsDistr((@view base.v[t, :]), inflfn.vspace)
        f_acum = cumsum!(f)

        # Computar percentiles de distribución f
        quantile!(q_f, f_acum, method.p)

        # Computar resumen intermensual basado en flpₜ
        mai_m[t] = renorm_f_flp_perf(f_acum, GLP, glp, q_f, q_flp, method.n)
    end

    mai_m
end

## Metodología MAI-FP

# Aplicación directa de las fórmulas de normalización de la
# metodología MAI-G, utilizando las distribuciones de ocurrencias

# Función de resumen intermensual para metodología MAI-FP
function (inflfn::InflationCoreMai{T})(cs::CountryStructure, ::CPIVarInterm, method::MaiFP) where T
    # Computar flp, tomando en cuenta observaciones de años completos en
    # la última base del CountryStructure
    V_star = _get_vstar(cs)
    flp = ObservationsDistr(V_star, inflfn.vspace)

    # Obtener distribuciones acumuladas y sus percentiles
    FLP = cumsum(flp)
    q_flp::Vector{T} = quantile(FLP, method.p)

    # Llamar al método de cómputo de inflación intermensual
    vm_fn = base -> inflfn(base, inflfn.method, flp, FLP, q_flp)
    vm = mapfoldl(vm_fn, vcat, cs.base)
    vm
end

# Variaciones intermensuales resumen con método de MAI-FP
function (inflfn::InflationCoreMai)(base::VarCPIBase{T}, method::MaiFP, flp, FLP, q_flp) where T

    mai_m = Vector{T}(undef, periods(base))
    q_f_list = [zeros(T, method.n+1) for _ in 1:Threads.nthreads()]

    # Utilizar la flp y la FLP para computar el resumen intermensual por
    # metodología de inflación subyacente MAI-G
    Threads.@threads for t in 1:periods(base)

        # Obtener lista de percentiles para el hilo
        j = Threads.threadid()
        q_f = q_f_list[j]

        # Computar distribución f y acumularla
        f = ObservationsDistr((@view base.v[t, :]), inflfn.vspace)
        f_acum = cumsum!(f)

        # Computar percentiles de distribución f
        quantile!(q_f, f_acum, method.p)

        # Computar resumen intermensual basado en flpₜ
        mai_m[t] = renorm_g_glp_perf(f_acum, FLP, flp, q_f, q_flp, method.n)
    end

    mai_m
end


text:
using OrdinaryDiffEq

f = function (du,u,p,t)
  du[1][1] = 0.85u[1][1]
end

u = [[1]]

prob = ODEProblem(f,u,(0.0,1.0))

sol = solve(prob,dt=1/2^4)


text:
using Dijkstra
using Test

graph = Graph{Symbol}()
add_edge!(graph, :S, :A, 4)
add_edge!(graph, :S, :B, 3)
add_edge!(graph, :S, :D, 7)
add_edge!(graph, :A, :C, 1)
add_edge!(graph, :B, :S, 3)
add_edge!(graph, :B, :D, 4)
add_edge!(graph, :C, :E, 1)
add_edge!(graph, :C, :D, 3)
add_edge!(graph, :D, :E, 1)
add_edge!(graph, :D, :T, 3)
add_edge!(graph, :D, :F, 5)
add_edge!(graph, :E, :G, 2)
add_edge!(graph, :G, :E, 2)
add_edge!(graph, :G, :T, 3)
add_edge!(graph, :T, :F, 5)

solution = ShortestPath(graph, :S)
@test get_distance(solution, :T) == 10


text:
# Deterministic GSGP
function dgsgp(cfg, X_train, y_train, X_test, y_test; rec=Record(), seed=1, verbose=true)
    Random.seed!(seed)
    init_id()
    init_ref_counter(cfg)
    init_entry_pool()

    data_train = Dataset(X_train, y_train)
    data_test  = Dataset(X_test, y_test)

    init_entries = cfg.create_population(cfg, data_train, data_test)
    add_entry!.(init_entries)

    pop = getfield.(init_entries, :index)
    new_pop = similar(pop)
    old_pop = similar(pop)

    avoid_overwrite.(pop)

    fitness = getfield.(init_entries, :train_fitness)

    register!(rec, "train_fitness", train_fitness, Float64)
    register!(rec, "test_fitness", test_fitness, Float64)
    add_record!(rec, 1, pop, cfg)

    verbose && print_latest_fitness(rec, 1)

    overwrite = false

    for g in 2:cfg.num_generations
        overwrite = g > 2
        elite_idx = argmin(fitness)
        new_pop[1] = pop[elite_idx]
        avoid_overwrite(pop[elite_idx])
        for i in 2:cfg.population_size
            rand_num = rand()

            if rand_num < cfg.p_crossover
                p1 = tournament_selection(pop, fitness, cfg)
                p2 = tournament_selection(pop, fitness, cfg)
                offspring = deterministic_crossover(p1, p2, data_train, data_test, cfg, overwrite)
                new_pop[i] = offspring.index
            else
                p = tournament_selection(pop, fitness, cfg)
                offspring = gsmutation(p, data_train, data_test, cfg, overwrite)
                new_pop[i] = offspring.index
            end
            # TODO: reproduction
        end

        update_fitness!(fitness, new_pop)
        add_record!(rec, 1, pop, cfg)

        g > 2 && update_entry_pool!(old_pop)

        old_pop = copy(pop)
        pop = copy(new_pop)

        verbose && print_latest_fitness(rec, g)
    end

    return pop[argmin(fitness)], rec
end

function deterministic_crossover(t1, t2, data_train, data_test, cfg, overwrite)
    increment_ref(entry_pool[t1])
    increment_ref(entry_pool[t2])

    local r, offspring
    if overwrite
        offspring = deterministic_crossover_entry(entry_pool[t1], entry_pool[t2], data_train, data_test, cfg)
    else
        offspring = deterministic_crossover_entry(entry_pool[t1], entry_pool[t2], data_train, data_test, cfg, index=get_counter_size()+1)
    end

    add_entry!(offspring)
    set_ref(offspring)

    # to avoid overwritten
    # need to decrement at the end of current generation
    increment_ref(offspring)

    return offspring
end

function deterministic_crossover_entry(t1, t2, data_train, data_test, cfg; index=get_usable_index())

    r = deterministic_coeff(t1.train_semantics, t2.train_semantics, data_train.target)

    train_semantics = Vector{Float64}(undef, length(data_train.target))
    @inbounds @simd for i in 1:length(train_semantics)
        train_semantics[i] = r * t1.train_semantics[i] + (1-r) * t2.train_semantics[i]
    end

    test_semantics = Vector{Float64}(undef, length(data_test.target))
    @inbounds @simd for i in 1:length(test_semantics)
        test_semantics[i] = r * t1.test_semantics[i] + (1-r) * t2.test_semantics[i]
    end

    train_fitness = cfg.fitness(train_semantics, data_train)
    test_fitness = cfg.fitness(test_semantics, data_test)

    id = get_new_id()
    #index = get_usable_index()

    return Entry(id, index, crossover, dummy_node(), (t1.id, t2.id, -1), (t1.index, t2.index, -1), train_semantics, test_semantics, train_fitness, test_fitness)
end

function dgsgp_opt_mate(cfg, data_train, data_test; seed=1, verbose=true)
    Random.seed!(seed)
    init_id()
    init_ref_counter(cfg)
    init_entry_pool()

    init_entries = cfg.create_population(cfg, data_train, data_test)
    add_entry!.(init_entries)

    pop = getfield.(init_entries, :index)
    new_pop = similar(pop)
    old_pop = similar(pop)

    fitness = getfield.(init_entries, :train_fitness)
    rec = Record()
    register!(rec, "train_fitness", train_fitness, Float64)
    register!(rec, "test_fitness", test_fitness, Float64)
    add_record!(rec, 1, pop, cfg)

    verbose && println(1, ": ", minimum(fitness))
    overwrite = false

    for g in 2:cfg.num_generations
        overwrite = g != 2
        elite_idx = argmin(fitness)
        new_pop[1] = pop[elite_idx]
        avoid_overwrite(elite_idx)
        for i in 2:cfg.population_size
            rand_num = rand()

            if rand_num < cfg.p_crossover
                p1 = tournament_selection(pop, fitness, cfg)
                local p2
                if rand() < cfg.p_mate
                    p2 = optimal_mate(p1, data_train.target, pop)
                else
                    p2 = tournament_selection(pop, fitness, cfg)
                end
                offspring = deterministic_crossover(p1, p2, data_train, data_test, cfg, overwrite)
                new_pop[i] = offspring.index
            else
                p = tournament_selection(pop, fitness, cfg)
                offspring = gsmutation(p, data_train, data_test, cfg, overwrite)
                new_pop[i] = offspring.index
            end
        end

        update_fitness!(fitness, new_pop)
        add_record!(rec, g, new_pop, cfg)

        g > 2 && update_entry_pool!(old_pop)

        old_pop = copy(pop)
        pop = copy(new_pop)

        verbose && println(g, ": ", minimum(fitness))
    end

    return pop[argmin(fitness)], rec
end

deterministic_coeff(t1, t2, p) = ifelse(t1 ≈ t2, 0.5, (norm(t1-t2) - norm(p-t1)*costheta(t2-t1, p-t1))/norm(t1-t2))

costheta(v1, v2) = dot(v1, v2)/(norm(v1)*norm(v2))

function optimal_mate(p1, target, pop)
    c = Vector{Float64}(undef, length(pop))
    idx = 0
    for i in 1:length(c)
        if pop[i] == p1
            idx = i
            continue
        end
        c[i] = costheta(target - sem(p1), sem(pop[i]) - sem(p1))
    end
    c[idx] = -Inf

    return pop[argmax(c)]
end

max_stars_repo_name: kogad/JuliaGSGP

max_issues_repo_name: kogad/JuliaGSGP

max_forks_repo_name: kogad/JuliaGSGP


text:
# L2 DistLoss
#Pkg.add("Latexify")

using LossFunctions
using LaTeXStrings, Measures, Latexify
gr(bg=:white)
l2fn = "L(r)=r^2"
l2dist = latexify(l2fn)
Plots.plot(L2DistLoss(), c=:olive, linewidth=2, annotation=(1, 3.3, L"L\left( r \right) = |r^{2}|", 18))

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
module Lehmann
using StaticArrays
using DelimitedFiles, LinearAlgebra
using Printf

include("utility/chebyshev.jl")

include("spectral.jl")
export Spectral

include("sample/sample.jl")
export Sample

include("discrete/builder.jl")
include("functional/builder.jl")

include("dlr.jl")
export DLRGrid

include("operation.jl")
export tau2dlr, dlr2tau, matfreq2dlr, dlr2matfreq, tau2matfreq, matfreq2tau, tau2tau, matfreq2matfreq

end


text:
module TestCapsule
using FactCheck
using SudokuSolver: solveSudoku

generateTestCase() = [0 0 0  0 0 0  0 0 0;
                      0 0 0  0 0 3  0 8 5;
                      0 0 1  0 2 0  0 0 0;

                      0 0 0  5 0 7  0 0 0;
                      0 0 4  0 0 0  1 0 0;
                      0 9 0  0 0 0  0 0 0;

                      5 0 0  0 0 0  0 7 3;
                      0 0 2  0 1 0  0 0 0;
                      0 0 0  0 4 0  0 0 9]

expectedSolution() = [9 8 7  6 5 4  3 2 1;
                      2 4 6  1 7 3  9 8 5;
                      3 5 1  9 2 8  7 4 6;

                      1 2 8  5 3 7  6 9 4;
                      6 3 4  8 9 2  1 5 7;
                      7 9 5  4 6 1  8 3 2;

                      5 1 9  2 8 6  4 7 3;
                      4 7 2  3 1 9  5 6 8;
                      8 6 3  7 4 5  2 1 9]

facts("System test.") do
    @time solutions = solveSudoku(generateTestCase())

    @fact solutions[1] => expectedSolution()
end

end #TestCapsule

max_stars_repo_name: Teo-ShaoWei/SudokuSolverSG

max_issues_repo_name: Teo-ShaoWei/SudokuSolverSG

max_forks_repo_name: Teo-ShaoWei/SudokuSolverSG


text:
"""Continuity equation's right-hand side with surface relaxation
-∂x(uh) - ∂y(vh) + γ*(η_ref - η)."""
function continuity_surf_relax!(η::Array{T,2},
                                Diag::DiagnosticVars{T,Tprog},
                                S::ModelSetup{T,Tprog},
                                t::Int) where {T,Tprog}

    @unpack dη = Diag.Tendencies
    @unpack dUdx,dVdy = Diag.VolumeFluxes
    @unpack η_ref = S.forcing
    @unpack γ = S.constants

    m,n = size(dη) .- (2,2)
    @boundscheck (m,n+2) == size(dUdx) || throw(BoundsError())
    @boundscheck (m+2,n) == size(dVdy) || throw(BoundsError())
    @boundscheck (m+2,n+2) == size(η) || throw(BoundsError())
    @boundscheck (m,n) == size(η_ref) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            dη[i+1,j+1] = -(Tprog(dUdx[i,j+1]) + Tprog(dVdy[i+1,j])) + Tprog(γ*(η_ref[i,j]-η[i+1,j+1]))
        end
    end
end

"""Continuity equation's right-hand side with time&space dependent forcing."""
function continuity_forcing!(   Diag::DiagnosticVars{T,Tprog},
                                S::ModelSetup{T,Tprog},
                                t::Int) where {T,Tprog}

    @unpack dη = Diag.Tendencies
    @unpack dUdx,dVdy = Diag.VolumeFluxes
    @unpack Fη = S.forcing

    m,n = size(dη) .- (2,2)         # cut off halo
    @boundscheck (m,n+2) == size(dUdx) || throw(BoundsError())
    @boundscheck (m+2,n) == size(dVdy) || throw(BoundsError())
    @boundscheck (m,n) == size(Fη) || throw(BoundsError())

    # avoid recomputation
    @unpack ωFη = S.constants
    Fηt = Ftime(T,t,ωFη)

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            dη[i+1,j+1] = -(Tprog(dUdx[i,j+1]) + Tprog(dVdy[i+1,j])) + Tprog(Fηt*Fη[i,j])
        end
    end
end

"""Continuity equation's right-hand side -∂x(uh) - ∂y(vh) without forcing."""
function continuity_itself!(Diag::DiagnosticVars{T,Tprog},
                            S::ModelSetup{T,Tprog},
                            t::Int) where {T,Tprog}

    @unpack dη = Diag.Tendencies
    @unpack dUdx,dVdy = Diag.VolumeFluxes

    m,n = size(dη) .- (2,2)     # cut off halo
    @boundscheck (m,n+2) == size(dUdx) || throw(BoundsError())
    @boundscheck (m+2,n) == size(dVdy) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            dη[i+1,j+1] = -(Tprog(dUdx[i,j+1]) + Tprog(dVdy[i+1,j]))
        end
    end
end

"""Transit function to call the specified continuity function."""
function continuity!(   η::AbstractMatrix,
                        Diag::DiagnosticVars,
                        S::ModelSetup,
                        t::Int)

    if S.parameters.surface_relax
        continuity_surf_relax!(η,Diag,S,t)
    elseif S.parameters.surface_forcing
        continuity_forcing!(Diag,S,t)
    else
        continuity_itself!(Diag,S,t)
    end
end


text:
struct NestedFilterSimulation{T1, T2, T3, T4, T5}
    hmodel::T1
    filter::T2
    hstate::T3
    fstate::T4
    tsteps::T5
end

function NestedFilterSimulation(
    N, p, q, σ, τ,
    Nrng, prng, qrng, σrng, τrng,
    m_out, m_in, width;
    timestep = RandomTimestep(Exponential(0.121))
)
    hmodel = ScalarBinomialModel(N, p, q, σ, τ)
    filter = NestedParticleFilter(width)
    hstate = ScalarBinomialState(N, 0)
    fstate = NestedParticleState(
                m_out, m_in,
                Nrng, prng, qrng, σrng, τrng
             )
    return NestedFilterSimulation(hmodel, filter, hstate, fstate, timestep)
end

function propagate!(sim::NestedFilterSimulation)
    obs = propagate_emit!(sim.hstate, sim.hmodel, sim.tsteps)
    update!(sim.fstate, obs, sim.filter)
    return obs
end

function run!(sim::NestedFilterSimulation; T::Int, plot_each_timestep = false)
    times = zeros(0)
    epsps = zeros(0)
    time = 0.
    for i in 1:T
        @time begin
            obs = propagate!(sim)
        end
        push!(times, time += obs.dt)
        push!(epsps, obs.EPSP)
        if plot_each_timestep
            posterior_plot(sim.fstate, times, epsps, truemodel = sim.hmodel)
        end
    end
    return times, epsps
end

MAP(sim::NestedFilterSimulation) = MAP(sim.fstate.model)


text:
using LinearAlgebra
using LoopVectorization
function RMP2{T}(refwfn::Fermi.HartreeFock.RHF,alg::Fermi.MollerPlesset.DF) where T <: AbstractFloat
    Fermi.MollerPlesset.print_header()
    ints = refwfn.ints
    nocc    = refwfn.ndocc
    nvir    = refwfn.nvir
    drop_occ = Fermi.CurrentOptions["drop_occ"]
    drop_vir = Fermi.CurrentOptions["drop_vir"]
    ints.orbs.frozencore = drop_occ
    ints.orbs.frozenvir = drop_vir
    eps = refwfn.eps
    ΔMP2 = 0.0
    @output "TEST\n"
    @output "\tComputing MP2 with DF algorithm\n\n"
    ttotal = @elapsed begin
    @output "\tComputing and transforming integrals ...\n"
    @output "\tComputing integrals ..."
    Fermi.Integrals.aux_ri!(ints)
    @output "\tBasis: {}\n" ints.bname["primary"]
    @output "\tDF basis: {}\n\n" ints.bname["aux"]
    t = @elapsed ints["B"]
    @output "\t done in {:5.2f} s\n" t

    ints.cache["B"] = convert(Array{T},ints.cache["B"])
    @output "\tTransforming integrals to MO basis ..."
    t = @elapsed Bov = ints["BOV"]
    @output "\t done in {:5.2f} s\n" t

    Bis = [zeros(T,size(Bov[:,1,:])) for i=1:Threads.nthreads()]
    Bjs = [zeros(T,size(Bis[1])) for i =1:Threads.nthreads()]
    BABs = [zeros(T,nvir,nvir) for i=1:Threads.nthreads()]
    BBAs = [zeros(T,nvir,nvir) for i=1:Threads.nthreads()]

    ΔMP2s = zeros(T,Threads.nthreads())
    epsd = eps[drop_occ+1:end] #occ - drop_occ
    epsv = eps[nocc+1:end] #virtuals
    _nvir = nvir-drop_vir
    _nocc = nocc-drop_occ
    t = @elapsed begin
        @sync for i in 1:_nocc
            Threads.@spawn begin
            id = Threads.threadid()
            Bi = Bis[id]
            Bj = Bjs[id]
            @views Bi[:,:] = Bov[:,i,:]
            BAB = BABs[id]
            BBA = BBAs[id]
            for j in i:_nocc
                if i != j
                    fac = 2
                else
                    fac = 1
                end
                @views Bj[:,:] .= Bov[:,j,:]
                @tensoropt begin
                    BAB[a,b] = Bi[Q,a]*Bj[Q,b]
                end
                transpose!(BBA,BAB)
                eij = epsd[i] + epsd[j]
                dmp2_1 = 0.0
                dmp2_2 = 0.0
                @fastmath for b in 1:_nvir
                    @inbounds eij_b = eij - epsv[b]
                    for a in 1:_nvir
                        @inbounds iajb = BAB[a,b]
                        #@inbounds ibja = BBA[a,b]
                        d = iajb/(eij_b - epsv[a])
                        @inbounds dmp2_1 += d*iajb
                        @inbounds dmp2_2 += d*BBA[a,b]
                    end
                end
                dmp2 = 2*dmp2_1 - dmp2_2
                ΔMP2s[id] += fac*dmp2
            end
            end
        end
        end
        ΔMP2 = sum(ΔMP2s)
    end
    @output "\t Energy computed in {:5.2f} s\n" t

    # generator function for T2 amplitudes, given ijab indices
    function gen(g::Fermi.GeneratedTensor,i,j,a,b)
        naux,nocc,nvir = size(g.data["Bov"])
        sum = 0.0
        for Q in 1:naux
            sum += g.data["Bov"][Q,i,a]*g.data["Bov"][Q,j,b]
        end
        eps = g.data["eps"]
        sum /= eps[i] + eps[j] - eps[a+nocc] - eps[b+nocc]
        sum
    end

    # data required to compute a T2 amplitude
    data = Dict{String,Any}(
                            "Bov" => Bov,
                            "eps" => eps
                           )

    no = convert(UInt16,nocc)
    nv = convert(UInt16,nvir)

    G = Fermi.GeneratedTensor{T}(gen,data,4,[no,no,nv,nv])

    @output "\n\tDF-MP2 energy is {:>14.10f}\n" ΔMP2
    @output "\tDF-MP2 done in {:>5.2f} s\n" ttotal
    @output repeat("-",80)*"\n"

    RMP2{T}(ΔMP2[],G)
end


text:
using Test
using ProxAL
using DelimitedFiles, Printf
using LinearAlgebra, JuMP, Ipopt
using CatViews
using CUDA
using MPI

MPI.Init()
DATA_DIR = joinpath(dirname(@__FILE__), "..", "data")
case = "case9"
T = 2
K = 0
ramp_scale = 0.5
load_scale = 1.0
rtol = 1e-2

# Load case
case_file = joinpath(DATA_DIR, "$(case).m")
load_file = joinpath(DATA_DIR, "mp_demand", "$(case)_oneweek_168")
# ctgs_arr = deepcopy(rawdata.ctgs_arr)

# Model/formulation settings
modelinfo = ModelInfo()
modelinfo.num_time_periods = T
modelinfo.load_scale = load_scale
modelinfo.ramp_scale = ramp_scale
modelinfo.allow_obj_gencost = true
modelinfo.allow_constr_infeas = false
modelinfo.time_link_constr_type = :penalty
modelinfo.case_name = case
modelinfo.num_ctgs = K

# Algorithm settings
algparams = AlgParams()
algparams.verbose = 0

# For JuMP
algparams.optimizer = optimizer_with_attributes(
    Ipopt.Optimizer,
    "print_level" => Int64(algparams.verbose > 0)*5,
)
# For ExaPF
algparams.gpu_optimizer = optimizer_with_attributes(
    Ipopt.Optimizer,
    "print_level" => 0,
    "limited_memory_max_history" => 50,
    "hessian_approximation" => "limited-memory",
    "tol" => 1e-5,
)

OPTIMAL_OBJVALUE = round(11258.316096599736*modelinfo.obj_scale, digits = 6)
OPTIMAL_PG = round.([0.8979870694509675, 1.3432060120295906, 0.9418738103137331, 0.9840203268625166, 1.448040098924617, 1.0149638876964715], digits = 5)
@testset "Test ProxAL on $(case) with $T-period, $K-ctgs, time_link=penalty and Ipopt" begin
    algparams.mode = :coldstart
    nlp = ProxALEvaluator(case_file, load_file, modelinfo, algparams, ProxAL.JuMPBackend())
    info = ProxAL.optimize!(nlp)
    @test isapprox(info.objvalue[end], OPTIMAL_OBJVALUE, rtol = rtol)
    @test isapprox(info.x.Pg[:], OPTIMAL_PG, rtol = rtol)
    @test isapprox(info.maxviol_c[end], 0.0)
    @test isapprox(info.maxviol_c_actual[end], 0.0)
    @test info.maxviol_t[end] <= algparams.tol
    @test info.maxviol_t_actual[end] <= algparams.tol
    @test info.maxviol_d[end] <= algparams.tol
    @test info.iter <= algparams.iterlim
end
MPI.Finalize()


text:
function initExternalEntities()
    externalEntities = findExternalModules("entities")
    append!(loadedEntities, externalEntities)
    loadModule.(externalEntities)
    loadExternalModules!(loadedModules, loadedEntities, "entities")
    registerPlacements!(entityPlacements, loadedEntities)
end

function initExternalTriggers()
    externalTriggers = findExternalModules("triggers")
    append!(loadedTriggers, externalTriggers)
    loadModule.(externalTriggers)
    loadExternalModules!(loadedModules, loadedTriggers, "triggers")
    registerPlacements!(triggerPlacements, loadedTriggers)
end

function initExternalEffects()
    externalEffects = findExternalModules("ffects")
    append!(loadedEffects, externalEffects)
    loadModule.(externalEffects)
    loadExternalModules!(loadedModules, loadedEffects, "effects")
    registerPlacements!(effectPlacements, loadedEffects)
end

function initExternalTools()
    externalTools = findExternalModules("tools")
    append!(loadedTools, externalTools)
    loadModule.(externalTools)
    loadExternalModules!(loadedModules, loadedTools, "tools")
end

function initExternalModules()
    initExternalEntities()
    initExternalTriggers()
    initExternalEffects()
    initExternalTools()
end


text:
@doc """

Calculates dfgls point optimal unit-root test statistic of Elliott, Rothenberg, and Stock (1996) value of the largest autorregressive coefficient under the alternative hypothesis given by ρ = 1+c/T, where c=-7 if trend="constant" and c=-13.5 if trend="linear")

"""->
function dfgls_test(y::Vector{Float64}; nar::Integer=5, trend::String="constant")
  #rewrite of mw's dfgls.prc
  t = length(y)
  dftau = NaN*zeros(1, 2)

  # @ -- GLS estimate of constant term -- @
  c = -7
  rho = (1+c/t)
  yrho = [y[1]; y[2:t]-rho*y[1:t-1]] #@ Quasi Differenced Y series @
  xcrho = [1.0; (1-rho)*ones(t-1)] #@ Quasi Difference Constant @
  mu = ((xcrho'*xcrho).^(-1)).*(xcrho'yrho)

  c = -13.5
  rho = (1+c/t)
  yrho = [y[1]; y[2:t]-rho*y[1:t-1]] #@ Quasi Differenced Y series @
  xcrho = [1.0; (1-rho)*ones(t-1)] #@ Quasi Difference Constant @
  trnd = collect(1:t)
  xtrho = [trnd[1]; trnd[2:t]-rho*trnd[1:t-1]]
  x = [xcrho xtrho]
  mut = inv(x'x)*(x'yrho)
  #  @ -- detrend y -- @
  y1 = y - mu.*ones(t)                #@ GLS Demeaned @
  y2 = y - [ones(t,1) trnd]*mut        #@ GLS Detrended @
  dy1 = diff(y1)
  dy2 = diff(y2)
  #compute adf
  x1 = ar_lagmatrix(dy1, nar)
  x2 = ar_lagmatrix(dy2, nar)

  z = [y1[1:size(x1,1)] x1]
  y = z[:, 1]
  x = z[:, 2:end]
  xxi = inv(x'x)
  xy = x'y
  bet = xxi*xy
  ssr = (y'y)-(xy)'bet
  s2 = ssr/(size(y,1)-size(x,2))
  b = bet[1]
  vb = s2.*xxi[1,1]
  dfmu = (b-1)./sqrt(vb)

  z = [y1[1:size(x1,1)] x1]
  y = z[:, 1]
  x = z[:, 2:end]
  xxi = inv(x'x)
  xy = x'y
  bet = xxi*xy
  ssr = (y'y)-(xy)'bet
  s2 = ssr/(size(y,1)-size(x,2))
  b = bet[1]
  vb = s2.*xxi[1,1]
  dftau = (b-1)./sqrt(vb)
  (trend=="constant") && return dfmu[1]
  (trend=="linear") && return dftau[1]
  (trend=="both") && return dfmu[1], dftau[1]
  end


text:
# d for integers:
@printf "%d\n" 1e5 #> 100000
x = 7.35679
# f = float format, rounded if needed:
@printf "x = %0.3f\n" x #> 7.357
aa = 1.5231071779744345
bb = 33.976886930000695
@printf("%.2f %.2f\n", aa, bb) #=> 1.52 33.98
# or to create another string:
str = @sprintf("%0.3f", x)
show(str) #> "7.357"
println()
# e = scientific format with e:
@printf "%0.6e\n" x #> 7.356790e+00
# c = for characters:
@printf "output: %c\n" 'α' #> output: α
# s for strings:
@printf "%s\n" "I like Julia"
# right justify:
@printf "%50s\n" "width 25, text right justified!"
# p for pointers:
@printf "a pointer: %p\n" 1e10 #> a pointer: 0x00000002540be400

max_stars_repo_path: Module 1/Chapter02/Old/formatting.jl

max_issues_repo_path: Module 1/Chapter02/Old/formatting.jl

max_forks_repo_path: Module 1/Chapter02/Old/formatting.jl


text:
# Author: Pedro Henrique Pereira Braga
# Date: 2019-01-25

# Add packages using Pkg.add("PackageName")
Pkg.add("Diversity")
Pkg.add("Phylo")
Pkg.add("SpatialEcology")
Pkg.add("Plots")
Pkg.add("JLD2")
Pkg.add("SparseArrays")
Pkg.add("DataFrames")
Pkg.add("Clustering")
Pkg.add("Distances")

# Load packages running
# using PackageName
using Diversity, Phylo
using SpatialEcology
using Plots
# You can load multiple pkgs separating them by commas
using JLD2, SparseArrays, DataFrames
using Distances

# Set working directory by: cd("../mydir/someotherdir/")

cd("C:\\Users\\phper\\Documents\\GitHub\\pedrohbraga\\ExploringJulia")

# Creating species
species = ["Dog", "Cat", "Human", "Monkey", "Alien"];

# Creating communities
communities = DataFrame(Dog = [1, 1, 1],
                        Cat = [0, 1, 1],
                        Human = [1, 0, 0],
                        Monkey = [1, 0, 1],
                        Alien = [0, 0, 1])
communities

# Get column names for the communitiesDataFrame
names(communities)

# Create a random phylogenetic tree on species
nt = rand(Nonultrametric(species))
nt

# Allows you to see the branch length information
collect(getbranches(nt))
getbranches(nt)

# Create a metaphylo object
# Before, since metaphylo works with Arrays, I had to convert communities to SparseArrays
communitiesArray = convert(Array, communities)
# metaphylo puts species in rows and sites in columns. Therefore, I needed to
# use the ' conjugate transpose operator
communitiesArray = communitiesArray'

# Using the stack and unstack functions from DataFrame, I was able to do the Same
# In summary, this creates an id column:
communities[:id] = 1:size(communities, 1)

# Stacks the information in a sort of long format
communitiesStack = stack(communities, names(communities)[1:length(names(communities))-1])
communitiesStack

# Unstack the stacked data frame, to then obtain a transposed one
communitiesUnstack = unstack(communitiesStack, :variable, :id, :value)

# Converts it to an Array, also taking away the id row
communitiesArray = convert(Array, communitiesUnstack[:,2:size(communitiesUnstack)[2]])

# Finally, create the Metaphylo object
metaphylo = Metacommunity(communitiesArray, PhyloTypes(nt));

# One can access what is inside the metaphylo object like This
?metaphylo
metaphylo.types # Returns the PhyloTypes within it
metaphylo.rawabundances # The raw community data
metaphylo.rawabundances

# Returns the names of the tips of the phylogeny
leafnames = gettypenames(metaphylo, true)
# Same as "getleafnames(nt)"

# Calculates gamma diversity
meta_gamma(metaphylo, 0)

# Can I calculate pairwise distances on the phylogeny?
?pairwise
distances(nt) # This seems to work!

max_stars_repo_name: pedrohbraga/ExploringJulia

max_issues_repo_name: pedrohbraga/ExploringJulia

max_forks_repo_name: pedrohbraga/ExploringJulia


text:
function jacobian(m::RepDSGEGovDebt)
    endo = augment_model_states(m.endogenous_states, n_model_states(m))
    eq   = m.equilibrium_conditions

    # Load in parameters/steady-state parameters
    c_scalar   = m[:cstar].value
    ell_scalar = m[:lstar].value
    ω          = m[:ωstar].value
    H          = m[:H].value
    ϕh         = m[:ϕh].value
    lam        = 1/c_scalar
    ϕ          = lam*ω/(H^ϕh)
    βrank      = exp(m[:γ].value)/(1+m[:r])
    y          = m[:ystar].value
    g          = m[:g].value
    Rk         = m[:Rkstar].value
    T          = m[:Tstar].value
    k          = m[:kstar].value
    x          = m[:xstar].value
    Tg         = m[:Tg].value
    spp        = m[:spp].value
    R          = m[:r].scaledvalue + 1
    γ          = m[:γ].scaledvalue
    δ          = m[:δ].value
    κw         = m[:κ_w].value
    κp         = m[:κ_p].value
    α          = m[:α].value
    ρR         = m[:ρR].value
    ψπ         = m[:ψπ].value
    ψy         = m[:ψy].value
    δb         = m[:δb].value
    bg         = m[:bg].value
    ρB         = m[:ρ_b].value
    ρG         = m[:ρ_g].value
    ρZ         = m[:ρ_z].value
    ρμ         = m[:ρ_μ].value
    ρlamw      = m[:ρ_λ_w].value
    ρlamf      = m[:ρ_λ_f].value
    ρmon       = m[:ρ_rm].value

    # make the Jacobian
    JJ = zeros(n_model_states(m), 2*n_model_states(m))

    # Euler equation
    JJ[eq[:eq_euler],endo[:l′_t]] = 1.
    JJ[eq[:eq_euler],endo[:z′_t]] = -1.
    JJ[eq[:eq_euler],endo[:b_t]]  = ell_scalar
    JJ[eq[:eq_euler],endo[:l_t]]  = -1.
    JJ[eq[:eq_euler],endo[:R_t]]  = 1.

    # mkt clearing
    JJ[eq[:eq_market_clearing],endo[:y_t]]   = y/g
    JJ[eq[:eq_market_clearing],endo[:g_t]]   = -y/g
    JJ[eq[:eq_market_clearing],endo[:I_t]]   = -x
    JJ[eq[:eq_market_clearing],endo[:l_t]] = c_scalar

    # lambda = average marginal utility
    JJ[eq[:eq_lambda],endo[:margutil_t]] = lam
    JJ[eq[:eq_lambda],endo[:l_t]] = -1/c_scalar

    # transfer
    JJ[eq[:eq_transfers],endo[:t_t]]  = T
    JJ[eq[:eq_transfers],endo[:capreturn_t]]  = -Rk*k
    JJ[eq[:eq_transfers],endo[:k_t]]  = -Rk*k
    JJ[eq[:eq_transfers],endo[:z_t]]   = Rk*k
    JJ[eq[:eq_transfers],endo[:I_t]]   = x
    JJ[eq[:eq_transfers],endo[:mc_t]]  = y
    JJ[eq[:eq_transfers],endo[:tg_t]]  = Tg

    # investment
    JJ[eq[:eq_investment],endo[:Q_t]]  = 1.
    JJ[eq[:eq_investment],endo[:μ_t]]  = 1.
    JJ[eq[:eq_investment],endo[:I′_t]] = spp*(exp(3*γ))/R
    JJ[eq[:eq_investment],endo[:z′_t]] = spp*(exp(3*γ))/R
    JJ[eq[:eq_investment],endo[:I_t]]  = -spp*(exp(3*γ))/R - spp*exp(2*γ)
    JJ[eq[:eq_investment],endo[:I_t1]] = spp*exp(2*γ)
    JJ[eq[:eq_investment],endo[:z_t]]  = -spp*exp(2*γ)

    # tobin's q
    JJ[eq[:eq_tobin_q],endo[:margutil_t]]  = 1.
    JJ[eq[:eq_tobin_q],endo[:margutil′_t]] = -1.
    JJ[eq[:eq_tobin_q],endo[:Q_t]]    = 1.
    JJ[eq[:eq_tobin_q],endo[:z′_t]]   = 1.
    JJ[eq[:eq_tobin_q],endo[:capreturn′_t]]  = -Rk/R
    JJ[eq[:eq_tobin_q],endo[:Q′_t]]   = -(1-δ)/R

    # capital accumulation
    JJ[eq[:eq_capital_accumulation],endo[:k′_t]] = 1.
    JJ[eq[:eq_capital_accumulation],endo[:k_t]]  = -(1-δ)
    JJ[eq[:eq_capital_accumulation],endo[:z_t]]   = (1-δ)
    JJ[eq[:eq_capital_accumulation],endo[:μ_t]]  = -x/k
    JJ[eq[:eq_capital_accumulation],endo[:I_t]]   = -x/k

    # wage phillips curve
    JJ[eq[:eq_wage_phillips],endo[:π_w_t]]  = -1.
    JJ[eq[:eq_wage_phillips],endo[:λ_w_t]] = 1.
    JJ[eq[:eq_wage_phillips],endo[:L_t]]   = κw*ϕh
    JJ[eq[:eq_wage_phillips],endo[:w_t]]    = -κw
    JJ[eq[:eq_wage_phillips],endo[:π_w′_t]]  = βrank
    JJ[eq[:eq_wage_phillips],endo[:margutil_t]]  = -κw

    # price phillips curve
    JJ[eq[:eq_price_phillips],endo[:π_t]]   = -1.
    JJ[eq[:eq_price_phillips],endo[:mc_t]]   = κp
    JJ[eq[:eq_price_phillips],endo[:λ_f_t]] = 1.
    JJ[eq[:eq_price_phillips],endo[:π′_t]]  = 1/R

    # marginal cost
    JJ[eq[:eq_marginal_cost],endo[:mc_t]] = 1.
    JJ[eq[:eq_marginal_cost],endo[:w_t]]  = -(1-α)
    JJ[eq[:eq_marginal_cost],endo[:capreturn_t]] = -α

    # gdp
    JJ[eq[:eq_gdp],endo[:y_t]]  = 1.
    JJ[eq[:eq_gdp],endo[:z_t]]  = α
    JJ[eq[:eq_gdp],endo[:k_t]] = -α
    JJ[eq[:eq_gdp],endo[:L_t]] = -(1-α)

    # optimal k/l ratio
    JJ[eq[:eq_optimal_kl],endo[:capreturn_t]] = 1.
    JJ[eq[:eq_optimal_kl],endo[:w_t]]  = -1.
    JJ[eq[:eq_optimal_kl],endo[:L_t]] = -1.
    JJ[eq[:eq_optimal_kl],endo[:k_t]] = 1.
    JJ[eq[:eq_optimal_kl],endo[:z_t]]  = -1.

    # taylor rule
    JJ[eq[:eq_taylor],endo[:i_t]]   = -1.
    JJ[eq[:eq_taylor],endo[:i_t1]]  = ρR
    JJ[eq[:eq_taylor],endo[:π_t]]   = (1-ρR)*ψπ
    JJ[eq[:eq_taylor],endo[:y_t]]    = (1-ρR)*ψy
    JJ[eq[:eq_taylor],endo[:y_t1]]  = -(1-ρR)*ψy
    JJ[eq[:eq_taylor],endo[:z_t]]    = (1-ρR)*ψy
    JJ[eq[:eq_taylor],endo[:rm_t]] = 1.

    # fisher eqn
    JJ[eq[:eq_fisher],endo[:R_t]]  = 1.
    JJ[eq[:eq_fisher],endo[:π′_t]] = 1.
    JJ[eq[:eq_fisher],endo[:i_t]]  = -1.

    # wage inflation
    JJ[eq[:eq_nominal_wage_inflation],endo[:π_w_t]] = 1.
    JJ[eq[:eq_nominal_wage_inflation],endo[:π_t]]  = -1.
    JJ[eq[:eq_nominal_wage_inflation],endo[:z_t]]   = -1.
    JJ[eq[:eq_nominal_wage_inflation],endo[:w_t]]   = -1.
    JJ[eq[:eq_nominal_wage_inflation],endo[:w_t1]]  = 1.

    # fiscal rule
    JJ[eq[:eq_fiscal_rule],endo[:tg_t]]  = -Tg
    JJ[eq[:eq_fiscal_rule],endo[:R_t]]   = δb*bg/R
    JJ[eq[:eq_fiscal_rule],endo[:bg_t]]  = δb*bg*ℯ^(-γ)
    JJ[eq[:eq_fiscal_rule],endo[:z_t]]   = -δb*bg*ℯ^(-γ)
    JJ[eq[:eq_fiscal_rule],endo[:y_t]]   = δb*(1-(1/g))*y
    JJ[eq[:eq_fiscal_rule],endo[:g_t]]   = δb*y/g

    # govt budget constraint
    JJ[eq[:eq_g_budget_constraint],endo[:bg′_t]] = -(bg/R)
    JJ[eq[:eq_g_budget_constraint],endo[:R_t]]   = bg/R
    JJ[eq[:eq_g_budget_constraint],endo[:bg_t]]  = bg*ℯ^(-γ)
    JJ[eq[:eq_g_budget_constraint],endo[:z_t]]   = -bg*ℯ^(-γ)
    JJ[eq[:eq_g_budget_constraint],endo[:y_t]]   = (1-(1/g))*y
    JJ[eq[:eq_g_budget_constraint],endo[:g_t]]   = y/g
    JJ[eq[:eq_g_budget_constraint],endo[:tg_t]]  = -Tg

    # update lagged variables
    JJ[eq[:LR],endo[:R′_t1]] = 1.
    JJ[eq[:LR],endo[:R_t]]   = -1.

    JJ[eq[:LI],endo[:i′_t1]] = 1.
    JJ[eq[:LI],endo[:i_t]]   = -1.

    JJ[eq[:LY],endo[:y′_t1]] = 1.
    JJ[eq[:LY],endo[:y_t]]   = -1.

    JJ[eq[:LW],endo[:w′_t1]] = 1.
    JJ[eq[:LW],endo[:w_t]]   = -1.

    JJ[eq[:LX],endo[:I′_t1]] = 1.
    JJ[eq[:LX],endo[:I_t]]   = -1.

    # discount factor shock
    JJ[eq[:eq_b],endo[:b′_t]] = 1.
    JJ[eq[:eq_b],endo[:b_t]]  = -ρB

    # g/y shock
    JJ[eq[:eq_g],endo[:g′_t]] = 1.
    JJ[eq[:eq_g],endo[:g_t]]  = -ρG

    # tfp growth shock
    JJ[eq[:eq_z],endo[:z′_t]] = 1.
    JJ[eq[:eq_z],endo[:z_t]]  = -ρZ

    # investment shock
    JJ[eq[:eq_μ],endo[:μ′_t]] = 1.
    JJ[eq[:eq_μ],endo[:μ_t]]  = -ρμ

    # wage mkup shock
    JJ[eq[:eq_λ_w],endo[:λ_w′_t]] = 1.
    JJ[eq[:eq_λ_w],endo[:λ_w_t]]  = -ρlamw

    # price mkup shock
    JJ[eq[:eq_λ_f],endo[:λ_f′_t]] = 1.
    JJ[eq[:eq_λ_f],endo[:λ_f_t]] = -ρlamf

    # monetary policy shock
    JJ[eq[:eq_rm],endo[:rm′_t]] = 1.
    JJ[eq[:eq_rm],endo[:rm_t]]  = -ρmon

    if !m.testing
        JJ  = normalize(m, JJ)
    end

    return JJ
end

function compose_normalization_matrices(m::RepDSGEGovDebt)
    Qleft     = eye(n_model_states(m))
    Qx        = eye(n_backward_looking_states(m))
    Qy        = eye(n_jumps(m))

    Qright    = cat(Qx', Qy', Qx', Qy', dims = [1,2])

    return Qx, Qy, Qleft, Qright
end

function normalize(m::RepDSGEGovDebt, JJ::Matrix{Float64})
    Qx, _, Qleft, Qright = compose_normalization_matrices(m)

    m <= Setting(:n_predetermined_variables, size(Qx, 1))

    Jac1 = Qleft*JJ*Qright

    return Jac1
end


text:
import Dagger: @par, @spawn, spawn
import Dagger: Chunk

@everywhere begin
    checkwid() = myid()==1

    function dynamic_fib(n)
        n <= 1 && return n
        t = Dagger.spawn(dynamic_fib, n-1)
        y = dynamic_fib(n-2)
        return (fetch(t)::Int) + y
    end

    struct ProcessLockedStruct
        x::Ptr{Int} # Zero'd during serialization
    end
    (pls::ProcessLockedStruct)(x) = x+UInt(pls.x)

    struct MulProc <: Dagger.Processor
        owner::Int
    end
    MulProc() = MulProc(myid())
    Dagger.get_parent(mp::MulProc) = OSProc(mp.owner)
    Dagger.move(src::MulProc, dest::Dagger.ThreadProc, x::Function) = Base.:*
end

@testset "@par" begin
    @testset "per-call" begin
        x = 2
        a = @par x + x
        @test a isa Dagger.Thunk
        b = @par sum([x,1,2])
        c = @par a * b
        @test collect(c) == 20
    end
    @testset "block" begin
        c = @par begin
            x = 2
            a = x + x
            b = sum([x,1,2])
            c = a * b
        end
        @test x isa Int
        @test a isa Dagger.Thunk
        @test c isa Dagger.Thunk
        @test collect(c) == 20
    end
end

@testset "@spawn" begin
    @test !isassigned(Dagger.Sch.EAGER_CONTEXT)
    @testset "per-call" begin
        x = 2
        a = @spawn x + x
        @test a isa Dagger.EagerThunk
        b = @spawn sum([x,1,2])
        c = spawn(*, a, b)
        @test c isa Dagger.EagerThunk
        @test fetch(a) == 4
        @test fetch(b) == 5
        @test fetch(c) == 20
    end
    @test Dagger.Sch.EAGER_CONTEXT[] isa Context
    @testset "waiting" begin
        a = @spawn sleep(1)
        @test !isready(a)
        wait(a)
        @test isready(a)
    end
    @testset "options" begin
        s = 1
        m = true
        a = @spawn single=s checkwid()
        b = @spawn meta=m ((_a)->_a isa Dagger.Chunk)(a)
        @test fetch(a)
        @test fetch(b)
    end
    @testset "errors" begin
        @testset "independent" begin
            a = @spawn error("Test")
            wait(a)
            @test isready(a)
            @test_throws_unwrap Dagger.ThunkFailedException fetch(a)
            b = @spawn 1+2
            @test fetch(b) == 3
        end
        @testset "direct vs indirect" begin
            a = @spawn error("Test")
            b = @spawn a+1

            ex = try
                fetch(a)
            catch err
                err
            end
            ex = Dagger.Sch.unwrap_nested_exception(ex)
            ex_str = sprint(io->Base.showerror(io,ex))
            @test occursin(r"^ThunkFailedException \(Thunk.*failure\):", ex_str)
            @test occursin("Test", ex_str)
            @test !occursin("due to a failure in", ex_str)

            ex = try
                fetch(b)
            catch err
                err
            end
            ex = Dagger.Sch.unwrap_nested_exception(ex)
            ex_str = sprint(io->Base.showerror(io,ex))
            @test occursin(r"Thunk.*failure due to a failure in", ex_str)
            @test occursin("Test", ex_str)
        end
        @testset "single dependent" begin
            a = @spawn error("Test")
            b = @spawn a+2
            @test_throws_unwrap Dagger.ThunkFailedException fetch(a)
        end
        @testset "multi dependent" begin
            a = @spawn error("Test")
            b = @spawn a+2
            c = @spawn a*2
            @test_throws_unwrap Dagger.ThunkFailedException fetch(b)
            @test_throws_unwrap Dagger.ThunkFailedException fetch(c)
        end
        @testset "dependent chain" begin
            a = @spawn error("Test")
            @test_throws_unwrap Dagger.ThunkFailedException fetch(a)
            b = @spawn a+1
            @test_throws_unwrap Dagger.ThunkFailedException fetch(b)
            c = @spawn b+2
            @test_throws_unwrap Dagger.ThunkFailedException fetch(c)
        end
        @testset "single input" begin
            a = @spawn 1+1
            b = @spawn (a->error("Test"))(a)
            @test fetch(a) == 2
            @test_throws_unwrap Dagger.ThunkFailedException fetch(b)
        end
        @testset "multi input" begin
            a = @spawn 1+1
            b = @spawn 2*2
            c = @spawn ((a,b)->error("Test"))(a,b)
            @test fetch(a) == 2
            @test fetch(b) == 4
            @test_throws_unwrap Dagger.ThunkFailedException fetch(c)
        end
        @testset "diamond" begin
            a = @spawn 1+1
            b = @spawn a+1
            c = @spawn a*2
            d = @spawn ((b,c)->error("Test"))(b,c)
            @test fetch(a) == 2
            @test fetch(b) == 3
            @test fetch(c) == 4
            @test_throws_unwrap Dagger.ThunkFailedException fetch(d)
        end
    end
    @testset "remote spawn" begin
        a = fetch(Distributed.@spawnat 2 Dagger.spawn(+, 1, 2))
        @test Dagger.Sch.EAGER_INIT[]
        @test fetch(Distributed.@spawnat 2 !(Dagger.Sch.EAGER_INIT[]))
        @test a isa Dagger.EagerThunk
        @test fetch(a) == 3

        # Mild stress-test
        @test dynamic_fib(10) == 55
    end
    @testset "undefined function" begin
        # Issues #254, #255

        # only defined on head node
        @eval evil_f(x) = x

        eager_thunks = map(1:10) do i
            single = isodd(i) ? 1 : first(workers())
            Dagger.@spawn single=single evil_f(i)
        end

        errored(t) = try
            fetch(t)
            false
        catch
            true
        end
        @test any(t->errored(t), eager_thunks)
        @test any(t->!errored(t), eager_thunks)
    end
    @testset "function chunks" begin
        @testset "lazy API" begin
            a = delayed(+)(1,2)
            @test !(a.f isa Chunk)

            a = delayed(+; scope=NodeScope())(1,2)
            @test a.f isa Chunk
            @test a.f.processor isa OSProc
            @test a.f.scope isa NodeScope

            a = delayed(+; processor=Dagger.ThreadProc(1,1))(1,2)
            @test a.f isa Chunk
            @test a.f.processor isa Dagger.ThreadProc
            @test a.f.scope isa AnyScope

            a = delayed(+; processor=Dagger.ThreadProc(1,1), scope=NodeScope())(1,2)
            @test a.f isa Chunk
            @test a.f.processor isa Dagger.ThreadProc
            @test a.f.scope isa NodeScope

            @testset "Scope Restrictions" begin
                pls = ProcessLockedStruct(Ptr{Int}(42))
                ctx = Context([1, workers()...])

                # Negative test
                @test_skip !all(x->x==43, collect(ctx, delayed(vcat)([delayed(pls)(1) for i in 1:10]...)))
                # Positive tests (no serialization)
                @test all(x->x==43, collect(ctx, delayed(vcat)([delayed(pls; scope=ProcessScope())(1) for i in 1:10]...)))
                @test all(x->x==1, collect(ctx, delayed(vcat)([delayed(pls; scope=ProcessScope(first(workers())))(1) for i in 1:10]...)))
            end
            @testset "Processor Data Movement" begin
                @everywhere Dagger.add_processor_callback!(()->MulProc(), :mulproc)
                @test collect(delayed(+; processor=MulProc())(3,4)) == 12
                @everywhere Dagger.delete_processor_callback!(:mulproc)
            end
        end
        @testset "eager API" begin
            _a = Dagger.spawn(+, 1, 2; scope=NodeScope())
            a = Dagger.Sch._find_thunk(_a)
            @test a.f isa Chunk
            @test a.f.processor isa OSProc
            @test a.f.scope isa NodeScope

            _a = Dagger.spawn(+, 1, 2; processor=Dagger.ThreadProc(1,1))
            a = Dagger.Sch._find_thunk(_a)
            @test a.f isa Chunk
            @test a.f.processor isa Dagger.ThreadProc
            @test a.f.scope isa AnyScope

            _a = Dagger.spawn(+, 1, 2; processor=Dagger.ThreadProc(1,1), scope=NodeScope())
            a = Dagger.Sch._find_thunk(_a)
            @test a.f isa Chunk
            @test a.f.processor isa Dagger.ThreadProc
            @test a.f.scope isa NodeScope
        end
    end
    @testset "parent fetch child, one thread" begin
        # Issue #282

        s = p -> p == Dagger.ThreadProc(1, 1)
        f = (x) -> 10 + x
        g = (x) -> fetch(Dagger.spawn(f, x; proclist=s))
        fetch(Dagger.spawn(g, 10; proclist=s))
    end
end


text:
# makes /opt/julia a shared depot path so that users can directly find and load
# pre-installed packages
#
# For non-root users, /opt/julia is read-only, trying to add packages would be impossible.
# A workaround for this is to start julia with a project folder that is writable, e.g.,
# `julia --project=$HOME/.julia/env/1.5"`

SHARE_DIR = ENV["JULIA_DEPOT_PATH"]

empty!(DEPOT_PATH)
push!(DEPOT_PATH, joinpath(homedir(), ".julia"), SHARE_DIR)
push!(LOAD_PATH, SHARE_DIR)

max_stars_repo_name: ksqsf/ustcmirror-images

max_issues_repo_name: ksqsf/ustcmirror-images

max_forks_repo_name: ksqsf/ustcmirror-images


text:
# ---
# title: 1551. Minimum Operations to Make Array Equal
# id: problem1551
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/minimum-operations-to-make-array-equal/description/>
# hidden: true
# ---
#
# You have an array `arr` of length `n` where `arr[i] = (2 * i) + 1` for all
# valid values of `i` (i.e. `0 <= i < n`).
#
# In one operation, you can select two indices `x` and `y` where `0 <= x, y < n`
# and subtract `1` from `arr[x]` and add `1` to `arr[y]` (i.e. perform `arr[x]
# -=1 `and `arr[y] += 1`). The goal is to make all the elements of the array
# **equal**. It is **guaranteed** that all the elements of the array can be made
# equal using some operations.
#
# Given an integer `n`, the length of the array. Return _the minimum number of
# operations_ needed to make all the elements of arr equal.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 3
#     Output: 2
#     Explanation: arr = [1, 3, 5]
#     First operation choose x = 2 and y = 0, this leads arr to be [2, 3, 4]
#     In the second operation choose x = 2 and y = 0 again, thus arr = [3, 3, 3].
#
#
# **Example 2:**
#
#
#
#     Input: n = 6
#     Output: 9
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 10^4`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
abstract LossFunction

type LogisticLoss <: LossFunction
    obj::Function
    gradient::Function
    hessian::Function

    function obj(y::Array, y_pred::Array)
        logistic_pred = sigmoid(y_pred)
        l = - y .* log(logistic_pred)
        r = - (1-y) .* log(1-logistic_pred)
        loss = 1/size(y,1) * (l + r)
        return loss
    end
    function gradient(y::Array, y_pred::Array)
        return (sigmoid(y_pred) - y)
    end

    function hessian(y::Array, y_pred::Array)
        logitsc_pred = sigmoid(y_pred)
        val = logitsc_pred .* ( 1 - logitsc_pred)
        return val
    end

    function LogisticLoss()
        loss = new(obj, gradient, hessian)
        return loss
    end

end

type SquareLoss <: LossFunction
    obj::Function
    gradient::Function
    hessian::Function

    function obj(y::Array, y_pred::Array)
        return 0.5 * sumabs2(y - y_pred)
    end
    function gradient(y::Array, y_pred::Array)
        return (y_pred - y)
    end

    function hessian(y::Array, y_pred::Array)
        n_sample = size(y,1)
        return ones(n_sample)
    end

    function SquareLoss()
        loss = new(obj, gradient, hessian)
        return loss
    end
end

type DecisionNode
    label::Union{Array, Int64, String, Float64}
    feature_index::Integer
    threshold::Features
    true_branch::Union{DecisionNode,String}
    false_branch::Union{DecisionNode, String}
    y_num::Int64
end

function DecisionNode(;
                      label = "nothing",
                      feature_index = 0,
                      threshold = Inf,
                      true_branch = "nothing",
                      false_branch = "nothing",
                      y_num = 1)
    return DecisionNode(label, feature_index, threshold,true_branch, false_branch, y_num)
end

abstract DecisionTree

type RegressionTree <: DecisionTree
    root::Union{DecisionNode,String}
    max_depth::Integer
    min_gain::Float64
    min_samples_split::Integer
    current_depth::Integer
    y_num::Integer
end

type ClassificationTree <: DecisionTree
    root::Union{DecisionNode,String}
    max_depth::Integer
    min_gain::Float64
    min_samples_split::Integer
    current_depth::Integer
    y_num::Integer
end

type XGBoostRegressionTree <: DecisionTree
    root::Union{DecisionNode,String}
    max_depth::Int64
    min_gain::Float64
    min_samples_split::Int64
    current_depth::Int64
    loss::LogisticLoss
    y_num::Integer
end

function XGBoostRegressionTree(;
                      root = "nothing",
                      min_samples_split = 2,
                      min_gain = 1e-7,
                      max_depth = 1e7,
                      current_depth::Int64 = 0,
                      loss::LogisticLoss = LogisticLoss(),
                      y_num = 1)
    return XGBoostRegressionTree(root, max_depth, min_gain,
        min_samples_split, current_depth, loss, y_num)
end

function ClassificationTree(;
                      root = "nothing",
                      min_samples_split = 2,
                      min_gain = 1e-7,
                      max_depth = 1e7,
                      current_depth = 0,
                      y_num = 1)
    return ClassificationTree(root, max_depth, min_gain,
        min_samples_split, current_depth, y_num)
end

function RegressionTree(;
                      root = "nothing",
                      min_samples_split = 2,
                      min_gain = 1e-7,
                      max_depth = 1e7,
                      current_depth = 0,
                      y_num = 1)
    return RegressionTree(root, max_depth, min_gain,
        min_samples_split, current_depth, y_num)
end

function train!(model::DecisionTree, X::Matrix, y::Array)
    model.current_depth = 0
    #Normalize
    #X_y = [X y]
    #X_y = normalize_(X_y)
    #X = X_y[:, 1:(end-1)]
    #y = X_y[:, end]
    model.root = build_tree(model, X, y)
end

function build_tree(model::DecisionTree, X::Matrix, y::Array)
    if typeof(model) <: XGBoostRegressionTree
        model.y_num = trunc(Int,size(y,2)/2)
    else
        model.y_num = size(y,2)
    end
    entropy = calc_entropy(y)
    largest_impurity = 0
    best_criteria = 0
    best_sets = 0
    n_features = size(X, 2)
    X_y = hcat(X, y)
    n_samples = size(X_y,1)
    if n_samples >= model.min_samples_split
        for i = 1:n_features
            feature_values = X_y[:, i]
            unique_values = unique(feature_values)
            #For large regression problem
            if typeof(model) == XGBoostRegressionTree
                num = 8
                if length(unique_values) >= num
                    num_ = length(unique_values)
                    indd = randperm(num_)[1:num]
                    unique_values = unique_values[indd]
                end
            end
            for threshold in unique_values
                Xy_1, Xy_2 = split_at_feature(X_y, i, threshold)
                if size(Xy_1,1) > 0 && size(Xy_2,1) > 0
                    y1 = Xy_1[:, (n_features+1):end]
                    y2 = Xy_2[:, (n_features+1):end]
                    impurity = impurity_calc(model, y, y1, y2)
                    if impurity > largest_impurity
                        #println("$(impurity)")
                        largest_impurity = impurity
                        best_criteria = Dict("feature_i" => i, "threshold" => threshold)
                        best_sets = Dict("left_branch" => Xy_1, "right_branch" => Xy_2)
                    end
                end
            end
        end
    end
    if model.current_depth < model.max_depth && largest_impurity > model.min_gain
        leftX, leftY = best_sets["left_branch"][:, 1:n_features], best_sets["left_branch"][:, (n_features+1):end]
        rightX, rightY = best_sets["right_branch"][:, 1:n_features], best_sets["right_branch"][:, (n_features+1):end]
        true_branch = build_tree(model, leftX, leftY)
        false_branch = build_tree(model, rightX, rightY)
        model.current_depth += 1
        return DecisionNode(feature_index = best_criteria["feature_i"],
         threshold = best_criteria["threshold"], true_branch = true_branch,
         false_branch = false_branch)
    end

    ## if not constructed

    leaf_value = leaf_value_calc(model, y)
    return DecisionNode(label = leaf_value)
end

function leaf_value_calc(model::XGBoostRegressionTree, y::Array)
    y, y_pred = split_(y)
    nominator = sum(model.loss.gradient(y, y_pred),1)
    denominator = sum(model.loss.hessian(y, y_pred),1)

    return nominator ./ denominator
end

function leaf_value_calc(model::RegressionTree, y::Array)
    return mean(y)
end

function leaf_value_calc(model::ClassificationTree, y::Array)
    labels = size(y,2)
    if labels > 1
        y = unhot(y)
    end
    feature = unique(y)
    most_common = nothing
    count_max = 0
    for i in feature
        count = sum(y .== i)
        if count > count_max
            count_max = count
            most_common = i
        end
    end
    if labels > 1
        most_common = trunc(Int64, most_common[1])
        most_common = eye(labels)[most_common,:]
    end
    return most_common
end

function split_(y)
    col = trunc(Int,size(y,2)/2)
    return y[:,1:col],y[:,(col+1):end]
end

function _gain_(model::XGBoostRegressionTree, y, y_pred)
    grad = sum(model.loss.gradient(y, y_pred))
    nominator = grad .* grad
    denominator = sum(model.loss.hessian(y, y_pred))
    return 0.5 * nominator / denominator
end

function impurity_calc(model::XGBoostRegressionTree, y, y1, y2)
    y, y_pred = split_(y)
    y1, y1_pred = split_(y1)
    y2, y2_pred = split_(y2)
    true_gain = _gain_(model, y1,y1_pred)
    false_gain = _gain_(model, y2,y2_pred)
    gain = _gain_(model, y, y_pred)
    return true_gain + false_gain - gain
end

function impurity_calc(model::RegressionTree, y, y1, y2)
    var_total = var(y)
    var_y1 = var(y1)
    var_y2 = var(y2)
    frac_1 = length(y1) / length(y)
    frac_2 = length(y2) / length(y)

    variance_reduction = var_total - (frac_1 * var_y1 + frac_2 * var_y2)
    return variance_reduction
end

function impurity_calc(model::ClassificationTree, y, y1, y2)
    p = size(y1,1)/size(y,1)
    entro = calc_entropy(y)
    #println("entro: $(entro), p: $(p), res: $(entro - (p * calc_entropy(y1) + (1-p) * calc_entropy(y2)))")
    return entro - (p * calc_entropy(y1) + (1-p) * calc_entropy(y2))
end

function predict(model::DecisionTree,
                 x::Matrix)
    n = size(x,1)
    res = zeros(n)
    if model.y_num == 1
        for i = 1:n
            res[i] = predict(model.root, x[i,:])
        end
    else
        res = zeros(n, model.y_num)
        for i = 1:n
            res[i,:] = predict(model.root, x[i,:])
        end
    end
    return res
end

function predict(model::DecisionNode,
                 x::Vector)
    if model.label != "nothing"
        return model.label
    end
    feature_current = x[model.feature_index]
    if typeof(feature_current) <:String
        if feature_current == model.threshold
            return predict(model.true_branch, x)
        else
            return predict(model.false_branch, x)
        end
    elseif typeof(feature_current) <: Real
        if feature_current <= model.threshold
            return predict(model.true_branch, x)
        else
            return predict(model.false_branch, x)
        end
    end
end

function print_tree(model::DecisionTree)
    if model.lable != "nothing"
        println($(model.label))
    else
        println("$(model.feature_index):$(model.threshold)?")
        println(" T->")
        print_tree(model.true_branch)
        println(" F->")
        print_tree(model.true_branch)
    end
end

function split_at_feature(X, feature_i, threshold)
    if typeof(threshold) <: Real
        ind1 = find(X[:,feature_i] .<= threshold)
        ind2 = find(X[:,feature_i] .> threshold)
    elseif typeof(threshhold) <: String
        ind1 = find(X[:,feature_i] .== threshold)
        ind2 = find(X[:,feature_i] .!= threshold)
    end
    X_1 = X[ind1,:]
    X_2 = X[ind2,:]
    return X_1, X_2
end

function test_ClassificationTree()
    X_train, X_test, y_train, y_test = make_iris()
    y_train = one_hot(y_train)
    y_test = one_hot(y_test)
    model = ClassificationTree()
    train!(model,X_train, y_train)
    predictions = predict(model,X_test)
    y_test = unhot(y_test)
    predictions = unhot(predictions)
    print("classification accuracy", accuracy(y_test, predictions))
    #PCA

    #pca_model = PCA()
    #train!(pca_model, X_test)
    #plot_in_2d(pca_model, X_test, predictions, "ClassificationTree")
end

function test_RegressionTree()
    X_train, X_test, y_train, y_test = make_reg(n_features = 1)
    model = RegressionTree()
    train!(model,X_train, y_train)
    predictions = predict(model,X_test)
    print("regression msea", mean_squared_error(y_test, predictions))
    PyPlot.scatter(X_test, y_test, color = "black")
    PyPlot.scatter(X_test, predictions, color = "green")
    legend(loc="upper right",fancybox="true")
end


text:
"""
for constructing boundary layers and evaluating [`PML`](@ref)/[`cPML`](@ref) conductivity
"""
module BoundaryLayers

export PML
export cPML
export noBL
export AbstractBL
export AbstractRealBL
export AbstractComplexBL
export getside

# PML params used in Boundaries
import ..EXTINCTION
import ..SCALING_ANGLE
import ..BL_DEPTH

import ..PRINTED_COLOR_NUMBER
import ..PRINTED_COLOR_DARK

using ..Points
using Formatting

"""
	AbstractBL{SIDE}
"""
abstract type AbstractBL{SIDE} end
"""
	AbstractRealBL{SIDE} <: AbstractBL{SIDE}
"""
abstract type AbstractRealBL{SIDE} <: AbstractBL{SIDE} end
"""
	AbstractComplexBL{SIDE} <: AbstractBL{SIDE}
"""
abstract type AbstractComplexBL{SIDE} <: AbstractBL{SIDE} end

"""
	getside(::AbstractBL) = side
"""
getside(::AbstractBL{SIDE}) where SIDE = SIDE

for (TBL,ATBL) ∈ zip((:PML,:cPML,:noBL),
					 (:AbstractComplexBL,:AbstractComplexBL,:AbstractRealBL))
	@eval begin
		"""
			$($TBL){SIDE}

		Fields: `depth`, `start`, `stop`
		"""
		mutable struct $TBL{SIDE} <: $ATBL{SIDE}
			depth::Float64
		    start::Float64
		    stop::Float64

			function $TBL{SIDE}(depth,start,stop) where SIDE
				if !($TBL <: noBL)
					depth>0 || throw(ArgumentError("cannot have $($TBL) with 0 depth"))
				end
				return new{SIDE}(depth,start,stop)
			end
		end

		if $TBL == :noBL
			noBL{SIDE}(depth::Real=BL_DEPTH) where SIDE = noBL{SIDE}(0,NaN,NaN)
		else
			"""
				$($TBL){SIDE}(depth) -> bl

			$($TBL) along dimension `SIDE` with `depth`. The location of the layer determined
			later by associating it with a shape in a call to `Boundary`.

			see also: [`Boundary`](@ref), [`PML`](@ref), [`cPML`](@ref), [`noBL`](@ref)
			"""
			$TBL{SIDE}(depth::Real=BL_DEPTH) where SIDE = $TBL{SIDE}(depth,NaN,NaN)
		end

		function Base.show(io::IO,tbl::$TBL{SIDE}) where SIDE
			printstyled(io,$TBL,color=PRINTED_COLOR_DARK)
			!get(io,:compact,false) ? print(io," SIDE=",SIDE) : nothing
			if $TBL!==noBL
				print(io," depth=")
				printstyled(io,fmt("1.1f",tbl.depth),color=PRINTED_COLOR_NUMBER)
			end
		end

		(tbl::$TBL)(x::Real,y...) = tbl(Point(x,y...))

		(bl::$TBL{SIDE})(;depth=bl.depth, start=bl.start, stop=bl.stop) where SIDE = $TBL{SIDE}(depth,start,stop)
	end
end
# noBL{SIDE}(depth::Real=BL_DEPTH) where SIDE = noBL{SIDE}(0,NaN,NaN)

"""
	(::PML)(p::Point) -> σ

evaluate the PML conductivity `σ` at the point `p`.
"""
(pml::PML)(p::Point) = _conductivity_profile(p[1+(getside(pml)-1)÷2],pml.start,pml.stop)

"""
	(::cPML)(p::Point) -> σ

evaluate the cPML conductivity `σ` at the point `p`.
"""
(cpml::cPML)(p::Point) = _conj_conductivity_profile(p[1+(getside(cpml)-1)÷2],cpml.start,cpml.stop)

"""
	(::noBL)(p::Point) -> 0.0
"""
(nbl::noBL)(p::Point) = 0.0

"""
	conj(pml) -> cpml
"""
Base.conj(pml::PML{SIDE}) where SIDE = cPML{SIDE}(pml.depth,pml.start,pml.stop)
"""
	conj(cpml) -> pml
"""
Base.conj(cpml::cPML{SIDE}) where SIDE = PML{SIDE}(cpml.depth,cpml.start,cpml.stop)
"""
	conj(nobl) -> nobl
"""
Base.conj(nbl::noBL) = nbl

# used inside PML. This is the place to modify the PML profile
"""
	_conductivity_profile(z,start,stop) -> σ

PML conductivity `σ` at point `z`
"""
function _conductivity_profile(z::Real,start::Real,stop::Real)
	if (start≤stop && start≤z) || (start≥stop && start≥z)
		depth = abs(stop-start)
	    β = 2*sqrt(depth)
	    α = -(1/4)*exp(complex(0,SCALING_ANGLE))*log(EXTINCTION)/depth/((2exp(β)-β)/(2β*expm1(β))-1/β^2)
        u = abs(z-start)/depth
        return α*u*expm1(β*u)/expm1(β)
    else
        return complex(0.0)
    end
end

"""
	_integrated_conductivity_profile(z,start,stop) -> σ

integrated PML conductivity `σ` at point `z`
"""
function _integrated_conductivity_profile(z::Real,start::Real,stop::Real)
	if (start≤stop && start≤z) || (start≥stop && start≥z)
		depth = abs(stop-start)
		β = 2*sqrt(depth)
		α = -(1/4)*exp(complex(0,SCALING_ANGLE))*log(EXTINCTION)/depth/((2exp(β)-β)/(2β*expm1(β))-1/β^2)
		u = abs(z-start)/depth
		return (α*depth/expm1(β))*(u*exp(β*u)-(1+1/β)*expm1(β*u)/β)
	else
		return complex(0.0)
	end
end

# used inside cPML
"""
	_conj_conductivity_profile(z,start,stop) -> σ

conjugate PML conductivity `σ` at point `z`
"""
_conj_conductivity_profile(z,start,stop) = -conj(_conductivity_profile(z,start,stop))

"""
	_conj_integrated_conductivity_profile(z,start,stop) -> σ

integrated conjugate PML conductivity `σ` at point `z`
"""
_conj_integrated_conductivity_profile(z,start,stop) = -conj(_integrated_conductivity_profile(z,start,stop))

end # module


text:
using ForwardDiff
using ReverseDiff
using Measures
include("helper_functions.jl")

begin
    λ1 = 0.4
    λ2 = 1.5
    α = 1
    μ(α) = 1 / (1.5 * (α))
    ν(α) = 1 / (0.7 * (α))
    parameters = (λ1, λ2, α, μ, ν)
    opponent_parameters = (λ2, λ1, α, μ, ν)
    p, q = equilibrium_probability(parameters...)
    step = 0.0005
    xs = 0:step:1
    ys = 0:step:1
end

begin
    threeD_cost = (p, q) -> total_cost(p, q, parameters...)
    heatcost = plot(xs, ys, threeD_cost, st = :contourf)
    scattersize = 5
    plot!([p], [q], st = :scatter, shape = :circle, msc = :black, ms = scattersize, label = "Equilibrium")
    plot!(xlabel = "Player 1 probability", ylabel = "Player 2 probability")

    begin
        gtape = ReverseDiff.GradientTape(vec -> threeD_cost(vec...), [p, q])
        out = zeros(2)
        dtotal(ϕ, ψ) = ReverseDiff.gradient!(out, gtape, [ϕ, ψ])
    end

    difference =  1
    p1, q1 = p, q
    while difference > 1e-6
        p0, q0 = p1, q1
        proots = first(find_valid_roots(a -> first(dtotal(a, q0))))
        p1 = proots[findmin(abs.(proots .- p0))[2]]
        qroots = first(find_valid_roots(b -> last(dtotal(p1, b))))
        q1 = qroots[findmin(abs.(qroots .- q0))[2]]
        difference = maximum([abs(p1 - p0), abs(q1 - q0)])
    end

    multiplier = 0.1
    plot_xlims = p1 > p ? ((1-multiplier) * p, (1+multiplier) * p1) : ((1-multiplier) * p1, (1+multiplier) * p) #(0.4, 0.45)
    plot_ylims = q1 > q ? ((1-multiplier) * q, (1+multiplier) * q1) : ((1-multiplier) * q1, (1+multiplier) * q) #(0.58, 0.62)
    plot_clims = (threeD_cost(p1, q1)*(1-multiplier), threeD_cost(p, q)*(1+multiplier))

plot!(heatcost, [p1], [q1], st = :scatter, ms = scattersize, label = "Minimum Response Time")
plot!(colorbar_title = "Sum of Average Response Time")
plot!(size = (600, 300), margin = 1.5mm)

heatcost
begin
    possibleroots = first(find_valid_roots(q -> ForwardDiff.derivative(ψ -> total_cost(pstar(ψ, parameters...), ψ, parameters...), q)))
    roots = possibleroots[ForwardDiff.derivative.(q -> ForwardDiff.derivative(ψ -> total_cost(pstar(ψ, parameters...), ψ, parameters...), q), possibleroots) .> 0]
    q_p1br = roots[findmin(abs.(roots .- q))[2]]
    p_p1br = pstar(q_p1br, parameters...)
    heatcost = plot(heatcost, [p_p1br], [q_p1br], st = :scatter, ms = scattersize, label = "Min where A1 best responds")
end

begin
    possibleroots = first(find_valid_roots(p -> ForwardDiff.derivative(ψ -> total_cost(ψ, qstar(ψ, parameters...), parameters...), p)))
    roots = possibleroots[ForwardDiff.derivative.(p -> ForwardDiff.derivative(ψ -> total_cost(ψ, qstar(ψ, parameters...), parameters...), p), possibleroots) .> 0]
    p_p2br = roots[findmin(abs.(roots .- p))[2]]
    q_p2br = qstar(p_p2br, parameters...)
    plot!(heatcost, [p_p2br], [q_p2br], st = :scatter, ms = scattersize, label = "Min where A2 best responds")
end

savefig("outputs/total_cost_contour_poster.pdf")

probability_path = he_equilibrium_probability.(λ1, λ2, α, μ, ν, 0:0.001:1)

begin
    heatcost_path = plot(heatcost, legend = :topright, size = (500, 300), xlims = plot_xlims, ylims = plot_ylims)
    n = 5
    for i in 1:n
        interval = Int(floor(length(probability_path) / n))
        slice = 1 + (i - 1) * interval:i * interval
        plot!(heatcost_path,
            first.(probability_path[slice]),
            last.(probability_path[slice]),
            lc = :black,
            lw = 1,
            ls = :solid,
            arrow = i == n ? false : true,
            label = i == 1 ? "Homo-Egualis path" : "")
        # if i != n
        #     annotate!(
        #         heatcost_path,
        #         probability_path[i*interval + 1]...,
        #         text(
        #             round(collect(0:0.001:1)[i*interval], digits = 2),
        #             :bottom,
        #             :green,
        #             11)
        #         )
        # end
    end
    heatcost_path
end

# savefig("outputs/homo-egualis_path_contour.pdf")