text:
using Test, StaticArrays
using CLIMA.VariableTemplates: varsindex, @vars
using StaticArrays

struct MoistureModel{FT} end
struct AtmosModel{FT}
    moisture::MoistureModel{FT}
end

function vars_state(::MoistureModel, FT)
    @vars begin
        ρq_tot::FT
        ρq_x::SVector{5, FT}
        ρq_liq::FT
        ρq_vap::FT
    end
end
function vars_state(m::AtmosModel, FT)
    @vars begin
        ρ::FT
        ρu::SVector{3, FT}
        ρe::FT
        moisture::vars_state(m.moisture, FT)
    end
end
@testset "Varsindex" begin
    FT = Float64
    m = AtmosModel(MoistureModel{FT}())

    @test 1:1 === varsindex(vars_state(m, FT), :ρ)
    @test 2:4 === varsindex(vars_state(m, FT), :ρu)
    @test 5:5 === varsindex(vars_state(m, FT), :ρe)

    # Since moisture is defined recusively this will get all the fields
    moist = varsindex(vars_state(m, FT), :moisture)
    @test 6:13 === moist

    # To get the specific onnes we need something like
    @test 6:6 === moist[varsindex(vars_state(m.moisture, FT), :ρq_tot)]
    @test 7:11 === moist[varsindex(vars_state(m.moisture, FT), :ρq_x)]
    @test 12:12 === moist[varsindex(vars_state(m.moisture, FT), :ρq_liq)]
    @test 13:13 === moist[varsindex(vars_state(m.moisture, FT), :ρq_vap)]
end


text:
# FixedPenalization penalty
# ------------------------------------------------------------
"""
    FixedPenalization

Functor object for using fixed penalty in a progressive-hedging algorithm. Create by supplying a [`Fixed`](@ref) object through `penalty` in the `ProgressiveHedgingSolver` factory function and then pass to a `StochasticPrograms.jl` model.

...
# Parameters
- `r::T = 1.00`: Fixed penalty
...
"""
mutable struct FixedPenalization{T <: AbstractFloat} <: AbstractPenalization
    r::T

    function FixedPenalization(r::AbstractFloat)
        T = typeof(r)
        return new{T}(r)
    end
end
function penalty(::AbstractProgressiveHedging, penalty::FixedPenalization)
    return penalty.r
end
function initialize_penalty!(::AbstractProgressiveHedging, ::FixedPenalization)
    nothing
end
function update_penalty!(::AbstractProgressiveHedging, ::FixedPenalization)
    nothing
end

# API
# ------------------------------------------------------------
"""
    Fixed

Factory object for [`FixedPenalization`](@ref). Pass to `penalty` in the `ProgressiveHedgingSolver` factory function. See ?FixedPenalization for parameter descriptions.

"""
struct Fixed{T <: AbstractFloat} <: AbstractPenalizer
    r::T

    function Fixed(; r::AbstractFloat = 1.0)
        T = typeof(r)
        return new{T}(r)
    end
end

function (fixed::Fixed)()
    return FixedPenalization(fixed.r)
end

function str(::Fixed)
    return "fixed penalty"
end

ext: jl


text:
export AbstractProbabilityVector, AbstractUnitIntervalThresholds, AbstractRealLineThresholds
export probabilities, thresholds, deformation, outcome, deformation_value, outcome_value, meanvalue
export ProbabilityVector, UnitIntervalThresholds, RealLineThresholds, DeformedRealLineThresholds

"""
    AbstractProbabilityVector{T}

Abstract type for probability vectors with entries of type `T` summing to one.
"""
abstract type AbstractProbabilityVector{T} end

"""
    probabilities(apv::AbstractProbabilityVector{T})

Get the probabilities of `apv` in a raw format, typically an `AbstractVector{T}`.
"""
function probabilities(apv::AbstractProbabilityVector) end

"""
    AbstractUnitIntervalThresholds{T}

Abstract type for a weakly sorted vector of thresholds of type `T` where values are between zero and one.
"""
abstract type AbstractUnitIntervalThresholds{T} end

"""
    thresholds(auit::AbstractUnitIntervalThresholds{T}) where {T}

Get the thresholds of `auit` in a raw format, typically an `AbstractVector{T}`.
"""
function thresholds(auit::AbstractUnitIntervalThresholds) end

"""
    AbstractRealLineThresholds{T}

Abstract type for a weakly sorted vector of thresholds of type `T` where values are unconstrained
"""
abstract type AbstractRealLineThresholds{T} end

"""
    thresholds(arlt::AbstractRealLineThresholds{T}) where {T}

Get the thresholds of `arlt` in a raw format, typically an `AbstractVector{T}`.
"""
function thresholds(arlt::AbstractRealLineThresholds{T}) where {T} end

"""
    ProbabilityVector{T, TVec <: AbstractVector{T}} <: AbstractProbabilityVector{T}

Trusted vector of probabilities of type `T` that sum to one.

    ProbabilityVector(ps::TVec) where {T, TVec<:AbstractVector{T}}

Convert a vector of probabilities to a `ProbabilityVector`, thereby ensuring that probabilities
are non-negative and sum to one (trusted). Regarding the last requirement, this constructor is very unforgiving.

    ProbabilityVector(apv::AbstractProbabilityVector{T}) where {T}
    ProbabilityVector(auit::AbstractUnitIntervalThresholds{T}) where {T}

Convert to `ProbabilityVector` other types that are trusted.
"""
struct ProbabilityVector{T, TVec <: AbstractVector{T}} <: AbstractProbabilityVector{T}
    probabilities::TVec

    function ProbabilityVector(ps::TVec) where {T, TVec<:AbstractVector{T}}
        if length(ps) == 0
            throw(DomainError(ps, "Vector of probabilities was empty so could not possibly sum to one."))
        end
        if sum(ps) != one(T)
            throw(DomainError(ps, "Vector of probabilities did not sum to one."))
        end
        for p in ps
            if p < zero(T)
                throw(DomainError(ps, "At least one probability was negative"))
            end
        end
        return new{T, TVec}(ps)
    end

    function ProbabilityVector(apv::AbstractProbabilityVector{T}) where {T}
        p = probabilities(apv)
        return new{T, typeof(p)}(p)
    end

    function ProbabilityVector(auit::AbstractUnitIntervalThresholds{T}) where {T}
        thrs = thresholds(auit)
        n = length(thrs) + 1
        probabilities = Vector{T}(undef, n)
        @inbounds begin
            if n == 1
                probabilities[1] = one(T) - zero(T)
            else
                probabilities[1] = thrs[1] - zero(T)
                for i in 2:(n-1)
                    probabilities[i] = thrs[i] - thrs[i-1]
                end
                probabilities[n] = one(T) - thrs[n-1]
            end
        end
        return new{T, Vector{T}}(probabilities)
    end
end

function probabilities(pv::ProbabilityVector)
    return pv.probabilities
end

"""
    UnitIntervalThresholds{T, TVec <: AbstractVector{T}} <: AbstractUnitIntervalThresholds{T}

A weakly sorted vector of values on the interval from zero and one.

The `UnitIntervalThresholds` struct is used to check that values are in interval and weakly sorted, and using
the type system to express that.

    UnitIntervalThresholds(thr::TVec) where {T, TVec <: AbstractVector{T}}

Convert a vector `thr` of thresholds of type `T` to a UnitIntervalThresholds

    UnitIntervalThresholds(auit::AbstractUnitIntervalThresholds{T})

Convert an AbstractUnitIntervalThresholds `auit` to a UnitIntervalThresholds by extracting its thresholds

    UnitIntervalThresholds(apv::AbstractProbabilityVector{T}) where {T}

Convert an `AbstractProbabilityVector` to a `UnitIntervalThresholds`

    UnitIntervalThresholds(arlt::AbstractRealLineThresholds{TReal}, atl::AbstractThresholdLink{TUnit, TReal}) where {TUnit, TReal}

Convert an AbstractRealLineThresholds `arlt` to a UnitIntervalThresholds using the threshold link `atl`
"""
struct UnitIntervalThresholds{T, TVec <: AbstractVector{T}} <: AbstractUnitIntervalThresholds{T}
    thresholds::TVec

    function UnitIntervalThresholds(thr::TVec) where {T, TVec <: AbstractVector{T}}
        if length(thr) == 0
            return new{T, TVec}(thr)
        end
        if first(thr) < zero(T)
            throw(DomainError(thr, "Unit interval thresholds had a value less than zero."))
        end
        if last(thr) > one(T)
            throw(DomainError(thr, "Unit interval thresholds had a value greater than one."))
        end
        if !issorted(thr)
            throw(DomainError(thr, "Unit interval thresholds were not weakly sorted."))
        end
        return new{T, TVec}(thr)
    end

    function UnitIntervalThresholds(auit::AbstractUnitIntervalThresholds{T}) where {T}
        thrs = thresholds(auit)
        return new{T, typeof(thrs)}(thrs)
    end

    function UnitIntervalThresholds(apv::AbstractProbabilityVector{T}) where {T}
        probs = probabilities(apv)
        n = length(probs)
        thresholds = Vector{T}(undef, n - 1)

        accumulated_probability = zero(T)
        @inbounds for i in 1:(n-1)
            accumulated_probability += probs[i]
            thresholds[i] = accumulated_probability
        end
        return new{T, Vector{T}}(thresholds)
    end

    function UnitIntervalThresholds(arlt::AbstractRealLineThresholds{TReal}, atl::AbstractThresholdLink{TUnit, TReal}) where {TUnit, TReal}
        return new{TUnit, Vector{TUnit}}(cdf(atl).(thresholds(arlt)))
    end
end

function thresholds(uit::UnitIntervalThresholds)
    return uit.thresholds
end

"""
    RealLineThresholds{T, TVec <: AbstractVector{T}} <: AbstractRealLineThresholds{T}

A weakly sorted vector of values on the real line.

The UnitIntervalThresholds struct is used to check that values are weakly sorted, and using
    the type system to express that.

    RealLineThresholds(rlt::TVec) where {T, TVec <: AbstractVector{T}}

Convert a raw vector of thresholds `rlt` of type `T` to a RealLineThresholds

    RealLineThresholds(arlt::AbstractRealLineThresholds{T}) where {T}

Convert an AbstractRealLineThresholds `arlt` to a RealLineThresholds by extracting the raw thresholds

    RealLineThresholds(auit::AbstractUnitIntervalThresholds{TUnit}, atl::AbstractThresholdLink{TUnit, TReal}) where {TUnit, TReal}

Convert an AbstractUnitIntervalThresholds `auit` to a RealLineThresholds using the link `atl`
"""
struct RealLineThresholds{T, TVec <: AbstractVector{T}} <: AbstractRealLineThresholds{T}
    thresholds::TVec

    function RealLineThresholds(rlt::TVec) where {T, TVec <: AbstractVector{T}}
        if (!issorted(rlt, lt=(<=)))
            throw(DomainError(ps, "Thresholds were not weakly sorted"))
        end
        return new{T, TVec}(rlt)
    end

    function RealLineThresholds(arlt::AbstractRealLineThresholds{T}) where {T}
        thrs = thresholds(arlt)
        return new{T, typeof(thrs)}(thrs)
    end

    function RealLineThresholds(auit::AbstractUnitIntervalThresholds{TUnit}, atl::AbstractThresholdLink{TUnit, TReal}) where {TUnit, TReal}
        return new{TReal, Vector{TReal}}(invcdf(atl).(thresholds(auit)))
    end
end

function thresholds(rlt::RealLineThresholds)
    return rlt.thresholds
end

"""
    DeformedRealLineThresholds{T, TBase <: AbstractRealLineThresholds{T}} <: AbstractRealLineThresholds{T}

A weakly sorted vector of values on the real line that has been affinely deformed.

DeformedRealLineThresholds holds the specification of the affine deformation and defers calculation
of the deformed thresholds until they are needed.
"""
struct DeformedRealLineThresholds{T, TBase <: AbstractRealLineThresholds{T}} <: AbstractRealLineThresholds{T}
    undeformed::TBase
    push::T
    norm::T
end

function thresholds(drlt::DeformedRealLineThresholds{T}) where {T}
    return (thresholds(drlt.undeformed) .+ drlt.push) ./ drlt.norm
end

# Operations on the above types

"""
Make a deformation of probabilities or thresholds

    deformation(arlt::AbstractRealLineThresholds, push, norm)
    deformation(auit::AbstractUnitIntervalThresholds, atl::AbstractThresholdLink, push, norm)
    deformation(apv::AbstractProbabilityVector, atl::AbstractThresholdLink, push, norm)

Deformations are made on the thresholds on the real line. They consist of a `push` that moves the thresholds
and a `norm` with which the pushed thresholds are divided for proper normalization.

Methods on probabilities and unit interval thresholds are provided for ease of use and also depend on a threshold link
`atl`.
"""
function deformation end

function deformation(arlt::AbstractRealLineThresholds, push, norm)
    return DeformedRealLineThresholds(arlt, push, norm)
end

function deformation(auit::AbstractUnitIntervalThresholds, atl::AbstractThresholdLink, push, norm)
    return UnitIntervalThresholds(deformation(RealLineThresholds(auit, atl), push, norm), atl)
end

function deformation(apv::AbstractProbabilityVector, atl::AbstractThresholdLink, push, norm)
    return ProbabilityVector(deformation(UnitIntervalThresholds(apv), atl, push, norm))
end

"""
Compute the outcome of a particular realization of scenario
"""
function outcome end

function outcome(ths, z)
    return searchsortedfirst(thresholds(ths), z)
end

function outcome(ths, zs::AbstractVector)
    return [outcome(ths, z) for z in zs]
end

"""
Compute the mean value when given the value for each outcome.
"""
function meanvalue end

function meanvalue(apv::AbstractProbabilityVector, values::AbstractVector)
    ps = probabilities(apv)
    n = length(ps)
    if n != length(values)
        throw(DomainError(values, "Values did not have the same length as probablity vector"))
    end
    return sum(ps .* values)
end

