text:
"""
SeaPearl Provides several pre-implemented search strategy to explore the search tree in a certain way. The specific strategy has to
be specified in the search! function. The search strategy choose how the _toCall_ Stack will be filled and empty ( ie. which branch will
be explored at a certain node ).
"""
abstract type SearchStrategy end

struct DFSearch <: SearchStrategy end

"""
    struct ILDSearch <: SearchStrategy
implements the basic version of the Iterative Limited Discrepancy Search, ddue to the fact that the depth of the search tree is unknow
before the search we cannot use the improved version of the Iterate Limited Discrepancy Search by Richard E. Korf. d is the max-discrepancy.
"""
struct ILDSearch <: SearchStrategy
    d::Int64
end

abstract type ExpandCriteria end

abstract type RBSearch{C} <: SearchStrategy where C <: ExpandCriteria end

"""
    struct staticRBSearch <: RBSearch
implements the static Restart-Based strategy where the stopping criteria `L` remains the same at each restart.
"""
struct staticRBSearch{C} <: RBSearch{C}
    L::Int64
    n::Int64
end

"""
    struct geometricRBSearch <: RBSearch
implements the geometric Restart-Based strategy where the stopping criteria `L` is increased by the geometric factor `α` at each restart.
"""
struct geometricRBSearch{C}<: RBSearch{C}
    L::Int64
    n::Int64
    α::Float32
end
"""
    struct lubyRBSearch <: RBSearch
implements the Luby Restart-Based strategy where the stopping criteria `L` is multiplied by the factor `Luby[i]` for the i-th restart.
The Luby sequence is a sequence of the following form: 1,1,2,1,1,2,4,1,1,2,1,1,2,4,8, . .and gives theoretical improvement on the search
in the general case.
"""
struct lubyRBSearch{C} <: RBSearch{C}
    L::Int64
    n::Int64
end

abstract type VisitedNodeCriteria <: ExpandCriteria end

"""
    function (criteria<:ExpandCriteria)(model::CPModel, limit::Int64)::Bool

Return a boolean, with `true` meaning that the stopping criteria has been reached
"""
function (criteria::RBSearch{VisitedNodeCriteria})(model::CPModel, limit::Int64)::Bool
    return model.statistics.numberOfNodesBeforeRestart <= limit
    # The stopping criteria is the number of Visited nodes in the search tree.
    # Here the inequality is large because we don't consider the root node of the tree as a visited node.
end

abstract type InfeasibleNodeCriteria <: ExpandCriteria end

function (criteria::RBSearch{InfeasibleNodeCriteria})(model::CPModel, limit::Int64)::Bool
    return model.statistics.numberOfInfeasibleSolutionsBeforeRestart < limit
    # The stopping criteria is the number of visited nodes where the partial solution led to a case where all constraints were not respected.
    # (ie. leading to the :Infeasible case in the expandRbs! function)
end

abstract type SolutionFoundCriteria <: ExpandCriteria end

function (criteria::RBSearch{SolutionFoundCriteria})(model::CPModel, limit::Int64)::Bool
    return model.statistics.numberOfSolutionsBeforeRestart < limit
    # The stopping criteria is the number of Solution Found during the search.
    # As long as L solutions have been found, the search restart at the top of the tree.
end


text:
#
# Copyright (c) 2021 Frederic Bruder, Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using Test
using FMIZoo

@testset "FMIZoo.jl" begin
    list_models()

    # get_model_filename#1
    path = get_model_filename("SpringDamperPendulum1D", "Dymola", "2022x")
    @test length(path) > 0
    split = splitpath(path)

    @test split[end] ==  "SpringDamperPendulum1D.fmu"
    @test split[end-1] ==  "2.0"
    @test split[end-2] ==  "2022x"
    @test split[end-3] ==  "Dymola"
    @test split[end-4] ==  "bin"
    @test split[end-5] ==  "models"
    @test split[end-6] ==  "FMIZoo.jl"

    path = get_model_filename("BouncingBall", "ModelicaReferenceFMUs", "0.0.14", "3.0")
    @test length(path) > 0
    split = splitpath(path)

    @test split[end] ==  "BouncingBall.fmu"
    @test split[end-1] ==  "BouncingBall"

    # generate_mos_scripts
    generate_mos_scripts(verbose=false)
    @test isfile(joinpath(FMIZoo.p_mos_scripts, "Dymola2022x.mos"))
end

ext: jl


text:
struct UniformPeriodicFluxDiffDisc2D{BalanceLaw, T, Fvol, Fnum, Fluxes, UseThreads<:Union{Val{true},Val{false}}} <: AbstractSemidiscretisation
    balance_law::BalanceLaw
    meshx::UniformPeriodicMesh1D{T}
    meshy::UniformPeriodicMesh1D{T}
    basis::LobattoLegendre{T}
    fvol::Fvol
    fnum::Fnum

    fluxes::Fluxes
    usethreads::UseThreads

    function UniformPeriodicFluxDiffDisc2D(balance_law::BalanceLaw, meshx::UniformPeriodicMesh1D{T}, meshy::UniformPeriodicMesh1D{T}, basis::LobattoLegendre{T}, fvol::Fvol, fnum::Fnum, fluxes::Fluxes, usethreads::UseThreads) where {BalanceLaw, T, Fvol, Fnum, Fluxes, UseThreads}
        @assert typeof(fluxes) <: AbstractArray{variables(balance_law), 4}
        @assert size(fluxes, 1) == length(basis.nodes)
        @assert size(fluxes, 2) == 2
        @assert size(fluxes, 3) == numedges(meshx)
        @assert size(fluxes, 4) == numedges(meshy)
        new{BalanceLaw, T, Fvol, Fnum, Fluxes, UseThreads}(balance_law, meshx, meshy, basis, fvol, fnum, fluxes, usethreads)
    end
end

function UniformPeriodicFluxDiffDisc2D(balance_law, meshx, meshy, basis, fvol, fnum, usethreads::Bool=false)
    fluxes = Array{variables(balance_law)}(undef, length(basis.nodes), 2, numedges(meshx), numedges(meshy))
    UniformPeriodicFluxDiffDisc2D(balance_law, meshx, meshy, basis, fvol, fnum, fluxes, Val{usethreads}())
end

@noinline function (semidisc::UniformPeriodicFluxDiffDisc2D)(du, u, p, t)
  @boundscheck begin
    if size(u) != size(du)
      error("size(u) = $(size(u)) != $(size(du)) = size(du)")
    end
    @assert size(u,1) == size(u,2) == length(semidisc.basis.nodes)
    @assert size(u,3) == numcells(semidisc.meshx)
    @assert size(u,4) == numcells(semidisc.meshy)
    if eltype(u) != variables(semidisc.balance_law)
      error("eltype(u) == $(eltype(u)) != $(variables(semidisc.balance_law)) == variables(semidisc.balance_law)")
    end
  end

  fill!(du, zero(eltype(du)))
  add_flux_differences!(du, u, semidisc)
  add_numerical_fluxes!(du, u, semidisc)

  nothing
end

function add_flux_differences!(du, u, semidisc::UniformPeriodicFluxDiffDisc2D)
    Nx = size(u, 3)
    Ny = size(u, 4)
    Pp1 = size(u, 1)
    @boundscheck @assert Pp1 == size(u, 2)

    @unpack balance_law, fvol, basis, usethreads = semidisc
    @unpack D = basis
    @boundscheck @assert Pp1 == size(D,1) == size(D,2)
    jacx = 2 / semidisc.meshx.Δx
    jacy = 2 / semidisc.meshy.Δx

    if Pp1 <= 1
        return nothing
    end

    add_flux_differences_inner_loop!(du, u, balance_law, fvol, Nx, Ny, Pp1, D, jacx, jacy, usethreads)
end

@inline function add_flux_differences_inner_loop!(du, u, balance_law, fvol, Nx, Ny, Pp1,
                                                    D, jacx, jacy, ::Val{true})
    dirx = Val{:x}()
    diry = Val{:y}()
    NxNy = Nx*Ny
    dims = (Pp1,Pp1,Nx*Ny)
    @inbounds Threads.@threads for ixy in Base.OneTo(NxNy)
        for ny in Base.OneTo(Pp1), nx in Base.OneTo(Pp1)
            idx = LinearIndices(dims)[nx, ny, ixy]
            u_idx = u[idx]
            # compute x derivative
            # at first for different indices
            for k in (nx+1):Pp1
                idxk = LinearIndices(dims)[k, ny, ixy]
                f = fvol(u_idx, u[idxk], balance_law, dirx)
                du[idx]  -= 2*jacx*D[nx,k] * f
                du[idxk] -= 2*jacx*D[k,nx] * f
            end
            # then for the diagonal element, using the consistency of the flux
            du[idx] -= 2*jacx*D[nx,nx] * flux(u_idx, balance_law, dirx)

            # compute y derivative
            # at first for different indices
            for k in (ny+1):Pp1
                idxk = LinearIndices(dims)[nx, k, ixy]
                f = fvol(u_idx, u[idxk], balance_law, diry)
                du[idx]  -= 2*jacy*D[ny,k] * f
                du[idxk] -= 2*jacy*D[k,ny] * f
            end
            # then for the diagonal element, using the consistency of the flux
            du[idx] -= 2*jacy*D[ny,ny] * flux(u_idx, balance_law, diry)
        end
    end
    nothing
end

@inline function add_flux_differences_inner_loop!(du, u, balance_law, fvol, Nx, Ny, Pp1,
                                                    D, jacx, jacy, ::Val{false})
    dirx = Val{:x}()
    diry = Val{:y}()
    NxNy = Nx*Ny
    dims = (Pp1,Pp1,Nx*Ny)
    @inbounds for ixy in Base.OneTo(NxNy)
        for ny in Base.OneTo(Pp1), nx in Base.OneTo(Pp1)
            idx = LinearIndices(dims)[nx, ny, ixy]
            u_idx = u[idx]
            # compute x derivative
            # at first for different indices
            for k in (nx+1):Pp1
                idxk = LinearIndices(dims)[k, ny, ixy]
                f = fvol(u_idx, u[idxk], balance_law, dirx)
                du[idx]  -= 2*jacx*D[nx,k] * f
                du[idxk] -= 2*jacx*D[k,nx] * f
            end
            # then for the diagonal element, using the consistency of the flux
            du[idx] -= 2*jacx*D[nx,nx] * flux(u_idx, balance_law, dirx)

            # compute y derivative
            # at first for different indices
            for k in (ny+1):Pp1
                idxk = LinearIndices(dims)[nx, k, ixy]
                f = fvol(u_idx, u[idxk], balance_law, diry)
                du[idx]  -= 2*jacy*D[ny,k] * f
                du[idxk] -= 2*jacy*D[k,ny] * f
            end
            # then for the diagonal element, using the consistency of the flux
            du[idx] -= 2*jacy*D[ny,ny] * flux(u_idx, balance_law, diry)
        end
    end
    nothing
end

function add_numerical_fluxes!(du, u, semidisc::UniformPeriodicFluxDiffDisc2D)
    Nx = size(u, 3)
    Ny = size(u, 4)
    Pp1 = size(u, 1)
    @boundscheck @assert Pp1 == size(u, 2)

    @unpack balance_law, meshx, meshy, fnum, fluxes, usethreads = semidisc
    ω = semidisc.basis.weights
    @boundscheck @assert Pp1 == length(ω)

    jacx = 2 / semidisc.meshx.Δx
    jacy = 2 / semidisc.meshy.Δx

    add_numerical_fluxes_inner_loop1!(du, fluxes, u, balance_law, fnum, Nx, Ny, Pp1, ω, jacx, jacy, usethreads)
    add_numerical_fluxes_inner_loop2!(du, fluxes, u, balance_law, fnum, Nx, Ny, Pp1, ω, jacx, jacy, usethreads)
end

@inline function add_numerical_fluxes_inner_loop1!(du, fluxes, u, balance_law, fnum,
                                                    Nx, Ny, Pp1, ω, jacx, jacy, ::Val{true})
    dirx = Val{:x}()
    diry = Val{:y}()
    dims = (Pp1, Pp1, Nx, Ny)
    # calculate numerical fluxes
    @inbounds Threads.@threads for ixy in Base.OneTo(Nx*Ny)
        ix, iy = Tuple(CartesianIndices((Nx,Ny))[ixy])
        ixm1 = ix ==  1 ? Nx : ix-1
        iym1 = iy ==  1 ? Ny : iy-1

        # flux x
        for ny in 1:Pp1
            # flux x - left
            fluxes[ny,1,ix,iy] = fnum(u[end,ny,ixm1,iy], u[1,ny,ix,iy], balance_law, dirx)
        end

        # flux y
        for nx in 1:Pp1
            # flux y - bottom
            fluxes[nx,2,ix,iy] = fnum(u[nx,end,ix,iym1], u[nx,1,ix,iy], balance_law, diry)
        end
    end
    nothing
end
@inline function add_numerical_fluxes_inner_loop2!(du, fluxes, u, balance_law, fnum,
                                                    Nx, Ny, Pp1, ω, jacx, jacy, ::Val{true})
    dirx = Val{:x}()
    diry = Val{:y}()
    @inbounds i_ω1 = 1 / ω[1]
    @inbounds i_ωend = 1 / ω[end]
    dims = (Pp1, Pp1, Nx, Ny)
    # add numerical fluxes
    @inbounds Threads.@threads for ixy in Base.OneTo(Nx*Ny)
        ix, iy = Tuple(CartesianIndices((Nx,Ny))[ixy])
        ixp1 = ix == Nx ?  1 : ix+1
        iyp1 = iy == Ny ?  1 : iy+1

        # flux x
        for ny in 1:Pp1
            # flux x - left
            idx = LinearIndices(dims)[1, ny, ix, iy]
            du[idx] += ( fluxes[ny,1,ix,iy] - flux(u[idx], balance_law, dirx) ) * jacx * i_ω1

            # flux x - right
            idx = LinearIndices(dims)[Pp1, ny, ix, iy]
            du[idx] -= ( fluxes[ny,1,ixp1,iy] - flux(u[idx], balance_law, dirx) ) * jacx * i_ωend
        end

        # flux y
        for nx in 1:Pp1
            # flux y - bottom
            idx = LinearIndices(dims)[nx, 1, ix, iy]
            du[idx] += ( fluxes[nx,2,ix,iy] - flux(u[idx], balance_law, diry) ) * jacy * i_ω1

            # flux y - top
            idx = LinearIndices(dims)[nx, Pp1, ix, iy]
            du[idx] -= ( fluxes[nx,2,ix,iyp1] - flux(u[idx], balance_law, diry) ) * jacy * i_ωend
        end
    end
    nothing
end

@inline function add_numerical_fluxes_inner_loop1!(du, fluxes, u, balance_law, fnum,
                                                    Nx, Ny, Pp1, ω, jacx, jacy, ::Val{false})
    dirx = Val{:x}()
    diry = Val{:y}()
    dims = (Pp1, Pp1, Nx, Ny)
    # calculate numerical fluxes
    @inbounds for ixy in Base.OneTo(Nx*Ny)
        ix, iy = Tuple(CartesianIndices((Nx,Ny))[ixy])
        ixm1 = ix ==  1 ? Nx : ix-1
        iym1 = iy ==  1 ? Ny : iy-1

        # flux x
        for ny in 1:Pp1
            # flux x - left
            fluxes[ny,1,ix,iy] = fnum(u[end,ny,ixm1,iy], u[1,ny,ix,iy], balance_law, dirx)
        end

        # flux y
        for nx in 1:Pp1
            # flux y - bottom
            fluxes[nx,2,ix,iy] = fnum(u[nx,end,ix,iym1], u[nx,1,ix,iy], balance_law, diry)
        end
    end
    nothing
