text:
"""
    convert_units!(p::Particles, h::AbstractGadgetHeader, units::Symbol=:full)

Converts the units of the particles properties (see [`Particles`](@ref)) from a snapshot or
subfind header in-place.
Use `:full` for units as `Unitful` quantities, `:physical` for values converted
to physical units (kpc, km/s, solar metallicities etc.), or `:sim` for values in simulation units.
The individual metallicities in "Zs" are converted to a one-dimensional `Vector` of values in solar
metallicities in any case.
"""
function convert_units!(p::Particles, h::AbstractGadgetHeader, units::Symbol=:full)
    if units === :physical
        for (key, vals) in pairs(p)
            convert_units_physical!(vals, key, h)
        end
    elseif units === :full
        for (key, vals) in pairs(p)
            p[key] = convert_units_full(vals, key, h)
        end
    end

    props = keys(p)
    if :zs in props && length(p.zs) > 0 && (:im in props || :mass in props)
        p.zs = convert_units_solar_metallicity(p.zs, :im in props ? p.iM : p.mass)
    end

    return p
end

function convert_to_unitful!(
    a::AbstractArray{T},
    prop::AbstractString,
    h::AbstractGadgetHeader,
    simunits::Dict,
) where {T<:Real}
    haskey(simunits, prop) || return a

    factor_func, u = simunits[prop]

    if !isnothing(factor_func)
        a .*= factor_func(h)
    end

    if !isnothing(u)
        Q = Quantity{T,dimension(u),typeof(u)}
        return reinterpret(Q, a)
    end

    return a
end

function convert_to_physical!(
    a::AbstractArray{<:Real},
    prop::AbstractString,
    h::AbstractGadgetHeader,
    simunits::Dict,
)
    haskey(simunits, prop) || return a

    factor_func, _ = simunits[prop]

    if !isnothing(factor_func)
        a .*= factor_func(h)
    end

    return a
end

function convert_to_sim!(
    a::AbstractArray{<:Real},
    prop::AbstractString,
    h::AbstractGadgetHeader,
    simunits::Dict,
)
    haskey(simunits, prop) || return a

    factor_func, _ = simunits[prop]

    if !isnothing(factor_func)
        a ./= factor_func(h)
    end

    return a
end

function convert_to_sim!(
    a::AbstractArray{Quantity},
    prop::AbstractString,
    h::AbstractGadgetHeader,
    simunits::Dict,
)
    astrip = ustrip(a)
    haskey(simunits, prop) || return astrip

    factor_func, _ = simunits[prop]

    if !isnothing(factor_func)
        astrip ./= factor_func(h)
    end

    return astrip
end

function convert_to_unitful(
    a::AbstractArray{T},
    prop::AbstractString,
    h::AbstractGadgetHeader,
    simunits::Dict,
) where {T<:Real}
    haskey(simunits, prop) || return copy(a)

    factor_func, u = simunits[prop]

    if isnothing(factor_func)
        if isnothing(u)
            return copy(a)
        else
            return a .* u
        end
    else
        if isnothing(u)
            return a .* factor_func(h)
        else
            return a .* factor_func(h) .* u
        end
    end
end

function convert_to_physical(
    a::AbstractArray{<:Real},
    prop::AbstractString,
    h::AbstractGadgetHeader,
    simunits::Dict,
)
    haskey(simunits, prop) || return copy(a)

    factor_func, _ = simunits[prop]

    if !isnothing(factor_func)
        return a .* factor_func(h)
    end

    return copy(a)
end

function convert_to_sim(
    a::AbstractArray{<:Real},
    prop::AbstractString,
    h::AbstractGadgetHeader,
    simunits::Dict,
)
    haskey(simunits, prop) || return copy(a)

    factor_func, _ = simunits[prop]

    if !isnothing(factor_func)
        return a ./ factor_func(h)
    end

    return copy(a)
end

function convert_to_sim(
    a::AbstractArray{Quantity},
    prop::AbstractString,
    h::AbstractGadgetHeader,
    simunits::Dict,
)
    astrip = ustrip(a)
    haskey(simunits, prop) || return copy(astrip)

    factor_func, _ = simunits[prop]

    if !isnothing(factor_func)
        return astrip ./ factor_func(h)
    end

    return copy(astrip)
