text:
"""
    forwardsub(L,b)

Solve the lower-triangular linear system with matrix `L` and right-hand side
vector `b`.
"""
function forwardsub(L,b)

    n = size(L,1)
    x = zeros(n)
    for i = 1:n
        s = i > 1 ? sum( L[i,j]*x[j] for j=1:i-1 ) : 0
        x[i] = ( b[i] - s ) / L[i,i]
    end

    return x
end

"""
    backsub(U,b)

Solve the upper-triangular linear system with matrix `U` and right-hand side
vector `b`.
"""
function backsub(U,b)

    n = size(U,1)
    x = zeros(n)
    for i = n:-1:1
        s = i < n ? sum( U[i,j]*x[j] for j=i+1:n ) : 0
        x[i] = ( b[i] - s ) / U[i,i]
    end

    return x
end

"""
    lufact(A)

Compute the LU factorization of square matrix `A`, returning the factors.
"""
function lufact(A)

    n = size(A,1)
    L = Matrix(Diagonal(ones(n)))    # puts ones on diagonal
    U = copy(A)

    # Gaussian elimination
    for j = 1:n-1
        for i = j+1:n
            L[i,j] = U[i,j] / U[j,j]   # row multiplier
            U[i,j:n] = U[i,j:n] - L[i,j]*U[j,j:n]
        end
    end

    return L,triu(U)
end

ext: jl


text:
using FdeSolver
using Test
using Statistics

@testset "FdeSolver.jl" begin

    tSpan = [0, 25]
    y0 = [34, 6]
    β = [0.98, 0.99]
    par = [0.55, 0.028, 0.80, 0.024]

    function F(t, n, β, y, par)

        F1 = par[1] .* y[n, 1] .- par[2] .* y[n, 1] .* y[n, 2]
        F2 = - par[3] .* y[n, 2] .+ par[4] .* y[n, 1] .* y[n, 2]

        [F1, F2]

    end

    t, Yapp = FDEsolver(F, tSpan, y0, β, nothing, par)

    @test @isdefined(t)
    @test @isdefined(Yapp)
    @test abs((mean(Yapp) - 28.19203)) < 0.05

end


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

"""
    slice(instance, range)

Creates a new instance, with only a subset of the time periods.
This function does not modify the provided instance. The initial
conditions are also not modified.

Example
-------

    # Build a 2-hour UC instance
    instance = UnitCommitment.read_benchmark("test/case14")
    modified = UnitCommitment.slice(instance, 1:2)

"""
function slice(
    instance::UnitCommitmentInstance,
    range::UnitRange{Int},
)::UnitCommitmentInstance
    modified = deepcopy(instance)
    modified.time = length(range)
    modified.power_balance_penalty = modified.power_balance_penalty[range]
    modified.reserves.spinning = modified.reserves.spinning[range]
    for u in modified.units
        u.max_power = u.max_power[range]
        u.min_power = u.min_power[range]
        u.must_run = u.must_run[range]
        u.min_power_cost = u.min_power_cost[range]
        u.provides_spinning_reserves = u.provides_spinning_reserves[range]
        for s in u.cost_segments
            s.mw = s.mw[range]
            s.cost = s.cost[range]
        end
    end
    for b in modified.buses
        b.load = b.load[range]
    end
    for l in modified.lines
        l.normal_flow_limit = l.normal_flow_limit[range]
        l.emergency_flow_limit = l.emergency_flow_limit[range]
        l.flow_limit_penalty = l.flow_limit_penalty[range]
    end
    for ps in modified.price_sensitive_loads
        ps.demand = ps.demand[range]
        ps.revenue = ps.revenue[range]
    end
    return modified
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   ModiaMath.DAE (ModiaMath/DAE/_module.jl)
#

# Utility functions used below
"""
    defaultVariableName(model, vcat, vindex)

Return default names for the variables (e.g. `x[1]`)
"""
defaultVariableName(model::Any, vcat::VariableCategory, vindex::Int) = vcat == Category_X ? "x[" * string(vindex) * "]" :
                                                                      (vcat == Category_W ? "w[" * string(vindex) * "]" :
                                                                                "der(x[" * string(vindex) * "])" )

#nameVector(name, nNames::Int) = [Symbol(name, "[", i, "]") for i = 1:nNames]
#fcNameVector(fc, names::AbstractVector) = [ Symbol(fc, "(", names[i], ")") for i in eachindex(names)]

nameVector(name, nNames::Int) = [string(name, "[", i, "]") for i = 1:nNames]
fcNameVector(fc, names::AbstractVector) = [ string(fc, "(", names[i], ")") for i in eachindex(names)]

"""
    getVariableName(model, vcat, vindex, nx=0;
                    xNames   =nameVector("x", nx),
                    derxNames=fcNameVector("der", xNames),
                    wNames   =String[])

Given category `vcat` and index `vindex`of category,
return the full path name of the respective variable.
"""
function getVariableName(model::Any, vcat::VariableCategory, vindex::Int, nx=0;
                         xNames::Vector{String}=nameVector("x", nx),
                         derxNames::Vector{String}=fcNameVector("der", xNames),
                         wNames::Vector{String}=String[])
    if vcat == Category_X
        return xNames[vindex]
    elseif vcat == Category_DERX
        return derxNames[vindex]
    elseif vcat == Category_W
        return wNames[vindex]
    else
        error("Wrong argument vcat = ", vcat)
    end
end

"""
   getResultNames(model)

Return a vector of result names.
"""
function getResultNames(model::Any)
    sim::SimulationState = model.simulationState
    resultNames = ["time"; String[sim.getVariableName(model, Category_X, i)  for i = 1:sim.nx];
                          String[sim.getVariableName(sim, Category_DERX, i) for i = 1:sim.nx];
                          String[sim.getVariableName(sim, Category_W, i)    for i = 1:sim.nw] ]
    return resultNames
end

"""
    @enum StructureOfDAE
          DAE_LinearDerivativesAndConstraints
          DAE_ExplicitDerivatives
          DAE_NoSpecialStructure

Enumeration defining the structure of the DAE of the simulation model.
The following DAE structures are supported
(function `getModelResidues!(model, t, x, derx, r, w; simulation=true)` returns the residues `r`):

# ModiaMath.DAE_LinearDerivativesAndConstraints (default)

```math
\\begin{align}
     z &= f_z(x,t) \\\\
 0 = r &= \\left[ \\begin{array}{l}
                    M_d(x,t,z_i>0) \\cdot \\dot{x} + b_d(x,t,z_i>0) \\;\\; (= f_d) \\\\
                    f_c(x,t,z_i>0)
                  \\end{array} \\right] \\\\
     J &= \\left[ M_d;
                  \\frac{\\partial f_c}{\\partial x} \\right] \\; \\text{is regular (matrix is invertible)}
\\end{align}
```

where

```math
\\lim_{\\epsilon \\rightarrow 0} x(t_0 - \\epsilon) = x_0^{-}
```

Equations ``f_d`` are linear in the derivatives ``\\dot{x}``.
It is required that the Jacobian ``J`` is **regular**, that is the DAE has an index 1
(= by differentiating ``f_c`` once, the system can be transformed to an ODE).

Equations ``z=z(t)`` are zero-crossing functions. Whenever a ``z_i(t)`` crosses zero,
an event is triggered and simulation is halted. During an event, ``z_i > 0`` can
change its value. The equations above are solved with a fixed-point iteration scheme (= *event iteration*)
until ``z_i > 0`` does not change anymore. Afterwards, integration is
restarted and the Boolean variable ``z_{pos} = z_{i,ev} > 0`` keeps its value until
the next event occurs.

At an event instant some ``f_c`` equations might become ``f_d`` equations and vice versa.
The constraint equations ``f_c`` can be at any position of the residue vector `r` and at an event
instant ``f_c`` equations might become ``f_d`` equations and vice versa.
When instantiating a [`SimulationState`](@ref), the initial definition
of the constraint equations is provided with vector argument `is_constraint`.
Note, it is also possible to define time events, so triggering events at pre-defined time
instants, for example to model sampled data systems
(see [`ModiaMath.setNextEvent!`](@ref)).

Initial conditions ``x_{ev}^{-}`` must be provided before simulation can start (``x_{ev}^{-} = x_0^{-}``) or at
an event restart. They need
**not** to fulfill the constraint equations, so ``f_c (x_{ev}^{-},t_{ev} ) \\neq 0 `` is allowed.
If this is the case, initialization/re-initialization will simulate for an infinitesimal small time instant
so that ``x_{ev}^{-}`` changes discontinuously to ``x_{ev}^{+}`` with ``f_c (x_{ev}^{+},t_{ev} )=0``.
This is performed by *analytically* integrating over the initial time or the event time and
might imply to integrate over Dirac impulses (in case ``x`` is discontinuous at this time instant).
Under certain conditions a numerical approximation of the
mathematical (exact) solution is computed.

The derivative of the constraint equations ``f_c`` can be provided at
event restart. In this case ``\\frac{\\partial f_c}{\\partial x}`` is computed
with the provided ``\\dot{f}_c``, instead of computing it with finite differences
(which is numerically less reliable). In both cases ``\\dot{x}_{ev}^{+}`` is then
computed by solving a linear equation system with Jacobian ``J``.

# ModiaMath.DAE_ExplicitDerivatives

```math
\\begin{align}
     z &= f_z(x, t) \\\\
 0 = r &= - \\dot{x} + b_d(x,t,z_i>0) \\;\\; (= f_d)
\\end{align}
```
where

```math
x(t_0) = x_0
```

This is a special case of the first form, where no constraint equations are present
and the derivatives are explicit. This form is also called ODE (Ordinary Differential Equations
in state space form). It has the advantage that many ODE integration methods can be used
to solve this system.

Initialization and re-initialization is trivial, because ``x_{ev}^{+}`` is provided
as initial value or at event restart from the model and then:

```math
der(x_{ev})^{+} := f_d(0, x_{ev}^{+}, t_{ev})
```

Note, if a Dirac impulse occurs in the model, then this property has to be
handled inside the model and the result of analytically integrating over the event
instant must be returned as ``x_{ev}^{+}`` from the model.

# ModiaMath.DAE_NoSpecialStructure

```math
\\begin{align}
     z &= f_z(x,t) \\\\
 0 = r &= f(\\dot{x},x,t,z_i>0)
 \\end{align}
```

This is a general DAE. Initialization and re-initialization is performed by using an implicit Euler step.
When appropriately scaling `r`, using a step size that tends to zero and under further
assumptions, then this solution can be interpreted as analytically integrating over
the time instant. This might mean to integrate over
Dirac impulses (in case `x` is discontinuous at this time instant).
Since the selected step size is not close to zero, the implicit Euler step
will give a very rough approximation of the analytical integral.
A much better approximation is achieved with option
`DAE_LinearDerivativesAndConstraints` above where a step size of zero is used.

This structure is only provided for backwards compatibility.
It is numerically not reliable and should no longer be used.
"""
@enum StructureOfDAE DAE_LinearDerivativesAndConstraints DAE_ExplicitDerivatives DAE_NoSpecialStructure