end
@inline function add_numerical_fluxes_inner_loop2!(du, fluxes, u, balance_law, fnum,
                                                    Nx, Ny, Pp1, ω, jacx, jacy, ::Val{false})
    dirx = Val{:x}()
    diry = Val{:y}()
    @inbounds i_ω1 = 1 / ω[1]
    @inbounds i_ωend = 1 / ω[end]
    dims = (Pp1, Pp1, Nx, Ny)
    # add numerical fluxes
    @inbounds for ixy in Base.OneTo(Nx*Ny)
        ix, iy = Tuple(CartesianIndices((Nx,Ny))[ixy])
        ixp1 = ix == Nx ?  1 : ix+1
        iyp1 = iy == Ny ?  1 : iy+1

        # flux x
        for ny in 1:Pp1
            # flux x - left
            idx = LinearIndices(dims)[1, ny, ix, iy]
            du[idx] += ( fluxes[ny,1,ix,iy] - flux(u[idx], balance_law, dirx) ) * jacx * i_ω1

            # flux x - right
            idx = LinearIndices(dims)[Pp1, ny, ix, iy]
            du[idx] -= ( fluxes[ny,1,ixp1,iy] - flux(u[idx], balance_law, dirx) ) * jacx * i_ωend
        end

        # flux y
        for nx in 1:Pp1
            # flux y - bottom
            idx = LinearIndices(dims)[nx, 1, ix, iy]
            du[idx] += ( fluxes[nx,2,ix,iy] - flux(u[idx], balance_law, diry) ) * jacy * i_ω1

            # flux y - top
            idx = LinearIndices(dims)[nx, Pp1, ix, iy]
            du[idx] -= ( fluxes[nx,2,ix,iyp1] - flux(u[idx], balance_law, diry) ) * jacy * i_ωend
        end
    end
    nothing
end

function semidiscretise(semidisc::UniformPeriodicFluxDiffDisc2D, u₀func, tspan)
    @unpack meshx, meshy, basis = semidisc
    u₀ = compute_coefficients(u₀func, meshx, meshy, basis)

    ODEProblem(semidisc, u₀, tspan)
end

"""
    max_dt(t, u, semidisc::UniformPeriodicFluxDiffDisc2D, cfl)

Compute the maximal time step `dt` satisfying the CFL condition `dt <= cfl * dx / ((2p+1)*λ)`,
where `dx` is the length of a cell, `λ` the greatest absolute value of the
speed in this cell, and `p` the polynomial degree.
"""
function max_dt(t, u, semidisc::UniformPeriodicFluxDiffDisc2D, cfl=0.5)
    @unpack balance_law, meshx, meshy, basis = semidisc
    Δx = min(meshx.Δx, meshy.Δx)
    p = length(basis.nodes)-1
    factor = Δx / (2p+1)

    dt = mapreduce(u->factor/max_abs_speed(u, balance_law), min, u, init=typemax(t))
    if dt == Inf
        dt = factor
    end
    cfl * dt
end

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

struct UniformPeriodicFluxDiffDisc3D{BalanceLaw, T, Fvol, Fnum, UseThreads<:Union{Val{true},Val{false}}} <: AbstractSemidiscretisation
    balance_law::BalanceLaw
    meshx::UniformPeriodicMesh1D{T}
    meshy::UniformPeriodicMesh1D{T}
    meshz::UniformPeriodicMesh1D{T}
    basis::LobattoLegendre{T}
    fvol::Fvol
    fnum::Fnum

    usethreads::UseThreads
end

function UniformPeriodicFluxDiffDisc3D(balance_law, meshx, meshy, meshz, basis, fvol, fnum, usethreads::Bool=false)
    UniformPeriodicFluxDiffDisc3D(balance_law, meshx, meshy, meshz, basis, fvol, fnum, Val{usethreads}())
end

@noinline function (semidisc::UniformPeriodicFluxDiffDisc3D)(du, u, p, t)
  @boundscheck begin
    if size(u) != size(du)
      error("size(u) = $(size(u)) != $(size(du)) = size(du)")
    end
    @assert size(u,1) == size(u,2) == size(u,3) == length(semidisc.basis.nodes)
    @assert size(u,4) == numcells(semidisc.meshx)
    @assert size(u,5) == numcells(semidisc.meshy)
    @assert size(u,6) == numcells(semidisc.meshz)
    if eltype(u) != variables(semidisc.balance_law)
      error("eltype(u) == $(eltype(u)) != $(variables(semidisc.balance_law)) == variables(semidisc.balance_law)")
    end
  end

  fill!(du, zero(eltype(du)))
  add_flux_differences!(du, u, semidisc)
  add_numerical_fluxes!(du, u, semidisc)

  nothing
end

function add_flux_differences!(du, u, semidisc::UniformPeriodicFluxDiffDisc3D)
    Nx = size(u, 4)
    Ny = size(u, 5)
    Nz = size(u, 6)
    Pp1 = size(u, 1)
    @boundscheck @assert Pp1 == size(u, 2) == size(u, 3)

    @unpack balance_law, fvol, basis, usethreads = semidisc
    @unpack D = basis
    @boundscheck @assert Pp1 == size(D,1) == size(D,2)
    jacx = 2 / semidisc.meshx.Δx
    jacy = 2 / semidisc.meshy.Δx
    jacz = 2 / semidisc.meshz.Δx

    if Pp1 <= 1
        return nothing
    end

    add_flux_differences_inner_loop!(du, u, balance_law, fvol, Nx, Ny, Nz, Pp1, D, jacx, jacy, jacz, usethreads)
end

@inline function add_flux_differences_inner_loop!(du, u, balance_law, fvol, Nx, Ny, Nz,
                                                    Pp1, D, jacx, jacy, jacz, ::Val{true})
    dirx = Val{:x}()
    diry = Val{:y}()
    dirz = Val{:z}()
    NxNyNz = Nx*Ny*Nz
    dims = (Pp1,Pp1,Pp1,NxNyNz)
    @inbounds Threads.@threads for ixyz in Base.OneTo(NxNyNz)
        for nz in Base.OneTo(Pp1), ny in Base.OneTo(Pp1), nx in Base.OneTo(Pp1)
            idx = LinearIndices(dims)[nx, ny, nz, ixyz]

            # compute x derivative
            # at first for different indices
            for k in (nx+1):Pp1
                idxk = LinearIndices(dims)[k, ny, nz, ixyz]
                f = fvol(u_idx, u[idxk], balance_law, dirx)
                du[idx]  -= 2*jacx*D[nx,k] * f
                du[idxk] -= 2*jacx*D[k,nx] * f
            end
            # then for the diagonal element, using the consistency of the flux
            du[idx] -= 2*jacx*D[nx,nx] * flux(u_idx, balance_law, dirx)

            # compute y derivative
            # at first for different indices
            for k in (ny+1):Pp1
                idxk = LinearIndices(dims)[nx, k, nz, ixyz]
                f = fvol(u_idx, u[idxk], balance_law, diry)
                du[idx]  -= 2*jacy*D[ny,k] * f
                du[idxk] -= 2*jacy*D[k,ny] * f
            end
            # then for the diagonal element, using the consistency of the flux
            du[idx] -= 2*jacy*D[ny,ny] * flux(u_idx, balance_law, diry)

            # compute z derivative
            # at first for different indices
            for k in (nz+1):Pp1
                idxk = LinearIndices(dims)[nx, ny, k, ixyz]
                f = fvol(u_idx, u[idxk], balance_law, dirz)
                du[idx]  -= 2*jacz*D[nz,k] * f
                du[idxk] -= 2*jacz*D[k,nz] * f
            end
            # then for the diagonal element, using the consistency of the flux
            du[idx] -= 2*jacz*D[nz,nz] * flux(u_idx, balance_law, dirz)
        end
    end
    nothing
end

@inline function add_flux_differences_inner_loop!(du, u, balance_law, fvol, Nx, Ny, Nz,
                                                    Pp1, D, jacx, jacy, jacz, ::Val{false})
    dirx = Val{:x}()
    diry = Val{:y}()
    dirz = Val{:z}()
    NxNyNz = Nx*Ny*Nz
    dims = (Pp1,Pp1,Pp1,NxNyNz)
    @inbounds for ixyz in Base.OneTo(NxNyNz)
        for nz in Base.OneTo(Pp1), ny in Base.OneTo(Pp1), nx in Base.OneTo(Pp1)
            idx = LinearIndices(dims)[nx, ny, nz, ixyz]
            u_idx = u[idx]

            # compute x derivative
            # at first for different indices
            for k in (nx+1):Pp1
                idxk = LinearIndices(dims)[k, ny, nz, ixyz]
                f = fvol(u_idx, u[idxk], balance_law, dirx)
                du[idx]  -= 2*jacx*D[nx,k] * f
                du[idxk] -= 2*jacx*D[k,nx] * f
            end
            # then for the diagonal element, using the consistency of the flux
            du[idx] -= 2*jacx*D[nx,nx] * flux(u_idx, balance_law, dirx)

            # compute y derivative
            # at first for different indices
            for k in (ny+1):Pp1
                idxk = LinearIndices(dims)[nx, k, nz, ixyz]
                f = fvol(u_idx, u[idxk], balance_law, diry)
                du[idx]  -= 2*jacy*D[ny,k] * f
                du[idxk] -= 2*jacy*D[k,ny] * f
            end
            # then for the diagonal element, using the consistency of the flux
            du[idx] -= 2*jacy*D[ny,ny] * flux(u_idx, balance_law, diry)

            # compute z derivative
            # at first for different indices
            for k in (nz+1):Pp1
                idxk = LinearIndices(dims)[nx, ny, k, ixyz]
                f = fvol(u_idx, u[idxk], balance_law, dirz)
                du[idx]  -= 2*jacz*D[nz,k] * f
                du[idxk] -= 2*jacz*D[k,nz] * f
            end
            # then for the diagonal element, using the consistency of the flux
            du[idx] -= 2*jacz*D[nz,nz] * flux(u_idx, balance_law, dirz)
        end
    end
    nothing
end

function add_numerical_fluxes!(du, u, semidisc::UniformPeriodicFluxDiffDisc3D)
    Nx = size(u, 4)
    Ny = size(u, 5)
    Nz = size(u, 6)
    Pp1 = size(u, 1)
    @boundscheck @assert Pp1 == size(u, 2) == size(u, 3)

    @unpack balance_law, meshx, meshy, fnum, usethreads = semidisc
    ω = semidisc.basis.weights
    @boundscheck @assert Pp1 == length(ω)

    jacx = 2 / semidisc.meshx.Δx
    jacy = 2 / semidisc.meshy.Δx
    jacz = 2 / semidisc.meshz.Δx

    add_numerical_fluxes_inner_loop!(du, u, balance_law, fnum, Nx, Ny, Nz, Pp1, ω, jacx, jacy, jacz, usethreads)
end

@inline function add_numerical_fluxes_inner_loop!(du, u, balance_law, fnum, Nx, Ny, Nz,
                                                    Pp1, ω, jacx, jacy, jacz, ::Val{true})
    dirx = Val{:x}()
    diry = Val{:y}()
    dirz = Val{:y}()
    @inbounds Threads.@threads for ixyz in Base.OneTo(Nx*Ny*Nz)
        ix, iy, iz = Tuple(CartesianIndices((Nx,Ny,Nz))[ixyz])
        ixm1 = ix ==  1 ? Nx : ix-1
        ixp1 = ix == Nx ?  1 : ix+1
        iym1 = iy ==  1 ? Ny : iy-1
        iyp1 = iy == Ny ?  1 : iy+1
        izm1 = iz ==  1 ? Nz : iz-1
        izp1 = iz == Nz ?  1 : iz+1

        # flux x
        for nz in Base.OneTo(Pp1), ny in Base.OneTo(Pp1)
            # flux x - left
            du[1,ny,nz,ix,iy,iz] += (
                fnum(u[end,ny,nz,ixm1,iy,iz], u[1,ny,nz,ix,iy,iz], balance_law, dirx)
                - flux(u[1,ny,nz,ix,iy,iz], balance_law, dirx) ) * jacx / ω[1]

            # flux x - right
            du[end,ny,nz,ix,iy,iz] -= (
                fnum(u[end,ny,nz,ix,iy,iz], u[1,ny,nz,ixp1,iy,iz], balance_law, dirx)
                - flux(u[end,ny,nz,ix,iy,iz], balance_law, dirx) ) * jacx / ω[end]
        end

        # flux y
        for nz in Base.OneTo(Pp1), nx in Base.OneTo(Pp1)
            # flux y - bottom
            du[nx,1,nz,ix,iy,iz] += (
                fnum(u[nx,end,nz,ix,iym1,iz], u[nx,1,nz,ix,iy,iz], balance_law, diry)
                - flux(u[nx,1,nz,ix,iy,iz], balance_law, diry) ) * jacy / ω[1]

            # flux y - top
            du[nx,end,nz,ix,iy,iz] -= (
                fnum(u[nx,end,nz,ix,iy,iz], u[nx,1,nz,ix,iyp1,iz], balance_law, diry)
                - flux(u[nx,end,nz,ix,iy,iz], balance_law, diry) ) * jacy / ω[end]
        end

        # flux z
        for ny in Base.OneTo(Pp1), nx in Base.OneTo(Pp1)
            # flux z - back
            du[nx,ny,1,ix,iy,iz] += (
                fnum(u[nx,ny,end,ix,iy,izm1], u[nx,ny,1,ix,iy,iz], balance_law, dirz)
                - flux(u[nx,ny,1,ix,iy,iz], balance_law, dirz) ) * jacz / ω[1]

            # flux z - front
            du[nx,ny,end,ix,iy,iz] -= (
                fnum(u[nx,ny,end,ix,iy,iz], u[nx,ny,1,ix,iy,izp1], balance_law, dirz)
                - flux(u[nx,ny,end,ix,iy,iz], balance_law, dirz) ) * jacz / ω[end]
        end
    end
    nothing
end

