text:
export PointSymmetry

export group_order,
       element, elements,
       element_name, element_names,
       group_multiplication_table,
       character_table,
       irrep, irreps, num_irreps, irrep_dimension,
       generator_elements, generator_indices,
       symmetry_product

export project

export little_group_elements
export little_group
export little_symmetry, little_symmetry_iso

export get_irrep_iterator
export read_point_symmetry
export symmetry_name

import LinearAlgebra
import MathExpr
using GroupTools

simplify_name(name::AbstractString) = replace(replace(name, r".*?"=>""), r".*?"=>"")

## PointSymmetry, constructor and related functions

"""
    PointSymmetry

# Arguments
* `elements::Vector{PointOperation{Int}}`
* `group::FiniteGroup`
* `generators::Vector{Int}`
* `conjugacy_classes::Vector{Vector{Int}}`
* `character_table::Matrix{ComplexF64}`
* `irreps::Vector{Vector{Matrix{ComplexF64}}}`
* `element_names::Vector{String}`
* `matrix_representations::Vector{Matrix{Int}}`
* `hermann_mauguin::String`
* `schoenflies::String`
"""
struct PointSymmetry <: AbstractSymmetry{PointOperation{Int}}
    elements::Vector{PointOperation{Int}}
    group::FiniteGroup

    generators::Vector{Int}

    conjugacy_classes::Vector{Vector{Int}}
    character_table::Matrix{ComplexF64}
    irreps::Vector{Vector{Matrix{ComplexF64}}}

    element_names::Vector{String}

    matrix_representations::Vector{Matrix{Int}}
    hermann_mauguin::String
    schoenflies::String

    function PointSymmetry(
        group::FiniteGroup,
        generators::AbstractVector{<:Integer},
        conjugacy_classes::AbstractVector{<:AbstractVector{<:Integer}},
        character_table::AbstractMatrix{<:Number},
        irreps::AbstractVector{<:AbstractVector{<:AbstractMatrix{<:Number}}},
        element_names::AbstractVector{<:AbstractString},
        matrix_representations::AbstractVector{<:AbstractMatrix{<:Integer}},
        hermann_mauguin::AbstractString,
        schoenflies::AbstractString;
        tol::Real=Base.rtoldefault(Float64)
    )
        let dim = size(first(matrix_representations), 1)
            if any(size(m) != (dim, dim) for m in matrix_representations)
                throw(DimensionMismatch("irrep matrix representations should be square matrices of the same dimension"))
            end
        end
        # number counting check
        if length(element_names) != group_order(group)
            throw(ArgumentError("number of elements different from order of group"))
        elseif length(matrix_representations) != group_order(group)
            throw(ArgumentError("number of matrix representations different from order of group"))
        elseif generate_subgroup(group, generators) != BitSet(1:group_order(group))
            throw(ArgumentError("Generators $(generators) does not generate the group"))
        elseif !ishomomorphic(matrix_representations, group)
            msg = """
            matrix representation not homomorphic to group:
            $matrix_representations
            $group
            """
            throw(ArgumentError(msg))
        end

        let prev_set = BitSet()
            for elements in conjugacy_classes
                !isempty(intersect(elements, prev_set)) && throw(ArgumentError("Same element in multiple conjugacy classes"))
                union!(prev_set, elements)
            end
            prev_set != BitSet(1:group_order(group)) && throw(ArgumentError("every element must belong to a conjugacy class"))
        end
        let n = length(conjugacy_classes)
            if size(character_table) != (n, n)
                throw(ArgumentError("size of character table does not match the number of conjugacy classes"))
            end
            if length(irreps) != n
                throw(ArgumentError("number of irreps match the number of conjugacy classes"))
            end
        end
        for rep in irreps
            if length(rep) != group_order(group)
                throw(ArgumentError("wrong number of matrices in irrep $rep"))
            end
            d = size(rep[1], 1)
            if !isapprox(rep[1], Matrix(LinearAlgebra.I, (d, d)); atol=tol)
                throw(ArgumentError("irrep matrix representation of identity should be identity"))
            end
            for m in rep
                if size(m) != (d,d)
                    throw(ArgumentError("irrep matrix representation should all have the same dimension"))
                end
            end
            if !ishomomorphic(rep, group; equal=(x,y) -> isapprox(x,y;atol=tol))
                throw(ArgumentError("irrep matrix representation not homomorphic to group"))
            end
        end

        if !allunique(matrix_representations)
            @warn "matrix representations not all unique"
        end

        elements = PointOperation.(matrix_representations)

        return new(
            elements,
            group,
            generators,
            conjugacy_classes,
            character_table,
            irreps,
            element_names,
            matrix_representations,
            hermann_mauguin,
            schoenflies,
        )
    end
end

"""
    read_point_symmetry(data::AbstractDict)
"""
function read_point_symmetry(data::AbstractDict)
    tol = Base.rtoldefault(Float64)
    read_matrix(obj) = collect(transpose(hcat(MathExpr.parseexpr(obj)...)))

    multiplication_table = read_matrix(data["MultiplicationTable"])
    group = FiniteGroup(multiplication_table)
    ord_group = group_order(group)
    generators = data["Generators"]
    conjugacy_classes = [item for item in data["ConjugacyClasses"]]
    character_table = MathExpr.cleanupnumber(read_matrix(data["CharacterTable"]), tol)
    irreps = [
        Matrix{ComplexF64}[MathExpr.cleanupnumber(read_matrix(elem), tol) for elem in item]
        for item in data["IrreducibleRepresentations"]
    ]
    element_names = data["ElementNames"]
    matrix_representations = Matrix{Int}[read_matrix(x) for x in data["MatrixRepresentations"]]
    hermann_mauguin = data["HermannMauguin"]
    schoenflies = data["Schoenflies"]
    return PointSymmetry(
        group, generators,
        conjugacy_classes, character_table, irreps,
        element_names, matrix_representations, hermann_mauguin, schoenflies,
    )
end

"""
    project(psym, projection; tol=√ϵ)

Project `psym` to a subspace using `projection`

# Arguments
* `psym::PointSymmetry`
* `projection::AbstractMatrix{<:Integer}`: n × d matrix, where d ≤ n.

# Optional Arguments
* `tol::Real`: singular values of projection must be 1 up to tolerance.
"""
function project(
    psym::PointSymmetry,
    projection::AbstractMatrix{<:Integer};
    tol::Real=Base.rtoldefault(Float64),
)
    dim = size(psym.matrix_representations[1], 1)
    if size(projection, 2) != dim
        throw(ArgumentError("projection does not match matrix_representations dimension"))
    end

    vals = LinearAlgebra.svdvals(projection)
    if any(x -> abs(x - 1) > tol, vals)
        throw(ArgumentError("projection has non-unity singular value: $vals"))
    end

    new_matrix_representations = [
        projection * m * transpose(projection)
        for m in psym.matrix_representations
    ]

    return PointSymmetry(
        psym.group,
        psym.generators,
        psym.conjugacy_classes,
        psym.character_table,
        psym.irreps,
        psym.element_names,
        new_matrix_representations,
        psym.hermann_mauguin,
        psym.schoenflies,
    )
end

## Basic properties

"""
    dimension(sym::PointSymmetry)

Spatial dimension of the point symmetry
"""
dimension(sym::PointSymmetry) = size(sym.matrix_representations[1], 1)

"""
    elements(sym::PointSymmetry)

Return the elements of the point symmetry.
"""
GroupTools.elements(sym::PointSymmetry) = sym.elements

"""
    element(sym::PointSymmetry, i)

Return the `i`th element of the point symmetry.
"""
GroupTools.element(sym::PointSymmetry, g) = sym.elements[g]

"""
    element_names(sym::PointSymmetry)

Return the names of the elements of the point symmetry.
"""
GroupTools.element_names(sym::PointSymmetry) = sym.element_names

"""
    element_name(sym::PointSymmetry, i)

Return the name of the `i`th element of the point symmetry.
"""
GroupTools.element_name(sym::PointSymmetry, g) = sym.element_names[g]

"""
    group(sym::PointSymmetry)

Group structure of the point symmetry
"""
GroupTools.group(psym::PointSymmetry) = psym.group

"""
    group_order(sym::PointSymmetry, g...)

Group order of the point symmetry. Calls `group_order(group(sym), g...)`
"""
GroupTools.group_order(psym::PointSymmetry, g...) = group_order(group(psym), g...)

"""
    group_multiplication_table(sym::PointSymmetry)

Group multiplication table of the point symmetry.
Calls `group_multiplication_table(group(sym))`
"""
GroupTools.group_multiplication_table(sym::PointSymmetry) = group_multiplication_table(group(sym))

"""
    generator_indices(sym::PointSymmetry)

Return indices of the generating translations.
"""
generator_indices(sym::PointSymmetry) = sym.generators

"""
    generator_elements(sym::PointSymmetry)

Return the generating point operations.
"""
generator_elements(sym::PointSymmetry) = element(sym, sym.generators)

"""
    symmetry_name(sym::FiniteTranslationSymmetry)

Name of the point symmetry. Return `FiniteTranslationSymmetry[<Hermann Mauguin>]`.
"""
symmetry_name(sym::PointSymmetry) = "PointSymmetry[$(sym.hermann_mauguin)]"

"""
    symmetry_product(psym::PointSymmetry)

Return a function which takes two point operations and combines them.
Not really necessary, but is defined for consistency with translation symmetry.
"""
function symmetry_product(sym::PointSymmetry)
    function product(lhs::PointOperation, rhs::PointOperation)
        return PointSymmetry(lhs.matrix * rhs.matrix)
    end
    return product
end

function symmetry_canonize(sym::PointSymmetry)
    canonize(arg::TranslationOperation{<:Integer}) = arg
    canonize(arg::PointOperation{<:Integer}) = arg
    canonize(arg::SpaceOperation{<:Integer, <:Integer}) = arg
    return canonize
end

Base.eltype(::PointSymmetry) = PointOperation{Int}
Base.valtype(::PointSymmetry) = PointOperation{Int}
Base.eltype(::Type{<:PointSymmetry}) = PointOperation{Int}
Base.valtype(::Type{<:PointSymmetry}) = PointOperation{Int}

Base.in(item::Any, sym::PointSymmetry) = false
# Base.in(item::IdentityOperation, sym::PointSymmetry) = dimension(item) == dimension(sym)
function Base.in(item::TranslationOperation, sym::PointSymmetry)
    return dimension(item) == dimension(sym) && ispoint(item)
end
Base.in(item::PointOperation{<:Integer}, sym::PointSymmetry) = in(item, elements(sym))
function Base.in(item::SpaceOperation{<:Integer, Tt}, sym::PointSymmetry) where {Tt}
    return ispoint(item) && Base.in(PointOperation(item.matrix), sym)
end

Base.iterate(sym::PointSymmetry) = iterate(elements(sym))
Base.iterate(sym::PointSymmetry, i) = iterate(elements(sym), i)

Base.length(sym::PointSymmetry) = length(elements(sym))

# Irreducible Representations

"""
    character_table(sym::PointSymmetry)

Return the character table of the point symmetry.
"""
character_table(sym::PointSymmetry) = sym.character_table

"""
    irreps(sym::PointSymmetry)

Return the irreducible representations of the point symmetry.
"""
irreps(sym::PointSymmetry) = sym.irreps

"""
    irrep(sym::PointSymmetry, idx::Integer)

Return the `idx`th irreducible representation of the point symmetry.
"""
irrep(sym::PointSymmetry, idx::Integer) = sym.irreps[idx]

"""
    num_irreps(sym::PointSymmetry)

Return the number of irreducible representations of the point symmetry.
Aliases: [`num_irreps(::PointSymmetry)`](@ref),
[`numirreps(::PointSymmetry)`](@ref), [`irrepcount(::PointSymmetry)`](@ref)
"""
num_irreps(sym::PointSymmetry) = length(sym.irreps)

"""
    numirreps(sym::PointSymmetry)

Return the number of irreducible representations of the point symmetry.
Aliases: [`num_irreps(::PointSymmetry)`](@ref),
[`numirreps(::PointSymmetry)`](@ref), [`irrepcount(::PointSymmetry)`](@ref)
"""
numirreps(sym::PointSymmetry) = length(sym.irreps)

"""
    irrepcount(sym::PointSymmetry)

Return the number of irreducible representations of the point symmetry.
Aliases: [`num_irreps(::PointSymmetry)`](@ref),
[`numirreps(::PointSymmetry)`](@ref), [`irrepcount(::PointSymmetry)`](@ref)
"""
irrepcount(sym::PointSymmetry) = length(sym.irreps)

"""
    irrep_dimension(sym::PointSymmetry, idx::Integer)

Dimension of the `idx`th irrep of the point symmetry.
"""
irrep_dimension(sym::PointSymmetry, idx::Integer) = size(first(irrep(sym, idx)), 2)

# """
#     get_site_permutation(lattice, matrix_representation, site_map)
#
# Get a list of `Permutation` which represents the element of the point group
# specified by the `matrix_representation` and `site_map`, which respectively
# contain information about how the Bravais lattice transforms, and how the
# basis sites transforms.
# """
# function get_site_permutation(
#             lattice::Lattice,
#             matrix_representation::AbstractMatrix{<:Integer},
#             site_map::AbstractVector{<:Tuple{<:Integer, <:AbstractVector{<:Integer}}})
#     p = zeros(Int, numsite(lattice.supercell))
#     for (i, (j, dR)) in enumerate(site_map)
#         namei = getsitename(lattice.unitcell, i)
#         namej = getsitename(lattice.unitcell, j)
#         for Ri in lattice.hypercube.coordinates
#             _, Rj = lattice.hypercube.wrap(matrix_representation * Ri + dR)
#             i_super = lattice.supercell.siteindices[(namei, Ri)]
#             j_super = lattice.supercell.siteindices[(namej, Rj)]
#             p[i_super] = j_super
#         end
#     end
#     return Permutation(p)
# end