function meanvalue(auit::AbstractUnitIntervalThresholds{T}, values::AbstractVector) where {T}
    thrs = thresholds(auit)
    n = length(thrs) + 1
    if n != length(values)
        throw(DomainError(values, "Values did not have length one greater than thresholds"))
    end
    @inbounds begin
        if n == 1
            return (one(T) - zero(T)) * values[1]
        else
            rv = (first(thrs) - zero(T)) * values[1]
            for i in 2:(n-1)
                rv += (thrs[i] - thrs[i-1]) * values[i]
            end
            rv += (one(T) - last(thrs)) * values[n]
            return rv
        end
    end
end

"""
Compute the mean value after a deformation when the values for each outcome is given.
"""
function deformation_value end

function deformation_value(auit::AbstractUnitIntervalThresholds, values::AbstractVector, atl::AbstractThresholdLink, push, norm)
    return meanvalue(deformation(auit, atl, push, norm), values)
end

function deformation_value(arlt::AbstractRealLineThresholds, values::AbstractVector, atl::AbstractThresholdLink, push, norm)
    return meanvalue(UnitIntervalThresholds(deformation(arlt, push, norm), atl), values)
end

function deformation_value(arlt::AbstractRealLineThresholds, values::AbstractVector, atl::AbstractThresholdLink, pushs::AbstractVector, norm)
    return [deformation_value(arlt, values, atl, push, norm) for push in pushs]
end

"""
Compute the values for one or more outcomes.
"""
function outcome_value end

function outcome_value(auit::AbstractUnitIntervalThresholds, values::AbstractVector, p)
    return values[outcome(auit, p)]
end

function outcome_value(auit::AbstractUnitIntervalThresholds, values::AbstractVector, ps::AbstractArray)
    return [outcome_value(auit, values, p) for p in ps]
end

function outcome_value(auit::AbstractUnitIntervalThresholds, values::AbstractVector, atl::AbstractThresholdLink, z)
    return outcome_value(auit, values, cdf(atl)(z))
end

function outcome_value(auit::AbstractUnitIntervalThresholds, values::AbstractVector, atl::AbstractThresholdLink, zs::AbstractArray)
    return [outcome_value(auit, values, cdf(atl)(z)) for z in zs]
end

function outcome_value(apv::AbstractProbabilityVector, values::AbstractVector, p)
    return outcome_value(UnitIntervalThresholds(apv), values, p)
end

function outcome_value(arlt::AbstractRealLineThresholds, values::AbstractVector, z)
    return values[outcome(arlt, z)]
end

function outcome_value(arlt::AbstractRealLineThresholds, values::AbstractVector, zs::AbstractVector)
    return [outcome_value(arlt, values, z) for z in zs]
end


text:
using ImageInpainting
using ColorTypes
using TestImages
using Plots; gr(size=(800,400),yflip=true,colorbar=false,ticks=false)
using VisualRegressionTests
using Test, Pkg

# list of maintainers
maintainers = ["juliohm"]

# environment settings
istravis = "TRAVIS" ∈ keys(ENV)
ismaintainer = "USER" ∈ keys(ENV) && ENV["USER"] ∈ maintainers
datadir = joinpath(@__DIR__,"data")

if ismaintainer
  Pkg.add("Gtk")
  using Gtk
end

# test images
blobs = testimage("blobs")
lighthouse = testimage("lighthouse")

# helper functions
function plot_crimisini_on_array(fname, img, inds, psize)
  fimg = Float64.(Gray.(img))
  mask = falses(size(fimg))
  mask[inds...] .= true
  fimg[mask] .= NaN
  fimg2 = inpaint(fimg, mask, Crimisini(psize))
  plt1 = heatmap(fimg, title="before inpainting")
  plt2 = heatmap(fimg2, title="after inpainting")
  plot(plt1, plt2)
  png(fname)
end

@testset "ImageInpainting.jl" begin
  @testset "Plain arrays" begin
    plot_blobs(fname) = plot_crimisini_on_array(fname, blobs, (50:150,50:150), (11,11))
    refimg = joinpath(datadir,"Blobs.png")
    @test test_images(VisualTest(plot_blobs, refimg), popup=!istravis) |> success

    plot_lighthouse(fname) = plot_crimisini_on_array(fname, lighthouse, (50:350,300:400), (30,30))
    refimg = joinpath(datadir,"LightHouse.png")
    @test test_images(VisualTest(plot_lighthouse, refimg), popup=!istravis) |> success
  end
end


text: if true else end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
include("../../example_cases/new node types/plotting.jl")
include("../../example_cases/new node types/ThirdOrderEq.jl")
include("../../example_cases/new node types/ThirdOrderEqGovernorIEEEG1.jl")
include("../../example_cases/new node types/ThirdOrderEqExciterIEEEDC1A.jl")

line_list = []
append!(line_list, [StaticLine(from = 1, to = 2, Y = -1im / 0.02)])
append!(line_list, [StaticLine(from = 2, to = 3, Y = -1im / 0.02)])
append!(line_list, [StaticLine(from = 3, to = 4, Y = -1im / 0.02)])
append!(line_list, [StaticLine(from = 4, to = 5, Y = -1im / 0.02)])
append!(line_list, [StaticLine(from = 5, to = 6, Y = -1im / 0.02)])
append!(line_list, [StaticLine(from = 1, to = 6, Y = -1im / 0.02)])

#=
append!(line_list, [StaticLine(from = 6, to = 7, Y = -1im / 0.02)])
append!(line_list, [StaticLine(from = 7, to = 8, Y = -1im / 0.02)])
append!(line_list, [StaticLine(from = 8, to = 9, Y = -1im / 0.02)])
append!(line_list, [StaticLine(from = 9, to = 10, Y = -1im / 0.02)])
append!(line_list, [StaticLine(from = 10, to = 11, Y = -1im / 0.02)])
append!(line_list, [StaticLine(from = 1, to = 11, Y = -1im / 0.02)])
append!(line_list, [StaticLine(from = 2, to = 7, Y = -1im / 0.02)])
append!(line_list, [StaticLine(from = 1, to = 4, Y = -1im / 0.02)])
append!(line_list, [StaticLine(from = 2, to = 11, Y = -1im / 0.02)])
append!(line_list, [StaticLine(from = 2, to = 9, Y = -1im / 0.02)])
=#

node_list = []

append!(node_list, [SlackAlgebraic(U = 1, Y_n = 0)])
#append!(node_list, [PQAlgebraic(S = 0.9, Y_n = 1)])
append!(node_list, [ThirdOrderEq(H = 3.318, P = -0.6337, D = 0.1, Ω = 50, E_f = 0.5, T_d_dash = 8.690, X_q_dash = 0.103, X_d_dash = 0.111, X_d = 0.1)])
#append!(node_list, [PQAlgebraic(S = -0.5, Y_n = 1)])
append!(node_list, [ThirdOrderEq(H = 3.318, P = -0.6337, D = 0.1, Ω = 50, E_f = 0.5, T_d_dash = 8.690, X_q_dash = 0.103, X_d_dash = 0.111, X_d = 0.1)])
#append!(node_list, [PQAlgebraic(S = 0.9, Y_n = 1)])
append!(node_list, [ThirdOrderEq(H = 3.318, P = 0.6337, D = 0.1, Ω = 50, E_f = 0.5, T_d_dash = 8.690, X_q_dash = 0.103, X_d_dash = 0.111, X_d = 0.1)])
#append!(node_list, [PQAlgebraic(S = 0.9, Y_n = 1)])
append!(node_list, [ThirdOrderEq(H = 3.318, P = 0.6337, D = 0.1, Ω = 50, E_f = 0.5, T_d_dash = 8.690, X_q_dash = 0.103, X_d_dash = 0.111, X_d = 0.1)])
#append!(node_list, [PQAlgebraic(S = 0.9, Y_n = 1)])
append!(node_list, [ThirdOrderEq(H = 3.318, P = 0.6337, D = 0.1, Ω = 50, E_f = 0.5, T_d_dash = 8.690, X_q_dash = 0.103, X_d_dash = 0.111, X_d = 0.1)])

#=
append!(node_list, [SlackAlgebraic(U = 1, Y_n = 1)])
append!(node_list, [ThirdOrderEqGovernorIEEEGone(H = 5.0, D = 0.1, Ω = 50, E_f = 0.2, T_d_dash = 0.1, X_d_dash = 0.5,  X_q_dash = 0.103, X_d = 0.5, P0 = 0.5, Pmax = 1.03, Pmin = 0, Pup = 1.38, Pdown = -1.25, T1 = 0.82, T2 = 0, T3 = 0.97, K = 20.0)])
append!(node_list, [ThirdOrderEqGovernorIEEEGone(H = 5.0, D = 0.1, Ω = 50, E_f = 0.2, T_d_dash = 0.1, X_d_dash = 0.5,  X_q_dash = 0.103, X_d = 0.5, P0 = 0.5, Pmax = 1.03, Pmin = 0, Pup = 1.38, Pdown = -1.25, T1 = 0.82, T2 = 0, T3 = 0.97, K = 20.0)])
append!(node_list, [ThirdOrderEqGovernorIEEEGone(H = 5.0, D = 0.1, Ω = 50, E_f = 0.2, T_d_dash = 0.1, X_d_dash = 0.5,  X_q_dash = 0.103, X_d = 0.5, P0 = 0.5, Pmax = 1.03, Pmin = 0, Pup = 1.38, Pdown = -1.25, T1 = 0.82, T2 = 0, T3 = 0.97, K = 20.0)])
append!(node_list, [ThirdOrderEqGovernorIEEEGone(H = 5.0, D = 0.1, Ω = 50, E_f = 0.2, T_d_dash = 0.1, X_d_dash = 0.5,  X_q_dash = 0.103, X_d = 0.5, P0 = 0.5, Pmax = 1.03, Pmin = 0, Pup = 1.38, Pdown = -1.25, T1 = 0.82, T2 = 0, T3 = 0.97, K = 20.0)])
append!(node_list, [ThirdOrderEqGovernorIEEEGone(H = 5.0, D = 0.1, Ω = 50, E_f = 0.2, T_d_dash = 0.1, X_d_dash = 0.5,  X_q_dash = 0.103, X_d = 0.5, P0 = 0.5, Pmax = 1.03, Pmin = 0, Pup = 1.38, Pdown = -1.25, T1 = 0.82, T2 = 0, T3 = 0.97, K = 20.0)])

=#

powergrid = PowerGrid(node_list, line_list)

return powergrid

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


text:
#unit-constants.jl
#testing each of the unum constant generators in four unum environments:
# Unum{0,0}, Unum{1,1}, Unum{4,6}, Unum{4,8}

#it's best to do this with a macro
macro uctest(T, f, ex, fr, fl, es, fs)
  creationstmt = :($f($T))

  arraytest = T.args[3] > 6 ? (:(x.fraction.a == $fr)) : (:(x.fraction == $fr))
  quote
    x = $creationstmt

    @test x.fsize == $fs
    @test x.esize == $es
    @test x.flags == $fl
    @test $arraytest
    @test x.exponent == $ex
  end
end

nan(T) = T(NaN)
inf(T) = Unums.inf(T)

@uctest Unum{0,0} zero z64 z64 z16 z16 z16
@uctest Unum{1,1} zero z64 z64 z16 z16 z16
@uctest Unum{4,6} zero z64 z64 z16 z16 z16
@uctest Unum{4,8} zero z64 [z64, z64, z64, z64] z16 z16 z16

@uctest Unum{0,0} one z64 t64 z16 z16 z16
@uctest Unum{1,1} one o64 z64 z16 o16 z16
@uctest Unum{4,6} one o64 z64 z16 o16 z16
@uctest Unum{4,8} one o64 [z64, z64, z64, z64] z16 o16 z16
@uctest Unum{0,8} one z64 [t64, z64, z64, z64] z16 z16 z16

@uctest Unum{4,6} neg_one o64 z64 0x0002 o16 z16

@uctest Unum{0,0} nan o64 t64                                     0x0001 z16 z16
@uctest Unum{1,1} nan 0x0000_0000_0000_0003 0xC000_0000_0000_0000 0x0001 0x0001 0x0001
@uctest Unum{4,6} nan 0x0000_0000_0000_FFFF 0xFFFF_FFFF_FFFF_FFFF 0x0001 0x000F 0x003F
@uctest Unum{4,8} nan 0x0000_0000_0000_FFFF [f64, f64, f64, f64]  0x0001 0x000F 0x00FF

@uctest Unum{0,0} inf o64 t64                                     0x0000 z16 z16
@uctest Unum{1,1} inf 0x0000_0000_0000_0003 0xC000_0000_0000_0000 0x0000 0x0001 0x0001
@uctest Unum{4,6} inf 0x0000_0000_0000_FFFF 0xFFFF_FFFF_FFFF_FFFF 0x0000 0x000F 0x003F
@uctest Unum{4,8} inf 0x0000_0000_0000_FFFF [f64, f64, f64, f64]  0x0000 0x000F 0x00FF

@uctest Unum{0,0} neg_inf o64 t64                                     0x0002 z16 z16
@uctest Unum{1,1} neg_inf 0x0000_0000_0000_0003 0xC000_0000_0000_0000 0x0002 0x0001 0x0001
@uctest Unum{4,6} neg_inf 0x0000_0000_0000_FFFF 0xFFFF_FFFF_FFFF_FFFF 0x0002 0x000F 0x003F
@uctest Unum{4,8} neg_inf 0x0000_0000_0000_FFFF [f64, f64, f64, f64]  0x0002 0x000F 0x00FF

@uctest Unum{0,0} mmr 0x0000_0000_0000_0001 0x0000_0000_0000_0000                  0x0001    z16    z16
@uctest Unum{1,1} mmr 0x0000_0000_0000_0003 0x8000_0000_0000_0000                  0x0001 0x0001 0x0001
@uctest Unum{4,6} mmr 0x0000_0000_0000_FFFF 0xFFFF_FFFF_FFFF_FFFE                  0x0001 0x000F 0x003F
@uctest Unum{4,8} mmr 0x0000_0000_0000_FFFF [f64, f64, f64, 0xFFFF_FFFF_FFFF_FFFE] 0x0001 0x000F 0x00FF