@inline function add_numerical_fluxes_inner_loop!(du, u, balance_law, fnum, Nx, Ny, Nz,
                                                    Pp1, ω, jacx, jacy, jacz, ::Val{false})
    dirx = Val{:x}()
    diry = Val{:y}()
    dirz = Val{:y}()
    @inbounds for ixyz in Base.OneTo(Nx*Ny*Nz)
        ix, iy, iz = Tuple(CartesianIndices((Nx,Ny,Nz))[ixyz])
        ixm1 = ix ==  1 ? Nx : ix-1
        ixp1 = ix == Nx ?  1 : ix+1
        iym1 = iy ==  1 ? Ny : iy-1
        iyp1 = iy == Ny ?  1 : iy+1
        izm1 = iz ==  1 ? Nz : iz-1
        izp1 = iz == Nz ?  1 : iz+1

        # flux x
        for nz in Base.OneTo(Pp1), ny in Base.OneTo(Pp1)
            # flux x - left
            du[1,ny,nz,ix,iy,iz] += (
                fnum(u[end,ny,nz,ixm1,iy,iz], u[1,ny,nz,ix,iy,iz], balance_law, dirx)
                - flux(u[1,ny,nz,ix,iy,iz], balance_law, dirx) ) * jacx / ω[1]

            # flux x - right
            du[end,ny,nz,ix,iy,iz] -= (
                fnum(u[end,ny,nz,ix,iy,iz], u[1,ny,nz,ixp1,iy,iz], balance_law, dirx)
                - flux(u[end,ny,nz,ix,iy,iz], balance_law, dirx) ) * jacx / ω[end]
        end

        # flux y
        for nz in Base.OneTo(Pp1), nx in Base.OneTo(Pp1)
            # flux y - bottom
            du[nx,1,nz,ix,iy,iz] += (
                fnum(u[nx,end,nz,ix,iym1,iz], u[nx,1,nz,ix,iy,iz], balance_law, diry)
                - flux(u[nx,1,nz,ix,iy,iz], balance_law, diry) ) * jacy / ω[1]

            # flux y - top
            du[nx,end,nz,ix,iy,iz] -= (
                fnum(u[nx,end,nz,ix,iy,iz], u[nx,1,nz,ix,iyp1,iz], balance_law, diry)
                - flux(u[nx,end,nz,ix,iy,iz], balance_law, diry) ) * jacy / ω[end]
        end

        # flux z
        for ny in Base.OneTo(Pp1), nx in Base.OneTo(Pp1)
            # flux z - back
            du[nx,ny,1,ix,iy,iz] += (
                fnum(u[nx,ny,end,ix,iy,izm1], u[nx,ny,1,ix,iy,iz], balance_law, dirz)
                - flux(u[nx,ny,1,ix,iy,iz], balance_law, dirz) ) * jacz / ω[1]

            # flux z - front
            du[nx,ny,end,ix,iy,iz] -= (
                fnum(u[nx,ny,end,ix,iy,iz], u[nx,ny,1,ix,iy,izp1], balance_law, dirz)
                - flux(u[nx,ny,end,ix,iy,iz], balance_law, dirz) ) * jacz / ω[end]
        end
    end
    nothing
end

function semidiscretise(semidisc::UniformPeriodicFluxDiffDisc3D, u₀func, tspan)
    @unpack meshx, meshy, meshz, basis = semidisc
    u₀ = compute_coefficients(u₀func, meshx, meshy, meshz, basis)

    ODEProblem(semidisc, u₀, tspan)
end


text:
using Test
using Domains, JSONSchema,JSON

schema = JSONSchema.Schema(read(Domains.jsonschemafile, String), parent_dir = dirname(Domains.jsonschemafile))

# For now exclude non lambert domains from tests
for domain in filter(d->d ∉ ["SCANDINAVIA_ROTM","NORWAY_POLAR", "RCR_POLAR"], Domains.domains)

    d = JSON.parsefile(joinpath(Domains.jsondir, "$domain.json"))

    @testset "$domain" begin

        # Validate against JSON schema
        @test JSONSchema.isvalid(d, schema)

        # Test that the North Pole [0, 90] is outside the domain with ezone.
        @test !in([0, 90], Domain(d), with_ezone = true)

    end # @testset "$domain" begin

end


text:
# problem 1a
N = 10
A = zeros(N,N)
for i in 1:N, j in 1:N
    abs(i-j)<=1 ? A[i,j]+=1 : nothing
    i==j ? A[i,j]-=3 : nothing
end
A

# problem 1b
sum(sum(1/j for j in 1:i) for i in 1:25)

#### Prepare Data

X = rand(1000, 3)               # feature matrix
a0 = rand(3)                    # ground truths
y = X * a0 + 0.1 * randn(1000);  # generate response

X2 = hcat(ones(1000),X)
lreg = X2\y

# problem 3
using DataFrames
data = DataFrame(X)
data[:y] = y
lm1 = fit(LinearModel, @formula(y ~ x1 + x2 + x3), data)
@test coef(lm1) ≈ lreg atol=1e-16

# problem 4
r = 2.9:.00005:4; numAttract = 150
steady = ones(length(r),1)*.25
for i=1:400 ## Get to steady state
  steady .= r.*steady.*(1-steady)
end
x = zeros(length(steady),numAttract)
x[:,1] = steady
@inbounds for i=2:numAttract ## Grab values at the attractor
  x[:,i] = r.*x[:,i-1].*(1-x[:,i-1])
end
using Plots
plot(collect(r),x,seriestype=:scatter,markersize=.002,legend=false)

macro ~(y,ex)
  new_ex = Meta.quot(ex)
  quote
    inner_ex = $(esc(new_ex))
    data_name = Symbol(string(inner_ex.args[end])[1])
    eval_ex = Expr(:(=),:data,data_name)
    eval(Main,eval_ex)
    new_X = Matrix{Float64}(size(data,1),length(inner_ex.args)-1)
    cur_spot = 0
    for i in 2:length(inner_ex.args)
      if inner_ex.args[i] == 1
        new_X[:,i-1] = ones(size(data,1))
      else
        col = parse(Int,string(string(inner_ex.args[i])[2]))
        new_X[:,i-1] = data[:,col]
      end
    end
    $(esc(y)),new_X
  end
end

y = rand(10)
X = rand(10,4)
y~1+X1+X2+X4

function solve_least_squares(y,X)
  X\y
end
solve_least_squares(tup::Tuple) = solve_least_squares(tup...)
solve_least_squares(y~1+X1+X2+X4)

function myquantile(d::UnivariateDistribution, q::Number)
    θ = mean(d)
    tol = Inf
    while tol > 1e-5
        θold = θ
        θ = θ - (cdf(d, θ) - q) / pdf(d, θ)
        tol = abs(θold - θ)
    end
    θ
end

for dist in [Gamma(5, 1), Normal(0, 1), Beta(2, 4)]
    @show myquantile(dist, .75)
    @show quantile(dist, .75)
    println()
end

using LightGraphs, Distributions
function mkTree(maxdepth::Int = 10, p::Float64 = 0.8, g::SimpleGraph = Graph(1), currhead::Int = 1)
    if (maxdepth <= 1) g
    else
       b = Binomial(2, p)
       nEdges = max(1, rand(b))
        for leaves in 1:nEdges
            add_vertex!(g)
            newnode = nv(g)
            add_edge!(g, currhead, newnode)
            mkTree(maxdepth-1, p, g, newnode)
        end
    end
    g
end

using Roots
f(x) = 10 - x + e*sin(x)
fzero(f,BigFloat(2.0))

f! = function (x,dx)
  dx[1] = x[1]   + x[2]   + x[3]^2 -12
  dx[2] = x[1]^2 - x[2]   + x[3]   - 2
  dx[3] = 2x[1]  - x[2]^2 + x[3]   - 1
end
using NLsolve
res = nlsolve(f!,[1.0;1.0;1.0])
res.zero
res = nlsolve(f!,[1.0;1.0;1.0],autodiff=true)
res.zero

max_stars_repo_name: snowdj/CoursePack

max_issues_repo_name: snowdj/CoursePack

max_forks_repo_name: snowdj/CoursePack


text:
struct Parameter{T} end
isparameter(::Variable) = false
isparameter(::Variable{<:Parameter}) = true

"""
$(SIGNATURES)

Define one or more known variables.
"""
macro parameters(xs...)
    esc(_parse_vars(:parameters, Parameter{Number}, xs))
end


text:
function gradient_st(tree::Matrix{Float64},x::Matrix{Float64},var::Int64)::Vector{Float64}
    N::Int64 = size(x)[1]
    terminalnodes::Vector{Int64} = find(x->(x==1),tree[:,6])
    terminal::Matrix{Float64} = tree[terminalnodes,:]
    Nt::Int64 = length(terminalnodes)
    logimat::Matrix{Float64} = Matrix{Float64}(N,Nt)
    node::Vector{Float64} = Vector{Float64}(9)
    logit::Vector{Float64} = Vector{Float64}(N)
    dlogit::Vector{Float64} = Vector{Float64}(N)
    x_split::Vector{Float64} = Vector{Float64}(N)
    variable::Int64 = 1::Int64
    parent::Int64 = 1::Int64
    logitaux::Vector{Float64}=Vector{Float64}(N)
    dlogitaux::Vector{Float64}=Vector{Float64}(N)
    for i in 1:Nt
        node=terminal[i,:]
        variable = floor(Int,node[7])
        x_split = x[:,variable]
        [logit[j] = 1/(1+exp(-node[4]*(x_split[j]-node[3]))) for j in 1:N]
        if variable == var
            dlogit = node[4]*logit.*(1-logit)
        else
            dlogit = zeros(N)
        end
        if node[1] == 2.0
            logit = 1.0-logit
            dlogit = -dlogit
        end

        parent = floor(Int,node[5])
        while parent!=0
            node = tree[parent,:]
            variable = floor(Int,node[7])
            x_split = x[:,variable]
            [logitaux[j] = 1/(1+exp(-node[4]*(x_split[j]-node[3]))) for j in 1:N]
            if variable == var
                dlogitaux = node[4]*logitaux.*(1-logitaux)
            else
                dlogitaux = zeros(N)
            end
            if node[1] == 2.0
                logitaux = 1.0-logitaux
                dlogitaux = -dlogitaux
            end
            [dlogit[j] = logit[j]*dlogitaux[j] + dlogit[j]*logitaux[j] for j in 1:N]
            [logit[j] = logit[j]*logitaux[j] for j in 1:N]
            parent = floor(Int,node[5])
        end
        logimat[:,i] = dlogit
    end
    fitted::Vector{Float64} = logimat*terminal[:,2]
    return fitted
end

function fastlm(x::Matrix{Float64},y::Vector{Float64})::Vector{Float64}
    K::Int64 = size(x,2)
    N::Int64 = size(x,1)
    C::Matrix{Float64} = Matrix{Float64}(K,K)
    D::Matrix{Float64} = Matrix{Float64}(K,N)
    b::Vector{Float64} = Vector{Float64}(K)
    ssr::Float64 = 0
    res::Vector{Float64} = Vector{Float64}(K+1)
    fit::Vector{Float64} = Vector{Float64}(N)
    gemm!('T','N',1.0,x,x, 0.0, C)
    gemm!('N','T',1.0,inv(C),x, 0.0, D)
    gemv!('N',1.0,D,y, 0.0, b)
    gemv!('N',1.0,x,b,0.0,fit)
    ssr = sum( [(y[i]-fit[i])^2 for i in 1:N])
    res[1] = ssr
    res[2:end] = b
    return res
end

function eval_tree(x::Matrix{Float64},tree::Matrix{Float64})::Vector{Float64}
    N::Int64 = size(x)[1]
    terminalnodes::Vector{Int64} = find(x->(x==1),tree[:,6])
    terminal::Matrix{Float64} = tree[terminalnodes,:]
    Nt::Int64 = length(terminalnodes)
    logimat::Matrix{Float64} = Matrix{Float64}(N,Nt)
    node::Vector{Float64} = Vector{Float64}(9)
    logit::Vector{Float64} = Vector{Float64}(N)
    x_split::Vector{Float64} = Vector{Float64}(N)
    variable::Int64 = 1::Int64
    parent::Int64 = 1::Int64
    logitaux::Vector{Float64}=Vector{Float64}(N)
    for i in 1:Nt
        node=terminal[i,:]
        variable = floor(Int,node[7])
        x_split = x[:,variable]
        [logit[j] = 1/(1+exp(-node[4]*(x_split[j]-node[3]))) for j in 1:N]
        if node[1] == 2.0
            logit = 1.0-logit
        end
        parent = floor(Int,node[5])

        while parent!=0
            node = tree[parent,:]
            variable = floor(Int,node[7])
            x_split = x[:,variable]
            [logitaux[j] = 1/(1+exp(-node[4]*(x_split[j]-node[3]))) for j in 1:N]
            if node[1] == 2.0
                logitaux = 1.0-logitaux
            end
            [logit[j] = logit[j]*logitaux[j] for j in 1:N]
            parent = floor(Int,node[5])
        end
        logimat[:,i] = logit
    end
    fitted::Vector{Float64} = logimat*terminal[:,2]

    return fitted
end

function grow_tree(x::Matrix{Float64},y::Vector{Float64},p::Float64,d_max::Int64,gamma::Vector{Float64},node_obs::Int64)::Tuple{Array{Float64,2},Array{Float64,1}}

    bf::Int64 = 0
    K::Int64 = size(x,2)
    N::Int64 = length(y)
    (tree::Matrix{Float64}, pmat::Matrix{Float64}) = grow0(x,y,p,N,K,gamma,node_obs)::Tuple{Array{Float64,2},Array{Float64,2}}
    iter::Int64 = 1
    while iter <= d_max
       (tree, pmat, bf, iter) = grow(x,y,p,N,K,gamma,node_obs,tree,pmat,d_max,iter,bf)::Tuple{Array{Float64,2},Array{Float64,2},Int64,Int64}
    end

    middlenodes::Vector{Int64} = find(x->(isnan(x)),pmat[1,:])
    [pmat[i,j] = 0 for i in 1:N, j in middlenodes]
    fitted::Vector{Float64} = pmat*tree[:,2]

    return tree, fitted
end

function grow(x::Matrix{Float64},y::Vector{Float64},p::Float64,N::Int64,K::Int64,gamma::Vector{Float64},node_obs::Int64,tree::Matrix{Float64},pmat::Matrix{Float64},d_max::Int64,iter::Int64,bf::Int64)::Tuple{Array{Float64,2},Array{Float64,2},Int64,Int64}
    terminals::Vector{Int64} = find(x->(x==1),tree[:,6])
    variables::Vector{Int64} = StatsBase.sample(1:K,Int64(round(p*K)),replace = false)
    test::Array{Tuple{Int64,Int64},2} = collect(Base.product(variables,terminals))
    gammai::Float64 = gamma[rand(1:length(gamma),1)[1]]
    fit::Dict{Float64,Vector{Float64}} = Dict{Float64,Vector{Float64}}()
    sizehint!(fit,length(test))

    terminal::Int64 = 0
    variable::Int64 = 0
    Nt::Int64 = min(20,N)
    deep::Int64 = 0
    xtest0::Vector{Float64} = Vector{Float64}(N)
    xtest::Vector{Float64} = Vector{Float64}(Nt)
    pr::Vector{Float64} = Vector{Float64}(N)
    gammascale::Float64 = 0.0
    ssr::Matrix{Float64} =  Matrix{Float64}(Nt,3+iter)
    bestres::Int64 = 1
    res0::Vector{Float64} = Vector{Float64}(5+iter)
    middlenodes::Vector{Int64} =  find(x->(isnan(x)),pmat[1,:])
    for i in 1:length(test)
        terminal = test[i][2]
        variable = test[i][1]
        deep = tree[terminal,9]
        xtest0 = x[:,variable]
        pr = pmat[:,terminal]
        xtest =  StatsBase.wsample(xtest0, pr+0.01, min(20,N),replace = false)
        gammascale = max(std(xtest0),0.1)
        for j in 1:length(xtest)
            ssr[j,:] = nvt(xtest0,y,xtest[j],gammai/gammascale,node_obs,pmat,deep+1,terminal,middlenodes)::Vector{Float64}
        end
        bestres = find(x->(x==minimum(ssr[:,1])),ssr[:,1])[1]
        res = vcat(xtest[bestres],gammai/gammascale,ssr[bestres,:])
        fit[i] = res
    end

    # = 3 e o best val = #
    fits::Vector{Float64} = [fit[i][3] for i in 1:length(fit)]
    best::Int64 = find(x->(x==minimum(fits)),fits)[1]
    node::Vector{Float64} = fit[best]

    if bf==5
        iter = d_max +1
        return tree, pmat, bf, iter
    end
    if node[3]==Inf
        bf = bf + 1
        return tree, pmat, bf, iter
    end

    bterminal::Int64 = test[best][2]

    ### 1-side, 2-b, 3-c0, 4-gamma, 5-parent, 6-terminal,7-variable, 8-id, 9-deep
    nodes::Matrix{Float64}=Matrix{Float64}(2,9)
    nodes[1,:] = [1, node[4], node[1], node[2], tree[bterminal,8], 1, test[best][1], size(tree,1)+1, tree[bterminal,9]+1]
    nodes[2,:] = [2, node[5], node[1], node[2], tree[bterminal,8], 1, test[best][1], size(tree,1)+2, tree[bterminal,9]+1]

    # adjust old tree terminal nodes and updated bs
    tree[bterminal,6] = 0
    rem::Vector{Int64} = find(x->x==0,tree[:,6])
    terminalst::Vector{Bool} = trues(size(tree,1))
    terminalst[rem] = false
    tree[terminalst,2] = node[6:end]

    tree = vcat(tree,nodes)

    xb::Vector{Float64} = x[:,test[best][1]]
    probs::Vector{Float64}=Vector{Float64}(N)
    [probs[i] = 1/(1+exp(-nodes[1,4]*(xb[i]-nodes[1,3]))) for i in 1:length(xb)]
    pmataux::Matrix{Float64} = hcat(probs,1-probs)
    [pmataux[i,j] = pmataux[i,j]*pmat[i,bterminal] for i in 1:length(xb), j in 1:2]

    pmat::Matrix{Float64} = hcat(pmat,pmataux)::Matrix{Float64}
    [pmat[i,bterminal] = NaN for i in 1:N]
    iter = iter + 1
    return tree, pmat, bf, iter
