text:
export draw_basin!, basins_of_attraction, automatic_Δt_basins

"""
    basins_of_attraction(grid::Tuple, ds::DynamicalSystem; kwargs...) -> basins, attractors
Compute an estimate of the basins of attraction of a dynamical system `ds` on
a partitioning of the state space given by `grid`. The method has been
inspired by the 2D grid approach developed by Nusse & Yorke [^Yorke1997].
It works _without_ knowledge of where attractors are; it identifies them automatically.

The dynamical system can be:
* An actual `DiscreteDynamicalSystem` or `ContinuousDynamicalSystem`.
* A Poincaré map of `ContinuousDynamicalSystem`: [`poincaremap`](@ref).
* A stroboscopic map, i.e. a periodically forced `ContinuousDynamicalSystem` (see examples
  for this particular application).

`grid` is a tuple of ranges defining the grid of initial conditions, for example
`grid = (xg, yg)` where `xg = yg = range(-5, 5; length = 100)`. The grid is not necessarilly
of the same dimension as the state space, attractors can be found in lower dimensional
projections.

The output `basins` is an integer-valued array on the `grid`, with its entries labelling
which basin of attraction the given grid point belongs to.
The output `attractors` is a dictionary whose keys correspond to the attractor number and
the values contains the points of the attractors found.
Notice that for some attractors this list may be incomplete.

See also [`match_attractors!`](@ref), [`basin_fractions`](@ref), [`tipping_probabilities`](@ref).

[^Yorke1997]: H. E. Nusse and J. A. Yorke, Dynamics: numerical explorations Ch. 7, Springer, New York, 1997

## Keyword Arguments
* `Δt`: Approximate time step of the integrator, which is `1` for discrete systems.
  For continuous systems, an automatic value is calculated using [`automatic_Δt_basins`](@ref).
  See that function for more info.
* `T` : Period of the stroboscopic map, in case of a continuous dynamical system with periodic
  time forcing. This argument is incompatible with `Δt`.
* `idxs = 1:length(grid)`: This vector selects the variables of the system that will define the
  subspace the dynamics will be projected into.
* `complete_state = zeros(D-Dg)`: This argument allows setting the _remaining_ variables
  of the dynamical system state on each initial condition `u`, with `Dg` the dimension
  of the grid. It can be either a vector of length `D-Dg`, or a function `f(y)` that
  returns a vector of length `D-Dg` given the _projected_ initial condition on the grid `y`.
* `diffeq = NamedTuple()`: Keyword arguments propagated to [`integrator`](@ref). Only
  useful for continuous systems. It is **strongly recommended** to choose high accuracy
  solvers for this application, e.g. `diffeq = (alg = Vern9(), reltol = 1e-9, abstol = 1e-9)`.
* `mx_chk_att = 2`: A parameter that sets the maximum checks of consecutives hits of an attractor
  before deciding the basin of the initial condition.
* `mx_chk_hit_bas = 10` : Maximum check of consecutive visits of the same basin of attraction.
  This number can be increased for higher accuracy.
* `mx_chk_fnd_att = 100` : Maximum check of unnumbered cell before considering we have an attractor.
  This number can be increased for higher accuracy.
* `mx_chk_loc_att = 100` : Maximum check of consecutive cells marked as an attractor before considering
  that we have all the available pieces of the attractor.
* `mx_chk_lost` : Maximum check of iterations outside the defined grid before we consider the orbit
  lost outside. This number can be increased for higher accuracy. It defaults to `20` if no
  attractors are given (see discussion on refining basins), and to `1000` if attractors are given.
* `horizon_limit = 1e6` : If the norm of the integrator state reaches this limit we consider that the
  orbit diverges.
* `show_progress = true` : By default a progress bar is shown using ProgressMeter.jl.
* `attractors, ε`: See discussion on refining basins below.

## Description
`basins` has the following organization:
* The basins are coded in sequential order from 1 up to the number of attractors.
* If the trajectory diverges or converges to an attractor outside the defined grid it is
  numbered `-1`

`attractors` has the following organization:
* The keys of the dictionary correspond to the number of the attractor.
* The value associated to a key is a [`Dataset`](@ref) with the *guessed* location of the
  attractor on the state space.

The method starts by picking the first available initial condition on the grid not yet
numbered. The dynamical system is then iterated until one of the following conditions happens:
1. The trajectory hits a known attractor already numbered `mx_chk_att` consecutive times: the
   initial condition is numbered with the corresponding number.
1. The trajectory spends `mx_chk_lost` steps outside the defiend grid: the initial
   condition is set to -1.
1. The trajectory hits a known basin `mx_chk_hit_bas` times in a row: the initial condition
   belongs to that basin and is numbered accordingly.
1. The trajectory hits `mx_chk_fnd_att` times in a row an unnumbered cell: it is considered
   an attractor and is labelled with a new number.

Regarding performace, this method is at worst as fast as tracking the attractors.
In most cases there is a signicative improvement in speed.

Notice that in the case we have to project the dynamics on a lower dimensional space,
there are edge cases where the system may have two attractors
that are close on the projected space but are far apart in another dimension. They could
be collapsed or confused into the same attractor. This is a drawback of this method.

## Refining basins of attraction
Sometimes one would like to be able to refine the found basins of attraction by recomputing
`basins_of_attraction` on a smaller, and more fine-grained, `grid`. If however this
new `grid` does not contain the attractors, `basins_of_attraction` would (by default)
attribute the value `-1` to all grid points. For these cases, an extra search clause can
be provided by setting the keywords `attractors, ε`. The `attractors` is a dictionary
mapping attractor IDs to `Dataset`s (i.e., the same as the return value of
`basins_of_attraction`). The algorithm checks at each step whether the system state is
`ε`-close (Euclidean norm) to any of the given attractors, and if so it attributes the stating grid point
to the basin of the close attractor. By default `ε` is equal to the mean grid spacing.

A word of advice while using this method: in order to work properly, `ε` should be
about the size of a grid cell that has been used to compute the given `attractors`. It is
recomended to keep the same step size (i.e., use the same integrator) since it may have an
influence in some cases. This algorithm is usually slower than the method with the
attractors on the grid.
"""
function basins_of_attraction(grid::Tuple, ds;
        Δt=nothing, T=nothing, idxs = 1:length(grid),
        complete_state = zeros(eltype(get_state(ds)), length(get_state(ds)) - length(grid)),
        diffeq = NamedTuple(), kwargs...
        # `kwargs` tunes the basin finding algorithm, e.g. `mx_chk_att`.
        # these keywords are actually expanded in `_identify_basin_of_cell!`
    )

    @assert length(idxs) == length(grid)
    integ = ds isa PoincareMap ? ds : integrator(ds; diffeq)
    idxs = SVector(idxs...)

    D = length(get_state(integ))
    if complete_state isa AbstractVector && (length(complete_state) ≠ D-length(idxs))
        error("Vector `complete_state` must have length D-Dg!")
    end
    # Check logic for automatic Δt computation:
    fixed_solver = haskey(diffeq, :dt) && haskey(diffeq, :adaptive)
    if ds isa ContinuousDynamicalSystem && isnothing(Δt) && isnothing(T) && !fixed_solver
        Δt = automatic_Δt_basins(ds, grid; idxs, complete_state, diffeq)
        if get(kwargs, :show_progress, true)
            @info "Automatic Δt estimation yielded Δt = $(Δt)"
        end
    end
    return basins_of_attraction(grid, integ, Δt, T, idxs, complete_state, fixed_solver; kwargs...)
end

function basins_of_attraction(
        grid, integ, Δt, T, idxs::SVector, complete_state, fixed_solver; kwargs...
    )

    complete_and_reinit! = CompleteAndReinit(complete_state, idxs, length(get_state(integ)))
    get_projected_state = (integ) -> view(get_state(integ), idxs)
    MDI = DynamicalSystemsBase.MinimalDiscreteIntegrator
    if !isnothing(T)
        iter_f! = (integ) -> step!(integ, T, true)
    elseif (integ isa PoincareMap) || (integ isa MDI) || fixed_solver
        iter_f! = step!
    else # generic case
        iter_f! = (integ) -> step!(integ, Δt) # we don't have to step _exactly_ `Δt` here
    end
    bsn_nfo = draw_basin!(
        grid, integ, iter_f!, complete_and_reinit!, get_projected_state; kwargs...
    )
    return bsn_nfo.basin, bsn_nfo.attractors
end

# Estimate Δt
using LinearAlgebra

"""
    automatic_Δt_basins(ds, grid; kwargs...) → Δt
Calculate an optimal `Δt` value for [`basins_of_attraction`](@ref).
This is done by evaluating the dynamic rule `f` (vector field) at `N` randomly chosen
points of the grid. The average `f` is then compared with the diagonal length of a grid
cell and their ratio provides `Δt`.

Keywords `idxs, complete_state, diffeq` are exactly as in [`basins_of_attraction`](@ref),
and the keyword `N` is `5000` by default.

Notice that `Δt` should not be too small which happens typically if the grid resolution
is high. It is okay for [`basins_of_attraction`](@ref) if the trajectory skips a few cells.
But if `Δt` is too small the default values for all other keywords such
as `mx_chk_hit_bas` need to be increased drastically.

Also, `Δt` that is smaller than the internal step size of the integrator will cause
a performance drop.
"""
function automatic_Δt_basins(ds, grid;
        idxs = 1:length(grid), N = 5000, diffeq = NamedTuple(),
        complete_state = zeros(eltype(get_state(ds)), length(get_state(ds)) - length(grid))
    )

    if ds isa Union{PoincareMap, DiscreteDynamicalSystem}
        return 1
    end
    integ = integrator(ds; diffeq...)
    complete_and_reinit! = CompleteAndReinit(complete_state, idxs, length(get_state(integ)))

    steps = step.(grid)
    s = sqrt(sum(x^2 for x in steps)) # diagonal length of a cell
    indices = CartesianIndices(length.(grid))
    random_points = [generate_ic_on_grid(grid, ind) for ind in rand(indices, N)]
    dudt = 0.0
    udummy = copy(integ.u)
    for p in random_points
        complete_and_reinit!(integ, p)
        deriv = if integ.u isa SVector
            integ.f(integ.u, integ.p, 0.0)
        else
            integ.f(udummy, integ.u, integ.p, 0.0)
            udummy
        end
        dudt += norm(deriv)
    end
    return Δt = 10*s*N/dudt
end


text:
using Distributions