end

function convert_to_unitful(a::Real, prop::AbstractString, h::AbstractGadgetHeader, simunits::Dict)
    haskey(simunits, prop) || return a

    factor_func, u = simunits[prop]

    if !isnothing(factor_func)
        a *= factor_func(h)
    end

    if isnothing(u)
        return a
    else
        return a * u
    end
end

function convert_to_physical(
    a::Real,
    prop::AbstractString,
    h::AbstractGadgetHeader,
    simunits::Dict,
)
    haskey(simunits, prop) || return a

    factor_func, _ = simunits[prop]

    if isnothing(factor_func)
        return a
    else
        return a * factor_func(h)
    end
end

function convert_to_sim(
    a::Real,
    prop::AbstractString,
    h::AbstractGadgetHeader,
    simunits::Dict,
)
    haskey(simunits, prop) || return a

    factor_func, _ = simunits[prop]

    if isnothing(factor_func)
        return a
    else
        return a / factor_func(h)
    end
end

function convert_to_sim(
    a::Quantity,
    prop::AbstractString,
    h::AbstractGadgetHeader,
    simunits::Dict,
)
    astrip = ustrip(a)
    haskey(simunits, prop) || return astrip

    factor_func, _ = simunits[prop]

    if isnothing(factor_func)
        return astrip
    else
        return astrip / factor_func(h)
    end
end

function convert_units_subfind_prop(
    val::Union{Real,AbstractArray{<:Real}},
    prop::AbstractString,
    h::AbstractGadgetHeader,
    units::Symbol=:full;
    verbose::Bool=false,
)
    # TODO: SPIN (units: angular momentum), DSUB (vel. dispersion units?),
    # SLUM, SLAT, SLOB, DUST, SZ (units?), SSFR (units really M⊙/yr?)
    if prop[1] === 'R' || prop in ["GPOS", "BGPO", "BGRA", "SPOS", "SCM", "SHMR"]
        return convert_units_pos(val, h, units)
    elseif prop[1] === 'V' || prop in ["SVEL"]
        return convert_units_vel(val, h, units)
    elseif (prop[1] === 'M' && prop != "MBID") || prop in ["BGMA", "SMST"]
        return convert_units_mass(val, h, units)
    elseif prop in ["SAGE"]
        return convert_units_age(val, h, units)
        # elseif prop in ["TGAS"] # TODO: unclear what KeV unit is
        # return convert_units_temp(val, h, units)
        # elseif prop in ["LGAS"] # TODO: energy unit conversion (10^44 erg/s)
        # return convert_units_temp(val, h, units)
    elseif prop in ["SSFR"] && units === :full
        return val * u"Msun" / Unitful.yr
    else
        verbose && @warn "The quantity $prop is returned in simulation units"
        return val
    end
end

@doc raw"""
    convert_units_age(val::Union{Real,AbstractArray{<:Real}}, h::AbstractGadgetHeader, units::Symbol=:full)
    convert_units_physical_age(val::Real, h)
    convert_units_physical_age(vals::AbstractArray{<:Real}, h)
    convert_units_physical_age!(vals::AbstractArray{<:Real}, h)
    convert_units_full_age(val::Real, h)
    convert_units_full_age(vals::AbstractArray{<:Real}, h)

Converts ages from the scale factor ``a`` via the lookback time according to the cosmology
defined by the snapshot header.
Full returns values in `Unitful` quantities, whereas physical returns the physical value in Gyr.
"""
function convert_units_age(
    val::Union{Real,AbstractArray{<:Real}},
    h::AbstractGadgetHeader,
    units::Symbol=:full,
)
    if units === :physical
        convert_units_physical_age(val, h)
    elseif units === :full
        convert_units_full_age(val, h)
    else
        val
    end
end
function convert_units_physical_age(val::Real, h)
    convert_units_full_age(val, h) |> ustrip
end
function convert_units_full_age(val::T, h) where {T<:Real}
    c = cosmology(h=h.h0, OmegaM=h.omega_0)
    lookback_time(c, 1 / val - 1) - lookback_time(c, h.z) # in Gyr