end

function grow0(x::Matrix{Float64},y::Vector{Float64},p::Float64,N::Int64,K::Int64,gamma::Vector{Float64},node_obs::Int64)::Tuple{Array{Float64,2},Array{Float64,2}}
    variables::Vector{Int64} = StatsBase.sample(1:K,Int64(round(p*K)),replace = false)
    gammai::Float64 = gamma[rand(1:length(gamma),1)[1]]
    fit0::Dict{Float64,Vector{Float64}} = Dict{Float64,Vector{Float64}}()
    sizehint!(fit0,length(variables))
    Nt::Int64 = min(20,N)
    ssr::Matrix{Float64} =  Matrix{Float64}(Nt,3)
    xtest0::Vector{Float64} = Vector{Float64}(N)
    xtest::Vector{Float64} = Vector{Float64}(Nt)
    gammascale::Float64 = 1.0::Float64
    bestres::Int64 = 1::Int64
    res0::Vector{Float64} = Vector{Float64}(5)

    for i in 1:length(variables)
        xtest0 = x[:,variables[i]]
        xtest = StatsBase.sample(xtest0, Nt,replace = false)
        gammascale = max(std(xtest0),0.1)
        for j in 1:length(xtest)
            ssr[j,:] = invt(xtest0,y,xtest[j],gammai/gammascale,node_obs)
        end

        bestres= find(x->(x==minimum(ssr[:,1])),ssr[:,1])[1]
        res0 = vcat(xtest[bestres],gammai/gammascale,ssr[bestres,:])
        fit0[i] = res0
    end

    # = 3 e o best val = #
    fit0s::Vector{Float64} = [fit0[i][3] for i in 1:length(fit0)]
    best0::Int64 = find(x->(x==minimum(fit0s)),fit0s)[1]
    node0::Vector{Float64} = fit0[best0]

    ### 1-side, 2-b, 3-c0, 4-gamma, 5-parent, 6-terminal,7-variable, 8-id, 9-deep
    node0left::Vector{Float64} = [1, node0[4], node0[1], node0[2], 0, 1, variables[best0], 1, 1]
    node0right::Vector{Float64} = [2, node0[5], node0[1], node0[2], 0, 1, variables[best0], 2, 1]

    #tree::Matrix{Float64} = (hcat(node0left,node0right))'
    tree::Matrix{Float64} = Matrix{Float64}(2,9)
    tree[1,:] = node0left
    tree[2,:] = node0right

    xb::Vector{Float64} = x[:,variables[best0]]
    probs::Vector{Float64}=Vector{Float64}(N)
    [probs[i] = 1/(1+exp(-node0left[4]*(xb[i]-node0left[3]))) for i in 1:length(xb)]
    pmat::Matrix{Float64} = hcat(probs,1-probs)

    return tree, pmat
end

function nvt(x::Vector{Float64},y::Vector{Float64},c0::Float64,gamma::Float64,node_obs::Int64,pmat::Matrix{Float64},deep::Int64,terminal::Int64,middlenodes::Vector{Int64})::Vector{Float64}
    N::Int64 = length(x)
    logit::Vector{Float64} = Vector{Float64}(N)
    [logit[i] = 1/(1+exp(-gamma*(x[i]-c0))) for i in 1:N]

    terminal_logit::Vector{Float64} = pmat[:,terminal]
    bs::Matrix{Float64} = Matrix{Float64}(N,2)
    for i in 1:N
        bs[i,1] = logit[i] .* terminal_logit[i]
        bs[i,2] = (1-logit[i]) .* terminal_logit[i]
    end
    l0::Float64 = sum(sign.(bs[:,1]-0.5^deep)+1)/2
    l1::Float64 = sum(sign.(bs[:,2]-0.5^deep)+1)/2

    rem::Vector{Int64} = vcat(terminal,middlenodes)
    Bools::Vector{Bool} = trues(size(pmat,2))
    Bools[rem] = false
    X::Matrix{Float64} = hcat(bs,pmat[:,Bools])

    if l0 < node_obs || l1 < node_obs
        e::Vector{Float64} = fill(Inf64,size(X,2)+1)
        return e
    end

    res::Vector{Float64} = try
        fastlm(X,y)::Vector{Float64}
    catch
        fill(Inf64,size(X,2)+1)::Vector{Float64}
    end

    return res
end

function invt(x::Vector{Float64},y::Vector{Float64},c0::Float64,gamma::Float64,node_obs::Int64)::Vector{Float64}
    N::Int64 = length(x)
    b0 = Vector{Float64}(N)
    [b0[i] = 1/(1+exp(-gamma*(x[i]-c0))) for i in 1:N]
    b1::Vector{Float64} = 1 - b0

    l0::Float64 =  sum(sign.(b0-0.5)+1)/2
    l1::Float64 =  sum(sign.(b1-0.5)+1)/2
    #l0 = length(find(x->(x>=0.5),b0))
    #l1 = length(find(x->(x>=0.5),b1))

    X = Matrix{Float64}(N,2)
    X[:,1] = b0
    X[:,2] = b1

    if l0 < node_obs || l1 < node_obs
        e::Vector{Float64} = fill(Inf64,size(X,2)+1)
        return e
    end

    res::Vector{Float64} = try
        fastlm(X,y)
    catch
        fill(Inf64,size(X,2)+1)
    end

    return res
end

