text:
#  main learn parameters
#  runs simulation to learn parameters and equilibrium levels of capital

ENV["GKS_ENCODING"] = "utf-8"

# Julia libraries used
using LinearAlgebra, Statistics,  Distributions, Random,  Plots, JLD2
# project functions used
include("setup_inital.jl")
include("sim_functions.jl")
include("show_functions.jl")

println()
println()
println("********** Setup and Initial Values *************************")
Random.seed!(123)  # set random seed to

# define parameters for simulation --- change parameters in that function (file setup_functions.jl)
sim_set, decide_set, econ_set = q_set_para_initial()

# defines functions used in simulation -- e.g. production function, functional form of decsion rule, ...
include("setup_functions.jl")
econ_set["price_para"] = price_para    # adjust price para to fit market type

# set up accounting arrays that keeps track of variables over time
account_econ, account_para = q_set_accounting(sim_set)

# create initial values -- state of the world and parameters for decsions function
initial_set = q_set_initial_values(sim_set["type"], sim_set["n_agents"])
state = initial_set["state"]
s_para = initial_set["decide_para"]

println("   Setup complete  " )

println()
println()
println("********** Run Learn_shock simulation *************************")

# simulated the economy for specified number of gnerations
account_econ, account_para, s_para, state = q_one_run(state, s_para, econ_set, decide_set, sim_set, account_econ, account_para)

println(" run learn shock complete ")

# show results
q_show_results(state, s_para, econ_set, decide_set, sim_set, account_econ, account_para)


text:
function DefaultTestSolver(;nl_solver=optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0, "sb" =>"yes"), solver_args...)
    solver_args_result = Vector{Pair{String, Any}}()
    push!(solver_args_result, "log_levels" => Symbol[])
    push!(solver_args_result, "nl_solver" => nl_solver)
    for v in solver_args
        push!(solver_args_result, string(v[1]) => v[2])
    end
    return solver_args_result
end

ext: jl


text:
using Romeo, GLFW, GLAbstraction, Reactive, ModernGL, GLWindow
obj = visualize(Float32[rand(Float32)  for i=0:10, j=0:10], color = rgba(1,0,0,1))
println(obj)
push!(Romeo.ROOT_SCREEN.renderlist, obj)

while Romeo.ROOT_SCREEN.inputs[:open].value
    Romeo.renderloop(Romeo.ROOT_SCREEN)
    sleep(0.0001)
end
GLFW.Terminate()


text:
## Interface to the StenoGraphs package

############################################################################################
### Define Modifiers
############################################################################################

# fixed parameter values
struct Fixed{N} <: EdgeModifier
    value::N
end
fixed(args...) = Fixed(args)
Fixed(value::Int) = Fixed(Float64(value))

# start values
struct Start{N} <: EdgeModifier
    value::N
end
start(args...) = Start(args)
Start(value::Int) = Start(Float64(value))

# labels for equality constraints
struct Label{N} <: EdgeModifier
    value::N
end
label(args...) = Label(args)

############################################################################################
### constructor for parameter table from graph
############################################################################################

function ParameterTable(;graph, observed_vars, latent_vars, g = 1, parname = :θ)
    graph = unique(graph)
    n = length(graph)
    from = Vector{Symbol}(undef, n)
    parameter_type = Vector{Symbol}(undef, n)
    to = Vector{Symbol}(undef, n)
    free = ones(Bool, n)
    value_fixed = zeros(n)
    start = zeros(n)
    estimate = zeros(n)
    identifier = Vector{Symbol}(undef, n); identifier .= Symbol("")
    # group = Vector{Symbol}(undef, n)
    # start_partable = zeros(Bool, n)

    sorted_vars = Vector{Symbol}()

    for (i, element) in enumerate(graph)
        if element isa DirectedEdge
            from[i] =  element.src.node
            to[i] =  element.dst.node
            parameter_type[i] = :→
        elseif element isa UndirectedEdge
            from[i] =  element.src.node
            to[i] =  element.dst.node
            parameter_type[i] = :↔
        elseif element isa ModifiedEdge
            if element.edge isa DirectedEdge
                from[i] =  element.edge.src.node
                to[i] =  element.edge.dst.node
                parameter_type[i] = :→
            elseif element.edge isa UndirectedEdge
                from[i] =  element.edge.src.node
                to[i] =  element.edge.dst.node
                parameter_type[i] = :↔
            end
            for modifier in values(element.modifiers)
                if modifier isa Fixed
                    if modifier.value[g] == :NaN
                        free[i] = true
                        value_fixed[i] = 0.0
                    else
                        free[i] = false
                        value_fixed[i] = modifier.value[g]
                    end
                elseif modifier isa Start
                    start_partable[i] = modifier.value[g] == :NaN
                    start[i] = modifier.value[g]
                elseif modifier isa Label
                    if modifier.value[g] == :NaN
                        throw(DomainError(NaN, "NaN is not allowed as a parameter label."))
                    end
                    identifier[i] = modifier.value[g]
                end
            end
        end
    end

    # make identifiers for parameters that are not labeled
    current_id = 1
    for i in 1:length(identifier)
        if (identifier[i] == Symbol("")) & free[i]
            identifier[i] = Symbol(parname, :_, current_id)
            current_id += 1
        elseif (identifier[i] == Symbol("")) & !free[i]
            identifier[i] = :const
        elseif (identifier[i] != Symbol("")) & !free[i]
            @warn "You labeled a constant. Please check if the labels of your graph are correct."
        end
    end

    return StructuralEquationModels.ParameterTable(
        Dict(
            :from => from,
            :parameter_type => parameter_type,
            :to => to,
            :free => free,
            :value_fixed => value_fixed,
            :start => start,
            :estimate => estimate,
            :identifier => identifier),
        Dict(
            :latent_vars => latent_vars,
            :observed_vars => observed_vars,
            :sorted_vars => sorted_vars)
    )
end

############################################################################################
### constructor for EnsembleParameterTable from graph
############################################################################################

function EnsembleParameterTable(;graph, observed_vars, latent_vars, groups)

    graph = unique(graph)

    partable = EnsembleParameterTable(nothing)

    for (i, group) in enumerate(groups)
        push!(
            partable.tables,
            Symbol(group) =>
                ParameterTable(;
                graph = graph,
                observed_vars = observed_vars,
                latent_vars = latent_vars,
                g = i,
                parname = Symbol(:g, i)))
    end

        return partable
end