"""
    simulationState = SimulationState(name, getModelResidues!, x_start,
                                      getVariableName; kwargs...)

Return a `simulationState` object that is described by a DAE
with one of the supported structures defined with enumeration [`StructureOfDAE`](@ref).

A `model` that shall be simulated with function
[`ModiaMath.simulate!`](@ref)`(model, ...)` is required to be defined as:

```julia
mutable struct ModelName <: ModiaMath.AbstractSimulationModel
    simulationState::ModiaMath.SimulationState

    # other definitions (e.g. parameters of model)
end
```

| Keyword arguments           | defaults                                           |
|:----------------------------|:---------------------------------------------------|
| structureOfDAE              | DAE_LinearDerivativesAndConstraints (see [`StructureOfDAE`](@ref)) |
| is_constraint               | fill(false, length(x_start))                       |
| has_constraintDerivatives   | false                                              |
| nz                          | 0                                                  |
| nw                          | 0                                                  |
| zDir                        | fill(0, nz)                                        |
| x_fixed                     | fill(false, length(x_start))                       |
| x_nominal                   | fill(NaN, length(x_start))                         |
| x_errorControl              | fill(true, length(x_start))                        |
| jac                         | nothing (not yet supported)                        |
| maxSparsity                 | 0.1                                                |
| hev                         | 1e-8  (only for DAE_NoSpecialStructure)            |
| scaleConstraintsAtEvents    | true  (not used; only for backwards compatibility) |
| nc                          | 0 (not used; only for backwards compatibility)     |
| getResultNames              | [`ModiaMath.getResultNames`](@ref)                 |
| storeResult!                | [`ModiaMath.storeRawResult!`](@ref)                |
| getResult                   | [`ModiaMath.getStringDictResult`](@ref)            |
| defaultTolerance            | 1e-4                                               |
| defaultStartTime            | 0.0                                                |
| defaultStopTime             | 1.0                                                |
| defaultInterval             | NaN                                                |
| defaultMaxStepSize          | NaN                                                |
| defaultMaxNumberOfSteps     | missing                                            |

# Required arguments

- `name::Union{AbstractString,Symbol}`: Name of model

- `getModelResidues!::Function`: Function with arguments `(model,t,x,derx,r,w)` to
  compute the residues `r` and auxiliary variables `w` from time `t`, vector `x` and
  its time derivative `derx`.

- `x_start::Vector{Float64}`: Start values of `x`.

- `getVariableName::Function=`[`ModiaMath.defaultVariableName`](@ref): Function that returns the
  name of a variable, given its type and its index.

# Optional (keyword) arguments:

- `structureOfDAE::`[`StructureOfDAE`](@ref): Structure of DAE.

- `is_constraint::Vector{Bool}`: = true, `residue[i]` is a constraint equation `fc` or its derivative.
                                 = false, `residue[i]` is not a constraint equation.
                                 (is only used for `structureOfDAE = DAE_LinearDerivativesAndConstraints`)

- `has_constraintDerivatives::Vector{Bool}`: if [`ModiaMath.compute_der_fc`](@ref) returns true
                        and `is_constraint[i] = true`:
                         = true , `residue[i]` is the derivative of an `fc` equation
                         = false, `residue[i]` is an `fc` equation
                        (is only used for `structureOfDAE = DAE_LinearDerivativesAndConstraints`)

- `nz::Int`: Number of event indicators

- `nw::Int`: Number of auxiliary variables (Float64 variables that are additionally computed
             and stored at communication points, and where start values can be provided
             for initialization)

- `zDir::Vector{Int}`: Interpretation of event indictors:
   zDir[i] = 0: Root is reported for both crossing directions,
           = 1: Root is reported when crossing from negative to positive direction
           = -1: Root is reported when crossing from positive to negative direction

- `x_fixed::Vector{Bool}`: = true, `x_start[i]` is fixed during initialization.
  = false, `x_start[i]` might be changed, e.g., due to an initial impulse.

- `x_nominal::Vector{Float64}`: Nominal values of `x`. if `x_nominal[i]=NaN` no nominal value
  is defined for `x[i]` and a nominal value is computed
  (`x_nominal[i] = abs(x_start[i]) > 1e-7 ? abs(x_start[i]) : 1.0`).

- `x_errorControl::Vector{Bool}`: = true, the absolute error tolerance is set to
  `0.1 * relativeTolerance * x_nominal[i]`. = false, the absolute error tolerance is
  switched off (is set to a large value). This is recommended for variables that are
  basically not limited (for example the angle of a shaft that is permantently rotating
  in the same direction and therefore becomes larger and larger).

- `hev::Float64`: Stepsize used during initialization and at event restart
  if `structureOfDAE = ModiaMath.DAE_NoSpecialStructure`.
  Otherwise `hev` is ignored.

- `scaleConstraintsAtEvents::Bool`: Dummy argument. Only kept for backwards compatibility.

- `nc::Int`: Dummy argument. Only kept for backwards compatibility.

- `jac`: Sparse Jacobian datastructure (currently not supported).

- `maxSparsity::Float64`: A sparse Jacobian is only used during simulation if
  sparseness of jac < maxSparsity (currently not supported)

- `getResultNames::Function`: Function that returns the names of the variables to be
  stored in the result data structure.

- `storeResult!::Function`: Function that stores the raw results.

- `getResult::Function`: Function that resturns the result data structure after the simulation.

- `defaultTolerance::Float64`: Model specific default relative tolerance, if not redefined in the call to
  [`ModiaMath.simulate!`](@ref).

- `defaultStartTime::Float64`: Model specific default start time in [s], if not redefined in the call to
  [`ModiaMath.simulate!`](@ref).

- `defaultStopTime::Float64`: Model specific default stop time in [s], if not redefined in the call to
  [`ModiaMath.simulate!`](@ref).

- `defaultInterval::Float64`: Model specific default interval in [s], if not redefined in the call to
  [`ModiaMath.simulate!`](@ref). Result data is stored every `defaultInterval` seconds.
  If `defaultInterval=NaN`, the default interval is computed as
  `interval = (stopTime - startTime)/500.0`.

- `defaultMaxStepSize::Float64`: Model specific default of the maximum absolute value of the step size.
  If `defaultMaxStepSize=NaN`, the default maximum step size of the integrator is used.

- `defaultMaxNumberOfSteps::Union{Int,Missing}`: Model specific default of the maximum number of
  steps to be taken by the solver in its attempt to reach the next output time.
  If `defaultMaxNumberOfSteps=missing`, the default value of the integrator is used (= 500).
"""
mutable struct SimulationState
    name::Symbol                      # Name of model
    model::Any                        # Model

    getModelResidues!::Function
    getVariableName::Function
    getResultNames::Function
    storeResult!::Function
    getResult::Function
    eventHandler::EventHandler
    structureOfDAE::StructureOfDAE
    is_constraint::Vector{Bool}
    has_constraintDerivatives::Bool

    nx::Int   # Length of x-vector
    nd::Int   # fd(der(x),x,t) = 0; nd = nx - nc
    nc::Int   #        fc(x,t) = 0; nc = count(is_constraint)
    nw::Int   # Number of auxiliary variables

    # Other model information
    nz::Int                           # Number of event indicators
    zDir::Vector{Cint}                # zDir[i] =  0: Root is reported for both crossing directions
                                     #         =  1: Root is reported when crossing from negative to positive direction
                                     #         = -1: Root is reported when crossing from positive to negative direction
    sparse::Bool                      # = true, if (sparse) jac present and shall be used
    # jac::Union{SparseMatrixCSC{Float64,Cint},Void}
                                     # Optional sparse Jacobian of DAE: der(f!,y) + cr*der(f!, yp)
                                     # (cr is a constant provided by the integrator)
    jac::Nothing
    #cg::Union{ModiaMath.SparseJacobian.ColumnGroups,Void} # if sparse, column groups of Jacobian jac
    cg::Nothing # if sparse, column groups of Jacobian jac

    # Model specific information not used by ModiaMath (can be used as default setting for ModiaMath)
    defaultTolerance::Float64    # default relative integration tolerance
    defaultStartTime::Float64    # default start time
    defaultStopTime::Float64     # default stop time
    defaultInterval::Float64     # default integration interval
    defaultMaxStepSize::Float64
    defaultMaxNumberOfSteps::Union{Int,Missing}

    # Run-time information
    time::ModiaMath.Time                        # Actual simulation time
    logger::ModiaMath.Logger                    # Logger object
    statistics::ModiaMath.SimulationStatistics  # Statistics object (complete after a full simulation run)
    tolerance::Float64                          # Relative integration tolerance
    startTime::Float64                          # Start time of the simulation
    stopTime::Float64                           # Stop time of the simulation
    interval::Float64                           # Interval of the simulation
    maxStepSize::Float64
    maxNumberOfSteps::Union{Int,Missing}

    x_start::Vector{Float64}
    x_fixed::Vector{Bool}
    x_nominal::Vector{Float64}
    x_errorControl::Vector{Bool}

    w::Vector{Float64}  # length(w) = nw

    # Auxiliary storage needed during initialization and at events
    eqInfo::ModiaMath.NonlinearEquationsInfo
    xev_old::Vector{Float64}
    xev_beforeEvent::Vector{Float64}
    xev::Vector{Float64}
    yNonlinearSolver::Vector{Float64}
    derxev::Vector{Float64}
    residues::Vector{Float64}
    residues2::Vector{Float64}
    derx_zero::Vector{Float64}
    jac2::Union{Matrix{Float64}, Nothing}
    yScale::Vector{Float64}
    rScale::Vector{Float64}
    tev::Float64
    hev::Float64
    FTOL::Float64
    scaleConstraintsAtEvents::Bool
    initialization::Bool

    # Information updated during simulation
    compute_der_fc::Bool          # = true, if derivative of constraint equations has to be returned in the residue
    storeResult::Bool             # = true, if DAE is called to store results of variables
    rawResult::ModiaMath.RawResult
    result::Any

    function SimulationState(name::Union{AbstractString,Symbol},
                             getModelResidues!::Function,
                             x_start::Vector{Float64},
                             getVariableName::Function=defaultVariableName;
                             structureOfDAE::StructureOfDAE = DAE_LinearDerivativesAndConstraints,
                             is_constraint::Vector{Bool} = fill(false, length(x_start)),
                             has_constraintDerivatives::Bool = false,
                             nc=0,
                             nz=0,
                             nw=0,
                             zDir::Vector{Int}=fill(0, nz),
                             x_fixed::Vector{Bool}=fill(false, length(x_start)),
                             x_nominal::Vector{Float64}=fill(NaN, length(x_start)),
                             x_errorControl::Vector{Bool}=fill(true, length(x_start)),
                             hev=1e-8,
                             scaleConstraintsAtEvents::Bool=true,
                             jac=nothing,
                             maxSparsity::Float64=0.1,
                             getResultNames::Function=ModiaMath.getResultNames,
                             storeResult!::Function=ModiaMath.storeRawResult!,
                             getResult::Function=ModiaMath.getStringDictResult,
                             defaultTolerance=1e-4,
                             defaultStartTime=0.0,
                             defaultStopTime=1.0,
                             defaultInterval=NaN,
                             defaultMaxStepSize=NaN,
                             defaultMaxNumberOfSteps=missing)

        # Check input arguments
        @assert(length(is_constraint) == length(x_start))
        @assert(nz >= 0)
        @assert(nw >= 0)
        @assert(length(x_fixed) == length(x_start))
        @assert(length(zDir) == nz)
        @assert(hev > 0.0)
        @assert(length(x_errorControl) == length(x_start))
        nx = length(x_start)

        if typeof(jac) != Nothing
            @assert(size(jac, 1) == nx)
            @assert(size(jac, 2) == nx)
        end

        @assert(0.0 <= maxSparsity <= 1.0)
        @assert(defaultTolerance > 0.0)
        @assert(defaultStopTime >= defaultStartTime)
        @assert(isnan(defaultInterval) || defaultInterval > 0.0)

        # Compute nominal values
        x_nominal2 = ones(nx)
        for i in 1:nx
            if isnan( x_nominal[i] )
                if abs(x_start[i]) > 1e-7
                    x_nominal2[i] = abs(x_start[i])
                end
            else
                x_nominal2[i] = x_nominal[i]
            end
        end

        yScale = ones(nx)
        for i = 1:nx
            yScale[i] = 1/x_nominal2[i]
        end

        # Compute utility elements
        nc2 = count(is_constraint)
        nd  = nx - nc2
        eventHandler = EventHandler(nz=nz)

        if structureOfDAE == DAE_LinearDerivativesAndConstraints
            eqInfo = ModiaMath.NonlinearEquationsInfo("???", nx, getEventResidues_DAE_LinearDerivativesAndConstraints!)
            jac2   = zeros(nx,nx)
            nc3    = nc2
        else
            eqInfo = ModiaMath.NonlinearEquationsInfo("???", nx, getEventResidues_DAE_NoSpecialStructure!)
            jac2   = nothing
            nc3    = missing
        end

        # If jac is provided and sparse enough, use sparse methods
        #sparse = typeof(jac) == Void ? false : nnz(jac)/(nx*nx) < maxSparsity
        #jac    = sparse ? convert(SparseMatrixCSC{Float64,Cint}, jac) : nothing
        #cg     = sparse ? ColumnGroups(jac) : nothing
        #if sparse
        #   @assert( nnz(jac) >= nx )   # otherwise jac is singular
        #end
        sparse = false
        jac = nothing
        cg = nothing

        new(Symbol(name), nothing, getModelResidues!, getVariableName, getResultNames,
            storeResult!, getResult, eventHandler, structureOfDAE, is_constraint,
            has_constraintDerivatives, nx, nd, nc2, nw, nz, zDir,
            sparse, jac, cg, defaultTolerance, defaultStartTime, defaultStopTime,
            defaultInterval, defaultMaxStepSize, defaultMaxNumberOfSteps, NaN, ModiaMath.Logger(),
            ModiaMath.SimulationStatistics(structureOfDAE, nx, nc3, sparse, sparse ? cg.ngroups : 0),
            NaN, NaN, NaN, NaN, NaN, missing,
            x_start, x_fixed, x_nominal2, x_errorControl, zeros(nw),
            eqInfo, zeros(nx), zeros(nx), zeros(nx),
            zeros(nx), zeros(nx), zeros(nx), zeros(nx), zeros(nx), jac2,
            yScale, ones(nx), 0.0, hev, 0.0,
            scaleConstraintsAtEvents, false, false, false)
    end
end

"""
    result = isNearlyEqual(x1, x2, x_nominal, tolerance)

The function returns true if x1 and x2 are nearly equal
"""
isNearlyEqual(x1::Float64, x2::Float64, x_nominal::Float64, tolerance::Float64) = abs(x1 - x2) <= max(tolerance*max(abs(x1),abs(x2)), 0.1*x_nominal*tolerance)

function getEventResidues_DAE_NoSpecialStructure!(eqInfo::ModiaMath.NonlinearEquationsInfo, y::Vector{Float64}, r::Vector{Float64})

    model = eqInfo.extraInfo
    sim::SimulationState = model.simulationState
    @assert(length(y) == eqInfo.ny)
    @assert(length(r) == eqInfo.ny)
    @assert(length(y) == sim.nx)

    residues = sim.residues
    hev      = sim.hev

    # Copy unknowns (y) to xev and derxev
    #if sim.nc == 0
    #    # Unknowns are at the left limit
    #    for i in eachindex(y)
    #        sim.xev[i]    = sim.xev_old[i]
    #        sim.derxev[i] = (sim.xev_old[i] - y[i]) / hev
    #    end
    #else # sim.nc > 0
        for i in eachindex(y)
            if sim.initialization && sim.x_fixed[i]
                # Unknowns are at the left limit
                sim.xev[i]    = sim.xev_old[i]
                sim.derxev[i] = (sim.xev_old[i] - y[i]) / hev
            else
                # Unknowns are at the right limit
                sim.xev[i]    = y[i]
                sim.derxev[i] = (y[i] - sim.xev_old[i]) / hev
            end
        end
    #end

    # Compute residues
    sim.time = sim.tev
    Base.invokelatest(sim.getModelResidues!, model, sim.tev, sim.xev, sim.derxev, residues, sim.w)

   # Copy to eq-residues and scale with h
   #=
    if sim.scaleConstraintsAtEvents
        for i = 1:sim.nd
            r[i] = residues[i]
        end

        for i = sim.nd + 1:sim.nx
            r[i] = residues[i] / hev
        end
    else
    =#
        for i = 1:sim.nx
            r[i] = residues[i]
        end
    #end
    eqInfo.lastNorm_r = norm(r, Inf)
    eqInfo.lastrScaledNorm_r = norm(sim.rScale .* r, Inf)

    # println("maxabs(r) = ", maxabs(r))
    # println("residue = "),display(r)
    # println("hev = ", hev, ", x_old = ", sim.xev_old[1], ", x = ", y[1], ", der(x) = ",sim.derxev[1], ", r = ", r[1])

    return nothing