function fastlm1(x::Matrix{Float64},y::Vector{Float64})::Vector{Float64}
    bhat::Vector{Float64} = inv(x'*x)*x'*y
    error::Vector{Float64} = y-x*bhat
    ssr::Float64 = sum( [error[i]^2 for i in 1:length(y)])
    res::Vector{Float64} = vcat(ssr,bhat)
    return res
end

max_stars_repo_name: gabrielrvsc/BooSTjl

max_issues_repo_name: gabrielrvsc/BooSTjl

max_forks_repo_name: gabrielrvsc/BooSTjl


text:
# # The Category of Elements
# A very useful construction in applied category theory is the Category of Elements, which is also called the Grothendieck construction. This is a very general technique in category theory, but we will look at how you can use it to explain why graphs are so useful in computer science. We have already seen that C-Sets are a model of relational databases that can be used to store data as a collection of interlocking tables. Relational databases are the bread and butter of the computing industry. Every company on earth uses software that is backed by a relational database. Most data that is not stored in a relational DB is often stored in some kind of graph data structure. This sketch will show how these approaches are interchangeable via the category of elements, which associates to every database instance a graph and a graph homomorphism into the schema of the graph.
using Catlab, Catlab.CategoricalAlgebra, Catlab.Graphs, Catlab.Graphics
using Catlab.Graphics.Graphviz
using Colors

# Let's tell Catlab how to draw categories of elements.
function graph(el::Elements)
  F = FinFunctor(
    Dict(:V => :El, :E => :Arr),
    Dict(:src => :src, :tgt => :tgt),
    BasicGraphs.TheoryGraph, CatElements.ThElements
  )
  ΔF = DeltaMigration(F, Elements{Symbol}, BasicGraphs.Graph)
  return ΔF(el)
end

safecolors(start, stop, len) = if len > 1
  return hex.(range(start, stop=stop, length=len))
else
  return [hex.(start)]
end

GraphvizGraphs.to_graphviz(f::Elements; kw...) =
  to_graphviz(GraphvizGraphs.to_graphviz_property_graph(f; kw...))

function GraphvizGraphs.to_graphviz_property_graph(f::Elements; kw...)
  pg = GraphvizGraphs.to_graphviz_property_graph(graph(f); kw...)
  vcolors = safecolors(colorant"#0021A5", colorant"#FA4616", nparts(f, :Ob))
  ecolors = safecolors(colorant"#6C9AC3", colorant"#E28F41", nparts(f, :Hom))
  for v in parts(f, :El)
    fv = f[v, :πₑ]
    set_vprops!(pg, v, Dict(:color => "#$(vcolors[fv])", :label=>"$v:$(f[v,[:πₑ, :nameo]])"))
  end
  for e in parts(f, :Arr)
    fe = f[e, :πₐ]
    set_eprops!(pg, e, Dict(:color => "#$(ecolors[fe])"))
  end
  pg
end

draw(g) = to_graphviz(g, node_labels=true, edge_labels=true, prog="neato")

# ## The simplest schema
# First we will look at discrete dynamical systems. The set S is our state space and the funct nxt associates to every state, the next state in the system. This is a deterministic dynamical system with finitely many states and discrete time.
@present TheoryDDS(FreeSchema) begin
  S::Ob
  nxt::Hom(S, S)
end

@acset_type DDS(TheoryDDS, index=[:nxt])

fₓ = @acset DDS begin
  S = 3
  nxt = [2,3,1]
end

# Now if you want to draw a DDS, you might think to draw each state as a vertex and the next relationship as an arrow from one state to the next state. If your intuition selected that representation, you have already discovered a special case of the category of elements

elᶠ = elements(fₓ)
draw(graph(elᶠ))

# As you are quite far in your category theoretic training, you are surely suspicious of the idea that you would construct the category of elements of a DDS and get a graph. We have a categorical construction, so your intuition should be saying "there should be some kind of morphism that goes with this object". And there is, the Elements of a C-Set X is a graph whose vertices are rows in the tables of X and whose edges are the foreign key relationships in X, along with a graph homomorphism into the schema C. We can draw this using the same color coding convention we used to draw graph homomorphisms.
draw(elᶠ)

# We can scale up our DDS drawing too.

Fₓ = @acset DDS begin
  S = 7
  nxt = [2,3,1, 1,7,7,4]
end
draw(elements(Fₓ))

# A category of elements derived from a C-Set is stored as a C-Set of on a different schema. You can see that it is the data of a graph homomorphism where the codomain graph has vertex and edge labels. The two projections πₑ and πₐ are the components of a natural transformation between graphs viewed as functors into Set. The names are attributes usually of type Symbol.
to_graphviz(ThElements)

# In the case of a DDS, we have only one object and one morphism in the schema. Since the graph with one edge and one vertices is terminal in *Graph*, there is only one vertex and edge color being used.

# ## The Elements of a Graph are its Vertices and Edges
# In what might appear as primordial ooze, we can examine the category of elements of a graph. We will look at the commuting triangle graph. Notice how there are 3 vertices and three edges with each edge incident to 2 vertices.
g = @acset BasicGraphs.Graph begin
  V = 3
  E = 3
  src = [1,2,1]
  tgt = [2,3,3]
end
draw(g)

# The category of elements has 6 vertices and 6 edges. The 6 vertices come from 3 vertices of type V and 3 vertices of type E. The 6 edges in this graph are 3 src relationships and 3 target relationships.
elᵍ = elements(g)
draw(elᵍ)

g = @acset BasicGraphs.Graph begin
  V = 6
  E = 7
  src = [1,2,1,3,5,6,4]
  tgt = [2,3,3,5,6,3,4]
end
draw(g)
draw(elements(g))

# Notice that the 3 vertices of type E are interspersed between the 3 vertices of type V and vice versa. You can see that the normal visual syntax for a graph is more compact for the special case of graphs. However the category of elements works for any schema C.

# ## Generality of the construction
# Discrete Dynamical Systems and Graphs are clearly data structures from mathematics and so it would make sense that they have a clean representation in the categorical language. But how about a database schema that comes not from mathematics, but from software engineering. We turn to everyone's favorite database example, the HR database at a fictitious company.
@present ThCompany(FreeSchema) begin
  (P, D, S)::Ob # Person, Department, Salary
  worksin::Hom(P, D)    # Every Person works in a Department
  makes::Hom(P, S)      # Every Person makes a Salary
  reportsto::Hom(P, P)  # Every Person reports to a Person
  managedby::Hom(D, P)  # Every Department is managed by a Person
  leq::Hom(S,S)         # Salaries are a finite total order
end

@acset_type Company(ThCompany, index=[])

# We can draw a company that has 4 people, 2 departments, and 3 distinct salaries.
cmpy = @acset Company begin
  P = 4
  D = 2
  S = 3
  worksin = [1,1,2,2]
  makes = [3,1,2,1]
  reportsto = [1, 1,1,3]
  managedby = [1,3]
  leq = [2,3,3]
end

# A visualization of the elements of this functor is commonly described as a knowledge graph.
elᶜ = elements(cmpy)
draw(elᶜ)

# The visualization of the category of elements isn't the most compact representation of data, but it is always an option. Once you know how to draw a graph homomorphism as a color-coded graph, you know how to draw any C-Set. The vertices and edges of the domain graph are the rows and fields of the C-Set, and the colors are from the vertices and edges of the schema C, viewed as a graph presenting the category. One thing you lose is the ability to represent path equations in C.

# ## Reconstructing the C-Set from its category of elemnts
# When looking at the graphical representation of the category of elements, you can see that the vertices are numbered sequentially within their vertex type. This order can be computed with the incident function in Catlab. There are inclusions of all the tables in the database instance into the global vertex set of the knowledge graph representation.
incident(elᶜ, :, :πₑ)

# The same goes for the edges (arrows) in the category of elements.
incident(elᶜ, :, :πₐ)

# From this information, you can reassemble the knowledge graph into a database. You use the codomain graph as the schema and then these two projections πₑ and πₐ to recover the instance data.

cset(T::Type, el::Elements) = begin
  X = T()
  ℓ = Dict{Int, Int}() # this map reverses the inclusions converting global element numbers into per-table numbers
  for ob in parts(el, :Ob)
    obname = el[ob, :nameo]
    eltsob = incident(el, ob, :πₑ)
    pts = add_parts!(X, obname, length(eltsob))
    map(zip(pts, eltsob)) do (i,x)
      ℓ[x] = i
    end
  end
  for h in parts(el, :Hom)
    nameh = el[h, :nameh]
    arrₕ = incident(el, h, :πₐ)
    doms = map(arrₕ) do e
      ℓ[el[e, :src]]
    end
    codoms = map(arrₕ) do e
      ℓ[el[e, :tgt]]
    end
    set_subpart!(X, doms, nameh, codoms)
  end
  return X
end

# We can apply this algorithm to our company knowledge graph to recover our original DB.
cset(Company, elᶜ)
cmpy == cset(Company, elᶜ)

# Hopefully you see one reason why graphs are ubiquitous in computer science. Every relational database is isomorphic to a typed graph. You can generate this graph from the database and recover the database from the typed graph. So while knowledge graphs are great for ingesting data when you aren't sure about the structure and useful for some types of computation that involve long paths. You can always build your software on relational databases knowing that you can easily get access to all the graph algorithms by applying the Grothendieck construction.

# ## The Slice-Elements Transform
# An amazing fact about presheaf toposes is that they are closed under taking slices. In the graphs section of this documentation, you can find a description of bipartite and k-partite graphs as a morphisms into a clique. That definition is very mathematically pleasing because it gives you a category of partitioned graphs that are derived from commuting triangles in Graph. However, for application oriented practitioners, the definition of a bipartite graph as "a graph with two sets of vertices, where all the edges go between the groups with no edges within either group" is probably more explicit. For example a classic way to get a bipartite graph would be to look at the graph of authors and papers that those authors wrote. People write papers, people do not write people and papers do not write papers so the authorship graph is bipartite. These two equivalent definitions of a bipartite graph are related via an isomorphism you can find on the [nlab](https://ncatlab.org/nlab/show/category+of+presheaves#RelWithOvercategories). It shows that a slice category of [C,Set]/X is isomorphic to [El(X), Set] which is a cateogory of presheaves on a different schema. Catlab knows how to use this idea to turn a category of elements into a schema for a new category of presheaves. The two directions of the isomorphism are not yet implemented.

e = @acset BasicGraphs.Graph begin
  V = 2
  E = 1
  src = 1
  tgt = 2
end

draw(elements(e))
ThBipartite = CatElements.presentation(elements(e))[1]
to_graphviz(ThBipartite)
@acset_type BipartiteGraph(ThBipartite)
b = @acset BipartiteGraph begin
  V_1 = 3
  V_2 = 2
  E_1 = 5
  src_E_1 = [1,2,3,3,2]
  tgt_E_1 = [1,2,2,1,1]
end
draw(elements(b))

# As a very advanced exercise, you could try to implement one or both directions of the isomorphism above.


text:
# # An unstable Bickley jet in Shallow Water model
#
# This example shows how to use `Oceananigans.ShallowWaterModel` to simulate
# the evolution of an unstable, geostrophically balanced, Bickley jet.
# The model solves the governing equations for the shallow water model in
# conservative form.  The geometry is that of a periodic channel
# in the ``x``-direction with a flat bottom and a free-surface. The initial
# conditions are that of a Bickley jet with small-amplitude perturbations.
# The interested reader can see ["The nonlinear evolution of barotropically unstable jets," J. Phys. Oceanogr. (2003)](https://doi.org/10.1175/1520-0485(2003)033<2173:TNEOBU>2.0.CO;2)
# for more details on this specific problem.
#
# Unlike the other models, the fields that are simulated are the mass transports,
# ``uh`` and ``vh`` in the ``x`` and ``y`` directions, respectively,
# and the height ``h``.  Note that the velocities ``u`` and ``v`` are not state
# variables but can be easily computed when needed, e.g., via `u = uh / h`.
#
# ## Install dependencies
#
# First we make sure that we have all of the packages that are required to
# run the simulation.
# ```julia
# using Pkg
# pkg"add Oceananigans, NCDatasets, Plots, Printf, Polynomials"
# ```

using Oceananigans
using Oceananigans.Models: ShallowWaterModel

# ## Two-dimensional domain
#
# The shallow water model is a two-dimensional model and thus the number of vertical
# points `Nz` must be set to one.  Note that ``L_z`` is the mean depth of the fluid.

Lx, Ly, Lz = 2π, 20, 1
Nx, Ny = 128, 128

grid = RegularRectilinearGrid(size = (Nx, Ny),
                              x = (0, Lx), y = (-Ly/2, Ly/2),
                              topology = (Periodic, Bounded, Flat))

# ## Physical parameters
#
# This is a toy problem and we choose the parameters so the jet idealizes a relatively narrow
# mesoscale jet. The physical parameters are
#
#   * ``f``: Coriolis parameter
#   * ``g``: Acceleration due to gravity
#   * ``U``: Maximum jet speed
#   * ``\Delta \eta``: Maximum free-surface deformation as dictated by geostrophy

 const f = 1
 const g = 9.8
 const U = 1.0
const Δη = f * U / g
nothing # hide

# ## Building a `ShallowWaterModel`
#
# We use `grid`, `coriolis` and `gravitational_acceleration` to build the model.
# Furthermore, we specify `RungeKutta3` for time-stepping and `WENO5` for advection.

model = ShallowWaterModel(
    timestepper=:RungeKutta3,
    advection=WENO5(),
    grid=grid,
    gravitational_acceleration=g,
    coriolis=FPlane(f=f),
    )

# ## Background state and perturbation
#
# The background velocity ``ū`` and free-surface ``η̄`` are chosen to represent a
# geostrophically balanced Bickely jet with maximum speed of ``U`` and maximum
# free-surface deformation of ``Δη``, i.e.,
#
# ```math
# \begin{align}
# η̄(y) & = - Δη \tanh(y) , \\
# ū(y) & = U \mathrm{sech}^2(y) .
# \end{align}
# ```
#
# The total height of the fluid is ``h = L_z + \eta``. Linear stability theory predicts that
# for the parameters we consider here, the growth rate for the most unstable mode that fits
# our domain is approximately ``0.139``.
#
# We also specify `ω̄` as the vorticity of the background state,
# ``ω̄ = - ∂_y ū = 2 U \mathrm{sech}^2(y) \tanh(y)``.

h̄(x, y, z) = Lz - Δη * tanh(y)
ū(x, y, z) = U * sech(y)^2
ω̄(x, y, z) = 2 * U * sech(y)^2 * tanh(y)
nothing # hide

# The initial conditions include a small-amplitude perturbation that decays away from the
# center of the jet.

 small_amplitude = 1e-4

 uⁱ(x, y, z) = ū(x, y, z) + small_amplitude * exp(-y^2) * randn()
 hⁱ(x, y, z) = h̄(x, y, z)
uhⁱ(x, y, z) = uⁱ(x, y, z) * hⁱ(x, y, z)
nothing # hide

# We set the initial conditions for the zonal mass transport `uhⁱ` and the fluid height `hⁱ`.

set!(model, uh = uhⁱ, h = hⁱ)

# We compute the total vorticity and the perturbation vorticity.

uh, vh, h = model.solution

        u = ComputedField(uh / h)
        v = ComputedField(vh / h)
        ω = ComputedField(∂x(v) - ∂y(u))
   ω_pert = ComputedField(ω - ω̄)
nothing #hide

# ## Running a `Simulation`
#
# We pick the time-step so that we make sure we resolve the surface gravity waves, which
# propagate with speed of the order ``\sqrt{g L_z}``. That is, with `Δt = 1e-2` we ensure
# that `` \sqrt{g L_z} Δt / Δx,  \sqrt{g L_z} Δt / Δy < 0.7``.

simulation = Simulation(model, Δt = 1e-2, stop_time = 150)

# ## Prepare output files
#
# Define a function to compute the norm of the perturbation on the cross channel velocity.
# We obtain the `norm` function from `LinearAlgebra`.

using LinearAlgebra: norm

function perturbation_norm(model)
    compute!(v)
    return norm(interior(v))
end
nothing # hide

# Choose the two fields to be output to be the total and perturbation vorticity.

outputs = (ω_total = ω, ω_pert = ω_pert)

# Build the `output_writer` for the two-dimensional fields to be output.
# Output every `t = 1.0`.

simulation.output_writers[:fields] =
    NetCDFOutputWriter(
        model,
        (ω = ω, ω_pert = ω_pert),
          filepath = joinpath(@__DIR__, "shallow_water_Bickley_jet.nc"),
          schedule = TimeInterval(1.0),
              mode = "c")

# Build the `output_writer` for the growth rate, which is a scalar field.
# Output every time step.

simulation.output_writers[:growth] =
    NetCDFOutputWriter(
        model,
        (perturbation_norm = perturbation_norm,),
          filepath = joinpath(@__DIR__, "perturbation_norm_shallow_water.nc"),
          schedule = IterationInterval(1),
        dimensions = (perturbation_norm=(),),
              mode = "c")

# And finally run the simulation.

run!(simulation)

# ## Visualize the results

# Load required packages to read output and plot.

using NCDatasets, Plots, Printf
nothing # hide

# Define the coordinates for plotting.

x, y = xnodes(ω), ynodes(ω)
nothing # hide

# Define keyword arguments for plotting the contours.

kwargs = (
         xlabel = "x",
         ylabel = "y",
         aspect = 1,
           fill = true,
         levels = 20,
      linewidth = 0,
          color = :balance,
       colorbar = true,
           ylim = (-Ly/2, Ly/2),
           xlim = (0, Lx)
)
nothing # hide

# Read in the `output_writer` for the two-dimensional fields and then create an animation
# showing both the total and perturbation vorticities.

ds = NCDataset(simulation.output_writers[:fields].filepath, "r")

iterations = keys(ds["time"])

anim = @animate for (iter, t) in enumerate(ds["time"])
     ω = ds["ω"][:, :, 1, iter]
    ωp = ds["ω_pert"][:, :, 1, iter]

    ωp_max = maximum(abs, ωp)

     plot_ω = contour(x, y, ω',
                       clim = (-1, 1),
                      title = @sprintf("Total vorticity, ω, at t = %.1f", t); kwargs...)

    plot_ωp = contour(x, y, ωp',
                       clim = (-ωp_max, ωp_max),
                      title = @sprintf("Perturbation vorticity, ω - ω̄, at t = %.1f", t); kwargs...)

    plot(plot_ω, plot_ωp, layout = (1, 2), size = (800, 440))
end

close(ds)

mp4(anim, "Bickley_Jet_ShallowWater.mp4", fps=15)

# Read in the `output_writer` for the scalar field (the norm of ``v``-velocity).

ds2 = NCDataset(simulation.output_writers[:growth].filepath, "r")

iterations = keys(ds2["time"])

     t = ds2["time"][:]
norm_v = ds2["perturbation_norm"][:]

close(ds2)
nothing # hide

# We import the `fit` function from `Polynomials.jl` to compute the best-fit slope of the
# perturbation norm on a logarithmic plot. This slope corresponds to the growth rate.

using Polynomials: fit

I = 6000:7000

degree = 1
linear_fit_polynomial = fit(t[I], log.(norm_v[I]), degree, var = :t)

# We can get the coefficient of the ``n``-th power from the fitted polynomial by using `n`
# as an index, e.g.,

constant, slope = linear_fit_polynomial[0], linear_fit_polynomial[1]

# We then use the computed linear fit coefficients to construct the best fit and plot it
# together with the time-series for the perturbation norm for comparison.

best_fit = @. exp(constant + slope * t)

plot(t, norm_v,
        yaxis = :log,
        ylims = (1e-3, 30),
           lw = 4,
        label = "norm(v)",
       xlabel = "time",
       ylabel = "norm(v)",
        title = "growth of perturbation norm",
       legend = :bottomright)

plot!(t[I], 2 * best_fit[I], # factor 2 offsets fit from curve for better visualization
           lw = 4,
        label = "best fit")

# The slope of the best-fit curve on a logarithmic scale approximates the rate at which instability
# grows in the simulation. Let's see how this compares with the theoretical growth rate.

println("Numerical growth rate is approximated to be ", round(slope, digits=3), ",\n",
        "which is very close to the theoretical value of 0.139.")


text:
"""
    ERA5Region

Structure that imports relevant [GeoRegion](https://github.com/JuliaClimate/GeoRegions.jl) properties used in the handling of the gridded ERA5 datasets.

All `ERA5Region` Types contain the following fields:
- `geo`   : The `GeoRegion` containing the geographical information
- `geoID` : The ID used to specify the `GeoRegion`
- `gres`  : The resolution of the gridded data to be downloaded/analysed
- `fstr`  : String, for specification of folder and file name
- `isglb` : A Bool, true if spans the globe, false if no
"""
struct ERA5Region{ST<:AbstractString, FT<:Real}
    geo   :: GeoRegion
    geoID :: ST
    gres  :: FT
    fstr  :: ST
    isglb :: Bool
end

"""
    ERA5Region(
        geo  :: GeoRegion;
        gres :: Real = 0,
        ST = String,
        FT = Float64
    ) -> egeo :: ERA5Region

Argument
========

- `geo`  : A `GeoRegion` structure type

Keyword Argument
================

- `gres` : The spatial resolution that ERA5 reanalysis data will be downloaded/analyzed, and 360 must be a multiple of `gres`
"""
function ERA5Region(
    geo  :: GeoRegion;
    gres :: Real = 0,
    ST = String,
    FT = Float64
)

    @info "$(modulelog()) - Creating an ERA5Region based on the GeoRegion \"$(geo.regID)\""
    gres = regionstep(geo.regID,gres)
    if geo.regID == "GLB"; isglb = true; else; isglb = false end

    return ERA5Region{ST,FT}(
        geo, geo.regID, gres,
        "$(geo.regID)x$(@sprintf("%.2f",gres))", isglb
    )

end

function regionstep(
    geoID :: AbstractString,
    gres  :: Real = 0
)

    @debug "$(modulelog()) - Determining spacing between grid points in the GeoRegion ..."
    if gres == 0
        @info "$(modulelog()) - No grid resolution specified, defaulting to the module default (1.0º for global GeoRegion, 0.25º for all others)"
        if geoID == "GLB";
              gres = 1.0;
        else; gres = 0.25;
        end
    else
        if !checkegrid(gres)
            error("$(modulelog()) - The grid resolution $(gres)º is not valid as it does not divide 360º without remainder")
        end
    end

    return gres

end

function checkegrid(gres::Real)

    if rem(360,gres) == 0
          return true
    else; return false
    end

end

function show(io::IO, ereg::ERA5Region)
    geo = ereg.geo

    if typeof(geo) <: PolyRegion
        print(
            io,
            "The ERA5Region wrapper for the \"$(ereg.geoID)\" GeoRegion has the following properties:\n",
            "    Region ID (regID) : ", ereg.geoID, '\n',
            "    Name       (name) : ", ereg.geo.name,  '\n',
            "    Resolution (gres) : ", ereg.gres,  '\n',
            "    Bounds  (N,S,E,W) : ",[geo.N,geo.S,geo.E,geo.W], '\n',
            "    Shape     (shape) : ", geo.shape, '\n',
            "        (is180,is360) : ",(geo.is180,geo.is360),"\n",
        )
    else
        print(
            io,
            "The ERA5Region wrapper for the \"$(ereg.geoID)\" GeoRegion has the following properties:\n",
            "    Region ID (regID) : ", ereg.geoID, '\n',
            "    Name       (name) : ", ereg.geo.name,  '\n',
            "    Resolution (gres) : ", ereg.gres,  '\n',
            "    Bounds  (N,S,E,W) : ",[geo.N,geo.S,geo.E,geo.W], '\n',
            "        (is180,is360) : ",(geo.is180,geo.is360),"\n",
        )
    end

end


text:
using Krypto
using Base.Test

#include("rsa.jl")
#include("rlwe.jl")
include("salsa20.jl")


text:
using CryptoMiniSat

#=
Basic sanity check that we can interface properly
with the library.
=#
function test_nvars()
    solver = SATSolver()

    new_vars(solver, 7)

    @assert nvars(solver) == 7
end

function test_satisfiable()
    # Example taken from:
    # https://github.com/msoos/cryptominisat#rust-binding

    # Create new instance of SAT solver
    solver = SATSolver()

    # Set number of threads
    set_num_threads(solver, 2)

    new_vars(solver, 3)

    # Create clauses
    @assert add_clause(solver, [Lit(0)])
    @assert add_clause(solver, [~Lit(1)])
    @assert add_clause(solver, [~Lit(0), Lit(1), Lit(2)])

    # Solve system and verify
    sat, mdl = solve(solver)
    @assert sat == true
    @assert mdl == [true, false, true]

    # no need to free the SAT solver,
    # as it will be finalized with gc
end

function test_unsatisfiable()
    solver = SATSolver()

    new_vars(solver, 3)

    # Create clauses
    @assert add_clause(solver, [Lit(0)])
    @assert add_clause(solver, [~Lit(1)])
    @assert add_clause(solver, [~Lit(0), Lit(1), Lit(2)])

    # Solve system with a conflicting assumption.
    sat, mdl = solve(solver, [Lit(1), ~Lit(2)])
    @assert sat == false

    # Verify that the conflict has been identified.
    conflict = get_conflict(solver)
    @assert conflict == [~Lit(1)]
end

function test_xor_clause()
    solver = SATSolver()

    new_vars(solver, 3)

    # Create clauses
    @assert add_clause(solver, [Lit(0)])
    @assert add_xor_clause(solver, Unsigned[0, 1], true)
    @assert add_xor_clause(solver, Unsigned[0, 1, 2], true)

    # Solve system with a conflicting assumption.
    sat, mdl = solve(solver)
    @assert sat == true

    # Verify that the solution has been identified.
    @assert sat == true
    @assert mdl == [true, false, false]

    # Verify that this is in fact the only solution
    @assert add_clause(solver, [~Lit(0), Lit(1), Lit(2)]) == false
end

test_nvars()
test_satisfiable()
test_unsatisfiable()
test_xor_clause()


text:
module FreqTables
    using CategoricalArrays
    using DataFrames
    using NamedArrays

    include("freqtable.jl")

    export freqtable
end # module


text:
# Julia wrapper for header: tiff.h
# Automatically generated using Clang.jl

# Julia wrapper for header: tiffconf.h
# Automatically generated using Clang.jl

# Julia wrapper for header: tiffio.h
# Automatically generated using Clang.jl

function TIFFGetVersion()
    ccall((:TIFFGetVersion, libtiff), Cstring, ())
end

function TIFFFindCODEC(arg1)
    ccall((:TIFFFindCODEC, libtiff), Ptr{TIFFCodec}, (uint16,), arg1)
end

function TIFFRegisterCODEC(arg1, arg2, arg3)
    ccall((:TIFFRegisterCODEC, libtiff), Ptr{TIFFCodec}, (uint16, Cstring, TIFFInitMethod), arg1, arg2, arg3)
end

function TIFFUnRegisterCODEC(arg1)
    ccall((:TIFFUnRegisterCODEC, libtiff), Cvoid, (Ptr{TIFFCodec},), arg1)
end

function TIFFIsCODECConfigured(arg1)
    ccall((:TIFFIsCODECConfigured, libtiff), Cint, (uint16,), arg1)
end

function TIFFGetConfiguredCODECs()
    ccall((:TIFFGetConfiguredCODECs, libtiff), Ptr{TIFFCodec}, ())
end

function _TIFFmalloc(s)
    ccall((:_TIFFmalloc, libtiff), Ptr{Cvoid}, (tmsize_t,), s)
end

function _TIFFcalloc(nmemb, siz)
    ccall((:_TIFFcalloc, libtiff), Ptr{Cvoid}, (tmsize_t, tmsize_t), nmemb, siz)
end

function _TIFFrealloc(p, s)
    ccall((:_TIFFrealloc, libtiff), Ptr{Cvoid}, (Ptr{Cvoid}, tmsize_t), p, s)
end

function _TIFFmemset(p, v, c)
    ccall((:_TIFFmemset, libtiff), Cvoid, (Ptr{Cvoid}, Cint, tmsize_t), p, v, c)
end

function _TIFFmemcpy(d, s, c)
    ccall((:_TIFFmemcpy, libtiff), Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}, tmsize_t), d, s, c)