function residuals_subperiods(dfm::DynamicFactorModel, break_period, variable_index)  # gets residuals from OLS regression of x on factors for both subsamples for variable_index (i.e. column of the x matrix)
    #dfm1 = DynamicFactorModel(dfm.y[1:break_period], reshape(ones(break_period), (break_period, 1)), dfm.x[1:break_period, :], "ICp1")  # select number of factors according to ICp1
    #dfm2 = DynamicFactorModel(dfm.y[break_period+1: end], reshape(ones(length(dfm.y)-break_period), (length(dfm.y)-break_period, 1)), dfm.x[break_period+1:end, :], "ICp1")  # select number of factors according to ICp1
    x1, y1 = dfm.factors[1:break_period, 1:dfm.number_of_factors], dfm.x[1:break_period, variable_index]
    x2, y2 = dfm.factors[break_period+1:end, 1:dfm.number_of_factors], dfm.x[break_period+1:end, variable_index]
    #factors1, loadings1, number_of_factors1 = calculate_factors(dfm.x[1:break_period, :])
    #factors2, loadings2, number_of_factors2 = calculate_factors(dfm.x[break_period+1:end, :])
    #return( (dfm.x[1:break_period, variable_index] - (factors1[:, 1:dfm.number_of_factors]*loadings1[:, dfm.number_of_factors]')[:, variable_index]), (dfm.x[break_period+1:end, variable_index] - (factors2[:, 1:dfm.number_of_factors]*loadings2[:, dfm.number_of_factors]')[:, variable_index]) )  # TODO: number of factors is likely to be higher if a structural break is present. --> number of factors should be reestimated for subperiods
    return( ((eye(break_period) - x1*inv(x1'x1)x1')y1), ((eye(length(y2)) - x2*inv(x2'x2)x2')y2) )   # use OLS for subperiods as in Breitung Eickmeier 2011 or Factor Model as in the 2009 version of the same paper?
end

# static factor Chow-tests Breitung, Eickmeier (2011)

# LR test
function LR_test(dfm::DynamicFactorModel, break_period::Int64, variable_index::Int64)  # TODO: in 2009 version of their paper Breitung and Eickmeier note that we can also use residuals from OLS regressions because they are asymptotically the same as PCA
    T = size(dfm.x, 1)
    likelihood_ratio = T * (log(sum(dfm.factor_residuals[:, variable_index].^2)) - log(apply(+, map(res->sum(res.^2), residuals_subperiods(dfm, break_period, variable_index)))))
    return likelihood_ratio
end

function Wald_test(dfm::DynamicFactorModel, break_period::Int64, variable_index::Int64)
    design_matrix = hcat(dfm.factors[:, 1:dfm.number_of_factors], dfm.factors[:, 1:dfm.number_of_factors].*[zeros(break_period), ones(length(dfm.y)-break_period)])
    ols_estimates = inv(design_matrix'design_matrix)design_matrix'dfm.x[:, variable_index]
    residuals = dfm.x[:, variable_index] - design_matrix * ols_estimates
    # calculate variance-covariance matrix according to White(1980) TODO: maybe this should be maximum likelihood estimation of var-cov matrix (i.e. negative inverse of information matrix)
    coefficient_covariance = inv(design_matrix'design_matrix)*(design_matrix'*diagm(residuals.^2)*design_matrix)*inv(design_matrix'design_matrix)
    wald_stat = (ols_estimates[dfm.number_of_factors+1:end]'inv(coefficient_covariance[dfm.number_of_factors+1:end, dfm.number_of_factors+1:end])ols_estimates[dfm.number_of_factors+1:end])[1]
    return wald_stat
end

function LM_test(dfm::DynamicFactorModel, break_period::Int64, variable_index::Int64)
    design_matrix = hcat(dfm.factors[:, 1:dfm.number_of_factors], dfm.factors[:, 1:dfm.number_of_factors].*[zeros(break_period), ones(length(dfm.y)-break_period)])
    ols_estimates = inv(design_matrix'design_matrix)design_matrix'dfm.factor_residuals[:, variable_index]
    residuals = dfm.factor_residuals[:, variable_index] - design_matrix * ols_estimates
    R_squared = 1 - sum(residuals.^2) / sum(dfm.factor_residuals[:, variable_index].^2)
    lagrange_multiplier = size(dfm.x, 1) * R_squared
    return lagrange_multiplier
end

ext: jl


text:
using DFTK
using Test

if mpi_nprocs() == 1  # Direct minimisation does not yet support MPI
@testset "External potential from Fourier coefficients" begin
    lattice = [[10 0 0.]; [0 0 0]; [0 0 0]]

    pot(G) = G == 0 ? zero(G) : 1 / abs(G)
    C = 1
    α = 2
    terms = [Kinetic(),
             ExternalFromFourier(G -> pot(G[1])),
             LocalNonlinearity(ρ -> C * ρ^α)]
    model = Model(lattice; n_electrons=1, terms, spin_polarization=:spinless)

    Ecut = 15
    basis = PlaneWaveBasis(model; Ecut, kgrid=(1, 1, 1))
    scfres_dm = direct_minimization(basis, tol=1e-10)
    scfres_scf = self_consistent_field(basis, tol=1e-10)
    @test abs(scfres_scf.energies.total - scfres_dm.energies.total) < 1e-6
end
end


text:
@doc raw"""
    NegativeMultinomial(x₀::Real, p::AbstractVector)

Negative Multinomial distribution defined as
```math
    p(\boldsymbol{x}|x_0, \boldsymbol{p}) = \Gamma\left(\sum_{i=0}^M x_i \right)\frac{p_0^{x_0}}{\Gamma(x_0)}\prod_{i=1}^M \frac{p_i^{x_i}}{x_i!}
```
where $p_0= 1-\sum_{i=1}^M p_i$.

For a detailed understanding of this distribution, see "Negative multinomial distribution" - Sibuya et al. - 1964
"""
struct NegativeMultinomial{Tx₀<:Real,Tp<:AbstractVector} <:
       Distributions.DiscreteMultivariateDistribution
    x₀::Tx₀
    p::Tp
    function NegativeMultinomial(x₀::Real, p::AbstractVector)
        x₀ > 0 || throw(ArgumentError("x₀ has to be positive"))
        (all(>=(0), p) && sum(p) < 1) || throw(
            ArgumentError(
                "All p should be positive and their sum should be strictly smaller than 1",
            ),
        )
        return new{typeof(x₀),typeof(p)}(x₀, p)
    end
end

_p₀(d::NegativeMultinomial) = 1 - sum(d.p)

Distributions.params(d::NegativeMultinomial) = (d.x₀, d.p)

Base.eltype(::NegativeMultinomial) = Int

Base.length(d::NegativeMultinomial) = length(d.p)

function Distributions._rand!(
    rng::AbstractRNG, d::NegativeMultinomial, x::AbstractVector{<:Real}
)
    p₀ = _p₀(d)
    θ = rand(rng, Gamma(d.x₀, inv(p₀) - 1))
    λ = d.p * θ / (1 - p₀) # convert parameters to the scaled Poisson ones
    for i in eachindex(x)
        x[i] = rand(rng, Poisson(λ[i]))
    end
    return x
end

function Distributions._logpdf(d::NegativeMultinomial, x::AbstractVector)
    return loggamma(d.x₀ + sum(x)) + d.x₀ * log(_p₀(d)) - loggamma(d.x₀) +
           mapreduce(+, d.p, x) do pᵢ, xᵢ
        xᵢ * log(pᵢ) - logfactorial(xᵢ)
    end
end

Distributions.mean(d::NegativeMultinomial) = d.x₀ / _p₀(d) * d.p

function Distributions.var(d::NegativeMultinomial)
    p₀ = _p₀(d)
    x₀ = d.x₀
    return x₀ / p₀^2 * abs2.(d.p) + x₀ / p₀ * p
end

function Distributions.cov(d::NegativeMultinomial)
    p₀ = _p₀(d)
    x₀ = d.x₀
    return x₀ / p₀^2 * d.p * d.p' + x₀ / p₀ * Diagonal(d.p)
end

function Distributions.mgf(d::NegativeMultinomial, t::AbstractVector)
    return (_p₀(d) / (1 - dot(d.p, exp.(t))))^d.x₀
end

function Distributions.kldivergence(p::NegativeMultinomial, q::NegativeMultinomial)
    p₀ = _p₀(p)
    x₀ = p.x₀
    x₀ == q.x₀ || error(
        "KL divergence cannot be computed with different x₍", "p: $(p.x₀), q: $(q.x₀)"
    )
    return x₀ * log(p₀) - x₀ * log(_p₀(q)) +
           x₀ / p₀ * sum(1:length(p)) do i
        p.p[i] * (log(p.p[i]) - log(q.p[i]))
    end
end


text:
using Compat

### Finite elements used to approximate fin_el ###

"""
## FiniteElement

Abstract finite element type.

### Type
```julia
abstract type FiniteElement end
```

### Subtypes
```julia
* Line::FiniteElement          : 1D Line(nod, nodof)
* Triangle::FiniteElement      : 2D Triangle(nod, nodof)
* Quadrilateral::FiniteElement : 2D Quadrilateral(nod, nodof)
* Hexahedron::FiniteElement    : 3D Hexahedron(nod, nodof)
* Tetrahedron::FiniteElement   : 3D Tetrahedron(nod, nodof)
```
"""
abstract type FiniteElement end

"""
## Line (Interval)

1D type finite element

### Constructor
```julia
Line(nod, nodof)
Line(nodof)
```

### Arguments
```julia
* nod::Int       : Number of nodes for finite element, defaults to 2
* nodof::Int     : Number of degrees of freedom per node
```

### Related help
```julia
?FiniteElement      : Help on finite element types
```

"""
struct Line <: FiniteElement
  nod::Int
  nodof::Int
end
Line(nodof::Int) = Line(2, nodof)

"""
## Triangle

2D type finite element

### Constructor
```julia
Triangle(nod, nodof)
```

### Arguments
```julia
* nod::Int       : Number of nodes for finite element (3, 6, 10, 15)
* nodof::Int     : Number of degrees of freedom per node
```

### Related help
```julia
?FiniteElement      : Help on finite element types
```

"""
struct Triangle <: FiniteElement
  nod::Int
  nodof::Int
end

"""
## Quadrilateral

2D type finite element

### Constructor
```julia
Quadrilateral(nod, nodof)
```

### Arguments
```julia
* nod::Int       : Number of nodes for finite element (4, 8, 9)
* nodof::Int     : Number of degrees of freedom per node
```

### Related help
```julia
?FiniteElement      : Help on finite element types
```

"""
struct Quadrilateral <: FiniteElement
  nod::Int
  nodof::Int
end

"""
## Tetrahedron

3D type finite element

### Constructor
```julia
Tetrahedron(nod, nodof)
Tetrahedron(nodof)
```

### Arguments
```julia
* nod::Int       : Number of nodes for finite element (defaults to 4)
* nodof::Int     : Number of degrees of freedom per node
```

### Related help
```julia
?FiniteElement      : Help on finite element types
```

"""
struct Tetrahedron <: FiniteElement
  nod::Int
  nodof::Int
end
Tetrahedron(nodof::Int) = Tetrahedron(4, nodof)

"""
## hexahedron

3D type finite element

### Constructor
```julia
Hexahedron(nod, nodof)
```

### Arguments
```julia
* nod::Int       : Number of nodes for finite element (8, 14, 20)
* nodof::Int     : Number of degrees of freedom per node
```

### Related help
```julia
?FiniteElement      : Help on finite element types
```

"""
struct Hexahedron <: FiniteElement
  nod::Int
  nodof::Int
end


text:
function general_header(op::Symbol, n::Integer, es::Integer, arity::Integer, status = false)
  prefix = "  "
  output = status ? "int" : p(n,es)
  ident  = string(p(n,es,op), status ? "" : "_e")
  statusarg = status ? ["$(p(n,es)) * res"] : []

  args   = join(vcat(statusarg, ["const $(p(n,es)) $(sym[idx])" for idx = 1:arity]),", ")
  "$prefix $output $ident($args)"
end

const binops = Dict(:add => :+, :sub => :-, :mul => :*, :div => :/, :lt => :<, :lte => :(<=) )

function cplusplus_header(op::Symbol, n::Integer, es::Integer, arity::Integer, status = false)
  if op in keys(binops)
    "extern $(c(n,es)) operator $(binops[op])($(c(n,es)) lhs, $(c(n,es)) rhs)"
  elseif op == :addinv
    "extern $(c(n,es)) operator -($(c(n,es)) arg)"
  else
    arglist = join(["$(c(n,es)) $(sym[n])" for n = 1:arity], ", ")
    "extern $(c(n,es)) $op($arglist)"
  end
end

header_set(op::Symbol, arity::Integer) = [(n, es) -> general_header(op, n, es, arity),
                                          (n, es) -> general_header(op, n, es, arity, true)]
#headers for special functions
boolean_header_c(op::Symbol) =   (n, es) -> "  bool $(p(n,es,op))(const $(p(n,es)) a,const  $(p(n,es)) b)"
boolean_header_cpp(op::Symbol) = (n, es) -> "extern bool operator $(binops[op])(const $(c(n,es)) a, const $(c(n,es)) b)"

fdp_header_c(op::Symbol) =       (n, es) -> "  $(p(n,es)) $(p(n,es,op))(fdp_cache$(n)e$(es)_t *fc, const $(p(n,es)) a, const $(p(n,es)) b)"
fdp_header_cpp(op::Symbol) =     (n, es) -> "extern $(c(n,es)) fdp(fdp_cache$(n)e$(es) *fc, const $(c(n,es)) a, const $(c(n,es)) b)"

dpi_header_c(op::Symbol) =       (n, es) -> "  void $(p(n,es,op))(fdp_cache$(n)e$(es)_t *fc)"
dpi_header_cpp(op::Symbol) =     (n, es) -> "extern void dpi(fdp_cache$(n)e$(es) *fc)"

newc_header(op::Symbol) =        (n, es) -> "  fdp_cache$(n)e$(es)_t* newcache$(n)e$(es)()"

const ordered_oplist =  [:add, :mul, :sub, :addinv, :lt, :lte, :gt, :gte, :eq, :div,
                         :mulinv, :log2, :exp2, :fma, :fms, :nfma, :nfms, :fas, :fcp, :fdp, :dpi, :newc,
                         :sqrt, :log1p, :log10, :log, :exp, :pow, :logx, :sin, :cos, :atan, :atan2]

#the ops dictionary lists

const ops = Dict(#required for all posit implementations
                 :add    => header_set(:add, 2),
                 :mul    => header_set(:mul, 2),
                 :sub    => header_set(:sub, 2),
                 :addinv => header_set(:addinv, 1),
                 :lt     => [nothing, boolean_header_c(:lt)],
                 :lte    => [nothing, boolean_header_c(:lte)],
                 :gt     => [nothing, boolean_header_c(:gt)],
                 :gte    => [nothing, boolean_header_c(:gte)],
                 :eq     => [nothing, boolean_header_c(:eq)],
                 #full posit required functions
                 :div    => header_set(:div, 2),
                 :mulinv => header_set(:mulinv, 1),
                 :log2   => header_set(:log2, 1),
                 :exp2   => header_set(:exp2, 1),
                 :fma    => header_set(:fma, 3),
                 :fms    => header_set(:fms, 3),
                 :nfma   => header_set(:nfma, 3),
                 :nfms   => header_set(:nfms, 3),
                 :fas    => header_set(:fas, 3),
                 :fcp    => header_set(:fcp, 4),
                 :fdp    => [nothing, fdp_header_c(:fdp)],
                 :dpi    => [nothing, dpi_header_c(:dpi)],
                 :newc   => [nothing, newc_header(:newc)],
                 #extended functions,
                 :sqrt   => header_set(:sqrt, 1),
                 :log1p  => header_set(:log1p, 1),
                 :log    => header_set(:log, 1),
                 :log10  => header_set(:log10, 1),
                 :exp    => header_set(:exp, 1),
                 :pow    => header_set(:pow, 2),
                 :logx   => header_set(:logx, 2),
                 :sin    => header_set(:sin, 1),
                 :cos    => header_set(:cos, 1),
                 :atan   => header_set(:atan, 1),
                 :atan2  => header_set(:atan2, 2))

function writefndef(io, op, n, es, version)
  op_fn = ops[op][version]
  write(io, op_fn(n, es))
end

function generate_posit_ops_status_h(io, posit_defs)
  write(io, """

    #ifndef __POSIT_OPS_STATUS_H
    #define __POSIT_OPS_STATUS_H

    #include "posit.h"

    extern "C"{

    """)
  for op in ordered_oplist
    for n in sort(collect(keys(posit_defs)))
      for es in posit_defs[n]
        writefndef(io, op, n, es, 2)
        write(io, ";\n")
      end
    end

    write(io, "\n")
  end
  write(io, """
  }

  #endif
  """)
end

function generate_posit_ops_errno_h(io, posit_defs)
  write(io, """

    #ifndef __POSIT_OPS_ERRNO_H
    #define __POSIT_OPS_ERRNO_H

    #include "posit.h"

    extern "C"{

    """)
  for op in ordered_oplist
    for n in sort(collect(keys(posit_defs)))
      for es in posit_defs[n]
        if ops[op][1] != nothing
          writefndef(io, op, n, es, 1)
          write(io, ";\n")
        end
      end
    end

    write(io, "\n")
  end
  write(io, """
  }

  #endif
  """)
end


text:
using SaveREPL
using Test

printREPL()

saveREPL("hist.jl")

copyREPL()


text:
ARSim(x::Vector{Float64},e::AR,s=randn(size(x))) = clamp((1-e.ρ)*e.μ + e.ρ*x+s*e.σ,e.x[1],e.x[end])
function MarkovSim!(ID::Vector{Int},s::Stochastic,r::Vector{Float64}=rand(length(ID)))
    csT=  cumsum(s.T,2)
    for i = 1:length(ID)
        j = 1
        while csT[ID[i],j]<r[i]
            j+=1
        end
        ID[i] = j
    end
    return ID
end

(s::AR)(ID::Vector{Int},r::Vector{Float64}=rand(length(ID))) = MarkovSim!(ID,s,r)
(s::AR)(ID::Int,r::Vector{Float64}=rand(length(ID))) = MarkovSim!([ID],s,r)[1]
(s::Markov)(ID::Vector{Int},r::Vector{Float64}=rand(length(ID))) = MarkovSim!(ID,s,r)

max_stars_repo_name: ZacLN/HADSGE

max_issues_repo_name: ZacLN/HADSGE

max_forks_repo_name: ZacLN/HADSGE


text:
function initialize(obj::ModeHandler, arg0::jboolean, arg1::AbsoluteDate)
    return jcall(obj, "initialize", void, (jboolean, AbsoluteDate), arg0, arg1)
end


text:
function ca_add_b!(c::Array{Float64,1}, a::Array{Float64,1}, b::Array{Float64,1}; nthreads::Int64=4)
    n = length(a)
    p = ccall((:ca_add_b_rr_, spmatveclib), Int64, (Ptr{Int64}, Ptr{Int64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
                                                    &nthreads , &n        , c           , a           , b             )
end

max_stars_repo_name: WneleiGao/ElasticWave


text:
#*******************************************************************************
# PACKAGES
#*******************************************************************************
using LinearAlgebra

#*******************************************************************************
# OPTIMIZATION FUNCTION DEFNITIONS
#*******************************************************************************
# Create an abstract type for first-order descent methods
abstract type FirstOrderMethod end

#*******************************************************************************
# Gradient Descent Implementation
mutable struct GradientDescent <: FirstOrderMethod
    α0; γ
end

gd_init!(GD::GradientDescent) = GD

function gd_step!(GD::GradientDescent, ∇f, x, k)
    α = GD.α0*GD.γ^(k-1)
    return x - α*∇f(x)
end

#*******************************************************************************
# Nesterov Momentum Implementation
mutable struct NesterovMomentum <: FirstOrderMethod
    α; β; v
end

nm_init!(NM::NesterovMomentum) = NM

function nm_step!(NM::NesterovMomentum, ∇f, x)
    α, β, v = NM.α, NM.β, NM.v
    v[:] = β*v - α*∇f(x + β*v)
    return x + v
end

#*******************************************************************************
# Adam Implementation
mutable struct Adam <: FirstOrderMethod
    α; γv; γs; ϵ; k; v; s
end

adam_init!(A::Adam) = A

function adam_step!(A::Adam, ∇f, x)
    α, γv, γs, ϵ, k, s, v, g = A.α , A.γv, A.γs, A.ϵ, A.k, A.s, A.v, ∇f(x)
    v[:] = γv*v + (1-γv)*g
    s[:] = γs*s + (1-γs)*g.*g
    A.k = k += 1
    v̂ = v ./ (1 - γv^k)
    ŝ = s ./ (1 - γs^k)
    return x - α*v̂ ./ (sqrt.(ŝ) .+ ϵ)
end

#*******************************************************************************
# Function to optimize the first simple problem
function optimize_simple1(f, g, x_best, n)
    #GD = GradientDescent(0.001, 0.995)
    #NM = NesterovMomentum(0.0006, 0.7, [0.0, 0.0])
    A = Adam(0.6, 0.4, 0.999, 1e-8, 0.0, [0.0,0.0], [0.0,0.0])
    k = 1
    while count(f, g) < n
        #x_best = gd_step!(GD, g, x_best, k)
        #x_best = nm_step!(NM, g, x_best)
        x_best = adam_step!(A, g, x_best)
        k = k + 1
    end
    return x_best
end

#*******************************************************************************
# Function to optimize the second simple problem
function optimize_simple2(f, g, x_best, n)
    #GD = GradientDescent(0.01, 0.995)
    #NM = NesterovMomentum(0.001, 0.9, [0.0, 0.0])
    A = Adam(0.5, 0.4, 0.999, 1e-8, 0.0, [0.0,0.0], [0.0,0.0])
    k = 1
    while count(f, g) < n
        #x_best = gd_step!(GD, g, x_best, k)
        #x_best = nm_step!(NM, g, x_best)
        x_best = adam_step!(A, g, x_best)
        k = k + 1
    end
    return x_best
end

#*******************************************************************************
# Function to optimize the third simple problem
function optimize_simple3(f, g, x_best, n)
    #GD = GradientDescent(0.001, 0.995)
    NM = NesterovMomentum(0.0009, 0.9, [0.0, 0.0, 0.0, 0.0])
    #A = Adam(0.4, 0.6, 0.999, 1e-8, 0.0, [0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0])
    k = 1
    while count(f, g) < n
        #x_best = gd_step!(GD, g, x_best, k)
        x_best = nm_step!(NM, g, x_best)
        #x_best = adam_step!(A, g, x_best)
        k = k + 1
    end
    return x_best
end

#*******************************************************************************
# Function to optimize the first secret problem
function optimize_secret1(f, g, x_best, n, len)
    v = [0.0 for x = 1:len]; s = [0.0 for x = 1:len]
    σ = [0.1 for x = 1:len]
    #GD = GradientDescent(0.12, 0.999)
    A = Adam(0.16, 0.75, 0.9, 1e-8, 0.0, v, s)
    k = 1
    while count(f, g) < n
        #x_best = gd_step!(GD, g, x_best, k)
        x_best = adam_step!(A, g, x_best)
        k = k + 1
    end
    return x_best
end

#*******************************************************************************
# Function to optimize the second secret problem
function optimize_secret2(f, g, x_best, n, len)
    v = [0.0 for x = 1:len]; s = [0.0 for x = 1:len]
    A = Adam(0.2, 0.6, 0.999, 1e-8, 0.0, v, s)
    while count(f, g) < n
        x_best = adam_step!(A, g, x_best)
    end
    return x_best
end

#*******************************************************************************
# MAIN OPTIMIZATION FUNCTION
#*******************************************************************************
function optimize(f, g, x0, n, prob)
    if prob == "simple1"
        x_best = optimize_simple1(f, g, x0, n)
    elseif prob == "simple2"
        x_best = optimize_simple2(f, g, x0, n)
    elseif prob == "simple3"
        x_best = optimize_simple3(f, g, x0, n)
    else
        len = length(x0)
        if len == 50
            x_best = optimize_secret1(f, g, x0, n, len)
        else
            x_best = optimize_secret2(f, g, x0, n, len)
        end
    end

    return x_best
end


text:
module testUtils

using Base.Test
using DataFrames
using TimeData
include(string(Pkg.dir("AssetMgmt"), "/src/AssetMgmt.jl"))

println("\n Running utils tests\n")

testMus = [0.000399905;
           0.000883977;
           0.000967358;
           0.00215034 ;
           0.00132068 ;
           0.000903217;
           0.00130242 ;
           0.00116732 ;
           0.00245227 ;
           0.00213858 ;
           0.00135902 ;
           0.00142933 ;
           0.00117371 ;
           0.00228266 ;
           0.00116749 ];

testSigmas = [0.000317419;
              0.00363872 ;
              0.00512444 ;
              0.0175377  ;
              0.0165166  ;
              0.00427322 ;
              0.0159598  ;
              0.0265965  ;
              0.0255393  ;
              0.0304155  ;
              0.0288717  ;
              0.0336391  ;
              0.0311551  ;
              0.0259919  ;
              0.023314];

############################
## grossMomentsLogMoments ##
############################

## transformation from gross to log to gross should be identity

nAss = size(testMus, 1)
for ii=1:nAss
    muGross = 1 + testMus[ii]
    sigma = sqrt(testSigmas[ii])

    ## to log moments
    muLog, sigmaLog =
        AssetMgmt.grossRetMomentsToLogRetMoments(muGross, sigma)

    ## and back
    muGrossOut, sigmaOut =
        AssetMgmt.logRetMomentsToGrossRetMoments(muLog, sigmaLog)

    @test muGross == muGrossOut
    @test_approx_eq_eps sigma sigmaOut 1e-14
end

#######################
## scaling functions ##
#######################

## at least shouldn't throw errors
scaledMus, scaledSigmas =
    AssetMgmt.defaultMuSigmaScaling(testMus, testSigmas)

##############################
## symbol string conversion ##
##############################

testSymbs = [:hello; :world]
expOut = ["hello"; "world"]
actOut = AssetMgmt.symbToStr(testSymbs)
@test expOut == actOut

testStr = UTF8String["hello"; "world"]
expOut = [:hello; :world]
actOut = AssetMgmt.strToSymb(testStr)
@test actOut == expOut

## EMACS_STOPPER_EMACS_STOPPER_EMACS_STOPPER

##############
## old code ##
##############

## ##################
## ## column means ##
## ##################

## df = DataFrame(rand(4, 2))
## tmp = AssetMgmt.mean(df, 1)
## @test isa(tmp, DataFrame)

## df = DataFrame(rand(4, 1))
## AssetMgmt.mean(df, 1)

## df = DataFrame(rand(1, 8))
## AssetMgmt.mean(df, 1)

## ## throw error for row sums
## @test_throws AssetMgmt.mean(df, 2)

## df = DataFrame(rand(50, 4))
## tmp = AssetMgmt.cov(df)
## @test isa(tmp, DataFrame)

## #########
## ## cov ##
## #########

## covMatr = AssetMgmt.cov(df)

## #############
## ## corrcov ##
## #############

## AssetMgmt.corrcov(covMatr)

## ##############
## ## randWgts ##
## ##############

## wgts = AssetMgmt.randWgts(10, 8)
## AssetMgmt.chkEqualsOne(wgts)

## ######################
## ## composeDataFrame ##
## ######################

## df = DataFrame(rand(8, 4))
## vals = array(df)
## nams = names(df)
## rename!(df, names(df), [:a1 :a2 :a3 :a4])
## df2 = AssetMgmt.composeDataFrame(vals, [:a1 :a2 :a3 :a4])
## @test (df == df2)

## ###############################################
## ## check matching investment and return data ##
## ###############################################

## vals = AssetMgmt.randWgts(10, 5)
## inv = AssetMgmt.Investments(DataFrame(vals))
## rets = Timematr(DataFrame(rand(10, 5)))
## AssetMgmt.chkMatchInvData(inv, rets)

## ## test invalid names
## rets = Timematr(rand(10, 5), [:a1 :a2 :a3 :a4 :a5])
## @test_throws AssetMgmt.chkMatchInvData(inv, rets)

## ## test invalid index
## vals = AssetMgmt.randWgts(10, 5)
## inv = AssetMgmt.Investments(DataFrame(vals), [20:29])
## rets = Timematr(rand(10, 5))
## @test_throws AssetMgmt.chkMatchInvData(inv, rets)

end


text:
# Section 4 of
# M. Philippe, R. Essick, G. E. Dullerud and R. M. Jungers.
# "Stability of discrete-time switching systems with constrained switching sequences."
# Automatica, 72:242-250, 2016
# The JSR is 0.9748171979372074

using LightGraphs

const expected_ub = [1.0035568400762171 0.9820958128927598
                     0.9863261446338583 0.9753566329765574
                     0.9769402921375085 0.9751185378641121
                     0.9749706065652288 0.9751185378641121
                     0.9749659545861287 0.9749374081223462
                     0.9749138871770363 0.9748393637263851]
const ρAs = 2.618033988749896
const ratio = [2, ρAs, ρAs, ρAs, ρAs, ρAs] .^ (1 ./ (2:2:12))
const expected_lb = expected_ub ./ ratio

@testset "[PEDJ] Section 4" begin
    include(joinpath(dirname(dirname(pathof(HybridSystems))), "examples", "PEDJ16s4.jl"))

    @testset "p-radius with $algo" for algo in (VeroneseLift(), KroneckerLift())
        pρlb = [0.87907456299, 0.89372555286, 0.90488811163]
        pρub = [1.42237252156, 1.13683646451, 1.06232506675]
        for d in 1:3
            lb, ub = pradiusb(hs, 2d, algo)
            @test isapprox(lb, pρlb[d])
            @test isapprox(ub, pρub[d])
        end
    end

    function t(pair::Pair)
        edge = HybridSystems.edge_object(hs.automaton, pair.first, pair.second)
        id = first(keys(hs.automaton.Σ[edge]))
        return HybridSystems.LightTransition(edge, id)
    end
    sbp = periodicswitching(hs, t.([3 => 3]))
    @test sbp.growthrate == 0.9392550239418471
    snp = periodicswitching(hs, t.([3 => 1, 1 => 3, 3 => 1, 1 => 3, 3 => 3, 3 => 3, 3 => 3, 3 => 3]))
    @test snp.growthrate == 0.9728940109399586
    smp = periodicswitching(hs, t.([3 => 1, 1 => 3, 3 => 1, 1 => 2, 2 => 3, 3 => 3, 3 => 3, 3 => 3]))
    @test smp.growthrate == 0.9748171979372074

    hsm = mdependentlift(hs, 2)
    function t(pair::Pair)
        edge = Edge(pair)
        id = first(keys(hsm.automaton.Σ[edge]))
        return HybridSystems.LightTransition(edge, id)
    end
    msbp = periodicswitching(hsm, t.([5 => 5]))
    @test msbp.growthrate == 0.9392550239418471
    msnp = periodicswitching(hsm, t.([5 => 5, 5 => 3, 3 => 8, 8 => 5, 5 => 5, 5 => 5, 5 => 5, 5 => 5, 5 => 3, 3 => 8, 8 => 5, 5 => 5, 5 => 5, 5 => 5, 5 => 5, 5 => 3, 3 => 8, 8 => 5, 5 => 5, 5 => 5, 5 => 5, 5 => 5, 5 => 3, 3 => 8, 8 => 5, 5 => 5, 5 => 5, 5 => 5, 5 => 5]))
    @test msnp.growthrate == 0.9653214971459174
    msmp = periodicswitching(hsm, t.([5 => 3, 3 => 8, 8 => 3, 3 => 7, 7 => 2, 2 => 5, 5 => 5, 5 => 5]))
    @test msmp.growthrate == 0.9748171979372074

    @testset "CJSR with $factory" for factory in sdp_factories
        @testset "With $(s == hs ? "original system" : "2-dependent lift")" for s in (hs, hsm)
            m = s === hs ? 1 : 2
            @testset "SOS d=$d" for d in 1:(7-m)
                tol = ismosek(factory) ? 6e-4 : 1e-3
                lb, ub = soslyapb(s, d, factory=factory, tol=tol)
                @test log(lb) ≈ log(expected_lb[d, m]) atol=tol
                @test log(ub) ≈ log(expected_ub[d, m]) atol=tol
            end
            @testset "EXT d=$d" for d in 1:(7-m)
                @testset "l=$l" for l in 1:2
                    @testset "v_0=$v_0" for v_0 in states(s)
                        seq = sosbuildsequence(s, d, p_0=:Primal, v_0=v_0, niter=100)
                        psw = findsmp(seq)
                        @test psw !== nothing
                        if s === hsm
                            if d == 1
                                if v_0 == 5
                                    @test psw == msbp
                                else
                                    @test psw == msnp
                                end
                            else
                                # I initially got only msmp but Mosek now also gets msbp
                                @test psw == msmp || psw == msbp
                            end
                        elseif d <= 3
                            if 2 <= d && v_0 == 3
                                @test psw == sbp
                            else
                                @test psw == snp
                            end
                        else
                            @test psw == smp
                        end
                    end
                end
            end
        end
    end
end


text:
##########################################

"""
$(TYPEDEF)

Structure containing the segments of a polygonal body.
To create an instance a set of indices have to be passed on.

# Fields

$(TYPEDFIELDS)

"""
struct BodySegments2D
    "(x,y) for first set of vertices (beginning of segments)"
    ver1::SubArray
    "(x,y) for second set of vertices (end of segments)"
    ver2::SubArray
    "total number of segments"
    nsegm::Integer

    function BodySegments2D(idx1::Union{Vector{<:Integer},UnitRange{<:Integer}},
                            vertices::Array{<:Real,2})
        # @assert ndims(ver1)==2
        @assert size(vertices,2)==2
        # circular shift to get second set of indices
        idx2 = circshift(idx1,-1)
        ## first set of vertices
        ver1 = view(vertices,idx1,:)
        ## second set of vertices
        ver2 = view(vertices,idx2,:)
        nsegm = size(ver1,1)
        return new(ver1,ver2,nsegm)
    end
end

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

"""
$(TYPEDEF)

Structure containing a set of polygonal bodies described by their segments and all vertices.
To create an instance, input an array of vectors of indices
  (of vertices) for each body and the array of all the vertices.

# Fields

$(TYPEDFIELDS)

"""
struct MagPolygBodies2D
    "array of bodies defined by their vertices"
    bo::Vector{BodySegments2D}
    "array of all vertices for all bodies"
    allvert::Array{<:Real,2}

    function MagPolygBodies2D(bodyindices::Vector{<:Vector{<:Integer}},allvert::Array{<:Real,2})
        @assert size(allvert,2)==2
        N=length(bodyindices)
        bodies = Vector{BodySegments2D}(undef,N)
        for i=1:N
            bodies[i] = BodySegments2D(bodyindices[i],allvert)
        end
        return new(bodies,allvert)
    end
end

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

"""
$(TYPEDEF)

Structure containing the components of a magnetization vector,
  i.e., module, inclination and declination angles.

# Fields

$(TYPEDFIELDS)

"""
Base.@kwdef struct MagnetizVector
    "modulus "
    mod::Real
    "inclination in degrees"
    Ideg::Real
    "declination in degrees"
    Ddeg::Real
end

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


text:
# ---
# title: 1642. Furthest Building You Can Reach
# id: problem1642
# author: AquaIndigo
# date: 2020-11-26
# difficulty: Medium
# categories: Binary Search, Heap
# link: <https://leetcode.com/problems/furthest-building-you-can-reach/description/>
# hidden: true
# ---
#
# You are given an integer array `heights` representing the heights of
# buildings, some `bricks`, and some `ladders`.
#
# You start your journey from building `0` and move to the next building by
# possibly using bricks or ladders.
#
# While moving from building `i` to building `i+1` ( **0-indexed** ),
#
#   * If the current building's height is **greater than or equal** to the next building's height, you do **not** need a ladder or bricks.
#   * If the current building's height is **less than** the next building's height, you can either use **one ladder** or `(h[i+1] - h[i])` **bricks**.
#
# _Return the furthest building index (0-indexed) you can reach if you use the
# given ladders and bricks optimally._
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/27/q4.gif)
#
#
#
#     Input: heights = [4,2,7,6,9,14,12], bricks = 5, ladders = 1
#     Output: 4
#     Explanation: Starting at building 0, you can follow these steps:
#     - Go to building 1 without using ladders nor bricks since 4 >= 2.
#     - Go to building 2 using 5 bricks. You must use either bricks or ladders because 2 < 7.
#     - Go to building 3 without using ladders nor bricks since 7 >= 6.
#     - Go to building 4 using your only ladder. You must use either bricks or ladders because 6 < 9.
#     It is impossible to go beyond building 4 because you do not have any more bricks or ladders.
#
#
# **Example 2:**
#
#
#
#     Input: heights = [4,12,2,7,3,18,20,3,19], bricks = 10, ladders = 2
#     Output: 7
#
#
# **Example 3:**
#
#
#
#     Input: heights = [14,3,19,3], bricks = 17, ladders = 0
#     Output: 3
#
#
#
#
# **Constraints:**
#
#   * `1 <= heights.length <= 105`
#   * `1 <= heights[i] <= 106`
#   * `0 <= bricks <= 109`
#   * `0 <= ladders <= heights.length`
#
#
## @lc code=start
using LeetCode
using DataStructures

function furthest_building(heights::Vector{Int}, bricks::Int, ladders::Int)
    q = Int[]
    for i in 1:(length(heights) - 1)
        diff = heights[i + 1] - heights[i]
        if diff > 0
            if diff > bricks && ladders == 0
                return i - 1
            end
            if bricks >= diff
                heappush!(q, diff)
                bricks -= diff
            else
                if !isempty(q) && q[1] > diff
                    bricks += heappop!(q) - diff
                    heappush!(q, diff)
                end
                ladders -= 1
            end
        end
    end
    return length(heights) - 1
end
## @lc code=end


text:
# THIS IS NOT A FULL JAVA WRAPPER!
# Even though the name seems to imply it, it only provides some convenience
# definitions for use in the MOI wrapper.

## Type definitions
# These correspond to Java types (which are directly returned by the functions
# in this binder). In particular, the inheritance defined in Java is not brought
# back to Julia.
const Store = JavaObject{Symbol("org.jacop.core.Store")}

const BooleanVar = JavaObject{Symbol("org.jacop.core.BooleanVar")}
const IntVar = JavaObject{Symbol("org.jacop.core.IntVar")}
const SetVar = JavaObject{Symbol("org.jacop.set.core.SetVar")}
const FloatVar = JavaObject{Symbol("org.jacop.floats.core.FloatVar")}
const CircuitVar = JavaObject{Symbol("org.jacop.constraints.CircuitVar")}

const Constraint = JavaObject{Symbol("org.jacop.constraints.Constraint")}
const LinearInt = JavaObject{Symbol("org.jacop.constraints.LinearInt")}
const Alldifferent = JavaObject{Symbol("org.jacop.constraints.Alldifferent")}
const In = JavaObject{Symbol("org.jacop.constraints.In")}
const XeqC = JavaObject{Symbol("org.jacop.constraints.XeqC")}

# Unions of types to model Java type hierarchy.

const Variable = Union{
    BooleanVar,
    IntVar,
    SetVar,
    FloatVar,
    CircuitVar,
}

# Add a constraint to a store.

function jacop_add_constraint_to_store(store::Store, constraint::Constraint)
    jcall(store, "impose", Nothing, (Constraint,), constraint)
    return
end


text:
const DOC_DISTRIBUTIONS =
"""
In the case the predictions `ŷ` are continuous probability
distributions, such as `Distributions.Normal`, replace the above sum
with an integral, and interpret `p` as the probablity density
function. In case of discrete distributions over the integers, such as
`Distributions.Poisson`, sum over all integers instead of `C`.
"""
const WITH_L2NORM_CONTINUOUS =
    [@eval(Distributions.$d) for d in [
        :Chisq,
        :Gamma,
        :Beta,
        :Chi,
        :Cauchy,
        :Normal,
        :Uniform,
        :Logistic,
        :Exponential]]

const WITH_L2NORM_COUNT =
    [@eval(Distributions.$d) for d in [
        :Poisson,
        :DiscreteUniform,
        :DiscreteNonParametric]]

const WITH_L2NORM = vcat([UnivariateFinite, ],
                         WITH_L2NORM_CONTINUOUS,
                         WITH_L2NORM_COUNT)

const UD = Distributions.UnivariateDistribution

# ========================================================
# AGGREGATED MEASURES

# ---------------------------------------------------------
# AreaUnderCurve

# Implementation drawn from
# https://www.kaggle.com/c/microsoft-malware-prediction/discussion/76013.

struct AreaUnderCurve <: Aggregated end

metadata_measure(AreaUnderCurve;
                 human_name = "area under the ROC",
                 instances = ["area_under_curve", "auc"],
                 target_scitype           = FiniteArrMissing{2},
                 prediction_type          = :probabilistic,
                 orientation              = :score,
                 supports_weights         = false,
                 distribution_type        = UnivariateFinite)

const AUC = AreaUnderCurve
@create_aliases AreaUnderCurve

@create_docs(AreaUnderCurve,
body=
"""
Returns the area under the ROC ([receiver operator
characteristic](https://en.wikipedia.org/wiki/Receiver_operating_characteristic))

If `missing` or `NaN` values are present, use `auc(skipinvalid(yhat, y)...)`.

$INVARIANT_LABEL
""",
scitpye = DOC_FINITE_BINARY)

# core algorithm:
function _auc(::Type{P}, ŷ, y) where P<:Real # type of probabilities
    lab_pos = classes(ŷ)[2] # 'positive' label
    scores  = pdf.(ŷ, lab_pos)     # associated scores
    y_sort  = y[sortperm(scores)]  # sort by scores
    n       = length(y)
    n_neg   = 0  # to keep of the number of negative preds
    auc     = P(0)
    @inbounds for i in 1:n
        # y[i] == lab_p --> it's a positive label in the ground truth
        # in that case increase the auc by the cumulative sum
        # otherwise increase the number of negatives by 1
        δ_auc, δ_neg = ifelse(y_sort[i] == lab_pos, (n_neg, 0), (0, 1))
        auc   += δ_auc
        n_neg += δ_neg
    end
    n_pos = n - n_neg
    return auc / (n_neg * n_pos)
end

# Missing values not supported, but allow `Missing` in eltype, because
# `skipinvalid(yhat, y)` does not tighten the type. See doc string above.

call(::AUC, ŷ::ArrMissing{UnivariateFinite{S,V,R,P}}, y) where {S,V,R,P} =
    _auc(P, ŷ, y)

# corner case of UnivariateFinite's of mixed type
call(::AUC, ŷ::ArrMissing{UnivariateFinite}, y) where {S,V,R,P} =
    _auc(Float64, ŷ, y)

# ========================================================
# UNAGGREGATED MEASURES

# ---------------------------------------------------------------------
# LogScore

struct LogScore{R <: Real} <: Unaggregated
    tol::R
end
LogScore(;eps=eps(), tol=eps) = LogScore(tol)

metadata_measure(LogScore;
                 instances                = ["log_score", ],
                 target_scitype           = Union{
                     Arr{<:Union{Missing,Multiclass}},
                     Arr{<:Union{Missing,OrderedFactor}},
                     Arr{<:Union{Missing,Continuous}},
                     Arr{<:Union{Missing,Count}}},
                 prediction_type          = :probabilistic,
                 orientation              = :score,
                 distribution_type        = Union{WITH_L2NORM...})

@create_aliases LogScore

@create_docs(LogScore,
body=
"""
Since the score is undefined in the case that the true observation is
predicted to occur with probability zero, probablities are clamped
between `tol` and `1-tol`, where `tol` is a constructor key-word
argument.

If `p` is the predicted probability mass or density function
corresponding to a *single* ground truth observation `η`, then the
score for that example is

    log(clamp(p(η), tol), 1 - tol)

For example, for a binary target with "yes"/"no" labels, and
predicted probability of "yes" equal to 0.8, an observation of "no"
scores `log(0.2)`.

The predictions `ŷ` should be an array of `UnivariateFinite`
distributions in the case of `Finite` target `y`, and otherwise a
supported `Distributions.UnivariateDistribution` such as `Normal` or
`Poisson`.

See also [`LogLoss`](@ref), which differs only in sign.
""")

# for single finite observation:
single(c::LogScore, d::UnivariateFinite{S,V,R,P}, η::Label) where {S,V,R,P} =
    log(clamp(pdf(d, η), P(c.tol), P(1) - P(c.tol)))

# for a single infinite observation:
single(c::LogScore, d::Distributions.UnivariateDistribution, η::Real) =
    log(clamp(pdf(d, η), c.tol, 1 - c.tol))

# performant broadasting in case of UnivariateFiniteArray:
function call(c::LogScore,
              ŷ::UnivariateFiniteArray{S,V,R,P,N},
              y::ArrMissing{V,N},
              w::Union{Nothing,Arr{<:Real,N}}=nothing) where {S,V,R,P<:Real,N}
    unweighted = log.(clamp.(broadcast(pdf, ŷ, y), P(c.tol), P(1) - P(c.tol)))
    if w === nothing
        return unweighted
    end
    return w .* unweighted
end

# ---------------------------------------------------------------------
# LogLoss

struct LogLoss{R <: Real} <: Unaggregated
    tol::R
end
LogLoss(;eps=eps(), tol=eps) = LogLoss(tol)

metadata_measure(LogLoss;
                 instances                = ["log_loss", "cross_entropy"],
                 target_scitype           = Union{
                     Arr{<:Union{Missing,Multiclass}},
                     Arr{<:Union{Missing,OrderedFactor}},
                     Arr{<:Union{Missing,Continuous}},
                     Arr{<:Union{Missing,Count}}},
                 prediction_type          = :probabilistic,
                 orientation              = :loss,
                 distribution_type        = Union{WITH_L2NORM...})

const CrossEntropy = LogLoss
@create_aliases LogLoss

@create_docs(LogLoss,
body=
"""
For details, see [`LogScore`](@ref), which differs only by a sign.
""")

# for single finite observation:
single(c::LogLoss, d::UnivariateFinite{S,V,R,P}, η::Label) where {S,V,R,P} =
    -single(LogScore(tol=c.tol), d, η)

# for a single infinite observation:
single(c::LogLoss, d::Distributions.UnivariateDistribution, η::Real) =
    -single(LogScore(tol=c.tol), d, η)

# performant broadasting in case of UnivariateFiniteArray:
call(c::LogLoss,
     ŷ::UnivariateFiniteArray{S,V,R,P,N},
     y::ArrMissing{V,N},
     w::Union{Nothing,Arr{<:Real,N}}=nothing) where {S,V,R,P<:Real,N} =
    -call(LogScore(tol=c.tol), ŷ, y, w)

# -----------------------------------------------------
# BrierScore

struct BrierScore <: Unaggregated end

metadata_measure(BrierScore;
                 human_name = "Brier score (a.k.a. quadratic score)",
                 instances                = ["brier_score",],
                 target_scitype           = Union{
                     Arr{<:Union{Missing,Multiclass}},
                     Arr{<:Union{Missing,OrderedFactor}},
                     Arr{<:Union{Missing,Continuous}},
                     Arr{<:Union{Missing,Count}}},
                 prediction_type          = :probabilistic,
                 orientation              = :score,
                 distribution_type        = Union{WITH_L2NORM...})

@create_aliases BrierScore

@create_docs(BrierScore,
body=
"""
Convention as in $PROPER_SCORING_RULES

*Finite case.* If `p` is the predicted probability mass function for a
*single* observation `η`, and `C` all possible classes, then the
corresponding score for that observation is given by

``2p(η) - \\left(\\sum_{c ∈ C} p(c)^2\\right) - 1``

*Warning.* `BrierScore()` is a "score" in the sense that bigger is
better (with `0` optimal, and all other values negative). In Brier's
original 1950 paper, and many other places, it has the opposite sign,
despite the name. Moreover, the present implementation does not treat
the binary case as special, so that the score may differ in the binary
case by a factor of two from usage elsewhere.

*Infinite case.* Replacing the sum above with an integral does *not*
lead to the formula adopted here in the case of `Continuous` or
`Count` target `y`. Rather the convention in the paper cited above is
adopted, which means returning a score of

``2p(η) - ∫ p(t)^2 dt``

in the `Continuous` case (`p` the probablity density function) or

``2p(η) - ∑_t p(t)^2``

in the `Count` cae (`p` the probablity mass function).
""",
scitype=DOC_FINITE)

# calling on single finite observation:
function single(::BrierScore,
                d::UnivariateFinite{S,V,R,P},
                η::Label) where {S,V,R,P}
    levels = classes(d)
    pvec = broadcast(pdf, d, levels)
    offset = P(1) + sum(pvec.^2)
    return P(2) * pdf(d, η) - offset
end

# calling on a single infinite observation:
single(::BrierScore, d::Distributions.UnivariateDistribution, η::Real) =
    2*pdf(d, η) - Distributions.pdfsquaredL2norm(d)

# Performant broadcasted version in case of UnivariateFiniteArray:
function call(::BrierScore,
              ŷ::UnivariateFiniteArray{S,V,R,P,N},
              y::ArrMissing{V,N},
              w::Union{Nothing,Arr{<:Real,N}}=nothing) where {S,V,R,P<:Real,N}

    probs = pdf(ŷ, classes(first(ŷ)))
    offset = P(1) .+ vec(sum(probs.^2, dims=2))

    unweighted = P(2) .* broadcast(pdf, ŷ, y) .- offset

    if w === nothing
        return unweighted
    end
    return w.*unweighted
end

# -----------------------------------------------------
# BrierLoss

struct BrierLoss <: Unaggregated end

metadata_measure(BrierLoss;
                 human_name = "Brier loss (a.k.a. quadratic loss)",
                 instances                = ["brier_loss",],
                 target_scitype           = Union{
                     Arr{<:Union{Missing,Multiclass}},
                     Arr{<:Union{Missing,OrderedFactor}},
                     Arr{<:Union{Missing,Continuous}},
                     Arr{<:Union{Missing,Count}}},
                 prediction_type          = :probabilistic,
                 orientation              = :loss,
                 distribution_type        = Union{WITH_L2NORM...})

@create_aliases BrierLoss

@create_docs(BrierLoss,
body=
"""
For details, see [`BrierScore`](@ref), which differs only by a sign.
""")

# calling on single finite observations:
single(::BrierLoss, d::UnivariateFinite{S,V,R,P}, η::Label) where {S,V,R,P} =
    - single(BrierScore(), d, η)

# calling on single infinite observations:
single(::BrierLoss, d::Distributions.UnivariateDistribution, η::Real) =
    -single(BrierScore(), d, η)

# to get performant broadcasting in case of UnivariateFiniteArray
call(m::BrierLoss,
     ŷ::UnivariateFiniteArray{S,V,R,P,N},
     y::ArrMissing{V,N},
     w::Union{Nothing,Arr{<:Real,N}}=nothing) where {S,V,R,P<:Real,N} =
         - call(BrierScore(), ŷ, y, w)

# -----------------------------------------------------
# SphericalScore

struct SphericalScore{T<:Real} <: Unaggregated
    alpha::T
end
SphericalScore(; alpha=2) = SphericalScore(alpha)

metadata_measure(SphericalScore;
                 human_name               = "Spherical score",
                 instances                = ["spherical_score",],
                 target_scitype           = Union{
                     Arr{<:Union{Missing,Multiclass}},
                     Arr{<:Union{Missing,OrderedFactor}},
                     Arr{<:Union{Missing,Continuous}},
                     Arr{<:Union{Missing,Count}}},
                 prediction_type          = :probabilistic,
                 orientation              = :score,
                 distribution_type        = Union{WITH_L2NORM...})

@create_aliases SphericalScore

@create_docs(SphericalScore,
body=
"""
Convention as in $PROPER_SCORING_RULES: If `η` takes on a finite
number of classes `C` and ``p(η)` is the predicted probability for a
*single* observation `η`, then the corresponding score for that
observation is given by

``p(y)^α / \\left(\\sum_{η ∈ C} p(η)^α\\right)^{1-α} - 1``

where `α` is the measure parameter `alpha`.

$DOC_DISTRIBUTIONS

""")

# calling on single observations:
function single(s::SphericalScore,
                d::UnivariateFinite{S,V,R,P}, η::Label) where {S,V,R,P}
    α = s.alpha
    levels = classes(d)
    pvec = broadcast(pdf, d, levels)
    return (pdf(d, η)/norm(pvec, α))^(α - 1)
end

single(s::SphericalScore, d::Distributions.UnivariateDistribution, η::Real) =
    pdf(d, η)/sqrt(Distributions.pdfsquaredL2norm(d))

# to compute the α-norm along last dimension:
_norm(A::AbstractArray{<:Any,N}, α) where N =
    sum(x -> x^α, A, dims=N)^(1/α)

# Performant version in case of UnivariateFiniteArray:
function call(s::SphericalScore,
              ŷ::UnivariateFiniteArray{S,V,R,P,N},
              y::ArrMissing{V,N},
              w::Union{Nothing,Arr{<:Real,N}}=nothing) where {S,V,R,P<:Real,N}

    α = s.alpha
    alphanorm(A) = _norm(A, α)

    predicted_probs = pdf(ŷ, classes(first(ŷ)))

    unweighted = (broadcast(pdf, ŷ, y) ./ alphanorm(predicted_probs)).^(α - 1)

    if w === nothing
        return unweighted
    end
    return w.*unweighted
end

# ---------------------------------------------------------------------------
# Extra check for L2 norm based proper scoring rules

err_l2_norm(m) = ArgumentError(
    "Distribution not supported by $m. "*
    "Supported distributions are "*
    join(string.(map(s->"`$s`", WITH_L2NORM)), ", ", ", and "))

const ERR_UNSUPPORTED_ALPHA = ArgumentError(
    "Only `alpha = 2` is supported, unless scoring a `Finite` target. ")

# not for export:
const L2ProperScoringRules = Union{LogScore,
                                   LogLoss,
                                   BrierScore,
                                   BrierLoss,
                                   SphericalScore}

function extra_check(measure::L2ProperScoringRules, yhat, args...)

    D = nonmissing(eltype(yhat))
    D <: Distributions.Distribution || D <: UnivariateFinite ||
        (D = typeof(findfirst(x->!isinvalid(x), yhat)))
    D <: Union{Nothing, WITH_L2NORM...} ||
        throw(err_l2_norm(measure))

    if measure isa SphericalScore
        measure.alpha == 2 || throw(ERR_UNSUPPORTED_ALPHA)
    end

    return nothing
end


text:
int_rules_4_3_7 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.3.7*(d*trig)^m*(a+b*(c*tan)^n)^p =#
    @apply_utils Antiderivative(~(u') * (~a + ~(b') * tan(~(e') + ~(f') * ~x) ^ 2) ^ ~p, ~x) => Antiderivative(ActivateTrig(~u * (~a * sec(~e + ~f * ~x) ^ 2) ^ ~p), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && EqQ(~a, ~b)
    @apply_utils Antiderivative(~(u') * (~(b') * tan(~(e') + ~(f') * ~x) ^ ~n) ^ ~p, ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], (((~b * ff ^ ~n) ^ IntPart(~p) * (~b * tan(~e + ~f * ~x) ^ ~n) ^ FracPart(~p)) / (tan(~e + ~f * ~x) / ff) ^ (~n * FracPart(~p))) * Antiderivative(ActivateTrig(~u) * (tan(~e + ~f * ~x) / ff) ^ (~n * ~p), ~x)) <-- FreeQ([~b, ~e, ~f, ~n, ~p], ~x) && (Not(IntegerQ(~p)) && (IntegerQ(~n) && (EqQ(~u, 1) || MatchQ(~u, (~(d') * (~trig)(~e + ~f * ~x)) ^ ~(m') <-- FreeQ([d, m], ~x) && MemberQ([sin, cos, tan, cot, sec, csc], trig)))))
    @apply_utils Antiderivative(~(u') * (~(b') * (~(c') * tan(~(e') + ~(f') * ~x)) ^ ~n) ^ ~p, ~x) => (((~b) ^ IntPart(~p) * (~b * (~c * tan(~e + ~f * ~x)) ^ ~n) ^ FracPart(~p)) / (~c * tan(~e + ~f * ~x)) ^ (~n * FracPart(~p))) * Antiderivative(ActivateTrig(~u) * (~c * tan(~e + ~f * ~x)) ^ (~n * ~p), ~x) <-- FreeQ([~b, ~c, ~e, ~f, ~n, ~p], ~x) && (Not(IntegerQ(~p)) && (Not(IntegerQ(~n)) && (EqQ(~u, 1) || MatchQ(~u, (~(d') * (~trig)(~e + ~f * ~x)) ^ ~(m') <-- FreeQ([d, m], ~x) && MemberQ([sin, cos, tan, cot, sec, csc], trig)))))
    @apply_utils Antiderivative(1 / (~a + ~(b') * tan(~(e') + ~(f') * ~x) ^ 2), ~x) => ~x / (~a - ~b) - (~b / (~a - ~b)) * Antiderivative(sec(~e + ~f * ~x) ^ 2 / (~a + ~b * tan(~e + ~f * ~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && NeQ(~a, ~b)
    @apply_utils Antiderivative((~a + ~(b') * (~(c') * tan(~(e') + ~(f') * ~x)) ^ ~n) ^ ~p, ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], ((~c * ff) / ~f) * Subst(Antiderivative((~a + ~b * (ff * ~x) ^ ~n) ^ ~p / ((~c) ^ 2 + ff ^ 2 * (~x) ^ 2), ~x), ~x, (~c * tan(~e + ~f * ~x)) / ff)) <-- FreeQ([~a, ~b, ~c, ~e, ~f, ~n, ~p], ~x) && (IntegersQ(~n, ~p) || (IGtQ(~p, 0) || (EqQ((~n) ^ 2, 4) || EqQ((~n) ^ 2, 16))))
    @apply_utils Antiderivative((~a + ~(b') * (~(c') * tan(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => Unintegrable((~a + ~b * (~c * tan(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~e, ~f, ~n, ~p], ~x)
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) ^ ~m * (~a + ~(b') * (~(c') * tan(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], ((~c * ff ^ (~m + 1)) / ~f) * Subst(Antiderivative(((~x) ^ ~m * (~a + ~b * (ff * ~x) ^ ~n) ^ ~p) / ((~c) ^ 2 + ff ^ 2 * (~x) ^ 2) ^ (~m / 2 + 1), ~x), ~x, (~c * tan(~e + ~f * ~x)) / ff)) <-- FreeQ([~a, ~b, ~c, ~e, ~f, ~n, ~p], ~x) && IntegerQ(~m / 2)
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * tan(~(e') + ~(f') * ~x) ^ 2) ^ ~(p'), ~x) => With([ff = FreeFactors(sec(~e + ~f * ~x), ~x)], (1 / (~f * ff ^ ~m)) * Subst(Antiderivative(((-1 + ff ^ 2 * (~x) ^ 2) ^ ((~m - 1) / 2) * ((~a - ~b) + ~b * ff ^ 2 * (~x) ^ 2) ^ ~p) / (~x) ^ (~m + 1), ~x), ~x, sec(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && IntegerQ((~m - 1) / 2)
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * tan(~(e') + ~(f') * ~x) ^ ~n) ^ ~(p'), ~x) => With([ff = FreeFactors(sec(~e + ~f * ~x), ~x)], (1 / (~f * ff ^ ~m)) * Subst(Antiderivative(((-1 + ff ^ 2 * (~x) ^ 2) ^ ((~m - 1) / 2) * (~a + ~b * (-1 + ff ^ 2 * (~x) ^ 2) ^ (~n / 2)) ^ ~p) / (~x) ^ (~m + 1), ~x), ~x, sec(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && (IntegerQ((~m - 1) / 2) && IntegerQ(~n / 2))
    @apply_utils Antiderivative((~(d') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~a + ~(b') * (~(c') * tan(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => Antiderivative(ExpandTrig((~d * sin(~e + ~f * ~x)) ^ ~m * (~a + ~b * (~c * tan(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n], ~x) && IGtQ(~p, 0)
    @apply_utils Antiderivative((~(d') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~a + ~(b') * tan(~(e') + ~(f') * ~x) ^ 2) ^ ~p, ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], ((ff * (~d * sin(~e + ~f * ~x)) ^ ~m * (sec(~e + ~f * ~x) ^ 2) ^ (~m / 2)) / (~f * tan(~e + ~f * ~x) ^ ~m)) * Subst(Antiderivative(((ff * ~x) ^ ~m * (~a + ~b * ff ^ 2 * (~x) ^ 2) ^ ~p) / (1 + ff ^ 2 * (~x) ^ 2) ^ (~m / 2 + 1), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~d, ~e, ~f, ~m, ~p], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative((~(d') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~a + ~(b') * (~(c') * tan(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => Unintegrable((~d * sin(~e + ~f * ~x)) ^ ~m * (~a + ~b * (~c * tan(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(d') * cos(~(e') + ~(f') * ~x)) ^ ~m * (~a + ~(b') * (~(c') * tan(~(e') + ~(f') * ~x)) ^ ~n) ^ ~p, ~x) => (~d * cos(~e + ~f * ~x)) ^ FracPart(~m) * (sec(~e + ~f * ~x) / ~d) ^ FracPart(~m) * Antiderivative((sec(~e + ~f * ~x) / ~d) ^ -(~m) * (~a + ~b * (~c * tan(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative((~(d') * tan(~(e') + ~(f') * ~x)) ^ ~(m') * (~a + ~(b') * (~(c') * tan(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], ((~c * ff) / ~f) * Subst(Antiderivative((((~d * ff * ~x) / ~c) ^ ~m * (~a + ~b * (ff * ~x) ^ ~n) ^ ~p) / ((~c) ^ 2 + ff ^ 2 * (~x) ^ 2), ~x), ~x, (~c * tan(~e + ~f * ~x)) / ff)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x) && (IGtQ(~p, 0) || (EqQ(~n, 2) || (EqQ(~n, 4) || IntegerQ(~p) && RationalQ(~n))))
    @apply_utils Antiderivative((~(d') * tan(~(e') + ~(f') * ~x)) ^ ~(m') * (~a + ~(b') * (~(c') * tan(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => Antiderivative(ExpandTrig((~d * tan(~e + ~f * ~x)) ^ ~m * (~a + ~b * (~c * tan(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n], ~x) && IGtQ(~p, 0)
    @apply_utils Antiderivative((~(d') * tan(~(e') + ~(f') * ~x)) ^ ~(m') * (~a + ~(b') * (~(c') * tan(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => Unintegrable((~d * tan(~e + ~f * ~x)) ^ ~m * (~a + ~b * (~c * tan(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(d') * cot(~(e') + ~(f') * ~x)) ^ ~m * (~a + ~(b') * tan(~(e') + ~(f') * ~x) ^ ~(n')) ^ ~(p'), ~x) => (~d) ^ (~n * ~p) * Antiderivative((~d * cot(~e + ~f * ~x)) ^ (~m - ~n * ~p) * (~b + ~a * cot(~e + ~f * ~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~d, ~e, ~f, ~m, ~n, ~p], ~x) && (Not(IntegerQ(~m)) && IntegersQ(~n, ~p))
    @apply_utils Antiderivative((~(d') * cot(~(e') + ~(f') * ~x)) ^ ~m * (~a + ~(b') * (~(c') * tan(~(e') + ~(f') * ~x)) ^ ~n) ^ ~p, ~x) => (~d * cot(~e + ~f * ~x)) ^ FracPart(~m) * (tan(~e + ~f * ~x) / ~d) ^ FracPart(~m) * Antiderivative((tan(~e + ~f * ~x) / ~d) ^ -(~m) * (~a + ~b * (~c * tan(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative(sec(~(e') + ~(f') * ~x) ^ ~m * (~a + ~(b') * (~(c') * tan(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], (ff / ((~c) ^ (~m - 1) * ~f)) * Subst(Antiderivative(((~c) ^ 2 + ff ^ 2 * (~x) ^ 2) ^ (~m / 2 - 1) * (~a + ~b * (ff * ~x) ^ ~n) ^ ~p, ~x), ~x, (~c * tan(~e + ~f * ~x)) / ff)) <-- FreeQ([~a, ~b, ~c, ~e, ~f, ~n, ~p], ~x) && (IntegerQ(~m / 2) && (IntegersQ(~n, ~p) || (IGtQ(~m, 0) || (IGtQ(~p, 0) || (EqQ((~n) ^ 2, 4) || EqQ((~n) ^ 2, 16))))))
    @apply_utils Antiderivative(sec(~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * tan(~(e') + ~(f') * ~x) ^ ~n) ^ ~(p'), ~x) => With([ff = FreeFactors(sin(~e + ~f * ~x), ~x)], (ff / ~f) * Subst(Antiderivative(ExpandToSum(~b * (ff * ~x) ^ ~n + ~a * (1 - ff ^ 2 * (~x) ^ 2) ^ (~n / 2), ~x) ^ ~p / (1 - ff ^ 2 * (~x) ^ 2) ^ ((~m + ~n * ~p + 1) / 2), ~x), ~x, sin(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (IntegerQ((~m - 1) / 2) && (IntegerQ(~n / 2) && IntegerQ(~p)))
    @apply_utils Antiderivative(sec(~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * tan(~(e') + ~(f') * ~x) ^ ~n) ^ ~(p'), ~x) => With([ff = FreeFactors(sin(~e + ~f * ~x), ~x)], (ff / ~f) * Subst(Antiderivative((1 / (1 - ff ^ 2 * (~x) ^ 2) ^ ((~m + 1) / 2)) * ((~b * (ff * ~x) ^ ~n + ~a * (1 - ff ^ 2 * (~x) ^ 2) ^ (~n / 2)) / (1 - ff ^ 2 * (~x) ^ 2) ^ (~n / 2)) ^ ~p, ~x), ~x, sin(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && (IntegerQ((~m - 1) / 2) && (IntegerQ(~n / 2) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative((~(d') * sec(~(e') + ~(f') * ~x)) ^ ~(m') * (~a + ~(b') * (~(c') * tan(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => Antiderivative(ExpandTrig((~d * sec(~e + ~f * ~x)) ^ ~m * (~a + ~b * (~c * tan(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n], ~x) && IGtQ(~p, 0)
    @apply_utils Antiderivative((~(d') * sec(~(e') + ~(f') * ~x)) ^ ~m * (~a + ~(b') * tan(~(e') + ~(f') * ~x) ^ 2) ^ ~p, ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], ((ff * (~d * sec(~e + ~f * ~x)) ^ ~m) / (~f * (sec(~e + ~f * ~x) ^ 2) ^ (~m / 2))) * Subst(Antiderivative((1 + ff ^ 2 * (~x) ^ 2) ^ (~m / 2 - 1) * (~a + ~b * ff ^ 2 * (~x) ^ 2) ^ ~p, ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~d, ~e, ~f, ~m, ~p], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative((~(d') * sec(~(e') + ~(f') * ~x)) ^ ~(m') * (~a + ~(b') * (~(c') * tan(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => Unintegrable((~d * sec(~e + ~f * ~x)) ^ ~m * (~a + ~b * (~c * tan(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(d') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~a + ~(b') * (~(c') * tan(~(e') + ~(f') * ~x)) ^ ~n) ^ ~p, ~x) => (~d * csc(~e + ~f * ~x)) ^ FracPart(~m) * (sin(~e + ~f * ~x) / ~d) ^ FracPart(~m) * Antiderivative((sin(~e + ~f * ~x) / ~d) ^ -(~m) * (~a + ~b * (~c * tan(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x) && Not(IntegerQ(~m))
end

max_stars_repo_path: src/rules/4 Trig functions/4.3 Tangent/4.3.7 (d trig)^m (a+b (c tan)^n)^p.jl


text:
module FinEtoolsVibInFluidsTutorials

end # module


text:
module NorthernLights

using ..Ahorn, Maple

placements = Maple.NorthernLights

function Ahorn.canFgBg(effect::Maple.NorthernLights)
    return true, true
end

end


text:
# dependencies
using MPI
using Compat
using Printf
using LinearAlgebra: BlasFloat, BlasReal
using DistributedArrays, Distributed
using DistributedArrays: DArray, defaultdist

# enicode encoded array to input ccall
f_pchar(scope::Char) = transcode(UInt8, string(scope))
f_pstring(scope::String) = transcode(UInt8, string(scope))

# path to scalapack dynamic library
const libscalapack = "/usr/local/Cellar/scalapack/2.0.2_15/lib/libscalapack.dylib"
# depend on envionment
# const ScaInt = Int32
const ScaInt = Int64
const zero = convert(ScaInt, 0)
const one = convert(ScaInt, 1)

# array descriptor index
struct descriptor_index
    # params
    dtype::ScaInt
    ctxt::ScaInt
    m::ScaInt; n::ScaInt;
    mb::ScaInt; nb::ScaInt;
    rsrc::ScaInt; csrc::ScaInt;
    lld::ScaInt;
    # constructor
    descriptor_index() = new(1, 2, 3, 4, 5, 6, 7, 8, 9)
end
const desc_idx = descriptor_index()


text:
"""
    R-LWE is performed over rings of the form ℛ = 𝔽q[x]/Φₘ(x)
    where Φₘ (Φ_m) is the m-th cyclotomic polynomial (the security of RLWE
    relies on cyclotomicity).

    This module consists of various standard parameter choices to match
    other popular libraries and standards.
"""
module CryptParameters

    using GaloisFields
    import ..NegacyclicRing
    export StdDistribution, HEStd_uniform, HEStd_error, HEStd_ternary,
        StdSecurity, HEStd_128_classic, HEStd_192_classic, HEStd_256_classic

    """
        These parameters match PALISADE's BGV defaults. See PALISADE's
        elementfactory.cpp.
    """
    const palisade_parameters = Dict(
        # m => ℛ
        16 => NegacyclicRing{GaloisField(1099511627873), 8}(108163207722),
        1024 => NegacyclicRing{GaloisField(525313), 512}(513496),
        2048 => NegacyclicRing{GaloisField(34359724033), 1024}(7225104974),
        4096 => NegacyclicRing{GaloisField(1152921504606830593), 2048}(811032584449645127)
    )

    @enum(StdDistribution,
        HEStd_uniform,
        HEStd_error,
        HeStd_ternary)

    @enum(StdSecurity,
        HEStd_128_classic = 1,
        HEStd_192_classic = 2,
        HEStd_256_classic = 3
    )

    struct StdParam
        dist::StdDistribution
        n::Int
        security::StdSecurity
        ringDim::Int
    end
    StdParam((a,b,c,d),) = StdParam(a,b,c,d)

    # Table from http://homomorphicencryption.org/wp-content/uploads/2018/11/HomomorphicEncryptionStandardv1.1.pdf
    const std_n = [1024, 2048, 4096, 8192, 16384, 32768]
    const std_params = Dict(
        HEStd_uniform => Dict(
             1024 => ( 29,  21,  16),
             2048 => ( 56,  39,  31),
             4096 => (111,  77,  60),
             8192 => (220, 154, 120),
            16384 => (440, 307, 239),
            32768 => (880, 612, 478)
        ),
        HEStd_error => Dict(
             1024 => ( 29,  21,  16),
             2048 => ( 56,  39,  31),
             4096 => (111,  77,  60),
             8192 => (220, 154, 120),
            16384 => (440, 307, 239),
            32768 => (883, 613, 478)
        ),
        HeStd_ternary => Dict(
            1024 => ( 27,  19,  14),
            2048 => ( 54,  37,  29),
            4096 => (109,  75,  58),
            8192 => (218, 152, 118),
           16384 => (438, 305, 237),
           32768 => (881, 611, 476)
        )
    )

    function std_ring_dim(dist_type, minSecLevel, curLogQ)
        std_n[findfirst(n->std_params[dist_type][n][Int(minSecLevel)] >= curLogQ, std_n)]
    end
end


text:
using Corona
using LinearAlgebra
using Plots
using FileIO, JLD2
using TimeSeries
using LaTeXStrings
using Formatting
using CSV
using FFTW

## Decide region and title
region = "Laputa";
plt_title = "Laputa";

## What to do
display_fig = true;
save_fig = true;
show_window = false;
show_actions = true;

### Load configuration
dataconfig = Corona.DataConfig();
include(joinpath(pwd(),"config/data.jl"));
if isfile(joinpath(pwd(),"config/paths.jl"))
    include(joinpath(pwd(),"config/paths.jl"))
end;

## Select paths
fname = Corona.build_filename(dataconfig,region);
figdir = joinpath(dataconfig.paths["figs"],fname);
resdir = joinpath(dataconfig.paths["results"],fname);
mkpath(figdir);
fig_ext = ["pdf", "png"];

## Load actions
act = CSV.read(joinpath(dataconfig.paths["raw"],"GOV",fname,"actions.csv"));
act_colors =  [:darkolivegreen, :palegreen2, :khaki, :darkslategray4, :sienna4];

## Load data
data, da = load(joinpath(resdir,"da.jld2"), "data", "da");

## Compute residual
J_min = Corona.residual(Corona.baseline(da), relative=true);

## Set title
# title_str = format(plt_title*", J = {:.3f} %",J_min*100);
title_str = plt_title*", $(data.lastdate)";
# title_str = plt_title;

## Determine datarange
datarange = data.outbreakdate:Day(1):data.lastdate;
datarange_dt = DateTime(data.outbreakdate):step(da.time):DateTime(data.lastdate);

## Determine result range
resrange = data.outbreakdate:Day(1):data.lastdate+Day(14);
resrange_dt = DateTime(resrange.start):step(da.time):DateTime(resrange.stop);

## Plot results
let plt = Plots.Plot()

        scatter!(plt, DateTime.(timestamp(data.cases[datarange])), values(data.cases.Confirmed[datarange]),
                    color=:orange, label="C (Confirmed) - data", Atickfontsize=12, legend=:topleft,
                    minorticks=true);
        scatter!(plt, DateTime.(timestamp(data.cases[datarange])), values(data.cases.Deaths[datarange]),
                    color=:black, label="D (Deaths) - data");

        plot!(plt, da.u.I[resrange_dt] .+ da.u.R[resrange_dt] .+ da.u.D[resrange_dt], color=:orange, lw=3, label="C (Confirmed)");
        plot!(plt, da.u.I[resrange_dt], color=:red, lw=3, label="I (Infectives)");
        plot!(plt, da.u.R[resrange_dt], color=:green, lw=3, label="R (Recovered)");
        plot!(plt, da.u.D[resrange_dt], color=:black, lw=3, label="D (Deaths)");
        title!(plt, title_str);

        if show_actions
            for i=1:size(act,1)
                    vline!(plt, [DateTime(act[i,:Date])], label=act[i,:Action] , lw=2, ls=:dash, color = act_colors[i])
            end
        end

        show_window && plot!(plt, ylims(plt)[2] .* da.σ.Confirmed, fillrange=ylims(plt)[1],
                                α=0.2, color=:grey, label="control");

        display_fig && display(plt)
        if save_fig
            for ext in fig_ext
                savefig(plt, joinpath(figdir,"da."*ext))
            end
        end

        plt
end;

## Plot results in log scale
let plt = Plots.Plot()

        scatter!(plt, DateTime.(timestamp(data.cases[datarange])), values(data.cases.Confirmed[datarange]) .+ 1,
                    color=:orange, label="C (data)", Atickfontsize=12, legend=:none, yaxis=:log10);
        scatter!(plt, DateTime.(timestamp(data.cases[datarange])), values(data.cases.Deaths[datarange]) .+ 1,
                    color=:black, label="D (data)");

        plot!(plt, da.u.I[datarange_dt] .+ da.u.R[datarange_dt] .+ da.u.D[datarange_dt] .+ 1.0, color=:orange, lw=3, label="C");
        plot!(plt, da.u.I[datarange_dt] .+ 1.0, color=:red, lw=3);
        plot!(plt, da.u.R[datarange_dt] .+ 1.0, color=:green, lw=3);
        plot!(plt, da.u.D[datarange_dt] .+ 1.0, color=:black, lw=3);
        title!(plt, title_str);

        if show_actions
            for i=1:size(act,1)
                    vline!(plt, [DateTime(act[i,:Date])], label=act[i,:Action] , lw=2, ls=:dash, color = act_colors[i])
            end
        end

        show_window && plot!(plt, ylims(plt)[2] .* da.σ.Confirmed[datarange_dt] .+ ylims(plt)[1], fillrange=ylims(plt)[1],
                                α=0.2, color=:grey, label="control");

        display_fig && display(plt)
        if save_fig
            for ext in fig_ext
                savefig(plt, joinpath(figdir,"da_log."*ext))
            end
        end

        plt
end

## Plot β
let plt = Plots.Plot()

        plot!(plt, da.p.β[datarange_dt], color=:orange, lw=3, label="\\beta (Contact rate)", legend=:topright);

        if show_actions
            for i=1:size(act,1)
                    vline!(plt, [DateTime(act[i,:Date])], label=act[i,:Action] , lw=2, ls=:dash, color = act_colors[i])
            end
        end

        show_window && plot!(plt, ylims(plt)[2] .* da.σ.Confirmed, fillrange=ylims(plt)[1],
                                α=0.2, color=:grey, label="control");

        title!(plt, title_str);

        display_fig && display(plt)
        if save_fig
            for ext in fig_ext
                savefig(plt, joinpath(figdir,"da_beta."*ext))
            end
        end

        plt
end

## Plot δ and γ
let plt = Plots.Plot()

        plot!(plt, da.p.γ[datarange_dt], color=:blue, lw=3, label="\\gamma", legend=:topleft);
        plot!(plt, da.p.δ[datarange_dt], color=:black, lw=3, label="\\delta");

        show_window && plot!(plt, ylims(plt)[2] .* da.σ.Confirmed[datarange_dt], fillrange=ylims(plt)[1],
                                α=0.2, color=:grey, label="control");

        if show_actions
            for i=1:size(act,1)
                    vline!(plt, [DateTime(act[i,:Date])], label=act[i,:Action], lw=2, ls=:dash, color = act_colors[i])
            end
        end

        title!(plt, title_str);

        display_fig && display(plt)
        if save_fig
            for ext in fig_ext
                savefig(plt, joinpath(figdir,"da_delta_gamma."*ext))
            end
        end

        plt
end

## Plot sigma
let plt = Plots.Plot()

        plot!(plt, da.p.β[datarange_dt]./(da.p.γ[datarange_dt] .+ da.p.δ[datarange_dt]), color=:red, lw=3, label="\\sigma (Contact number)" , legend=:topright)
        ylims!(plt,0,3)

        plot!(plt, timestamp(da.p.β[datarange_dt]), ones(size(da.p.β[datarange_dt],1)) , fillrange=ylims(plt)[1],
                                α=0.2, color=:green, label=:none)

        show_window && plot!(plt, ylims(plt)[2] .* da.σ.Confirmed, fillrange=ylims(plt)[1],
                                α=0.2, color=:grey, label="control");

        if show_actions
            for i=1:size(act,1)
                    vline!(plt, [DateTime(act[i,:Date])], label=act[i,:Action], lw=2, ls=:dash, color = act_colors[i])
            end
        end

        title!(plt, title_str);

        display_fig && display(plt)
        if save_fig
            for ext in fig_ext
                savefig(plt, joinpath(figdir,"da_contact_number."*ext))
            end
        end

        plt
end;

## Compute eigenvalues of A and extract A22
Λ = TimeArray(timestamp(da.u), complex(zeros(size(da.u))), [:λ₁, :λ₂, :λ₃, :λ₄]);
A22 = TimeArray(timestamp(da.u), zeros(size(da.u,1)), [:A22]);

for i=1:size(da.u,1)
        A = Corona.dfdu(values(da.u)[i,:], values(da.p)[i,:])
        values(A22)[i] = A[2,2]
        values(Λ)[i,:] = eigvals(A)
end;

## Plot eigenvalues of A
let plt = Plots.Plot()
        plot!(plt, real.(Λ.λ₁)[datarange_dt], label="\\lambda 1", lw=3, color=:orange)
        plot!(plt, real.(Λ.λ₂)[datarange_dt], label="\\lambda 2", lw=3, color=:blue)
        plot!(plt, real.(Λ.λ₃)[datarange_dt], label="\\lambda 3", lw=3, color=:green)
        plot!(plt, real.(Λ.λ₄)[datarange_dt], label="\\lambda 4", lw=3, color=:red)

        ylims!(plt, -.1, 0.4)

        show_window && plot!(plt, ylims(plt)[2] .* da.σ.Confirmed, fillrange=ylims(plt)[1],
                                α=0.2, color=:grey, label="control");

        if show_actions
            for i=1:size(act,1)
                    vline!(plt, [DateTime(act[i,:Date])], label=act[i,:Action] , lw=2, ls=:dash, color = act_colors[i])
            end
        end

        title!(plt, title_str);

        display_fig && display(plt)
        if save_fig
            for ext in fig_ext
                savefig(plt, joinpath(figdir,"da_lambda."*ext))
            end
        end

        plt
end;

## Plot A22
let plt = Plots.Plot()
        plot!(plt, A22[datarange_dt], label="A22", lw=3, color=:blue)
        # ylims!(plt, -.01, 0.4)

        plot!(plt, timestamp(A22[datarange_dt]), ylims(plt)[1]*ones(size(A22[datarange_dt],1)) , fillrange=0.0,
                                α=0.2, color=:green, label=:none)

        show_window && plot!(plt, ylims(plt)[2] .* da.σ.Confirmed, fillrange=ylims(plt)[1],
                                α=0.2, color=:grey, label="control");

        if show_actions
            for i=1:size(act,1)
                    vline!(plt, [DateTime(act[i,:Date])], label=act[i,:Action] ,  lw=2, ls=:dash, color = act_colors[i])
            end
        end

        title!(plt, title_str);

        display_fig && display(plt)
        if save_fig
            for ext in fig_ext
                savefig(plt, joinpath(figdir,"da_A22."*ext))
            end
        end

        plt
end;

## Compute FFT of A22
freq = 24*(0:(size(A22,1)/2))/size(A22,1);
P_A22 = abs.(fft(values(A22)));

let plt = Plots.Plot()
        sticks!(plt, freq.^(-1), P_A22, lw=3, label="FFT(A22)",
                xaxis="T (days)", legend=:top, color=:blue)
        # ylims!(plt, -.01, 0.4)

        xlims!(plt,2,9)
        ylims!(plt,0,10)

        title!(plt, title_str);

        display_fig && display(plt)
        if save_fig
            for ext in fig_ext
                savefig(plt, joinpath(figdir,"da_FFT_A22."*ext))
            end
        end

        plt
end;

## Export simulation result to CSV
let tah = when(da.u, hour, 0)
    tad = TimeArray(Date.(timestamp(tah)), values(tah), colnames(tah))
    writetimearray(tad, joinpath(resdir,"da_SIRD.csv"))
end;

## Export parameters to CSV
let tah = when(da.p, hour, 0)
    tad = TimeArray(Date.(timestamp(tah)), values(tah), colnames(tah))
    writetimearray(tad, joinpath(resdir,"da_param.csv"))
end;


text:
module TEST_NRES
using NRES
using Test

#=
#   NRES.jl : Eit fokus: Representere NRES som container av conditional-celler av type node.
#
#   -> dimensjonalitet, type av conditionals, osv. er vidare definert av under-typer.
=#

@testset "test test" begin

end

end#module TEST_NRES

max_stars_repo_name: leikanger/NRES

max_issues_repo_name: leikanger/NRES

max_forks_repo_name: leikanger/NRES


text:
# Exercise 20-1
function avoids(word, banned)
    return length(intersect(word, banned)) == 0
end

# println(avoids("abcdefghi", "wxyz"))
# println(avoids("abcdefghi", "d"))


text:
@blueprint struct Readme <: MustacheBlueprint
    template::TemplateFile = TemplateFile(file="README.md")
    inline_badges::Bool = false
end

function blueprint_view(p::Readme, ctx::Context)
    strings = String[]
    return Dict(
        "HAS_CITATION" => has_blueprint(ctx, Citation) && get_blueprint(ctx, Citation).readme,
        "HAS_INLINE_BADGES" => !isempty(strings) && p.inline_badges,
        "PKG" => ctx.kwargs.name,
    )
end

badge_order() = [
    Documenter{GitHubActions},
    Documenter{GitLabCI},
    Documenter{TravisCI},
    GitHubActions,
    GitLabCI,
    TravisCI,
    AppVeyor,
    DroneCI,
    CirrusCI,
    Codecov,
    Coveralls,
]


text:
immutable PINningSynapseParameter
end

@withkw type PINningSynapse
    param::PINningSynapseParameter = PINningSynapseParameter()
    colptr::Vector{Int} # column pointer of sparse W
    I::Vector{Int}      # postsynaptic index of W
    W::Vector{Float}  # synaptic weight
    rI::Vector{Float} # postsynaptic rate
    rJ::Vector{Float} # presynaptic rate
    g::Vector{Float}  # postsynaptic conductance
    P::Vector{Float}  # <rᵢrⱼ>⁻¹
    q::Vector{Float}  # P * r
    f::Vector{Float}  # postsynaptic traget
    records::Dict = Dict()
end

function PINningSynapse(pre, post; σ = 1.5, p = 0.0, α = 1)
    w = σ / √(p * pre.N) * sprandn(post.N, pre.N, p)
    rowptr, colptr, I, J, index, W = dsparse(w)
    rI, rJ, g = post.r, pre.r, post.g
    P = α .* (I .== J)
    f, q = zeros(post.N), zeros(post.N)
    PINningSynapse(;@symdict(colptr, I, W, rI, rJ, g, P, q, f)...)
end

@replace function forward!(c::PINningSynapse, param::PINningSynapseParameter)
    fill!(q, zero(Float))
    fill!(g, zero(Float))
    @inbounds for j in 1:(length(colptr) - 1)
        rJj = rJ[j]
        for s = colptr[j]:(colptr[j+1] - 1)
            i = I[s]
            g[i] += W[s] * rJj
            q[i] += P[s] * rJj
        end
    end
end

@replace function plasticity!(c::PINningSynapse, param::PINningSynapseParameter, dt::Float, t::Float)
    C = 1 / (1 + dot(q, rI))
    @inbounds for j in 1:(length(colptr) - 1)
        for s in colptr[j]:(colptr[j+1] - 1)
            i = I[s]
            P[s] += -C * q[i] * q[j]
            W[s] += C * (f[i] - g[i]) * q[j]
        end
    end
end