# function get_site_permutations(lattice::Lattice, point_symmetry::PointSymmetry)
#     # sitemap contains how sites of the "unitcell" transform under
#     # every point group operations. Uses heuristics
#     sitemap = findsitemap(lattice.unitcell, point_symmetry)

#     permutations = Permutation[]
#     sizehint!(permutations, group_order(point_symmetry))

#     for (i_elem, (matrep, orbmap)) in enumerate(zip(point_symmetry.matrix_representations, sitemap))
#         p = zeros(Int, numsite(lattice.supercell))
#         for (i, (j, dR)) in enumerate(orbmap)
#             namei = getsitename(lattice.unitcell, i)
#             namej = getsitename(lattice.unitcell, j)
#             for Ri in lattice.hypercube.coordinates
#                 _, Rj = lattice.hypercube.wrap(matrep * Ri + dR)
#                 i_super = lattice.supercell.siteindices[(namei, Ri)]
#                 j_super = lattice.supercell.siteindices[(namej, Rj)]
#                 p[i_super] = j_super
#             end
#         end
#         push!(permutations, Permutation(p))
#     end
#     return permutations
# end

# function get_irrep_iterator(lattice::Lattice,
#                             tsym::FiniteTranslationSymmetry,
#                             psym::PointSymmetry,
#                             tsym_irrep_index::Integer,
#                             psym_irrep_index::Integer,
#                             tsym_irrep_compo::Integer=1,
#                             psym_irrep_compo::Integer=1,
#                             tol::Real=Base.rtoldefault(Float64)
#                             )
#
#     if tsym_irrep_compo != 1 || psym_irrep_compo != 1
#         @warn "Currently only supports Gamma point, trivial point irrep"
#     end
#     #@assert iscompatible(tsym, psym)
#
#     tsym_permutations = get_site_permutations(lattice, tsym)
#     tsym_irrep = irrep(tsym, tsym_irrep_index)
#     tsym_irrep_components = [m[tsym_irrep_compo, tsym_irrep_compo] for m in tsym_irrep]
#
#     psym_permutations = get_site_permutations(lattice, psym)
#
#     psym_irrep = irrep(psym, psym_irrep_index)
#     psym_irrep_components = [m[psym_irrep_compo, psym_irrep_compo] for m in psym_irrep]
#
#     return (
#         (psym_perm * tsym_perm ,  psym_phase * tsym_phase)
#         for (tsym_perm, tsym_phase) in zip(tsym_permutations,
#                                            tsym_irrep_components)
#         for (psym_perm, psym_phase) in zip(psym_permutations,
#                                            psym_irrep_components)
#         #if abs(tsym_phase)>tol && abs(psym_phase) > tol
#     )
# end

# function parse_seitz(name::AbstractString)
#     m = match(r"([12346m])(([-+]))?((.*)?)?", name)
#     main_str, chiral_str, axis_str = m[1], m[3], m[5]
#     axis = Int[]
#     while !isempty(axis_str)
#         m2 = match(r"^(-?\d)(.*)$", axis_str)
#         push!(axis, parse(Int, m2[1]))
#         axis_str = m2[2]
#     end
#     if !isempty(axis) && first(axis[axis!=0]) < 0
#         axis[:] = -axis
#     end
#     return (main_str, isnothing(chiral_str) ? "" : chiral_str, axis)
# end

# export get_hermann_mauguin
# function get_hermann_mauguin(group::FiniteGroup,
#                               element_names::AbstractVector{<:AbstractString})
#     # 1. find principal
#     is_rotation = [startswith(name, r"-?[2346]") for name in element_names]
#     is_reflection = [startswith(name, r"m") for name in element_names]
#     elements = parse_seitz.(element_names)
#
#     if any(is_rotation)
#         _, principal_rotation_index = minimum((-pl, i) for (i, pl) in enumerate(group.period_lengths) if is_rotation[i])
#
#
#     else
#         if any(is_reflection)
#             conjugacy_elements = [x[1] for x in elements[ [first(cc) for cc in group.conjugacy_classes[2:end]] ]]
#             return join(conjugacy_elements)
#         else
#             return "1"
#         end
#     end
# end


text:
using Clang.cindex
using Test

@testset "function" begin
    top = cindex.parse_header(joinpath(@__DIR__, "cxx/cxxbasic.h");
                              cplusplus = true)

    funcs = cindex.search(top, "func")
    @test length(funcs) == 1
    f = funcs[1]

    #@test map(spelling, cindex.function_args(f)) == ["Int", "Int"]
    @test map(spelling, cindex.function_args(f)) == ["x", "y"]
    @test isa(return_type(f), IntType)

    defs = cindex.function_arg_defaults(f)
    @test defs == (nothing, -10)

    funcs = cindex.search(top, "func2")
    @test length(funcs) == 1
    f = funcs[1]

    @test all(map(x->x.text, function_return_modifiers(f)) .== ["const", ])
end

ext: jl


text:
ordered_keys = [:format, :style, :color, :cellcolor]

function apply(row::CellArray)

    string(join(apply.(row), " & "), " \\\\ \n")
end

function apply(M::CellMatrix)
    rows = map(enumerate(eachrow(M))) do (i, row)
        string(join(apply.(row), " & "), " \\\\ \n")
    end
    return reduce(string, rows)
end

apply(c::Cell) = apply(c.style, c.val)

function apply(styles::Dict{Symbol, LatexStyle}, x)
    val = x
    for key in ordered_keys
        haskey(styles, key) || continue
        val = apply(styles[key], val)
    end
    return val
end

# Color
mutable struct Color <: LatexStyle
    val::Symbol
end

apply(s::Color, x) = string("\\color{", s.val,"}{", x, "}")

# Style
mutable struct Style <: LatexStyle
    val::Symbol
end

function apply(s::Style, x)
    stl = s.val
    if !(stl in [:normal, :bold, :italic])
        throw(ArgumentError("Style $(stl) not supported. Use one of the following [:normal, :bold, :italic]"))
    end

    if stl == :normal
        return x
    elseif stl == :bold
        l_style = "\\textbf"
    else
        l_style = "\\textit"
    end

    i1 = findfirst('\$', x)
    i2 = findlast('\$', x)
    if !isnothing(i1) && !isnothing(i2)
        return x[1:i1] * l_style * "{" * x[i1+1:i2-1] * "}" * x[i2:end]
    else
        return l_style * "{" *  x * "}"
    end
end

# CellColor
mutable struct CellColor <: LatexStyle
    val::Symbol
end

apply(s::CellColor, x) = string("\\cellcolor{", s.val,"}", x)

# Format
mutable struct Format <: LatexStyle
    type::Char
    precision::Int
end

Format(s::Char) = Format(string(s))

function Format(s::String)

    s_len = length(s)
    if s_len > 2
        throw(ArgumentError("Format $s too long. It should be \"3f\""))
    elseif s_len == 2
        if !(s[2] in ['d', 'e', 'f', 'p', 's'])
            throw(ArgumentError("Format $s has wrong second char. Use one of the following ['d', 'e', 'f', 'p', 's']"))
        else
            return Format(s[2], parse(Int, s[1]))
        end
    elseif s_len == 1
        if !(s[1] in ['d', 's'])
            throw(ArgumentError("Format $s has one char only. Use one of the following ['d', 's'] or use two chars for remaining."))
        else
            return Format(s[1], 0)
        end
    end
end

string_dig(x, digits::Int) = format(x, precision=digits)

apply(s::Format, x) = x

function apply(s::Format, x::Real)

    prec = s.precision
    type = s.type
    if type == 'f'
        r = string_dig(x, prec)
    elseif type == 'p'
        r = string_dig(100*x, prec) * "\\%"
    elseif type == 'd'
        r = string_dig(x, 0)
    elseif type == 'e'
        if x == 0
            r = "0"
        else
            order_of_mag = floor(log10(abs(x)))
            value_print  = x*10^(-order_of_mag)
            r = string_dig(value_print, prec) * "\\cdot 10^{" * string_dig(order_of_mag, 0) * "}"
        end
    elseif type == 's'
        return string(x)
    else
        error("Type not defined.")
    end

    return "\$" * r * "\$"
end


text:
import Term: Dendogram, Panel, Tree, inspect, theme, expressiontree
using Term.dendogram

# -------------------------------- playground -------------------------------- #
# WORKS
# e = :(2x + 3y + 2)
# e = :(2x + 3y + 2z)
# e = :(2x + 3 + 2 + 2y)
# e = :(2x^(3+y))
# e = :(x^(3+y))
# e = :(1 + 1 - 2x^2)
# e = :(mod(22, 6))
# e = :(2x^(3+y) + 2z)
# e = :(mod(x, 2) + sum(map(y->2y, 1:10)))
# e = :(2x + 3 * √(3x^2))
# e = :(print(x))
# e = :(print(lstrip("test")))

# NOWORKS
# e = :(for i in 1:10; println(i); end)

e = :(2x + 3y + 2z^2)

theme.tree_max_width = 100

print("\n"^3)
expressiontree(e)
inspect(e)

# TODO make Tree work for vector of Pair/Dict


text:
using TravelerProblem;
using Test;

@testset "TravelerProblem.jl" begin

        @testset "LP vs Bruteforce" begin
                #Compiling
                tm = TravelerProblem.Test.trcmatrix(4)[1];
                TravelerProblem.solve_with_linear_programming(tm);
                TravelerProblem.solve_with_very_brute_force(tm);
                println();
                printstyled("Testing LP vs Bruteforce",color = :green, bold = true);
                for i in 5:8
                        cm,rr = TravelerProblem.Test.trcmatrix(i);
                        println();
                        printstyled("Matrix dimension ($i,$i)\n", color = :blue)
                        println("Solving with linear programming:")
                        @time lpr = TravelerProblem.solve_with_linear_programming(cm)
                        println("Solving with very brute force:")
                        @time vbf = TravelerProblem.solve_with_very_brute_force(cm);
                        @test lpr.rute == rr || lpr.rute == TravelerProblem.reverserute(rr);
                        @test vbf.rute == rr || vbf.rute == TravelerProblem.reverserute(rr);
                end
                println();
        end

        @testset "LP with bigger matrices" begin

                println();
                printstyled("Testing LP for bigger matrices!!!",color = :green, bold = true);
                for i in 10:10:150
                        cm,rr = TravelerProblem.Test.trcmatrix(i, Int32);
                        println();
                        printstyled("Matrix dimension ($i,$i)\n", color = :blue)
                        println("Solving with linear programming:")
                        @time lpr = TravelerProblem.solve_with_linear_programming(cm)
                        @test lpr.rute == rr || lpr.rute == TravelerProblem.reverserute(rr);
                end
                println();

        end

end


text:
# This is modified from ForwardDiff in order to allow
# partials to be a tuple of heterogeneous types.
# It is basically unrolling a loop over all partial derivatives
# and performing the basic operations that required

# Switch to the right generated function depending on type
function add_partials(a::T1, b::T2)  where {T1, T2}
    if isnamed(T1)
        add_named_partials(a,b)
    else
        add_unnamed_partials(a,b)
    end
end
function sub_partials(a::T1, b::T2)  where {T1, T2}
    if isnamed(T1)
        sub_named_partials(a,b)
    else
        sub_unnamed_partials(a,b)
    end
end
function neg_partials(a::T1)  where T1
    if isnamed(T1)
        neg_named_partials(a)
    else
        neg_unnamed_partials(a)
    end
end
function mul_partials(a::T1, b::T2, afactor, bfactor) where {T1, T2}
    if isnamed(T1)
        mul_named_partials(a,b, afactor, bfactor)
    else
        mul_unnamed_partials(a,b, afactor, bfactor)
    end
end
function scale_partials(a::T1, x)  where T1
    if isnamed(T1)
        scale_named_partials(a, x)
    else
        scale_unnamed_partials(a, x)
    end
end

# Operations on partial derivatives that are not named
function partialexpr(f, N)
    ex = Expr(:tuple, [f(i) for i=1:N]...)
    return quote
        $(Expr(:meta, :inline))
        @inbounds return $ex
    end
end

@generated function add_unnamed_partials(a::UnitPartials{N,T1}, b::UnitPartials{N,T2})  where {N, T1 <: Tuple, T2 <: Tuple}
    return partialexpr(i -> :(a.partials[$i] + b.partials[$i]), N)
end
@generated function sub_unnamed_partials(a::UnitPartials{N,T1}, b::UnitPartials{N,T2})  where {N, T1 <: Tuple, T2 <: Tuple}
    return partialexpr(i -> :(a.partials[$i] - b.partials[$i]), N)
end
@generated function neg_unnamed_partials(a::UnitPartials{N,T1})  where {N, T1 <: Tuple, T2 <: Tuple}
    return partialexpr(i -> :(-a.partials[$i]), N)
end
@generated function mul_unnamed_partials(a::UnitPartials{N,T1}, b::UnitPartials{N,T2}, afactor, bfactor)  where {N, T1 <: Tuple, T2 <: Tuple}
    return partialexpr(i -> :( (afactor*a.partials[$i]) + (bfactor*b.partials[$i])), N)