end

function _TIFFmemcmp(p1, p2, c)
    ccall((:_TIFFmemcmp, libtiff), Cint, (Ptr{Cvoid}, Ptr{Cvoid}, tmsize_t), p1, p2, c)
end

function _TIFFfree(p)
    ccall((:_TIFFfree, libtiff), Cvoid, (Ptr{Cvoid},), p)
end

function TIFFGetTagListCount(arg1)
    ccall((:TIFFGetTagListCount, libtiff), Cint, (Ptr{TIFF},), arg1)
end

function TIFFGetTagListEntry(arg1, tag_index)
    ccall((:TIFFGetTagListEntry, libtiff), uint32, (Ptr{TIFF}, Cint), arg1, tag_index)
end

function TIFFFindField(arg1, arg2, arg3)
    ccall((:TIFFFindField, libtiff), Ptr{TIFFField}, (Ptr{TIFF}, uint32, TIFFDataType), arg1, arg2, arg3)
end

function TIFFFieldWithTag(arg1, arg2)
    ccall((:TIFFFieldWithTag, libtiff), Ptr{TIFFField}, (Ptr{TIFF}, uint32), arg1, arg2)
end

function TIFFFieldWithName(arg1, arg2)
    ccall((:TIFFFieldWithName, libtiff), Ptr{TIFFField}, (Ptr{TIFF}, Cstring), arg1, arg2)
end

function TIFFFieldTag(arg1)
    ccall((:TIFFFieldTag, libtiff), uint32, (Ptr{TIFFField},), arg1)
end

function TIFFFieldName(arg1)
    ccall((:TIFFFieldName, libtiff), Cstring, (Ptr{TIFFField},), arg1)
end

function TIFFFieldDataType(arg1)
    ccall((:TIFFFieldDataType, libtiff), TIFFDataType, (Ptr{TIFFField},), arg1)
end

function TIFFFieldPassCount(arg1)
    ccall((:TIFFFieldPassCount, libtiff), Cint, (Ptr{TIFFField},), arg1)
end

function TIFFFieldReadCount(arg1)
    ccall((:TIFFFieldReadCount, libtiff), Cint, (Ptr{TIFFField},), arg1)
end

function TIFFFieldWriteCount(arg1)
    ccall((:TIFFFieldWriteCount, libtiff), Cint, (Ptr{TIFFField},), arg1)
end

function TIFFAccessTagMethods(arg1)
    ccall((:TIFFAccessTagMethods, libtiff), Ptr{TIFFTagMethods}, (Ptr{TIFF},), arg1)
end

function TIFFGetClientInfo(arg1, arg2)
    ccall((:TIFFGetClientInfo, libtiff), Ptr{Cvoid}, (Ptr{TIFF}, Cstring), arg1, arg2)
end

function TIFFSetClientInfo(arg1, arg2, arg3)
    ccall((:TIFFSetClientInfo, libtiff), Cvoid, (Ptr{TIFF}, Ptr{Cvoid}, Cstring), arg1, arg2, arg3)
end

function TIFFCleanup(tif)
    ccall((:TIFFCleanup, libtiff), Cvoid, (Ptr{TIFF},), tif)
end

function TIFFClose(tif)
    ccall((:TIFFClose, libtiff), Cvoid, (Ptr{TIFF},), tif)
end

function TIFFFlush(tif)
    ccall((:TIFFFlush, libtiff), Cint, (Ptr{TIFF},), tif)
end

function TIFFFlushData(tif)
    ccall((:TIFFFlushData, libtiff), Cint, (Ptr{TIFF},), tif)
end

function TIFFReadDirectory(tif)
    ccall((:TIFFReadDirectory, libtiff), Cint, (Ptr{TIFF},), tif)
end

function TIFFReadCustomDirectory(tif, diroff, infoarray)
    ccall((:TIFFReadCustomDirectory, libtiff), Cint, (Ptr{TIFF}, toff_t, Ptr{TIFFFieldArray}), tif, diroff, infoarray)
end

function TIFFReadEXIFDirectory(tif, diroff)
    ccall((:TIFFReadEXIFDirectory, libtiff), Cint, (Ptr{TIFF}, toff_t), tif, diroff)
end

function TIFFScanlineSize64(tif)
    ccall((:TIFFScanlineSize64, libtiff), uint64, (Ptr{TIFF},), tif)
end

function TIFFScanlineSize(tif)
    ccall((:TIFFScanlineSize, libtiff), tmsize_t, (Ptr{TIFF},), tif)
end

function TIFFRasterScanlineSize64(tif)
    ccall((:TIFFRasterScanlineSize64, libtiff), uint64, (Ptr{TIFF},), tif)
end

function TIFFRasterScanlineSize(tif)
    ccall((:TIFFRasterScanlineSize, libtiff), tmsize_t, (Ptr{TIFF},), tif)
end

function TIFFStripSize64(tif)
    ccall((:TIFFStripSize64, libtiff), uint64, (Ptr{TIFF},), tif)
end

function TIFFStripSize(tif)
    ccall((:TIFFStripSize, libtiff), tmsize_t, (Ptr{TIFF},), tif)
end

function TIFFRawStripSize64(tif, strip)
    ccall((:TIFFRawStripSize64, libtiff), uint64, (Ptr{TIFF}, uint32), tif, strip)
end

function TIFFRawStripSize(tif, strip)
    ccall((:TIFFRawStripSize, libtiff), tmsize_t, (Ptr{TIFF}, uint32), tif, strip)
end

function TIFFVStripSize64(tif, nrows)
    ccall((:TIFFVStripSize64, libtiff), uint64, (Ptr{TIFF}, uint32), tif, nrows)
end

function TIFFVStripSize(tif, nrows)
    ccall((:TIFFVStripSize, libtiff), tmsize_t, (Ptr{TIFF}, uint32), tif, nrows)
end

function TIFFTileRowSize64(tif)
    ccall((:TIFFTileRowSize64, libtiff), uint64, (Ptr{TIFF},), tif)
end

function TIFFTileRowSize(tif)
    ccall((:TIFFTileRowSize, libtiff), tmsize_t, (Ptr{TIFF},), tif)
end

function TIFFTileSize64(tif)
    ccall((:TIFFTileSize64, libtiff), uint64, (Ptr{TIFF},), tif)
end

function TIFFTileSize(tif)
    ccall((:TIFFTileSize, libtiff), tmsize_t, (Ptr{TIFF},), tif)
end

function TIFFVTileSize64(tif, nrows)
    ccall((:TIFFVTileSize64, libtiff), uint64, (Ptr{TIFF}, uint32), tif, nrows)
end

function TIFFVTileSize(tif, nrows)
    ccall((:TIFFVTileSize, libtiff), tmsize_t, (Ptr{TIFF}, uint32), tif, nrows)
end

function TIFFDefaultStripSize(tif, request)
    ccall((:TIFFDefaultStripSize, libtiff), uint32, (Ptr{TIFF}, uint32), tif, request)
end

function TIFFDefaultTileSize(arg1, arg2, arg3)
    ccall((:TIFFDefaultTileSize, libtiff), Cvoid, (Ptr{TIFF}, Ptr{uint32}, Ptr{uint32}), arg1, arg2, arg3)
end

function TIFFFileno(arg1)
    ccall((:TIFFFileno, libtiff), Cint, (Ptr{TIFF},), arg1)
end

function TIFFSetFileno(arg1, arg2)
    ccall((:TIFFSetFileno, libtiff), Cint, (Ptr{TIFF}, Cint), arg1, arg2)
end

function TIFFClientdata(arg1)
    ccall((:TIFFClientdata, libtiff), thandle_t, (Ptr{TIFF},), arg1)
end

function TIFFSetClientdata(arg1, arg2)
    ccall((:TIFFSetClientdata, libtiff), thandle_t, (Ptr{TIFF}, thandle_t), arg1, arg2)
end

function TIFFGetMode(arg1)
    ccall((:TIFFGetMode, libtiff), Cint, (Ptr{TIFF},), arg1)
end

function TIFFSetMode(arg1, arg2)
    ccall((:TIFFSetMode, libtiff), Cint, (Ptr{TIFF}, Cint), arg1, arg2)
end

function TIFFIsTiled(arg1)
    ccall((:TIFFIsTiled, libtiff), Cint, (Ptr{TIFF},), arg1)
end

function TIFFIsByteSwapped(arg1)
    ccall((:TIFFIsByteSwapped, libtiff), Cint, (Ptr{TIFF},), arg1)
end

function TIFFIsUpSampled(arg1)
    ccall((:TIFFIsUpSampled, libtiff), Cint, (Ptr{TIFF},), arg1)
end

function TIFFIsMSB2LSB(arg1)
    ccall((:TIFFIsMSB2LSB, libtiff), Cint, (Ptr{TIFF},), arg1)
end

function TIFFIsBigEndian(arg1)
    ccall((:TIFFIsBigEndian, libtiff), Cint, (Ptr{TIFF},), arg1)
end

function TIFFGetReadProc(arg1)
    ccall((:TIFFGetReadProc, libtiff), TIFFReadWriteProc, (Ptr{TIFF},), arg1)
end

function TIFFGetWriteProc(arg1)
    ccall((:TIFFGetWriteProc, libtiff), TIFFReadWriteProc, (Ptr{TIFF},), arg1)
end

function TIFFGetSeekProc(arg1)
    ccall((:TIFFGetSeekProc, libtiff), TIFFSeekProc, (Ptr{TIFF},), arg1)
end

function TIFFGetCloseProc(arg1)
    ccall((:TIFFGetCloseProc, libtiff), TIFFCloseProc, (Ptr{TIFF},), arg1)
end

function TIFFGetSizeProc(arg1)
    ccall((:TIFFGetSizeProc, libtiff), TIFFSizeProc, (Ptr{TIFF},), arg1)
end

function TIFFGetMapFileProc(arg1)
    ccall((:TIFFGetMapFileProc, libtiff), TIFFMapFileProc, (Ptr{TIFF},), arg1)
end

function TIFFGetUnmapFileProc(arg1)
    ccall((:TIFFGetUnmapFileProc, libtiff), TIFFUnmapFileProc, (Ptr{TIFF},), arg1)
end

function TIFFCurrentRow(arg1)
    ccall((:TIFFCurrentRow, libtiff), uint32, (Ptr{TIFF},), arg1)
end

function TIFFCurrentDirectory(arg1)
    ccall((:TIFFCurrentDirectory, libtiff), uint16, (Ptr{TIFF},), arg1)
end

function TIFFNumberOfDirectories(arg1)
    ccall((:TIFFNumberOfDirectories, libtiff), uint16, (Ptr{TIFF},), arg1)
end

function TIFFCurrentDirOffset(arg1)
    ccall((:TIFFCurrentDirOffset, libtiff), uint64, (Ptr{TIFF},), arg1)
end

function TIFFCurrentStrip(arg1)
    ccall((:TIFFCurrentStrip, libtiff), uint32, (Ptr{TIFF},), arg1)
end

function TIFFCurrentTile(tif)
    ccall((:TIFFCurrentTile, libtiff), uint32, (Ptr{TIFF},), tif)
end

function TIFFReadBufferSetup(tif, bp, size)
    ccall((:TIFFReadBufferSetup, libtiff), Cint, (Ptr{TIFF}, Ptr{Cvoid}, tmsize_t), tif, bp, size)
end

function TIFFWriteBufferSetup(tif, bp, size)
    ccall((:TIFFWriteBufferSetup, libtiff), Cint, (Ptr{TIFF}, Ptr{Cvoid}, tmsize_t), tif, bp, size)
end