@uctest Unum{0,0} neg_mmr 0x0000_0000_0000_0001 0x0000_0000_0000_0000                  0x0003    z16    z16
@uctest Unum{1,1} neg_mmr 0x0000_0000_0000_0003 0x8000_0000_0000_0000                  0x0003 0x0001 0x0001
@uctest Unum{4,6} neg_mmr 0x0000_0000_0000_FFFF 0xFFFF_FFFF_FFFF_FFFE                  0x0003 0x000F 0x003F
@uctest Unum{4,8} neg_mmr 0x0000_0000_0000_FFFF [f64, f64, f64, 0xFFFF_FFFF_FFFF_FFFE] 0x0003 0x000F 0x00FF

@uctest Unum{0,0} big_exact 0x0000_0000_0000_0001 0x0000_0000_0000_0000                  0x0000    z16    z16
@uctest Unum{1,1} big_exact 0x0000_0000_0000_0003 0x8000_0000_0000_0000                  0x0000 0x0001 0x0001
@uctest Unum{4,6} big_exact 0x0000_0000_0000_FFFF 0xFFFF_FFFF_FFFF_FFFE                  0x0000 0x000F 0x003F
@uctest Unum{4,8} big_exact 0x0000_0000_0000_FFFF [f64, f64, f64, 0xFFFF_FFFF_FFFF_FFFE] 0x0000 0x000F 0x00FF

@uctest Unum{0,0} neg_big_exact 0x0000_0000_0000_0001 0x0000_0000_0000_0000                  0x0002    z16    z16
@uctest Unum{1,1} neg_big_exact 0x0000_0000_0000_0003 0x8000_0000_0000_0000                  0x0002 0x0001 0x0001
@uctest Unum{4,6} neg_big_exact 0x0000_0000_0000_FFFF 0xFFFF_FFFF_FFFF_FFFE                  0x0002 0x000F 0x003F
@uctest Unum{4,8} neg_big_exact 0x0000_0000_0000_FFFF [f64, f64, f64, 0xFFFF_FFFF_FFFF_FFFE] 0x0002 0x000F 0x00FF

@uctest Unum{0,0} sss z64 z64                   0x0001    z16    z16
@uctest Unum{1,1} sss z64 z64                   0x0001 0x0001 0x0001
@uctest Unum{4,6} sss z64 z64                   0x0001 0x000F 0x003F
@uctest Unum{4,8} sss z64 [z64, z64, z64, z64]  0x0001 0x000F 0x00FF

@uctest Unum{0,0} neg_sss z64 z64                    0x0003    z16    z16
@uctest Unum{1,1} neg_sss z64 z64                    0x0003 0x0001 0x0001
@uctest Unum{4,6} neg_sss z64 z64                    0x0003 0x000F 0x003F
@uctest Unum{4,8} neg_sss z64 [z64, z64, z64, z64]   0x0003 0x000F 0x00FF

@uctest Unum{0,0} small_exact z64 t64                   0x0000    z16    z16
@uctest Unum{1,1} small_exact z64 0x4000_0000_0000_0000 0x0000 0x0001 0x0001
@uctest Unum{4,6} small_exact z64 o64                   0x0000 0x000F 0x003F
@uctest Unum{4,8} small_exact z64 [z64, z64, z64, o64]  0x0000 0x000F 0x00FF

@uctest Unum{0,0} neg_small_exact z64 t64                   0x0002    z16    z16
@uctest Unum{1,1} neg_small_exact z64 0x4000_0000_0000_0000 0x0002 0x0001 0x0001
@uctest Unum{4,6} neg_small_exact z64 o64                   0x0002 0x000F 0x003F
@uctest Unum{4,8} neg_small_exact z64 [z64, z64, z64, o64]  0x0002 0x000F 0x00FF

#this doesn't seem to always work.
@test is_mmr(mmr(Unum{0,0}))
@test is_mmr(mmr(Unum{1,1}))
@test is_mmr(mmr(Unum{2,2}))
@test is_mmr(mmr(Unum{3,3}))
@test is_mmr(mmr(Unum{3,5}))
@test is_mmr(mmr(Unum{4,6}))
@test is_mmr(mmr(Unum{4,7}))
#turns out that the generation of the top_mask implicitly converted an UInt16 to
#Int64 by adding one, which has a different top_mask definition.

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
using Underscores
using Test

@testset "Underscores Examples" begin
    @test [1,"a",2.0] == @_ map(_, [1,"a",2.0])

    strs = ["ab", "ca", "ad"]
    @test [true, false, true] == @_ map(startswith(_, "a"), strs)

    data = [(x="a", y=1),
            (x="b", y=2),
            (x="c", y=3)]

    @test data[1:1] == @_ filter(startswith(_.x, "a"), data)
    @test data[2:3] == @_ filter(_.y >= 2, data)

    # Use with indexing
    @test data[1] == @_ filter(startswith(_.x, "a"), data)[end]
    @test data[2:3] == @_ filter(_.y >= 2, data)[1:2]
    @test data[2] == @_ data[findfirst(_.y >= 2, data)]

    # Operators
    @test -2 == @_ -findfirst(_.x == "b", data)
    @test [14] == @_ [1 2 3] * map(_.y, data)
    @test (1:3)./3 == @_ data |> map(_.y, __) ./ length(__)

    # Multiple args
    @test [0,0] == @_ map(_-_, [1,2])

    # Ternary
    @test ["a","b","c"] == @_ data |>
        (any(_.y>2, __) ? map(_.x, __) : map(_.y, __))

    # Use with piping and __
    @test [1] == @_ data |>
                    filter(startswith(_.x, "a"), __) |>
                    map(_.y, __)

    @test [3,2,1] == @_ data |>
                        sort(__, by=_.x, rev=true) |>
                        map(_.y, __)

    @test [1] == (@_ map(_.y, __) ∘
                     filter(startswith(_.x, "a"), __))(data)

    @test [0,0,0] == @_ data |> map(_1.y + _2, __, [-1,-2,-3])

    @test 3 == @_ [[1], [1,2], [1,2,3]] |>
                                map(_[_[end]], __) |>
                                __[end]

    # Use with piping and lazy versions of map and filter
    Filter(f) = x->filter(f,x)
    Map(f) = x->map(f,x)

    @test [1] == @_ data |>
                    Filter(startswith(_.x, "a")) |>
                    Map(_.y)

    @test [1] == @_(Map(_.y) ∘ Filter(startswith(_.x, "a")))(data)

    # Getproperty
    @test "a" == @_ data |> __[1].x
    @test "b" == @_ (x="a", y=(z="b", t="c")) |> __.y.z
    @test "c" == @_ (x="a", y=(z="b", t="c")) |> __.y[end]
    @test -3 == @_ [1+2im,3,4+5im] |> sum(_^2, __).re

    # Interpolation
    @test ["1","a","2.0"] == @_ map("$_", [1,"a",2.0])

    # Broadcasting
    @test [[1],[2],[3]] == @_ data |> filter.(_ isa Int, collect.(__))
    @test [9,25] == @_ 2 |>  __ .+ (1:3)  .|>  __^2 |> filter(isodd(_), __)
    @test sqrt.([3,2,1]) == @_ data .|> identity |>
                        sort(__, by=_.x, rev=true) |>
                        map(_.y, __) .|> sqrt

    # Comprehensions
    @test [[],[],[3]] == @_ [[1], [1,2], [1,2,3]] |>
                                [filter(_>2, x) for x in __]
    @test fill(:y,3) == @_ data |>
                        Symbol[findfirst(_ isa Int, x) for x in __]

    # Matrix construction
    @test 4*ones(2,2) == @_ ones(4) |> [ sum(__)      sum(_^2, __)
                                         sum(_^3, __) sum(_^4, __) ]
end

@testset "Underscores lowering" begin
    cleanup! = Base.remove_linenums!
    lower(ex) = cleanup!(Underscores.lower_underscores(ex))

    # Simple cases
    # _
    @test lower(:(f(_))) == cleanup!(:(f((_1,)->_1)))
    @test lower(:(f(g(h(_))))) == cleanup!(:(f(((_1,)->g(h(_1))))))
    # __
    @test lower(:(f(__))) == cleanup!(:((__1,)->f(__1)))
    @test lower(:(f(g(h(__))))) == cleanup!(:((__1,)->f(g(h(__1)))))
    # Tricky case when function name is itself an underscore expression.
    # We want the function name to be also expanded here.
    @test lower(:(g(_)(x))) == cleanup!(:(((_1,)->g(_1))(x)))

    # Multiple arguments
    # _
    @test lower(:(f(_,a))) == cleanup!(:(f(((_1,)->_1), a)))
    @test lower(:(f(a,_))) == cleanup!(:(f(a, ((_1,)->_1))))
    @test lower(:(f(_,_))) == cleanup!(:(f(((_1,)->_1), ((_1,)->_1))))
    @test lower(:(f(g(_,_)))) == cleanup!(:(f(((_1,)->g(_1,_1)))))
    # __
    @test lower(:(f(__,__))) == cleanup!(:((__1,)->f(__1, __1)))

    # Numbered arguments
    # _
    @test lower(:(f(_1))) == cleanup!(:(f((_1,)->_1)))
    @test lower(:(f(_2))) == cleanup!(:(f((_1,_2)->_2)))
    @test lower(:(f(_2+_1))) == cleanup!(:(f((_1,_2)->_2+_1)))
    # __
    @test lower(:(f(__1))) == cleanup!(:((__1,)->f(__1)))
    @test lower(:(f(__2))) == cleanup!(:((__1,__2)->f(__2)))
    @test lower(:(f(__2+__1))) == cleanup!(:((__1,__2)->f(__2+__1)))

    # Cute subscript-numbered arguments
    @test lower(:(f(_₁))) == cleanup!(:(f((_1,)->_1)))
    @test lower(:(f(_₂))) == cleanup!(:(f((_1,_2)->_2)))
    @test lower(:(f(__₁))) == cleanup!(:((__1,)->f(__1)))

    # Can't mix numbered and non-numbered placeholders
    @test_throws ArgumentError lower(:(f(_+_1)))
    @test_throws ArgumentError lower(:(f(__+__1)))

    # piping and composition
    # _
    @test lower(:(f(_) |> g(_) |> h(_))) ==
          cleanup!(:(f((_1,)->_1) |>
                     g((_1,)->_1) |>
                     h((_1,)->_1)))
    @test lower(:(f(_) ∘ g(_) ∘ h(_))) ==
          cleanup!(:(f((_1,)->_1) ∘
                     g((_1,)->_1) ∘
                     h((_1,)->_1)))
    @test lower(:(f(_) <| g(_) <| h(_))) ==
          cleanup!(:(f((_1,)->_1) <|
                     g((_1,)->_1) <|
                     h((_1,)->_1)))
    @test lower(:(f(_) .|> g(_) .|> h(_))) ==
          cleanup!(:(f((_1,)->_1) .|>
                     g((_1,)->_1) .|>
                     h((_1,)->_1)))
    # __
    @test lower(:(f(_, __) |> g(_, __))) ==
        cleanup!(:(((__1,)->f((_1,)->_1, __1)) |>
                   ((__1,)->g((_1,)->_1, __1))))
    @test lower(:(f(_, __) .|> g(__))) ==
        cleanup!(:(((__1,)->f((_1,)->_1, __1)) .|> ((__1,)->g(__1))))

    # Keyword arguments
    @test lower(:(f(x, k=_+1))) == cleanup!(:(f(x, k=((_1,)->_1+1))))
    @test lower(:(f(x; k=_+1))) == cleanup!(:(f(x; k=((_1,)->_1+1))))

    # Broadcast notation
    @test lower(:(f.(_)))   == cleanup!(:(f.((_1,)->_1)))
    @test lower(:(f.(_,_))) == cleanup!(:(f.((_1,)->_1, (_1,)->_1)))
    @test lower(:(f.(__)))  == cleanup!(:((__1,)->f.(__1)))
    @test lower(:((_).(x)))  == cleanup!(:(((_1,)->_1).(x)))

    # Indexing
    @test lower(:(f(_)[2])) == cleanup!(:(f((_1,)->_1)[2]))
    @test lower(:(f(g(_[3]),h)[4])) == cleanup!(:(f((_1,)->g(_1[3]),h)[4]))
    @test lower(:(f(__)[5])) == cleanup!(:((__1,)->f(__1)[5]))

    # Quote
    @test lower(:(map(:(a[$_]), 1:4))) == cleanup!(:(map((_1,)->:(a[$_1]), 1:4)))

    # Random sample of other syntax
    @test lower(:([_]))  == cleanup!(:([(_1,)->_1]))
    @test lower(:((f(_), g(_))))  == cleanup!(:(((_1,)->f(_1)), ((_1,)->g(_1))))
    @test lower(:([__])) == cleanup!(:((__1,)->[__1]))
    @test lower(:(__.x)) == cleanup!(:((__1,)->__1.x))

    # do syntax is disabled for now because desired behaviour is not entirely
    # clear. See #4
    @test_throws ErrorException lower(:(f() do
                                            body
                                        end))
end


text:
"""
    connect_mesh(EToV,fv)

Initialize element connectivity matrices, element to element and element to face
connectivity.

Inputs:
- `EToV` is a `K` by `Nv` matrix whose rows identify the `Nv` vertices
which make up an element.
- `fv` (an array of arrays containing unordered indices of face vertices).

Output: `FToF`, `length(fv)` by `K` index array containing face-to-face connectivity.
"""
function connect_mesh(EToV, fv)
    Nfaces = length(fv)
    K = size(EToV,1)

    # sort and find matches
    fnodes = [[sort(EToV[e,ids]) for ids in fv, e = 1:K]...]
    p = sortperm(fnodes) # sorts by lexicographic ordering by default
    fnodes = fnodes[p,:]

    FToF = reshape(collect(1:Nfaces*K),Nfaces,K)
    for f = 1:size(fnodes,1)-1
        if fnodes[f,:]==fnodes[f+1,:]
            f1 = FToF[p[f]]
            f2 = FToF[p[f+1]]
            FToF[p[f]] = f2
            FToF[p[f+1]] = f1
        end
    end
    return FToF
end