end
function convert_units_physical_age(vals::AbstractArray{<:Real}, h)
    convert_units_full_age(vals, h) .|> ustrip
end
function convert_units_physical_age!(vals::AbstractArray{<:Real}, h)
    vals .= convert_units_physical_age(vals, h)
end
function convert_units_full_age(vals::AbstractArray{T}, h) where {T<:Real}
    c = cosmology(h=h.h0, OmegaM=h.omega_0)
    lookback_time.((c,), 1 ./ vals .- 1) .- lookback_time(c, h.z) # in Gyr
end

"""
    convert_units_solar_metallicity(zs::AbstractMatrix{<:Number}, mass::AbstractVector{<:Number})

Converts an `11×n` mass matrix and an `n`-length mass vector to a metallicity vector of length `n`,
in solar metallicities.
"""
function convert_units_solar_metallicity(zs::AbstractMatrix{<:Number}, mass::AbstractVector{<:Number})
    dropdims(sum(zs[2:11, :]; dims=1); dims=1) ./ (mass .- dropdims(sum(zs; dims=1); dims=1)) ./ 0.0134 # in solar metallicities
end


text:
using CollisionDetection

abstract type LocalOperator <: Operator end

function allocatestorage(op::LocalOperator, test_functions, trial_functions,
    storage_trait::Type{Val{:bandedstorage}}, longdelays_trait)

    T = scalartype(op, test_functions, trial_functions)

    M = Int[]
    N = Int[]
    V = T[]

    function store(v,m,n)
        push!(M,m)
        push!(N,n)
        push!(V,v)
    end

    function freeze()
        nrows = numfunctions(test_functions)
        ncols = numfunctions(trial_functions)
        return sparse(M,N,V, nrows, ncols)
    end

    return freeze, store
end

function allocatestorage(op::LocalOperator, testfunctions, trialfunctions,
    storage_trait::Type{Val{:sparsedicts}}, longdelays_trait)

    T = scalartype(op, testfunctions, trialfunctions)

    m = numfunctions(testfunctions)
    n = numfunctions(trialfunctions)
    Z = SparseMatrixDict{T,Int}(m,n)

    store(v,m,n) = (Z[m,n] += v)
    freeze() = SparseArrays.SparseMatrixCSC(Z)

    return freeze, store
end

function allocatestorage(op::LocalOperator, test_functions, trial_functions,
    storage_trait::Type{Val{:densestorage}}, longdelays_trait)

    T = scalartype(op, test_functions, trial_functions)

    Z = zeros(T, numfunctions(test_functions), numfunctions(trial_functions))
    store(v,m,n) = (Z[m,n] += v)
    freeze() = Z

    return freeze, store
end

function assemble!(biop::LocalOperator, tfs::Space, bfs::Space, store,
        threading::Type{Threading{:multi}};
        quadstrat=defaultquadstrat(biop, tfs, bfs))

    if geometry(tfs) == geometry(bfs)
        return assemble_local_matched!(biop, tfs, bfs, store; quadstrat)
    end

    if CompScienceMeshes.refines(geometry(tfs), geometry(bfs))
        return assemble_local_refines!(biop, tfs, bfs, store; quadstrat)
    end

    return assemble_local_mixed!(biop, tfs, bfs, store; quadstrat)
end

function assemble!(biop::LocalOperator, tfs::Space, bfs::Space, store,
    threading::Type{Threading{:single}};
    quadstrat=defaultquadstrat(biop, tfs, bfs))

if geometry(tfs) == geometry(bfs)
    return assemble_local_matched!(biop, tfs, bfs, store; quadstrat)
end

if CompScienceMeshes.refines(geometry(tfs), geometry(bfs))
    return assemble_local_refines!(biop, tfs, bfs, store; quadstrat)
end

return assemble_local_mixed!(biop, tfs, bfs, store; quadstrat)
end