end
@generated function scale_unnamed_partials(a::UnitPartials{N,T1}, x)  where {N, T1 <: Tuple, T2 <: Tuple}
    return partialexpr(i -> :( x*a.partials[$i] ), N)
end

# Operations on partial derivatives that are named
function namedpartialexpr(f, ::Val{names}) where names
    ex = :(UnitPartials(;$((f(i) for i in names)...)))
    return quote
        $(Expr(:meta, :inline))
        @inbounds return $ex
    end
end

@generated function add_named_partials(a::UnitPartials{N,T1}, b::UnitPartials{N,T2})  where {N, T1 <: NamedTuple, T2 <: NamedTuple}
    return namedpartialexpr(i -> :($i = a.partials.$i + b.partials.$i), Val(T1.parameters[1]))
end
@generated function sub_named_partials(a::UnitPartials{N,T1}, b::UnitPartials{N,T2})  where {N, T1 <: NamedTuple, T2 <: NamedTuple}
    return namedpartialexpr(i -> :($i = a.partials.$i - b.partials.$i), Val(T1.parameters[1]))
end
@generated function neg_named_partials(a::UnitPartials{N,T1})  where {N, T1 <: NamedTuple, T2 <: NamedTuple}
    return namedpartialexpr(i -> :($i = -a.partials.$i), Val(T1.parameters[1]))
end
@generated function mul_named_partials(a::UnitPartials{N,T1}, b::UnitPartials{N,T2}, afactor, bfactor)  where {N, T1 <: NamedTuple, T2 <: NamedTuple}
    return namedpartialexpr(i -> :($i = (afactor*a.partials.$i) + (bfactor*b.partials.$i)), Val(T1.parameters[1]))
end
@generated function scale_named_partials(a::UnitPartials{N,T1}, x)  where {N, T1 <: NamedTuple, T2 <: NamedTuple}
    return namedpartialexpr(i -> :($i = x*a.partials.$i), Val(T1.parameters[1]))
end


text:
function ObservableSatellite(arg0::jint)
    return ObservableSatellite((jint,), arg0)
end

function get_clock_offset_driver(obj::ObservableSatellite)
    return jcall(obj, "getClockOffsetDriver", ParameterDriver, ())
end

function get_propagator_index(obj::ObservableSatellite)
    return jcall(obj, "getPropagatorIndex", jint, ())
end


text:
using JET, Documenter, Literate, Markdown

const DOC_SRC_DIR          = normpath(@__DIR__, "src")
const INDEX_FILENAME       = normpath(DOC_SRC_DIR, "index.md")
const PLUGIN_API_FILENAME  = normpath(DOC_SRC_DIR, "generated-plugin-api.md")
const PLUGIN_EXAMPLES_DIRS = (normpath(@__DIR__, "..", "examples"), normpath(DOC_SRC_DIR, "generated-plugin-examples"))

function generate_index!()
    isfile(INDEX_FILENAME) && rm(INDEX_FILENAME)
    open(INDEX_FILENAME, write=true) do io
        s = string(@doc JET)
        write(io, s)
    end

    return relpath(INDEX_FILENAME, DOC_SRC_DIR)
end

function generate_example_docs!(dir = PLUGIN_EXAMPLES_DIRS[1], outs = String[])
    # clean up first
    outdir = PLUGIN_EXAMPLES_DIRS[2]
    isdir(outdir) && rm(outdir; recursive = true)

    for (root, dirs, files) in walkdir(dir)
        for file in files
            endswith(file, ".jl") || continue
            push!(outs, Literate.markdown(normpath(root, file), outdir; documenter=true))
        end
        for dir in dirs
            gen_example_doc!(normpath(root, dir), outs)
        end
    end

    return relpath.(outs, DOC_SRC_DIR)
end

function codeblock(s, header = "@docs")
    return """
    ```$header
    $s
    ```
    """ |> Markdown.parse
end

function generate_api_doc(examples_pages)
    out = relpath(PLUGIN_API_FILENAME, DOC_SRC_DIR)

    isfile(PLUGIN_API_FILENAME) && rm(PLUGIN_API_FILENAME)
    open(PLUGIN_API_FILENAME, write=true) do io
        contents = codeblock("Pages = $(repr([out]))", "@contents")
        interface_docs = codeblock(join(JET.JETInterfaces.DOCUMENTED_NAMES, '\n'))
        examples_contents = codeblock("Pages = $(repr(examples_pages))", "@contents")

        s = md"""
        # [JET Analyzer Framework](@id JET-Analyzer-Framework)

        $contents

        JET offers an infrastructure to implement a "plugin" code analyzer.
        Actually, JET's default error analyzer that explained [in the usage section](@ref usages)
        is one specific instance of such a pluggin analyzer built on top of the framework.

        In this documentation we will try to elaborate the framework APIs and showcase example analyzers.

        !!! warning
            The APIs described in this page is _very_ experimental and subject to changes.
            And this documentation is also very WIP.

        ```@meta
        CurrentModule = JET
        ```

        ## Interfaces

        $interface_docs

        ## Examples

        $examples_contents

        """ |> string

        write(io, s)
    end
    @assert isfile(PLUGIN_API_FILENAME)

    return out
end

let
    examples = generate_example_docs!()
    makedocs(; modules = [JET],
               sitename="JET.jl",
               pages = [
                    "README" => generate_index!(),
                    "Usages" => "usages.md",
                    "Configurations" => "config.md",
                    "Internals" => "internals.md",
                    "JET Analyzer Framework" => Any[
                        "API"      => generate_api_doc(examples),
                        "Examples" => examples,
                    ]
               ],
               format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
               )
end

deploydocs(; repo = "github.com/aviatesk/JET.jl.git",
             push_preview = true,
             )


text:
using Test
using Dates
using ExtensibleScheduler
using ExtensibleScheduler: InstantTrigger, get_next_dt_fire
using Base: IteratorSize, HasLength, length

@testset "InstantTrigger" begin

    @testset "trigger" begin
        dt_fire_at, dt_previous_fire, dt_now, dt_next_fire = DateTime(2009, 7, 6), DateTime(0), DateTime(2008, 5, 4), DateTime(2009, 7, 6)

        trigger = Trigger(dt_fire_at)

        @test trigger isa InstantTrigger
        @test get_next_dt_fire(trigger, dt_previous_fire, dt_now) == dt_next_fire

        @test IteratorSize(trigger) == HasLength()
        @test length(trigger) == 1
    end

    @testset "sample" begin
        n_executed = 0

        function sample_action()
            n_executed += 1
        end

        dt_start = DateTime(2010, 1, 1, 0, 0)
        dt = dt_start
        clock = SimClock(dt)
        sched = BlockingScheduler(clock=clock, delayfunc=NoSleep)
        action = Action(sample_action)
        add(sched, action, Trigger(DateTime(2010, 1, 1, 0, 3)))
        add(sched, action, Trigger(DateTime(2010, 1, 1, 0, 5)))
        run_pending(sched)
        @test n_executed == 0

        dt = dt + Dates.Minute(1)  # 00:01
        println("Clock set to $dt")
        set(clock, dt)
        run_pending(sched)
        @test n_executed == 0

        dt = dt + Dates.Minute(1)  # 00:02
        println("Clock set to $dt")
        set(clock, dt)
        run_pending(sched)
        @test n_executed == 0

        dt = dt + Dates.Minute(1)  # 00:03
        println("Clock set to $dt")
        set(clock, dt)
        run_pending(sched)
        @test n_executed == 1

        dt = dt + Dates.Minute(1)  # 00:04
        println("Clock set to $dt")
        set(clock, dt)
        run_pending(sched)
        @test n_executed == 1

        dt = dt + Dates.Minute(1)  # 00:05
        println("Clock set to $dt")
        set(clock, dt)
        run_pending(sched)
        @test n_executed == 2

    end
end


text:
#module GraphicReasoning

using Calculus

#puts a semiCircle on the graph starting at x=a,y=fa, and going until x=b.
#integral is (fa+pi/2)*(b-a)
function semiCircle(a::Int,fa,b::Int)
	r = (b-a)/2.0
	x0 = (b+a)/2.0

	flip = rand([true false])
	if flip
		x -> fa - sqrt((r)^2-(x-(x0))^2)
	else
		x -> fa + sqrt((r)^2-(x-(x0))^2)
	end
end

function makeLineEnd(exc)
	y = range(-5,5,length=11)
	let cond = true
		while(cond)
		fb = Int(rand(y))
		if fb != exc
			return fb
		end
		end
	end
end

#integral is line(b)*(b-a)
function line(a::Int,fa,b::Int)
	y = range(-5,5,length=11)
	cond = true
	fb = makeLineEnd(fa)

	m = (fb-fa)/(b-a)

	x -> m*(x-a)+fa
end

function plotQuestion()
	p = plot()

	d = [2,3,4]
	y_start = 0
	pos = -5
	ys = []
	while pos<=5
		ds = rand(d)
		is_line = rand([true,false])
			if is_line
				f = line(pos,y_start,pos+ds)
				y_start = Int(f(pos+ds))
				txt = "($(pos+ds),$(y_start))"
				plot!(f,pos,pos+ds;annotation=(pos+ds,y_start,text(txt,8)))
				push!(ys,y_start)
			else
				f = semiCircle(pos,y_start,pos+ds)
				plot!(f,pos,pos+ds)
				annotate!([(pos+ds,y_start,Plots.text("($(pos+ds),$(y_start))",:left,8))])
				annotate!([(pos+(ds/2),f(pos+(ds/2)),Plots.text("($pos,$(f(pos+(ds/2))))",:bottom,8))])
				push!(ys,f(pos+(ds/2)))
			end
			if pos+ds>5
				push!(ys,pos+ds)
				break
			end
			pos = pos+ds
	end
	#going through x=-5 to x=5
	y0 = min(minimum(ys),-5)
	yf = max(maximum(ys),5)

	plot!(
	xlims=[y0,yf],
	ylims=[y0,yf],
	title="graph of f",
	minorgrid=true,
	size=(600,600),
	legend=false
	)

	savefig("question.png")

end

#end


text:
module ExtInt

push!(LOAD_PATH, ".")

using Error
using Lexer
export parse, calc, NumVal, ClosureVal

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

abstract OWL
abstract Environment
abstract RetVal
type Num <: OWL
	n::Real
end
type Binop <: OWL
	op::Function
	lhs::OWL
	rhs::OWL
end
type Unop <: OWL
	op::Function
	operand::OWL
end
type If0 <: OWL
     	condition::OWL
	zero_branch::OWL
	nonzero_branch::OWL
end
type Binder <: OWL
     	name::Symbol
	binding_expr::OWL
end
type With <: OWL
     	binders::Array{Binder}
	body::OWL
end
type Id <: OWL
     	name::Symbol
end
type FunDef <: OWL
     	formal_parameters::Array{Symbol}
	fun_body::OWL
end
type FunApp <: OWL
     	fun_expr::OWL
	arg_exprs::Array{OWL}
end
type NumVal <: RetVal
     	n::Real
end
type ClosureVal <: RetVal
     	params::Array{Symbol}
	body::OWL
	env::Environment  # this is the environment at definition time
end
type mtEnv <: Environment
end
type SymVal <: OWL
     	name::Symbol
	value::RetVal
end
type CEnvironment <: Environment
     	symvals::Array{SymVal}
	parent::Environment
end

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

function collatz(n::Real)
	return collatz_helper(n, 0)
end
function collatz_helper(n::Real, num_iters::Int)
	if n == 1
	      	return num_iters
	end
	if mod(n,2)==0
		return collatz_helper(n/2, num_iters+1)
	else
		return collatz_helper(3*n+1, num_iters+1)
	end
end

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

# All operators will be stored in this dicitonary
op_table = Dict([	(:+,+),
			(:-,-),
			(:*,*),
			(:/,/),
			(:mod,mod),
			(:collatz,collatz)
			])
reserved_words=[:+,:-,:*,:/,:mod,:collatz,:if0,:with,:lambda]

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

# This is the 'base case' of the parse tree.
function parse(expr::Real)
	return Num(expr)
end
# Parse things with type Id
function parse(expr::Symbol)
	if !(expr in reserved_words) # Make sure the word is not reserved
		return Id(expr)
	else
	   	throw(LispError("Symbol \'$expr\' is a reserved word"))
	end