"""
    build_node_maps(Xf, FToF)

Intialize the connectivity table along all edges and boundary node tables of all
elements. `mapM` - map minus (interior). `mapP` - map plus (exterior).

`Xf = (xf, yf, zf)` and `FToF` is size `(Nfaces * K)` and `FToF[face]` = face neighbor

`mapM`, `mapP` are size `Nfp` x `(Nfaces*K)`

# Examples
```julia
julia> mapM, mapP, mapB = build_node_maps((xf, yf), FToF)
```
"""

function build_node_maps(FToF, Xf...; tol = 1e-12)

    NODETOL = tol
    NfacesK = length(FToF)
    dims = length(Xf)

    # number nodes consecutively
    Nfp  = length(Xf[1]) ÷ NfacesK
    mapM = reshape(collect(1:length(Xf[1])), Nfp, NfacesK);
    mapP = copy(mapM);

    D = zeros(Nfp,Nfp)
    idM,idP = zeros(Int,Nfp),zeros(Int,Nfp)
    for (f1,f2) in enumerate(FToF)

        fill!(D,zero(eltype(D)))

        # find volume node numbers of left and right nodes
        for i = 1:dims
            Xfi = reshape(Xf[i],Nfp,NfacesK)
            for j = 1:Nfp, k = 1:Nfp
                D[j,k] += abs(Xfi[j,f1]-Xfi[k,f2])
            end
        end

        refd = maximum(D[:])
        map!(id->id[1], idM, findall(@. D < NODETOL*refd))
        map!(id->id[2], idP, findall(@. D < NODETOL*refd))
        @. mapP[idM,f1] = idP + (f2-1)*Nfp
    end

    mapB = map(x->x[1],findall(@. mapM[:]==mapP[:]))
    return mapM, mapP, mapB
end

# old deprecated interface
@deprecate make_periodic(rd::RefElemData, md::MeshData) make_periodic(md)
@deprecate make_periodic(rd::RefElemData, md::MeshData, args...) make_periodic(md,args...)
@deprecate make_periodic(md::MeshData, rd::RefElemData, args...) make_periodic(md,args...)

"""
    make_periodic(md::MeshData{Dim}, is_periodic...) where {Dim}
    make_periodic(md::MeshData{Dim}, is_periodic = ntuple(x->true,Dim)) where {Dim}
    make_periodic(md::MeshData, is_periodic = true)

Returns new MeshData such that the node maps `mapP` and face maps `FToF` are now periodic.
Here, `is_periodic` is a tuple of `Bool` indicating whether or not to impose periodic
BCs in the `x`,`y`, or `z` coordinate.
"""
make_periodic(md::MeshData{Dim}, is_periodic::Bool = true) where {Dim} = make_periodic(md, ntuple(_->is_periodic, Dim))

function make_periodic(md::MeshData{Dim}, is_periodic::NTuple{Dim, Bool}) where {Dim, Bool}

    @unpack mapM, mapP, mapB, xyzf, FToF = md
    NfacesTotal = length(FToF)
    FToF_periodic = copy(FToF)
    mapPB = build_periodic_boundary_maps!(xyzf...,is_periodic...,NfacesTotal,
                                          mapM, mapP, mapB, FToF_periodic)
    mapP_periodic = copy(mapP)
    mapP_periodic[mapB] = mapPB
    mapB_periodic = mapB[mapPB .== mapP[mapB]] # keep only non-periodic boundary nodes
    return setproperties(md,(mapB=mapB_periodic, mapP = mapP_periodic,
                             FToF = FToF_periodic, is_periodic = is_periodic)) # from Setfield.jl
end

# specializes to 1D - periodic = find min/max indices of xf and reverse their order
function make_periodic(md::MeshData{1, Tv, Ti}, is_periodic::Bool = true) where {Tv, Ti}

    if is_periodic == true
        @unpack mapP, mapB, xf, FToF = md
        mapPB = argmax(vec(xf)),argmin(vec(xf))
        mapP_periodic = copy(mapP)
        mapP_periodic[mapB] .= mapPB
        FToF_periodic = copy(FToF)
        FToF_periodic[[1,length(FToF)]] .= mapPB
        mapB_periodic = Ti[]
        return setproperties(md,(mapB = mapB_periodic, mapP = mapP_periodic,
                                 FToF = FToF_periodic, is_periodic = (true,))) # from Setfield.jl
    end
end

# Helper functions for `make_nodemaps_periodic!`, 2D version which modifies FToF.
function build_periodic_boundary_maps!(xf, yf, is_periodic_x, is_periodic_y,
                                       NfacesTotal, mapM, mapP, mapB, FToF)

        # find boundary faces (e.g., when FToF[f] = f)
        Flist = 1:length(FToF)
        Bfaces = findall(vec(FToF) .== Flist)

        xb = xf[mapB]
        yb = yf[mapB]
        Nfp = convert(Int,length(xf)/NfacesTotal)
        Nbfaces = convert(Int,length(xb)/Nfp)
        xb = reshape(xb,Nfp,Nbfaces)
        yb = reshape(yb,Nfp,Nbfaces)

        # compute centroids of faces
        xc = vec(sum(xb,dims=1)/Nfp)
        yc = vec(sum(yb,dims=1)/Nfp)
        mapMB = reshape(mapM[mapB],Nfp,Nbfaces)
        mapPB = reshape(mapP[mapB],Nfp,Nbfaces)

        xmin,xmax = extrema(xc)
        ymin,ymax = extrema(yc)

        NODETOL = 1e-12
        LX,LY = map((x->x[2]-x[1])∘extrema,(xf,yf))
        if abs(abs(xmax-xmin)-LX)>NODETOL && is_periodic_x
            error("periodicity requested in x, but max_dist(xf) != max_dist(xc)")
        end
        if abs(abs(ymax-ymin)-LY)>NODETOL && is_periodic_y
            error("periodicity requested in y, but max_dist(yf) != max_dist(yc)")
        end

        # determine which faces lie on x and y boundaries
        Xscale, Yscale = max(1,LX), max(1,LY)
        yfaces = map(x->x[1],findall(@. (@. abs(yc-ymax)<NODETOL*Yscale) | (@. abs(yc-ymin)<NODETOL*Yscale)))
        xfaces = map(x->x[1],findall(@. (@. abs(xc-xmax)<NODETOL*Xscale) | (@. abs(xc-xmin)<NODETOL*Xscale)))

        D = zeros(eltype(xb), size(xb,1), size(xb,1))
        ids = zeros(Int,size(xb,1))
        if is_periodic_y # find matches in y faces
            for i in yfaces, j in yfaces
                if i!=j
                    if abs(xc[i]-xc[j])<NODETOL*Xscale && abs(abs(yc[i]-yc[j])-LY)<NODETOL*Yscale
                        # create distance matrix
                        @. D = abs(xb[:,i] - xb[:,j]')
                        map!(x->x[1], ids, findall(@. D < NODETOL*Xscale))
                        @. mapPB[:,i] = mapMB[ids,j]

                        FToF[Bfaces[i]] = Bfaces[j]
                    end
                end
            end
        end

        if is_periodic_x # find matches in x faces
            for i in xfaces, j in xfaces
                if i!=j
                    if abs(yc[i]-yc[j])<NODETOL*Yscale && abs(abs(xc[i]-xc[j])-LX)<NODETOL*Xscale
                        @. D = abs(yb[:,i] - yb[:,j]')
                        map!(x->x[1], ids, findall(@. D < NODETOL*Yscale))
                        mapPB[:,i]=mapMB[ids,j]

                        FToF[Bfaces[i]] = Bfaces[j]
                    end
                end
            end
        end

        return mapPB[:]
end

# 3D version of build_periodic_boundary_maps, modifies FToF
function build_periodic_boundary_maps!(xf, yf, zf,
                                       is_periodic_x, is_periodic_y, is_periodic_z,
                                       NfacesTotal, mapM, mapP, mapB, FToF)

    # find boundary faces (e.g., when FToF[f] = f)
    Flist = 1:length(FToF)
    Bfaces = findall(vec(FToF) .== Flist)

    xb,yb,zb = xf[mapB],yf[mapB],zf[mapB]
    Nfp = length(xf) ÷ NfacesTotal
    Nbfaces = length(xb) ÷ Nfp
    xb, yb, zb = (x->reshape(x, Nfp, Nbfaces)).((xb, yb, zb))

    # compute centroids of faces
    xc = vec(sum(xb,dims=1)/Nfp)
    yc = vec(sum(yb,dims=1)/Nfp)
    zc = vec(sum(zb,dims=1)/Nfp)
    mapMB = reshape(mapM[mapB],Nfp,Nbfaces)
    mapPB = reshape(mapP[mapB],Nfp,Nbfaces)

    xmin,xmax = extrema(xc)
    ymin,ymax = extrema(yc)
    zmin,zmax = extrema(zc)

    NODETOL = 1e-12
    LX,LY,LZ = map((x->x[2]-x[1])∘extrema,(xf,yf,zf))
    if abs(abs(xmax-xmin)-LX)>NODETOL && is_periodic_x
        error("periodicity requested in x, but max_dist(xf) != max_dist(xc)")
    end
    if abs(abs(ymax-ymin)-LY)>NODETOL && is_periodic_y
        error("periodicity requested in y, but max_dist(yf) != max_dist(yc)")
    end
    if abs(abs(zmax-zmin)-LZ)>NODETOL && is_periodic_z
        error("periodicity requested in z, but max_dist(zf) != max_dist(zc)")
    end

    # determine which faces lie on x and y boundaries
    xfaces = map(x->x[1],findall(@. (@. abs(xc-xmax)<NODETOL*LX) | (@. abs(xc-xmin)<NODETOL*LX)))
    yfaces = map(x->x[1],findall(@. (@. abs(yc-ymax)<NODETOL*LY) | (@. abs(yc-ymin)<NODETOL*LY)))
    zfaces = map(x->x[1],findall(@. (@. abs(zc-zmax)<NODETOL*LZ) | (@. abs(zc-zmin)<NODETOL*LZ)))

    D = zeros(eltype(xb), size(xb,1), size(xb,1))
    ids = zeros(Int,size(xb,1))
    if is_periodic_x # find matches in x faces
        for i in xfaces, j in xfaces
            if i!=j
                if abs(yc[i]-yc[j]) < NODETOL*LY && abs(zc[i]-zc[j]) < NODETOL*LZ && abs(abs(xc[i]-xc[j])-LX) < NODETOL*LX
                    # create distance matrix
                    @. D = abs(yb[:,i] - yb[:,j]') + abs(zb[:,i] - zb[:,j]')
                    map!(x->x[1], ids, findall(@.D < NODETOL*LY))
                    @. mapPB[:,i] = mapMB[ids,j]

                    FToF[Bfaces[i]] = Bfaces[j]
                end
            end
        end
    end

    # find matches in y faces
    if is_periodic_y
        for i in yfaces, j = yfaces
            if i!=j
                if abs(xc[i]-xc[j]) < NODETOL*LX && abs(zc[i]-zc[j]) < NODETOL*LZ && abs(abs(yc[i]-yc[j])-LY) < NODETOL*LY
                    @. D = abs(xb[:,i] - xb[:,j]') + abs(zb[:,i] - zb[:,j]')
                    map!(x->x[1], ids, findall(@.D < NODETOL*LX))
                    @. mapPB[:,i] = mapMB[ids,j]

                    FToF[Bfaces[i]] = Bfaces[j]
                end
            end
        end
    end

    # find matches in z faces
    if is_periodic_z
        for i in zfaces, j in zfaces
            if i!=j
                if abs(xc[i]-xc[j]) < NODETOL*LX && abs(yc[i]-yc[j]) < NODETOL*LY && abs(abs(zc[i]-zc[j])-LZ) < NODETOL*LZ
                    @. D = abs(xb[:,i] - xb[:,j]') + abs(yb[:,i] - yb[:,j]')
                    map!(x->x[1], ids, findall(@.D < NODETOL*LX))
                    @. mapPB[:,i] = mapMB[ids,j]

                    FToF[Bfaces[i]] = Bfaces[j]
                end
            end
        end
    end

    return mapPB[:]
end


text:
# Julia wrapper for header: /usr/include/libavutil/avutil.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    avutil_version,
    av_version_info,
    avutil_configuration,
    avutil_license,
    av_get_media_type_string,
    av_get_picture_type_char,
    av_x_if_null,
    av_int_list_length_for_size,
    av_fopen_utf8,
    av_get_time_base_q,
    av_fourcc_make_string

function avutil_version()
    ccall((:avutil_version, libavutil), UInt32, ())
end

function av_version_info()
    ccall((:av_version_info, libavutil), Cstring, ())
end

function avutil_configuration()
    ccall((:avutil_configuration, libavutil), Cstring, ())
end

function avutil_license()
    ccall((:avutil_license, libavutil), Cstring, ())
end

function av_get_media_type_string(media_type::AVMediaType)
    ccall((:av_get_media_type_string, libavutil), Cstring, (AVMediaType,), media_type)
end

function av_get_picture_type_char(pict_type::AVPictureType)
    ccall((:av_get_picture_type_char, libavutil), UInt8, (AVPictureType,), pict_type)
end

function av_x_if_null(p, x)
    ccall((:av_x_if_null, libavutil), Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}), p, x)
end

function av_int_list_length_for_size(elsize::Integer, list, term::UInt64)
    ccall((:av_int_list_length_for_size, libavutil), UInt32, (UInt32, Ptr{Cvoid}, UInt64), elsize, list, term)
end

function av_fopen_utf8(path, mode)
    ccall((:av_fopen_utf8, libavutil), Ptr{FILE}, (Cstring, Cstring), path, mode)
end

function av_get_time_base_q()
    ccall((:av_get_time_base_q, libavutil), AVRational, ())
end

function av_fourcc_make_string(buf, fourcc::Integer)
    ccall((:av_fourcc_make_string, libavutil), Cstring, (Cstring, UInt32), buf, fourcc)
end