function assemble_local_matched!(biop::LocalOperator, tfs::Space, bfs::Space, store;
    quadstrat=defaultquadstrat(biop, tfs, bfs))

    tels, tad, ta2g = assemblydata(tfs)
    bels, bad, ba2g = assemblydata(bfs)

    bg2a = zeros(Int, length(geometry(bfs)))
    for (i,j) in enumerate(ba2g) bg2a[j] = i end

    trefs = refspace(tfs)
    brefs = refspace(bfs)

    qd = quaddata(biop, trefs, brefs, tels, bels, quadstrat)
    locmat = zeros(scalartype(biop, trefs, brefs), numfunctions(trefs), numfunctions(brefs))
    for (p,cell) in enumerate(tels)
        P = ta2g[p]
        q = bg2a[P]
        q == 0 && continue

        qr = quadrule(biop, trefs, brefs, cell, qd, quadstrat)
        fill!(locmat, 0)
        cellinteractions_matched!(locmat, biop, trefs, brefs, cell, qr)

        for i in 1 : size(locmat, 1), j in 1 : size(locmat, 2)
            for (m,a) in tad[p,i], (n,b) in bad[q,j]
                store(a * locmat[i,j] * b, m, n)
end end end end

function assemble_local_refines!(biop::LocalOperator, tfs::Space, bfs::Space, store;
    quadstrat=defaultquadstrat(biop, tfs, bfs))

    println("Using 'refines' algorithm for local assembly:")

    tgeo = geometry(tfs)
    bgeo = geometry(bfs)
    @assert CompScienceMeshes.refines(tgeo, bgeo)

    trefs = refspace(tfs)
    brefs = refspace(bfs)

    tels, tad, ta2g = assemblydata(tfs)
    bels, bad, ba2g = assemblydata(bfs)

    bg2a = zeros(Int, length(geometry(bfs)))
    for (i,j) in enumerate(ba2g) bg2a[j] = i end

    qd = quaddata(biop, trefs, brefs, tels, bels, quadstrat)

    print("dots out of 10: ")
    todo, done, pctg = length(tels), 0, 0
    for (p,tcell) in enumerate(tels)

        P = ta2g[p]
        Q = CompScienceMeshes.parent(tgeo, P)
        q = bg2a[Q]

        bcell = bels[q]
        @assert overlap(tcell, bcell)

        isct = intersection(tcell, bcell)
        for cell in isct

            P = restrict(brefs, bcell, cell)
            Q = restrict(trefs, tcell, cell)

            qr = quadrule(biop, trefs, brefs, cell, qd, quadstrat)
            zlocal = cellinteractions(biop, trefs, brefs, cell, qr)
            zlocal = Q * zlocal * P'

            for i in 1 : numfunctions(trefs)
                for j in 1 : numfunctions(brefs)
                    for (m,a) in tad[p,i]
                        for (n,b) in bad[q,j]
                            store(a * zlocal[i,j] * b, m, n)
                        end # next basis function this cell supports
                    end # next test function this cell supports
                end # next refshape on basis side
            end # next refshape on test side

        end # next cell in intersection

        done += 1
        new_pctg = round(Int, done / todo * 100)
        if new_pctg > pctg + 9
            print(".")
            pctg = new_pctg
        end
    end # next cell in the test geometry

    println("")

end

function assemble_local_matched!(biop::LocalOperator, tfs::subdBasis, bfs::subdBasis, store;
    quadstrat=defaultquadstrat(biop, tfs, bfs))

    tels, tad = assemblydata(tfs)
    bels, bad = assemblydata(bfs)

    trefs = refspace(tfs)
    brefs = refspace(bfs)

    qd = quaddata(biop, trefs, brefs, tels, bels, quadstrat)
    for (p,cell) in enumerate(tels)

        qr = quadrule(biop, trefs, brefs, cell, qd, quadstrat)
        locmat = cellinteractions(biop, trefs, brefs, cell, qr)

        for i in 1 : size(locmat, 1), j in 1 : size(locmat, 2)
            for (m,a) in tad[p][i], (n,b) in bad[p][j]
                store(a * locmat[i,j] * b, m, n)
end end end end