function TIFFSetupStrips(arg1)
    ccall((:TIFFSetupStrips, libtiff), Cint, (Ptr{TIFF},), arg1)
end

function TIFFWriteCheck(arg1, arg2, arg3)
    ccall((:TIFFWriteCheck, libtiff), Cint, (Ptr{TIFF}, Cint, Cstring), arg1, arg2, arg3)
end

function TIFFFreeDirectory(arg1)
    ccall((:TIFFFreeDirectory, libtiff), Cvoid, (Ptr{TIFF},), arg1)
end

function TIFFCreateDirectory(arg1)
    ccall((:TIFFCreateDirectory, libtiff), Cint, (Ptr{TIFF},), arg1)
end

function TIFFCreateCustomDirectory(arg1, arg2)
    ccall((:TIFFCreateCustomDirectory, libtiff), Cint, (Ptr{TIFF}, Ptr{TIFFFieldArray}), arg1, arg2)
end

function TIFFCreateEXIFDirectory(arg1)
    ccall((:TIFFCreateEXIFDirectory, libtiff), Cint, (Ptr{TIFF},), arg1)
end

function TIFFLastDirectory(arg1)
    ccall((:TIFFLastDirectory, libtiff), Cint, (Ptr{TIFF},), arg1)
end

function TIFFSetDirectory(arg1, arg2)
    ccall((:TIFFSetDirectory, libtiff), Cint, (Ptr{TIFF}, uint16), arg1, arg2)
end

function TIFFSetSubDirectory(arg1, arg2)
    ccall((:TIFFSetSubDirectory, libtiff), Cint, (Ptr{TIFF}, uint64), arg1, arg2)
end

function TIFFUnlinkDirectory(arg1, arg2)
    ccall((:TIFFUnlinkDirectory, libtiff), Cint, (Ptr{TIFF}, uint16), arg1, arg2)
end

function TIFFUnsetField(arg1, arg2)
    ccall((:TIFFUnsetField, libtiff), Cint, (Ptr{TIFF}, uint32), arg1, arg2)
end

function TIFFWriteDirectory(arg1)
    ccall((:TIFFWriteDirectory, libtiff), Cint, (Ptr{TIFF},), arg1)
end

function TIFFWriteCustomDirectory(arg1, arg2)
    ccall((:TIFFWriteCustomDirectory, libtiff), Cint, (Ptr{TIFF}, Ptr{uint64}), arg1, arg2)
end

function TIFFCheckpointDirectory(arg1)
    ccall((:TIFFCheckpointDirectory, libtiff), Cint, (Ptr{TIFF},), arg1)
end

function TIFFRewriteDirectory(arg1)
    ccall((:TIFFRewriteDirectory, libtiff), Cint, (Ptr{TIFF},), arg1)
end

function TIFFPrintDirectory(arg1, arg2, arg3)
    ccall((:TIFFPrintDirectory, libtiff), Cvoid, (Ptr{TIFF}, Ptr{FILE}, Clong), arg1, arg2, arg3)
end

function TIFFReadScanline(tif, buf, row, sample)
    ccall((:TIFFReadScanline, libtiff), Cint, (Ptr{TIFF}, Ptr{Cvoid}, uint32, uint16), tif, buf, row, sample)
end

function TIFFWriteScanline(tif, buf, row, sample)
    ccall((:TIFFWriteScanline, libtiff), Cint, (Ptr{TIFF}, Ptr{Cvoid}, uint32, uint16), tif, buf, row, sample)
end

function TIFFReadRGBAImage(arg1, arg2, arg3, arg4, arg5)
    ccall((:TIFFReadRGBAImage, libtiff), Cint, (Ptr{TIFF}, uint32, uint32, Ptr{uint32}, Cint), arg1, arg2, arg3, arg4, arg5)
end

function TIFFReadRGBAImageOriented(arg1, arg2, arg3, arg4, arg5, arg6)
    ccall((:TIFFReadRGBAImageOriented, libtiff), Cint, (Ptr{TIFF}, uint32, uint32, Ptr{uint32}, Cint, Cint), arg1, arg2, arg3, arg4, arg5, arg6)
end

function TIFFReadRGBAStrip(arg1, arg2, arg3)
    ccall((:TIFFReadRGBAStrip, libtiff), Cint, (Ptr{TIFF}, uint32, Ptr{uint32}), arg1, arg2, arg3)
end

function TIFFReadRGBATile(arg1, arg2, arg3, arg4)
    ccall((:TIFFReadRGBATile, libtiff), Cint, (Ptr{TIFF}, uint32, uint32, Ptr{uint32}), arg1, arg2, arg3, arg4)
end

function TIFFReadRGBAStripExt(arg1, arg2, arg3, stop_on_error)
    ccall((:TIFFReadRGBAStripExt, libtiff), Cint, (Ptr{TIFF}, uint32, Ptr{uint32}, Cint), arg1, arg2, arg3, stop_on_error)
end

function TIFFReadRGBATileExt(arg1, arg2, arg3, arg4, stop_on_error)
    ccall((:TIFFReadRGBATileExt, libtiff), Cint, (Ptr{TIFF}, uint32, uint32, Ptr{uint32}, Cint), arg1, arg2, arg3, arg4, stop_on_error)
end

function TIFFRGBAImageOK(arg1, arg2)
    ccall((:TIFFRGBAImageOK, libtiff), Cint, (Ptr{TIFF}, Ptr{UInt8}), arg1, arg2)
end

function TIFFRGBAImageBegin(arg1, arg2, arg3, arg4)
    ccall((:TIFFRGBAImageBegin, libtiff), Cint, (Ptr{TIFFRGBAImage}, Ptr{TIFF}, Cint, Ptr{UInt8}), arg1, arg2, arg3, arg4)
end

function TIFFRGBAImageGet(arg1, arg2, arg3, arg4)
    ccall((:TIFFRGBAImageGet, libtiff), Cint, (Ptr{TIFFRGBAImage}, Ptr{uint32}, uint32, uint32), arg1, arg2, arg3, arg4)
end

function TIFFRGBAImageEnd(arg1)
    ccall((:TIFFRGBAImageEnd, libtiff), Cvoid, (Ptr{TIFFRGBAImage},), arg1)
end

function TIFFOpen(arg1, arg2)
    ccall((:TIFFOpen, libtiff), Ptr{TIFF}, (Cstring, Cstring), arg1, arg2)
end

function TIFFFdOpen(arg1, arg2, arg3)
    ccall((:TIFFFdOpen, libtiff), Ptr{TIFF}, (Cint, Cstring, Cstring), arg1, arg2, arg3)
end

function TIFFClientOpen(arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9, arg10)
    ccall((:TIFFClientOpen, libtiff), Ptr{TIFF}, (Cstring, Cstring, thandle_t, TIFFReadWriteProc, TIFFReadWriteProc, TIFFSeekProc, TIFFCloseProc, TIFFSizeProc, TIFFMapFileProc, TIFFUnmapFileProc), arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9, arg10)
end

function TIFFFileName(arg1)
    ccall((:TIFFFileName, libtiff), Cstring, (Ptr{TIFF},), arg1)
end

function TIFFSetFileName(arg1, arg2)
    ccall((:TIFFSetFileName, libtiff), Cstring, (Ptr{TIFF}, Cstring), arg1, arg2)
end

function TIFFSetErrorHandler(arg1)
    ccall((:TIFFSetErrorHandler, libtiff), TIFFErrorHandler, (TIFFErrorHandler,), arg1)
end

function TIFFSetErrorHandlerExt(arg1)
    ccall((:TIFFSetErrorHandlerExt, libtiff), TIFFErrorHandlerExt, (TIFFErrorHandlerExt,), arg1)
end

function TIFFSetWarningHandler(arg1)
    ccall((:TIFFSetWarningHandler, libtiff), TIFFErrorHandler, (TIFFErrorHandler,), arg1)
end

function TIFFSetWarningHandlerExt(arg1)
    ccall((:TIFFSetWarningHandlerExt, libtiff), TIFFErrorHandlerExt, (TIFFErrorHandlerExt,), arg1)
end

function TIFFSetTagExtender(arg1)
    ccall((:TIFFSetTagExtender, libtiff), TIFFExtendProc, (TIFFExtendProc,), arg1)
end

function TIFFComputeTile(tif, x, y, z, s)
    ccall((:TIFFComputeTile, libtiff), uint32, (Ptr{TIFF}, uint32, uint32, uint32, uint16), tif, x, y, z, s)
end

function TIFFCheckTile(tif, x, y, z, s)
    ccall((:TIFFCheckTile, libtiff), Cint, (Ptr{TIFF}, uint32, uint32, uint32, uint16), tif, x, y, z, s)
end

function TIFFNumberOfTiles(arg1)
    ccall((:TIFFNumberOfTiles, libtiff), uint32, (Ptr{TIFF},), arg1)
end

function TIFFReadTile(tif, buf, x, y, z, s)
    ccall((:TIFFReadTile, libtiff), tmsize_t, (Ptr{TIFF}, Ptr{Cvoid}, uint32, uint32, uint32, uint16), tif, buf, x, y, z, s)
end

function TIFFWriteTile(tif, buf, x, y, z, s)
    ccall((:TIFFWriteTile, libtiff), tmsize_t, (Ptr{TIFF}, Ptr{Cvoid}, uint32, uint32, uint32, uint16), tif, buf, x, y, z, s)
end

function TIFFComputeStrip(arg1, arg2, arg3)
    ccall((:TIFFComputeStrip, libtiff), uint32, (Ptr{TIFF}, uint32, uint16), arg1, arg2, arg3)
end

function TIFFNumberOfStrips(arg1)
    ccall((:TIFFNumberOfStrips, libtiff), uint32, (Ptr{TIFF},), arg1)
end

function TIFFReadEncodedStrip(tif, strip, buf, size)
    ccall((:TIFFReadEncodedStrip, libtiff), tmsize_t, (Ptr{TIFF}, uint32, Ptr{Cvoid}, tmsize_t), tif, strip, buf, size)
end

function TIFFReadRawStrip(tif, strip, buf, size)
    ccall((:TIFFReadRawStrip, libtiff), tmsize_t, (Ptr{TIFF}, uint32, Ptr{Cvoid}, tmsize_t), tif, strip, buf, size)
end

function TIFFReadEncodedTile(tif, tile, buf, size)
    ccall((:TIFFReadEncodedTile, libtiff), tmsize_t, (Ptr{TIFF}, uint32, Ptr{Cvoid}, tmsize_t), tif, tile, buf, size)
end

function TIFFReadRawTile(tif, tile, buf, size)
    ccall((:TIFFReadRawTile, libtiff), tmsize_t, (Ptr{TIFF}, uint32, Ptr{Cvoid}, tmsize_t), tif, tile, buf, size)
end

function TIFFWriteEncodedStrip(tif, strip, data, cc)
    ccall((:TIFFWriteEncodedStrip, libtiff), tmsize_t, (Ptr{TIFF}, uint32, Ptr{Cvoid}, tmsize_t), tif, strip, data, cc)
end

function TIFFWriteRawStrip(tif, strip, data, cc)
    ccall((:TIFFWriteRawStrip, libtiff), tmsize_t, (Ptr{TIFF}, uint32, Ptr{Cvoid}, tmsize_t), tif, strip, data, cc)
end

function TIFFWriteEncodedTile(tif, tile, data, cc)
    ccall((:TIFFWriteEncodedTile, libtiff), tmsize_t, (Ptr{TIFF}, uint32, Ptr{Cvoid}, tmsize_t), tif, tile, data, cc)
end

function TIFFWriteRawTile(tif, tile, data, cc)
    ccall((:TIFFWriteRawTile, libtiff), tmsize_t, (Ptr{TIFF}, uint32, Ptr{Cvoid}, tmsize_t), tif, tile, data, cc)
end

function TIFFDataWidth(arg1)
    ccall((:TIFFDataWidth, libtiff), Cint, (TIFFDataType,), arg1)
end

function TIFFSetWriteOffset(tif, off)
    ccall((:TIFFSetWriteOffset, libtiff), Cvoid, (Ptr{TIFF}, toff_t), tif, off)
end

function TIFFSwabShort(arg1)
    ccall((:TIFFSwabShort, libtiff), Cvoid, (Ptr{uint16},), arg1)
end

function TIFFSwabLong(arg1)
    ccall((:TIFFSwabLong, libtiff), Cvoid, (Ptr{uint32},), arg1)
end

function TIFFSwabLong8(arg1)
    ccall((:TIFFSwabLong8, libtiff), Cvoid, (Ptr{uint64},), arg1)
end

function TIFFSwabFloat(arg1)
    ccall((:TIFFSwabFloat, libtiff), Cvoid, (Ptr{Cfloat},), arg1)
end

function TIFFSwabDouble(arg1)
    ccall((:TIFFSwabDouble, libtiff), Cvoid, (Ptr{Cdouble},), arg1)
end

function TIFFSwabArrayOfShort(wp, n)
    ccall((:TIFFSwabArrayOfShort, libtiff), Cvoid, (Ptr{uint16}, tmsize_t), wp, n)
end

function TIFFSwabArrayOfTriples(tp, n)
    ccall((:TIFFSwabArrayOfTriples, libtiff), Cvoid, (Ptr{uint8}, tmsize_t), tp, n)
end

function TIFFSwabArrayOfLong(lp, n)
    ccall((:TIFFSwabArrayOfLong, libtiff), Cvoid, (Ptr{uint32}, tmsize_t), lp, n)
end

function TIFFSwabArrayOfLong8(lp, n)
    ccall((:TIFFSwabArrayOfLong8, libtiff), Cvoid, (Ptr{uint64}, tmsize_t), lp, n)
end

function TIFFSwabArrayOfFloat(fp, n)
    ccall((:TIFFSwabArrayOfFloat, libtiff), Cvoid, (Ptr{Cfloat}, tmsize_t), fp, n)
end

function TIFFSwabArrayOfDouble(dp, n)
    ccall((:TIFFSwabArrayOfDouble, libtiff), Cvoid, (Ptr{Cdouble}, tmsize_t), dp, n)
end

function TIFFReverseBits(cp, n)
    ccall((:TIFFReverseBits, libtiff), Cvoid, (Ptr{uint8}, tmsize_t), cp, n)
end

function TIFFGetBitRevTable(arg1)
    ccall((:TIFFGetBitRevTable, libtiff), Ptr{Cuchar}, (Cint,), arg1)
end

function LogL16toY(arg1)
    ccall((:LogL16toY, libtiff), Cdouble, (Cint,), arg1)
end

function LogL10toY(arg1)
    ccall((:LogL10toY, libtiff), Cdouble, (Cint,), arg1)
end

function XYZtoRGB24(arg1, arg2)
    ccall((:XYZtoRGB24, libtiff), Cvoid, (Ptr{Cfloat}, Ptr{uint8}), arg1, arg2)
end

function uv_decode(arg1, arg2, arg3)
    ccall((:uv_decode, libtiff), Cint, (Ptr{Cdouble}, Ptr{Cdouble}, Cint), arg1, arg2, arg3)
end

function LogLuv24toXYZ(arg1, arg2)
    ccall((:LogLuv24toXYZ, libtiff), Cvoid, (uint32, Ptr{Cfloat}), arg1, arg2)
end

function LogLuv32toXYZ(arg1, arg2)
    ccall((:LogLuv32toXYZ, libtiff), Cvoid, (uint32, Ptr{Cfloat}), arg1, arg2)
end

function LogL16fromY(arg1, arg2)
    ccall((:LogL16fromY, libtiff), Cint, (Cdouble, Cint), arg1, arg2)
end