end

function getEventResidues_DAE_LinearDerivativesAndConstraints!(eqInfo::ModiaMath.NonlinearEquationsInfo, y::Vector{Float64}, r::Vector{Float64})

    model = eqInfo.extraInfo
    sim::SimulationState = model.simulationState
    @assert(length(y) == eqInfo.ny)
    @assert(length(r) == eqInfo.ny)
    @assert(length(y) == sim.nx)

    residues  = sim.residues
    residues2 = sim.residues2

    # Copy unknowns (y) to xev and derxev
    if sim.nc == 0
        # Unknowns are at the left limit
        for i in eachindex(y)
            sim.xev[i]    = sim.xev_old[i]
            sim.derxev[i] = sim.xev_old[i] - y[i]
        end
    else # sim.nc > 0
        for i in eachindex(y)
            if sim.initialization && sim.x_fixed[i]
                # Unknowns are at the left limit
                sim.xev[i]    = sim.xev_old[i]
                sim.derxev[i] = sim.xev_old[i] - y[i]
            else
                # Unknowns are at the right limit
                sim.xev[i]    = y[i]
                sim.derxev[i] = y[i] - sim.xev_old[i]
            end
        end
    end

    # Compute residues
    sim.time = sim.tev
    Base.invokelatest(sim.getModelResidues!, model, sim.tev, sim.xev, sim.derx_zero, residues , sim.w)
    Base.invokelatest(sim.getModelResidues!, model, sim.tev, sim.xev, sim.derxev   , residues2, sim.w)

    for i = 1:sim.nx
        if sim.is_constraint[i]
            r[i] = residues[i]
        else
            r[i] = residues[i] - residues2[i]
        end
    end

    eqInfo.lastNorm_r        = norm(r, Inf)
    eqInfo.lastrScaledNorm_r = norm(sim.rScale .* r, Inf)

    # println("maxabs(r) = ", maxabs(r))
    # println("residue = "),display(r)

    return nothing
end

const small = sqrt( eps(Float64) )

function reinitialize!(model, sim::SimulationState, tev::Float64)
    nx = sim.nx
    nd = sim.nd
    nc = sim.nc

    # Determine consistent xev, der(xev)
    if ModiaMath.isLogEvents(sim)
        if sim.structureOfDAE == DAE_NoSpecialStructure
            println("        determine consistent DAE variables x,der(x) (with implicit Euler step; step size = ", sim.hev, ")")
        elseif sim.structureOfDAE == DAE_LinearDerivativesAndConstraints
            if nc == 0
                println("        for given x, determine consistent DAE variables der(x) (solving a linear equation system)")
            else
                println("        determine consistent DAE variables x,der(x) (with analytical integral over time instant)")
            end
        elseif sim.structureOfDAE == DAE_ExplicitDerivatives
            println("        for given x, compute der(x)")
        else
            error("... Cannot occur. Variable structureOfDAE = ", sim.structureOfDAE, " is wrong.")
        end
    end

    # Compute consistent xev(t+)
    if  sim.structureOfDAE == DAE_NoSpecialStructure ||
       (sim.structureOfDAE == DAE_LinearDerivativesAndConstraints && nc > 0)

        for i in eachindex(sim.xev)
            sim.xev_old[i] = sim.xev[i]
        end

        sim.tev = tev
        sim.eqInfo.lastNorm_r        = 1.0
        sim.eqInfo.lastrScaledNorm_r = 1.0

        for i in eachindex(sim.yNonlinearSolver)
            sim.yNonlinearSolver[i] = sim.xev[i]
        end

        ModiaMath.solveNonlinearEquations!(sim.eqInfo, sim.yNonlinearSolver; yScale=sim.yScale, rScale=sim.rScale, FTOL=sim.FTOL)

        # Print log message
        if ModiaMath.isLogInfos(sim)
            for i = 1:nx
                if !isNearlyEqual(sim.xev_beforeEvent[i], sim.xev[i], sim.x_nominal[i], sim.tolerance)
                    xname = sim.getVariableName(model, Category_X, i)
                    println("            ", xname, " = ", sim.xev_beforeEvent[i], " changed to ", sim.xev[i])
                end
            end
        end
    end

    # Compute consistent derivatives der(xev)(t+)
    if sim.structureOfDAE == DAE_ExplicitDerivatives
        Base.invokelatest(sim.getModelResidues!, model, sim.tev, sim.xev, sim.derx_zero, sim.derxev, sim.w)

    elseif sim.structureOfDAE == DAE_LinearDerivativesAndConstraints
        # Determine linear factors
        jac2 = sim.jac2
        sim.compute_der_fc = true
        if ModiaMath.isLogEvents(sim) && nc > 0 && sim.has_constraintDerivatives
            println("        compute der(x) with Jacobian that is constructed with model provided constraint derivatives (der(fc))")
        end
        Base.invokelatest(sim.getModelResidues!, model, sim.tev, sim.xev, sim.derx_zero, sim.residues, sim.w)

        # Determine der(fd, der(x)) and der(fc, x) using the fact that fd and der(fc) are linear in der(x)
        for j = 1:nx
            sim.derx_zero[j] = 1.0
            Base.invokelatest(sim.getModelResidues!, model, sim.tev, sim.xev, sim.derx_zero, sim.residues2, sim.w)
            for i = 1:nx
                if !sim.is_constraint[i] || sim.has_constraintDerivatives
                    # println("... 1: jac2[$i,$j]")
                    jac2[i,j] = sim.residues2[i] - sim.residues[i]
                end
            end
            sim.derx_zero[j] = 0.0
        end
        sim.compute_der_fc = false

        # Determine der(fc, x) with forward difference, if not yet computed
        if !sim.has_constraintDerivatives && nc > 0
        	inc::Float64 = 0.0
            xinc = sim.yNonlinearSolver
            for j = 1:nx
                xinc[j] = 0.0
            end
            for j = 1:nx
                xj  = sim.xev[j]
                inc = small*max( abs(xj), sim.x_nominal[j] )
                inc = (xj + inc) - xj
                sim.xev[j] += inc
                Base.invokelatest(sim.getModelResidues!, model, sim.tev, sim.xev, sim.derx_zero, sim.residues2, sim.w)
                for i = 1:nx
                    if sim.is_constraint[i]
                        # println("... 2: jac2[$i,$j]")
                        jac2[i,j] = sim.residues2[i]/inc - sim.residues[i]/inc
                    end
                end
                sim.xev[j] = xj
            end
        end

        # Solve linear system of equation: jac2*derxev = -residue
        for j = 1:nx
            sim.derxev[j] = -sim.residues[j]
        end
        ldiv!(lu!(jac2),sim.derxev)
    end
end

function eventIteration!(model, sim::SimulationState, tev::Float64)
    eh::EventHandler = sim.eventHandler

    # Initialize event iteration
    initEventIteration!(eh, tev)

    # Perform event iteration
	iter_max = 20
	success  = false
    for iter = 1:iter_max
        # Determine event branches
        for i = 1:sim.nx
            sim.xev_beforeEvent[i] = sim.xev[i]
        end
        eh.event = true
        Base.invokelatest(sim.getModelResidues!, model, tev, sim.xev, sim.derxev, sim.residues, sim.w)

        # Fix event branches and determine new consistent sim.derxev_start (and sim.xev_start during initialization)
        eh.event   = false
        eh.initial = false
        reinitialize!(model, sim, tev)

		if terminateEventIteration!(eh)
            eh.event = false
			success  = true
            break
        end
    end

	if !success
		error("Maximum number of event iterations (= $iter_max) reached")
    end
end

function initialize!(model, sim::SimulationState, t0::Float64, nt::Int, tolerance::Float64)
    eh::EventHandler = sim.eventHandler
    eh.zDir    = sim.zDir
    eh.logger  = sim.logger
    eh.initial = true
    eh.afterSimulationStart = false
    sim.eqInfo.extraInfo = model
    sim.rawResult = ModiaMath.RawResult(nt, sim.getResultNames(model))
    sim.model     = model
    sim.FTOL      = sim.tolerance    # max(sim.tolerance, eps(Float64)^(1 / 3))   # eps(Float64)^(1 / 3)   # residue tolerance for nonlinear solver
    # println("... FTOL = ", sim.FTOL)
    sim.initialization = true

    # Initialize auxiliary arrays for event iteration
    # println("... initialize!: scaleConstraintsAtEvents = ", sim.scaleConstraintsAtEvents, ", nc = ", sim.nc, ", nd = ", sim.nd)
    for i in eachindex(sim.xev),
        sim.xev[i]    = sim.x_start[i]
        sim.derxev[i] = 0.0
    end

    #if sim.scaleConstraintsAtEvents
    #    for i = 1:sim.nd
    #        sim.rScale[i] = 1.0
    #    end
    #    for i = sim.nd + 1:sim.nx
    #        sim.rScale[i] = 1.0 / sim.hev
    #    end
    #else
        for i = 1:sim.nx
            sim.rScale[i] = 1.0
        end
    #end

    # Print initial values
    if ModiaMath.isLogEvents(sim)
        println("        initial values:")
        x_table = DataFrames.DataFrame(name=String[sim.getVariableName(model, Category_X, i)  for i = 1:sim.nx],
                                       start=sim.x_start, fixed=sim.x_fixed, nominal=sim.x_nominal)

        # Print x_table, but indented and type-information removed from the heading
        io = IOBuffer()
        show(io, x_table, summary=false, rowlabel=:x)
        str=String(take!(io))
        newline = isequal('\n')
        i1=findfirst(newline, str)
        i2=findnext(newline, str, i1+1)
        i3=findnext(newline, str, i2+1)
        str2="          " * str[i1+1:i2] * str[i3+1:end]
        str3=replace(str2, "\n" => "\n          ")
        println(str3, "\n")
    end

    # Perform initial event iteration
    eventIteration!(model, sim, t0)
    sim.initialization = false

    # Compute auxiliary variables w and store all results
    computeAndStoreResult!(model, sim, t0, sim.xev, sim.derxev)
    eh.afterSimulationStart = true

    return InitInfo(sim.xev, sim.derxev;
                    y_nominal=sim.x_nominal,
                    y_errorControl=sim.x_errorControl,
                    maxTime=eh.maxTime,
                    nextEventTime=eh.nextEventTime,
                    integrateToEvent=eh.integrateToEvent,
                    terminate=eh.restart == Terminate)
end

getResidues!(model, sim::SimulationState, t::Float64, y::Vector{Float64}, yp::Vector{Float64}, residues::Vector{Float64}, hcur) =
   Base.invokelatest(sim.getModelResidues!, model, t, y, yp, residues, sim.w)

function computeAndStoreResult!(model, sim::SimulationState, t::Float64, y::Vector{Float64}, yp::Vector{Float64})
    sim.storeResult = true
    sim.time        = t
    Base.invokelatest(sim.getModelResidues!, model, t, y, yp, sim.residues, sim.w)
    sim.storeResult = false
    sim.storeResult!(sim.rawResult, model, t, y, yp, sim.w)
    return nothing
end

function terminate!(model, sim::SimulationState, t::Float64, y::Vector{Float64}, yp::Vector{Float64})
    eh::EventHandler = sim.eventHandler
    eh.terminal = true
    sim.time    = t
    Base.invokelatest(sim.getModelResidues!, model, t, y, yp, sim.residues, sim.w)
    eh.terminal = false
    #sim.result,nt = Result.getDictResult(sim.rawResult)
    sim.result, nt  = sim.getResult(model, sim.rawResult)
    ModiaMath.set_nResultsForSimulationStatistics!(sim.statistics, nt)
    return sim.result
end

function processEvent!(model, sim::SimulationState, t::Float64, y::Vector{Float64}, yp::Vector{Float64}, eventInfo::EventInfo)
    eh::EventHandler = sim.eventHandler

    # Event iteration
    for i in eachindex(sim.xev)
        sim.xev[i]    = y[i]
        sim.derxev[i] = yp[i]
    end
    eventIteration!(model, sim, t)
    for i in eachindex(sim.xev)
        y[i]  = sim.xev[i]
        yp[i] = sim.derxev[i]
    end

    # Compute auxiliary variables w and store all results
    computeAndStoreResult!(model, sim, t, y, yp)

    eventInfo.restart          = eh.restart
    eventInfo.maxTime          = eh.maxTime
    eventInfo.nextEventTime    = eh.nextEventTime
    eventInfo.integrateToEvent = eh.integrateToEvent

    return nothing
end

function getEventIndicators!(model, sim::SimulationState, t::Float64, y::Vector{Float64}, yp::Vector{Float64}, z::Vector{Float64})
    eh::EventHandler = sim.eventHandler
    eh.crossing = true
    sim.time    = t
    Base.invokelatest(sim.getModelResidues!, model, t, y, yp, sim.residues, sim.w)
    eh.crossing = false

    for i = 1:eh.nz
        z[i] = eh.z[i]
    end
    return nothing
end