text:
#agent declaration
approximator_GNN = GeometricFlux.GraphConv(64 => 64, Flux.leakyrelu)
target_approximator_GNN = GeometricFlux.GraphConv(64 => 64, Flux.leakyrelu)
gnnlayers = 10
numInFeatures = 3
approximator_model = SeaPearl.CPNN(
    graphChain = Flux.Chain(
        GeometricFlux.GraphConv(numInFeatures => 64, Flux.leakyrelu),
        [approximator_GNN for i = 1:gnnlayers]...
    ),
    nodeChain = Flux.Chain(
        Flux.Dense(64, 32, Flux.leakyrelu),
        Flux.Dense(32, 32, Flux.leakyrelu),
        Flux.Dense(32, 16, Flux.leakyrelu),
    ),
    outputChain = Flux.Dense(16, 3),
)
target_approximator_model = SeaPearl.CPNN(
    graphChain = Flux.Chain(
        GeometricFlux.GraphConv(numInFeatures => 64, Flux.leakyrelu),
        [target_approximator_GNN for i = 1:gnnlayers]...
    ),
    nodeChain = Flux.Chain(
        Flux.Dense(64, 32, Flux.leakyrelu),
        Flux.Dense(32, 32, Flux.leakyrelu),
        Flux.Dense(32, 16, Flux.leakyrelu),
    ),
    outputChain = Flux.Dense(16, 3),
)

agent = RL.Agent(
    policy = RL.QBasedPolicy(
        learner = RL.DQNLearner(
            approximator = RL.NeuralNetworkApproximator(
                model = approximator_model,
                optimizer = ADAM(0.0005f0)
            ),
            target_approximator = RL.NeuralNetworkApproximator(
                model = target_approximator_model,
                optimizer = ADAM(0.0005f0)
            ),
            loss_func = Flux.Losses.huber_loss,
            stack_size = nothing,
            γ = 0.9999f0,
            batch_size = 1, #32,
            update_horizon = 25,
            min_replay_history = 1,
            update_freq = 10,
            target_update_freq = 200,
        ),
        explorer = RL.EpsilonGreedyExplorer(
            ϵ_stable = 0.001,
            kind = :exp,
            ϵ_init = 1.0,
            warmup_steps = 0,
            decay_steps = 5000,
            step = 1,
            is_break_tie = false,
            #is_training = true,
            rng = MersenneTwister(33)
        )
    ),
    trajectory = RL.CircularArraySARTTrajectory(
        capacity = 10,
        state = SeaPearl.DefaultTrajectoryState[] => (),
    )
)

@testset "basicmetrics.jl" begin

    @testset "constructor basicmetrics.jl" begin
        trailer = SeaPearl.Trailer()
        cpmodel = SeaPearl.CPModel(trailer)
        basicheuristic = SeaPearl.BasicHeuristic()

        metrics  = SeaPearl.BasicMetrics(cpmodel, basicheuristic)

        @test isempty(metrics.nodeVisited)== true
        @test isempty(metrics.meanNodeVisitedUntilEnd) == true
        @test isempty(metrics.timeneeded) == true
        @test isnothing(metrics.scores) == true
        @test isnothing(metrics.totalReward) == true
        @test isnothing(metrics.loss) == true
        @test metrics.nbEpisodes == 0

        metrics  = SeaPearl.BasicMetrics(cpmodel, basicheuristic;meanOver=100)
        @test metrics.meanOver == 100

        trailer = SeaPearl.Trailer()
        cpmodel = SeaPearl.CPModel(trailer)
        x = SeaPearl.IntVar(1, 2, "x", trailer)
        y = SeaPearl.IntVar(1, 2, "y", trailer)

        SeaPearl.addVariable!(cpmodel, x)
        SeaPearl.addVariable!(cpmodel, y)

        basicheuristic = SeaPearl.BasicHeuristic()
        learnedheuristic = SeaPearl.LearnedHeuristic(agent)

        metrics  = SeaPearl.BasicMetrics(cpmodel, basicheuristic)
        @test typeof(metrics) ==  SeaPearl.BasicMetrics{SeaPearl.DontTakeObjective, SeaPearl.BasicHeuristic}

        metrics  = SeaPearl.BasicMetrics(cpmodel, learnedheuristic)
        @test typeof(metrics) ==  SeaPearl.BasicMetrics{SeaPearl.DontTakeObjective, SeaPearl.LearnedHeuristic{SeaPearl.DefaultStateRepresentation{SeaPearl.DefaultFeaturization,SeaPearl.DefaultTrajectoryState}, SeaPearl.DefaultReward, SeaPearl.FixedOutput}}
        @test isnothing(metrics.totalReward) == false
        @test isnothing(metrics.loss) == false

        SeaPearl.addObjective!(cpmodel,y)

        metrics  = SeaPearl.BasicMetrics(cpmodel, basicheuristic)
        @test typeof(metrics) ==  SeaPearl.BasicMetrics{SeaPearl.TakeObjective, SeaPearl.BasicHeuristic}
        @test isnothing(metrics.scores) == false

        metrics  = SeaPearl.BasicMetrics(cpmodel, learnedheuristic)
        @test typeof(metrics) ==  SeaPearl.BasicMetrics{SeaPearl.TakeObjective, SeaPearl.LearnedHeuristic{SeaPearl.DefaultStateRepresentation{SeaPearl.DefaultFeaturization,SeaPearl.DefaultTrajectoryState}, SeaPearl.DefaultReward, SeaPearl.FixedOutput}}

    end

    @testset "BasicMetrics{DontTakeObjective, BasicHeuristic} " begin

        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)
        basicheuristic = SeaPearl.BasicHeuristic()

        x = SeaPearl.IntVar(2, 3, "x", trailer)
        y = SeaPearl.IntVar(2, 3, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))

        metrics  = SeaPearl.BasicMetrics(model, basicheuristic)
        dt = @elapsed SeaPearl.search!(model, SeaPearl.DFSearch(), SeaPearl.MinDomainVariableSelection(), basicheuristic)
        metrics(model,dt)

        @test metrics.nodeVisited[1] == [2,3]
        @test size(metrics.timeneeded,1) == 1
        @test metrics.nbEpisodes == 1

    end

    @testset "BasicMetrics{DontTakeObjective, LearnedHeuristic}" begin

        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)
        learnedheuristic = SeaPearl.LearnedHeuristic(agent)
        x = SeaPearl.IntVar(2, 3, "x", trailer)
        y = SeaPearl.IntVar(2, 3, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))

        metrics  = SeaPearl.BasicMetrics(model, learnedheuristic)
        dt = @elapsed SeaPearl.search!(model, SeaPearl.DFSearch(), SeaPearl.MinDomainVariableSelection(), learnedheuristic)

        metrics(model,dt)

        @test isempty(metrics.nodeVisited) == false #the number of solution found and their relative score depend on the stochastic heuristic...
        @test metrics.nbEpisodes == 1
        @test size(metrics.timeneeded,1) == 1
        @test size(metrics.totalReward,1) == 1
        @test size(metrics.loss,1) == 1

    end

    @testset "BasicMetrics{TakeObjective, BasicHeuristic}" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)
        basicheuristic = SeaPearl.BasicHeuristic()

        x = SeaPearl.IntVar(2, 3, "x", trailer)
        y = SeaPearl.IntVar(2, 3, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))
        SeaPearl.addObjective!(model,y)
        metrics  = SeaPearl.BasicMetrics(model, basicheuristic)
        dt = @elapsed SeaPearl.search!(model, SeaPearl.DFSearch(), SeaPearl.MinDomainVariableSelection(), basicheuristic)
        metrics(model,dt)

        @test metrics.nodeVisited[1] == [2, 3] #only one soluion found due to Objective prunning
        @test metrics.scores[1] == [3, 2]
        @test size(metrics.timeneeded,1) == 1
        @test metrics.nbEpisodes == 1
    end

    @testset "BasicMetrics{TakeObjective, LearnedHeuristic}" begin

        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)
        learnedheuristic = SeaPearl.LearnedHeuristic(agent)
        x = SeaPearl.IntVar(2, 3, "x", trailer)
        y = SeaPearl.IntVar(2, 3, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))
        SeaPearl.addObjective!(model,y)

        metrics  = SeaPearl.BasicMetrics(model, learnedheuristic)
        dt = @elapsed SeaPearl.search!(model, SeaPearl.DFSearch(), SeaPearl.MinDomainVariableSelection(), learnedheuristic)

        metrics(model,dt)

        @test isempty(metrics.nodeVisited) == false #the number of solution found and their relative score depend on the stochastic heuristic...
        @test isempty(metrics.scores) == false
        @test metrics.nbEpisodes == 1
        @test size(metrics.timeneeded,1) == 1
        @test size(metrics.totalReward,1) == 1
        @test size(metrics.loss,1) == 1
    end

    @testset "test !" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)
        basicheuristic = SeaPearl.BasicHeuristic()
        metrics  = SeaPearl.BasicMetrics(model, basicheuristic; meanOver=4)

        push!(metrics.nodeVisited, [10,0])
        push!(metrics.nodeVisited, [10,0])
        push!(metrics.nodeVisited, [10,0])
        push!(metrics.nodeVisited, [10,0])
        push!(metrics.nodeVisited, [0,10])
        push!(metrics.nodeVisited, [0,10])
        push!(metrics.nodeVisited, [0,10])
        push!(metrics.nodeVisited, [0,10])

        push!(metrics.meanNodeVisitedUntilEnd, 0)
        push!(metrics.meanNodeVisitedUntilEnd, 0)
        push!(metrics.meanNodeVisitedUntilEnd, 0)
        push!(metrics.meanNodeVisitedUntilEnd, 0)
        push!(metrics.meanNodeVisitedUntilEnd, 10)
        push!(metrics.meanNodeVisitedUntilEnd, 10)
        push!(metrics.meanNodeVisitedUntilEnd, 10)
        push!(metrics.meanNodeVisitedUntilEnd, 10)

        push!(metrics.meanNodeVisitedUntilfirstSolFound, 10)
        push!(metrics.meanNodeVisitedUntilfirstSolFound, 10)
        push!(metrics.meanNodeVisitedUntilfirstSolFound, 10)
        push!(metrics.meanNodeVisitedUntilfirstSolFound, 10)
        push!(metrics.meanNodeVisitedUntilfirstSolFound, 0)
        push!(metrics.meanNodeVisitedUntilfirstSolFound, 0)
        push!(metrics.meanNodeVisitedUntilfirstSolFound, 0)
        push!(metrics.meanNodeVisitedUntilfirstSolFound, 0)

        metrics.nbEpisodes = 8
        SeaPearl.computemean!(metrics)
        @test metrics.meanNodeVisitedUntilfirstSolFound == [10.0, 7.5, 5.0, 2.5, 0.0]
        @test metrics.meanNodeVisitedUntilEnd == [0.0, 2.5, 5.0, 7.5, 10.0]

    end
end


text:
using MathPhysicalConstants, Measurements, Unitful

using Test

using MathPhysicalConstants

@testset begin
    @test MathPhysicalConstants.MKS.PlancksConstantH == 6.62606896e-34
    @test iszero(measurement(MathPhysicalConstants.IS.SanchezElectrConstant) - measurement(MathPhysicalConstants.IS.SanchezElectrConstant))
    # \\\ @test setprecision(BigFloat, 768) do; precision(ustrip(big(c))) end == 768
    #\\\@test measurement(h) === measurement(h)
    #\\\@test iszero(measurement(α) - measurement(α))
    #\\\@test isone(measurement(BigFloat, atm) / measurement(BigFloat, atm))
    #\\\@test iszero(measurement(BigFloat, ħ) - (measurement(BigFloat, h) / 2big(pi)))
    #\\\@test isone(measurement(BigFloat, ħ) / (measurement(BigFloat, h) / 2big(pi)))
end


text:
using SimilarityMetrics

# Set

a, b = Set{Int}(), Set{Int}()

add!(a, 1); add!(a, 2)
add!(b, 2); add!(b, 3)
expl2(a, b)

add!(a, 4)
expl2(a, b)

add!(b, 1)
expl2(a, b)

# IntSet

a, b = IntSet(), IntSet()

add!(a, 1); add!(a, 2)
add!(b, 2); add!(b, 3)
expl2(a, b)

add!(a, 4)
expl2(a, b)

add!(b, 1)
expl2(a, b)

# AbstractVector
a, b = falses(4), falses(4)

a[1] = true; a[2] = true
b[2] = true; b[3] = true
expl2(a, b)

a[4] = true
expl2(a, b)

b[1] = true
expl2(a, b)


text:
module SimplexTypes
    using Requires
    function __init__()
        @require Simplices="d5428e67-3037-59ba-9ab1-57a04f0a3b6a" begin
            import .Simplices

            include("AbstractSimplex.jl")
            include("ImmutableSimplex.jl")
            include("MutableSimplex.jl")
            include("simplex_subsampling.jl")
        end
    end
end


text:
using Test

using Pandoc
const Markdown = Pandoc.Markdown
using Markdown

BASE_DIR = abspath(joinpath(dirname(Base.find_package("Pandoc")), ".."))

# Copied @includetests from https://github.com/ssfrr/TestSetExtensions.jl.
# Ideally, we could import and use TestSetExtensions.  Its functionality was broken by changes
# in Julia v0.7.  Refer to https://github.com/ssfrr/TestSetExtensions.jl/pull/7.

"""
Includes the given test files, given as a list without their ".jl" extensions.
If none are given it will scan the directory of the calling file and include all
the julia files.
"""
macro includetests(testarg...)
    if length(testarg) == 0
        tests = []
    elseif length(testarg) == 1
        tests = testarg[1]
    else
        error("@includetests takes zero or one argument")
    end

    quote
        tests = $tests
        rootfile = @__FILE__
        if length(tests) == 0
            tests = readdir(dirname(rootfile))
            tests = filter(f->endswith(f, ".jl") && f != basename(rootfile), tests)
        else
            tests = map(f->string(f, ".jl"), tests)
        end
        println()
        for test in tests
            print(splitext(test)[1], ": ")
            include(test)
            println()
        end
    end
end

function run_tests()
    # Testing Topological components of the schema
    @time @testset "Begin tests" begin
        @includetests ARGS
    end
end

run_tests()


text:
export Grid2D

"""
    Grid2D(;xlims=(-1,1), ylim=(-1, 1), res=(21, 21))

Create a 2D grid.
"""
function Grid2D(;xlims=(-1, 1), ylims=(-1, 1), res::NTuple{2, Int}=(21, 21))
    Grid2D(
        central_range(xlims..., res[1]),
        central_range(ylims..., res[2]))