function LogL10fromY(arg1, arg2)
    ccall((:LogL10fromY, libtiff), Cint, (Cdouble, Cint), arg1, arg2)
end

function uv_encode(arg1, arg2, arg3)
    ccall((:uv_encode, libtiff), Cint, (Cdouble, Cdouble, Cint), arg1, arg2, arg3)
end

function LogLuv24fromXYZ(arg1, arg2)
    ccall((:LogLuv24fromXYZ, libtiff), uint32, (Ptr{Cfloat}, Cint), arg1, arg2)
end

function LogLuv32fromXYZ(arg1, arg2)
    ccall((:LogLuv32fromXYZ, libtiff), uint32, (Ptr{Cfloat}, Cint), arg1, arg2)
end

function TIFFCIELabToRGBInit(arg1, arg2, arg3)
    ccall((:TIFFCIELabToRGBInit, libtiff), Cint, (Ptr{TIFFCIELabToRGB}, Ptr{TIFFDisplay}, Ptr{Cfloat}), arg1, arg2, arg3)
end

function TIFFCIELabToXYZ(arg1, arg2, arg3, arg4, arg5, arg6, arg7)
    ccall((:TIFFCIELabToXYZ, libtiff), Cvoid, (Ptr{TIFFCIELabToRGB}, uint32, int32, int32, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}), arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function TIFFXYZToRGB(arg1, arg2, arg3, arg4, arg5, arg6, arg7)
    ccall((:TIFFXYZToRGB, libtiff), Cvoid, (Ptr{TIFFCIELabToRGB}, Cfloat, Cfloat, Cfloat, Ptr{uint32}, Ptr{uint32}, Ptr{uint32}), arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function TIFFYCbCrToRGBInit(arg1, arg2, arg3)
    ccall((:TIFFYCbCrToRGBInit, libtiff), Cint, (Ptr{TIFFYCbCrToRGB}, Ptr{Cfloat}, Ptr{Cfloat}), arg1, arg2, arg3)
end

function TIFFYCbCrtoRGB(arg1, arg2, arg3, arg4, arg5, arg6, arg7)
    ccall((:TIFFYCbCrtoRGB, libtiff), Cvoid, (Ptr{TIFFYCbCrToRGB}, uint32, int32, int32, Ptr{uint32}, Ptr{uint32}, Ptr{uint32}), arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function TIFFMergeFieldInfo(arg1, arg2, arg3)
    ccall((:TIFFMergeFieldInfo, libtiff), Cint, (Ptr{TIFF}, Ptr{TIFFFieldInfo}, uint32), arg1, arg2, arg3)
end
# Julia wrapper for header: tiffvers.h
# Automatically generated using Clang.jl


text:
export getConvMatPeriodic, getAverageMatrix, getCoupledConvMat

function getAverageMatrix(TYPE::Type,nImg,nc)

    A1 = spdiagm((fill(convert(TYPE,.5),nImg[1]),fill(convert(TYPE,.5),nImg[1]-1)),(0,1),nImg[1],nImg[1])
    A1 = A1[1:2:end,:]
    A2 = spdiagm((fill(convert(TYPE,.5),nImg[2]),fill(convert(TYPE,.5),nImg[2]-1)),(0,1),nImg[2],nImg[2])
    A2 = A2[1:2:end,:]
    Av    = kron(A2,A1);
    return kron(speye(TYPE,nc),Av);
end
function getCoupledConvMat(TYPE::Type,K::Array,nImg,sK;stride=[1;1])
    # build single 2D convolution operator
    nK  = prod(sK[1:2])
    K   = reshape(K,tuple(sK...)) # basic kernel
    K = flipdim(K,2)
    K = flipdim(K,1)
    # K   = permutedims(K,(2,1,3,4))
    # apply stride
    idx = falses(nImg[1],nImg[2])
    idx[1:stride[1]:end,:] = true
    idx[:,1:stride[2]:end] = true

    A = spzeros(TYPE,0,prod(nImg[1:2])*sK[3])
    for i=1:sK[4]
        Ak = spzeros(TYPE,prod(nImg[1:2]),0)
        for j=1:sK[3]
            Ak = [Ak getConvMatPeriodic(TYPE,K[:,:,j,i],[nImg[1:2]; 1])]
        end
        A = [A; Ak[vec(idx).==true,:]]
    end
    return A
end

function getConvMatPeriodic(TYPE::Type,K::Array,n::Array{Int64})

G  = reshape(collect(1:prod(n)),tuple(n...))

jj = zeros(Int64,0); ii = zeros(Int64,0); vv = zeros(TYPE,0);
I1,I2,I3 = ndgrid(1:n[1],1:n[2],1:n[3]);
for i=1:size(K,1)
    for j=1:size(K,2)
        for k=1:size(K,3)
            ofi = i - div(size(K,1)-1,2) - 1
            ofj = j - div(size(K,2)-1,2) - 1
            ofk = k - div(size(K,3)-1,2) - 1

            t1 = vec(I1.+ofi)
            t1[t1.<=0]   += n[1]
            t1[t1.>n[1]] -= n[1]

            t2 = vec(I2.+ofj)
            t2[t2.<=0]   += n[2]
            t2[t2.>n[2]] -= n[2]

            t3 = vec(I3.+ofk)
			t3[t3.<=0]   += n[3]
            t3[t3.>n[3]] -= n[3]

            jj = [jj; vec(G)];
            ii = [ii; vec(G[sub2ind(tuple(n...),t1, t2,t3)])];
            vv = [vv; ones(prod(n))*K[i,j,k]];
        end
    end
end
A = convert(SparseMatrixCSC{TYPE},sparse(ii, jj ,vv,prod(n),prod(n)));
return A
end


text:
using LevenshteinToolkit
using Test

tests = [ "matrix_strings", "row_strings",
          "matrix_arrays",  "row_arrays",
          "matrix_costs",   "row_costs",
          "matrix_stress",  "row_stress" ]

for test ∈ tests
    include("$test.jl")
end

max_stars_repo_name: jlegare/LevenshteinToolkit

max_issues_repo_name: jlegare/LevenshteinToolkit

max_forks_repo_name: jlegare/LevenshteinToolkit


text:
module dapr
  const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
  module proto
    const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
    module operator
      const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
      module v1
        const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
        include("operator_pb.jl")
      end
    end
  end
end


text:
"""
    FlatLens(N₁, N₂, ..., Nₙ)

# Examples
```jldoctest
julia> using Setfield, Kaleido

julia> l = MultiLens((
           (@lens _.x) ∘ IndexBatchLens(:a, :b, :c),
           (@lens _.y) ∘ IndexBatchLens(:d, :e),
       )) ∘ FlatLens(3, 2);

julia> get((x=(a=1, b=2, c=3), y=(d=4, e=5)), l)
(1, 2, 3, 4, 5)

julia> set((x=(a=1, b=2, c=3), y=(d=4, e=5)), l, (10, 20, 30, 40, 50))
(x = (a = 10, b = 20, c = 30), y = (d = 40, e = 50))
```
"""
FlatLens

struct FlatLens{lengths} <: Lens end

Setfield.get(obj, ::FlatLens) = _cat(obj...)

Setfield.set(::Any, ::FlatLens{lengths}, val) where lengths =
    let val = Tuple(val)
        _foldl(lengths, ((), 1)) do (ys, start), n
            (ys..., val[start:start + n - 1]), start + n
        end[1]
    end

FlatLens(lengths::Vararg{Integer}) = FlatLens{UInt.(lengths)}()

Base.show(io::IO, ::FlatLens{lengths}) where lengths =
    print_apply(io, FlatLens, Int.(lengths))

"""
    SingletonLens()

Inverse of `FlatLens(1)`.
"""
struct SingletonLens <: Lens end

Setfield.get(obj, ::SingletonLens) = (obj,)
Setfield.set(::Any, ::SingletonLens, obj) = obj[1]


text:
include("sec_6_2.jl")

using PyPlot

r = 3
M = 2r
ppI = 40
store = Store(α, r+1, M, ppI)

N = 5
t, U = FODEdG!(λ, tmax, f, u0, N, r, M, store)
τ = FractionalTimeDG.reconstruction_pts(T, r)
pcwise_t, pcwise_U = evaluate_pcwise_poly!(U, t, τ[0:r], store)
pcwise_u = [ u(t, λ, u0, f, 10, store) for t in pcwise_t ]
pcwise_err = pcwise_U - pcwise_u
Uhat = reconstruction(U, u0, store)
pcwise_t_fine, pcwise_Uhat = evaluate_pcwise_poly!(Uhat, t, ppI, store)
pcwise_u_fine = [ u(t, λ, u0, f, 10, store) for t in pcwise_t_fine ]
pcwise_rec_err = pcwise_Uhat - pcwise_u_fine
pcwise_rec_err[1,1] = NaN
a = FractionalTimeDG.Fourier_Legendre_coefs(r, t, store) do s
    return u(s, λ, u0, f, 10, store)
end

figure(4)
trans_pcwise_t = copy(pcwise_t')
trans_pcwise_err = copy(pcwise_err')
line0 = semilogy(pcwise_t_fine, abs.(pcwise_rec_err), "k:")
line1 = semilogy(trans_pcwise_t[:,1], abs.(trans_pcwise_err[:,1]), "o",
                 markersize=3)
line2 = semilogy(trans_pcwise_t[:,2], abs.(trans_pcwise_err[:,2]), "v",
                 markersize=3)
line3 = semilogy(trans_pcwise_t[:,3], abs.(trans_pcwise_err[:,3]), "+",
                 markersize=4)
line4 = semilogy(trans_pcwise_t[:,4], abs.(trans_pcwise_err[:,4]), "x",
                 markersize=4)
legend((line0[1], line1[1], line2[1], line3[1], line4[1]),
       (L"$\widehat E(t)$", L"$j=0$", L"j=1", L"j=2", L"j=3"))
xlabel(L"$t$", size=12)
ylabel("Error")
grid(true)
savefig("fig4.eps")


text:
using GuessworkQuantumSideInfo
using SDPAFamily
using Test
using LinearAlgebra, GenericLinearAlgebra

# Fix a compatibility problem
LinearAlgebra.eigmin(A::Hermitian{Complex{BigFloat},Array{Complex{BigFloat},2}}) =
    minimum(real.(eigvals(A)))

function testPOVM(Es; tol = 1e-25)
    dB = size(first(Es), 1)
    @test sum(Es) ≈ complex(1.0) * I(dB) atol = tol

    for E in Es
        @test E ≈ E' atol = tol
    end
    @test all(isposdef, Hermitian(E) + tol * I(dB) for E in Es)
end

T = BigFloat
default_sdp_solver(T) = SDPAFamily.Optimizer{T}(presolve = true)

@testset "BB84" begin
    ρBs = GuessworkQuantumSideInfo.BB84_states(T)

    p = ones(T, 4) / 4
    output = guesswork(p, ρBs; solver = default_sdp_solver(T))
    testPOVM(output.Es)
    relaxed_output =
        guesswork_upper_bound(p, ρBs; num_constraints = 4, make_solver = () -> default_sdp_solver(T))

    true_val = (big(1) / big(4)) * (10 - sqrt(big(10)))
    @test output.optval ≈ true_val rtol = 1e-25

    @test output.optval ≈ relaxed_output.optval rtol = 1e-4

    lb = guesswork_lower_bound(p, ρBs; solver = default_sdp_solver(T)).optval
    @test lb <= output.optval + 1e-4
end


text:
using ImplicitFields
using Base.Test

type Foo
    x::Int
    u::Symbol
end

type Bar
    x::Float32
end

f = Foo(2, :result)
b = Bar(1.6f0)

@test resolve(
    :(u += x * x_post + a * dt),
    f, post=b, dt=1e-3) == :(result += 2 * 1.6f0 + a * 0.001)


text:
using ITensors
using ITensorGPU

# Set to identity to run on CPU
gpu = cu

N = 50
sites = siteinds("S=1", N)

ampo = AutoMPO()
for j in 1:(N - 1)
  ampo .+= 0.5, "S+", j, "S-", j + 1
  ampo .+= 0.5, "S-", j, "S+", j + 1
  ampo .+= "Sz", j, "Sz", j + 1
end
H = gpu(MPO(ampo, sites))

ψ₀ = gpu(randomMPS(sites))

sweeps = Sweeps(6)
maxdim!(sweeps, 10, 20, 40, 100)
mindim!(sweeps, 1, 10)
cutoff!(sweeps, 1e-11)
noise!(sweeps, 1e-10)
energy, ψ = @time dmrg(H, ψ₀, sweeps)
@show energy


text:
using LinearAlgebra
using RoME
# , IncrementalInference, TransformUtils, Distributions
using Test
import  IncrementalInference: getSample

##

mutable struct RotationTest <: IncrementalInference.AbstractRelativeRoots
  z::MvNormal
end

getSample(cfo::CalcFactor{<:RotationTest}, N::Int=1) = (reshape(rand(cfo.factor.z,N),3,:),)

# 3 dimensional line, z = [a b][x y]' + c
function (cfo::CalcFactor{<:RotationTest})( meas,
                                            var1,
                                            var2)
  #
  #FIXME JT - I'm createing new res for simplicity, it may not hold up well though
  res = Vector{eltype(var1)}(undef, 3)

  z = meas
  dq = convert(Quaternion, Euler(z...))
  @show var1
  @show s1 = TU.so3(var1)
  s2 = TU.so3(var2)
  q1 = convert(Quaternion, s1)
  q2 = convert(Quaternion, s2)
  q12 = q1*q_conj(q2)
  qq = dq*q_conj(q12)
  @show res
  vee!(res, convert(TU.so3, qq))
  return res
end

rr = RotationTest(MvNormal(zeros(3), 0.001*diagm(ones(3))))

##

@testset "Increased dimension root finding test" begin

## known rotations

eul = zeros(3)

R1 = zeros(3)
R2 = zeros(3)

res = randn(3)

# should be Sphere3
res = testFactorResidualBinary( rr,
                                ContinuousEuclid{3},
                                ContinuousEuclid{3},
                                R1,
                                R2,
                                (eul,))
# rr(res, nothing, 1, (zeros(3,1),), R1, R2)

@test norm(res) < 1e-10

##

# random rotations
eul = 0.25*randn(3)

R1 = rand(3)
R2 = rand(3)

# res = zeros(3)

res = testFactorResidualBinary( rr,
                                ContinuousEuclid{3},
                                ContinuousEuclid{3},
                                R1,
                                R2,
                                (eul,))
# rr(res, nothing, 1, (zeros(3),), R1, R2)

@test norm(res) > 1e-3

##

end

#


text:
" solution status (from ipopt) "
ApplicationReturnStatus = Dict(
0=>:Solve_Succeeded,
1=>:Solved_To_Acceptable_Level,
2=>:Infeasible_Problem_Detected,
3=>:Search_Direction_Becomes_Too_Small,
4=>:Diverging_Iterates,
5=>:User_Requested_Stop,
6=>:Feasible_Point_Found,
-1=>:Maximum_Iterations_Exceeded,
-2=>:Restoration_Failed,
-3=>:Error_In_Step_Computation,
-4=>:Maximum_CpuTime_Exceeded,
-5=>:Optimize_not_called,
-6=>:Method_not_defined,
-10=>:Not_Enough_Degrees_Of_Freedom,
-11=>:Invalid_Problem_Definition,
-12=>:Invalid_Option,
-13=>:Invalid_Number_Detected,
-100=>:Unrecoverable_Exception,
-102=>:Insufficient_Memory,
-199=>:Internal_Error)

max_issues_repo_name: exanauts/ActiveSetMethods