end
# When an array is encountered use this
function parse(expr::Array{Any})
	# Check to make sure array is not empty
	if length(expr) >= 1
	   	# The operator should be the first argument
		op_symbol = expr[1]
		if haskey(op_table, op_symbol) # Check to make sure this is a valid op.
		   	if length(expr) == 3 # All Operators have 3 or 2 length => 2 or 1 arguments respectively.
			   	if op_symbol != :collatz # Only the collatz operator needs 1 argument.
				   	return Binop(op_table[op_symbol], parse(expr[2]), parse(expr[3]))
				else
					throw(LispError("Too many parameters given for \'$op_symbol\'")) # Throws if 'collatz' is called with more than one argument.
				end
			elseif length(expr) == 2
			       	if op_symbol == :- || op_symbol == :collatz # These are the only two allowed 1-arg symbols
				   	return Unop(op_table[op_symbol], parse(expr[2]))
				else
					throw(LispError("Not enough parameters given for \'$op_symbol\'")) # Throw this if not uniary operator
				end
			else
				throw(LispError("Wrong amount of parameters given for \'$op_symbol\'")) # Throw this if not correct arg length
			end
			# Handle non-binary/uniary operators and AST elements seperately
		elseif op_symbol==:if0
		        if length(expr)==4 # Verify correct input length
				return If0(parse(expr[2]), parse(expr[3]), parse(expr[4])) # Good Input gives good return
				# In case it isn't clear, 'If0' is a node in the AST and has a type defined for it.
			else
				throw(LispError("Wrong amount of parameters given for \'If0\'"))
			end
		elseif op_symbol == :with
	       	        if length(expr) == 3 # Verify correct argument length
		   	   	binding_exprs=expr[2] # This should be an array of 'pairs'
				if no_dups(binding_exprs) # Check to make sure there's no duplicates in the binders
					binders=Binder[]
					for i in 1:size(binding_exprs,1) # Go through list and bind values
				      	       push!(binders, Binder(binding_exprs[i][1], parse(binding_exprs[i][2]))) # Evaluate/parse all bindings
					end
					return With(binders, parse(expr[3])) # The third element is the actual body
					# We'll handle the 'With' node in the AST in the calc statement
				else
					throw(LispError("Duplicate symbols in \'with\' statement"))
				end
			else
				throw(LispError("Wrong amount of parameters given for \'with\'"))
			end
		elseif op_symbol == :lambda
	       	        if length(expr) == 3 # Verify correct argument length
		   	   	if no_dups(expr[2]) # Make sure the arguments aren't duplicated
					return FunDef(expr[2],parse(expr[3])) # elements: 2- argument list. 3-function body
				else
					throw(LispError("Duplicate symbols in \'lambda\' expression"))
				end
	       		else
				throw(LispError("Wrong amount of parameters given for \'lambda\'"))
			end
		else # Everything else must evaluate to a function
			arg_exprs=OWL[]
			for i in 2:size(expr,1)
		      	      	push!(arg_exprs, parse(expr[i]))
			end
			return FunApp(parse(op_symbol), arg_exprs) # The op symbol isn't a symbol. It's a function expression
		end
	else
		throw(LispError("No parameters given: Empty array."));
	end
end
# Catch all other parse types.
function parse(expr::Any)
	throw(LispError("Invalid type \'$expr\'"))
end
# This is the helper function that checks for duplicate symbols in with and lambda statements
function no_dups(exprs::Array{Any})
	if length(exprs) > 0
  	   	syms = Any[]
		if typeof(exprs[1]) == Array{Any,1} # With nodes should have this
    	   		 for i in 1:size(exprs,1)
			       	if typeof(exprs[i][1]) != Symbol
				   	throw(LispError("Must use symbols for binding expressions."))
				elseif in(exprs[i][1], reserved_words) # Check for reserved words
				        throw(LispError("Reserved word \'$exprs[i][1]\' cannot be used in binding expression."))
      		      	 	elseif in(exprs[i][1], syms) # Check for duplicates of partial forward-list
        		    	      	return false # Duplicate encountered
				else
					push!(syms, exprs[i][1]) # Push good symbols
			 	end
			end
			return true
		elseif typeof(exprs[1]) == Symbol # Used in lambda nodes
    	       	        for i in 1:size(exprs,1)
      		      	      	if typeof(exprs[i]) != Symbol # Make sure we're using symbols
        		    	   	throw(LispError("Must use symbols for single-dim argument lists."))
				elseif in(exprs[i], reserved_words) # Check for reserved words
				        throw(LispError("Reserved word \'$exprs[i]\' cannot be used in argument lists."))
				elseif in(exprs[i], syms) # Check for duplicates
			    	        return false # Duplicate encountered
			 	else
					push!(syms, exprs[i]) # All is good in the hood.
			 	end
			end
			return true
		else
			throw(LispError("Improper format. Probably missing parentheses."))
		end
	else
    	   	return true # No need to throw error or even check. Empty lists are okay.
  	end
end
# Default behavior to catch incorrect format. Throw error if it's not a list
function no_dups(exprs::Any)
	throw(LispError("Improper format. Array expected."))
end

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

# Calculate the basic OWL/root abstract syntax tree
function calc(ast::OWL)
	 return calc(ast, mtEnv())
end
# Base case is a number or closure which we will return
function calc(owl::Num, env::Environment)
	 return NumVal(owl.n)
end
# Handle binary operations as in the last lab
function calc(owl::Binop, env::Environment) # The environment will just be used on calculating either side.
	 left = calc(owl.lhs, env)
	 right = calc(owl.rhs, env)
	 if typeof(left) == NumVal && typeof(right) == NumVal # Make sure things can evaluated in the calcs
	    	 if (owl.op == /) && right.n == 0 # Make sure that we can't devide by zero
		    	  throw(LispError("Cannot divide by zero in binary op"))
		 else
			  return NumVal(owl.op(left.n, right.n)) # Return correct calculation
		 end
	 else
		 throw(LispError("Invalid binary op types did not evaluate correctly"))
	 end
end
# Handle bnary operations in a similar fashion to binary operations
function calc(owl::Unop, env::Environment) # As in Binop, the environment is just used to calculate the expression recursively
	 unary = calc(owl.operand, env) # Calculate the expression
	 if typeof(unary) == NumVal # Make sure it evaluated to a value
		 if !(owl.op == collatz && unary.n <= 0) # Check to make sure the collatz function is called with a correct input.
			  return NumVal(owl.op(unary.n))
		 else
			  throw(LispError("Cannot perform collatz on \'$unary\'"))
		 end
	 else
		 throw(LispError("Invalid uninary op types did not evaluate correctly"))
	 end
end
# Handle If0 nodes
function calc(owl::If0, env::Environment) # Pass in the environment as always
	 cond = calc(owl.condition, env)
	 if typeof(cond) == NumVal # Make sure that the condition calculated correctly
	    	 if cond.n == 0
		    	  return calc(owl.zero_branch, env) # This is the 'then' branch
		 else
			  return calc(owl.nonzero_branch, env) # This is the 'else' branch
		 end
	 else
		 throw(LispError("Invalid \'If0\' type with value \'$cond\'."))
	 end
end
# Handle Id/Symbols in grammer
function calc(owl::Id, env::Environment)
	 if env == mtEnv() # Check to see if it's an mpty environment.
	    	 throw(LispError("Could not find symbol \'$owl\'"))
	 elseif hasSymVal(env.symvals, owl.name) # See helper function
	 	 return getValue(env.symvals, owl.name) # Return symbol value.
	 else
		 return calc(owl, env.parent) # I believe I need this for recursion
	 end
end
# Handle With nodes
function calc(owl::With, env::Environment)
	 symvals = SymVal[] # The With node has a bunch of sumbols we need to evaluate.
	 for i in 1:size(owl.binders,1)
	       	 push!(symvals, SymVal(owl.binders[i].name, calc(owl.binders[i].binding_expr, env)))
	 end
	 return calc(owl.body, CEnvironment(symvals, env)) # Extend the environment and pass it on.
end
# Define the FunDef Node and pass of the closure
function calc(owl::FunDef, env::Environment)
	 return ClosureVal(owl.formal_parameters, owl.fun_body, env) # super simple I know.
end
# This is the FunApp node
function calc( owl::FunApp, env::Environment )
	 # The function expression should result in a ClosureVal
	 the_closure = calc(owl.fun_expr, env)
	 if typeof(the_closure) == ClosureVal # Check for correct evaluation
	    	 if length(owl.arg_exprs) == length(the_closure.params) # These two lists should be equal.
	 	    	  # Extend the current environment by binding the actual parameters to the formal parameters
	 	 	  symvals = SymVal[] # Build array
	 	 	  for i in 1:size(owl.arg_exprs,1) # Iterate through array and push calcualted symbols
	       	       	      	   push!(symvals, SymVal(the_closure.params[i], calc(owl.arg_exprs[i], env)))
	 	 	  end
	      	 	  return calc(the_closure.body, CEnvironment(symvals, the_closure.env)) # Calculate new environment and return
		 else
		    	  throw(LispError("Number of parameters and arguments does not match in function call/declaration."))
	 	 end
	 else
		 throw(LispError("Invalid type did not return closure"))
	 end
end
# Default calc case
function calc(owl::Any)
	 throw(LispError("Cannot calculate unknown type"))
end
# This could also happen with environments.
function calc(owl::Any, env::Environment)
	 throw(LispError("Cannot calculate unknown type"))
end
# Helper function to check for symbol name in environment array
function hasSymVal(symvals::Array{SymVal}, name::Symbol)
	for i in 1:size(symvals,1) # Iterate through and find
	       	 if symvals[i].name == name
	       	    	 return true # Found
		 end
	end
	return false
	# return (name in symvals) # this won't work with different types fyi
end
# Helper function to get value associated with symbol name
function getValue(symvals::Array{SymVal}, name::Symbol)
	 for i in 1:size(symvals,1)
	       	 if symvals[i].name == name
		    	 return symvals[i].value
		 end
	 end
	 throw(LispError("Symbol with \'$name\' not found - Good job on getting here though"))
end

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

# This is for debugging
function calc(expr::AbstractString)
	 return calc(parse(Lexer.lex(expr)))
end
#
## TODO: check recursive function
#println(calc("(with ( (recur (lambda (x) (
#		    if0 x 0 (recur (- x))
#)) ) ) (recur 5))"))
#

#println(calc("(with ((with 1)) 1)"))

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

end # module

max_stars_repo_name: scottcrossen/JuliaDevelopment

max_issues_repo_name: scottcrossen/JuliaDevelopment

max_forks_repo_name: scottcrossen/JuliaDevelopment


text:
# Autogenerated wrapper script for CGAL_jll for any
using boost_jll
using GMP_jll
using MPFR_jll
JLLWrappers.@generate_wrapper_header("CGAL")
function __init__()
    JLLWrappers.@generate_init_header(boost_jll, GMP_jll, MPFR_jll)
    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
Function to calculate weights for new particles

# Arguments
* `newParticles`: particles from new time step of abcpmc algorithm
* `oldParticles`: particles from previous step of abcpmc algorithm
* `oldWeights::Weights`: weights from previous step of abcpmc algorithm
* `kernelbw`: variance used in transition kernel when drawing new particles
* `kernelDensity::Function`: density of transition kernel, takes a input value,
 old particle, and kernelbw
* `priorDensity::Function`: returns the prior density when given a particle