end
function Grid2D(xmin, xmax, ymin, ymax, xres, yres)
    Grid2D(
        central_range(xmin, xmax, xres),
        central_range(ymin, ymax, yres))
end
# Grid2D(;lims=(-1, 1), res::Int=21) = Grid2D(lims, res)
Grid2D(lims::Tuple{<:Real, <:Real}, res::Int) = Grid2D( central_range(lims..., res), central_range(lims..., res))
function Grid2D(f::MP.AbstractPolynomial, res=600; factor=1.5)
    Grid2D(amoeba_carcase_domain_heuristic(f, factor=factor)..., res, res)
end

Base.broadcastable(G::Grid2D) = Ref(G)
Base.size(G::Grid2D) = (length(G.xrange), length(G.yrange))
Base.size(G::Grid2D, i) = size(G)[i]
Base.length(G::Grid2D) = prod(size(G))
Base.ndims(G::Grid2D) = 3
@inline function Base.getindex(G::Grid2D, k::Integer)
    i, j = _ind2sub(size(G), k)
    G[i, j]
end
Base.getindex(G::Grid2D, I::CartesianIndex{2}) = G[Tuple(I)...]
@inline function Base.getindex(G::Grid2D, i::Integer, j::Integer)
    G.xrange[i], G.yrange[j]
end
Base.getindex(G::Grid2D, ij::Tuple{<:Integer, <:Integer}) = G[ij[1], ij[2]]
Base.getindex(G::Grid2D, ij::SVector{2, <:Integer}) = G[ij[1], ij[2]]
Base.getindex(G::Grid2D, xy::NTuple{2, Float64}) = gridpoint(G, xy[1], xy[2])

xcoordinates(G::Grid2D) = G.xrange
ycoordinates(G::Grid2D) = G.yrange
xcoordinate(G::Grid2D, i) = G.xrange[i]
ycoordinate(G::Grid2D, i) = G.yrange[i]

xgridpoint(G::Grid2D, x::Float64) = rangeindex(G.xrange, x)
ygridpoint(G::Grid2D, y::Float64) = rangeindex(G.yrange, y)
gridpoint(G::Grid2D, x, y) = (xgridpoint(G, x), ygridpoint(G, y))
function limits(G::Grid2D)
    xcoords, ycoords = xcoordinates(G), ycoordinates(G)
    xmin = xcoords[1] - 0.5 * step(xcoords)
    xmax = xcoords[end] + 0.5 * step(xcoords)
    ymin = ycoords[1] - 0.5 * step(ycoords)
    ymax = ycoords[end] + 0.5 * step(ycoords)

    xmin, xmax, ymin, ymax
end

"""
    line(f, G, p0, p1)

Map the line between `p0` and `p1` on the grid `G` and call function `f` for each
grid point on that line.
This is an implemenation of
[Bresenhams's line algorithm](https://en.wikipedia.org/wiki/Bresenham%27s_line_algorithm).
"""
function line(f::Function, G::Grid2D, p0, p1)
    x, y = p0
    x1, y1 = p1
    m, n = size(G)

    dx = round(Int, abs(x1-x))
    dy = round(Int, abs(y1-y))

    sx = x < x1 ? 1 : -1
    sy = y < y1 ? 1 : -1;

    err = 0.5 * (dx > dy ? dx : -dy)

    while x != x1 || y != y1
        if 1 ≤ x ≤ m && 1 ≤ y ≤ n
            f(x, y)
        end

        e2 = err;
        if e2 > -dx
            err -= dy
            x += sx
        end
        if e2 < dy
            err += dx
            y += sy
        end
    end
end


text:
function horner(coeff, x0)
	degree = length(coeff) - 1
	y = coeff[1]
	z = 0
	for j ∈ 2:(degree + 1)
		println(coeff[j])
		z = y + x0 * z
		y = x0 * y + coeff[j]
	end
	return (y, z)
end

println(horner([5, 2, 1], 5)[1])


text:
immutable ChebyshevJacobiConstants{D,T}
    α::T
    β::T
    M::Int
    N::Int
    nM₀::Int
    αN::T
    K::Int
end

function ChebyshevJacobiConstants{T}(c::AbstractVector{T},α::T,β::T;M::Int=7,D::Bool=FORWARD)
    N = length(c)-1
    if modαβ(α) == 0.5 && modαβ(β) == 0.5 return ChebyshevJacobiConstants{D,T}(α,β,M,N,0,zero(T),0) end
    if D == FORWARD
        nM₀ = floor(Int,min((eps(T)*2.0^(2M-1)*sqrtpi/absf(modαβ(α),modαβ(β),M,1/2))^(-1/(M+1/2)),N))
        αN = min(one(T)/log(N/nM₀),half(T))
        K = N > 0 ? ceil(Int,log(N/nM₀)/log(1/αN)) : 0
    else#if D == BACKWARD
        nM₀ = floor(Int,min((eps(T)*2.0^(2M-1)*sqrtpi/absf(modαβ(α),modαβ(β),M,1/2))^(-1/(M+1/2)),2N))
        αN = min(one(T)/log(2N/nM₀),half(T))
        K = N > 0 ? ceil(Int,log(2N/nM₀)/log(1/αN)) : 0
    end
    ChebyshevJacobiConstants{D,T}(α,β,M,N,nM₀,αN,K)
end

immutable ChebyshevJacobiIndices
    i₁::Vector{Int}
    i₂::Vector{Int}
    j₁::Vector{Int}
    j₂::Vector{Int}
end

function ChebyshevJacobiIndices{D,T}(α::T,β::T,CJC::ChebyshevJacobiConstants{D,T},tempmindices::Vector{T},cfs::Matrix{T},tempcos::Vector{T},tempsin::Vector{T},tempcosβsinα::Vector{T})
    M,N,nM₀,αN,K = getconstants(CJC)

    i₁,i₂ = zeros(Int,K+1),zeros(Int,K+1)
    j₁,j₂ = zeros(Int,K+1),zeros(Int,K+1)

    if D == FORWARD
        i₁[1],i₂[1],j₂[1] = 1,N+1,N+1
        for k=1:K j₁[k] = ceil(Int,αN^k*N) end
    else#if D == BACKWARD
        i₁[1],i₂[1],j₂[1] = 1,2N+1,2N+1
        for k=1:K j₁[k] = ceil(Int,αN^k*2N) end
    end

    for k=1:K
        if j₁[k] < nM₀ break end
        i₁[k+1],i₂[k+1] = findmindices!(tempmindices,cfs,α,β,j₁[k],M,tempcos,tempsin,tempcosβsinα) # 4 allocations from the Beta function
    end
    for k=1:K j₂[k+1] = j₁[k]-1 end

    ChebyshevJacobiIndices(i₁,i₂,j₁,j₂)
end

type ChebyshevJacobiPlan{D,T,DCT,DST,SA} <: FastTransformPlan{D,T}
    CJC::ChebyshevJacobiConstants{D,T}
    CJI::ChebyshevJacobiIndices
    p₁::DCT
    p₂::DST
    rp::RecurrencePlan{T}
    c₁::Vector{T}
    c₂::SA
    um::Vector{T}
    vm::Vector{T}
    cfs::Matrix{T}
    θ::Vector{T}
    tempcos::Vector{T}
    tempsin::Vector{T}
    tempcosβsinα::Vector{T}
    tempmindices::Vector{T}
    cnαβ::Vector{T}
    cnmαβ::Vector{T}
    w::Vector{T}
    anαβ::Vector{T}
    c_cheb2::Vector{T}
    pr::Vector{T}
    function ChebyshevJacobiPlan(CJC::ChebyshevJacobiConstants{D,T},CJI::ChebyshevJacobiIndices,p₁::DCT,p₂::DST,rp::RecurrencePlan{T},c₁::Vector{T},c₂::SA,um::Vector{T},vm::Vector{T},cfs::Matrix{T},θ::Vector{T},tempcos::Vector{T},tempsin::Vector{T},tempcosβsinα::Vector{T},tempmindices::Vector{T},cnαβ::Vector{T},cnmαβ::Vector{T})
        P = new()
        P.CJC = CJC
        P.CJI = CJI
        P.p₁ = p₁
        P.p₂ = p₂
        P.rp = rp
        P.c₁ = c₁
        P.c₂ = c₂
        P.um = um
        P.vm = vm
        P.cfs = cfs
        P.θ = θ
        P.tempcos = tempcos
        P.tempsin = tempsin
        P.tempcosβsinα = tempcosβsinα
        P.tempmindices = tempmindices
        P.cnαβ = cnαβ
        P.cnmαβ = cnmαβ
        P
    end
    function ChebyshevJacobiPlan(CJC::ChebyshevJacobiConstants{D,T},CJI::ChebyshevJacobiIndices,p₁::DCT,p₂::DST,rp::RecurrencePlan{T},c₁::Vector{T},c₂::SA,um::Vector{T},vm::Vector{T},cfs::Matrix{T},θ::Vector{T},tempcos::Vector{T},tempsin::Vector{T},tempcosβsinα::Vector{T},tempmindices::Vector{T},cnαβ::Vector{T},cnmαβ::Vector{T},w::Vector{T},anαβ::Vector{T},c_cheb2::Vector{T},pr::Vector{T})
        P = ChebyshevJacobiPlan{D,T,DCT,DST,SA}(CJC,CJI,p₁,p₂,rp,c₁,c₂,um,vm,cfs,θ,tempcos,tempsin,tempcosβsinα,tempmindices,cnαβ,cnmαβ)
        P.w = w
        P.anαβ = anαβ
        P.c_cheb2 = c_cheb2
        P.pr = pr
        P
    end
    function ChebyshevJacobiPlan(CJC::ChebyshevJacobiConstants{D,T})
        P = new()
        P.CJC = CJC
        P
    end
end

function ForwardChebyshevJacobiPlan{T}(c_jac::AbstractVector{T},α::T,β::T,M::Int)
    # Initialize constants
    CJC = ChebyshevJacobiConstants(c_jac,α,β;M=M,D=FORWARD)
    if modαβ(α) == 0.5 && modαβ(β) == 0.5 return ChebyshevJacobiPlan{FORWARD,T,Any,Any,Any}(CJC) end
    M,N,nM₀,αN,K = getconstants(CJC)
    α,β = modαβ(α),modαβ(β)

    # Initialize recurrence coefficients
    rp = RecurrencePlan(α,β,N+1)

    # Initialize temporary arrays
    c₁,c₂,um,vm = zero(c_jac),view(zero(c_jac),2:N),zero(c_jac),zero(c_jac)

    # Initialize DCT-I and DST-I plans
    p₁,p₂ = applyTN_plan(c₁),applyUN_plan(c₂)

    # Initialize coefficients of the asymptotic formula
    cfs = init_cfs(α,β,M)

    # Clenshaw-Curtis points
    θ = N > 0 ? T[k/N for k=zero(T):N] : T[0]

    # Initialize sines and cosines
    tempsin = sinpi(θ/2)
    tempcos = reverse(tempsin)
    tempcosβsinα,tempmindices = zero(c_jac),zero(c_jac)
    @inbounds for i=1:N+1 tempcosβsinα[i] = tempcos[i]^(β+1/2)*tempsin[i]^(α+1/2) end

    # Initialize normalizing constant
    cnαβ = Cnαβ(0:N,α,β)
    cnmαβ = zero(cnαβ)

    # Get indices
    CJI = ChebyshevJacobiIndices(α,β,CJC,tempmindices,cfs,tempcos,tempsin,tempcosβsinα)

    ChebyshevJacobiPlan{FORWARD,T,typeof(p₁),typeof(p₂),typeof(c₂)}(CJC,CJI,p₁,p₂,rp,c₁,c₂,um,vm,cfs,θ,tempcos,tempsin,tempcosβsinα,tempmindices,cnαβ,cnmαβ)
end

function BackwardChebyshevJacobiPlan{T}(c_cheb::AbstractVector{T},α::T,β::T,M::Int)
    # Initialize constants
    CJC = ChebyshevJacobiConstants(c_cheb,α,β;M=M,D=BACKWARD)
    if modαβ(α) == 0.5 && modαβ(β) == 0.5 return ChebyshevJacobiPlan{BACKWARD,T,Any,Any,Any}(CJC) end
    M,N,nM₀,αN,K = getconstants(CJC)
    α,β = modαβ(α),modαβ(β)

    # Array of almost double the size of the coefficients
    c_cheb2 = N > 0 ? zeros(T,2N+1) : T[0]

    # Initialize recurrence coefficients and temporary array
    rp,pr = RecurrencePlan(α,β,2N+1),zero(c_cheb2)

    # Initialize temporary arrays
    c₁,c₂,um,vm = zero(c_cheb2),view(zero(c_cheb2),2:2N),zero(c_cheb2),zero(c_cheb2)

    # Initialize DCT-I and DST-I plans
    p₁,p₂ = applyTN_plan(c₁),applyUN_plan(c₂)

    # Initialize coefficients of the asymptotic formula
    cfs = init_cfs(α,β,M)

    # Clenshaw-Curtis nodes and weights
    θ = N > 0 ? T[k/2N for k=zero(T):2N] : T[0]
    w = N > 0 ? clenshawcurtisweights(2N+1,α,β,p₁) : T[0]

    # Initialize sines and cosines
    tempsin = sinpi(θ/2)
    tempcos = reverse(tempsin)
    tempcosβsinα,tempmindices = zero(c_cheb2),zero(c_cheb2)
    @inbounds for i=1:2N+1 tempcosβsinα[i] = tempcos[i]^(β+1/2)*tempsin[i]^(α+1/2) end

    # Initialize normalizing constant
    cnαβ = Cnαβ(0:2N,α,β)
    cnmαβ = zero(cnαβ)

    # Initialize orthonormality constants
    anαβ = Anαβ(0:N,α,β)

    # Get indices
    CJI = ChebyshevJacobiIndices(α,β,CJC,tempmindices,cfs,tempcos,tempsin,tempcosβsinα)

    ChebyshevJacobiPlan{BACKWARD,T,typeof(p₁),typeof(p₂),typeof(c₂)}(CJC,CJI,p₁,p₂,rp,c₁,c₂,um,vm,cfs,θ,tempcos,tempsin,tempcosβsinα,tempmindices,cnαβ,cnmαβ,w,anαβ,c_cheb2,pr)