text:
# Julia wrapper for header: /usr/include/scip/cons_varbound.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeConshdlrVarbound(scip)
    ccall((:SCIPincludeConshdlrVarbound, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPcreateConsVarbound(scip, cons, name, var, vbdvar, vbdcoef, lhs, rhs, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateConsVarbound, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{SCIP_VAR}, Ptr{SCIP_VAR}, Cdouble, Cdouble, Cdouble, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, var, vbdvar, vbdcoef, lhs, rhs, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
end

function SCIPcreateConsBasicVarbound(scip, cons, name, var, vbdvar, vbdcoef, lhs, rhs)
    ccall((:SCIPcreateConsBasicVarbound, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{SCIP_VAR}, Ptr{SCIP_VAR}, Cdouble, Cdouble, Cdouble), scip, cons, name, var, vbdvar, vbdcoef, lhs, rhs)
end

function SCIPgetLhsVarbound(scip, cons)
    ccall((:SCIPgetLhsVarbound, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetRhsVarbound(scip, cons)
    ccall((:SCIPgetRhsVarbound, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetVarVarbound(scip, cons)
    ccall((:SCIPgetVarVarbound, libscip), Ptr{SCIP_VAR}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetVbdvarVarbound(scip, cons)
    ccall((:SCIPgetVbdvarVarbound, libscip), Ptr{SCIP_VAR}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetVbdcoefVarbound(scip, cons)
    ccall((:SCIPgetVbdcoefVarbound, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetDualsolVarbound(scip, cons)
    ccall((:SCIPgetDualsolVarbound, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetDualfarkasVarbound(scip, cons)
    ccall((:SCIPgetDualfarkasVarbound, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetRowVarbound(scip, cons)
    ccall((:SCIPgetRowVarbound, libscip), Ptr{SCIP_ROW}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end


text:
using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = test_page_model()
    include("../src/components/N2Ocycle.jl")

    add_comp!(m, n2ocycle)

    set_param!(m, :n2ocycle, :e_globalN2Oemissions, readpagedata(m, "test/validationdata/$valdir/e_globalN2Oemissions.csv"))
    set_param!(m, :n2ocycle, :y_year, [2020.,2030.,2040.,2050.,2075.,2100.,2150.,2200.,2250.,2300.]) # real value
    set_param!(m, :n2ocycle, :y_year_0, 2015.) # real value
    set_param!(m, :n2ocycle, :rtl_g_landtemperature, readpagedata(m, "test/validationdata/$valdir/rtl_g_landtemperature.csv"))

    ##running Model
    run(m)

    conc = m[:n2ocycle,  :c_N2Oconcentration]
    conc_compare = readpagedata(m, "test/validationdata/$valdir/c_n2oconcentration.csv")

    @test conc ≈ conc_compare rtol = 1e-4
end


text:
using Printf
using FiniteDiff
using ForwardDiff
using KernelAbstractions
using LinearAlgebra
using Random
using SparseArrays
using Test
using TimerOutputs
using ExaPF
import ExaPF: PowerSystem, AutoDiff

const PS = PowerSystem

@testset "Compute reduced gradient on CPU" begin
    @testset "Case $case" for case in ["case9.m", "case30.m"]
        datafile = joinpath(dirname(@__FILE__), "..", "data", case)
        tolerance = 1e-8
        pf = PS.PowerNetwork(datafile, 1)

        polar = PolarForm(pf, CPU())
        cache = ExaPF.get(polar, ExaPF.PhysicalState())

        xk = ExaPF.initial(polar, State())
        u = ExaPF.initial(polar, Control())
        p = ExaPF.initial(polar, Parameters())

        jx, ju, ∂obj = ExaPF.init_autodiff_factory(polar, cache)

        # solve power flow
        conv = powerflow(polar, jx, cache, tol=1e-12)
        ExaPF.get!(polar, State(), xk, cache)
        # No need to recompute ∇gₓ
        ∇gₓ = jx.J
        ∇gᵤ = ExaPF.jacobian(polar, ju, cache)
        # test jacobian wrt x
        ∇gᵥ = ExaPF.jacobian(polar, jx, cache)
        @test isapprox(∇gₓ, ∇gᵥ)

        # Test with Matpower's Jacobian
        V = cache.vmag .* exp.(im * cache.vang)
        Ybus = pf.Ybus
        J = ExaPF.residual_jacobian(V, Ybus, pf.pv, pf.pq)
        @test isapprox(∇gₓ, J)

        # Test gradients
        @testset "Reduced gradient" begin
            # Refresh cache with new values of vmag and vang
            ExaPF.refresh!(polar, PS.Generator(), PS.ActivePower(), cache)
            # We need uk here for the closure
            uk = copy(u)
            ExaPF.∂cost(polar, ∂obj, cache)
            ∇fₓ = ∂obj.∇fₓ
            ∇fᵤ = ∂obj.∇fᵤ

            ## ADJOINT
            # lamba calculation
            λk  = -(∇gₓ') \ ∇fₓ
            grad_adjoint = ∇fᵤ + ∇gᵤ' * λk
            # ## DIRECT
            S = - inv(Array(∇gₓ)) * ∇gᵤ
            grad_direct = ∇fᵤ + S' * ∇fₓ
            @test isapprox(grad_adjoint, grad_direct)

            # Compare with finite difference
            function reduced_cost(u_)
                # Ensure we remain in the manifold
                ExaPF.transfer!(polar, cache, xk, u_, p)
                convergence = powerflow(polar, jx, cache, tol=1e-14)
                ExaPF.refresh!(polar, PS.Generator(), PS.ActivePower(), cache)
                return ExaPF.cost_production(polar, cache.pg)
            end

            grad_fd = FiniteDiff.finite_difference_gradient(reduced_cost, uk)
            @test isapprox(grad_fd, grad_adjoint, rtol=1e-4)
        end
        @testset "Reduced Jacobian" begin
            for cons in [ExaPF.state_constraint, ExaPF.power_constraints]
                m = ExaPF.size_constraint(polar, cons)
                for icons in 1:m
                    ExaPF.jacobian(polar, cons, icons, ∂obj, cache)
                end
            end
        end
    end
end


text:
function test_random(AT)
    @testset "Random" begin
        @testset "rand" begin  # uniform
            for T in (Float32, Float64, Int64, Int32)
                @test length(rand(AT{T,1}, (4,))) == 4
                @test length(rand(AT{T}, (4,))) == 4
                @test length(rand(AT{T}, 4)) == 4
                @test eltype(rand(AT, 4)) == Float32
                @test length(rand(AT, T, 4)) == 4
                @test length(rand(AT{T,2}, (4,5))) == 20
                @test length(rand(AT, T, 4, 5)) == 20
                A = rand(AT{T,2}, (2,2))
                B = copy(A)
                @test all(A .== B)
                rand!(B)
                @test !any(A .== B)
            end
        end
    end
end


text:
using ICI_simulateur
using Test

Principal() ;


text:
using Documenter, VCFTools

ENV["DOCUMENTER_DEBUG"] = "true"

makedocs(
    format = Documenter.HTML(),
    sitename = "VCFTools",
    modules = [VCFTools],
    clean = true,
    pages = [
        "Home" => "index.md",
        "VCF Summary" => "man/summaryinfo.md",
        "Filter" => "man/filter.md",
        "Convert" => "man/convert.md",
        "Match markers in two VCF" => "man/conformgt.md",
        "API"  => "man/api.md"
    ]
)

deploydocs(
    repo   = "github.com/biona001/VCFTools.jl.git",
    target = "build"
)


text:
"""
Betina Robust Sampled based Formulation

Further information:
  - Fernandes, B., Street, A., ValladA˜ £o, D., e Fernandes, C. (2016). An adaptive robust portfolio optimization model with loss constraints based on data-driven polyhedral uncertainty sets. European Journal of Operational Research, 255(3):961 – 970. ISSN 0377-2217. URL.

"""
struct RobustBetina <: AbstractSampleBased
    sampled_returns::Array{Float64,2}
    number_of_assets::Int
    number_of_samples::Int
end

function RobustBetina(;
    sampled_returns::Array{Float64,2}
)
    number_of_samples, number_of_assets = size(sampled_returns)

    return RobustBetina(
        sampled_returns, number_of_assets, number_of_samples
    )
end

function _portfolio_return_latex_RobustBetina_dual()
    return """
        ```math
        \\max_{\\theta} \\quad  \\theta \\\\
        s.t.  \\quad \\theta \\leq r_s ' w \\quad \\forall s = 1:\\mathcal{S} \\\\
        ```
        """
end

"""
    portfolio_return!(model::JuMP.Model, w, formulation::RobustBetina)

Returns worst case return in Betina's uncertainty set, defined by the following dual problem:

$(_portfolio_return_latex_RobustBetina_dual())
"""
function portfolio_return!(model::JuMP.Model, w, formulation::RobustBetina)
    # auxilary variables
    θ = @variable(model, θ)
    # convex hull
    @constraint(model, sum(formulation.sampled_returns * w , dims=2) .>= θ)

    return θ
end


text:
struct EvalFunc{F} <: Function
  f::F
end
(f::EvalFunc)(args...) = f.f(args...)

NO_TSPAN_PROBS = Union{AbstractLinearProblem,AbstractNonlinearProblem,
  AbstractQuadratureProblem,
  AbstractSteadyStateProblem,AbstractJumpProblem}

has_kwargs(_prob::DEProblem) = has_kwargs(typeof(_prob))
Base.@pure has_kwargs(::Type{T}) where {T} = :kwargs ∈ fieldnames(T)

const allowedkeywords = (
  :dense,
  :saveat,
  :save_idxs,
  :tstops,
  :d_discontinuities,
  :save_everystep,
  :save_on,
  :save_start,
  :save_end,
  :initialize_save,
  :adaptive,
  :abstol,
  :reltol,
  :dt,
  :dtmax,
  :dtmin,
  :force_dtmin,
  :internalnorm,
  :controller,
  :gamma,
  :beta1,
  :beta2,
  :qmax,
  :qmin,
  :qsteady_min,
  :qsteady_max,
  :qoldinit,
  :failfactor,
  :calck,
  :alias_u0,
  :maxiters,
  :callback,
  :isoutofdomain,
  :unstable_check,
  :verbose,
  :merge_callbacks,
  :progress,
  :progress_steps,
  :progress_name,
  :progress_message,
  :timeseries_errors,
  :dense_errors,
  :calculate_errors,
  :initializealg,
  :alg,
  :save_noise,
  :delta,
  :seed,
  :alg_hints,
  :kwargshandle,
  :trajectories,
  :batch_size,
  :sensealg,
  :advance_to_tstop,
  :stop_at_next_tstop,
  # These two are from the default algorithm handling
  :default_set,
  :second_time
)

const KWARGWARN_MESSAGE =
"""
Unrecognized keyword arguments found. Future versions will error.
The only allowed keyword arguments to `solve` are:
$allowedkeywords
See https://diffeq.sciml.ai/stable/basics/common_solver_opts/ for more details.

Set kwargshandle=KeywordArgError for an error message and more details.
Set kwargshandle=KeywordArgSilent to ignore this message.
"""

const KWARGERROR_MESSAGE =
"""
Unrecognized keyword arguments found.
The only allowed keyword arguments to `solve` are:
$allowedkeywords
See https://diffeq.sciml.ai/stable/basics/common_solver_opts/ for more details.
"""

struct CommonKwargError <: Exception
  kwargs
end

function Base.showerror(io::IO, e::CommonKwargError)
  println(io, KWARGERROR_MESSAGE)
  notin = collect(map(x -> x ∉ allowedkeywords, keys(e.kwargs)))
  unrecognized = collect(keys(e.kwargs))[notin]
  print(io, "Unrecognized keyword arguments: $unrecognized")
end

@enum KeywordArgError KeywordArgWarn KeywordArgSilent

const INCOMPATIBLE_U0_MESSAGE =
"""
Initial condition incompatible with functional form.
Detected an in-place function with an initial condition of type Number or SArray.
This is incompatible because Numbers cannot be mutated, i.e.
`x = 2.0; y = 2.0; x .= y` will error.

If using a immutable initial condition type, please use the out-of-place form.
I.e. define the function `du=f(u,p,t)` instead of attempting to "mutate" the immutable `du`.

If your differential equation function was defined with multiple dispatches and one is
in-place, then the automatic detection will choose in-place. In this case, override the
choice in the problem constructor, i.e. `ODEProblem{false}(f,u0,tspan,p,kwargs...)`.

For a longer discussion on mutability vs immutability and in-place vs out-of-place, see:
https://diffeq.sciml.ai/stable/tutorials/faster_ode_example/#Example-Accelerating-a-Non-Stiff-Equation:-The-Lorenz-Equation
"""

struct IncompatibleInitialConditionError <: Exception end

function Base.showerror(io::IO, e::IncompatibleInitialConditionError)
  print(io, INCOMPATIBLE_U0_MESSAGE)
end

const NO_DEFAULT_ALGORITHM_MESSAGE =
"""
Default algorithm choices require DifferentialEquations.jl.
Please specify an algorithm (e.g., `solve(prob, Tsit5())` for an ODE)
or import DifferentialEquations directly.

You can find the list of available solvers at https://diffeq.sciml.ai/stable/solvers/ode_solve/
and its associated pages.
"""

struct NoDefaultAlgorithmError <: Exception end

function Base.showerror(io::IO, e::NoDefaultAlgorithmError)
  print(io, NO_DEFAULT_ALGORITHM_MESSAGE)
end

const NO_TSPAN_MESSAGE =
"""
No tspan is set in the problem or chosen in the init/solve call
"""

struct NoTspanError <: Exception end

function Base.showerror(io::IO, e::NoTspanError)
  print(io, NO_TSPAN_MESSAGE)
end

const NON_SOLVER_MESSAGE =
"""
The arguments to solve are incorrect.
The second argument must be a solver choice, `solve(prob,alg)`
where `alg` is a `<: DEAlgorithm`, i.e. `Tsti5()`.

Please double check the arguments being sent to the solver.
"""

struct NonSolverError <: Exception end

function Base.showerror(io::IO, e::NonSolverError)
  print(io, NON_SOLVER_MESSAGE)
end

const PROBSOLVER_PAIRING_MESSAGE =
"""
Incompatible problem+solver pairing.
For example, this can occur if an ODE solver is passed with an SDEProblem.
Solvers are only capable of handling specific problem types. Please double
check that the chosen pairing is capable for handling the given problems.
"""

struct ProblemSolverPairingError <: Exception
  prob
  alg
end

function Base.showerror(io::IO, e::ProblemSolverPairingError)
  println(io, PROBSOLVER_PAIRING_MESSAGE)
  println(io, "Problem type: $(SciMLBase.__parameterless_type(typeof(e.prob)))")
  println(io, "Solver type: $(SciMLBase.__parameterless_type(typeof(e.alg)))")
  println(io, "Problem types compatible with the chosen solver: $(compatible_problem_types(e.prob,e.alg))")
end

function compatible_problem_types(prob,alg)
  if alg isa AbstractODEAlgorithm
    ODEProblem
  elseif alg isa AbstractSDEAlgorithm
    (SDEProblem, SDDEProblem)
  elseif alg isa AbstractDDEAlgorithm # StochasticDelayDiffEq.jl just uses the SDE alg
    DDEProblem
  elseif alg isa AbstractDAEAlgorithm
    DAEProblem
  elseif alg isa AbstractSteadyStateAlgorithm
    SteadyStateProblem
  end
end

const DIRECT_AUTODIFF_INCOMPATABILITY_MESSAGE =
"""
Incompatible solver + automatic differentiation pairing.
The chosen automatic differentiation algorithm requires the ability
for compiler transforms on the code which is only possible on pure-Julia
solvers such as those from OrdinaryDiffEq.jl. Direct differentiation methods
which require this ability include:

- Direct use of ForwardDiff.jl on the solver
- `ForwardDiffSensitivity`, `ReverseDiffAdjoint`, `TrackerAdjoint`, and `ZygoteAdjoint`
  sensealg choices for adjoint differentiation.

Either switch the choice of solver to a pure Julia method, or change the automatic
differentiation method to one that does not require such transformations.

For more details on automatic differentiation, adjoint, and sensitivity analysis
of differential equations, see the documentation page:

https://diffeq.sciml.ai/stable/analysis/sensitivity/
"""

struct DirectAutodiffError <: Exception end

function Base.showerror(io::IO, e::DirectAutodiffError)
  println(io, DIRECT_AUTODIFF_INCOMPATABILITY_MESSAGE)
end

function init_call(_prob, args...; merge_callbacks=true, kwargshandle=KeywordArgWarn, kwargs...)

  if has_kwargs(_prob)
    if merge_callbacks && haskey(_prob.kwargs, :callback) && haskey(kwargs, :callback)
      kwargs_temp = NamedTuple{Base.diff_names(Base._nt_names(
          values(kwargs)), (:callback,))}(values(kwargs))
      callbacks = NamedTuple{(:callback,)}((DiffEqBase.CallbackSet(_prob.kwargs[:callback], values(kwargs).callback),))
      kwargs = merge(kwargs_temp, callbacks)
    end
    kwargs = isempty(_prob.kwargs) ? kwargs : merge(values(_prob.kwargs), kwargs)
  end

  checkkwargs(kwargshandle; kwargs...)

  if hasfield(typeof(_prob), :f) && hasfield(typeof(_prob.f), :f) && typeof(_prob.f.f) <: EvalFunc
    Base.invokelatest(__init, _prob, args...; kwargs...)#::T
  else
    __init(_prob, args...; kwargs...)#::T
  end
end

function init(prob::DEProblem, args...; kwargs...)
  if haskey(kwargs, :alg) && (isempty(args) || args[1] === nothing)
    alg = kwargs[:alg]
    _prob = get_concrete_problem(prob, isadaptive(alg); kwargs...)
    check_prob_alg_pairing(_prob, alg)
    init_call(_prob, alg, args...; kwargs...)
  elseif !isempty(args) && typeof(args[1]) <: DEAlgorithm
    alg = args[1]
    _prob = get_concrete_problem(prob, isadaptive(alg); kwargs...)
    check_prob_alg_pairing(_prob, alg)
    init_call(_prob, args...; kwargs...)
  else
    _prob = get_concrete_problem(prob, !(typeof(prob) <: DiscreteProblem); kwargs...)
    init_call(_prob, args...; kwargs...)
  end
end

function solve_call(_prob, args...; merge_callbacks=true, kwargshandle=KeywordArgWarn, kwargs...)
  if has_kwargs(_prob)
    if merge_callbacks && haskey(_prob.kwargs, :callback) && haskey(kwargs, :callback)
      kwargs_temp = NamedTuple{Base.diff_names(Base._nt_names(
          values(kwargs)), (:callback,))}(values(kwargs))
      callbacks = NamedTuple{(:callback,)}((DiffEqBase.CallbackSet(_prob.kwargs[:callback], values(kwargs).callback),))
      kwargs = merge(kwargs_temp, callbacks)
    end
    kwargs = isempty(_prob.kwargs) ? kwargs : merge(values(_prob.kwargs), kwargs)
  end

  checkkwargs(kwargshandle; kwargs...)

  if hasfield(typeof(_prob), :f) && hasfield(typeof(_prob.f), :f) && typeof(_prob.f.f) <: EvalFunc
    Base.invokelatest(__solve, _prob, args...; kwargs...)#::T
  else
    __solve(_prob, args...; kwargs...)#::T
  end
end

# save_idxs and saveat are here due to https://github.com/FluxML/Zygote.jl/issues/664
function solve(prob::DEProblem, args...; sensealg=nothing,
  u0=nothing, p=nothing, kwargs...)
  u0 = u0 !== nothing ? u0 : prob.u0
  p = p !== nothing ? p : prob.p
  if sensealg === nothing && haskey(prob.kwargs, :sensealg)
    sensealg = prob.kwargs[:sensealg]
  end
  solve_up(prob, sensealg, u0, p, args...; kwargs...)
end

function solve_up(prob::DEProblem, sensealg, u0, p, args...; kwargs...)

  if haskey(kwargs, :alg) && (isempty(args) || args[1] === nothing)
    alg = kwargs[:alg]
    _alg = prepare_alg(alg, u0, p, prob)
    _prob = get_concrete_problem(prob, isadaptive(_alg); u0=u0, p=p, kwargs...)
    check_prob_alg_pairing(_prob, alg) # use alg for improved inference
    solve_call(_prob, _alg, args...; kwargs...)
  elseif !isempty(args) && typeof(args[1]) <: DEAlgorithm
    alg = args[1]
    _alg = prepare_alg(alg, u0, p, prob)
    _prob = get_concrete_problem(prob, isadaptive(_alg); u0=u0, p=p, kwargs...)
    check_prob_alg_pairing(_prob, alg)
    solve_call(_prob, _alg, Base.tail(args)...; kwargs...)
  elseif isempty(args) # Default algorithm handling
    _prob = get_concrete_problem(prob, !(typeof(prob) <: DiscreteProblem); u0=u0, p=p, kwargs...)
    solve_call(_prob, args...; kwargs...)
  else
    _prob = get_concrete_problem(prob, !(typeof(prob) <: DiscreteProblem); u0=u0, p=p, kwargs...)
    solve_call(_prob, args...; kwargs...)
  end
end

function solve(prob::EnsembleProblem, args...; kwargs...)
  if isempty(args) || length(args) == 1 && typeof(args[1]) <: EnsembleAlgorithm
    __solve(prob, nothing, args...; kwargs...)
  else
    __solve(prob, args...; kwargs...)
  end
end

function solve(prob::AbstractNoiseProblem, args...; kwargs...)
  __solve(prob, args...; kwargs...)
end

function solve(prob::AbstractJumpProblem, args...; kwargs...)
  __solve(prob, args...; kwargs...)
end

function checkkwargs(kwargshandle; kwargs...)
  if any(x -> x ∉ allowedkeywords, keys(kwargs))
    if kwargshandle == KeywordArgError
      throw(CommonKwargError(kwargs))
    elseif kwargshandle == KeywordArgWarn
      @warn KWARGWARN_MESSAGE
    else
      @assert kwargshandle == KeywordArgSilent
    end
  end
end

@non_differentiable checkkwargs(kwargshandle)

function get_concrete_problem(prob::AbstractJumpProblem, isadapt; kwargs...)
  prob
end

function get_concrete_problem(prob::SteadyStateProblem, isadapt; kwargs...)
  u0 = get_concrete_u0(prob, isadapt, Inf, kwargs)
  u0 = promote_u0(u0, prob.p, nothing)
  remake(prob; u0=u0)
end

function get_concrete_problem(prob::NonlinearProblem, isadapt; kwargs...)
  u0 = get_concrete_u0(prob, isadapt, nothing, kwargs)
  u0 = promote_u0(u0, prob.p, nothing)
  remake(prob; u0=u0)
end

function get_concrete_problem(prob::AbstractEnsembleProblem, isadapt; kwargs...)
  prob
end

function solve(prob::PDEProblem, alg::DiffEqBase.DEAlgorithm, args...;
  kwargs...)
  solve(prob.prob, alg, args...; kwargs...)
end

function init(prob::PDEProblem, alg::DiffEqBase.DEAlgorithm, args...;
  kwargs...)
  init(prob.prob, alg, args...; kwargs...)
end

function get_concrete_problem(prob, isadapt; kwargs...)
  p = get_concrete_p(prob, kwargs)
  tspan = get_concrete_tspan(prob, isadapt, kwargs, p)
  u0 = get_concrete_u0(prob, isadapt, tspan[1], kwargs)
  u0_promote = promote_u0(u0, p, tspan[1])
  f_promote = promote_f(prob.f, u0_promote)
  tspan_promote = promote_tspan(u0_promote, p, tspan, prob, kwargs)
  if isconcreteu0(prob, tspan[1], kwargs) && typeof(u0_promote) === typeof(prob.u0) &&
     prob.tspan == tspan && typeof(prob.tspan) === typeof(tspan_promote) &&
     p === prob.p && f_promote === prob.f
    return prob
  else
    return remake(prob; f=f_promote, u0=u0_promote, p=p, tspan=tspan_promote)
  end
end

function get_concrete_problem(prob::DAEProblem, isadapt; kwargs...)

  p = get_concrete_p(prob, kwargs)
  tspan = get_concrete_tspan(prob, isadapt, kwargs, p)
  u0 = get_concrete_u0(prob, isadapt, tspan[1], kwargs)
  du0 = get_concrete_du0(prob, isadapt, tspan[1], kwargs)

  u0_promote = promote_u0(u0, p, tspan[1])
  du0_promote = promote_u0(du0, p, tspan[1])

  f_promote = promote_f(prob.f, u0_promote)
  tspan_promote = promote_tspan(u0_promote, p, tspan, prob, kwargs)
  if isconcreteu0(prob, tspan[1], kwargs) && typeof(u0_promote) === typeof(prob.u0) &&
     isconcretedu0(prob, tspan[1], kwargs) && typeof(du0_promote) === typeof(prob.du0) &&
     prob.tspan == tspan && typeof(prob.tspan) === typeof(tspan_promote) &&
     p === prob.p && f_promote === prob.f
    return prob
  else
    return remake(prob; f=f_promote, du0=du0_promote, u0=u0_promote, p=p, tspan=tspan_promote)
  end
end

function get_concrete_problem(prob::DDEProblem, isadapt; kwargs...)
  p = get_concrete_p(prob, kwargs)
  tspan = get_concrete_tspan(prob, isadapt, kwargs, p)
  u0 = get_concrete_u0(prob, isadapt, tspan[1], kwargs)

  if prob.constant_lags isa Function
    constant_lags = prob.constant_lags(p)
  else
    constant_lags = prob.constant_lags
  end

  u0 = promote_u0(u0, p, tspan[1])
  tspan = promote_tspan(u0, p, tspan, prob, kwargs)

  remake(prob; u0=u0, tspan=tspan, p=p, constant_lags=constant_lags)
end

function promote_f(f::F, u0) where {F}
  # Ensure our jacobian will be of the same type as u0
  uElType = u0 === nothing ? Float64 : eltype(u0)
  if isdefined(f, :jac_prototype) && f.jac_prototype isa AbstractArray
    f = @set f.jac_prototype = similar(f.jac_prototype, uElType)
  end
  return f
end

promote_f(f::SplitFunction, u0) = typeof(f.cache) === typeof(u0) && isinplace(f) ? f : remake(f, cache=zero(u0))
prepare_alg(alg, u0, p, f) = alg

function get_concrete_tspan(prob, isadapt, kwargs, p)
  if prob.tspan isa Function
    tspan = prob.tspan(p)
  elseif haskey(kwargs, :tspan)
    tspan = kwargs[:tspan]
  elseif prob.tspan === (nothing, nothing)
    throw(NoTspanError())
  else
    tspan = prob.tspan
  end

  isadapt && eltype(tspan) <: Integer && (tspan = float.(tspan))

  tspan
end

function isconcreteu0(prob, t0, kwargs)
  !eval_u0(prob.u0) && prob.u0 !== nothing && !isdistribution(prob.u0)
end

function isconcretedu0(prob, t0, kwargs)
  !eval_u0(prob.u0) && prob.du0 !== nothing && !isdistribution(prob.du0)
end

function get_concrete_u0(prob, isadapt, t0, kwargs)
  if eval_u0(prob.u0)
    u0 = prob.u0(prob.p, t0)
  elseif haskey(kwargs, :u0)
    u0 = kwargs[:u0]
  else
    u0 = prob.u0
  end

  isadapt && eltype(u0) <: Integer && (u0 = float.(u0))

  _u0 = handle_distribution_u0(u0)

  if isinplace(prob) && (_u0 isa Number || _u0 isa SArray)
    throw(IncompatibleInitialConditionError())
  end

  _u0
end

function get_concrete_du0(prob, isadapt, t0, kwargs)
  if eval_u0(prob.du0)
    du0 = prob.du0(prob.p, t0)
  elseif haskey(kwargs, :du0)
    du0 = kwargs[:du0]
  else
    du0 = prob.du0
  end

  isadapt && eltype(du0) <: Integer && (du0 = float.(du0))

  _du0 = handle_distribution_u0(du0)

  if isinplace(prob) && (_du0 isa Number || _du0 isa SArray)
    throw(IncompatibleInitialConditionError())
  end

  _du0
end

function get_concrete_p(prob, kwargs)
  if haskey(kwargs, :p)
    p = kwargs[:p]
  else
    p = prob.p
  end
end

handle_distribution_u0(_u0) = _u0
handle_distribution_u0(_u0::Distributions.Sampleable) = rand(_u0)
isdistribution(_u0::Distributions.Sampleable) = true

eval_u0(u0::Function) = true
eval_u0(u0) = false

function __solve(prob::DEProblem, args...; default_set=false, second_time=false, kwargs...)
  if second_time
    throw(NoDefaultAlgorithmError())
  elseif length(args) > 0 && !(typeof(args[1]) <: Union{Nothing,DEAlgorithm})
    throw(NonSolverError())
  else
    __solve(prob::DEProblem, nothing, args...; default_set=false, second_time=true, kwargs...)
  end
end

function check_prob_alg_pairing(prob, alg)
  if prob isa ODEProblem && !(alg isa AbstractODEAlgorithm) ||
     prob isa SDEProblem && !(alg isa AbstractSDEAlgorithm) ||
     prob isa SDDEProblem && !(alg isa AbstractSDEAlgorithm) ||
     prob isa DDEProblem && !(alg isa AbstractDDEAlgorithm) ||
     prob isa DAEProblem && !(alg isa AbstractDAEAlgorithm) ||
     prob isa SteadyStateProblem && !(alg isa AbstractSteadyStateAlgorithm)

     throw(ProblemSolverPairingError(prob, alg))
  end

  if isdefined(prob, :u0) && eltype(prob.u0) <: ForwardDiff.Dual && !SciMLBase.isautodifferentiable(alg)
    throw(DirectAutodiffError())
  end
end

################### Differentiation

struct SensitivityADPassThrough <: SciMLBase.DEAlgorithm end

function ChainRulesCore.frule(::typeof(solve_up), prob,
  sensealg::Union{Nothing,AbstractSensitivityAlgorithm},
  u0, p, args...;
  kwargs...)
  _solve_forward(prob, sensealg, u0, p, args...; kwargs...)
end

function ChainRulesCore.rrule(::typeof(solve_up), prob::SciMLBase.DEProblem,
  sensealg::Union{Nothing,AbstractSensitivityAlgorithm},
  u0, p, args...;
  kwargs...)
  _solve_adjoint(prob, sensealg, u0, p, args...; kwargs...)
end

###
### Legacy Dispatches to be Non-Breaking
###

@deprecate concrete_solve(prob::SciMLBase.DEProblem, alg::Union{SciMLBase.DEAlgorithm,Nothing},
  u0=prob.u0, p=prob.p, args...; kwargs...) solve(prob, alg, args...; u0=u0, p=p, kwargs...)

function _solve_adjoint(prob, sensealg, u0, p, args...; merge_callbacks=true, kwargs...)

  _prob = if haskey(kwargs, :alg) && (isempty(args) || args[1] === nothing)
    alg = kwargs[:alg]
    get_concrete_problem(prob, isadaptive(alg); u0=u0, p=p, kwargs...)
  elseif !isempty(args) && typeof(args[1]) <: DEAlgorithm
    alg = args[1]
    get_concrete_problem(prob, isadaptive(alg); u0=u0, p=p, kwargs...)
  elseif isempty(args) # Default algorithm handling
    get_concrete_problem(prob, !(typeof(prob) <: DiscreteProblem); u0=u0, p=p, kwargs...)
  else
    get_concrete_problem(prob, !(typeof(prob) <: DiscreteProblem); u0=u0, p=p, kwargs...)
  end

  if has_kwargs(_prob)
    if merge_callbacks && haskey(_prob.kwargs, :callback) && haskey(kwargs, :callback)
      kwargs_temp = NamedTuple{Base.diff_names(Base._nt_names(
          values(kwargs)), (:callback,))}(values(kwargs))
      callbacks = NamedTuple{(:callback,)}((DiffEqBase.CallbackSet(_prob.kwargs[:callback], values(kwargs).callback),))
      kwargs = merge(kwargs_temp, callbacks)
    end
    kwargs = isempty(_prob.kwargs) ? kwargs : merge(values(_prob.kwargs), kwargs)
  end

  if isempty(args)
    _concrete_solve_adjoint(_prob, nothing, sensealg, u0, p; kwargs...)
  else
    _concrete_solve_adjoint(_prob, args[1], sensealg, u0, p, Base.tail(args)...; kwargs...)
  end
end

function _solve_forward(prob, sensealg, u0, p, args...; merge_callbacks=true, kwargs...)

  _prob = if haskey(kwargs, :alg) && (isempty(args) || args[1] === nothing)
    alg = kwargs[:alg]
    get_concrete_problem(prob, isadaptive(alg); u0=u0, p=p, kwargs...)
  elseif !isempty(args) && typeof(args[1]) <: DEAlgorithm
    alg = args[1]
    get_concrete_problem(prob, isadaptive(alg); u0=u0, p=p, kwargs...)
  elseif isempty(args) # Default algorithm handling
    get_concrete_problem(prob, !(typeof(prob) <: DiscreteProblem); u0=u0, p=p, kwargs...)
  else
    get_concrete_problem(prob, !(typeof(prob) <: DiscreteProblem); u0=u0, p=p, kwargs...)
  end

  if has_kwargs(_prob)
    if merge_callbacks && haskey(_prob.kwargs, :callback) && haskey(kwargs, :callback)
      kwargs_temp = NamedTuple{Base.diff_names(Base._nt_names(
          values(kwargs)), (:callback,))}(values(kwargs))
      callbacks = NamedTuple{(:callback,)}((DiffEqBase.CallbackSet(_prob.kwargs[:callback], values(kwargs).callback),))
      kwargs = merge(kwargs_temp, callbacks)
    end
    kwargs = isempty(_prob.kwargs) ? kwargs : merge(values(_prob.kwargs), kwargs)
  end

  if isempty(args)
    _concrete_solve_forward(prob, nothing, sensealg, u0, p; kwargs...)
  else
    _concrete_solve_forward(prob, args[1], sensealg, u0, p, Base.tail(args)...; kwargs...)
  end
end

function _concrete_solve_adjoint(args...; kwargs...)
  error("No adjoint rules exist. Check that you added `using DiffEqSensitivity`")
end

function _concrete_solve_forward(args...; kwargs...)
  error("No sensitivity rules exist. Check that you added `using DiffEqSensitivity`")
end

max_stars_repo_name: JuliaDiffEq/DifferentialEquationsBase

max_issues_repo_name: JuliaDiffEq/DifferentialEquationsBase

max_forks_repo_name: JuliaDiffEq/DifferentialEquationsBase


text:
#MAIN ENTRY POINT
#final columns for the decompressed file

#merges crsp and compustat
# note refresh merge will be automatically run if any of the other two are run
function loadccm(;
  datapath::String=DATA_PATH,
  ccmname::String = CCM_NAME,
  incsvstream::Function = IN_CSV_STREAM)

  local ccm::DataFrame

  #makes a serialization object so we don't have to keep parsing the CSV
  ccm = incsvstream("$datapath\\$ccmname.csv.lz4") |> CSV.read |> DataFrame

  return ccm

end

function mergecompccm!(comp::DataFrame;
    testoutput::Bool = TEST_OUTPUT)

  #first
  local ccm::DataFrame = loadccm()
  ccm = preprocessccm!(ccm)

  testoutput && CSV.write("output\\ccm.csv", ccm)

  #println("Comp rows before merge: $(size(comp, 1))")
  comp = join(comp, ccm, on=:gvkey, kind=:inner)

  #fmatch the dates
  #println("Comp rows after merge, pre-filter: $(size(comp, 1))")
  comp.keep = trues(size(comp,1))

  comp.keep = (r::DataFrameRow->
    r.adateq ∈ r.linkeffdate:Day(1):r.linkenddate).(eachrow(comp))

  comp = comp[comp.keep,:]
  select!(comp, Not(:keep))
  println("Comp rows after merge, post-filter: $(size(comp, 1))")

  comp = reconcilecompccmintervals!(comp)

  return comp
end

#formats the ccm table in a reasonable way
function preprocessccm!(ccm::DataFrame;
  retainedcolsccm::Vector{Symbol} = RETAINED_COLS_CCM,
  lastdatestring::String = LAST_DATE_STRING)

  #names ot lower case
  names!(ccm, (s::Symbol->Symbol(lowercase(string(s)))).(names(ccm)))

  #fix the end date
  ccm.linkenddt .= (s::MString->
    coalesce(s,"")=="E" ? lastdatestring : s).(ccm.linkenddt)

  filter!(r::DataFrameRow->r.linkprim =='P' || r.linkprim =='C', ccm)

  #parse the dates
  ccm.linkeffdate = (s->parseccm(Date, s)).(ccm.linkdt)
  ccm.linkenddate = (s->parseccm(Date, s)).(ccm.linkenddt)

  #only need some of the fields
  ccm = select!(ccm, retainedcolsccm)

  return ccm
end

#this function reconciles the link intervals with the announcement intervals
#see
function reconcilecompccmintervals!(comp::DataFrame)
  #subtract one day since the effective interval is inclusive of the
  #first date while the begin-end date is exclusve of the first date
  comp.begindate .= (max).(comp.begindate, comp.linkeffdate .- Day(1))
  comp.enddate .= (min).(comp.enddate, comp.linkenddate)

  #this gets rid of one day records
  #(this assumes the previous effective record is the current record)
  #not many of these ~140
  #println("records before interval reconciliation: $(size(comp,1))")
  #println(comp[comp.begindate .≥ comp.enddate,
  #  [:gvkey, :adateq, :begindate, :enddate, :linkeffdate, :linkenddate]])
  comp = comp[comp.begindate .< comp.enddate, :]
  println("records after interval reconciliation: $(size(comp,1))")

  return comp
end

function parseccm(::Type{T}, s::String;
    #=parsedmissings = CCM_PARSED_MISSINGS=#) where T

  #if it doesn't parse to the right type, set to missing
  v::Union{T,Missing} = something(tryparse(T, s), missing)

  #check against the list of missing codes
  #(!ismissing(v)) && (v ∈ parsedmissings) && (v=missing)

  return v
end

#the date case
parseccm(::Type{Date}, s::String;
    ccmdateformat::DateFormat = CCM_DATE_FORMAT) = Dates.Date(s, ccmdateformat)

#helper methods and special cases
parseccm(::Type{<:Any}, v::Missing) = missing
parseccm(::Type{Date}, i::Int) = parseccm(Date, "$i")

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using NamedTypeParameters
using Test

# Test types
struct NoTypes end
struct AnyTypes1{A} end
struct AnyTypes2{A,B} end
struct AnyTypes3{A,B,C} end
struct SuperTypes1{A<:Real} end
struct SuperTypes2{A<:Real,B<:Array} end
struct SuperTypes3{A<:Real,B<:Array,C<:AbstractDict} end

@testset begin

### Notes on parameter types

# specifying a type with <: equals adding an unnamed where parameter, e.g.
@test SuperTypes1{<:Real} == SuperTypes1{A} where A<:Real
# note: using invalid where supertypes (invalid wrt the definition of the parameteric type)
# does not throw
SuperTypes1{A} where A<:String # works, although SuperTypes1 is defined with A<:Real
# but inserting an invalid type as a concrete type parameter does throw
@test_throws TypeError SuperTypes1{String}

### Pkg tests

# return type itself if no parameters are given
types = [ NoTypes,
          AnyTypes1, AnyTypes2, AnyTypes2,
          SuperTypes1, SuperTypes2, SuperTypes3 ]
for t in types
  @test @parameterize(t) == t
end

# using the parameter type name alone inserts its default supertype
@test @parameterize(AnyTypes1{A}) == AnyTypes1{<:Any}
@test @parameterize(AnyTypes2{A,B}) == AnyTypes2{<:Any, <:Any}
@test @parameterize(AnyTypes3{A,B,C}) == AnyTypes3{<:Any, <:Any, <:Any}
@test @parameterize(SuperTypes1{A}) == SuperTypes1{<:Real}
@test @parameterize(SuperTypes2{A,B}) == SuperTypes2{<:Real, <:Array}
@test @parameterize(SuperTypes3{A,B,C}) == SuperTypes3{<:Real, <:Array, <:AbstractDict}

# override type parameter defaults
@test @parameterize(AnyTypes1{A<:Float64}) == AnyTypes1{<:Float64}
@test @parameterize(AnyTypes2{A<:Float64,B<:String}) == AnyTypes2{<:Float64,<:String}
@test @parameterize(AnyTypes3{A<:Float64,B<:String,C<:Symbol}) == AnyTypes3{<:Float64,<:String,<:Symbol}
@test @parameterize(SuperTypes1{A<:Float64}) == SuperTypes1{<:Float64}
@test @parameterize(SuperTypes2{A<:Float64,B<:Vector}) == SuperTypes2{<:Float64,<:Vector}
@test @parameterize(SuperTypes3{A<:Float64,B<:Vector,C<:Dict}) == SuperTypes3{<:Float64,<:Vector,<:Dict}

# skip parameters
@test @parameterize(AnyTypes2{A<:Float64}) == AnyTypes2{<:Float64,<:Any}
@test @parameterize(AnyTypes2{B<:String}) == AnyTypes2{<:Any,<:String}
@test @parameterize(AnyTypes3{A<:Float64}) == AnyTypes3{<:Float64,<:Any,<:Any}
@test @parameterize(AnyTypes3{B<:Vector}) == AnyTypes3{<:Any,<:Vector,<:Any}
@test @parameterize(AnyTypes3{C<:Dict}) == AnyTypes3{<:Any,<:Any,<:Dict}
@test @parameterize(AnyTypes3{A<:Float64,B<:String}) == AnyTypes3{<:Float64,<:String,<:Any}
@test @parameterize(AnyTypes3{A<:Float64,C<:Symbol}) == AnyTypes3{<:Float64,<:Any,<:Symbol}
@test @parameterize(AnyTypes3{B<:String,C<:Symbol}) == AnyTypes3{<:Any,<:String,<:Symbol}
@test @parameterize(SuperTypes2{A<:Float64}) == SuperTypes2{<:Float64,<:Array}
@test @parameterize(SuperTypes2{B<:Vector}) == SuperTypes2{<:Real,<:Vector}
@test @parameterize(SuperTypes3{A<:Float64}) == SuperTypes3{<:Float64,<:Array,<:AbstractDict}
@test @parameterize(SuperTypes3{B<:Vector}) == SuperTypes3{<:Real,<:Vector,<:AbstractDict}
@test @parameterize(SuperTypes3{C<:Dict}) == SuperTypes3{<:Real,<:Array,<:Dict}
@test @parameterize(SuperTypes3{A<:Float64,B<:Vector}) == SuperTypes3{<:Float64,<:Vector,<:AbstractDict}
@test @parameterize(SuperTypes3{A<:Float64,C<:Dict}) == SuperTypes3{<:Float64,<:Array,<:Dict}
@test @parameterize(SuperTypes3{B<:Vector,C<:Dict}) == SuperTypes3{<:Real,<:Vector,<:Dict}

# throw if more overrides than type parameters are given
@test_throws LoadError @eval @parameterize(NoTypes{A})
@test_throws LoadError @eval @parameterize(AnyTypes1{A,B})
@test_throws LoadError @eval @parameterize(AnyTypes2{A,B,C})
@test_throws LoadError @eval @parameterize(AnyTypes3{A,B,C,D})
@test_throws LoadError @eval @parameterize(SuperTypes1{A,B})
@test_throws LoadError @eval @parameterize(SuperTypes2{A,B,C})
@test_throws LoadError @eval @parameterize(SuperTypes3{A,B,C,D})

# set a parameter type with =
@test @parameterize(AnyTypes1{A=Float64}) == AnyTypes1{Float64}
@test @parameterize(AnyTypes2{A=Float64,B=String}) == AnyTypes2{Float64,String}
@test @parameterize(AnyTypes3{A=Float64,B=String,C=Symbol}) == AnyTypes3{Float64,String,Symbol}
@test @parameterize(SuperTypes1{A=Float64}) == SuperTypes1{Float64}
@test @parameterize(SuperTypes2{A=Float64,B=Vector}) == SuperTypes2{Float64,Vector}
@test @parameterize(SuperTypes3{A=Float64,B=Vector,C=Dict}) == SuperTypes3{Float64,Vector,Dict}
# can also be used with where syntax
# note that AnyTypes1{<:Any} == AnyTypes1{T} where <: Any
@test (@parameterize(AnyTypes1{A=T}) where {T<:Real}) == AnyTypes1{<:Real}
@test (@parameterize(AnyTypes2{A=T,B=S}) where {T<:Real,S<:Array}) == AnyTypes2{<:Real, <:Array}
@test (@parameterize(AnyTypes3{A=T,B=S,C=R}) where {T<:Real,S<:Array,R<:Symbol}) == AnyTypes3{<:Real, <:Array, <:Symbol}
@test (@parameterize(SuperTypes1{A=T}) where {T<:Float64}) == SuperTypes1{<:Float64}
@test (@parameterize(SuperTypes2{A=T,B=S}) where {T<:Float64,S<:Vector}) == SuperTypes2{<:Float64, <:Vector}
@test (@parameterize(SuperTypes3{A=T,B=S,C=R}) where {T<:Float64,S<:Vector,R<:Dict}) == SuperTypes3{<:Float64, <:Vector, <:Dict}

# duplicated parameter names
@test_throws LoadError @eval @parameterize(AnyTypes2{A,A})
@test_throws LoadError @eval @parameterize(AnyTypes3{A,A,A})
@test_throws LoadError @eval @parameterize(AnyTypes3{A,A,C})
@test_throws LoadError @eval @parameterize(AnyTypes3{A,B,A})
@test_throws LoadError @eval @parameterize(AnyTypes3{A,B,B})
@test_throws LoadError @eval @parameterize(SuperTypes2{A,A})
@test_throws LoadError @eval @parameterize(SuperTypes3{A,A,A})
@test_throws LoadError @eval @parameterize(SuperTypes3{A,A,C})
@test_throws LoadError @eval @parameterize(SuperTypes3{A,B,A})
@test_throws LoadError @eval @parameterize(SuperTypes3{A,B,B})

# unknown parameter names
@test_throws LoadError @eval @parameterize(AnyTypes1{X})
@test_throws LoadError @eval @parameterize(AnyTypes2{X})
@test_throws LoadError @eval @parameterize(AnyTypes2{X,Y})
@test_throws LoadError @eval @parameterize(AnyTypes3{X})
@test_throws LoadError @eval @parameterize(AnyTypes3{X,Y})
@test_throws LoadError @eval @parameterize(AnyTypes3{X,Y,Z})
@test_throws LoadError @eval @parameterize(AnyTypes3{A,Y})
@test_throws LoadError @eval @parameterize(AnyTypes3{A,B,Z})
@test_throws LoadError @eval @parameterize(SuperTypes1{X})
@test_throws LoadError @eval @parameterize(SuperTypes2{X})
@test_throws LoadError @eval @parameterize(SuperTypes2{X,Y})
@test_throws LoadError @eval @parameterize(SuperTypes3{X})
@test_throws LoadError @eval @parameterize(SuperTypes3{X,Y})
@test_throws LoadError @eval @parameterize(SuperTypes3{X,Y,Z})
@test_throws LoadError @eval @parameterize(SuperTypes3{A,Y})
@test_throws LoadError @eval @parameterize(SuperTypes3{A,B,Z})

end # testset


text:
using Revise
using PyCall
using PyPlot
using Bem2d

"""
    plotlocal(els, idx, dispfault, dispfreesurfflat, dispfreesurftopo, stressfault, stressfreesurfflat, stressfreesurftopo, xobs, yobs, npts)

Generate four panel plots to shw displacements and stresses for the
flat and topographic cases.
"""
function plotlocal(els, idx, dispfault, dispfreesurfflat, dispfreesurftopo, stressfault, stressfreesurfflat, stressfreesurftopo, xobs, yobs, npts)
    # Pretty picture of displacements and stresses
    xfreesurf = unique([els.x1[idx["freesurfflat"]] ; els.x2[idx["freesurfflat"]]])
    xfill = [xfreesurf ; [10e3 ; -10e3 ; -10e3]]
    yfreesurf = unique([els.y1[idx["freesurfflat"]] ; els.y2[idx["freesurfflat"]]])
    yfill = [yfreesurf ; [5e3 ; 5e3 ; minimum(yfreesurf)]]

    # Combine two fields for total displacement and stress fields
    dispfield = sqrt.((dispfault - dispfreesurfflat)[:, 1].^2 + (dispfault - dispfreesurfflat)[:, 2].^2)
    stressxx = (stressfault - stressfreesurfflat)[:, 1]
    stressyy = (stressfault - stressfreesurfflat)[:, 2]
    stressxy = (stressfault - stressfreesurfflat)[:, 3]
    I1 = stressxx + stressyy  # 1st invariant
    I2 = stressxx .* stressyy - stressxy.^2  # 2nd invariant
    J2 = (I1.^2) ./ 3.0 - I2  # 2nd invariant (deviatoric)
    stressfield = log10.(abs.(J2))

    ncontours = 20
    dispcontours = collect(LinRange(0.0, 1.0, ncontours))
    stresscontours = collect(LinRange(10.0, 16.0, ncontours))

    figure(figsize = (15, 10))
    subplot(2, 2, 1)
    contourf(reshape(xobs, npts, npts), reshape(yobs, npts, npts), reshape(dispfield, npts, npts), dispcontours, cmap = get_cmap("plasma"))
    colorbar(fraction = 0.020, pad = 0.05, extend = "both", label = L"$||u||$ (m)")
    contour(reshape(xobs, npts, npts), reshape(yobs, npts, npts), reshape(dispfield, npts, npts), dispcontours, linewidths = 0.5, colors = "gray")
    fill(xfill, yfill, "w", zorder = 30)
    for i in 1:els.endidx
        if els.name[i] == "fault" || els.name[i] == "freesurfflat"
            plot([els.x1[i], els.x2[i]], [els.y1[i], els.y2[i]], "-k", linewidth = 2.0, zorder=40)
        end
    end
    xlim([minimum(xobs), maximum(xobs)])
    ylim([minimum(yobs), maximum(yobs)])
    xticks([minimum(xobs), 0, maximum(xobs)])
    yticks([minimum(yobs), 0, maximum(yobs)])
    gca().set_aspect("equal")
    xlabel(L"$x$ (m)")
    ylabel(L"$y$ (m)")

    subplot(2, 2, 3)
    contourf(reshape(xobs, npts, npts), reshape(yobs, npts, npts), reshape(stressfield, npts, npts), stresscontours, cmap = get_cmap("hot_r"))
    colorbar(fraction = 0.020, pad = 0.05, extend = "both", label = L"$\log_{10}|\mathrm{J}_2|$ (Pa$^2$)")
    contour(reshape(xobs, npts, npts), reshape(yobs, npts, npts), reshape(stressfield, npts, npts), stresscontours, linewidths = 0.5, colors = "gray")
    fill(xfill, yfill, "w", zorder = 30)
    for i in 1:els.endidx
        if els.name[i] == "fault" || els.name[i] == "freesurfflat"
            plot([els.x1[i], els.x2[i]], [els.y1[i], els.y2[i]], "-k", linewidth = 2.0, zorder=40)
        end
    end
    xlim([minimum(xobs), maximum(xobs)])
    ylim([minimum(yobs), maximum(yobs)])
    xticks([minimum(xobs), 0, maximum(xobs)])
    yticks([minimum(yobs), 0, maximum(yobs)])
    gca().set_aspect("equal")
    xlabel(L"$x$ (m)")
    ylabel(L"$y$ (m)")

    # Topography
    xfreesurf = unique([els.x1[idx["freesurftopo"]] ; els.x2[idx["freesurftopo"]]])
    xfill = [xfreesurf ; [10e3 ; -10e3 ; -10e3]]
    yfreesurf = unique([els.y1[idx["freesurftopo"]] ; els.y2[idx["freesurftopo"]]])
    yfill = [yfreesurf ; [5e3 ; 5e3 ; minimum(yfreesurf)]]

    # Combine two fields for total displacement and stress fields
    dispfield = sqrt.((dispfault - dispfreesurftopo)[:, 1].^2 + (dispfault - dispfreesurftopo)[:, 2].^2)
    stressxx = (stressfault - stressfreesurftopo)[:, 1]
    stressyy = (stressfault - stressfreesurftopo)[:, 2]
    stressxy = (stressfault - stressfreesurftopo)[:, 3]
    I1 = stressxx + stressyy  # 1st invariant
    I2 = stressxx .* stressyy - stressxy.^2  # 2nd invariant
    J2 = (I1.^2) ./ 3.0 - I2  # 2nd invariant (deviatoric)
    stressfield = log10.(abs.(J2))

    subplot(2, 2, 2)
    contourf(reshape(xobs, npts, npts), reshape(yobs, npts, npts), reshape(dispfield, npts, npts), dispcontours, cmap = get_cmap("plasma"))
    colorbar(fraction = 0.020, pad = 0.05, extend = "both", label = L"$||u||$ (m)")
    contour(reshape(xobs, npts, npts), reshape(yobs, npts, npts), reshape(dispfield, npts, npts), dispcontours, linewidths = 0.5, colors = "gray")
    fill(xfill, yfill, "w", zorder = 30)
    for i in 1:els.endidx
        if els.name[i] == "fault" || els.name[i] == "freesurftopo"
            plot([els.x1[i], els.x2[i]], [els.y1[i], els.y2[i]], "-k", linewidth = 2.0, zorder=40)
        end
    end
    xlim([minimum(xobs), maximum(xobs)])
    ylim([minimum(yobs), maximum(yobs)])
    xticks([minimum(xobs), 0, maximum(xobs)])
    yticks([minimum(yobs), 0, maximum(yobs)])
    gca().set_aspect("equal")
    xlabel(L"$x$ (m)")
    ylabel(L"$y$ (m)")

    subplot(2, 2, 4)
    contourf(reshape(xobs, npts, npts), reshape(yobs, npts, npts), reshape(stressfield, npts, npts), stresscontours, cmap = get_cmap("hot_r"))
    colorbar(fraction = 0.020, pad = 0.05, extend = "both", label = L"$\log_{10}|\mathrm{J}_2|$ (Pa$^2$)")
    contour(reshape(xobs, npts, npts), reshape(yobs, npts, npts), reshape(stressfield, npts, npts), stresscontours, linewidths = 0.5, colors = "gray")
    fill(xfill, yfill, "w", zorder = 30)
    for i in 1:els.endidx
        if els.name[i] == "fault" || els.name[i] == "freesurftopo"
            plot([els.x1[i], els.x2[i]], [els.y1[i], els.y2[i]], "-k", linewidth = 2.0, zorder=40)
        end
    end
    xlim([minimum(xobs), maximum(xobs)])
    ylim([minimum(yobs), maximum(yobs)])
    xticks([minimum(xobs), 0, maximum(xobs)])
    yticks([minimum(yobs), 0, maximum(yobs)])
    gca().set_aspect("equal")
    xlabel(L"$x$ (m)")
    ylabel(L"$y$ (m)")
    tight_layout()
    show()
    return nothing
end

"""
    thrustfaulttopography()

Comparison of volume displacemnts and stresses from a fault with
and without topography.  Includes both constant and quadratic
element cases.
"""
function thrustfaulttopography()
    close("all")
    mu = 30e9
    nu = 0.25
    els = Elements(Int(1e5))

    # Observation points for internal evaluation and visualization
    npts = 50
    xobs, yobs = obsgrid(-10e3, -5e3, 10e3, 5e3, npts)
    nfreesurf = 20
    nfault = 10

    # Flat free surface
    x1, y1, x2, y2 = discretizedline(-20e3, 0, 20e3, 0, nfreesurf)
    y1 = -1e-3 * atan.(x1 / 1e3)
    y2 = -1e-3 * atan.(x2 / 1e3)
    for i in 1:length(x1)
        els.x1[els.endidx + i] = x1[i]
        els.y1[els.endidx + i] = y1[i]
        els.x2[els.endidx + i] = x2[i]
        els.y2[els.endidx + i] = y2[i]
        els.name[els.endidx + i] = "freesurfflat"
    end
    standardize_elements!(els)

    # Topographic free surface
    x1, y1, x2, y2 = discretizedline(-20e3, 0, 20e3, 0, nfreesurf)
    y1 = -1e3 * atan.(x1 / 1e3)
    y2 = -1e3 * atan.(x2 / 1e3)
    for i in 1:length(x1)
        els.x1[els.endidx + i] = x1[i]
        els.y1[els.endidx + i] = y1[i]
        els.x2[els.endidx + i] = x2[i]
        els.y2[els.endidx + i] = y2[i]
        els.name[els.endidx + i] = "freesurftopo"
    end
    standardize_elements!(els)

    # Curved fault
    x1, y1, x2, y2 = discretizedline(-7e3, 0e3, 0, 0, nfault)
    y1 = 3e3 * atan.(x1 / 1e3)
    y2 = 3e3 * atan.(x2 / 1e3)
    for i in 1:length(x1)
        els.x1[els.endidx + i] = x1[i]
        els.y1[els.endidx + i] = y1[i]
        els.x2[els.endidx + i] = x2[i]
        els.y2[els.endidx + i] = y2[i]
        els.name[els.endidx + i] = "fault"
    end
    standardize_elements!(els)

    # Create convience tools
    idx = getidxdict(els)
    partialsconst = initpartials(els)
    partialsquad = initpartials(els)

    #! CS elements
    # Partial derivatves for solving BEM problem
    _, _, partialsconst["trac"]["fault"]["freesurfflat"] = partialsconstdispstress(slip2dispstress, els, idx["fault"], idx["freesurfflat"], mu, nu)
    _, _, partialsconst["trac"]["freesurfflat"]["freesurfflat"] = partialsconstdispstress(slip2dispstress, els, idx["freesurfflat"], idx["freesurfflat"], mu, nu)
    _, _, partialsconst["trac"]["fault"]["freesurftopo"] = partialsconstdispstress(slip2dispstress, els, idx["fault"], idx["freesurftopo"], mu, nu)
    _, _, partialsconst["trac"]["freesurftopo"]["freesurftopo"] = partialsconstdispstress(slip2dispstress, els, idx["freesurftopo"], idx["freesurftopo"], mu, nu)

    # Solve the BEM problem for unit slip in the x-direction
    faultslip = zeros(2 * nfault)
    faultslip[1:2:end] .= 1.0 # Global coordinate system
    dispfreesurfflatslip = inv(partialsconst["trac"]["freesurfflat"]["freesurfflat"]) * partialsconst["trac"]["fault"]["freesurfflat"] * faultslip
    dispfreesurftoposlip = inv(partialsconst["trac"]["freesurftopo"]["freesurftopo"]) * partialsconst["trac"]["fault"]["freesurftopo"] * faultslip

    # Forward model for volumetric displacements and stresses
    dispfault, stressfault = constdispstress(slip2dispstress, xobs, yobs, els, idx["fault"], ones(size(idx["fault"])), zeros(size(idx["fault"])), mu, nu)
    dispfreesurfflat, stressfreesurfflat = constdispstress(slip2dispstress, xobs, yobs, els, idx["freesurfflat"], dispfreesurfflatslip[1:2:end], dispfreesurfflatslip[2:2:end], mu, nu)
    dispfreesurftopo, stressfreesurftopo = constdispstress(slip2dispstress, xobs, yobs, els, idx["freesurftopo"], dispfreesurftoposlip[1:2:end], dispfreesurftoposlip[2:2:end], mu, nu)

    # Pretty plotting
    plotlocal(els, idx, dispfault, dispfreesurfflat, dispfreesurftopo, stressfault, stressfreesurfflat, stressfreesurftopo, xobs, yobs, npts)

    #! 3QN elements
    # Partial derivatves for solving BEM problem
    _, _, partialsquad["trac"]["fault"]["freesurfflat"] = partialsquaddispstress(slip2dispstress, els, idx["fault"], idx["freesurfflat"], mu, nu)
    _, _, partialsquad["trac"]["freesurfflat"]["freesurfflat"] = partialsquaddispstress(slip2dispstress, els, idx["freesurfflat"], idx["freesurfflat"], mu, nu)
    _, _, partialsquad["trac"]["fault"]["freesurftopo"] = partialsquaddispstress(slip2dispstress, els, idx["fault"], idx["freesurftopo"], mu, nu)
    _, _, partialsquad["trac"]["freesurftopo"]["freesurftopo"] = partialsquaddispstress(slip2dispstress, els, idx["freesurftopo"], idx["freesurftopo"], mu, nu)

    # Solve the BEM problem for unit slip in the x-direction
    faultslip = zeros(6 * nfault)
    faultslip[1:2:end] .= 1.0 # Global coordinate system
    dispfreesurfflatslip = inv(partialsquad["trac"]["freesurfflat"]["freesurfflat"]) * partialsquad["trac"]["fault"]["freesurfflat"] * faultslip
    dispfreesurftoposlip = inv(partialsquad["trac"]["freesurftopo"]["freesurftopo"]) * partialsquad["trac"]["fault"]["freesurftopo"] * faultslip

    # Forward model for volumetric displacements and stresses
    dispfault, stressfault = quaddispstress(slip2dispstress, xobs, yobs, els, idx["fault"], quadstack(ones(3 * length(idx["fault"]))), quadstack(zeros(3 * length(idx["fault"]))), mu, nu)
    dispfreesurfflat, stressfreesurfflat = quaddispstress(slip2dispstress, xobs, yobs, els, idx["freesurfflat"], quadstack(dispfreesurfflatslip[1:2:end]), quadstack(dispfreesurfflatslip[2:2:end]), mu, nu)
    dispfreesurftopo, stressfreesurftopo = quaddispstress(slip2dispstress, xobs, yobs, els, idx["freesurftopo"], quadstack(dispfreesurftoposlip[1:2:end]), quadstack(dispfreesurftoposlip[2:2:end]), mu, nu)

    # Pretty plotting
    plotlocal(els, idx, dispfault, dispfreesurfflat, dispfreesurftopo, stressfault, stressfreesurfflat, stressfreesurftopo, xobs, yobs, npts)

end
thrustfaulttopography()


text:
using Test
using DistributedData, Distributed, Random

@testset "DistributedData tests" begin
    include("base.jl")
    include("tools.jl")
end


text:
## Destructive matrix exponential using algorithm from Higham, 2008,
## "Functions of Matrices: Theory and Computation", SIAM
##
## Non-allocating version of `LinearAlgebra.exp!`. Modifies `A` to
## become (approximately) `exp(A)`.
function _exp!(A::StridedMatrix{T}; caches=nothing) where T <: LinearAlgebra.BlasFloat
    X = A
    n = LinearAlgebra.checksquare(A)
    # if ishermitian(A)
        # return copytri!(parent(exp(Hermitian(A))), 'U', true)
    # end

    if caches == nothing
        A2   = Matrix{T}(undef, n, n)
        P    = Matrix{T}(undef, n, n)
        U    = Matrix{T}(undef, n, n)
        V    = Matrix{T}(undef, n, n)
        temp = Matrix{T}(undef, n, n)
    else
        A2, P, U, V, temp = caches
    end
    fill!(P, zero(T)); fill!(@diagview(P), one(T)) # P = Inn

    ilo, ihi, scale = LAPACK.gebal!('B', A)    # modifies A
    nA = opnorm(A, 1)
    ## For sufficiently small nA, use lower order Padé-Approximations
    if (nA <= 2.1)
        if nA > 0.95
            C = T[17643225600.,8821612800.,2075673600.,302702400.,
                     30270240.,   2162160.,    110880.,     3960.,
                           90.,         1.]
        elseif nA > 0.25
            C = T[17297280.,8648640.,1995840.,277200.,
                     25200.,   1512.,     56.,     1.]
        elseif nA > 0.015
            C = T[30240.,15120.,3360.,
                    420.,   30.,   1.]
        else
            C = T[120.,60.,12.,1.]
        end
        mul!(A2, A, A)
        @. U = C[2] * P
        @. V = C[1] * P
        for k in 1:(div(size(C, 1), 2) - 1)
            k2 = 2 * k
            mul!(temp, P, A2); P, temp = temp, P # equivalent to P *= A2
            @. U += C[k2 + 2] * P
            @. V += C[k2 + 1] * P
        end
        mul!(temp, A, U); U, temp = temp, U # equivalent to U = A * U
        @. X = V + U
        @. temp = V - U
        LAPACK.gesv!(temp, X)
    else
        s  = log2(nA/5.4)               # power of 2 later reversed by squaring
        if s > 0
            si = ceil(Int,s)
            A ./= convert(T,2^si)
        end
        C  = T[64764752532480000.,32382376266240000.,7771770303897600.,
                1187353796428800.,  129060195264000.,  10559470521600.,
                    670442572800.,      33522128640.,      1323241920.,
                        40840800.,           960960.,           16380.,
                             182.,                1.]
        mul!(A2, A, A)
        @. U = C[2] * P
        @. V = C[1] * P
        for k in 1:6
            k2 = 2 * k
            mul!(temp, P, A2); P, temp = temp, P # equivalent to P *= A2
            @. U += C[k2 + 2] * P
            @. V += C[k2 + 1] * P
        end
        mul!(temp, A, U); U, temp = temp, U # equivalent to U = A * U
        @. X = V + U
        @. temp = V - U
        LAPACK.gesv!(temp, X)

        if s > 0            # squaring to reverse dividing by power of 2
            for t=1:si
                mul!(temp, X, X)
                X .= temp
            end
        end
    end

    # Undo the balancing
    for j = ilo:ihi
        scj = scale[j]
        for i = 1:n
            X[j,i] *= scj
        end
        for i = 1:n
            X[i,j] /= scj
        end
    end

    if ilo > 1       # apply lower permutations in reverse order
        for j in (ilo-1):-1:1; LinearAlgebra.rcswap!(j, Int(scale[j]), X) end
    end
    if ihi < n       # apply upper permutations in forward order
        for j in (ihi+1):n;    LinearAlgebra.rcswap!(j, Int(scale[j]), X) end
    end
    return X
end


text:
# This file is a part of SimilaritySearch.jl

using Test, SimilaritySearch, LinearAlgebra

@testset "closestpair" begin
    dist = CosineDistance()
    dim, mindist = 2, 1e-4
    db = MatrixDatabase(rand(Float32, dim, 1000))
    G = SearchGraph(; db, dist)
    tG = @elapsed index!(G)
    tG += @elapsed i, j, d = closestpair(G)
    @test i != j
    @test d < mindist
    @show i, j, d
    i, j, d = closestpair(G; parallel=true)
    @test i != j
    @test d < mindist
    @show i, j, d, :parallel
    tE = @elapsed i, j, d = closestpair(ExhaustiveSearch(; dist, db))
    @info "NOTE: the exact method will be faster on small datasets due to the preprocessing step of the approximation method"
    @info "closestpair computation time", :approx => tG, :exact => tE
end

max_stars_repo_name: sadit/SimilaritySearch

max_issues_repo_name: sadit/SimilaritySearch

max_forks_repo_name: sadit/SimilaritySearch


text:
function test_nzrows(A::RecombineMatrix)
    # Compare output of nzrows with similar functions for sparse arrays.
    # See ?nzrange.
    @testset "nzrows" begin
        S = sparse(A)
        rows = rowvals(S)
        vals = nonzeros(S)
        for j in axes(A, 2)
            ra = nzrows(A, j)
            rs = nzrange(S, j)
            @test length(ra) == length(rs)  # same number of non-zero rows
            for (ia, is) in zip(ra, rs)
                @test ia == rows[is]        # same non-zero rows
                @test A[ia, j] == vals[is]  # same values
            end
        end
    end
    nothing
end

function test_recombine_matrix(A::RecombineMatrix)
    bcleft, bcright = constraints(A)
    @assert bcleft === bcright  # only supported case for now
    @testset "Recombination matrix: $(bcleft)" begin
        test_nzrows(A)
        N, M = size(A)
        nl, nr = num_recombined(A)
        cl, cr = num_constraints(A)
        @assert nl == nr && cl == cr  # only supported case for now
        n = nl
        c = cl
        @test M == N - sum(num_constraints(A))

        # The sum of all rows is 2 in the recombined regions (because 2
        # recombined B-splines), and 1 in the centre.
        # Note that this is an arbitrary choice of the implementation.
        @test vec(sum(A, dims=1)) ≈ vcat(2 * ones(n), ones(M - 2n), 2 * ones(n))

        # Upper left corner, lower right corner, and centre.
        @test @views A[1:(n + c), 1:n] == A.ul
        @test @views A[(N - n - c + 1):N, (M - n + 1):M] == A.lr
        @test @views A[(n + c + 1):(N - n - c), (n + 1):(M - n)] == I

        u = rand(M)
        v1 = similar(u, N)
        v2 = similar(v1)

        # Test custom matrix-vector multiplication.
        # The second version uses the getindex() function to build a regular
        # sparse matrix, which is less efficient.
        mul!(v1, A, u)
        mul!(v2, sparse(A), u)
        @test v1 ≈ v2

        # Test left division.
        let v = copy(v1)
            u2 = A \ v
            @test u2 ≈ u

            # Modify coefficients of original basis near the border, so that the
            # resulting function has no representation in the recombined basis.
            rand!(v)
            @test_throws NoUniqueSolutionError A \ v
        end

        # Test 5-argument mul!.
        mul!(v2, A, u, 1, 0)  # equivalent to 3-argument mul!
        @test v1 ≈ v2

        randn!(u)  # use different u, otherwise A * u = v1
        mul!(v2, A, u, 2, -3)
        @test v2 ≈ 2 * (A * u) - 3 * v1
    end

    nothing
end

function test_boundary_conditions(R::RecombinedBSplineBasis)
    bl, br = constraints(R)
    @assert bl === br "different BCs on each boundary not yet supported"
    test_boundary_conditions(bl, R)
end

function test_boundary_conditions(ops::Tuple{Vararg{Derivative}},
                                  R::RecombinedBSplineBasis)
    # TODO adapt this for different BCs on the left and right boundaries
    @assert (ops, ops) === constraints(R)
    @test length(ops) == num_constraints(R)[1] == num_constraints(R)[2]
    @test length(R) == length(parent(R)) - 2 * length(ops)

    a, b = boundaries(R)
    N = length(R)
    k = order(R)

    # Verify that all basis functions satisfy the boundary conditions
    # for the derivative `D`, while they leave at least one degree of freedom
    # to set the other derivatives.

    # For instance, if D = 1 (Neumann BC), we want:
    #   (1) ϕⱼ'(a) = 0 ∀j,
    #   (2) ∑_j |ϕⱼ(a)| > 0,
    #   (3) ∑_j |ϕⱼ''(a)| > 0 (if k >= 3; otherwise the B-splines are not C²),
    # etc...

    for n = 0:(k - 1)
        # We cound use `dot` here, but it doesn't make a difference for
        # single derivatives. It does make a difference for linear combinations
        # of derivatives (e.g. Robin BCs).
        bsum = sum(1:N) do i
            f = R[i]
            fa = f(a, Derivative(n))
            fb = f(b, Derivative(n))
            abs(fa) + abs(fb)  # each of these must be 0 if n = D
        end
        # The sum must be zero if and only if n = D.
        # We consider that the result is zero, if it is negligible wrt the
        # derivative at the border of the first B-spline of the original
        # basis.
        B = parent(R)
        ε = 2 * num_recombined(R)[1] * eps(B[1](a, Derivative(n)))
        if Derivative(n) ∈ ops
            @test bsum ≤ ε
        else
            @test bsum > 1
        end
    end

    nothing
end

# More general BCs (Robin-like BCs and more combinations).
function test_boundary_conditions(ops, R::RecombinedBSplineBasis)
    @assert (ops, ops) === constraints(R)
    @test length(ops) == num_constraints(R)[1] == num_constraints(R)[2]
    @test length(R) == length(parent(R)) - 2 * length(ops)

    a, b = boundaries(R)
    N = length(R)

    for op in ops
        op_a, op_b = dot.(op, (LeftNormal(), RightNormal()))
        bsum = sum(1:N) do i
            f = R[i]
            fa = f(a, op_a)
            fb = f(b, op_b)
            abs(fa) + abs(fb)
        end
        B = parent(R)
        ε = 2 * num_recombined(R)[1] * eps(B[1](a, op_a))
        @test bsum ≤ ε
    end

    nothing
end

# Test relation between `support` and `nonzero_in_segment` (they're inverse of
# each other, in some sense...).
function test_support(R)
    for i in eachindex(R)
        sup = support(R, i)
        for n in sup[1:end-1]
            i ∈ nonzero_in_segment(R, n) || return false
        end
    end
    true
end

# Verify that basis recombination gives the right boundary conditions.
function test_basis_recombination()
    knots_base = gauss_lobatto_points(40)
    k = 4
    B = BSplineBasis(k, knots_base)
    @test constraints(B) === ((), ())
    @testset "Mixed derivatives" begin
        @inferred RecombineMatrix(Derivative.((0, 1)), B)
        @inferred RecombineMatrix(Derivative.((0, 2)), B)

        # Derivative orders must be in increasing order.
        @test_throws ArgumentError RecombineMatrix(Derivative.((1, 0)), B)

        # Unsupported combinations of derivatives.
        @test_throws ArgumentError RecombineMatrix(Derivative.((1, 2)), B)

        M = RecombineMatrix(Derivative.((0, 1)), B)
        test_recombine_matrix(M)
    end

    @test_throws ArgumentError RecombinedBSplineBasis(Derivative(k), B)

    @testset "Spline" begin
        R = RecombinedBSplineBasis(Derivative(1), B)
        coefs = rand(length(R))

        @testset "Equality" begin
            @test R != B
            @test R != RecombinedBSplineBasis(Derivative(0), B)
            @test R == RecombinedBSplineBasis(Derivative(1), B)
        end

        # This constructs a Spline using the recombined basis.
        S = @inferred Spline(R, coefs)
        @test S isa Spline

        @test coefficients(S) === coefs
        @test length(S) == length(R)
        @test basis(S) === R

        @test @allocated(Recombinations.parent_coefficients(R, coefs)) == 0

        # Construct a Spline in the original B-spline basis.
        @test @allocated(Splines.parent_spline(S)) == 0
        Sp = Splines.parent_spline(S)
        @test coefficients(Sp) == recombination_matrix(R) * coefs
        @test length(Sp) == length(B)
        @test basis(Sp) === B

        # Check that both splines are exactly the same
        @test Sp(0.4242) == S(0.4242)
        @test diff(Sp) == diff(S)  # this always returns a Spline in the original basis
        @test integral(Sp) == integral(S)  # same as above
    end

    @testset "Order $D" for D = 0:(k - 1)
        R = @inferred RecombinedBSplineBasis(Derivative(D), B)
        ops = (Derivative(D), )
        @test constraints(R) === (ops, ops)

        δl, δr = @inferred num_constraints(R)
        @test δl == δr == 1  # only one BC per boundary

        # In segment [t[k + 1], t[k + 2]], non-zero recombined functions are 1:k.
        @test @inferred(nonzero_in_segment(R, k + 1)) == 1:k

        # Near the boundaries, there are fewer recombined functions.
        @test nonzero_in_segment(R, k) == 1:(k - 1)
        @test nonzero_in_segment(R, k - 1) == 1:(k - 2)
        @test nonzero_in_segment(R, 2) == 1:1
        @test isempty(nonzero_in_segment(R, 1))

        N = length(parent(R))
        Nr = length(R)
        @test isempty(nonzero_in_segment(R, N + k - δr))
        @test nonzero_in_segment(R, N + k - 1 - δr) == Nr:Nr
        @test test_support(R)

        test_recombine_matrix(recombination_matrix(R))
        test_boundary_conditions(R)

        let
            N = length(R)
            a, b = boundaries(R)
            cl, cr = constraints(R)
            for f ∈ (repr, summary)
                @test startswith(
                    f(R), "$N-element RecombinedBSplineBasis of order $k",
                )
            end
        end

        @testset "Robin-like BCs" begin
            op = 1.1 * Derivative(D) + 4.2 * Derivative(D + 1)
            @test BSplineKit.max_order(op) == D + 1
            if D + 1 == k
                # "cannot resolve operators (...) with B-splines of order 4"
                @test_throws ArgumentError RecombinedBSplineBasis(op, B)
            else
                let Rs = RecombinedBSplineBasis(op, B)
                    @test constraints(Rs) === ((op, ), (op, ))
                    test_recombine_matrix(recombination_matrix(Rs))
                    test_boundary_conditions(Rs)
                end
                if D != 0
                    # Combine with Dirichlet BCs
                    ops = (Derivative(0), op)
                    Rs = RecombinedBSplineBasis(ops, B)
                    @test constraints(Rs) === (ops, ops)
                    test_recombine_matrix(recombination_matrix(Rs))
                    test_boundary_conditions(Rs)
                end
            end
        end

        # Simultaneously satisfies BCs of orders 0 to D.
        @testset "Mixed BCs" begin
            ops = ntuple(d -> Derivative(d - 1), D + 1)
            Rs = RecombinedBSplineBasis(ops, B)
            @test constraints(Rs) === (ops, ops)
            test_recombine_matrix(recombination_matrix(Rs))
            test_boundary_conditions(Rs)
        end
    end
    nothing
end

function test_natural_recombination()
    knots_base = gauss_lobatto_points(40)
    B = BSplineBasis(BSplineOrder(6), knots_base)
    R = RecombinedBSplineBasis(Natural(), B)
    M = recombination_matrix(R)
    test_nzrows(M)
    nothing
end

@testset "Basis recombination" begin
    test_basis_recombination()
    @testset "Natural BCs" begin
        test_natural_recombination()
    end
end

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines


text:
header("PARSE_LABELS")

header("UNCLAIMED LABEL")
@testset "unclaimed" begin
    # The funk are key, value pairs expected to be in the PnmlDict.
    for (node,funk) in [
        xml"""<declarations> </declarations>""" => [
            :content => ""
            ],
        xml"""<declarations atag="test1"> </declarations>""" => [
            :atag => "test1",
            :content => ""
            ],

            xml"""<declarations atag="test2">
                <something> some content </something>
                <something> other stuff </something>
                <something2 tag2="two"> <value/> <value tag3="three"/> </something2>
              </declarations>""" => [
            :atag => "test2",
            :something  => [Dict(:content => "some content"), Dict(:content => "other stuff")],
            :something2 => Dict(:value => [Dict(), Dict(:tag3 => "three")], :tag2 => "two")
            ],

        xml"""<foo><declarations> </declarations></foo>""" => [
            :declarations => Dict(:content => "")
            ],

        # no content, no attribute results in empty dict.
        xml"""<null></null>""" => [],
        xml"""<null2/>""" => [],
        # no content, with attribute
        xml"""<null at="null"></null>""" => [
            :at => "null"
            ],
        xml"""<null2 at="null2" />""" => [
            :at => "null2"
            ],

        # empty content, no attribute
        xml"""<empty> </empty>""" => [
            :content => ""
            ],
        # empty content, with attribute
        xml"""<empty at="empty"> </empty>""" => [
            :content => "",
            :at => "empty"
            ],

        xml"""<foo id="testid"/>""" => [
            :id => :testid
            ],
        ]
        reg1=PNML.IDRegistry()
        reg2=PNML.IDRegistry()

        u = PNML.unclaimed_label(node, reg=reg1)
        l = PNML.PnmlLabel(u, node)
        a = PNML.anyelement(node, reg=reg2)

        @test_call PNML.unclaimed_label(node, reg=reg1)
        @test_call PNML.PnmlLabel(u, node)
        @test_call PNML.anyelement(node, reg=reg2)

        @test !isnothing(u)
        @test !isnothing(l)
        @test !isnothing(a)

        @test u isa Pair{Symbol, Dict{Symbol, Any}}
        @test l isa PNML.PnmlLabel
        @test a isa PNML.AnyElement

        nn = Symbol(nodename(node))
        @test u.first === nn
        @test tag(l) === nn
        @test tag(a) === nn

        @test u.second isa PnmlDict
        @test l.dict isa PnmlDict
        @test a.dict isa PnmlDict

        for (key,val) in funk
            @test u.second[key] == val
            @test l.dict[key] == val
            @test a.dict[key] == val
        end

        haskey(u.second, :id) && @test PNML.isregistered(reg1, u.second[:id])
        haskey(l.dict, :id) && @test PNML.isregistered(reg1, l.dict[:id])
        haskey(a.dict, :id) && @test PNML.isregistered(reg2, a.dict[:id])

        @test_call  PNML.isregistered(reg2, :id)

        @show u
        @show l
        @show a
        println()
    end
end

#header("GET_LABEL")
@testset "get rate label" begin
    n = parse_node(xml"""<transition id ="birth">
        <rate> <text>0.3</text> </rate>
    </transition>""", reg=PNML.IDRegistry())
    printnode(n)
    @show l = PNML.labels(n)
    @test PNML.tag(first(l)) === :rate # only label
    @test PNML.get_label(n, :rate) === first(PNML.labels(n))
    @test PNML.rate(n) ≈ 0.3

    @test_call PNML.has_labels(n)
    @test_call PNML.labels(n)
    @test_call PNML.has_label(n, :rate)
    @test_call PNML.get_label(n, :rate)
    @test_call PNML.rate(n)
end

header("DECLARATION")
@testset "empty declarations" begin
    # The attribute should be ignored.
    n = parse_node(xml"""
        <declaration key="test">
          <structure>
           <declarations>
           </declarations>
          </structure>
        </declaration>
        """; reg = PNML.IDRegistry())

    @show typeof(n), fieldnames(typeof(n))
    printnode(n)

    @test typeof(n) <: PNML.Declaration
    @test xmlnode(n) isa Maybe{EzXML.Node}
    @test typeof(PNML.declarations(n)) <: Vector{PNML.AbstractDeclaration}
    @test length(PNML.declarations(n)) == 0

    @test typeof(n.com) <: PNML.ObjectCommon
    @test PNML.graphics(n) === nothing
    @test PNML.tools(n) === nothing
    @test PNML.labels(n) === nothing
    @test isempty(n.com)

    @test_call PNML.declarations(n)
    @test_call PNML.graphics(n)
    @test_call PNML.tools(n)
    @test_call PNML.labels(n)
end

@testset "declaration tree" begin
        node = xml"""
        <declaration>
            <structure>
                <declarations>
                    <namedsort id="LegalResident" name="LegalResident">
                        <cyclicenumeration>
                            <feconstant id="LegalResident0" name="0"/>
                            <feconstant id="LegalResident1" name="1"/>
                        </cyclicenumeration>
                    </namedsort>
                    <namedsort id="MICSystem" name="MICSystem">
                        <cyclicenumeration>
                            <feconstant id="MICSystem0" name="0"/>
                            <feconstant id="MICSystem1" name="1"/>
                        </cyclicenumeration>
                    </namedsort>
                    <namedsort id="CINFORMI" name="CINFORMI">
                        <cyclicenumeration>
                            <feconstant id="CINFORMI0" name="0"/>
                            <feconstant id="CINFORMI1" name="1"/>
                        </cyclicenumeration>
                    </namedsort>
            </declarations>
        </structure>
    </declaration>
    """
    reg = PNML.IDRegistry()
    n = parse_node(node; reg)
    printnode(n)

    @test typeof(n) <: PNML.Declaration
    @test xmlnode(n) isa Maybe{EzXML.Node}
    @test length(PNML.declarations(n)) == 3
    @test_call PNML.declarations(n)

    for d in PNML.declarations(n)
        @test typeof(d) <: PNML.AbstractDeclaration
        @test typeof(d) <: PNML.SortDeclaration
        @test typeof(d) <: PNML.NamedSort

        #@show d
        #@show fieldtypes(typeof(d))
        #@show fieldnames(typeof(d))
        #@show fieldtypes(typeof(d.def))
        #@show fieldnames(typeof(d.def))

        @test PNML.isregistered(reg, pid(d))
        @test_call PNML.isregistered(reg, pid(d))
        @test Symbol(PNML.name(d)) === pid(d) # name and id are the same.
        @test d.def isa PNML.AnyElement #TODO implement definitions?
        @test tag(d.def) === :cyclicenumeration
        @test haskey(d.def.dict, :feconstant)

        for x in d.def.dict[:feconstant]
            @test x isa PnmlDict
            @test PNML.isregistered(reg, pid(x))
            @test x[:name] isa String
            @test endswith(string(pid(x)), x[:name])
        end
        #println()
    end
end

header("PT initMarking")
@testset "PT initMarking" begin
    node = xml"""
 <initialMarking>
    <text>1.0</text>
    <toolspecific tool="org.pnml.tool" version="1.0">
        <tokengraphics>
            <tokenposition x="6" y="9"/>
        </tokengraphics>
    </toolspecific>
 </initialMarking>
 """

    n = parse_node(node; reg=PNML.IDRegistry())
    printnode(n)
    @test typeof(n) <: PNML.PTMarking
    @test xmlnode(n) isa Maybe{EzXML.Node}
    @test typeof(n.value) <: Number
    @test n.value == n()

    mark1 = PNML.PTMarking(2)
    @test_call PNML.PTMarking(2)
    @test typeof(mark1()) == typeof(2)
    @test mark1() == 2
    @test_call mark1()

    mark2 = PNML.PTMarking(3.5)
    @test_call PNML.PTMarking(3.5)
    @test typeof(mark2()) == typeof(3.5)
    @test mark2() ≈ 3.5
    @test_call mark2()

    mark3 = PNML.PTMarking()
    @test_call PNML.PTMarking()
    @test typeof(mark3()) == typeof(PNML.default_marking(PnmlCore())())
    @test mark3() == PNML.default_marking(PnmlCore())()
    @test_call mark3()
end

header("HL Marking")
@testset "HL initMarking" begin
    str = """
 <hlinitialMarking>
     <text><All,All></text>
     <structure>
            <tuple>
              <subterm>
                <all>
                  <usersort declaration="N1"/>
                </all>
              </subterm>
              <subterm>
                <all>
                  <usersort declaration="N2"/>
                </all>
              </subterm>
            </tuple>
     </structure>
 </hlinitialMarking>
    """
    n = parse_node(to_node(str); reg = PNML.IDRegistry())
    @show typeof(n), fieldnames(typeof(n))
    printnode(n)

    @test typeof(n) <: PNML.AbstractLabel
    @test typeof(n) <: PNML.HLMarking
    @test xmlnode(n) isa Maybe{EzXML.Node}
    @test n.text == "<All,All>"
    @test n.term !== nothing
    @test n.term isa PNML.AbstractTerm

    #@show typeof(n.term), fieldnames(typeof(n.term))
    #@show n.term
    @test tag(n.term) === :tuple
    @test n.term.dict[:subterm][1][:all] !== nothing
    @test n.term.dict[:subterm][1][:all][:usersort][:declaration] == "N1"
    @test n.term.dict[:subterm][2][:all][:usersort][:declaration] == "N2"
end

@testset "text" begin
    str1 = """
 <text>ready</text>
    """
    n = parse_node(to_node(str1); reg = PNML.IDRegistry())
    @test n == "ready"

    str2 = """
 <text>
ready
</text>
    """
    n = parse_node(to_node(str2); reg = PNML.IDRegistry())
    @test n == "ready"

    str3 = """
 <text>    ready  </text>
    """
    n = parse_node(to_node(str3); reg = PNML.IDRegistry())
    @test n == "ready"

    str4 = """
     <text>ready
to
go</text>
    """
    n = parse_node(to_node(str4); reg = PNML.IDRegistry())
    @test n == "ready\nto\ngo"
end

header("STRUCTURE")
@testset "structure" begin
    node = xml"""
     <structure>
         <tuple>
              <subterm>
                <all>
                  <usersort declaration="N1"/>
                </all>
              </subterm>
              <subterm>
                <all>
                  <usersort declaration="N2"/>
                </all>
              </subterm>
         </tuple>
     </structure>
    """

    n = parse_node(node; reg = PNML.IDRegistry())
    printnode(n)
    @test n isa PNML.Structure
    @test xmlnode(n) isa Maybe{EzXML.Node}
    @test tag(n) === :structure
    @test n.dict isa PnmlDict
    @test tag(n) === :structure
    @test n.dict[:tuple][:subterm][1][:all][:usersort][:declaration] == "N1"
    @test n.dict[:tuple][:subterm][2][:all][:usersort][:declaration] == "N2"
end

header("SORT TYPE")
@testset "type" begin
    n1 = xml"""
 <type>
     <text>N2</text>
     <structure> <usersort declaration="N2"/> </structure>
 </type>
    """
    @testset for node in [n1]
        n = parse_node(node; reg = PNML.IDRegistry())
        printnode(n)
        @test typeof(n) <: PNML.AnyElement
        @test tag(n) === :type
        @test n.dict[:text][:content] == "N2"
        @test n.dict[:structure][:usersort][:declaration] == "N2"
    end
end

header("CONDITION")
@testset "condition" begin
    n1 = xml"""
 <condition>
     <text>(x==1 and y==1 and d==1)</text>
     <structure> <or> #TODO </or> </structure>
 </condition>
    """
    @testset for node in [n1]
        n = parse_node(node; reg = PNML.IDRegistry())
        printnode(n)
        @test typeof(n) <: PNML.Condition
        @test n.text !== nothing
        @test n.term !== nothing
        @test tag(n.term) === :or
        @test n.com.graphics === nothing
        @test n.com.tools === nothing || isempty(n.com.tools)
        @test n.com.labels === nothing || isempty(n.com.labels)
    end
end

@testset "inscription" begin
    n1 = xml"<inscription> <text> 12 </text> </inscription>"
    @testset for node in [n1]
        n = parse_node(node; reg = PNML.IDRegistry())
        printnode(n)
        @test typeof(n) <: PNML.PTInscription
        @test xmlnode(n) isa Maybe{EzXML.Node}
        @test n.value == 12
        @test n.com.graphics === nothing
        @test n.com.tools === nothing || isempty(n.com.tools)
        @test n.com.labels === nothing || isempty(n.com.labels)
    end
end

@testset "hlinscription" begin
    n1 = xml"""
 <hlinscription>
     <text><x,v></text>
     <structure>
            <tuple>
              <subterm>
                <variable refvariable="x"/>
              </subterm>
              <subterm>
                <variable refvariable="v"/>
              </subterm>
            </tuple>
     </structure>
 </hlinscription>
 """
    @testset for node in [n1]
        n = parse_node(node; reg = PNML.IDRegistry())
        printnode(n)
        @test typeof(n) <: PNML.HLInscription
        @test xmlnode(n) isa Maybe{EzXML.Node}
        @test tag(n.term) === :tuple
        @test n.term.dict[:subterm][1][:variable][:refvariable] == "x"
        @test n.term.dict[:subterm][2][:variable][:refvariable] == "v"
        @test n.text == "<x,v>"
    end
end


text:
using ProjectiveDictionaryPairLearning
using MAT, Base.Test

###### Testing TrainDPL
# load DictMat, EncoderMat, TrData, TrLabel, DictSize, tau, lambda, gamma
# prepared using
    # matwrite("testDataForTrainDpl.mat", @compat Dict(
    #     "DictMat" => DictMat,
    #     "EncoderMat" => EncoderMat,
    #     "TrData" => TrData,
    #     "TrLabel" => TrLabel,
    #     "DictSize" => DictSize,
    #     "tau" => tau,
    #     "lambda" => lambda,
    #     "gamma" => gamma,
    # ))
matFileName = "testDataForTrainDpl.mat"
println("matFileName: $matFileName.")

tic()
data = matread(joinpath(dirname(@__FILE__), matFileName))
DictMat, EncoderMat, TrData, TrLabel = data["DictMat"], data["EncoderMat"], data["TrData"], data["TrLabel"]
DictSize, tau, lambda, gamma = data["DictSize"], data["tau"], data["lambda"], data["gamma"]
println("Data loaded in $(toq()) s")

tic()
DEMO = true
TestDictMat, TestEncoderMat = TrainDPL(TrData, TrLabel, DictSize, tau, lambda, gamma, DEMO)
println("TrainDPL passed in $(toq()) s")

# for i = 1:length(DictSize)
#     @test_approx_eq(DictMat[i], TestDictMat[i])
# end
# @test_approx_eq(EncoderMat, TestEncoderMat)

testMat = [i+j for i in 1.0:5, j in 1.0:5]
resultMat = [
 0.21081851067789195 0.25819888974716115 0.29019050004400465 0.31311214554257477 0.3302891295379082
 0.31622776601683794 0.34426518632954817 0.36273812505500586 0.3757345746510897 0.3853373177942262
 0.4216370213557839 0.43033148291193524 0.435285750066007 0.4383570037596047 0.4403855060505443
 0.5270462766947299 0.5163977794943223 0.5078333750770082 0.5009794328681196 0.4954336943068623
 0.6324555320336759 0.6024640760767093 0.5803810000880093 0.5636018619766345 0.5504818825631803
]
@test_approx_eq(resultMat, normcol_lessequal(testMat))

accuracy = dpldemo();
@test accuracy > 0.95

println("Finished unit tests!")


text:
using DFTShims
using AxisArrays
using Unitful

macro lintpragma(s) end
@lintpragma("Ignore use of undeclared D")
@lintpragma("Ignore use of undeclared U")

const DD = Dispatch.Dimensions
const DH = Dispatch.Hartree

spin = AxisArray(zeros((10, 2)), Axis{:radius}(0u"m":1u"m":9u"m"), Axis{:spin}((:+, :-)))
nospin = AxisArray(zeros((10, 2)), Axis{:radius}(0u"m":1u"m":9u"m"), Axis{:n}((1, 2)))
const polarized = ColinearSpin()
const unpolarized = SpinDegenerate()

@testset "has_axis" begin
    @test has_axis(spin, :spin)
    @test has_axis(spin, Axis{:spin})
    @test has_axis(typeof(spin), :spin)
    @test has_axis(typeof(spin), Axis{:spin})
    @test !has_axis(nospin, :spin)
    @test !has_axis(nospin, Axis{:spin})
    @test !has_axis(typeof(nospin), :spin)
    @test !has_axis(typeof(nospin), Axis{:spin})
end

@testset "spin polarization" begin
    @test @inferred is_spin_polarized(spin)
    @test ! @inferred is_spin_polarized(nospin)

    @test @inferred is_spin_polarized(typeof(spin))
    @test ! @inferred is_spin_polarized(typeof(nospin))

    @test @inferred(SpinCategory(spin)) === polarized
    @test @inferred(SpinCategory(nospin)) === unpolarized
end

@testset "components" begin
    @test @inferred(components(UnitfulHartree.ρ, unpolarized)) == (:ρ,)
    @test components(UnitfulHartree.ρ, polarized) == (:α, :β)
    @test components(UnitfulHartree.ϵ, polarized) == (:α, :β)
    @test length(components(DD.Scalars.∂ϵ_∂σ, unpolarized)) == 1
    @test length(@inferred(components(UnitfulHartree.∂ϵ_∂σ, polarized))) == 3
    @test length(@inferred(components(UnitfulHartree.∂²ϵ_∂ρ², polarized))) == 3
    @test length(@inferred(components(UnitfulHartree.∂²ϵ_∂ρ∂σ, unpolarized))) == 1
    @test length(@inferred(components(UnitfulHartree.∂²ϵ_∂ρ∂σ, polarized))) == 6
    @test length(@inferred(components(UnitfulHartree.∂²ϵ_∂σ², polarized))) == 6
    @test length(@inferred(components(UnitfulHartree.∂³ϵ_∂ρ³, polarized))) == 4
    @test length(@inferred(components(UnitfulHartree.∂³ϵ_∂ρ²∂σ, polarized))) == 9
    @test length(@inferred(components(UnitfulHartree.∂³ϵ_∂ρ∂σ², polarized))) == 12
    @test length(@inferred(components(UnitfulHartree.∂³ϵ_∂σ³, polarized))) == 10
end

@testset "functional category" begin
    @test @inferred(FunctionalCategory(DH.Scalars.ρ{Int64})) === LDA()
    @test @inferred(FunctionalCategory(DH.Scalars.∂ϵ_∂σ)) === GGA()
    @test @inferred(FunctionalCategory(DD.Scalars.∂²ϵ_∂ρ²)) === LDA()
    @test @inferred(FunctionalCategory(DD.Scalars.∂²ϵ_∂ρ∂σ)) === GGA()
    @test @inferred(FunctionalCategory(DD.Scalars.∂²ϵ_∂σ²)) === GGA()
    @test @inferred(FunctionalCategory(UnitfulHartree.∂³ϵ_∂ρ³)) === LDA()
    @test @inferred(FunctionalCategory(UnitfulHartree.∂³ϵ_∂ρ²∂σ)) === GGA()
    @test @inferred(FunctionalCategory(UnitfulHartree.∂³ϵ_∂ρ∂σ²)) === GGA()
    @test @inferred(FunctionalCategory(UnitfulHartree.∂³ϵ_∂σ³)) === GGA()
end

@testset "Type concretization" begin
    @test @inferred(concretize_type(DD.Scalars.ρ, Int16)) === DH.Scalars.ρ{Int16}
    @test @inferred(concretize_type(DH.Scalars.ρ, Int16)) === DH.Scalars.ρ{Int16}
    @test @inferred(concretize_type(typeof(1.0u"m"), Int16)) === typeof(1.0u"m")
    @test @inferred(concretize_type(DH.Scalars.ρ, Quantity{Int16, D, U} where {D, U})) ===
            DH.Scalars.ρ{Int16}
    @test @inferred(concretize_type(DD.Scalars.ρ, Quantity{Int16, D, U} where {D, U})) ===
            DH.Scalars.ρ{Int16}
    const Q = Quantity{Int16, typeof(dimension(1u"C")), U} where U
    @test @inferred(concretize_type(DD.Scalars.ρ, Q)) === DH.Scalars.ρ{Int16}
end

max_stars_repo_name: mdavezac/DFTShims

max_issues_repo_name: mdavezac/DFTShims

max_forks_repo_name: mdavezac/DFTShims


text:
# Autogenerated wrapper script for OpenCV_jll for x86_64-apple-darwin
export libopencv_calib3d, libopencv_core, libopencv_dnn, libopencv_features2d, libopencv_flann, libopencv_gapi, libopencv_highgui, libopencv_imgcodecs, libopencv_imgproc, libopencv_julia, libopencv_objdetect, libopencv_stitching, libopencv_video, libopencv_videoio

using Qt5Base_jll
using Libglvnd_jll
using libcxxwrap_julia_jll
JLLWrappers.@generate_wrapper_header("OpenCV")
JLLWrappers.@declare_library_product(libopencv_calib3d, "@rpath/libopencv_calib3d.4.5.dylib")
JLLWrappers.@declare_library_product(libopencv_core, "@rpath/libopencv_core.4.5.dylib")
JLLWrappers.@declare_library_product(libopencv_dnn, "@rpath/libopencv_dnn.4.5.dylib")
JLLWrappers.@declare_library_product(libopencv_features2d, "@rpath/libopencv_features2d.4.5.dylib")
JLLWrappers.@declare_library_product(libopencv_flann, "@rpath/libopencv_flann.4.5.dylib")
JLLWrappers.@declare_library_product(libopencv_gapi, "@rpath/libopencv_gapi.4.5.dylib")
JLLWrappers.@declare_library_product(libopencv_highgui, "@rpath/libopencv_highgui.4.5.dylib")
JLLWrappers.@declare_library_product(libopencv_imgcodecs, "@rpath/libopencv_imgcodecs.4.5.dylib")
JLLWrappers.@declare_library_product(libopencv_imgproc, "@rpath/libopencv_imgproc.4.5.dylib")
JLLWrappers.@declare_library_product(libopencv_julia, "@rpath/libopencv_julia.dylib")
JLLWrappers.@declare_library_product(libopencv_objdetect, "@rpath/libopencv_objdetect.4.5.dylib")
JLLWrappers.@declare_library_product(libopencv_stitching, "@rpath/libopencv_stitching.4.5.dylib")
JLLWrappers.@declare_library_product(libopencv_video, "@rpath/libopencv_video.4.5.dylib")
JLLWrappers.@declare_library_product(libopencv_videoio, "@rpath/libopencv_videoio.4.5.dylib")
function __init__()
    JLLWrappers.@generate_init_header(Qt5Base_jll, Libglvnd_jll, libcxxwrap_julia_jll)
    JLLWrappers.@init_library_product(
        libopencv_calib3d,
        "lib/libopencv_calib3d.4.5.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_core,
        "lib/libopencv_core.4.5.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_dnn,
        "lib/libopencv_dnn.4.5.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_features2d,
        "lib/libopencv_features2d.4.5.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_flann,
        "lib/libopencv_flann.4.5.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_gapi,
        "lib/libopencv_gapi.4.5.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_highgui,
        "lib/libopencv_highgui.4.5.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_imgcodecs,
        "lib/libopencv_imgcodecs.4.5.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_imgproc,
        "lib/libopencv_imgproc.4.5.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_julia,
        "lib/libopencv_julia.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_objdetect,
        "lib/libopencv_objdetect.4.5.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_stitching,
        "lib/libopencv_stitching.4.5.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_video,
        "lib/libopencv_video.4.5.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_videoio,
        "lib/libopencv_videoio.4.5.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()