function elementstree(elements)

    nverts = dimension(eltype(elements)) + 1
    ncells = length(elements)

    @assert !isempty(elements)
    P = eltype(elements[1].vertices)
    T = coordtype(eltype(elements))

    points = zeros(P, ncells)
    radii = zeros(T, ncells)

    for i in 1 : ncells

        verts = elements[i].vertices

        bary = verts[1]
        for j in 2:length(verts)
            bary += verts[j]
        end

        points[i] = bary / nverts
        for j in 1 : nverts
            radii[i] = max(radii[i], norm(verts[j]-points[i]))
        end
    end

    return Octree(points, radii)
end

"""
    assemble_local_mixed(biop::LocalOperator, tfs, bfs)

For use when basis and test functions are defined on different meshes
"""
function assemble_local_mixed!(biop::LocalOperator, tfs::Space, bfs::Space, store;
    quadstrat=defaultquadstrat(biop, tfs, bfs))

    tol = sqrt(eps(Float64))

    trefs = refspace(tfs)
    brefs = refspace(bfs)

    tels, tad = assemblydata(tfs)
    bels, bad = assemblydata(bfs)

    qd = quaddata(biop, trefs, brefs, tels, bels, quadstrat)

    # store the bcells in an octree
    tree = elementstree(bels)

    print("dots out of 10: ")
    todo, done, pctg = length(tels), 0, 0
    for (p,tcell) in enumerate(tels)

        tc, ts = boundingbox(tcell.vertices)
        pred = (c,s) -> boxesoverlap(c,s,tc,ts)

        for box in boxes(tree, pred)
            for q in box
                bcell = bels[q]

                if overlap(tcell, bcell)

                    isct = intersection(tcell, bcell)
                    for cell in isct
                        volume(cell) < tol && continue

                        P = restrict(brefs, bcell, cell)
                        Q = restrict(trefs, tcell, cell)

                        qr = quadrule(biop, trefs, brefs, cell, qd, quadstrat)
                        zlocal = cellinteractions(biop, trefs, brefs, cell, qr)
                        zlocal = Q * zlocal * P'

                        for i in 1 : numfunctions(trefs)
                            for j in 1 : numfunctions(brefs)
                                for (m,a) in tad[p,i]
                                    for (n,b) in bad[q,j]
                                        store(a * zlocal[i,j] * b, m, n)
                                    end # next basis function this cell supports
                                end # next test function this cell supports
                            end # next refshape on basis side
                        end # next refshape on test side

                    end # next cell in intersection
                end # if overlap
            end # next cell in the basis geometry
        end # next box in the octree

        done += 1
        new_pctg = round(Int, done / todo * 100)
        if new_pctg > pctg + 9
            print(".")
            pctg = new_pctg
        end
    end # next cell in the test geometry

    println("")
end

function cellinteractions_matched!(zlocal, biop, trefs, brefs, cell, qr)

    num_tshs = length(qr[1][3])
    num_bshs = length(qr[1][4])

    # zlocal = zeros(Float64, num_tshs, num_bshs)
    for q in qr

        w, mp, tvals, bvals = q[1], q[2], q[3], q[4]
        j = w * jacobian(mp)
        kernel = kernelvals(biop, mp)

        for n in 1 : num_bshs
            bval = bvals[n]
            for m in 1 : num_tshs
                tval = tvals[m]

                igd = integrand(biop, kernel, mp, tval, bval)
                zlocal[m,n] += j * igd
            end
        end
    end

    return zlocal
end

function cellinteractions(biop, trefs, brefs, cell, qr)

    num_tshs = length(qr[1][3])
    num_bshs = length(qr[1][4])

    zlocal = zeros(Float64, num_tshs, num_bshs)
    for q in qr

        w, mp, tvals, bvals = q[1], q[2], q[3], q[4]
        j = w * jacobian(mp)
        kernel = kernelvals(biop, mp)

        for m in 1 : num_tshs
            tval = tvals[m]

            for n in 1 : num_bshs
                bval = bvals[n]

                igd = integrand(biop, kernel, mp, tval, bval)
                zlocal[m,n] += j * igd

            end
        end
    end

    return zlocal
end

ext: jl


text:
"""
    Atoms

Atoms objects characterize structural representations of molecular structures.
These representations can either be standardized forms
(like graph representation using an adjacency matrix, point cloud representation, etc),
or be custom user-defined representations too.
"""
module Atoms

using ..ChemistryFeaturization.AbstractType: AbstractAtoms