# Value
weights for new particles
"""
function kernel_weights(newParticles::Array{T, 1},
                        previousStep::A,
                        kernelbw::Float64,
                        density_kernel::Function,
                        density_prior::Function) where {T <: Number, A <: UnivariateAbcPmc}
    if length(newParticles) != length(previousStep.particles)
        error("Length of new particles and old particles must match")
    end

    samplingDensity = zeros(Float64, size(newParticles, 1))
    for (p, w) in zip(previousStep.particles, previousStep.weights.values)
        samplingDensity += (density_kernel.(newParticles, p, kernelbw) .* w)
    end
    priorDensity = density_prior.(newParticles)
    newWeights = priorDensity ./ samplingDensity
    return StatsBase.AnalyticWeights(newWeights ./ sum(newWeights))
end

function kernel_weights(newParticles::Array{T, 2},
                        previousStep::A,
                        kernelbw::Array{Float64, 1},
                        density_kernel::Function,
                        density_prior::Function) where {T <: Number, A <: MultivariateAbcPmc}
    if size(newParticles) != size(previousStep.particles)
        error("Length of new particles and old particles must match")
    end
    n = size(newParticles, 1)
    newWeights = zeros(Float64, n)
    for ii in 1:n
        samplingDensity = 0.0
        for (jj, w) in enumerate(previousStep.weights.values)
            samplingDensity += density_kernel(newParticles[ii, :],
                                              previousStep.particles[jj, :],
                                              kernelbw) * w
        end
        newWeights[ii] = density_prior(newParticles[ii, :]) / samplingDensity
    end
    return StatsBase.AnalyticWeights(newWeights ./ sum(newWeights))
end

"""
    Initialize ABC PMC Algorithm
    """
function pmc_start(summaryStatistics::Any,
                   nparticles::Int64,
                   ninitial::Int64,
                   sample_prior::Function,
                   forward_model::Function,
                   compute_distance::Function,
                   rank_distances::Function,
                   initial_acceptbw::Function,
                   verbose::Bool = true)

    ##Draw first sample to determine typing
    proposal = sample_prior()
    simulatedData = forward_model(proposal)
    dist = compute_distance(simulatedData, summaryStatistics)

    ##Initialize Arrays based on draws
    if typeof(proposal) <: Array
        particles = Array{eltype(proposal)}(ninitial, length(proposal))
    else
        particles = Array{eltype(proposal)}(ninitial)
    end

    if typeof(dist) <: Array
        distances = Array{eltype(dist)}(ninitial, length(dist))
    else
        distances = Array{eltype(dist)}(ninitial)
    end

    particles[1, :] = proposal
    distances[1, :] = dist

    ##Calculate initial samples
    if verbose println("Sampling Initial Particles and Computing Distances...") end
    for ii in 2:ninitial
        proposal = sample_prior()
        simulatedData = forward_model(proposal)
        particles[ii, :] = proposal
        distances[ii, :] = compute_distance(simulatedData, summaryStatistics)
    end

    ##Rank Particles
    if verbose println("Selecting Particles...") end
    ranks = rank_distances(distances)
    ordinals = sortperm(ranks)

    acceptbw = initial_acceptbw(distances, ordinals, nparticles)

    ##Subset Particles
    idxs = ordinals[1:nparticles]
    particles = subset(particles, idxs)
    distances = subset(distances, idxs)
    #=
    if typeof(proposal) <: Array
    particles = particles[idxs, :]
    else
    particles = particles[idxs]
    end

    if typeof(dist) <: Array
    distances = distances[idxs, :]
    acceptbw = vec(maximum(distances, 1))
    else
    distances = distances[idxs]
    acceptbw = maximum(distances)
    end
    =#

    ##Return result
    weights = StatsBase.AnalyticWeights(fill(1.0 / nparticles, nparticles))
    nsampled = fill(floor(Int64, ninitial / nparticles), nparticles)
    return AbcPmc(particles, distances, weights, acceptbw, nsampled)
end

"""
    Function to new step in ABC Population Monte Carlo Algorithm

    # Arguments
    * `previousStep::ABCPMCStep` object containing information on previous step
    * `shrink_acceptbw::Function` rule for computing reduced acceptance acceptbw
    * `transition_kernel::Function` transition kernel for adding noise to previously
     accepted parameter values
    * `forward_model::Function` simulates new data given parameter value
    * `compute_distance::Function` calculates a distance measure between output of
     `forawrd_model` and supplied summary statistics
    * `density_kernel::Function` density value for transition kernel returning
    likelihood for new parameter values given old parameter values
    * `density_prior::Function` evaluates density of parameter value for prior function
    * `verbose::Bool = true` default value is true, should phase of step be reported
    # Value
    ABCPMCStep object containing new posterior distribution, parameters and weights, as well
    as computed distances, total number of parameters sampled, and acceptance acceptbw that
    was used.
    """
function pmc_step(previousStep::A, summaryStatistics::Any,
                  nparticles::Int64, kernel_bandwidth::Function,
                  shrink_acceptbw::Function,
                  sample_kernel::Function,
                  forward_model::Function,
                  compute_distance::Function,
                  density_kernel::Function,
                  density_prior::Function,
                  accept_reject::Function,
                  verbose::Bool = true) where A <: AbcPmcStep

    #Allocate variables
    newParticles = zeros(previousStep.particles)
    newDistances = zeros(previousStep.distances)
    nsampled = zeros(previousStep.nsampled)

    #Compute kernel variance / std
    if verbose println("Calculating Bandwidth...") end
    kernelbw = kernel_bandwidth(previousStep)

    #Compute new acceptbw
    if verbose println("Calculating Acceptbw...") end
    acceptbw = shrink_acceptbw(previousStep)

    #Find new particles
    if verbose println("Sampling Particles...") end
    numAccepted = 0
    while numAccepted < nparticles
        nsampled[numAccepted + 1] += 1

        #Draw Sample
        proposal = StatsBase.sample(previousStep)
        proposal = sample_kernel(proposal, kernelbw)

        #Simulate data
        simulatedData = forward_model(proposal)
        proposalDistance = compute_distance(simulatedData, summaryStatistics)

        #Accept / reject
        #if all(proposalDistance .< acceptbw)
        if accept_reject(proposalDistance, acceptbw)
            numAccepted += 1
            newParticles[numAccepted, :] = proposal
            newDistances[numAccepted, :] = proposalDistance
        end
    end

    if verbose  println("Calculating Weights...") end
    newWeights = kernel_weights(newParticles, previousStep, kernelbw,
                                density_kernel, density_prior)

    #Return result as AbcPmcStep
    return AbcPmc(newParticles, newDistances, newWeights, acceptbw, nsampled)
end

"""
    Formatting for logging run setup
    """
function initlog(filename::String, startTime::DateTime, ninitial::Int64, nparticles::Int64, nsteps::Int64)
    logfile = open(filename, "w")
    write(logfile, "Run Info\n")
    write(logfile, "————————————————————","\n")
    write(logfile, string("Start Time: ", Dates.format(startTime, "Y-mm-dd HH:MM:SS"), "\n"))
    write(logfile, string("Initial Samples: ", ninitial, "\n"))
    write(logfile, string("Particles: ", nparticles, "\n"))
    write(logfile, string("Steps: ", nsteps, "\n"))
    write(logfile, "————————————————————","\n\n")
    close(logfile)
    nothing
end

"""
    Formatting for logging step info
    """
function steplog(filename::String, nstep::Int64, totalTime::Base.Dates.Millisecond, stepTime::Base.Dates.Millisecond, stepacceptbw::G, nsampled::Int64) where G <: AbstractFloat
    logfile = open(filename, "a")
    write(logfile, string("Step ", nstep, " Info\n"))
    write(logfile, "————————————————————","\n")
    write(logfile, string("Total Time: ", duration_to_string(totalTime), "\n"))
    write(logfile, string("Step Time: ", duration_to_string(stepTime), "\n"))
    write(logfile, string("Acceptbw: ", stepacceptbw, "\n"))
    write(logfile, string("Particles Sampled: ", nsampled, "\n"))
    write(logfile, "————————————————————","\n\n")
    close(logfile)
    nothing
end

function steplog(filename::String, nstep::Int64, totalTime::Base.Dates.Millisecond, stepTime::Base.Dates.Millisecond, stepacceptbw::Array{G, 1}, nsampled::Int64) where G <: AbstractFloat
    logfile = open(filename, "a")
    write(logfile, string("Step ", nstep, " Info\n"))
    write(logfile, "————————————————————","\n")
    write(logfile, string("Total Time: ", duration_to_string(totalTime), "\n"))
    write(logfile, string("Step Time: ", duration_to_string(stepTime), "\n"))
    write(logfile, string("Acceptbw: ", stepacceptbw, "\n"))
    write(logfile, string("Particles Sampled: ", nsampled, "\n"))
    write(logfile, "————————————————————","\n\n")
    close(logfile)
    nothing
end

"""
Population Monte Carlo Approximate Bayesian Computation Algorithm

# Arguments
* `summaryStatistics` data for posterior to be conditioned on
* `steps::Int64` number of iterations to run
* `nparticles::Int64` number of posterior samples to return for each iteration
* `ninitial::Int64` number of samples to base initial iterationon
* `sample_prior::Function` generates a random draw from the prior
* `density_prior::Function` pdf of prior
* `compute_distance::Function` compute distance metric between output of
  `forward_model` and provided summary statistics
* `forward_model::Function` function to simulate summary statistics according to new parameter value
* `transition_kernel::Function` transition kernel or "noise" function applied to particles
* `density_kernel::Function` pdf functions for transition kernel
* `rank_distances::Function` method for ranking particles based on computed distances and
  choosing initial selection acceptbw
* `shrink_acceptbw::Function` method for shrinking acceptbw between iterations based on output
  from previous iteration
* `log::Bool = true` should progress be written to a log file
* `logFile::String = "log.txt"` file to write progress to, ignored if `log = false`
* `save::Bool = true` should results be written to a file
* `saveFile::String = "results.jld"` file to write results to
* `verbose::Bool = true` should progress be reported
"""
function abc_pmc(summaryStatistics::Any, nsteps::Int64,
                 nparticles::Int64, ninitial::Int64,
                 sample_prior::Function, density_prior::Function,
                 sample_kernel::Function, density_kernel::Function,
                 forward_model::Function, compute_distance::Function,
                 rank_distances::Function, kernel_bandwidth::Function,
                 shrink_acceptbw::Function;
                 initial_acceptbw::Function = maxdist_bandwidth,
                 accept_reject::Function = boxcar_accept,
                 verbose::Bool = true,
                 log::Bool = true, logFile::String = "log.txt",
                 save::Bool = true, saveFile::String = "results.jld")

    ##Write initial information is running log file
    if log
        startTime = now()
        initlog(logFile, startTime, ninitial, nparticles, nsteps)
    end

    ##Run first step
    results = [pmc_start(summaryStatistics,
                         nparticles, ninitial,
                         sample_prior, forward_model,
                         compute_distance, rank_distances,
                         initial_acceptbw, verbose)]
    ##Save results
    if verbose println("Saving Step...") end
    if save JLD.@save saveFile summaryStatistics results end

    ##Print to log file
    if verbose println(string("Step ", 1, " Complete")) end

    if log
        stepTime = now() - startTime
        totalTime = stepTime
        steplog(logFile, 1, totalTime, stepTime, results[1].acceptbw, ninitial)
    end

    #Subsequent steps
    for ii in 2:nsteps
        push!(results, pmc_step(results[ii - 1], summaryStatistics, nparticles,
                                kernel_bandwidth, shrink_acceptbw, sample_kernel,
                                forward_model, compute_distance, density_kernel,
                                density_prior, accept_reject, verbose))

        ##Save results
        if verbose println("Saving Step...") end
        if save JLD.@save saveFile summaryStatistics results end

        #Print to log file
        if verbose println(string("Step ", ii, " Complete")) end

        if log
            stepTime = now() - startTime - totalTime
            totalTime += stepTime
            steplog(logFile, ii, totalTime, stepTime, results[ii].acceptbw, sum(results[ii].nsampled))
        end

    end
    return results
end

function abc_pmc_warmstart(nsteps::Int64, nparticles::Int64,
                           sample_prior::Function, density_prior::Function,
                           sample_kernel::Function, density_kernel::Function,
                           forward_model::Function, compute_distance::Function,
                           kernel_bandwidth::Function, shrink_acceptbw::Function;
                           accept_reject::Function = boxcar_accept,
                           verbose::Bool = true,
                           log::Bool = true, logFile::String = "log.txt",
                           save::Bool = true, saveFile::String = "results.jld")

    ##Print to log file
    if verbose println("Resuming algorithm...") end

    ##Write initial information is running log file
    if log
        startTime = now()
        totalTime = startTime - now()
    end

    ##Load results
    summaryStatistics = JLD.load(saveFile, "summaryStatistics")
    results = JLD.load(saveFile, "results")

    ##Check Number of steps
    if nsteps <= length(results)
        error("$nsteps alreadty complete")
    end

    #Subsequent steps
    for ii in length(results) + 1:nsteps
        push!(results, pmc_step(results[ii - 1], summaryStatistics, nparticles,
                                kernel_bandwidth, shrink_acceptbw, sample_kernel,
                                forward_model, compute_distance, density_kernel,
                                density_prior, accept_reject, verbose))

        ##Save results
        if verbose println("Saving Step...") end
        if save JLD.@save saveFile results end

        #Print to log file
        if verbose println(string("Step ", ii, " Complete")) end

        if log
            stepTime = now() - startTime - totalTime
            totalTime += stepTime
            steplog(logFile, ii, totalTime, stepTime, results[ii].acceptbw, sum(results[ii].nsampled))
        end

    end
    return results
end


text:
import NNlib
using AutoGrad
using Statistics: mean
using LinearAlgebra: lmul!
export pool, mean

#= NOTES
- maxpool/meanpool vs op=maximum, op=mean (would require Statistics)
- adaptivexxx vs specifying output=1 or window=Inf.
- torch.nn: MaxPoolNd, MaxUnpoolNd, AvgPoolNd, FractionalMaxPoolNd, LPPoolNd, AdaptiveMaxPoolNd, AdaptiveAvgPoolNd
- torch.nn.functional: avg_poolNd, max_poolNd, max_unpoolNd, lp_poolNd, adaptive_max_poolNd, adaptive_avg_poolNd
- adaptive versions take output_size as argument, regular versions take kernel_size.
- torch also supports dilation? cudnn does not.
- torch supports maxunpool but not avgunpool.
- existence of lp may make the op kwarg more preferrable (other ops in the future?)
- keras has MaxPoolingNd, AveragePoolingNd, GlobalMaxPoolingNd, GlobalAveragePoolingNd: only output=1 supported
- keras GlobalPooling includes reshape.
=#

"""
    pool(x; kwargs...)

Compute pooling of input values (i.e., the maximum or average of several adjacent values) to
produce an output with smaller height and/or width.

If `x` has dimensions `(X1,X2,...,Cx,N)`, the result `y` will have dimensions
`(Y1,Y2,...,Cx,N)` where

    Yi=1+floor((Xi+2*padding[i]-window[i])/stride[i])

Here `Cx` is the number of input channels, `N` is the number of instances, and `Xi,Yi` are
spatial dimensions. `window`, `padding`, `stride` are keyword arguments that can be
specified as a single number (in which case they apply to all dimensions), or an array/tuple
with entries for each spatial dimension. If `window[i] > Xi+2*padding[i]` it is truncated to
fit the padded input, e.g. one can use `window=typemax(Int)` to implement global pooling.

# Keywords:

* `op=maximum`: pooling operation, can be `maximum` or `mean`
* `window=2`: the pooling window size for each dimension.
* `padding=0`: the number of extra zeros implicitly concatenated at the start and at the end of each dimension.
* `stride=window`: the number of elements to slide to reach the next pooling window.
* `propagateNaN=false`: if false do not propagate NaN values in max pooling.
* `includePadding=false`: if false do not include padded values in the count for average.
* `channelmajor = false`: assume channel-major format tensors if true
* `alpha=1`: can be used to scale the result.