end

getplanαβ(CJC::ChebyshevJacobiConstants) = CJC.α,CJC.β
getplanαβ(plan::ChebyshevJacobiPlan) = getplanαβ(plan.CJC)

getconstants(CJC::ChebyshevJacobiConstants) = CJC.M,CJC.N,CJC.nM₀,CJC.αN,CJC.K
getconstants(plan::ChebyshevJacobiPlan) = getconstants(plan.CJC)

getindices(CJI::ChebyshevJacobiIndices) = CJI.i₁,CJI.i₂,CJI.j₁,CJI.j₂
getindices(plan::ChebyshevJacobiPlan) = getindices(plan.CJI)

getplan(plan::ChebyshevJacobiPlan{FORWARD})  = plan.p₁,plan.p₂,plan.rp,plan.c₁,plan.c₂,plan.um,plan.vm,plan.cfs,plan.θ,plan.tempcos,plan.tempsin,plan.tempcosβsinα,plan.tempmindices,plan.cnαβ,plan.cnmαβ
getplan(plan::ChebyshevJacobiPlan{BACKWARD}) = plan.p₁,plan.p₂,plan.rp,plan.c₁,plan.c₂,plan.um,plan.vm,plan.cfs,plan.θ,plan.tempcos,plan.tempsin,plan.tempcosβsinα,plan.tempmindices,plan.cnαβ,plan.cnmαβ,plan.w,plan.anαβ,plan.c_cheb2,plan.pr


text:
# ------------------------------------------------------------------------------------------
# # Introduction to DataFrames
# **[Bogumił Kamiński](http://bogumilkaminski.pl/about/), May 23, 2018**
# ------------------------------------------------------------------------------------------

using DataFrames # load package

# ------------------------------------------------------------------------------------------
# ## Manipulating columns of DataFrame
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Renaming columns
#
# Let's start with a `DataFrame` of `Bool`s that has default column names.
# ------------------------------------------------------------------------------------------

x = DataFrame(Bool, 3, 4)

# ------------------------------------------------------------------------------------------
# With `rename`, we create new `DataFrame`; here we rename the column `:x1` to `:A`.
# (`rename` also accepts collections of Pairs.)
# ------------------------------------------------------------------------------------------

rename(x, :x1 => :A)

# ------------------------------------------------------------------------------------------
# With `rename!` we do an in place transformation.
#
# This time we've applied a function to every column name.
# ------------------------------------------------------------------------------------------

rename!(c -> Symbol(string(c)^2), x)

# ------------------------------------------------------------------------------------------
# We can also change the name of a particular column without knowing the original.
#
# Here we change the name of the third column, creating a new `DataFrame`.
# ------------------------------------------------------------------------------------------

rename(x, names(x)[3] => :third)

# ------------------------------------------------------------------------------------------
# With `names!`, we can change the names of all variables.
# ------------------------------------------------------------------------------------------

names!(x, [:a, :b, :c, :d])

# ------------------------------------------------------------------------------------------
# We get an error when we try to provide duplicate names
# ------------------------------------------------------------------------------------------

names!(x, fill(:a, 4))

# ------------------------------------------------------------------------------------------
#  unless we pass `makeunique=true`, which allows us to handle duplicates in passed names.
# ------------------------------------------------------------------------------------------

names!(x, fill(:a, 4), makeunique=true)

# ------------------------------------------------------------------------------------------
# ### Reordering columns
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We can reorder the names(x) vector as needed, creating a new DataFrame.
# ------------------------------------------------------------------------------------------

srand(1234)
x[shuffle(names(x))]

# ------------------------------------------------------------------------------------------
# also `permutecols!` will be introduced in next release of DataFrames
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Merging/adding columns
# ------------------------------------------------------------------------------------------

x = DataFrame([(i,j) for i in 1:3, j in 1:4])

# ------------------------------------------------------------------------------------------
# With `hcat` we can merge two `DataFrame`s. Also [x y] syntax is supported but only when
# DataFrames have unique column names.
# ------------------------------------------------------------------------------------------

hcat(x, x, makeunique=true)

# ------------------------------------------------------------------------------------------
# We can also use `hcat` to add a new column; a default name `:x1` will be used for this
# column, so `makeunique=true` is needed.
# ------------------------------------------------------------------------------------------

y = hcat(x, [1,2,3], makeunique=true)

# ------------------------------------------------------------------------------------------
# You can also prepend a vector with `hcat`.
# ------------------------------------------------------------------------------------------

hcat([1,2,3], x, makeunique=true)

# ------------------------------------------------------------------------------------------
# Alternatively you could append a vector with the following syntax. This is a bit more
# verbose but cleaner.
# ------------------------------------------------------------------------------------------

y = [x DataFrame(A=[1,2,3])]

# ------------------------------------------------------------------------------------------
# Here we do the same but add column `:A` to the front.
# ------------------------------------------------------------------------------------------

y = [DataFrame(A=[1,2,3]) x]

# ------------------------------------------------------------------------------------------
# A column can also be added in the middle. Here a brute-force method is used and a new
# DataFrame is created.
# ------------------------------------------------------------------------------------------

using BenchmarkTools
@btime [$x[1:2] DataFrame(A=[1,2,3]) $x[3:4]]

# ------------------------------------------------------------------------------------------
# We could also do this with a specialized in place method `insert!`. Let's add `:newcol` to
# the `DataFrame` `y`.
# ------------------------------------------------------------------------------------------

insert!(y, 2, [1,2,3], :newcol)

# ------------------------------------------------------------------------------------------
# If you want to insert the same column name several times `makeunique=true` is needed as
# usual.
# ------------------------------------------------------------------------------------------

insert!(y, 2, [1,2,3], :newcol, makeunique=true)

# ------------------------------------------------------------------------------------------
# We can see how much faster it is to insert a column with `insert!` than with `hcat` using
# `@btime`.
# ------------------------------------------------------------------------------------------

@btime insert!(copy($x), 3, [1,2,3], :A)

# ------------------------------------------------------------------------------------------
# Let's use `insert!` to append a column in place,
# ------------------------------------------------------------------------------------------

insert!(x, ncol(x)+1, [1,2,3], :A)

# ------------------------------------------------------------------------------------------
# and to in place prepend a column.
# ------------------------------------------------------------------------------------------

insert!(x, 1, [1,2,3], :B)

# ------------------------------------------------------------------------------------------
# With `merge!`, let's merge the second DataFrame into first, but overwriting duplicates.
# ------------------------------------------------------------------------------------------

df1 = DataFrame(x=1:3, y=4:6)
df2 = DataFrame(x='a':'c', z = 'd':'f', new=11:13)
df1, df2, merge!(df1, df2)

# ------------------------------------------------------------------------------------------
#  For comparison: merge two `DataFrames`s but renaming duplicate names via `hcat`.
# ------------------------------------------------------------------------------------------

df1 = DataFrame(x=1:3, y=4:6)
df2 = DataFrame(x='a':'c', z = 'd':'f', new=11:13)
hcat(df1, df2, makeunique=true)

# ------------------------------------------------------------------------------------------
# ### Subsetting/removing columns
#
# Let's create a new `DataFrame` `x` and show a few ways to create DataFrames with a subset
# of `x`'s columns.
# ------------------------------------------------------------------------------------------

x = DataFrame([(i,j) for i in 1:3, j in 1:5])

# ------------------------------------------------------------------------------------------
# First we could do this by index
# ------------------------------------------------------------------------------------------

x[[1,2,4,5]]

# ------------------------------------------------------------------------------------------
# or by column name.
# ------------------------------------------------------------------------------------------

x[[:x1, :x4]]

# ------------------------------------------------------------------------------------------
# We can also choose to keep or exclude columns by `Bool`. (We need a vector whose length is
# the number of columns in the original `DataFrame`.)
# ------------------------------------------------------------------------------------------

x[[true, false, true, false, true]]

# ------------------------------------------------------------------------------------------
# Here we create a single column `DataFrame`,
# ------------------------------------------------------------------------------------------

x[[:x1]]

# ------------------------------------------------------------------------------------------
# and here we access the vector contained in column `:x1`.
# ------------------------------------------------------------------------------------------

x[:x1]

# ------------------------------------------------------------------------------------------
# We could grab the same vector by column number
# ------------------------------------------------------------------------------------------

x[1]

# ------------------------------------------------------------------------------------------
# and remove everything from a `DataFrame` with `empty!`.
# ------------------------------------------------------------------------------------------

empty!(y)

# ------------------------------------------------------------------------------------------
# Here we create a copy of `x` and delete the 3rd column from the copy with `delete!`.
# ------------------------------------------------------------------------------------------

z = copy(x)
x, delete!(z, 3)

# ------------------------------------------------------------------------------------------
# ### Modify column by name
# ------------------------------------------------------------------------------------------

x = DataFrame([(i,j) for i in 1:3, j in 1:5])

# ------------------------------------------------------------------------------------------
# With the following syntax, the existing column is modified without performing any copying.
# ------------------------------------------------------------------------------------------

x[:x1] = x[:x2]
x

# ------------------------------------------------------------------------------------------
# We can also use the following syntax to add a new column at the end of a `DataFrame`.
# ------------------------------------------------------------------------------------------

x[:A] = [1,2,3]
x

# ------------------------------------------------------------------------------------------
# A new column name will be added to our `DataFrame` with the following syntax as well (7 is
# equal to `ncol(x)+1`).
# ------------------------------------------------------------------------------------------

x[7] = 11:13
x

# ------------------------------------------------------------------------------------------
# ### Find column name
# ------------------------------------------------------------------------------------------

x = DataFrame([(i,j) for i in 1:3, j in 1:5])

# ------------------------------------------------------------------------------------------
# We can check if a column with a given name exists via
# ------------------------------------------------------------------------------------------

:x1 in names(x)

# ------------------------------------------------------------------------------------------
# and determine its index via
# ------------------------------------------------------------------------------------------

findfirst(names(x), :x2)

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
include(joinpath("src", "simulation.jl"))
include(joinpath("src", "03runeval.jl"))

run_batch(batch_netsize, batch_name = "Netsize", stop_at = 1)
run_batch(batch_addfriends, batch_name = "Addfriends")
run_batch(batch_unfriend, batch_name = "Unfriend", resume_at=5)

raw = load("results/Addfriends_run1.jld2")

addfriends = raw[first(keys(raw))]

netsize
RunEval(netsize[51])

teest = rerun_single(netsize[51])

RunEval(teest)

batch_netsize = Config[]

for agent_count in 100:100:500

    push!(batch_netsize, Config(
            network = cfg_net(
                agent_count = agent_count,
                m0 = Int(agent_count/10),
                new_follows = 10
            ),
            simulation = cfg_sim(
                ticks = 1000,
                addfriends = "hybrid",
                repcount = 2,
                agent_logging = false
            ),
            opinion_threshs = cfg_ot(
                backfire = 0.4,
                befriend = 0.2,
                unfriend = 0.6
            ),
            agent_props = cfg_ag(
                own_opinion_weight = 0.95,
                unfriend_rate = 0.05,
                min_friends_count = 5
            )
        )
    )
end

batch_unfriend = Config[]

for unfriend in [0.4, 0.6, 0.8, 1.0, 1.2]

    push!(batch_unfriend, Config(
            network = cfg_net(
                agent_count = 500,
                m0 = Int(500/10),
                new_follows = 10
            ),
            simulation = cfg_sim(
                ticks = 1000,
                addfriends = "hybrid",
                repcount = 50,
                agent_logging = false
            ),
            opinion_threshs = cfg_ot(
                backfire = 0.4,
                befriend = 0.2,
                unfriend = unfriend,
            ),
            agent_props = cfg_ag(
                own_opinion_weight = 0.95,
                unfriend_rate = 0.05,
                min_friends_count = 5
            )
        )
    )
end

batch_addfriends = Config[]

for addfriends in ["neighborsofneighbors", "hybrid", "random"]

    push!(batch_addfriends, Config(
            network = cfg_net(
                agent_count = 500,
                m0 = Int(500/10),
                new_follows = 10
            ),
            simulation = cfg_sim(
                ticks = 1000,
                addfriends = addfriends,
                repcount = 50,
                agent_logging = false
            ),
            opinion_threshs = cfg_ot(
                backfire = 0.4,
                befriend = 0.2,
                unfriend = 0.6
            ),
            agent_props = cfg_ag(
                own_opinion_weight = 0.95,
                unfriend_rate = 0.05,
                min_friends_count = 5
            )
        )
    )
end


text:
# Conservative 7-color palette
# Wong, Bang. "Points of view: Color blindness." (2011): 441.
# https://www.nature.com/articles/nmeth.1618?WT.ec_id=NMETH-201106

function default_palettes()
    return (
        color=[
            RGB(0/255, 114/255, 178/255), # blue
            RGB(230/255, 159/255, 0/255), # orange
            RGB(0/255, 158/255, 115/255), # green
            RGB(204/255, 121/255, 167/255), # reddish purple
            RGB(86/255, 180/255, 233/255), # sky blue
            RGB(213/255, 94/255, 0/255), # vermillion
            RGB(240/255, 228/255, 66/255), # yellow
        ],
        marker=[:circle, :utriangle, :cross, :rect, :diamond, :dtriangle, :pentagon, :xcross],
        linestyle=[:solid, :dash, :dot, :dashdot, :dashdotdot],
        side=[:left, :right],
        layout=wrap,
    )
end

# Batlow colormap
# Crameri, Fabio, Grace E. Shephard, and Philip J. Heron. "The misuse of colour in science communication." Nature communications 11.1 (2020): 1-10.
# https://www.nature.com/articles/s41467-020-19160-7

function default_styles()
    return (
        color=:gray25,
        strokecolor=RGBA(0, 0, 0, 0),
        outlierstrokecolor=RGBA(0, 0, 0, 0),
        mediancolor=:white,
        marker=:circle,
        markersize=9,
        linewidth=1.5,
        medianlinewidth=1.5,
        colormap=:batlow,
    )
end

# axis defaults