import ..ChemistryFeaturization.elements
elements(a::AbstractAtoms) = throw(MethodError(elements, a))
export elements

include("atomgraph.jl")
export AtomGraph, visualize, elements

end


text:
using MacroTools: @capture, prewalk
using Calculus: differentiate

export AnalyticFunction, @analytic

import FunctionWrappers
import FunctionWrappers: FunctionWrapper
const ScalarFun{T} = FunctionWrapper{T, Tuple{T}}

"""
`F64fun`: `FunctionWrapper` to wrap many different potentials within a single
type. Can be used as  `F64fun(::Function)` or as
`F64fun(::AnalyticFunction)`
"""
const F64fun = ScalarFun{Float64}

"""
take an expression of the form `r -> f(r)` and return the expression
`r -> f'(r)`
"""
function fdiff( ex )
   @assert @capture(ex, var_ -> expr_)
   return :(  $var -> $(differentiate(expr.args[2], var)) )
end

"""
auxiliary function to allow expression substitution from `@Analytic`;
see `?@Analytic` for more detail
"""
function substitute(args)
    subs = Dict{Symbol,Union{Expr, Symbol}}()
    for i = 2:length(args)
        @assert @capture(args[i], var_ = sub_)
        subs[var] = sub
    end
    prewalk(expr -> get(subs, expr, expr), args[1])
end

# ================== Analytical Potentials ==========================

# documentation attached below
@pot struct AnalyticFunction{F0,F1,F2} <: PairPotential
   f::F0
   f_d::F1
   f_dd::F2
end

"""
`struct AnalyticFunction`: described an analytic function, allowing to
evaluate at least 2 derivatives.

Formally, `AnalyticFunction <: PairPotential`, which simplifies dispatch,
but it an `AnalyticFunction` should normally **not be used as a Calculator**!
This type hierarchy may need to be revisited.

### Usage:
```julia
lj = @analytic r -> r^(-12) - 2.0*r^(-6)
morse = let A = 4.0, r0 = 1.234
   @analytic r -> exp(-2.0*A*(r/r0-1.0)) - 2.0*exp(-A*(r/r0-1.0))
end

To create a "wrapped" AnalyticFunction{F64fun, ...} , use
```
lj_wrapped = F64fun(lj)
```

if an formula for an analytic function contains a sub-expression that
occurs multiple time, then this can be constructed as follows:
```
V = @analytic( r -> exp(s) * s, s = r^2 )
# is the same as
V = @analytic r -> exp(r^2) * r^2
```
"""
AnalyticFunction

const WrappedAnalyticFunction = AnalyticFunction{F64fun, F64fun, F64fun}

F64fun(p::AnalyticFunction) =
   AnalyticFunction(F64fun(p.f), F64fun(p.f_d), F64fun(p.f_dd))

evaluate(p::AnalyticFunction, r::Number) = p.f(r)
evaluate_d(p::AnalyticFunction, r::Number) = p.f_d(r)
evaluate_dd(p::AnalyticFunction, r::Number) = p.f_dd(r)
cutoff(V::AnalyticFunction) = Inf

"""
`@analytic`: generate C2 function from symbol
"""
macro analytic(args...)
   fexpr = substitute(args)
   quote
      AnalyticFunction(
        $(Base.FastMath.make_fastmath(esc(fexpr))),
        $(Base.FastMath.make_fastmath(esc(fdiff(fexpr)))),
        $(Base.FastMath.make_fastmath(esc(fdiff(fdiff(fexpr)))))
      )
   end
end

# TODO: this is a hack to make tight-binding work; but it should be reconsidered
# right now I am thinking it is actually ok as is!
evaluate(p::AnalyticFunction, r, R) = evaluate(p, r)
evaluate_d(p::AnalyticFunction, r, R) = evaluate_d(p, r)