"""
function pool(x; op=maximum, window=2, padding=0, stride=window, propagateNaN=false, includePadding=false, channelmajor=false, alpha=1)
    checkpoolopts(x, op, window, padding, stride, propagateNaN, includePadding, channelmajor, alpha)
    if channelmajor; x = permutedims(x,(2,3,1,4)); end
    N = ndims(x)
    window = NNlib.expand(Val(N-2), window)
    stride = NNlib.expand(Val(N-2), stride)
    padding = NNlib.expand(Val(N-2), padding)
    window = min.(window, size(x)[1:N-2] .+ 2 .* padding)
    pdims = NNlib.PoolDims(x, window; padding = padding, stride = stride)
    y = (op == maximum ? NNlib.maxpool(x, pdims) : NNlib.meanpool(x, pdims))
    if channelmajor; y = permutedims(y, (3,1,2,4)); end
    alpha == 1 ? y : lmul!(eltype(y)(alpha), y)
end

function ∇pool(x,y,dy; op=maximum, window=2, padding=0, stride=window, propagateNaN=false, includePadding=false, channelmajor=false, alpha=1)
    checkpoolopts(x, op, window, padding, stride, propagateNaN, includePadding, channelmajor, alpha)
    if alpha != 1; y = y ./ eltype(y)(alpha); end
    if channelmajor; x,y,dy = (a->permutedims(a,(2,3,1,4))).((x,y,dy)); end
    N = ndims(x)
    window = NNlib.expand(Val(N-2), window)
    stride = NNlib.expand(Val(N-2), stride)
    padding = NNlib.expand(Val(N-2), padding)
    window = min.(window, size(x)[1:N-2] .+ 2 .* padding)
    pdims = NNlib.PoolDims(x, window; padding = padding, stride = stride)
    dx = (op == maximum ? NNlib.∇maxpool(dy, y, x, pdims) : NNlib.∇meanpool(dy, y, x, pdims))
    if channelmajor; dx = permutedims(dx, (3,1,2,4)); end
    alpha == 1 ? dx : lmul!(eltype(dx)(alpha), dx)
end

@primitive1 pool(x;o...),dy,y  ∇pool(x,y,dy;o...)
@primitive1 ∇pool(x,y,dy;o...),ddx,dx  nothing  nothing  pool(ddx;o...)

function checkpoolopts(x, op, window, padding, stride, propagateNaN, includePadding, channelmajor, alpha)
    @assert op ∈ (maximum, mean) "Pooling op=$op should be one of (maximum,mean)"
    if op == maximum && !propagateNaN
        @warn "propagateNaN=false not yet implemented in NNlib, NaNs will be propagated. See https://github.com/FluxML/NNlib.jl/issues/218" maxlog=1
    end
    if op == mean && any(padding .> 0) && !includePadding
        @warn "includePadding=false not yet implemented in NNlib, pads will be included. See https://github.com/FluxML/NNlib.jl/issues/218" maxlog=1
    end
    if channelmajor
        @warn "channelmajor is not yet implemented in NNlib, using slow permutedims, see NNlib#267" maxlog=1
    end
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#    Tests of row name column.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Show row names" begin
    expected = """
┌───────┬─────┬───────┬───────┬─────┐
│       │  C1 │    C2 │    C3 │  C4 │
│       │ Int │  Bool │ Float │ Hex │
├───────┼─────┼───────┼───────┼─────┤
│ Row 1 │   1 │ false │   1.0 │   1 │
│ Row 2 │   2 │  true │   2.0 │   2 │
│ Row 3 │   3 │ false │   3.0 │   3 │
│ Row 4 │   4 │  true │   4.0 │   4 │
│ Row 5 │   5 │ false │   5.0 │   5 │
│ Row 6 │   6 │  true │   6.0 │   6 │
└───────┴─────┴───────┴───────┴─────┘
"""

    result = pretty_table(
        String, data;
        header = (
            ["C1",  "C2",   "C3",    "C4"],
            ["Int", "Bool", "Float", "Hex"]
        ),
        row_names = ["Row $i" for i in 1:6]
    )

    @test result == expected

    expected = """
┌───────────┬─────┬───────┬───────┬─────┐
│ Row names │  C1 │    C2 │    C3 │  C4 │
│           │ Int │  Bool │ Float │ Hex │
├───────────┼─────┼───────┼───────┼─────┤
│     Row 1 │   1 │ false │   1.0 │   1 │
│     Row 2 │   2 │  true │   2.0 │   2 │
│     Row 3 │   3 │ false │   3.0 │   3 │
│     Row 4 │   4 │  true │   4.0 │   4 │
│     Row 5 │   5 │ false │   5.0 │   5 │
│     Row 6 │   6 │  true │   6.0 │   6 │
└───────────┴─────┴───────┴───────┴─────┘
"""

    result = pretty_table(
        String,
        data;
        header = (
            ["C1",  "C2",   "C3",    "C4"],
            ["Int", "Bool", "Float", "Hex"]
        ),
        row_names = ["Row $i" for i in 1:6],
        row_name_column_title = "Row names"
    )

    @test result == expected

    expected = """
┌───────────┬─────┬───────┬───────┬─────┐
│ Row names │  C1 │    C2 │    C3 │  C4 │
│           │ Int │  Bool │ Float │ Hex │
├───────────┼─────┼───────┼───────┼─────┤
│   Row 1   │   1 │ false │   1.0 │   1 │
│   Row 2   │   2 │  true │   2.0 │   2 │
│   Row 3   │   3 │ false │   3.0 │   3 │
│   Row 4   │   4 │  true │   4.0 │   4 │
│   Row 5   │   5 │ false │   5.0 │   5 │
│   Row 6   │   6 │  true │   6.0 │   6 │
└───────────┴─────┴───────┴───────┴─────┘
"""

    result = pretty_table(
        String, data;
        header = (
            ["C1",  "C2",   "C3",    "C4"],
            ["Int", "Bool", "Float", "Hex"]
        ),
        row_names = ["Row $i" for i in 1:6],
        row_name_column_title = "Row names",
        row_name_alignment = :c
    )

    @test result == expected

    expected = """
┌───────────┬─────┬───────┬───────┬─────
│ Row names │  C1 │    C2 │    C3 │  C ⋯
│           │ Int │  Bool │ Float │ He ⋯
├───────────┼─────┼───────┼───────┼─────
│   Row 1   │   1 │ false │   1.0 │    ⋯
│   Row 2   │   2 │  true │   2.0 │    ⋯
│   Row 3   │   3 │ false │   3.0 │    ⋯
│     ⋮     │  ⋮  │   ⋮   │   ⋮   │  ⋮ ⋱
└───────────┴─────┴───────┴───────┴─────
             1 column and 3 rows omitted
"""

    result = pretty_table(
        String,
        data;
        header = (
            ["C1",  "C2",   "C3",    "C4"],
            ["Int", "Bool", "Float", "Hex"]
        ),
        crop = :both,
        display_size = (12, 40),
        row_names = ["Row $i" for i in 1:6],
        row_name_column_title = "Row names",
        row_name_alignment = :c
    )
    @test result == expected
end


text:
using DataDrivenDiffEq
using ModelingToolkit
using OrdinaryDiffEq
using LinearAlgebra
using Plots
gr()

# Create a test problem
function pendulum(u, p, t)
    x = u[2]
    y = -9.81sin(u[1]) - 0.1u[2]
    return [x;y]
end

u0 = [0.2π; -1.0]
tspan = (0.0, 40.0)
prob = ODEProblem(pendulum, u0, tspan)
sol = solve(prob)

plot(sol)

# Create the differential data
DX = similar(sol[:,:])
for (i, xi) in enumerate(eachcol(sol[:,:]))
    DX[:,i] = pendulum(xi, [], 0.0)
end

# Create a basis
@variables u[1:2]

# Lots of polynomials
polys = [u[1]^0]
for i ∈ 1:3
    for j ∈ 1:3
        push!(polys, u[1]^i*u[2]^j)
    end
end

# And some other stuff
h = [1u[1];1u[2]; cos(u[1]); sin(u[1]); u[1]*u[2]; u[1]*sin(u[2]); u[2]*cos(u[2]); polys...]

basis = Basis(h, u)

# Get the reduced basis via the sparse regression
Ψ = SInDy(sol[:,:], DX, basis, ϵ = 1e-10)

# Transform into ODE System
sys = ODESystem(Ψ)

# Simulate
estimator = ODEProblem(dynamics(Ψ), u0, tspan)
sol_ = solve(estimator, saveat = sol.t)

# Yeah! We got it right
plot(sol, vars = (1,2))
plot!(sol, vars = (1,2))

norm(sol-sol_) # ≈ 1.89e-13


text:
module CUTENSOR

using ..CuArrays
using ..CuArrays: active_context

using CUDAapi

using CUDAdrv
using CUDAdrv: CUstream

using CEnum
const cudaDataType_t = cudaDataType

const libcutensor = if Sys.iswindows()
    # no ccall by soname, we need the filename
    # NOTE: we discover the full path here, while only the wordsize and toolkit versions
    #       would have been enough to construct "cutensor64_10.dll"
    toolkit = find_toolkit()
    path = find_cuda_library("cutensor", toolkit)
    if path === nothing
        nothing
    else
        basename(path)
    end
else
    # ccall by soname; CuArrays.__init__ will have populated Libdl.DL_LOAD_PATH
    "libcutensor"
end

# core library
include("libcutensor_common.jl")
include("error.jl")
include("libcutensor.jl")

# low-level wrappers
include("tensor.jl")
include("wrappers.jl")

# high-level integrations
include("interfaces.jl")

const _handles = Dict{CuContext,cutensorHandle_t}()
const _handle = Ref{cutensorHandle_t}(C_NULL)

function handle()
    if _handle[] == C_NULL
        @assert isassigned(active_context) # some other call should have initialized CUDA
        _handle[] = get!(_handles, active_context[]) do
            context = active_context[]
            handle = cutensorCreate()
            atexit(()->CUDAdrv.isvalid(context))
            handle
        end
    end
    return _handle[]
end

end


text:
using CircGeometry, Test

center = Point(-0.9,0.0)
p1 = Point(-1.0,-1.0)
p2 = Point(1.0,1.0)

fo = CircGeometry.FillingCircle(0.1,center,10.0)
olist = [fo]
outline = CircGeometry.OutlineRectangle(p1,p2)

Fx,Fy = CircGeometry.compute_wall_repulsion(1,olist,outline.wlist)
bool1 = CircGeometry.is_intersecting_walls(outline,fo)
bool2 = CircGeometry.is_inside_outline(fo,outline)

max_stars_repo_name: pseastham/CircGeometry

max_issues_repo_name: pseastham/CircGeometry


text:
using CSyntax.CSwitch
using Test

include("vsSwitch.jl")
@testset "CSwitch" begin
    @testset "fall through" begin
        @test switch_with_fallthrough(1) == Int[1, 2, 3, 4]
        @test switch_with_fallthrough(3) == Int[3, 4]
        @test switch_with_fallthrough(100) == Int[4]
    end
    @testset "leading default" begin
        @test switch_with_leading_default(1) == Int[1, 2, 3]
        @test switch_with_leading_default(3) == Int[3]
        @test switch_with_leading_default(100) == Int[4, 1, 2, 3]
    end
    @testset "with expressions" begin
        @test switch_with_expressions(1) == Int[1, 2, 3, 4]
        @test switch_with_expressions(3) == Int[3, 4]
        @test switch_with_expressions(100) == Int[4]
    end
    @testset "with break" begin
        @test switch_with_break(1) == Int[1]
        @test switch_with_break(3) == Int[3]
        @test switch_with_break(100) == Int[4]
    end
    @testset "enums" begin
        @enum test t=1 f=2
        tester = t
        @cswitch tester begin
            @case t
            x = 1
            break
            @case f
            x = 2
            break
        end
        @test x == 1
    end
end


text:
using Test

include("../src/aes-modes.jl")

# Test cases from FIPS 197
@test gmul(0x57, 0x13) == 0xfe
@test gmul(0x00, 0x00) == 0x00
@test gmul(0x01, 0x01) == 0x01

# Test for property of gmul and gmulinv
for b=0x01:0xFF
	@test gmul(b, gmulinv(b)) == 0x01
end

# Test cases from FIPS 197
const keyunexp1 = [ 0x2b, 0x7e, 0x15, 0x16, 0x28, 0xae, 0xd2, 0xa6, 0xab, 0xf7, 0x15, 0x88, 0x09, 0xcf, 0x4f, 0x3c ]
const keyexp1 = KeyExpansion(keyunexp1, Nks[1], Nrs[1])
const keyexp1expect = "2b7e151628aed2a6abf7158809cf4f3ca0fafe1788542cb123a339392a6c7605f2c295f27a96b9435935807a7359f67f3d80477d4716fe3e1e237e446d7a883bef44a541a8525b7fb671253bdb0bad00d4d1c6f87c839d87caf2b8bc11f915bc6d88a37a110b3efddbf98641ca0093fd4e54f70e5f5fc9f384a64fb24ea6dc4fead27321b58dbad2312bf5607f8d292fac7766f319fadc2128d12941575c006ed014f9a8c9ee2589e13f0cc8b6630ca6"

@test bytes2hex(keyexp1) == keyexp1expect

# Tests from OpenDocument spreadsheet available at
# http://www.nayuki.io/page/aes-cipher-internals-in-excel
const key128 =    "2b7e151628aed2a6abf7158809cf4f3c"
const plain128 =  "3243f6a8885a308d313198a2e0370734"
const cipher128 = "3925841d02dc09fbdc118597196a0b32"

@test AESEncrypt(plain128, key128) == cipher128

# Test vectors from
# http://www.inconteam.com/software-development/41-encryption/55-aes-test-vectors

# AES-128
const key1 =    "2b7e151628aed2a6abf7158809cf4f3c"
const plain1 =  "6bc1bee22e409f96e93d7e117393172a"
const cipher1 = "3ad77bb40d7a3660a89ecaf32466ef97"

# AES-192
const key2 =    "8e73b0f7da0e6452c810f32b809079e562f8ead2522c6b7b"
const plain2 =  "6bc1bee22e409f96e93d7e117393172a"
const cipher2 = "bd334f1d6e45f25ff712a214571fa5cc"

# AES-256
const key3 =    "603deb1015ca71be2b73aef0857d77811f352c073b6108d72d9810a30914dff4"
const plain3 =  "6bc1bee22e409f96e93d7e117393172a"
const cipher3 = "f3eed1bdb5d2a03c064b5a7e3db181f8"

# Encryption tests
@test AESEncrypt(plain1, key1) == cipher1
@test AESEncrypt(plain2, key2) == cipher2
@test AESEncrypt(plain3, key3) == cipher3

# Decryption tests
@test AESDecrypt(cipher1, key1) == plain1
@test AESDecrypt(cipher2, key2) == plain2
@test AESDecrypt(cipher3, key3) == plain3

# AES ECB
const key4 =    "2b7e151628aed2a6abf7158809cf4f3c"
const plain4 =  "6bc1bee22e409f96e93d7e117393172a"
const cipher4 = "3ad77bb40d7a3660a89ecaf32466ef97"

@test AESECB(plain4, key4, true) == cipher4
@test AESECB(cipher4, key4, false) == plain4

# AES CBC
const iv5 =     "000102030405060708090a0b0c0d0e0f"
const key5 =    key4
const plain5 =  plain4
const cipher5 = "7649abac8119b246cee98e9b12e9197d"

@test AESCBC(plain5, key5, iv5, true) == cipher5
@test AESCBC(cipher5, key5, iv5, false) == plain5

# AES CFB
const iv6 =     iv5
const key6 =    key4
const plain6 =  plain4
const cipher6 = "3b3fd92eb72dad20333449f8e83cfb4a"

@test AESCFB(plain6, key6, iv6, true) == cipher6
@test AESCFB(cipher6, key6, iv6, false) == plain6

# AES OFB
const iv7 =     iv5
const key7 =    key4
const plain7 =  plain4
const cipher7 = cipher6

@test AESOFB(plain7, key7, iv7) == cipher7
@test AESOFB(cipher7, key7, iv7) == plain7

# AES CTR
# Test with four blocks for AES to test counter incrementation
# Source: https://github.com/servo/nss/blob/3d07e85a597ce9c5c1ec80f85983efcb26aa58e1/cmd/bltest/tests/aes_ctr/aes_ctr_tests_source.txt#L27-L54
const iv8 =     "f0f1f2f3f4f5f6f7f8f9fafbfcfdfeff" # first input block
#               "f0f1f2f3f4f5f6f7f8f9fafbfcfdff00" # second input block
#               "f0f1f2f3f4f5f6f7f8f9fafbfcfdff01" # third input block
#               "f0f1f2f3f4f5f6f7f8f9fafbfcfdff02" # fourth input block
const key8 =    key4
const plain8 =  string(plain4, "ae2d8a571e03ac9c9eb76fac45af8e51", "30c81c46a35ce411e5fbc1191a0a52ef", "f69f2445df4f9b17ad2b417be66c3710")
const cipher8 = string("874d6191b620e3261bef6864990db6ce", "9806f66b7970fdff8617187bb9fffdff", "5ae4df3edbd5d35e5b4f09020db03eab", "1e031dda2fbe03d1792170a0f3009cee")

@test AESCTR(plain8, key8, iv8) == cipher8
@test AESCTR(cipher8, key8, iv8) == plain8

# Test overflow of 8-byte big endian counter
const iv9 =     "aabbccddeeffaabbffffffffffffffff" # first input block
#               "aabbccddeeffaabb0000000000000000" # second input block
#               "aabbccddeeffaabb0000000000000001" # third input block
#               "aabbccddeeffaabb0000000000000002" # fourth input block
const key9 =    key8
const plain9 =  plain8
const cipher9 = string("036c5f64f290c0a3efed8bb8fc3a98a5", "40083d03ba096ab72543aa92c590817c", "d1ea82fc1aa999725dd4c4e627703a46", "375d1e844f3cf18b57f230d2def50441")

@test AESCTR(plain9, key9, iv9) == cipher9
@test AESCTR(cipher9, key9, iv9) == plain9

# Encrypt four random blocks using different modes of operation and
# check if decryption recovers original blocks
const BLOCK_BYTES = 16
const ivrand =     rand(UInt8, BLOCK_BYTES)
const keysrand =   rand(UInt8, div(128, 8))
const keymrand =   rand(UInt8, div(192, 8))
const keylrand =   rand(UInt8, div(256, 8))
const plainrand =  rand(UInt8, 4 * BLOCK_BYTES)
const plainrandl = rand(UInt8, 4 * BLOCK_BYTES + 1)

# AES ECB
for key in (keysrand, keymrand, keylrand)
	cipherrand = AESECB(plainrand, key, true)
	@test cipherrand != plainrand
	@test AESECB(cipherrand, key, false) == plainrand
end

# AES CBC
for key in (keysrand, keymrand, keylrand)
	cipherrand = AESCBC(plainrand, key, ivrand, true)
	@test cipherrand != plainrand
	@test AESCBC(cipherrand, key, ivrand, false) == plainrand
end

# AES CFB
for key in (keysrand, keymrand, keylrand)
	for plain in (plainrand, plainrandl)
		cipherrand = AESCFB(plain, key, ivrand, true)
		@test cipherrand != plain
		@test AESCFB(cipherrand, key, ivrand, false) == plain
	end
end

# AES OFB
for key in (keysrand, keymrand, keylrand)
	for plain in (plainrand, plainrandl)
		cipherrand = AESOFB(plain, key, ivrand)
		@test cipherrand != plain
		@test AESOFB(cipherrand, key, ivrand) == plain
	end
end

# AES CTR
for key in (keysrand, keymrand, keylrand)
	for plain in (plainrand, plainrandl)
		cipherrand = AESCTR(plain, key, ivrand)
		@test cipherrand != plain
		@test AESCTR(cipherrand, key, ivrand) == plain
	end
end


text:
##################################################################################
# This file is part of ModelBaseEcon.jl
# BSD 3-Clause License
# Copyright (c) 2020, Bank of Canada
# All rights reserved.
##################################################################################

export Transformation, NoTransform, LogTransform, NegLogTransform, transformation, inverse_transformation

"""
    transformation(::Type{<:Transformation})