#=
"""
    computeFullJacobian!(model, sim, t, y, yp, residues, cj, ewt)

Compute full Jacobian

```julia
jac = der(f,y) + cj*der(f,yp).
```

numerically by finite differences
"""
function computeJacobian!(model, sim::SimulationState, t::Float64, y::Vector{Float64}, yp::Vector{Float64}, r::Vector{Float64},
                          fulljac::Matrix{Float64}, hcur::Float64, cj::Float64, ewt::Vector{Float64})
    # Compute one column of the Jacobian once at a time
    inc::Float64 = 0.0
    yj::Float64  = 0.0
    ypj::Float64 = 0.0
    for j = 1:sim.nx
        # Determine increment for column j
        yj  = y[j]
        ypj = yp[j]
        inc = max( small*max( abs(yj), abs(hcur*ypj) ), 1.0/ewt[j] )
        if hcur*ypj < 0.0
            inc = -inc
        end
        inc = (yj + inc) - yj
        y[j]  += inc
        yp[j] += cj*inc

        # Compute residues
        Base.invokelatest(sim.getModelResidues!, model, t, y, yp, sim.residues, sim.w)

        # Compute finite differences and store them in the j column of the Jacobian
        for i = 1:sim.nx
            fulljac[i,j] = sim.residues[i]/inc - r[j]/inc
        end

        # Reset y and yp
        y[j]  = yj
        yp[j] = ypj
    end
end
=#


text:
module ParallelTest

	function hosts()
		hosts = []
		pids = []
		for i in workers()
			host, pid = fetch(@spawnat i (gethostname(), getpid()))
			push!(hosts, host)
			push!(pids, pid)
		end
		println("your hosts = $hosts")
		println("your pids = $pids")
	end

	function addprocs_ec2(machines::Array{AbstractString})
		pids = addprocs(machinces,sshflags=`-o UserKnownHostsFile=/dev/null -o StrictHostKeyChecking=no -o LogLevel=ERROR`,tunnel=true)
		return pids
	end

	# to be run on all nodes

	function sayhello()
	   println("hi I am worker number $(myid()), I live on $(gethostname())")
	end

	function domath(i::Integer)
		if i < 2
			error("need i>2")
		end
		println("Hi, I am worker number $(myid()) doing some math")
		x = rand(i,i)
		mean(x)
	end

	serial() = map( n -> sum(svd(rand(n,n))[1]) , [800 for i in 1:64])
	parallel() = pmap( n -> sum(svd(rand(n,n))[1]) , [800 for i in 1:64])

	function doBIGmath(n=30000,m=10000)
		println("allocating a $n by $m matrix ")
		x = rand(n,m);
		r = round(sizeof(x) /  1.074e+9 , 2)
		println("my matrix has $r GB")
		mean(x)
	end

	function machines(; verbose = false)
		if verbose sayhello() end
		# which scheduler system
		if haskey(ENV,"PBS_NODEFILE")
			sys = :pbs
		    # read PBS_NODEFILE
		    filestream = open(ENV["PBS_NODEFILE"])
		    seekstart(filestream)
		    node_file = readlines(filestream)

		    # strip eol
		    node_file = map(x->strip(x,['\n']),node_file)

		    # get number of workers on each node
		    procs = Dict{ASCIIString,Int}()
		    for n in node_file
		        procs[n] = get(procs,n,0) + 1
		    end

		    println("name of compute nodes and number of workers:")
		    println(procs)

		    # add processes on master itself
		    master = ENV["HOSTNAME"]

		    wrker = 0
		    while wrker < ppn
		        addprocs(1)
		        wrker += 1
		    end
		    println("added $wrker processes on master itself")

		    # add procs on other machines
		    for (k,v) in procs
		        wrker = 0
		        if k!=master
		            while wrker < ppn
		                addprocs([k], dir= JULIA_HOME)
		                # println("addprocs($k)")
		                wrker += 1
		            end
		            println("added $wrker processes on machine $k")
		        end
		    end
		elseif haskey(ENV,"SLURM_JOB_NODELIST")
			# magi[90-91]   salloc -N 2
			# magi[50,90-91]  salloc -N 3 -n 15

			# salloc -N 3 --ntasks-per-node=5
			# magi3 /home/dist/florian.oswald $ scontrol show hostnames $SLURM_JOB_NODELIST
			# magi50
			# magi90
			# magi91
			machine_file = readlines(`scontrol show hostnames`)
			if (gethostname() != machine_file[1])
				# println("shutting down this process because not master")
				exit()
			end
			println("ParallelTest :: machine file $machine_file")
			pids = readlines(`pgrep julia`)  # get process ids on surviving machine

			println("ParallelTest :: pids = $pids")

			# designate first process id as master who will call addprocs on the others.
			# so, if you are not master, exit
			if getpid() != parse(Int,pids[1])
				println("I am not the master, so goodbye")
				exit()
			end

			sleep(5)
			if haskey(ENV,"SLURM_NTASKS")
				tasks_per_cpu = floor(Int,parse(Int,ENV["SLURM_NTASKS"]) / length(machine_file))
			else
				tasks_per_cpu = length(machine_file)
			end
			println("ParallelTest :: tasks per cpu: $tasks_per_cpu")

			machines = [(i,tasks_per_cpu) for i in machine_file]
			println("ParallelTest :: array of machines: $machines")
			addprocs(machines)
		        println("ParallelTest :: done. added $(length(workers()))")
		        println("ParallelTest :: workers: $(workers())")
			return(workers())
		elseif haskey(ENV,"PE_HOSTFILE")
			sys = :sge
		else
			sys = :local
		end
	end

end # module


text:
using LazySets.Arrays: extend,
                       _vector_type, _matrix_type

for _dummy_ in 1:1 # avoid global variable warnings
    # reseeding with random seed
    rng = LazySets.GLOBAL_RNG
    seed = rand(1:10000)
    LazySets.reseed(rng, seed)
    n1 = rand(Int)
    LazySets.reseed(rng, seed)
    n2 = rand(Int)
    @test n1 == n2

    # StrictlyIncreasingIndices
    vectors = Vector{AbstractVector{Int}}()
    for v in LazySets.StrictlyIncreasingIndices(5, 4)
        push!(vectors, copy(v))
    end
    @test vectors == [[1, 2, 3, 4], [1, 2, 3, 5], [1, 2, 4, 5], [1, 3, 4, 5], [2, 3, 4, 5]]

    # square matrix
    @test LazySets.issquare([2 3; 0 0])
    @test LazySets.issquare(sparse([1], [1], [1], 3, 3))
    @test !LazySets.issquare(hcat([2 3]))
    @test !LazySets.issquare(sparse([1], [1], [1], 3, 4))

    # invertible matrix: (1) invertible, (2) non-invertible, (3) non-square
    # dense matrix
    @test LazySets.isinvertible([2 3; 1 2])
    @test !LazySets.isinvertible([2 3; 0 0])
    @test !LazySets.isinvertible(hcat([2 3]))
    # sparse matrix
    @test LazySets.isinvertible(sparse([2 3; 1 2]))
    @test !LazySets.isinvertible(sparse([2 3; 0 0]))
    @test !LazySets.isinvertible(sparse(hcat([2 3])))
    # diagonal matrix
    @test LazySets.isinvertible(Diagonal([2 0; 0 2]))
    @test !LazySets.isinvertible(Diagonal([2 0; 0 0]))

    # matrix rank
    A = sprandn(2, 10, 0.4)
    @test rank(Matrix(A)) == rank(A) == rank(view(A, :, :))

    for N in [Float64, Rational{Int}, Float32]
        # substitution
        x = N[1, 2, 3]
        substitution = Dict(1 => N(4), 3 => N(0))
        @test LazySets.substitute(substitution, x) == N[4, 2, 0]
        LazySets.substitute!(substitution, x)
        @test x == N[4, 2, 0]

        A = N[1 4; 2 5; 3 6]
        x1 = N[0, 2, 0]
        y1 = N[3, 0]
        x2 = SingleEntryVector(2, 3, N(2))
        y2 = SingleEntryVector(1, 2, N(3))
        @test -x2 == SingleEntryVector(2, 3, N(-2))
        @test inner(x1, A, y1) == dot(x1, A * y1) == inner(x2, A, y2) ==
              dot(x2, A * y2) == N(12)

        x = N[0, 1, -1]
        @test LazySets.rectify(x) == N[0, 1, 0]

        # approximate permutation check
        v1 = [N[1, 2], N[2, 3], N[3, 4]]
        tol = Base.rtoldefault(N)/2
        v2 = [N[3 + tol, 4 - tol], N[2 - tol, 3 + tol], N[1 + tol, 2 + tol]]
        v3 = [N[2, 3], N[1, 2], N[4, 4]]
        @test ispermutation(v1, v2) && !ispermutation(v2, v3)
        @test ispermutation(N[1], [N(1) + Base.rtoldefault(N)])

        # to_negative_vector
        for v in [N[-1, 0, 1], sparsevec([1, 3], N[-1, 1], 3), N(-1):N(1):N(1)]
            u = LazySets.Arrays.to_negative_vector(v)
            @test u isa Vector{N} && u == N[1, 0, -1]
        end

        # ============================================
        # Corresponding vector types and matrix types
        # ============================================

        # sparse
        vec = sparsevec([1, 3], N[1, 3], 3)
        mat = sparse([1, 3], [1, 3], N[1, 3], 3, 3)
        @test _vector_type(typeof(vec)) == SparseVector{N, Int}
        @test _matrix_type(typeof(vec)) == SparseMatrixCSC{N,Int64}
        @test _vector_type(typeof(mat)) == SparseVector{N, Int}
        @test _matrix_type(typeof(mat)) == SparseMatrixCSC{N,Int64}

        # regular
        vec = N[1, 0, 3]
        mat = N[1 0 0; 0 0 0; 0 0 3]
        @assert _vector_type(typeof(vec)) == Vector{N}
        @assert _matrix_type(typeof(vec)) == Matrix{N}
        @assert _vector_type(typeof(mat)) == Vector{N}
        @assert _matrix_type(typeof(mat)) == Matrix{N}

        # other: Diagonal
        mat = Diagonal(N[1, 2])
        @test _vector_type(typeof(mat)) == Vector{N}
        @assert _matrix_type(typeof(mat)) == Diagonal{N,Array{N,1}}
    end

    for N in [Float64, Float32]
        # modified dot product
        @test isnan(dot(N[1, 0], N[Inf, -Inf]))
        @test LazySets.dot_zero(N[1, 0], N[Inf, -Inf]) == N(Inf)

        # remove zero columns
        m = 3
        n = 10
        A = ones(N, m, n)
        A[:, 1] = A[:, 5] = A[:, n] = zeros(N, m)
        nzcol = LazySets.Arrays.nonzero_columns(A)
        B = A[:, nzcol]
        @test size(B) == (m, n-3)
        Asp = sparse(A)
        @test LazySets.Arrays.nonzero_columns(Asp) == [2, 3, 4, 6, 7, 8, 9]

        # extend by orthogonal complement
        M = N[1 1; 2 2; 3 4.]
        @assert rank(M) == 2
        Mext, inv_Mext = extend(M)
        @test rank(Mext) == 3
        @test Mext * inv_Mext ≈ Matrix(one(N)I, 3, 3)
        Md = N[1 1; 2 2; 3 3.]
        @test_throws ArgumentError extend(Md) # test default argument check
        @test_throws ArgumentError extend(Md, check_rank=true)
    end
end


text:
function MOI.eval_objective(p::Problem{T}, z::Vector{T}) where T
    trajectory!(p.trajopt.x, p.trajopt.u, z,
        p.trajopt.model.idx.x, p.trajopt.model.idx.u)
    eval_obj(p.trajopt.obj, p.trajopt.x, p.trajopt.u, p.trajopt.w)
end

function MOI.eval_objective_gradient(p::Problem{T}, grad, z) where T
    fill!(grad, 0.0)
    trajectory!(p.trajopt.x, p.trajopt.u, z,
        p.trajopt.model.idx.x, p.trajopt.model.idx.u)
    eval_obj_grad!(grad, p.trajopt.model.idx.xu, p.trajopt.obj, p.trajopt.x, p.trajopt.u, p.trajopt.w)
end

function MOI.eval_constraint(p::Problem{T}, con, z) where T
    fill!(con, 0.0)
    trajectory!(p.trajopt.x, p.trajopt.u, z,
        p.trajopt.model.idx.x, p.trajopt.model.idx.u)
    eval_con!(con, p.idx.dyn_con, p.trajopt.model.dyn, p.trajopt.x, p.trajopt.u, p.trajopt.w)
    !isempty(p.idx.stage_con) && eval_con!(con, p.idx.stage_con, p.trajopt.con.stage, p.trajopt.x, p.trajopt.u, p.trajopt.w)
end

function MOI.eval_constraint_jacobian(p::Problem{T}, jac, z) where T
    fill!(jac, 0.0)
    trajectory!(p.trajopt.x, p.trajopt.u, z,
        p.trajopt.model.idx.x, p.trajopt.model.idx.u)
    eval_jac!(jac, p.idx.dyn_jac, p.trajopt.model.dyn, p.trajopt.x, p.trajopt.u, p.trajopt.w)
    !isempty(p.idx.stage_jac) && eval_jac!(jac, p.idx.stage_jac, p.trajopt.con.stage, p.trajopt.x, p.trajopt.u, p.trajopt.w)
    return nothing
end

function MOI.eval_hessian_lagrangian(p::MOI.AbstractNLPEvaluator, hess, z, σ, λ)
    fill!(hess, 0.0)
    trajectory!(p.trajopt.x, p.trajopt.u, z,
        p.trajopt.model.idx.x, p.trajopt.model.idx.u)
    duals!(p.trajopt.λ_dyn, p.trajopt.λ_stage, λ, p.idx.dyn_con, p.idx.stage_con)
    eval_obj_hess!(hess, p.idx.obj_hess, p.trajopt.obj, p.trajopt.x, p.trajopt.u,p.trajopt.w, σ)
    eval_hess_lag!(hess, p.idx.dyn_hess, p.trajopt.model.dyn, p.trajopt.x, p.trajopt.u, p.trajopt.w, p.trajopt.λ_dyn)
    eval_hess_lag!(hess, p.idx.stage_hess, p.trajopt.con.stage, p.trajopt.x, p.trajopt.u, p.trajopt.w, p.trajopt.λ_stage)