text:
# TODO: Make everything work with Float64 as well as BigFloat
# Precompile the module
__precompile__(true)
"""
## References:
 1. P. Deuflhard:
     Newton Methods for Nonlinear Problems. -
     Affine Invariance and Adaptive Algorithms.
     Series Computational Mathematics 35, Springer (2004)
 2. U. Nowak, L. Weimann:
     A Family of Newton Codes for Systems of Highly Nonlinear
     Equations - Algorithm, Implementation, Application.
     ZIB, Technical Report TR 90-10 (December 1990)

Currently the following two solvers are implemented:
- NLEQ1: Damped Newton-algorithm for systems of highly nonlinear
    equations - damping strategy due to Ref. (1).
- NLEQ2: Damped Newton-algorithm with rank strategy for systems of
    highly nonlinear equations - damping strategy due to Ref.(1).

## Licence
The MIT License (MIT)

Copyright (c) 2016 Vishal Sontakke

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

## Warranty
This code has been tested up to a certain level. Defects and
weaknesses, which may be included in the code, do not establish
any warranties by the author. The author does not take over any liabilities
which may follow from acquisition or application of this code.
"""
module NewtonMethods

using ForwardDiff

# Export the required methods
export nleq1, nleq2, OptionsNLEQ, clearWorkspace, clearAllWorkspaces

# Export the required options
export OPT_RTOL, OPT_QSUCC, OPT_MODE, OPT_JACGEN, OPT_JACFCN, OPT_MSTOR, OPT_ML,
    OPT_MU, OPT_ISCAL, OPT_PRINTWARNING, OPT_PRINTITERATION, OPT_PRINTIOWARN,
    OPT_PRINTIOMON, OPT_PRINTIOSOL, OPT_PRINTSOLUTION, OPT_NONLIN, OPT_QRANK1,
    OPT_QORDI, OPT_QSIMPL, OPT_NOROWSCAL, OPT_BOUNDEDDAMP, OPT_IORMON, OPT_NITMAX,
    OPT_FCBAND, OPT_SIGMA, OPT_SIGMA2, OPT_AJDEL, OPT_AJMIN, OPT_ETADIF, OPT_ETAINI,
    OPT_NBROY, OPT_FCSTART, OPT_FCMIN, OPT_IRANK, OPT_COND, OPT_STORE

# Export the required statistics
export STATS_NITER, STATS_NCORR, STATS_NFCN, STATS_NFCNJ, STATS_NJAC, STATS_NREJR1,
    STATS_XSCAL, STATS_RTOL, STATS_XITER, STATS_NATLEVEL, STATS_SIMLEVEL,
    STATS_STDLEVEL, STATS_PRECISION, STATS_DAMPINGFC, STATS_NEW, STATS_ICONV,
    STATS_CONV, STATS_SUMX, STATS_DLEVF, STATS_SUBCOND, STATS_SENS, STATS_IFAIL
# Include common files
include("Options.jl")
include("Constants.jl")
include("Common.jl")
include("Jacobian.jl")
include("SolverSpecific.jl")

global wkNLEQ1 = OptionsNLEQ()
global wkNLEQ2 = OptionsNLEQ()

"""
function clearWorkspace(name::AbstractString)

Function used to clear a workspace.

## Input parameter
`name` is the string which signifies which workspace is to be cleared.

For clearing the workspace related to NLEQ1 the following three strings are valid:
`NLEQ1`, `nleq1`, and `Nleq1`.

For clearing the workspace related to NLEQ2 the following three strings are valid:
`NLEQ2`, `nleq2`, and `Nleq2`
"""
function clearWorkspace(name::AbstractString)
    if name == "NLEQ1" || name == "nleq1" || name == "Nleq1"
        empty!(wkNLEQ1.options)
    elseif name == "NLEQ2" || name == "nleq2" || name == "Nleq2"
        empty!(wkNLEQ2.options)
    else
        println("Invalid option. Please specify the correct argument.")
        println("\"Nleq1\" or \"NLEQ1\" or \"nleq1\" - to clear the workspace for nleq1 function")
        println("\"Nleq2\" or \"NLEQ2\" or \"nleq2\" - to clear the workspace for nleq2 function")
    end
end

"""
function clearWorkspace()

Function used to clear all workspaces.
"""
function clearWorkspace()
    empty!(wkNLEQ1.options)
    empty!(wkNLEQ2.options)
end

# Include the solver specific files
include("NLEQ1.jl")
include("NLEQ2.jl")

end