Return a `Function` that will be substituted into the model equations and will be
called to transform the input data before solving. See also
[`inverse_transformation`](@ref).

It is expected that `transformation(T) ∘ inverse_transformation(T) == identity`
and `inverse_transformation(T) ∘ transformation(T) == identity`, but these is
not verified.

"""
function transformation end

"""
    inverse_transformation(::Type{<:Transformation})

Return a `Function` that will be called to transform the simulation data after solving. See also
[`transformation`](@ref).

It is expected that `transformation(T) ∘ inverse_transformation(T) == identity`
and `inverse_transformation(T) ∘ transformation(T) == identity`, but these is
not verified.

"""
function inverse_transformation end

"""
    abstract type Transformation end

The base class for all variable transformations.
"""
abstract type Transformation end
transformation(T::Type{<:Transformation}) = error("Transformation of type $T is not defined.")
inverse_transformation(T::Type{<:Transformation}) = error("Inverse transformation of type $T is not defined.")

"""
    NoTransform <: Transformation

The identity transformation.
"""
struct NoTransform <: Transformation end
transformation(::Type{NoTransform}) = Base.identity
inverse_transformation(::Type{NoTransform}) = Base.identity

"""
    LogTransform <: Transformation

The `log` transformation. The inverse is of course `exp`. This is the default
for variables declared with `@log`.
"""
struct LogTransform <: Transformation end
transformation(::Type{LogTransform}) = Base.log
inverse_transformation(::Type{LogTransform}) = Base.exp

"""
    NegLogTransform <: Transformation