end

MOI.features_available(p::MOI.AbstractNLPEvaluator) = p.hess_lag ? [:Grad, :Jac, :Hess] : [:Grad, :Jac]
MOI.initialize(p::MOI.AbstractNLPEvaluator, features) = nothing
MOI.jacobian_structure(p::MOI.AbstractNLPEvaluator) = p.sp_jac
MOI.hessian_lagrangian_structure(p::MOI.AbstractNLPEvaluator) = p.sp_hess_lag


text:
abstract type FermiAction{Dim,Dirac,fermion,gauge} end

include("./StaggeredFermiAction.jl")
include("./WilsonFermiAction.jl")
include("./DomainwallFermiAction.jl")
include("./MobiusDomainwallFermiAction.jl")
include("./WilsontypeFermiAction.jl")

function FermiAction(D::Dirac_operator{Dim},parameters_action;covneuralnet = nothing) where {NC,Dim}
    diractype = typeof(D)
    if covneuralnet ==  nothing
        hascovnet = false
    else
        hascovnet = true
    end

    if diractype <: Staggered_Dirac_operator
        return StaggeredFermiAction(D,hascovnet,covneuralnet,parameters_action)
    elseif diractype <: Wilson_Dirac_operators
        return WilsonFermiAction(D,hascovnet,covneuralnet,parameters_action)
    elseif diractype <: Domainwall_Dirac_operator
        return DomainwallFermiAction(D,hascovnet,covneuralnet)
    elseif diractype <: MobiusDomainwall_Dirac_operator
        return MobiusDomainwallFermiAction(D,hascovnet,covneuralnet)
    elseif diractype <: Wilson_GeneralDirac_operator
        return Wilson_GeneralDirac_FermiAction(D,hascovnet,covneuralnet,parameters_action)
    else
        error("Action type $diractype is not supported")
    end

end

function evaluate_FermiAction(fermi_action::FermiAction,U,ϕ::AbstractFermionfields)
    error("evaluate_FermiAction(fermi_action,U,ϕ) is not implemented in type fermi_action:$(typeof(fermi_action)), U:$(typeof(U)), and ϕ:$(typeof(ϕ)),  ")
end

function gauss_sampling_in_action!(η::AbstractFermionfields,U,fermi_action::FermiAction)
    error("gauss_sampling_in_action!(η,fermi_action) is not implemented in type η:$(typeof(η)), fermi_action:$(typeof(fermi_action))")
end

#=
 Conventional case:
det(D)^Nf = det(D^+ D)^{Nf/2}
 = int dphi dphi^* exp[- phi^* (D^+ D)^{-1} phi]
 = int dphi dphi^* exp[- phi^* D^{-1} (D^+)^{-1} phi]

 RHMC case:
det(D)^Nf =
 = int dphi dphi^* exp[- phi^* D^{-Nf} phi]
 = int dphi dphi^* exp[- phi^* D^{-Nf/2} D^{-Nf/2} phi]
=#

function sample_pseudofermions!(ϕ::AbstractFermionfields,U,fermi_action::FermiAction,ξ)
    error("sample_pseudofermions!(ϕ,fermi_action,ξ) is not implemented in type ϕ:$(typeof(ϕ)), fermi_action:$(typeof(fermi_action)), ξ:$(typeof(ξ))")
end

function calc_UdSfdU(fermi_action::FermiAction{Dim,Dirac,fermion,gauge},U::Vector{<: AbstractGaugefields},ϕ::AbstractFermionfields) where {Dim,Dirac,fermion,gauge}
    x = U[1]
    UdSfdU = Array{typeof(x),1}(undef,Dim)
    for μ=1:Dim
        UdSfdU[μ] = similar(x)
    end
    calc_UdSfdU!(UdSfdU,fermi_action,U,ϕ)
end

function calc_UdSfdU!(UdSfdU::Vector{<: AbstractGaugefields},fermi_action::FermiAction,U::Vector{<: AbstractGaugefields},ϕ::AbstractFermionfields)
    error("cald_UdSfdU!(UdSfdU,fermi_action,U) is not implemented in type UdSfdU:$(typeof(UdSfdU)), fermi_action:$(typeof(fermi_action)), U:$(typeof(U)), ϕ:$(typeof(ϕ))")
end


text:
"""
```julia
rm(name::AbstractString)
```

Equivalent to `Base.rm`, but it doesn't error if file `name` is missing.
"""
rm(f::AbstractString) = isfile(f) && Base.rm(f)

"""
```julia
clean(name::AbstractString, filelist)
clean(filelist)
```

Create a task `name` (defaults to "clean") that will delete files given in
`filelist`.
"""
function clean(name::AbstractString, filelist)
    task(name) do
        for f in filelist
            rm(f)
        end
    end
end
clean(filelist) = clean("clean", filelist)


text:
__precompile__()

module CLaP

import MathOptInterface
const MOI = MathOptInterface

export MOI,
    NONLINEAR_CONE, POLYHEDRAL_CONE, MOI_CONES, SCALAR_SETS, SUPPORTED_SETS,
    StandardProblem, build_standard_form,
    extract_implied_integer,
    build_cgcp,
    strengthen,
    lift_and_project

include("standard_form.jl")
include("prepross.jl")
include("cgcp.jl")
include("strengthening.jl")

# Lift-and-project callback
include("lift_and_project.jl")

end  # module

max_stars_repo_name: mtanneau/CLaP

max_issues_repo_name: mtanneau/CLaP

max_forks_repo_name: mtanneau/CLaP


text:
export NFFTAlg

using NFFT
function NFFTAlg(u::AbstractArray, v::AbstractArray; padfac=1, m=10)
    uv = Matrix{eltype(u)}(undef, 2, length(u))
    uv[1,:] .= u
    uv[2,:] .= v
    return ObservedNUFT(NFFTAlg(;padfac, m), uv)
end

function NFFTAlg(ac::ArrayConfiguration; padfac=1, m=10)
    u, v = getuv(ac)
    return NFFTAlg(u, v; padfac, m)
end

function padimage(alg::NFFTAlg, img)
    padfac = alg.padfac
    # if no padding exit now
    (padfac == 1) && return convert(Matrix{Complex{eltype(img)}}, img)

    ny,nx = size(img)
    nnx = nextpow(2, padfac*nx)
    nny = nextpow(2, padfac*ny)
    nsx = nnx÷2-nx÷2
    nsy = nny÷2-ny÷2
    cimg = convert(Matrix{Complex{eltype(img)}}, img)
    return PaddedView(zero(eltype(cimg)), cimg,
                      (1:nnx, 1:nny),
                      (nsx+1:nsx+nx, nsy+1:nsy+ny)
                     )
end

function plan_nuft(alg::ObservedNUFT{<:NFFTAlg}, img, dx, dy)
    uv2 = similar(alg.uv)
    uv2[1,:] .= alg.uv[1,:]*dx
    uv2[2,:] .= alg.uv[2,:]*dy
    plan = plan_nfft(uv2, size(img); precompute=NFFT.POLYNOMIAL)
    return plan
end

function make_phases(alg::ObservedNUFT{<:NFFTAlg}, img)
    dx, dy = pixelsizes(img)
    u = @view alg.uv[1,:]
    v = @view alg.uv[2,:]
    return cispi.((u.*dx .+ v.*dy)).*visibility_point.(Ref(img.pulse), u.*dx, v.*dy).*dx.*dy
end

@inline function create_cache(alg::ObservedNUFT{<:NFFTAlg}, plan, phases, img)
    return NUFTCache(alg, plan, phases, transpose(img))
end


text:
"""
Submodule of `SirenSeq`.
Used for playing `SirenSeq` expressions.
Must have some synthesizer setup and must have the `pmidi` program installed.
"""
module Play

export setDefaultPlayPort, playMidi, stop, makeStopMidi

using SirenSeq
using SirenSeq.defaultBpm

## default ALSA midi port to play to
defaultPlayPort = "130:0"

"""
	setDefaultPlayPort(pp::AbstractString)

Sets the default ALSA midi port to which `playMidi` will play to. \\
For example, if 128:0 is the qsynth midi input port then use, `setDefaultPlayPort(\"130:0\")`.
"""
setDefaultPlayPort(pp::AbstractString) = ( global defaultPlayPort = pp ; )

"""
	playMidi([z]; port=defaultPlayPort, path=\"temp\", tpq=96, bpm=defaultBpm)

Creates a midi file caled \"`path`.mid\" (do not include the \".mid\" part) using audio expression `z`.
It then plays the file on ALSA port `port` (e.g. \"130:0\") using the program pmidi, which must be installed on the system.
`tpq` is the midi ticks per quarter note.  `bpm` is the intial beats per minute (may be changed within `z` at any time).
If `z` is not specified it will look for a file called \"temp.mid\" and play that.
"""
function playMidi(x; port::AbstractString=defaultPlayPort, path::AbstractString="temp", tpq::Int=96, bpm::Int=defaultBpm)
	makeMidi(x;path=path,tpq=tpq,bpm=bpm)
	sleep(0.2)
	playMidi(port=port,path=path)
end

function playMidi(; port::AbstractString=defaultPlayPort, path::AbstractString="temp")
	try
		println("playing: $(path).mid ... \t")
		p = spawn(`pmidi -d 0.5 -p $(port) $(path).mid`)
		if !process_running(p) ; throw(error("failed to run pmidi")) ; end
		sleep(0.5)
		return p
	catch err
		if isa(err,InterruptException)
			print("interrupted!\n")
			return
		else
			throw(err)
		end
	end
end

## used by stop
function playMidiQuick(; port::AbstractString=defaultPlayPort, path::AbstractString="stop")
	try
		run(`pmidi -d 0.1 -p $(port) $(path).mid`)
		sleep(0.1)
	catch err
		if isa(err,InterruptException)
			print("interrupted!\n")
			return
		else
			throw(err)
		end
	end
end

"""
	stop([p]; port=defaultPlayPort, path=\"stop\"

Stops the process `p`, which is assumed to be the one returned by `playMidi`.
If `p` is not given, it will simply kill all instances of the program pmidi.
Plays the file \"`path`.mid\" on ALSA port `port` using the program pmidi.
By default, `path` == \"stop.mid\", which is assumed to be a short midi file that tells a sequencers to be quiet.
Use the function `makeStopMidi` to generate this file in your project directory.
"""
function stop(p; port::AbstractString=defaultPlayPort, path::AbstractString="stop")
	if process_running(p) ; kill(p) ; end
	if !isfile("$(path).mid")
		makeStopMidi(path)
	end
	sleep(0.1)
	playMidiQuick(port=defaultPlayPort,path="stop")
	sleep(0.1)
	print("play stopped\n")
	return
end

function stop()
	try
		run(`killall pmidi`)
	catch err
	end
	if !isfile("stop.mid")
		makeStopMidi("stop")
	end
	sleep(0.1)
	playMidiQuick(port=defaultPlayPort,path="stop")
	sleep(0.1)
	print("play stopped\n")
	return
end