const font_folder = joinpath(dirname(@__DIR__), "assets", "fonts")

firasans(weight) = joinpath(font_folder, "FiraSans-$(weight).ttf")
opensans(weight) = joinpath(font_folder, "OpenSans-$(weight).ttf")

function aog_theme()
    Axis = (
        xgridvisible=false,
        ygridvisible=false,
        topspinevisible=false,
        rightspinevisible=false,
        bottomspinecolor=:darkgray,
        leftspinecolor=:darkgray,
        xtickcolor=:darkgray,
        ytickcolor=:darkgray,
        xticklabelfont=firasans("Light"),
        yticklabelfont=firasans("Light"),
        xlabelfont=firasans("Medium"),
        ylabelfont=firasans("Medium"),
        titlefont=firasans("Medium"),
    )
    Axis3 = (
        protrusions=55, # to include label on z axis, should be fixed in AbstractPlotting
        xgridvisible=false,
        ygridvisible=false,
        zgridvisible=false,
        xspinecolor=:darkgray,
        yspinecolor=:darkgray,
        zspinecolor=:darkgray,
        xtickcolor=:darkgray,
        ytickcolor=:darkgray,
        ztickcolor=:darkgray,
        xticklabelfont=firasans("Light"),
        yticklabelfont=firasans("Light"),
        zticklabelfont=firasans("Light"),
        xlabelfont=firasans("Medium"),
        ylabelfont=firasans("Medium"),
        zlabelfont=firasans("Medium"),
        titlefont=firasans("Medium"),
    )
    Legend = (
        framevisible=false,
        gridshalign=:left,
        padding=(0f0, 0f0, 0f0, 0f0),
        labelfont=firasans("Light"),
        titlefont=firasans("Medium"),
    )
    Colorbar = (
        colormap=:batlow,
        flip_vertical_label=true,
        spinewidth=0,
        ticklabelfont=firasans("Light"),
        labelfont=firasans("Medium"),
    )
    return (; Axis, Axis3, Legend, Colorbar)
end

function set_aog_theme!()
    theme = aog_theme()
    return set_theme!(; theme...)
end


text:
const _GW_REG_KEY = "_register"
const _GW_REG_RFILE_KEY = "_rfkey"
const _GW_REG_UPTIME_KEY = "_uptime"

function reg_sdat!(gw::GitWorker, fkey::String, dat)
    rfkey = worker_relpath(gw, fkey)

    # ram version
    keys = splitpath(rfkey)
    sdat = get_sreg(gw)
    for k in keys
        sdat = get!(sdat, k, Dict{String, Any}())
    end

    for (dk, val) in dat
        sdat[string(dk)] = val
    end
    sdat[_GW_REG_RFILE_KEY] = rfkey # cross-machine compatible
    sdat[_GW_REG_UPTIME_KEY] = time()

    return sdat
end

reg_sdat!(gw::GitWorker, fkey::String; dat...) = reg_sdat!(gw::GitWorker, fkey::String, dat)

get_sreg(gw::GitWorker) = get!(() -> Dict{String, Any}(), gw, _GW_REG_KEY)

function get_sdat(gw::GitWorker, rkey::String)
    rfkey = worker_relpath(gw, rkey)

    keys = splitpath(rfkey)
    dict = get_sreg(gw)
    for k in keys
        dict = get!(dict, k) do
            Dict{String, Any}()
        end
    end
    return dict
end

function empty_sreg!(gw::GitWorker, rkey::String)
    sdat = get_sdat(gw, rkey)
    (sdat isa Dict) && empty!(sdat)
    return sdat
end

function _find_by_khint(gw::GitWorker, rootk, khint)
    khint = string(khint)
    sdat = get_sdat(gw, rootk)
    for (k, dat) in sdat
        contains(k, khint) && return dat
    end
    return Dict{String, Any}()
end

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
using MultilevelEstimators
using DifferentialEquations
using Reporter
using FileIO
using DelimitedFiles

DIR=pwd()*"/MultilevelEstimators/data"

# DifferentialEquations SDE example partly taken from https://diffeq.sciml.ai/stable/tutorials/sde_example/

# define the problem
α=0.05
β=0.2
u₀=100.
f(u,p,t) = α*u
g(u,p,t) = β*u
dt = 1/2^(4)
T = 1.
tspan = (0.0, T)
P_func = (x) -> exp(-α*T)*maximum([0, x-100.])
prob = SDEProblem(f,g,u₀,(0.0,1.0), dt=dt)

sol = solve(prob,EM())
using Plots; plotly() # Using the Plotly backend
p1 = plot(sol)

MEAN = 0
P = 10^4
for p=1:P
    sol = solve(prob,EM())
    MEAN += sol[end]
end
MEAN /= P

# The mean of the solution at the end, should equal ???

# Multilevelusing FileIOEstimators.jl setup, adapted from documentation

LL = 20
SDE_probs = zeros(LL)
for i in 1:LL
    dt = 1//2^(i+1)
    SDE_probs[i] = dt
end

function sample_SDE(level, ω; TEST=false)
    # solve on finest grid
    SDE_prob_coarse = SDEProblem(f, g, u₀, (0.0, 1.0), dt=SDE_probs[level + one(level)])
    sol1 = solve(SDE_prob_coarse, EM(), save_noise = true) # TO DO, dt needs to be varied
    Qf = P_func(sol1.u[end])

    # compute difference when not on coarsest grid
    dQ = copy(Qf)
    if level != Level(0)
        SDE_prob_fine = SDEProblem(f, g, u₀, (0.0, 1.0), dt=SDE_probs[level], noise=NoiseWrapper(sol1.W))
        sol2 = solve(SDE_prob_fine, EM())
        Qc = P_func(sol2.u[end])
        dQ -= Qc
    end

    if TEST
        return dQ, Qf, Qc
    else
        return dQ, Qf
    end
end

distributions = [Normal()]
estimator = Estimator(ML(), MC(), sample_SDE, distributions, folder=DIR)

# run Estimator
h = run(estimator, 5e-2)

S = MultilevelEstimators.samples_diff(estimator)
name = estimator[:name]
sf_dir= string(name[1:end-5])

isdir(joinpath(DIR,sf_dir)) || mkdir(joinpath(DIR,sf_dir))
isdir(joinpath(DIR,sf_dir,"Samples")) || mkdir(joinpath(DIR,sf_dir,"Samples"))

for idx in CartesianIndices(S)
    idx_dir = joinpath(DIR,sf_dir,"Samples", join(idx.I,"_"))
    isdir(idx_dir) || mkdir(idx_dir)
    for k in keys(S[idx])
        writedlm(joinpath(idx_dir,string("samples_level_",k[1]-1,".txt")),S[idx][k])
    end
end

full_history = History(h)

filePath = DIR*"/UntitledEstimator.jld2"
history = load(filePath,"history")
report(history, folder=DIR*"/testReport", include_preamble=true)


text:
# Text provided under a Creative Commons Attribution license, CC-BY.  All code is made available under the FSF-approved BSD-3 license.  (c) Lorena A. Barba, Gilbert F. Forsyth 2017. Thanks to NSF for support via CAREER award #1149784.
# [@LorenaABarba](https://twitter.com/LorenaABarba)
# Translation to Julia by [@miguelraz_](https://twitter.com/miguelraz_)

# 12 steps to Navier–Stokes
# =====
# ***

# This lesson complements the first interactive module of the online [CFD Python](https://github.com/barbagroup/CFDPython) class,
# by Prof. Lorena A. Barba, called **12 Steps to Navier–Stokes.** The interactive module starts with
# simple exercises in 1D that at first use little of the power of Python.
# We now present some new ways of doing the same things that are more efficient and produce prettier code, but in Julia 🚀 ,
# which is hosted at the [CFD Julia](TODO LINK) repo

# This lesson was written with BU graduate student Gilbert Forsyth.

# Defining Functions in Julia
# ----

# In steps 1 through 8, we wrote Python code that is meant to run from top to bottom.
# We were able to reuse code (to great effect!) by copying and pasting, to incrementally build a solver for the Burgers' equation.
# But moving forward there are more efficient ways to write our Python codes.
# In this lesson, we are going to introduce *function definitions*,
# which will allow us more flexibility in reusing and also in organizing our code.

# We'll begin with a trivial example: a function which adds two numbers.

# To create a function in Julia, we start with the following function block:

function simpleadd(a,b)
    return a + b
end

# The `return` statement tells Julia what data to return in response to being called.
# In Julia, if you omit the `return` inside the body of a funciton, the last evaluated expression is implicitly returned.
# Now we can try calling our `simpleadd` function:

simpleadd(3, 4)

# For simple functions, we can use the inline function syntax:
simpleadd(a,b) = a + b

# Where both are equivalent.

# Of course, there can be much more happening between the `function` line and the `return` line.
# In this way, one can build code in a *modular* way.
# Let's try a function which returns the `n`-th number in the Fibonacci sequence.

function fibonacci(n)
    a, b = 0, 1
    for i in 1:n
        a, b = b, a + b
    end
    return a
end

fibonacci(7)

# Once defined, the function `fibonacci` can be called like any of the built-in Python functions that we've already used.
# For exmaple, we might want to print out the Fibonacci sequence up through the `n`-th value:

for n in 1:10
    println(fibonacci(n))
end

# We will use the capacity of defining our own functions in Julia to help us build code that is easier to reuse, easier to maintain, easier to share!

##### Exercise

# (Pending.)

# Learn more
# -----
# ***

# Remember our short detour on using [array operations with Julia](./07_Step_5.jl)?

# Well, there are a few more ways to make your scientific codes in Julia run faster.
# The original tutorial said...
# > We recommend the article on the Technical Discovery blog about [Speeding Up Python](http://technicaldiscovery.blogspot.com/2011/06/speeding-up-python-numpy-cython-and.html) (June 20, 2011), which talks about NumPy, Cython and Weave. It uses as example the Laplace equation (which we will solve in [Step 9](./12_Step_9.ipynb)) and makes neat use of defined functions.
#  We *encourage* instead that you read [the performance tips section of the Julia manual](https://docs.julialang.org/en/v1/manual/performance-tips/).

# From the original tutorial, it says...
# > But a recent new way to get fast Python codes is [Numba](http://numba.pydata.org).
# Which is actually a good comparison to how Julia's JIT operates internally, but unfortunately breaks down in many cases
# [that we are interested in](http://www.stochasticlifestyle.com/why-numba-and-cython-are-not-substitutes-for-julia/)

max_stars_repo_name: mkitti/CFDJulia

max_issues_repo_name: mkitti/CFDJulia

max_forks_repo_name: mkitti/CFDJulia


text:
"""
	SpinModule

Store informations about the spin model.

# Type list

* SpinInfo
* TrotterInfo for Quantum Spins
* SpinModel
* ClassicalSpinModel
* IsingModel
* PottsModel
* ClockModel
* QuantumSpinModel
* QuantumIsingModel

# Method list

* Methods for SpinInfo:
  * isZeroTemperature
  * isClassical
  * getmodelname
  * getStates
  * isSymmetricFactorization
  * getTemperature
  * getExternalfield
  * getEnvParameters
  * setTemperature!
  * setExternalfield!
  * setEnvParameters!
* Methods for TrotterInfo
  * getTrotterparameter
  * getTrotteriteration
  * getTrotterlayers
  * getBeta
  * setTrotterparameter!
  * setTrotteriteration!
  * iteration2layer
* Methods for SpinModel
  * buildSpinSystem
  * getMeasureOperator
  * getFactorW
  * isPottsModel
  * getHamiltonian
* Methods for IsingModel
  * initialize!
* Methods for QuantumSpinModel
  * checkInfo!
* Methods for QuantumIsingModel
  * getFactorWp

# Test list

* testSpinInfo
* testTrotterInfo
* testSpinModel
* testIsingModel
* testPottsModel
* testClockModel
* testQuantumIsingModel
"""
module SpinModule

export
## spin_info.jl
SpinInfo,
isZeroTemperature,
isClassical,
getModelname,
getStates,
isSymmetricFactorization,
getTemperature,
getExternalfield,
getEnvParameters,
setTemperature!,
setExternalfield!,
setEnvParameters!,
testSpinInfo,
## spin_trotter.jl
TrotterInfo,
getTrotterparameter,
getTrotteriteration,
getTrotterlayers,
getBeta,
setTrotterparameter!,
setTrotterIteration!,
iteration2layer,
testTrotterInfo,
	## spin_model.jl
SpinModel,
buildSpinSystem,
getMeasureOperator,
getFactorW,
isPottsModel,
getHamiltonian,
testSpinModel,
	## spin_classical.jl
ClassicalSpinModel,
getBeta,
	## spin_classical_ising.jl
IsingModel,
initialize!,
testIsingModel,
	## spin_classical_potts.jl
PottsModel,
testPottsModel,
	## spin_classical_clock.jl
ClockModel,
testClockModel,
	## spin_quantum.jl
QuantumSpinModel,
checkInfo!,
	## spin_quantum_ising.jl
QuantumIsingModel,
getFactorWp,
testQuantumIsingModel

include("spin_info.jl")
include("spin_trotter.jl")
include("spin_model.jl")
include("spin_classical.jl")
include("spin_classical_ising.jl")
include("spin_classical_potts.jl")
include("spin_classical_clock.jl")
include("spin_quantum.jl")
include("spin_quantum_ising.jl")

end


text:
using VectorSphericalWaves
using ProfileView

m, n, kr_r, kr_i, θ, ϕ = 1, 1, 1e5, 0, 0.2, 0.3
kr = complex(kr_r, kr_i)

@profiler VectorSphericalWaves.πₘₙ( m, n, θ)


text:
let i, n, A, B, C, θ
n = 5 # rand(1:5)

A = matrixdepot("rando", Int, n, 1)
B = matrixdepot("rando", Int, n, 2)
C = matrixdepot("rando", Int, n, 3)

for i in length(A)
    @test A[i] == 0 || A[i] == 1
end

@test all([bi == -1 || bi == 1 for bi in B])
@test all([bi == -1 || bi == 0 || bi == 1 for bi in C])

@test_throws ArgumentError matrixdepot("rando", 3, 3, 5)

θ = matrixdepot("rando", n)
end
println("'rando' passed test...")