The `log(-x)`, with the inverse being `-exp(x)`. Use this when the variable is
negative with exponential behaviour (toward -∞).
"""
struct NegLogTransform <: Transformation end
"logm(x) = log(-x)" @inline logm(x) = log(-x)
"mexp(x) = -exp(x)" @inline mexp(x) = -exp(x)
transformation(::Type{NegLogTransform}) = logm
inverse_transformation(::Type{NegLogTransform}) = mexp

export logm, mexp


text:
export SincMatrix

#
# A SincMatrix is a sum of products of Diagonal and Toeplitz matrices.
# S = D1*T1 + D2*T2 + ... + Dm*Tm
# Diagonal and Toeplitz matrices can be stored as vectors. Therefore,
# m diagonal or Toeplitz matrices can be stored as a single matrix.
#

type SincMatrix{T} <: AbstractMatrix{T}
    D::Matrix{T}
    TM::Matrix{T}
    TFFT::Matrix{Complex{T}}
end

function SincMatrix{T}(D::Matrix{T},TM::Matrix{T})
    TMn,TMm = size(TM)
    Dn,Dm = size(D)
    TMn == 2Dn-1 && TMm == Dm || throw(DimensionMismatch(""))
    n,m=int((Dn-1)/2),Dm
    pad = zeros(T,4n-1)
    TFFT = zeros(Complex{T},8n,m)
    for i =1:m
        TFFT[:,i] = fft([TM[:,i],pad])
    end
    SincMatrix(D,TM,TFFT)
end

SincMatrix{T}(D::Matrix{T},TM::Vector{T}) = SincMatrix(D,reshape(TM,length(TM),1))
SincMatrix{T}(D::Vector{T},TM::Matrix{T}) = SincMatrix(reshape(D,length(D),1),TM)
SincMatrix{T}(D::Vector{T},TM::Vector{T}) = SincMatrix(reshape(D,length(D),1),reshape(TM,length(TM),1))
SincMatrix{T}(TM::Vector{T}) = SincMatrix(ones(T,div(size(TM)[1]+1,2),1),reshape(TM,length(TM),1))
SincMatrix{T}(TM::Matrix{T}) = SincMatrix(vec(sum(TM,2)))

Base.ndims(A::SincMatrix) = 2
Base.size(A::SincMatrix) = size(A.D)[1],size(A.D)[1]
Base.size(A::SincMatrix,d::Integer) = d == 1 || d == 2 ? size(A.D)[1] : throw(DimensionMismatch(""))

function Base.getindex{T}(A::SincMatrix{T},k::Int64,j::Int64)
    n,m=size(A.D)
    sum(A.D[k,1:m]*A.TM[n+k-j,1:m])
end

function (*){T1<:Number,T2<:Number}(A::SincMatrix{T1},x::Vector{T2})
    nx = int((length(x)-1)/2)
    n,m = size(A.D)
    n = int((n-1)/2)
    n == nx || throw(DimensionMismatch(""))
    T = promote_type(T1,T2)
    xfft = fft([x,zeros(T,6n-1)])
    temp = zeros(T,2n+1)
    for i=1:m
        temp += A.D[:,i].*real(ifft(A.TFFT[:,i].*xfft))[2n+1:4n+1]
    end
    return temp
end


text:
## Copyright 2013 Nick Polson, James Scott, and Jesse Windle.

## This file is part of BayesLogit, distributed under the GNU General Public
## License version 3 or later and without ANY warranty, implied or otherwise.

################################################################################
## CALCULATE SADDLE POINT APPROXIMATION ##
################################################################################

mutable struct FD
    val::Float64
    der::Float64
    FD() = new(0.0, 0.0)
    FD(val, der) = new(val, der)
end

mutable struct Line
    slope::Float64
    icept::Float64
    Line() = new(0.0, 0.0)
    Line(slope, icept) = new(slope, icept)
end

struct PolyaGammaSPSampler <:
       Distributions.Sampleable{Distributions.Univariate,Distributions.Continuous}
    h::Float64
    z::Float64
    maxiter::Int64
    PolyaGammaSPSampler(h, z) = new(h, z, 1000)
    PolyaGammaSPSampler(h, z, maxiter) = new(h, z, maxiter)
end

function cos_rt_sp(v::Float64)
    y = 0.0
    r = sqrt(abs(v))
    if (v >= 0)
        y = cos(r)
    else
        y = cosh(r)
    end
    return y
end

function sp_approx_1(x::Float64, n::Float64, z::Float64)
    v = v_eval(x, 1e-9, 1000)
    u = v / 2
    t = u + z^2 / 2
    m = y_eval(-z^2)

    phi = _log_c(cosh(z)) - _log_c(cos_rt_sp(v)) - t * x

    K2 = x^2 + (1 - x) / (2 * u)
    if (u < 1e-5) && (u > -1e-5)
        K2 = x^2 - (1 / 3) - (2 / 15) * (2 * u)
    end
    spa = sqrt((n / (2 * __PI)) / K2) * _exp_c(phi * n)
    return spa
end

function delta_func_sp(x::Float64, mid::Float64)
    if (x >= mid)
        return FD(_log_c(x) - _log_c(mid), 1 / x)
    else
        return FD(0.5 * (1 - 1 / x) - 0.5 * (1 - 1 / mid), 0.5 / (x^2))
    end
end

function phi_func_sp(x::Float64, z::Float64)
    v = v_eval(x, 1e-9, 1000)
    #println("v =", v)
    u = v / 2
    t = u + (z^2 / 2)
    return FD(_log_c(cosh(abs(z))) - _log_c(cos_rt_sp(v)) - t * x, -t)
end

function tangent_lines_eta(x::Float64, z::Float64, mid::Float64)
    eta = FD()
    delta = FD()
    phi = FD()

    phi = phi_func_sp(x, z)
    delta = delta_func_sp(x, mid)

    eta.val = phi.val - delta.val
    eta.der = phi.der - delta.der

    return Line(eta.der, eta.val - eta.der * x)
end
################################################################################
## rrtigauss ##
################################################################################

## pigauss - cumulative distribution function for Inv-Gauss(mu, lambda).
## NOTE: note using mu = mean, using Z = 1/mu.
##------------------------------------------------------------------------------
#Distributions.pdf(InverseGaussian(Z, lambda), x)

function rrtinvch2_1(rng::Distributions.AbstractRNG, scale::Float64, trnc::Float64)
    R = trnc / scale
    #E = rtnorm(1, 0, 1, left=1/sqrt(R), right=Inf)
    E = Distributions.rand(rng, Distributions.TruncatedNormal(0.0, 1.0, 1 / sqrt(R), Inf))
    X = scale / E^2
    return X
end

function rigauss_1(rng::Distributions.AbstractRNG, mu::Float64, lambda::Float64)
    nu = Random.randn(rng)
    y = nu^2
    x =
        mu + 0.5 * mu^2 * y / lambda -
        0.5 * mu / lambda * sqrt(4 * mu * lambda * y + (mu * y)^2)
    if (Random.rand(rng) > mu / (mu + x))
        x = mu^2 / x
    end
    x
end

function rrtigauss_1(
    rng::Distributions.AbstractRNG,
    mu::Float64,
    lambda::Float64,
    trnc::Float64,
)
    ## trnc is truncation point
    accept = false
    X = trnc + 1
    if (trnc < mu)
        alpha = 0.0
        while (!accept)
            X = rrtinvch2_1(rng, lambda, trnc)
            l_alpha = -0.5 * lambda / mu^2 * X
            accept = _log_c(Random.rand(rng)) < l_alpha
        end
        ## cat("rtigauss.ch, part i:", X, "\n");
    else  ## trnc >= mu
        while (X > trnc)
            Y = Random.randn(rng)^2
            X =
                mu + 0.5 * mu^2 / lambda * Y -
                0.5 * mu / lambda * sqrt(4 * mu * lambda * Y + (mu * Y)^2)
            if (Random.rand(rng) > (mu / (mu + X)))
                X = mu^2 / X
            end
        end
    end
    return X
end
function rltgamma_dagpunar_1(
    rng::Distributions.AbstractRNG,
    shape::Float64,
    rate::Float64,
    trnc::Float64,
)
    ## y ~ Ga(shape, rate, trnc)
    ## x = y/t
    ## x ~ Ga(shape, rate t, 1)
    a = shape
    b = rate * trnc

    if (shape == 1)
        return (Random.randexp(rng) / rate + trnc)
    end

    d1 = b - a
    d3 = a - 1
    c0 = 0.5 * (d1 + sqrt(d1^2 + 4 * b)) / b

    x = 0
    accept = false

    while (!accept)
        x = b + Random.randexp(rng) / c0
        u = Random.rand(rng)

        l_rho = d3 * _log_c(x) - x * (1 - c0)
        l_M = d3 * _log_c(d3 / (1 - c0)) - d3

        accept = _log_c(u) <= (l_rho - l_M)
    end
    return trnc * x / b
end
function pigauss(x::Float64, Z::Float64, lambda::Float64)

    ## I believe this works when Z = 0
    ##Z = 1/mu
    b = sqrt(lambda / x) * (x * Z - 1)
    a = sqrt(lambda / x) * (x * Z + 1) * -1.0
    y =
        Distributions.cdf(Distributions.Normal(), b) +
        exp(2 * lambda * Z) * Distributions.cdf(Distributions.Normal(), a)
    # y2 = 2 * pnorm(-1.0 / sqrt(x));
    return y
end
##------------------------------------------------------------------------------
## Generate sample for J^*(n,z) using SP approx.
function sp_sampler_1(
    rng::Distributions.AbstractRNG,
    n::Float64,
    z::Float64;
    maxiter = 1000,
)

    xl = y_eval(-z^2)
    mid = xl * 1.1
    xr = xl * 1.2
    v_mid = v_eval(mid, 1e-9, maxiter)
    K2_mid = mid^2 + (1 - mid) / v_mid
    al = mid^3 / K2_mid
    ar = mid^2 / K2_mid
    ll = tangent_lines_eta(xl, z, mid)
    lr = tangent_lines_eta(xr, z, mid)
    #tl = tangent_lines_eta(xl, xr, z, mid)
    rl = -ll.slope
    rr = -lr.slope
    il = copy(ll.icept)
    ir = copy(lr.icept)

    wl =
        sqrt(al) *
        _exp_c(-n * sqrt(2 * rl) + n * il + 0.5 * n - 0.5 * n * (1 - 1 / mid)) *
        pigauss(mid, sqrt(2 * rl), n)
    #Distributions.cdf(Distributions.InverseGaussian(sqrt(2*rl), n), mid)
    wr =
        sqrt(ar * n / (__PI * 2)) *
        _exp_c(-n * _log_c(n * rr) + n * ir - n * _log_c(mid) + loggamma(n)) *
        (1 - Distributions.cdf(Distributions.Gamma(n, 1 / (n * rr)), mid))

    wt = wl + wr
    pl = wl / wt

    go = true
    iter = 0
    X = 2
    FX = 0

    while (go && iter <= maxiter)
        iter = iter + 1

        if (wt * Random.rand(rng) < wl)
            ## sample left
            X = rrtigauss_1(rng, 1 / sqrt(2 * rl), n, mid)
            ## while (X > mid) X = rigauss.1(1/sqrt(2*rl), n)
            phi_ev = n * (-rl * X + il) + 0.5 * n * ((1 - 1 / X) - (1 - 1 / mid))
            FX = sqrt(al * n / (2 * __PI)) * X^(-1.5) * _exp_c(phi_ev)

        else
            ## sample right
            X = rltgamma_dagpunar_1(rng, n, n * rr, mid)

            #X = Distributions.rand(rng, TruncatedGamma(n, 1 / (n * rr), mid, Inf))

            phi_ev = n * (-rr * X + ir) + n * (_log_c(X) - _log_c(mid))
            FX = sqrt(ar * 0.5 * n / __PI) * _exp_c(phi_ev) / X
        end

        spa = sp_approx_1(X, n, z)

        if (FX * Random.rand(rng) < spa)
            go = false
        end
    end

    return X
end

## draw with saddle point method
function Distributions.rand(rng::Distributions.AbstractRNG, s::PolyaGammaSPSampler)
    n = s.h
    z = 0.5 * s.z

    x = sp_sampler_1(rng, n, z)

    return 0.25 * n * x
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libdeflate_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libdeflate")
JLLWrappers.@generate_main_file("libdeflate", UUID("46979653-d7f6-5232-b59e-dd310c4598de"))
end  # module libdeflate_jll


text:
function ar1_cov(T, rho)
  acv = rho.^(0:T-1);
  acv = acv/(1-rho^2);
  covmat = acv2cov(acv)
  corr = acv./acv[1]
  return covmat, corr
end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

include("cirrange.jl")
export cirrange

include("rad2sec.jl")
export rad2sec

include("sec2rad.jl")
export sec2rad


text:
module RetroCap

export CapStrategy, UpperBound, MonotonicUpperBound
export LatestVersionOption, CapLatestVersion, ExcludeLatestVersion
export add_caps

include("types.jl")

include("Compress/Compress.jl")

include("add_caps.jl")
include("assert.jl")
include("compat_entries.jl")
include("parse_registry.jl")
include("utils.jl")

end # module


text:
include("four_state_function.jl")

parameter_list = []
filename_list = []

# local
parameters = (; U = 1e0, γ = 1e2, κ = 1e-2, ω = 1e2)
filename = "local.jld2"
push!(parameter_list, parameters)
push!(filename_list, filename)
# semilocal
parameters = (; U = 1e0, γ = 1e1, κ = 1e-1, ω = 1e1)
filename = "semilocal.jld2"
push!(parameter_list, parameters)
push!(filename_list, filename)
# nonlocal
parameters = (; U = 1e0, γ = 1e0, κ = 1e0, ω = 1e0)
filename = "nonlocal.jld2"
push!(parameter_list, parameters)
push!(filename_list, filename)
# nonlocal symmetric
parameters = (; U = 1e0, γ = 1e0, κ = 1e0, ω = 0.0)
filename = "nonlocal_symmetric.jld2"
push!(parameter_list, parameters)
push!(filename_list, filename)
# nonlocal more U
parameters = (; U = 1e0, γ = 1e-1, κ = 1e-1, ω = 1e-1)
filename = "nonlocal_more_velocity.jld2"
push!(parameter_list, parameters)
push!(filename_list, filename)
# nonlocal less diffusivity
parameters = (; U = 1e0, γ = 1e-0, κ = 1e-1, ω = 1e-0)
filename = "nonlocal_less_diffusivity.jld2"
push!(parameter_list, parameters)
push!(filename_list, filename)
# nonlocal less transition
parameters = (; U = 1e0, γ = 1e-1, κ = 1e-0, ω = 1e-1)
filename = "nonlocal_less_transition.jld2"
push!(parameter_list, parameters)
push!(filename_list, filename)
# nonlocal more diffisuvity less transition
parameters = (; U = 1e0, γ = 1e-1, κ = 1e1, ω = 1e-1)
filename = "nonlocal_more_diffusivity_less_transition.jld2"
push!(parameter_list, parameters)
push!(filename_list, filename)
# reverse transition
parameters = (; U = 1e0, γ = 1e-0, κ = 1e0, ω = -1e-0)
filename = "nonlocal_reverse.jld2"
push!(parameter_list, parameters)
push!(filename_list, filename)

total_time = 0
for (parameters, filename) in zip(parameter_list, filename_list)
    println("------")
    println("Currently on ", filename)
    tic = Base.time()
    four_state(parameters; M = 8 * 8, N = 16, filename = filename, dirichlet = false)
    toc = Base.time()
    println("Finished running in ", toc - tic, " seconds.")
    global total_time += toc - tic
    println("------")
end
println("Total time to reproduce data = ", total_time, " seconds")

max_stars_repo_name: sandreza/StatisticalNonlocality

max_issues_repo_name: sandreza/StatisticalNonlocality

max_forks_repo_name: sandreza/StatisticalNonlocality


text:
# This program creates a simple 3D vel model for seismic modeling
#
# by GRAM | 8 Apr 2016

using Seismic

# def size
nx = 100
ny = 100
nz = 200
dz = 10 # depth increment
dx = 20
dy = 20

vel = zeros(Float32,nz,nx,ny)

# build model
for i = 1:nz, j = 1:nx, k = 1:ny # overburden 1
	vel[i,j,k] = 2000
end

for i = 1:nz, j = 1:nx, k = 1:ny # overburden 2
	if i > 30
		vel[i,j,k] = 2500
	end
end

for i = 1:nz, j = 1:nx, k = 1:ny # domal pyramid 1
	if ((i > floor(100 - 1.6 * j + 0.016 * j^2)) & (i > floor(100 - 1.6 * k + 0.016 * k^2)))
		vel[i,j,k] = 3000
	end
end

for i = 1:nz, j = 1:nx, k = 1:ny # domal pyramid 2
	if ((i > floor(140 - 1.6 * j + 0.016 * j^2)) & (i > floor(140 - 1.6 * k + 0.016 * k^2)))
		vel[i,j,k] = 4000
	end
end

for i = 1:nz, j = 1:nx, k = 1:ny # basement
	if i >= 170
		vel[i,j,k] = 5000
	end
end

# Here we "pick" erroneously high vels below the target
# zone, hypothetically because the geophysicist didn't
# want to pick a velocity inversion.

for i = 1:nz, j = 1:nx, k = 1:ny # target
	if ((i < 170) & (i > 130) & (i > floor(140 - 1.6 * j + 0.016 * j^2)) & (i > floor(140 - 1.6 * k + 0.016 * k^2)))
		vel[i,j,k] = 4050
	end
end

#= display
print(size(vel), "\n")
print(vel,"\n")
=#

#= write out text file for QC
out = open("../dat/vel_1.txt","w")
for k = 1:ny
	for i = 1:nz
		for j = 1:nx
			showcompact(out, vel[i,j,k])
			write(out," ")
		end
		write(out, "\r\n")
	end
	write(out, "\r\n")
end
close(out)
=#

# write out vel model for Seismic.jl
h = Array(Header,nx*ny)
ex = Seismic.Extent(convert(Int32,nz), convert(Int32,nx), convert(Int32,ny),
	1, 1, 0, 0, 0, 0, 0, convert(Float32,dz), convert(Float32,dx),
	convert(Float32,dy), 1, 1, "Depth", "mx", "my", "", "", "", "", "",
	"", "", "")

for i = 1 : ny
	for j = 1 : nx
		h[j + nx * (i - 1)] = Seismic.InitSeisHeader()
		h[j + nx * (i - 1)].tracenum = convert(Int32, j + nx * (i - 1))
		h[j + nx * (i - 1)].o1 = convert(Float32, 0)
		h[j + nx * (i - 1)].n1 = convert(Int32, nz)
		h[j + nx * (i - 1)].d1 = convert(Float32, dz)
		h[j + nx * (i - 1)].mx = convert(Float32, j)
		h[j + nx * (i - 1)].my = convert(Float32, i)
	end
end

SeisWrite("../dat/vel_incorrect",vel,h,ex)