"""
	makeStopMidi(path::AbstractString=\"stop\")

Writes a midi file called \"`path`.mid\" where `path` defaults to \"stop\".
The file contains a midi messages to turn all notes off and stop all sounds.
This function should be run once to make the file in you projec directory so it can be used by the `stop` function.
"""
function makeStopMidi(path::AbstractString="stop")
	x = S(NotesOff(),B(1//32),SoundsOff(),B(1//32))
	makeMidi(x;path=path,bpm=-1)
end

end


text:
using NumericFuns
using Base.Test
using Compat

## auxiliary function

function check_rtype{T<:Number}(f::Functor{1}, ::Type{T})
    pt = result_type(f, T)
    rt = typeof(evaluate(f, one(T)))
    if !(pt == rt)
        error("result_type($f, $T) should be $rt (but predicted $pt)")
    end
end

function check_rtype{T1<:Number,T2<:Number}(f::Functor{2}, ::Type{T1}, ::Type{T2})
    pt = result_type(f, T1, T2)
    rt = typeof(evaluate(f, one(T1), one(T2)))
    if !(pt == rt)
        error("result_type($f, $T1, $T2) should be $rt (but predicted $pt)")
    end
end

## type categories

const inttypes = [Int8,
                  Int16,
                  Int32,
                  Int64,
                  Int128,
                  UInt8,
                  UInt16,
                  UInt32,
                  UInt64,
                  UInt128]

const integraltypes = [Bool, inttypes...]
const fptypes = [Float32, Float64]
const realtypes = [integraltypes..., fptypes...]

const complextypes = [Complex{Bool},
                      Complex{Int16},
                      Complex{Int32},
                      Complex{Int64},
                      Complex64,
                      Complex128]

const numerictypes_r = [realtypes..., Complex64, Complex128]
const numerictypes = [realtypes..., complextypes...]

## test fptype

println("    fptype")

for t in realtypes
    ft = typeof(convert(AbstractFloat, zero(t)))
    @test fptype(t) == ft
end

@test fptype(Complex{Int}) == Complex128
@test fptype(Complex64) == Complex64
@test fptype(Complex128) == Complex128

## arithmetics

println("    arithmetics (unary)")

for F in [Negate, AbsFun, Abs2Fun, RealFun, ImagFun, SqrFun, RcpFun]
    for T in numerictypes
        check_rtype(F(), T)
    end
end

for F in [SignFun, SignbitFun]
    for T in realtypes
        check_rtype(F(), T)
    end
end

println("    arithmetics (binary)")

for F in [Add, Subtract, Multiply, Divide, RDivide, Pow]
    println("        $F")
    for T1 in numerictypes, T2 in numerictypes
        check_rtype(F(), T1, T2)
    end
end

println("    arithmetics (quotient & modulo)")

for F in [DivFun, FldFun, RemFun, ModFun, MaxFun, MinFun]
    println("        $F")
    for T1 in realtypes, T2 in realtypes
        check_rtype(F(), T1, T2)
    end
end

## comparison

println("    comparison operators")

for F in [LT, GT, LE, GE, EQ, NE]
    println("        $F")
    for T1 in realtypes, T2 in realtypes
        check_rtype(F(), T1, T2)
    end
end

## logical & bitwise

println("    logical operators")

check_rtype(Not(), Bool)
check_rtype(And(), Bool, Bool)
check_rtype(Or(),  Bool, Bool)

println("    bitwise operators")

for T in integraltypes
    check_rtype(BitwiseNot(), T)
end

for F in [BitwiseAnd, BitwiseOr, BitwiseXor]
    for T1 in integraltypes, T2 in integraltypes
        check_rtype(F(), T1, T2)
    end
end

# rounding functions

println("    rounding functions")

for F in [FloorFun, CeilFun, TruncFun, RoundFun,
          IfloorFun, IceilFun, ItruncFun, IroundFun]
    for T in realtypes
        check_rtype(F(), T)
    end
end

# number classification functions

println("    number classification")

for F in [IsnanFun, IsinfFun, IsfiniteFun]
    for T in realtypes
        check_rtype(F(), T)
    end
end

# algebraic functions

println("    algebraic functions")

for F in [SqrtFun, RsqrtFun]
    for T in numerictypes
        check_rtype(F(), T)
    end
end

for F in [CbrtFun, RcbrtFun]
    for T in realtypes
        check_rtype(F(), T)
    end
end

for T1 in realtypes, T2 in realtypes
    check_rtype(HypotFun(), T1, T2)
end

println("    exponential & logarithm")

for F in [ExpFun, Exp2Fun, Exp10Fun,
          LogFun, Log2Fun, Log10Fun]
    for T in numerictypes_r
        check_rtype(F(), T)
    end
end

for F in [Expm1Fun, Log1pFun, XlogxFun,
          LogisticFun, LogitFun, SoftplusFun, InvsoftplusFun]
    for T in realtypes
        check_rtype(F(), T)
    end
end

for F in [XlogyFun, LogsumexpFun]
    for T1 in realtypes, T2 in realtypes
        check_rtype(F(), T1, T2)
    end
end

println("    trigonometric functions")

for F in [SinFun, CosFun, TanFun, CotFun, SecFun, CscFun,
          AsinFun, AcosFun, AtanFun, AcotFun, AsecFun, AcscFun,
          SincFun, CoscFun, SinpiFun, CospiFun]
    for T in realtypes
        check_rtype(F(), T)
    end
end

for T1 in realtypes, T2 in realtypes
    check_rtype(Atan2Fun(), T1, T2)
end

println("    hyperbolic functions")

for F in [SinhFun, CoshFun, TanhFun, CothFun, SechFun, CschFun,
          AsinhFun, AcoshFun, AtanhFun, AcothFun, AsechFun, AcschFun]
    for T in numerictypes_r
        check_rtype(F(), T)
    end
end

println("    error functions & friends")

for F in [ErfFun,ErfcFun,ErfiFun,ErfcxFun]
    for T in numerictypes_r
        check_rtype(F(), T)
    end
end

for F in [ErfinvFun, ErfcinvFun]
    for T in realtypes
        check_rtype(F(), T)
    end
end

println("    gamma & friends")

for F in [GammaFun,LgammaFun,EtaFun,ZetaFun]
    for T in numerictypes_r
        check_rtype(F(), T)
    end
end

for T in realtypes
    check_rtype(DigammaFun(), T)
end

println("    beta & lbeta")

for F in [BetaFun, LbetaFun]
    for T1 in realtypes, T2 in realtypes
        check_rtype(F(), T1, T2)
    end
end

println("    airy & friends")

for F in [AiryFun,AiryprimeFun,AiryaiFun,AiryaiprimeFun,AirybiFun,AirybiprimeFun]
    for T in numerictypes_r
        check_rtype(F(), T)
    end
end

println("    bessel & friends")

for F in [Besselj0Fun, Besselj1Fun, Bessely0Fun, Bessely1Fun]
    for T in numerictypes_r
        check_rtype(F(), T)
    end
end


text: include("CTLinearValueIterationIRL.jl")


text:
using JLD
import JLD: save, load

export WritePlan,StorePlan, initialize_storage, store_data!

"""
    WritePlan(file,write_Δt,varlist)

Create a plan for writing data to file. The filename `file` is specified,
to be written to every `write_Δt` time units. The variable names to be written are
specified as a vector of strings in `varlist`.
"""
struct WritePlan
    filen::String
    write_Δt::Float64
    varlist::Vector{String}
end

"""
    save(t,R::WritePlan,v)

Check if time `t` is appropriate for writing to file, according to the WritePlan `R`,
and if so, write the specified variables `v` (which may be separated by commas)
in `R` to the file in `R`.
"""
function save(t,R::WritePlan,v...)
    tol = 1e-8
    if (isapprox(mod(t,R.write_Δt),0,atol=tol) ||
        isapprox(mod(t,R.write_Δt),R.write_Δt,atol=tol))
        save(R.filen,_writelist(R,v)...)
    end
end

"""
    load(R::WritePlan[,var::String])

Load the data stored in the file specified by the WritePlan `R`. If
variable `var` is specified, outputs the specified variable.
"""
load(R::WritePlan) = load(R.filen)
load(R::WritePlan,v...) = load(R.filen,v...)

function _writelist(R::WritePlan,v)
    list = ()
    @assert length(R.varlist)==length(v)
    for (i,vname) in enumerate(R.varlist)
        list = (list...,vname,v[i])
    end
    return list
end

"""
    StorePlan(min_t,max_t,store_Δt,v[;htype=RegularHistory])

Create a plan for storing history data. The storage of data is specified to
start at time unit `min_t` and to proceed until (and including) `max_t`, and
is stored every `store_Δt` time units. The list of variables to be stored is
specified as a list of variables `v`. Tuple-type variables are unwrapped
into separate storage. The optional argument `htype` can be used to set the
history data to `PeriodicHistory` or `RegularHistory` (the default). In the
case of `PeriodicHistory`, the data is assumed to repeat with a period equal
to length(history)+1.
"""
struct StorePlan{H<:HistoryType}
    min_t::Float64
    max_t::Float64
    store_Δt::Float64
    varlist::Vector{DataType}
    StorePlan(min_t,max_t,store_Δt,varlist...;htype=RegularHistory) = new{htype}(min_t,max_t,store_Δt,_get_type(varlist))
end

#function StorePlan(min_t,max_t,store_Δt,v...)
#    return StorePlan(min_t,max_t,store_Δt,_get_type(v))
#end

"""
    initialize_storage(S::StorePlan) -> Vector

Initialize a storage data stack for the storage plan `S`. The output is
an empty vector of `History` vectors.
"""
function initialize_storage(S::StorePlan{H}) where {H}
    data = []
    for v in S.varlist
        push!(data,History(v,htype=H))
    end
    return data
end

"""
    store_data!(data,t,S::StorePlan,v)

Check whether time `t` is a time for saving for storage as described by plan `S`,
and if so, push the variables specified in `v` onto the `data` stack.
"""
function store_data!(data,t,S::StorePlan,v...)
  tol = 1e-8
  if t >= (S.min_t-tol) && t <= (S.max_t + tol) &&
      ((isapprox(mod(t,S.store_Δt),0,atol=tol) ||
        isapprox(mod(t,S.store_Δt),S.store_Δt,atol=tol)))
        store_data!(data,S::StorePlan,v)
  end
  return data
end

function _get_type(u)
    tlist = DataType[]
    _get_type!(tlist,u)
    return tlist
end
function _get_type!(tlist,u::Tuple)
    for ui in u
        _get_type!(tlist,ui)
    end
    return tlist
end
function _get_type!(tlist,u)
    push!(tlist,typeof(u))
    return tlist
end

function store_data!(data,S::StorePlan,u)
    cnt = [0]
    return _store_data!(data,cnt,S,u)
end
function _store_data!(data,cnt,S::StorePlan,u::Tuple)
    for ui in u
        _store_data!(data,cnt,S,ui)
    end
    return data
end
function _store_data!(data,cnt,S::StorePlan,u)
    cnt[1] += 1
    @assert S.varlist[cnt[1]] == typeof(u)
    push!(data[cnt[1]],u)
    return data
end


text:
using Discord
include("token.jl")
include("set.jl")
include("battle.jl")
include("reserve.jl")
include("attack.jl")
include("fin.jl")
include("la.jl")
include("adjust.jl")
include("cancel.jl")

# Create a client.
c = Client(TOKEN,prefix=".")

function set_command(c::Client,m::Message,argument)

    if m.author.bot === true
        return
    end
    #@show m.author.bot

    new_argument = argument
    #全角スペース，全角数字に対応
    replace_pair = ["　" => " ",
                    "０" => "0",
                    "１" => "1",
                    "２" => "2",
                    "３" => "3",
                    "４" => "4",
                    "５" => "5",
                    "６" => "6",
                    "７" => "7",
                    "８" => "8",
                    "９" => "9",
                    "．" => ".",
                    "万" => ""
                    ]
    for pair in replace_pair
        new_argument = replace(new_argument,pair)
    end

    argument_list = split(new_argument," ",keepempty=false)

    if argument_list[1] == ".set" && length(argument_list) == 3
        return set_channel_id(c,m,argument_list[2],argument_list[3])

    elseif argument_list[1] == ".start" || argument_list[1] == ".開始"
        @show now(),m.author.username,new_argument
        if length(argument_list) == 1
            return battle(c,m)
        elseif length(argument_list) == 3
            if isdecimal(argument_list[2]) && isdecimal(argument_list[3])
                return battle(c,m,lap=int(argument_list[2]),boss_number=int(argument_list[3]))
            end
        end

    elseif argument_list[1] == ".reserve" || argument_list[1] == ".予約"
        @show now(),m.author.username,new_argument
        if length(argument_list) == 3
            if isdecimal(argument_list[2]) && isdecimal(argument_list[3])
                return reserve(c,m,int(argument_list[2]),int(argument_list[3]))
            end
        elseif length(argument_list) == 4
            if isdecimal(argument_list[2]) && isdecimal(argument_list[3]) && (argument_list[4] == "over" || argument_list[4] == "持ち越し")
                return reserve(c,m,int(argument_list[2]),int(argument_list[3]),over=true)
            end
        end

    elseif argument_list[1] == ".attack" ||  argument_list[1] == ".凸"
        @show now(),m.author.username,new_argument
        if length(argument_list) == 1
            return attack(c,m)
        elseif length(argument_list) ==2
            if isdecimal(argument_list[2])
                #boss_numberは適当に入れてるだけのクソ実装，現在のボスに凸する場合
                return reserve(c,m,1,int(argument_list[2]),attack=true)
            end
        elseif length(argument_list) == 3
            if isdecimal(argument_list[2]) &&  (argument_list[3] == "over" || argument_list[3] == "持ち越し")
                return reserve(c,m,1,int(argument_list[2]),over=true,attack=true) #持ち越しで現在のボスに凸する場合
            end
        end
    elseif argument_list[1] == ".fin" || argument_list[1] == ".完了"
        @show now(),m.author.username,new_argument
        if length(argument_list) == 2
            if isdecimal(argument_list[2])
                return fin(c,m,int(argument_list[2]))
            end
        end

    elseif argument_list[1] == ".la" || argument_list[1] == ".討伐"
        @show now(),m.author.username,new_argument
        if length(argument_list) == 1
            return  la(c,m)
        elseif length(argument_list) == 2
            return  la(c,m,over_time = argument_list[2])
        end

    elseif argument_list[1] == ".adjust" || argument_list[1] == ".調整"
        @show now(),m.author.username,new_argument
        if length(argument_list) == 2
            if isdecimal(argument_list[2])
                return adjust(c,m,int(argument_list[2]))
            end
        end

    elseif argument_list[1] == ".cancel" || argument_list[1] == ".キャンセル"
        @show now(),m.author.username,new_argument
        if length(argument_list) == 2
            if isdecimal(argument_list[2])
                return cancel(c,m,int(argument_list[2]))
            end
        end
    end

end

add_handler!(
    c,MessageCreate, (c,e) -> set_command(c,e.message,e.message.content)
)

# Log in to the Discord gateway.
open(c)
# Wait for the client to disconnect.
wait(c)

max_stars_repo_name: hibibol/clan-battle

max_issues_repo_name: hibibol/clan-battle

max_forks_repo_name: hibibol/clan-battle


text:
macro lintpragma(s)
end

@lintpragma("Ignore undefined module AutoHashEquals.@auto_hash_equals")
@lintpragma("Ignore undefined module DataStructures.SortedDict")
@lintpragma("Ignore undefined module Lint.lintpkg")
@lintpragma("Ignore undefined module Lists.List")
@lintpragma("Ignore undefined module Lists.ListNode")
@lintpragma("Ignore undefined module Nemo.PowerSeriesRing")
@lintpragma("Ignore undefined module Nemo.QQ")
@lintpragma("Ignore undefined module Nemo.ZZ")
@lintpragma("Ignore undefined module Nemo.coeff")
@lintpragma("Ignore undefined module Nemo.divexact")
@lintpragma("Ignore undefined module Nemo.fac")
@lintpragma("Ignore undefined module Pipe.@pipe")
@lintpragma("Ignore undefined module Reexport.@reexport")
@lintpragma("Ignore undefined module Roots.fzeros")

@lintpragma("Ignore undefined module Brainstorm.dropwhile")
@lintpragma("Ignore undefined module Brainstorm.flatten")
@lintpragma("Ignore undefined module Brainstorm.takewhile")
@lintpragma("Ignore undefined module Brainstorm.Math.checked_add")
@lintpragma("Ignore undefined module Brainstorm.Math.checked_sub")

@lintpragma("Ignore use of undeclared variable _")


text:
using StaticArrays
using IntervalArithmetic
import IntervalArithmetic: interval

using LinearAlgebra

import Base: +, *, ^, -

"""
Affine form with center `c`, affine components `γ` and error `Δ`.

Variant where Δ is an interval
"""
struct AF{N,T<:AbstractFloat}
    c::T   # mid-point
    γ::SVector{N,T}  # affine terms
    Δ::Interval{T}   # error term
end

function Base.show(io::IO, C::AF{N,T}) where {N,T}
    print(io, "⟨", C.c, "; ", C.γ, "; ", C.Δ, "⟩")
end

# ==(C::Affine, D::Affine) = C.c == D.c && C.γ == D.γ

"""
Make an `AF` based on an interval, which is number `i` of `n` total variables.
"""
function AF(X::Interval{T}, n, i) where {T}
    c = mid(X)
    r = radius(X)

    γ = SVector(ntuple(j->i==j ? r : zero(r), n))

    return AF(c, γ, Interval{T}(0))
end

+(x::AF{N,T}, y::AF{N,T}) where {N,T} = AF(x.c + y.c, x.γ .+ y.γ, x.Δ + y.Δ)

-(x::AF{N,T}, y::AF{N,T}) where {N,T} = AF(x.c - y.c, x.γ .- y.γ, x.Δ - y.Δ)

interval(C::AF) = C.c + sum(abs.(C.γ))*(-1..1) + C.Δ

function *(x::AF{N,T}, y::AF{N,T}) where {N,T}
    c = x.c * y.c

    γ = x.c .* y.γ + y.c .* x.γ

    Δ = (x.γ ⋅ y.γ) * (0..1)  # ϵ_i^2

    if N > 1
        Δ += sum(x.γ[i] * y.γ[j] for i in 1:N, j in 1:N if i ≠ j) * (-1..1)  # ϵ_i * ϵ_j
    end

    Δ += (x.c + sum(abs.(x.γ))*(-1..1)) * y.Δ
    Δ += (y.c + sum(abs.(y.γ))*(-1..1)) * x.Δ

    Δ += x.Δ * y.Δ

    return AF(c, γ, Δ)

end

*(x::AF, α::Real) = AF(α*x.c, α.*x.γ, α*x.Δ)
*(α::Real, x::AF) = x * α

+(x::AF, α::Real) = AF(α+x.c, x.γ, x.Δ)
+(α::Real, x::AF) = x + α

-(x::AF) = AF(-x.c, .-(x.γ), -x.Δ)
-(x::AF, α::Real) = AF(x.c - α, x.γ, x.Δ)
-(α::Real, x::AF) = α + (-x)

/(x::AF, α::Real) = AF(x.c/α, x.γ/α, x.Δ/α)

function ^(x::AF, n::Integer)

    invert = false

    if n < 0
        invert = true
        n = -n
        @show n
    end

    result = Base.power_by_squaring(x, n)

    if invert
        result = inv(result)
    end

    return result
end

Base.literal_pow(::typeof(^), x::AF, ::Val{p}) where {T,p} = x^p

x = AF{2,Float64}(0.0, SVector(1.0, 0.0), 0..0)
y = AF{2,Float64}(0.0, SVector(0.0, 1.0), 0..0)

x = AF(3..5, 2, 1)
y = AF(2..4, 2, 2)
#
# 3-x
# interval(3-x)
#
# x + y
#
#
# interval(x+y)
#
# x * y
# interval(x * y)
#
# interval(x * y)
# interval(x) * interval(y)
#
# z = AF(-1..1, 1, 1)
# z^2
# interval(z^2)
#
# using Polynomials
#
# p = Poly([-3, 1])
# p2 = p^8
#
# x = 4 ± 1e-4
# y = AF(x, 1, 1)
#
# interval(y)
# interval(p2(x))
# interval(p2(y))
#
# @time interval(p2(y))
#
#
# f( (x, y) ) = [x^3 + y, (x - y)^2]
#
# X = IntervalBox(-1..1, -1..1)
#
# f(X)
#
# xx = AF(X[1], 2, 1)
# yy = AF(X[2], 2, 2)
#
# interval.(f((xx, yy)))
#
# f(X)
#
#
#
#
# x = AF(4..6, 1, 1)    # example from Messine
# f(x) = x * (10 - x)
#
# f(x)
# interval(f(x))
#
# interval(10*x - x^2)

"General formula for affine approximation of nonlinear functions"
function affine_approx(x::AF, α, ζ, δ)

    c = α * x.c + ζ
    γ = α .* x.γ
    δ += α * x.Δ  # interval

    return AF(c, γ, δ)
end

function Base.sqrt(x::AF, X=interval(x))

    a, b = X.lo, X.hi

    # @show a, b

    # min-range:  de Figuereido book, pg. 64
    α = 1 / (2*√b)
    ζ = (√b) / 2
    δ = ( (√(b) - √(a))^2 / (2*√b) ) * (-1..1)

    return affine_approx(x, α, ζ, δ)
end

function Base.inv(x::AF, X=interval(x))

    a, b = X.lo, X.hi

    # @show a, b

    # min-range:  de Figuereido book, pg. 70
    α = -1 / (b^2)
    d = interval(1/b - α*b, 1/a - α*a)
    ζ = mid(d)
    δ = radius(d) * (-1..1)

    if a < 0
        ζ = -ζ
    end

    return affine_approx(x, α, ζ, δ)
end


text:
import CompactBasisFunctions: vandermonde_matrix, vandermonde_matrix_inverse
import LinearAlgebra: I
import OffsetArrays: OffsetArray

@testset "$(rpad("Vandermonde Matrix",80))" begin

    n = 5
    x = rand(n)
    A = vandermonde_matrix(x)
    B = vandermonde_matrix_inverse(x)

    @test A * B ≈ Matrix(I, n, n)
    @test inv(A) ≈ B
    @test inv(B) ≈ A

    y = OffsetArray(x, 0:n-1)
    C = vandermonde_matrix(y)
    D = vandermonde_matrix_inverse(y)

    @test C * D ≈ Matrix(I, n, n)
    @test inv(C) ≈ D
    @test inv(D) ≈ C

    @test A == C
    @test B == D

end


text:
## Using cubic spline model to define earth's curvature.
## Data points are obtained from: https://www.ngdc.noaa.gov/mgg/topo/report/s6/s6A.html

"""
The area of each pixel is added to determine the total area of a mask. The res parameter is needed for the resolution in arc seconds. By default the resolution is 30 arc seconds.
```
mask = rand(0:1,8000,7100)
mask_area(INDIA_COORDINATE_SYSTEM, mask)
```
"""
function mask_area(geometry::CoordinateSystem, mask, res = 30)
    LAT =[-90.0, -89, -86, -82, -78, -74, -70, -60, -50, -40, -30, -20, -10, 0, 10, 20, 30, 40, 50, 60, 70, 74, 78, 82, 86, 89, 90]
    EW  =[0, 16.0, 64, 133, 193, 256, 318, 465, 598, 712, 804, 872, 914, 928, 914, 872, 804, 712, 598, 465, 318, 256, 193, 133, 64, 16, 0]/1000
    NS  =[931, 931.0, 931, 931, 930, 930, 930, 929, 927, 925, 924, 923, 922, 921, 922, 923, 924, 925, 927, 929, 930, 930, 930, 931, 931,931, 931]/1000 # We divide by thousand to convert into kilometer from meter
    EW = EW / (15 / res) # The values are for 30", our geometry is 15"
    NS = NS / (15 / res)

    MODEL_EW = CubicSpline(LAT, EW)
    MODEL_NS = CubicSpline(LAT, NS)
    area = 0
    for i in 1:geometry.height
        for j in 1:geometry.width
            if mask[i, j] ==1
                lat = row_to_lat(geometry, i)
                ns = MODEL_NS[lat]
                ew = MODEL_EW[lat]
                area = area + ew * ns
            end
        end
    end
    return area
end


text:
get_entropy(snapshot) = get_genotype_distribution(snapshot) |> get_entropy
get_entropy(distribution::Dict) = -sum([p * log2(p) for p in values(distribution)])

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
using DataDeps

register(DataDep(
  "credit-approval",
  "http://archive.ics.uci.edu/ml/datasets/Credit+Approval",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/credit-screening/crx.data",
  "fff49bc186cbddb3ace7371d40d9fbbb3af4f126019c13ff3f562249b1454f4d",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=16,
        feature_indices=1:15,
        categoric_indices=[1; 4:7; 9:10; 12:13],
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
const INT_STR_ID = "@int_str"
const REAL_STR_ID = "@real_str"
const STR_STR_ID = "@str_str"
const TRUE_STR_ID = "@true"
const FALSE_STR_ID = "@false"
const EMPTY_STR_ID = "@empty"

const STRING_TYPE_ID = "%string"
const ARRAY_TYPE_ID = "%array"
const INT_ARRAY_TYPE_ID = "%int_array"
const REAL_ARRAY_TYPE_ID = "%real_array"
const BOOL_ARRAY_TYPE_ID = "%boolean_array"
const STRING_ARRAY_TYPE_ID = "%string_array"

const CONCAT_ID = "@concat"
const COMB_LEN_ID = "@combined_length"
const WRITE_BOOL_ID = "@print_bool"
const COMP_STRS_ID = "@compare_strings"
const GET_ARRAY_SIZE_ID = "@get_array_size"
const GET_ARRAY_PTR_ID = "@get_array_ptr"

const write_bool = """
define void $WRITE_BOOL_ID(i1 %a) {
entry:
  br i1 %a, label %is_true, label %is_false
is_true:
  call i32 (i8*, ...) @printf(i8* getelementptr inbounds ([5 x i8], [5 x i8]* $(TRUE_STR_ID), i32 0, i32 0), i1 %a)
  br label %end
is_false:
  call i32 (i8*, ...) @printf(i8* getelementptr inbounds ([6 x i8], [6 x i8]* $(FALSE_STR_ID), i32 0, i32 0), i1 %a)
  br label %end
end:
  ret void
}
"""

const combined_length = """
define i32 $COMB_LEN_ID($STRING_TYPE_ID* %str1, $STRING_TYPE_ID* %str2) {
  %ptr1 = getelementptr $STRING_TYPE_ID, $STRING_TYPE_ID* %str1, i32 0, i32 1
  %ptr2 = getelementptr $STRING_TYPE_ID, $STRING_TYPE_ID* %str2, i32 0, i32 1
  %X.0 = load i32, i32* %ptr1
  %X.1 = load i32, i32* %ptr2
  %X.2 = add i32 %X.0, %X.1
  %X.3 = add i32 %X.2, 1
  ret i32 %X.3
}
"""

const concat = """
define void $CONCAT_ID($STRING_TYPE_ID* %dest, i8* %new_string, $STRING_TYPE_ID* %head, $STRING_TYPE_ID* %tail) {
  %dest.str.ptr = getelementptr $STRING_TYPE_ID, $STRING_TYPE_ID* %dest, i32 0, i32 0
  store i8* %new_string, i8** %dest.str.ptr
  %head.str.ptr = getelementptr $STRING_TYPE_ID, $STRING_TYPE_ID* %head, i32 0, i32 0
  %head.str = load i8*, i8** %head.str.ptr
  %tail.str.ptr = getelementptr $STRING_TYPE_ID, $STRING_TYPE_ID* %tail, i32 0, i32 0
  %tail.str = load i8*, i8** %tail.str.ptr
  call i8* @strcpy(i8* %new_string, i8* %head.str)
  call i8* @strcat(i8* %new_string, i8* %tail.str)
  ret void
}
"""

const compare_strings = """
define i32 $COMP_STRS_ID($STRING_TYPE_ID* %str1, $STRING_TYPE_ID* %str2) {
  %str1.ptr = getelementptr $STRING_TYPE_ID, $STRING_TYPE_ID* %str1, i32 0, i32 0
  %str1.str = load i8*, i8** %str1.ptr
  %str2.ptr = getelementptr $STRING_TYPE_ID, $STRING_TYPE_ID* %str2, i32 0, i32 0
  %str2.str = load i8*, i8** %str2.ptr
  %res = call i32 @strcmp(i8* %str1.str, i8* %str2.str)
  ret i32 %res
}
"""

const get_array_size = """
define i32 $GET_ARRAY_SIZE_ID($ARRAY_TYPE_ID* %arr) {
  %len.ptr = getelementptr $ARRAY_TYPE_ID, $ARRAY_TYPE_ID* %arr, i32 0, i32 1
  %len = load i32, i32* %len.ptr
  ret i32 %len
}
"""

const get_array_ptr = """
define i8* $GET_ARRAY_PTR_ID($ARRAY_TYPE_ID* %arr) {
  %ptr.ptr = getelementptr $ARRAY_TYPE_ID, $ARRAY_TYPE_ID* %arr, i32 0, i32 0
  %ptr = load i8*, i8** %ptr.ptr
  ret i8* %ptr
}
"""

HELPER_FUNCTIONS = [
  write_bool,
  combined_length,
  concat,
  compare_strings,
  get_array_size,
  get_array_ptr,
]

max_stars_repo_name: juhamyllari/minipascal-julia

max_issues_repo_name: juhamyllari/minipascal-julia

max_forks_repo_name: juhamyllari/minipascal-julia


text:
# ------------------------------------------------------------------------------------------
# # Empezando
#
# Temas:
# 1. Cómo conseguir documentación
# 2. Cómo imprimir a pantalla
# 3. Cómo asignar variables
# 4. Cómo poner comentarios
# 5. Sintáxis para matemáticas básicas
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Cómo conseguir documentación para funciones de Julia
#
# Para conseguir docs para funciones con las que uno no está familiarizado, pon un signo de
# interrogación antes. (¡También funciona en la terminal!)
# ------------------------------------------------------------------------------------------

?println

# ------------------------------------------------------------------------------------------
# ## Cómo imprimir
#
# En Julia usamos `println()` para imprimir texto a la pantalla.
# ------------------------------------------------------------------------------------------

println("¡Estoy emocionado por aprender Julia!")

# ------------------------------------------------------------------------------------------
# Si es tu primera vez usando los notebooks, toma nota que la última línea de cada casilla
# es la que imprime cuando ejecutas la casilla.
# ------------------------------------------------------------------------------------------

123
456

# ------------------------------------------------------------------------------------------
# ## Cómo asignar variables
#
# ¡Sólo necesitas un nombre y un signo de igualdad!<br>
# Julia se encargará de saber el tipo de los datos por nosotros.
# ------------------------------------------------------------------------------------------

my_answer = 42
typeof(my_answer)

my_pi = 3.14159
typeof(my_pi)

my_name = "Jane"
typeof(my_name)

# ------------------------------------------------------------------------------------------
# Después de asignar un valor a una variable, podemos reasignar un valor de un tipo
# diferente a esa variable sin ningún problema.
# ------------------------------------------------------------------------------------------

my_answer = my_name

typeof(my_answer)

# ------------------------------------------------------------------------------------------
# ## Cómo comentar
# ------------------------------------------------------------------------------------------

# Puedes dejar un comentario en una sola línea usando la tecla de gato

#=

Para comentarios de varias lineas,
usa la secuencia de '#= =#' .

=#

# ------------------------------------------------------------------------------------------
# ## Sintáxis para matemáticas básicas
# ------------------------------------------------------------------------------------------

suma = 3 + 7

resta = 10 - 3

producto = 20 * 5

cociente = 100 / 10

potencia = 10 ^ 2

módulo = 101 % 2

# ------------------------------------------------------------------------------------------
# ### Ejercicios
#
# 1.1 Busca los docs para las funciones `convert` y `parse`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 1.2 Asigna `365` a una variable llamada `days`. Convierte `days` a un número flotante.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 1.3 Fíjate que sucede cuando ejecutas
#
# ```julia
# convert(Int64, '1')
# ```
# y
#
# ```julia
# parse(Int64, '1')
# ```
#
# ¿Cuál es la diferencia?
# ------------------------------------------------------------------------------------------

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
# ---
# title: 969. Pancake Sorting
# id: problem969
# author: Indigo
# date: 2021-02-19
# difficulty: Medium
# categories: Array, Sort
# link: <https://leetcode.com/problems/pancake-sorting/description/>
# hidden: true
# ---
#
# Given an array of integers `arr`, sort the array by performing a series of
# **pancake flips**.
#
# In one pancake flip we do the following steps:
#
#   * Choose an integer `k` where `1 <= k <= arr.length`.
#   * Reverse the sub-array `arr[1...k]`.
#
# For example, if `arr = [3,2,1,4]` and we performed a pancake flip choosing `k
# = 3`, we reverse the sub-array `[3,2,1]`, so `arr = [ **1,2,3** ,4]` after the
# pancake flip at `k = 3`.
#
# Return _the`k`-values corresponding to a sequence of pancake flips that sort
# `arr`_. Any valid answer that sorts the array within `10 * arr.length` flips
# will be judged as correct.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [3,2,4,1]
#     Output: [4,2,4,3]
#     Explanation:
#     We perform 4 pancake flips, with k values 4, 2, 4, and 3.
#     Starting state: arr = [3, 2, 4, 1]
#     After 1st flip (k = 4): arr = [ **1, 4, 2, 3** ]
#     After 2nd flip (k = 2): arr = [ **4, 1** , 2, 3]
#     After 3rd flip (k = 4): arr = [ **3, 2, 1, 4** ]
#     After 4th flip (k = 3): arr = [ **1, 2, 3** , 4], which is sorted.
#     Notice that we return an array of the chosen k values of the pancake flips.
#
#
# **Example 2:**
#
#
#
#     Input: arr = [1,2,3]
#     Output: []
#     Explanation: The input is already sorted, so there is no need to flip anything.
#     Note that other answers, such as [3, 3], would also be accepted.
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 100`
#   * `1 <= arr[i] <= arr.length`
#   * All integers in `arr` are unique (i.e. `arr` is a permutation of the integers from `1` to `arr.length`).
#
#
## @lc code=start
using LeetCode

function pancake_sort(arr::Vector{Int})
    len = length(arr)
    res = Int[]
    for i in len:-1:1
        v = @view(arr[1:i])
        agm = argmax(v)
        if i != agm
            append!(res, (agm, i))
            pancake_sort!(arr, (agm, i))
        end
    end
    return res
end
function pancake_sort!(arr::Vector{Int}, ops)
    for op in ops
        reverse!(@view(arr[1:op]))
    end
    return arr
end

## @lc code=end


text:
using Base.Test
using CUDArt
using CUDNN
@show CUDNN_VERSION

import CUDNN.cudnnConvolutionForward
import CUDNN.cudnnConvolutionBackwardFilter
import CUDNN.cudnnConvolutionBackwardData
import CUDNN.cudnnGetConvolutionNdForwardOutputDim
import CUDNN.cudnnGetPoolingNdForwardOutputDim
import CUDNN.cudnnPoolingForward
import CUDNN.cudnnPoolingBackward
include(joinpath("..","src","util","conv_pool_cpu.jl"))

padding=0
stride=1
x = rand(Float32,64,64,3,128); tx = CudaArray(x)
w = rand(Float32,8,8,3,10); tw = CudaArray(w)

dw = zeros(Float32, size(w)); tdw = CudaArray(dw);
dx = zeros(Float32, size(x)); tdx = CudaArray(dx);

dw2 = zeros(Float32, size(w))
dx2 = zeros(Float32, size(x));

y = zeros(Float32, cudnnGetConvolutionNdForwardOutputDim(tx,tw; padding=padding,stride=stride)); ty = CudaArray(y)
y2 = zeros(Float32, size(y))
dy = rand(Float32,size(y)); tdy = CudaArray(dy)

cudnnConvolutionForward(tx,tw,ty; padding=padding, stride=stride); y = to_host(ty)
@time cudnnConvolutionForward(x,w,y2; padding=padding, stride=stride);
@test_approx_eq y y2

cudnnConvolutionBackwardFilter(tx,tdy,tdw); dw = to_host(tdw)
@time cudnnConvolutionBackwardFilter(x,dy,dw2)
@test_approx_eq dw dw2

cudnnConvolutionBackwardData(tw, tdy, tdx); dx = to_host(tdx)
@time cudnnConvolutionBackwardData(w, dy, dx2);
@test_approx_eq dx dx2

using CUDNN: PD, CUDNN_POOLING_MAX, cudnnGetPoolingNdForwardOutputDim
x = rand(Float32,64,64,3,128); tx = CudaArray(x);
psize, padding, stride = 5, 0, 5
pd1 = PD(2, psize, padding, stride, CUDNN_POOLING_MAX)
@assert cudnnGetPoolingNdForwardOutputDim(pd1, tx) == cudnnGetPoolingNdForwardOutputDim(x, window=psize, padding=padding, stride=stride, mode=0)
@show ydims = cudnnGetPoolingNdForwardOutputDim(x, window=psize, padding=padding, stride=stride, mode=0)
y = zeros(Float32, ydims); ty = CudaArray(y);
y2 = zeros(y)
cudnnPoolingForward(tx, ty; window=psize, padding=padding, stride=stride, mode=0); y = to_host(ty);
@time cudnnPoolingForward(x, y2; window=psize, padding=padding, stride=stride, mode=0);
@test_approx_eq y y2

dx = zeros(Float32, size(x)); tdx = CudaArray(dx);
dx2 = zeros(Float32, size(x));
dy = rand(Float32, size(y)); tdy = CudaArray(dy);
cudnnPoolingBackward(ty, tdy, tx, tdx; window=psize, padding=padding, stride=stride, mode=0); dx = to_host(tdx);
@time cudnnPoolingBackward(y, dy, x, dx2; window=psize, padding=padding, stride=stride, mode=0);
@test_approx_eq dx dx2

:ok


text:
#############################################################################
#############################################################################
#
# This file implements polynomial multiplication
#
#############################################################################
#############################################################################

"""
Multiply two polynomials.
"""

function *(p1::Polynomial, p2::Polynomial)::Polynomial
    p_out = Polynomial()
    for t in p1
        p_out = p_out + (t * p2)
    end
    return p_out
end

"""
Multiply two polynomials modulo P.
"""
function *(p1::PolynomialModP, p2::PolynomialModP)::Polynomial
    @assert p1.prime == p2.prime "Primes must be the same"
    p_out = Polynomial()
    for t in p1.terms
        p_out = p_out + (t * p2.terms)
    end
    return mod(p_out, p1.prime)
end

"""
Power of a polynomial.
"""

function ^(p::Polynomial, n::Int)
    n < 0 && error("No negative power")
    out = one(p)
    binary_arr = digits(n, base=2)
    for i in 1:length(binary_arr)
        binary_arr[i] == 1 ? out = out*p : out *= 1
        p = p*p
    end
    return out
end

# I decided to keep the ^ as a general function and not restrict it to only polynomials since this also works with raising integers and other types to a power n mod p.
function ^(p, n::Int, prime::Int)
    n < 0 && error("No negative power")
    out = one(p)
    binary_arr = digits(n, base=2)
    p = mod(p, prime)
    for i in 1:length(binary_arr)
        binary_arr[i] == 1 ? out = mod(out*p, prime) : out *= 1
        p = mod(p*p, prime)
    end
    return mod(out, prime)
end

^(p::PolynomialModP, n::Int) = ^(p.terms, n, p.prime)

function CRT(poly_arr, prime_arr)

    c = 0
    max_input_deg = max(degree(poly_arr[1]), degree(poly_arr[2]))+1
    exponent_arr_1 = zeros(Int64, 1, max_input_deg)
    exponent_arr_2 = zeros(Int64, 1, max_input_deg)
    for i in poly_arr[1]
        exponent_arr_1[i.degree+1] = i.coeff
    end
    for i in poly_arr[2]
        exponent_arr_2[i.degree+1] = i.coeff
    end
    """
    The lines above actually result in a marginally slower CRT implementation as opposed to the approach where you cycle through exponent arrays 1 and 2 and check if k matches the exponent value before sending those values off to the iCRT call. I chose to keep this approach because the computational time difference was only marginal and I felt like it shows that I *really* tried to tackle the task given.
    """

    for k in max_input_deg-1:-1:0
        ak = exponent_arr_1[k+1]
        bk = exponent_arr_2[k+1]
        ck = iCRT([ak, bk], prime_arr)
        k == 0 ? c = c + ck : c = c + ck * x^k # Can probably remove the turnery, its only there because I had a strange bug with Term(1,1)
    end
    return c
end

function mult_poly_with_crt(a, b)
    height_a = maximum(coeffs(a))
    height_b = maximum(coeffs(b))
    B = 2*height_a*height_b*min(degree(a)+1, degree(b)+1)
    M = 3
    c = poly_modP_multiplication(a, b, M)
    while M < B
        nextprime(M) == M ? p = nextprime(M, 2) : p = nextprime(M)
        d = poly_modP_multiplication(a, b, p)
        c = CRT([c, d], [M, p])
        M = M*p
    end
    return c
end

function iCRT(u, m)
    v = Vector{Int}(undef, 2)
    v[1] = u[1]
    v[2] = mod((u[2] - v[1])*inverse_mod(m[1], m[2]), m[2])
    output = u[1] + v[2]*m[1]
    smod(output, m[1]) == smod(u[1], m[1]) && smod(output, m[2]) == smod(u[2], m[2]) && return smod(output, m[1]*m[2])
    error("Output modulo prime is not consistent with inputs")
end

function poly_modP_multiplication(a::Polynomial, b::Polynomial, prime::Int)
    if length(a) == 0 || length(b) == 0
        return 0
    end
    if leading(a) == a.terms[1] && leading(b) == b.terms[1]
        return mod(mod(a, prime) * mod(b, prime), prime)
    end

    if length(a.terms) == 1
        a1 = a
        a2 = Polynomial()
    else
        a1 = Polynomial(a.terms[1:(length(a.terms)÷2)])
        a2 = Polynomial(a.terms[(length(a.terms)÷2)+1:length(a.terms)])
    end

    if length(b.terms) == 1
        b1 = b
        b2 = Polynomial()
    else
        b1 = Polynomial(b.terms[1:(length(b.terms)÷2)])
        b2 = Polynomial(b.terms[(length(b.terms)÷2)+1:length(b.terms)])
    end

    return mod(poly_modP_multiplication(a1, b1, prime) + poly_modP_multiplication(a1, b2, prime) + poly_modP_multiplication(a2, b1, prime) + poly_modP_multiplication(a2, b2, prime), prime)
end

function mult_vec_el_up_to(arr::Vector, el::Int)
    @assert length(arr) > 0 "Array must contain at least 1 element."
    el == 0 ? 1 : arr[el]*mult_vec_el_up_to(arr, el-1)
end

function i_ext_euclid_alg(a,b)
    a == 0 && return b, 0, 1
    g, s, t = i_ext_euclid_alg(b % a, a)
    return g, t - (b ÷ a)*s, s
end

smod(a::Int,m::Int)::Int = mod(a,m) > m ÷ 2 ?  mod(a,m) - m : mod(a,m)

inverse_mod(a,m) = mod(i_ext_euclid_alg(a,m)[2],m);