text:
# License for this file: MIT (expat)
# Copyright 2017-2021, DLR Institute of System Dynamics and Control
# Author: Martin Otter, DLR-SR

#=
This file contains functions to transform the ODAE generated by the Pantelides
algorithm to an ODE, if this is reasonably possible
(e.g. without solving nonlinear equation systems inside the model or introducing stabilizing terms.
The approach is based on the dummy derivative method from
[Mattsson and Söderlind (1992](https://ieeexplore.ieee.org/document/274429),
the transformation algorithm from
[Otter and Elmqvist (2017)](https://modelica.org/events/modelica2017/proceedings/html/submissions/ecp17132565_OtterElmqvist.pdf),
and on some new ideas that are not yet published.

# Functions exported by this module:

- [`getSortedAndSolvedAST`](@ref): Given the BLT
  (after Pantelides algorithm),
  this function returns the Abstract Syntax Tree (AST) of
  the sorted and solved equations and selects the ODE states.
=#

export getSortedAndSolvedAST
export StateSelectionFunctions
export showCodeWithoutComments

#using OrderedCollections
using  ModiaBase
using  ModiaBase.BLTandPantelides
import ModiaBase: revertAssociation
using  StaticArrays

const INFORMATION = 0
const WARNING     = 1
const ERROR       = 2

"""
    stateSelectionFunctions = StateSelection.StateSelectionFunctions(;
        var_name               = v          -> nothing,
        var_julia_name         = v          -> nothing,
        var_unit               = v          -> "",
        var_startInitFixed     = v          -> (nothing,false),   # (startOrInit, fixed)
        var_is_state           = v_original -> false,
        equation               = e          -> ""
        isSolvableEquation     = (e_original,v_original) -> false,
        isLinearEquation       = (e_original,v_original) -> (false, false),
        getSolvedEquationAST   = (e, v)                  -> nothing,
        getResidualEquationAST = (e, residualName)       -> nothing,
        showMessage            = (message; severity=0,from="???",details="",variables=Int[],equations=Int[]) -> println("...")
    )

Generate an instance of an immutable struct that holds the callback functions
as needed by [`getSortedAndSolvedAST`](@ref).

The functions need to have the following arguments:

- `var_name(v::Int)::String`:\\
   Return full (Modia) path name of `v` as String (e.g. `"a.b.c"` or `"a.b.der(c)"`).

- `var_julia_name(v::Int)::String`:\\
  Return full Julia path name of variable `v` as `Symbol`
  (e.g. Symbol("m.a.b.c") or Symbol("a_b_der_c")).

- `var_unit(v::Int)::String`:\\
  Return unit of variable `v` as `String` (for example "N*m") or `""`, if no unit is defined.

- `var_startInitFixed(v::Int)::Tuple(Any,Bool)`:\\
  Return `(startOrInit, fixed)`, where `startOrInit` is the `start` or `init` value or `nothing`
  (if neither `start` nor `init` defined) and `fixed` is `true`, if an `init` value is defined.

-  `var_is_state(v_original::Int)::Bool`:\\
   Return true, if variable `v_original` is defined to be a state.

- `equation(e::Int)`:\\
   Return equation as string or "", if equation is not provided.

- `isSolvableEquation(e_original::Int, v_original::Int)::Bool`:\\
   Return `true` if original, undifferentiated equation `e_original`
   can be uniquely solved for variable `v_original`.

- `(isLinear, hasConstantCoefficient) = isLinearEquation(e_original::Int, v_original::Int)`:\\
   Return the information whether the original, undifferentiated equation `e`
   is **linear** with respect to variable `v` (`isLinear=true`) and if this is the case whether
   the coefficient of this linear term is **constant** after initialization
   (`hasConstantCoefficient=true`). If `v` is an array, then the array equation `e` must
   be linear with respect to all elements of array `v` in order to return `isLinear=true`.

- `AST = getSolvedEquationAST(e:Int, v:Int)`:\\
  Return Abstract Syntax Tree `AST::Expr` of equation `e` solved for variable `v`
  as instance of `Expr`.

- `AST = getResidualEquationAST(e::Int, residualName::Symbol)`:\\
  Return Abstract Syntax Tree `AST::Expr` of eqation `e` as residual equation with ustrip(..) \\
  `residualName = ustrip( < residual form of equation e> )`

- `showMessage(message; severity=0,from="???",details="",variables=Int[],equations=Int[])`:\\
  Print information, warning or error message.
"""
struct StateSelectionFunctions
    var_name::Function
    var_julia_name::Function
    var_unit::Function
    var_startInitFixed::Function
    var_nominal::Function
    var_unbounded::Function
    var_is_state::Function
    equation::Function
    isSolvableEquation::Function
    isLinearEquation::Function
    getSolvedEquationAST::Function
    getResidualEquationAST::Function
    showMessage::Function

    StateSelectionFunctions(;
        var_name               = v          -> nothing,
        var_julia_name         = v          -> nothing,
        var_unit               = v          -> "",
        var_startInitFixed     = v          -> (nothing,false),   # (startOrInit, fixed)
        var_nominal            = v_original -> NaN,
        var_unbounded          = v_original -> false,
        var_is_state           = v_original -> false,
        equation               = e          -> "",
        isSolvableEquation     = (e_original,v_original) -> false,
        isLinearEquation       = (e_original,v_original) -> (false, false),
        getSolvedEquationAST   = (e, v)                  -> nothing,
        getResidualEquationAST = (e, residualName)       -> nothing,
        showMessage            = (message; severity=0,from="???",details="",variables=Int[],equations=Int[]) ->
                                       println("\n!!! Message from $from with severity=$severity:\n$message,\ndetails\nvariables=$variables,\nequations=$equations")
    ) = new(var_name,
            var_julia_name,
            var_unit,
            var_startInitFixed,
            var_nominal,
            var_unbounded,
            var_is_state,
            equation,
            isSolvableEquation,
            isLinearEquation,
            getSolvedEquationAST,
            getResidualEquationAST,
            showMessage)
end

function printVector(vec)
    for v in vec
        println("    ", v)
    end
end

"""
    showMessage(message, severity, from, details, variables, equations)

Print message.

# Arguments:

- `severity::Int`: Severity of message: 0=Information, 1=warning, 2=Error
- `from::String`: Name of function that throws error message
- `message::String`: One line description of issue
- `details::String`: Detailed description of issue
- `variables::Vector{String}`: List of involved variables
- `equations::Vector{String}`: List of involved equations
"""
function showMessage(message::String, severity::Int, from::String, details::String,
                     variables::AbstractVector, equations::AbstractVector)::Nothing
    if severity == 0
        println("\nInformation message from $(from):")
    elseif severity == 1
        println("\nWarning message from $(from):")
    else
        println("\nError message from $(from):")
    end
    println(message)
    if length(details) > 0
        println(details)
    end
    if length(variables) > 0
        println("Involved variables:")
        printVector(variables)
    end
    if length(equations) > 0
        println("Involved equations:"),
        printVector(equations)
    end
	println("")
    return nothing
end

showCodeWithoutComments(code) = println("code = ", replace(sprint(show,code), r"( *#= .*=#\n)|(#= .*=#)" => "") )

"""
    (eConstraints,vConstraints) = getConstraintSets(BLT_i, assignRev, Arev, Brev,
                                       var_is_state, showMessage)

Determines the set of equations/constraints and their unknowns that are associated with
the highest derivative level BLT block i (BLT_i), according to the dummy derivative
method of [Mattsson and Söderlind (1992](https://ieeexplore.ieee.org/document/274429),
(note, eConstraints[j] and vConstraints[j] are of type Vector{Int}):

- eConstraints[end] = BLT_i, so the equations of BLT block i (must be highest derivative equations)
- eConstraints[j]: Constraint equations that form the integral of eConstraints[j+1].
- eConstraints[1]: Lowest order constraint equations. All equations in eConstraints[1]
                   are original, non-differentiated equations.

vConstraints[j] are the unknowns of eConstraints[j].

# Arguments:

- `BLT_i::Vector{Int}`: Vector of equation indices of BLT block i.

- `assignRev::Vector{Int}`: Reverted assign vector: vi = assignRev[ei] is variable vi assigned to equation ei

- `Arev::Vector{Int}`: Reverted variable association: `Arev[i] = der(V[k]) == V[i] ? k : 0`.

- `Brev::Vector{Int}`: Reverted equation association: `Erev[i] = der(E[k]) == E[i] ? k : 0`.

- `var_is_state::Function`: var_is_state(`vi`) returns true, if `vi` is defined to be a definite state.

- showMessage: showMessage2 below to print error/warning message

"""
function getConstraintSets(BLT_i::Vector{Int},
                           assignRev::Vector{Int},
                           Arev::Vector{Int},
                           Brev::Vector{Int},
                           var_is_state::Function,
                           showMessage::Function)
    # Determine unknowns of BLT block i
    BLT_i_unknowns = fill(0, length(BLT_i))
    for j in eachindex(BLT_i)
        BLT_i_unknowns[j] = assignRev[BLT_i[j]]
    end

    # Initialize constraint sets
    eConstraints = [BLT_i]
    vConstraints = [BLT_i_unknowns]

    # Check that number of equations and number of unknowns agree
    if length(BLT_i) != length(BLT_i_unknowns)
        showMessage("Number of equations and number of unknowns do not agree for BLT block";
                    severity  = ERROR,
                    variables = BLT_i_unknowns,
                    equations = BLT_i)
        return (eConstraints, vConstraints, false)
    end

    while true
        # Determine constraints at one differentiation order less
        ceq = fill(0, 0)
        for eq in eConstraints[1]
            if Brev[eq] > 0
                push!(ceq, Brev[eq])
            end
        end
        if length(ceq) == 0; break; end

        pushfirst!(eConstraints, ceq)   # move ceq to the beginning of the constraints vector

        # Determine unknowns of ceq
        veq        = fill(0,0)
        veq_states = fill(0,0)
        for vc in vConstraints[1]
            if Arev[vc] > 0
                if var_is_state( Arev[vc] )
                    push!(veq_states, Arev[vc])
                else
                    push!(veq, Arev[vc])
                end
            end
        end
        if length(veq) < length(ceq)
            nRemoveStates = length(veq_states) - (length(ceq) - length(veq))
            showMessage("$nRemoveStates of the followig defined states cannot be states!";
                severity  = ERROR,
                variables = veq_states,
                equations = ceq)
            return (eConstraints, vConstraints, false)
        end
        pushfirst!(vConstraints, veq)   # move veq to the beginning of the constraints vector
    end

    @assert(length(eConstraints) == length(vConstraints))
    return (eConstraints, vConstraints, true)
end

"""
    names = getNames(var_name::Function, v::Vector{Int})
    names = getNames(eq::EquationGraph , v::Vector{Int})

Get names of variables `v` as a string.
"""
function getNames(var_name::Function, v)::String
    if length(v) > 0
        str = var_name(v[1])
        for i in 2:length(v)
            str = str * ", " * var_name(v[i])
        end
        return str
    else
        return ""
    end
end
getNames(eq, v) = getNames(eq.fc.var_name, v)

function printVariables(eq, v; indent=4)
    for vi in v
        println(lpad(vi,indent), ": ", eq.fc.var_name(vi))
    end
    return nothing
end

"""
    printEquations(equation::Function, e::Vector{Int}; indent=4)
    printEquations(eq::EquationGraph, e::Vector{Int}; indent=4)
    printEquations(eq::EquationGraph, expr; indent=4)

Print equations `e` as a string or as vector of equations.
"""
function printEquations(equation::Function, e::Vector{Int}; indent=4)::Nothing
    for ei in e
        println(lpad(ei,indent), ": ", equation(ei))
    end
    return nothing
end
printEquations(eq, e; indent=4) = printEquations(eq.fc.equation, e; indent=indent)

function printEquations(expr::Expr; indent=4)
    indentation = repeat(" ",indent)
    str = string(expr)
    str = replace(str, r".*#= .*\n" => "")
    println(indentation, str)
end

#=

"""
    Arev = revertAssociation(A::Vector{Int})

Revert the association Vector `A[i] = j`, such that
`Arev[j] = i` (`A[i]=0` is allowed and is ignored).
"""
function revertAssociation(A::Vector{Int})::Vector{Int}
    nArev = max(maximum(A), length(A))
    Arev  = fill(0, nArev)
    for i in eachindex(A)
        if A[i] != 0
            Arev[ A[i] ] = i
        end
    end
    return Arev
end
=#

"""
    tearedEquations = TearedEquations(..)

Return an immutable struct that contains the information about a teared system of equations.
"""
struct TearedEquations
    isLinear::Bool                 # = true if equation system is linear
    hasConstantCoefficients::Bool  # = true if linear equation system has a constant coefficient matrix A
    vSolved::Vector{Int}           # Variables that are explicitly solved by tearing
    vTear::Vector{Int}             # Iteration variables of tearing
    eSolved::Vector{Int}           # Equations that are explicitly solved by tearing
    eResiduals::Vector{Int}        # Residual equations of the teared equations
    AST::Expr                      # AST of equation system
end

"""
    eq = EquationGraph(G, BLT, assign, A, B, log, showMessage::Function)

Construct the representation of the equation graph as needed by state selection.

# Arguments

- `G`: G[i] is the Integer vector of variable indices of equation i.
   A variable i can be a scalar or an array.
   G is of type Vector{ Vector{Int} } or Vector{ Any }.

- `BLT::Vector{Vector{Int}}`: BLT[i] is the vector of equations belonging to BLT-block i.
  In principal, **BLT[i] must include only highest derivative equations**.
  In order to ease the generation of BLT[i], to not be forced to map indices forth and back
  when calling the function that computes the BLT, the following assumption is
  made: If BLT[j] contains only **one equation with B[j][1] > 0** then this BLT block is
  a lower-derivative (dummy) block and is ignored.

- `assign::Vector{Int}`: ei = assign[vi] is equation ei assigned to variable vi

- `A`: A-Vector of Pantelides algorithm:
       `A[i] = if der(v[i]) == v[k] then k else 0`
       where `v[i]` is variable `i`.

- `B`: B-Vector of Pantelides algorithm:
       `B[i] = if der(e[i]) == e[k] then k else 0`
       where `e[i]` is equation `i`.

- `log`: = true, if logging is on

- `showMessage`: Function used to print warning/error messages.
"""
mutable struct EquationGraph
    success::Bool                                  # = true, if operation was successful, otherwise false
    log::Bool                                      # = true, if logging is on
    logDetails::Bool                               # = true, if detailled logging is on
    showMessage::Function
    fc::StateSelectionFunctions                    # Callback functions
    G                                              # Structure of equations
    A::Vector{Int}
    B::Vector{Int}
    Gunknowns::Vector{Vector{Int}}
    vOriginal::Vector{Bool}                        # vOriginal[v] = true, if variable v appears in an original (non-differentiated) equation
    vActive::Vector{Bool}                          # vActive[v] = false, if v is a state
    vSolvedWithFixedTrue::Vector{Int}              # Explicitly solved variables with fixed=true

    # Reverted info vectors
    Arev::Vector{Int}                              # Reverted vector A
    Brev::Vector{Int}                              # Reverted vector B
    assignRev::Vector{Int}                         # Reverted vector assign; vi = assignRev[ei] if equation ei is assigned to variable vi

    # Equation/Constraint sets
    eConstraintsVec::Vector{Vector{Vector{Int}}}   # eConstraintsVec[i] is the constraint set of highest derivative BLT block i
    vConstraintsVec::Vector{Vector{Vector{Int}}}   # vConstraintsVec[i] are the assigned variables of eConstraintsVec[i]

    # Tearing
    td::TearingSetup                               # Data structure for equation tearing to traverse a set of equations
                                                   # that are represented as a DAG (Directed Acyclic Graph).
    vSolved::Vector{Int}                           # Variables that are explicitly solved by tearing
    vTear::Vector{Int}                             # Iteration variables of tearing
    eSolved::Vector{Int}                           # Equations that are explicitly solved by tearing
    eResiduals::Vector{Int}                        # Residual equations of the teared equations
    vTearWithoutStart::Vector{Int}                 # Summary of all vTear variables that have no start value (to print warning message at the end)

    der_vSolvedLower::Vector{Int}                  # Initialize tearing with variables that are already explicitly solved
                                                   # (this is known from the constraints one differentiation level less)
    der_vTearLower::Vector{Int}                    # Highest priority variables to solve for
                                                   # (= differentiated tearing variables from one differentiation level less)
                                                   # (this is known from the constraints one differentiation level less)
    der_eSolvedLower::Vector{Int}                  # Initialize tearing with equations that are already explicitly solved
                                                   # (this is known from the constraints one differentiation level less)
    vc::Vector{Int}                                # Variables of the current component that are unknowns of the tearing
    ec::Vector{Int}                                # Equations of the current component that shall be teared
    vcCandidates::Vector{Vector{Int}}              # Try to solve equations ec in the following order for vc:
                                                   # vcCandidates[1], ..., [6]
    vLinActive::Vector{Bool}                       # Utility vector (vLinActive[i]=true, if variable is inspected)
    equationInfo::Modia.EquationInfo               # Info about the structure of the equations
                                                   # (equationInfo.x_names, equationInfo.linearEquations)

    # Sorted equations
    fullAssignRev::Vector{Int}                     # fullAssignRev[e] = v: If e is a scalar equation system, it is solved for v.
                                                   # If v = -1, the scalar equation system is not solved directly, but
                                                   #            as linear equation system (with tearing).
    tearedEquations::Vector{TearedEquations}       # Info about teared equation system
    tearedEquations_indices::Vector{Int}           # teq = tearedEquations_indices[e] is tearedEquations[teq]
    eAST::Vector{Expr}                             # eAST[e] is the AST of equation e
    AST_aux::Vector{Expr}                          # Abstract Syntax Tree of some code (auxiliary vector)
    AST::Vector{Expr}                              # Abstract Syntax Tree of the sorted and solved equations as vector of Expr
                                                   # Statements are evaluted in the order: AST[1], AST[2], ....

    function EquationGraph(G, BLT, assign, A, B, stateSelectionFunctions, log, logDetails, showMessage::Function)

        # Check input arguments
        @assert(length(B) == length(G))

        # Revert association vectors
        Arev      = revertAssociation(A)
        Brev      = revertAssociation(B)
        assignRev = revertAssociation(assign)

        # Initialize vector of variables of the non-differentiated, original equations
        vOriginal = fill(false, length(A))

        # Initialize the active variables
        vActive = fill(true, length(A))
        for v = 1:length(vActive)
            if Arev[v] > 0
                vActive[ Arev[v] ] = false    # Arev[v] is a potential state
            end
        end

        # Define empty constraint sets
        eConstraintsVec = Vector{Vector{Int}}[]
        vConstraintsVec = Vector{Vector{Int}}[]

        # Gunknowns has the same dimension as G.
        # Gunknowns[i] is not empty, if it is an original, non-differentiated
        # equation. In this case it is a subset of G[i] containing only the
        # unknown variables of the corresponding constraint set as needed by tearing.
        # Note, tearing is only applied on non-differentiated equations
        # Tearing information of differentiated equations is
        # deduced from the tearing of the non-differentiated equations.
        Gunknowns = fill(Int[], length(G))

        # Construct constraint sets of the dummy derivative method from the BLT
        success = true
        for (i,BLT_i) in enumerate(BLT)
            # Check that BLT block contains only higher derivative equations
            # (with one exception due to an agreed convention with Modia:
            # If BLT[j] has only one equation and B( BLT[j][1] ) > 0, so it is
            # a lower derivative equation, then this BLT block is ignored).
            highest = true
            for eqj in BLT_i
                if B[eqj] != 0 # One equation of block BLT_i appears also differentiated
                    highest = false
                    if length(BLT_i) > 1
                        success = false
                        showMessage("Equation $eqj in BLT[$i]=$BLT_i is not a highest derivative equation.";
                        details   = "This means that there is an error in argument BLT.",
                        severity  = ERROR,
                        equations = [eqj])
                    end
                    break
                end
            end
            if !success
                break
            end

            if highest
                # Get all equation sets eConstraints and their corresponding unknowns vConstraints
                # from lowest to highest differentiation order (eConstraints[end] is c)
                (eConstraints, vConstraints, success2) = getConstraintSets(BLT_i, assignRev, Arev, Brev,
                                                         stateSelectionFunctions.var_is_state, showMessage)
                if !success2
                    success=false
                    break
                end

                # Initialize vActive
                #for vc in vConstraints
                #    for v in vc
                #        if Arev[v] > 0
                #            vActive[ Arev[v] ] = false    # Arev[v] is a potential state
                #        end
                #    end
                #end

                push!(eConstraintsVec, eConstraints)
                push!(vConstraintsVec, vConstraints)

                # Determine unknown variables of G as used by tearing
                for i in eachindex(eConstraints)
                    for eq in eConstraints[i]
                        # Gunknowns[eq] contains the variables of G[eq] that are unknowns of vConstraints[i]
                        Gunknowns[eq] = intersect(G[eq], vConstraints[i])

                        if Brev[eq] == 0
                            # Original, non-differentiated equation
                            for v in Gunknowns[eq]
                                # v is a variable from the original, non-differentiated equations
                                vOriginal[v] = true
                            end
                        end
                    end
                end
            end
        end

        new(success, log, logDetails, showMessage, stateSelectionFunctions, G, A, B, Gunknowns, vOriginal, vActive,
            Int[], Arev, Brev, assignRev,
            eConstraintsVec, vConstraintsVec,
            TearingSetup(Gunknowns,length(Arev)),
            Int[], Int[], Int[], Int[], Int[], Int[], Int[], Int[], Int[], Int[],
            Vector{Int}[Int[], Int[], Int[], Int[], Int[], Int[]], fill(false, length(Arev)),
            Modia.EquationInfo(status=Modia.CODE_GENERATION),
            fill(0, length(B)), TearedEquations[], fill(0, length(B)),
            Vector{Expr}(undef, length(B)), Expr[], Expr[])
    end
end

function undifferentiated(eq::EquationGraph, v::Int)::Int
    while eq.Arev[v] > 0
        v = eq.Arev[v]
    end
    return v
end

"""
    isLinearOriginalOrDifferentiatedEquation(eq, e, v)
"""
function isLinearOriginalOrDifferentiatedEquation(eq::EquationGraph, e::Int, v::Int)
    # Determine original, undifferentiated equation
    while eq.Brev[e] > 0
        e = eq.Brev[e]
        v = eq.Arev[v]
    end

    # Check whether original, undifferentiated equation is linear
    return eq.fc.isLinearEquation(e, v)
end

"""
    (isLinear, hasConstantCoefficients) = isLinearEquationSystem!(eq, es, vs)

Returns `isLinear = true`, if equations `es::Vector{Int}` are linear with respect to
variables `vs::Vector{Int}`. If **all** coefficients of the linear factors are additionally constant,
`hasConstantCoefficients=true` is returned, otherwise false.
`eq::EquationGraph` is the data structure holding all information about the graph.
"""
function isLinearEquationSystem!(eq::EquationGraph, es::Vector{Int}, vs::Vector{Int})::Tuple{Bool,Bool}
    # It is assumed, that on input, eq.vLinActive[i] = false
    # On return, eq.vLinActive[i] = false again

    # Mark that vs are active variables
    eq.vLinActive[vs] .= true

    # Check whether es is linear with respect to vs
    isLinear::Bool = true
    hasConstantCoefficients::Bool = true
    for e in es
        for v in eq.G[e]
            if eq.vLinActive[v]  # Only check variables vs in G[e]
                (e_isLinear, e_hasConstantCoefficient) = isLinearOriginalOrDifferentiatedEquation(eq, e, v)
                if !e_isLinear
                    isLinear = false
                    hasConstantCoefficients = false
                    break
                else
                    hasConstantCoefficients = hasConstantCoefficients && e_hasConstantCoefficient
                end
            end
        end
        if !isLinear
            break
        end
    end

    # Unmark that vs are active variables
    eq.vLinActive[vs] .= false

    return (isLinear, hasConstantCoefficients)
end

"""
    tearEquationsWithCandidates!(eq::EquationGraph)

Tear equations (eq.der_eSolvedLower, eq.ec) with respect to variables
(eq.der_vSolvedLower, eq.vc) and initialize tearing with the
Directed-Acylic-Graph (eq.der_eSolvedLower, eq.der_vSolvedLower).

It is tried to solve the equations eq.ec with respect to variables `eq.vc`
in the following order:

1. eq.der_vTearLower\\
   (These are differentiated variables from the previous differentiation level that
    have been used as tearing variables at the previous differentation level.
    In general, it is adviceable to use non-differentiated
    variables as states and by trying to eliminate der_vTearLower first, it is tried to solve
    for all differentiated variables first; note these variables have no start values).

2. Variable has start=nothing and the variable does not appear differentiated
   in the set of original equations.
   (Variables without a start value should be explicitly computed during initialization and
   therefore these variables should be solved for;
   however, if a variable appears also differentiated in the original equations,
   then it should be tried to use it as a state and therefore, it should be first tried
   to not solve for this variable).

3. Variable has start=nothing and the variable appears differentiated
   in the set of original equations.
   (Variables without a start value should be explicitly computed during initialization and
   therefore these variables should be solved for).

4. Variable has fixed=false\\
   (The start value is a guess value, so it is safe to solve for the variable and
   hereyby ignore the start value)

5. Variable has fixed=nothing\\
   (The intention is to respect the start value if possible, that is, it should be avoided
   to solve for such a variable)

6. Variable has fixed=true\\
   (The start value must be respected, that is, if any possible no equation should be solved
   for this variable; if it is solved for, an error is triggered during initialization if the
   start value and the computed value do not agree)

# Example

Assume the constraint `phi1 = phi2` and differentiated constraints
`der(phi1) = der(phi2), w1 = der(phi1), w2 = der(phi2)` and start values are defined
for `phi1` and `w1`.

- Since no start value is defined for `phi2` (but a start value is defined for `phi1`),
  equation `phi1 = phi2` is solved for `phi2`, so `phi2 := phi1` and `phi1` is the tearing variable.
  This means that `der_vSolvedLower = der(phi1); der_vTearLower = der(phi2)`.

- On the next differentiation level, the differentiated equation `der(phi1) = der(phi2)`
  is solved in the same way as on the previous level due to `der_eSolvedLower, der_vSolvedLower`, that is
  `der(phi2) := der(phi1)`.

- Since `der_vTearLower = der(phi1)`, it is first tried to solve for `der(phi1)`. There is only
  one equation left  with `der(phi1)` and therefore this equation is solved for it:
  `der(phi1) := w1`.

- Since a start value is defined for `w1` and no start value for `w2`, it is tried to solve
  for `w2`, that is: `w2 := der(phi2)` and the result is that `w1` is used as tearing variable.

Resulting equation set with states = `phi1, w1`:

```julia
phi2      := phi1
der(phi1) := w1
der(phi2) := der(phi1)
w2        := der(phi2)
```
"""
function tearEquationsWithCandidates!(eq::EquationGraph)::Nothing
    # Empty all candidate vectors
    for c in eq.vcCandidates
        empty!(c)
    end

    # Determine priorities to solve for variables eq.vc
    #   Assume equations:
    #     (1)   s = ...
    #     (1')  w = der(s)
    #     (1'') a = der(w)
    #   Then it is most natural to use "w" and "a" as tearing variables if "s" is a tearing variable
    #   and solve for "der(phi)" and "der(w)".
    #   This means, that eq.der_vTearLower should be tried first to be solved
    eq.vcCandidates[1] = eq.der_vTearLower

    for v in setdiff(eq.vc, eq.der_vTearLower)
        # Variable has no start value because:
        #   (a) is not from the original, non-differentiated equations
        #   (b) is a differentiated variable
        #   (c) no start or init value is defined
        (startOrInit, fixed) = eq.fc.var_startInitFixed(v)
        hasNoStart = !eq.vOriginal[v] || eq.Arev[v] > 0 || startOrInit isa Nothing
        if hasNoStart
            # Has no start value
            if eq.vOriginal[v] && (eq.A[v]>0 && eq.vOriginal[eq.A[v]])
                # Is from the original, non-differentiated equations and its derivative
                # is also in the original equations (so try to use it as state).
                push!(eq.vcCandidates[3], v)
            else
                # All other variables without a start value are tried
                push!(eq.vcCandidates[2], v)
            end

        else
            if startOrInit isa Nothing
                push!(eq.vcCandidates[5], v)
            elseif fixed
                push!(eq.vcCandidates[6], v)
            else
                push!(eq.vcCandidates[4], v)
            end
        end
    end

    if eq.log && eq.logDetails
        println("    eSolvedFixed = ", eq.der_eSolvedLower)
        println("    vSolvedFixed = ", eq.der_vSolvedLower)
        println("    vcCandidates = ", eq.vcCandidates)
    end

    # Tear equations (ec are non-differentiated, original equations)
    (eq.eSolved, eq.vSolved, eq.eResiduals, eq.vTear) =
        tearEquations!(eq.td, eq.fc.isSolvableEquation, eq.ec, eq.vcCandidates;
                       eSolvedFixed = eq.der_eSolvedLower,
                       vSolvedFixed = eq.der_vSolvedLower,
                       log = eq.log && eq.logDetails)
    return nothing
end

"""
    v = derivativeOfVlower(vLower::Vector{Int}, A::Vector{Int})

Given a set of variables vLower and the A-Vector of Pantelides,
return the derivative of vLower (so the differentiated variables of vLower).
"""
function derivativeOfVlower(vLower::Vector{Int}, A::Vector{Int})::Vector{Int}
    v = fill(0, length(vLower))
    for i in eachindex(vLower)
        derv = A[ vLower[i] ]
        @assert(derv > 0)
        v[i] = derv
    end
    return v
end

"""
    addFixedVariable!(eq::EquationGraph, v::Int)

If variable v is not differentiated AND has fixed=true store it in eq.equationInfo
(this is used to check during initialization, that the start value is respected)
"""
function addFixedVariable!(eq::EquationGraph, v::Int)::Nothing
    (v_startOrInit,v_fixed) = eq.fc.var_startInitFixed(v)
    if eq.Arev[v] == 0 && v_fixed
        push!(eq.vSolvedWithFixedTrue, v)
        push!(eq.equationInfo.vSolvedWithFixedTrue, eq.fc.var_name(v))
    end
    return nothing
end

"""
    addSolvedEquations!(eq::EquationGraph, eSolved::Vector{Int}, vSolved::Vector{Int})

Solve equation eSolved[i] for variable vSolved[i] and push this information to
eq.SortedEquations. Furthermore, generate the AST.
"""
function addSolvedEquations!(eq::EquationGraph, eSolved::Vector{Int}, vSolved::Vector{Int})::Nothing
    if eq.log
        println("Julia code:")
    end

    for i in eachindex(eSolved)
        # Push solved equation to eAST
        v = vSolved[i]
        e = eSolved[i]
        eq.eAST[e] = eq.fc.getSolvedEquationAST(e, v)   # filter_lineno( )
        if eq.log
            printEquations(eq.eAST[e])
        end
        eq.fullAssignRev[e] = v
        addFixedVariable!(eq,v)
    end

    return nothing
end

"""
    addLinearEquations!(eq::EquationGraph, hasConstantCoefficients::Bool, unitless::Bool)

A linear equation system is solved. Push this information to
eq.SortedEquations. Furthermore, generate the AST for a for-loop that builds
and solves a linear equation system A*x = b
from the solved equations and the residual equations,
and push this information to eq.SortedEquations.AST.
"""
function addLinearEquations!(eq::EquationGraph, hasConstantCoefficients::Bool, unitless::Bool)::Nothing
    # Construct body of for-loop
    empty!(eq.AST_aux)
    while_body = eq.AST_aux
    vTear_names   = String[]
    vTear_lengths = Int[]

    # Assign iteration variables
    #   unitless = false:
    #     v_i = leq.x[i]
    #
    #   unitless = true:
    #     v_i = leq.x[i]*@u_str($v_unit)
    #
    #   assign vector-valued variables vec_j at the end (storage is allocated inside LinearEquations)
    #     vec_j = leq.x_vec[j]
    #
    vAssigned_names   = Any[]
    v_vec_julia_names = Any[]
    v_vec             = Int[]
    v_vec_startOrInit = Any[]
    i1 = 0
    i2 = 0
    for (i,v) in enumerate(eq.vTear)
        (v_startOrInit, v_fixed) = eq.fc.var_startInitFixed(undifferentiated(eq,v))
        v_name = eq.fc.var_name(v)

        if isFixedLengthStartOrInit(v_startOrInit, v_name)
            # length(v_startOrInit) is fixed after compilation
            v_julia_name = eq.fc.var_julia_name(v)
            v_unit       = unitless ? "" : eq.fc.var_unit(v)
            v_length     = isnothing(v_startOrInit) ? 1 : length(v_startOrInit)
            i1 = i2 + 1
            i2 = i1 + v_length - 1
            if v_startOrInit isa AbstractVector
                # v is a vector
                x_elements = Expr[]
                for i in i1:i2
                    push!(x_elements, :( _leq_mode.x[$i] ))
                end
                if v_unit == ""
                    push!(while_body, :( $v_julia_name = Modia.SVector{$v_length,_FloatType}($(x_elements...))) )
                else
                    push!(while_body, :( $v_julia_name = Modia.SVector{$v_length,_FloatType}($(x_elements...))*@u_str($v_unit)) )
                end
            elseif v_startOrInit isa Number || v_startOrInit isa Nothing
                # v is a scalar or nothing (= assumed to be a scalar)
                if v_unit == ""
                    push!(while_body, :( $v_julia_name = _leq_mode.x[$i1] ) )
                else
                    push!(while_body, :( $v_julia_name = _leq_mode.x[$i1]*@u_str($v_unit) ) )
                end
            else
                error("Should not occur: v_name = $v_name, v_startOrInit = $v_startOrInit")
            end
            push!(vAssigned_names, v_julia_name)
            push!(vTear_lengths  , v_length)

            # Store iteration variables in linearEquations to be used in error messages
            push!(vTear_names, v_name)
        else
            # length(v_startOrInit) may change after compilation
            push!(v_vec, v)
            push!(v_vec_startOrInit, v_startOrInit)
        end
    end

    # Assignment code for vector-valued variables
    nx_fixedLength = length(vAssigned_names)
    for (i,v) in enumerate(v_vec)
        v_startOrInit = v_vec_startOrInit[i]
        v_julia_name  = eq.fc.var_julia_name(v)
        v_unit        = unitless ? "" : eq.fc.var_unit(v)
        if v_unit == ""
            push!(while_body, :( $v_julia_name = _leq_mode.x_vec[$i] ) )
        else
            push!(while_body, :( $v_julia_name = _leq_mode.x_vec[$i]*@u_str($v_unit) ))
        end
        push!(vAssigned_names  , v_julia_name)
        push!(vTear_lengths    , length(v_startOrInit))
        push!(v_vec_julia_names, v_julia_name)

        # Store iteration variables in linearEquations to be used in error messages
        push!(vTear_names, eq.fc.var_name(v))
    end

    # Add solved equations
    for (i,e) in enumerate(eq.eSolved)
        push!(while_body, eq.fc.getSolvedEquationAST(e, eq.vSolved[i]))
        push!(vAssigned_names, eq.fc.var_julia_name(eq.vSolved[i]))
    end

    # Add residual equations
    #   appendResidual!(leq.residuals, < equation in residual form >)
    for (i,e) in enumerate(eq.eResiduals)
        e_AST = eq.fc.getResidualEquationAST(e, :(_leq_mode.residuals[$i]) )
        if !isnothing(e_AST)
            push!(while_body, e_AST)
        end
    end

    # Generate LinearEquations data structure
    push!(eq.equationInfo.linearEquations, (vTear_names, v_vec_julia_names, vTear_lengths, nx_fixedLength, hasConstantCoefficients))

    # Construct for-loop
    leq_index = length(eq.equationInfo.linearEquations)
    while_loop = quote
        local $(vAssigned_names...)
        _leq_mode = initLinearEquationsIteration!(_m, $leq_index)
        Modia.TimerOutputs.@timeit _m.timer "Modia LinearEquationsIteration!" while Modia.LinearEquationsIteration!(_leq_mode, _m.isInitial, _m.solve_leq, _m.storeResult, _m.time, _m.timer)
            $(while_body...)
        end
        _leq_mode = nothing
    end

    if eq.log
        showCodeWithoutComments(while_loop)
    end

    # Store information about linear equation system
    teq = TearedEquations(true, hasConstantCoefficients,
                          copy(eq.vSolved), copy(eq.vTear),
                          copy(eq.eSolved), copy(eq.eResiduals),
                          while_loop)
    push!(eq.tearedEquations, teq)
    teq_index = length(eq.tearedEquations)

    # Reference teq_index for all equations
    for eSolved in eq.eSolved
        eq.tearedEquations_indices[eSolved] = teq_index
    end
    for eResiduals in eq.eResiduals
        eq.tearedEquations_indices[eResiduals] = teq_index
    end

    # Store fixed=true variables
    for vTear in eq.vTear
        addFixedVariable!(eq,vTear)
    end
    for vSolved in eq.vSolved
        addFixedVariable!(eq,vSolved)
    end

    # Store vTear without a start value
    for vTear in eq.vTear
        (vTear_startOrInit, vTear_fixed) = eq.fc.var_startInitFixed(vTear)
        if vTear_startOrInit isa Nothing
            push!(eq.vTearWithoutStart, vTear)
        end
    end

    return nothing
end

"""
    sortEquations!(eq:EquationGraph)

Generate the sorted equations
"""
function sortEquations!(eq::EquationGraph, Goriginal)::Nothing
    if eq.log
         println("Sort equations (BLT on all equations under the assumption that the ODE states are known).")
    #    println("\nSorted equations (BLT on all equations under\nthe assumption that the ODE states are known):")
    #
    #    code = quote
    #        $(eq.AST...)
    #    end
    #    showCodeWithoutComments(code)
    end

    # Matching of all equations
    assign    = matching(Goriginal, length(eq.A), eq.vActive)
    assignRev = revertAssociation(assign)

    # Sort equations and find strong components
    blt = BLT(Goriginal,assign)

    # Generate AST of sorted equations
    for blt_i in blt
        @assert(length(blt_i) >= 1)

        solveLinearEquation = true
        if length(blt_i) == 1 # One equation in one unknown
            e = blt_i[1]

            # Check
            v = eq.fullAssignRev[e]
            if v != -1
                @assert(e == assign[v])

                # Push AST
                push!(eq.AST, eq.eAST[e])
                solveLinearEquation = false
            end
        end

        if solveLinearEquation # Equation system
            # Check that all equations reference the same equation system
            e1 = blt_i[1]
            teq_index = eq.tearedEquations_indices[e1]
            @assert(teq_index != 0)
            for j=2:length(blt_i)
                e = blt_i[j]
                if eq.tearedEquations_indices[e] != teq_index
                    error("... should not occur: teq_index=$teq_index, e1=$e1, e_teq_index=", eq.tearedEquations_indices[e])
                end
            end

            # Check that the assigned variables of the BLT block are the variables computed from the equation system
            teq = eq.tearedEquations[teq_index]
            vTearComputed = union(teq.vTear, teq.vSolved)
            vBltComputed = Int[]
            for e in blt_i
                push!(vBltComputed, assignRev[e])
            end
            vdiff = setdiff(vBltComputed, vTearComputed)
            @assert(length(vdiff) == 0)

            # Push AST
            push!(eq.AST, teq.AST)
        end
    end

    return nothing
end

"""
    (success, AST, equationInfo) = getSortedAndSolvedAST(
                                        G, BLT, assign, A, B, stateSelectionFunctions;
                                        log = false, logDetails=false, modelName = "???",
                                        defaultParameterAndStartValues = nothing,
                                        unitless=false)

From the BLT on highest derivative level `(G, BLT, assign, A, B)`,
a set of callback functions `stateSelectionFunctions`, this function returns the
AST (Abstract Syntax Tree) of the sorted and solved equations as
a vector of `Expr` and information about the ODE equations
(especially names of the ODE states).

# Input arguments

- `G`: G[i] is the Integer vector of variable indices of equation i.
   A variable i can be a scalar or an array.
   is of type Vector{ Vector{Int} } or Vector{ Any }.

- `BLT::Vector{Vector{Int}}`: BLT[i] is the vector of equations belonging to BLT-block i.
  In principal, **BLT[i] must include only highest derivative equations**.
  In order to ease the generation of BLT[i], to not be forced to map indices forth and back
  when calling the function that computes the BLT, the following assumption is
  made: If BLT[j] contains only **one equation with B[j][1] > 0** then this BLT block is
  a lower-derivative (dummy) block and is ignored.

- `assign::Vector{Int}`: ei = assign[vi] is equation ei assigned to variable vi

- `A`: A-Vector of Pantelides algorithm:
       `A[i] = if der(v[i]) == v[k] then k else 0`
       where `v[i]` is variable `i`.

- `B`: B-Vector of Pantelides algorithm:
       `B[i] = if der(e[i]) == e[k] then k else 0`
       where `e[i]` is equation `i`.

- `stateSelectionFunctions`:
  An instance of immutable struct [`StateSelectionFunctions`](@ref)
  that holds the information about the callback functions needed by
  `getSortedAndSolvedAST`.

- `log`: = true: Print debug information (to find bugs in the code).

- `logDetails`: = true: Print detailed debug information, if `log = true`.

- `modelName`: Name of model used in messages (otherwise not used).

# Output arguments

A tuple with the following values:

- `success::Bool`: = true, if generation was successful. = false, if
  an error occured. Arguments `AST, equationInfo` hold the information that
  have been collected until the error was detected.

- `AST::Vector{Expr}`: The AST (Abstract Syntax Tree) of the sorted and solved equations
  as Vector of `Expr`.

- `equationInfo::`[`Modia.EquationInfo`](@ref): Object that defines
  the states of the ODE.

When variable `v` has neither a `start` nor `init` attribute defined,
it is assumed to have `start=0.0`.

| defined   | Description                                                |
|:----------|:-----------------------------------------------------------|
| `start`   | `start` is allowed to be changed during initialization     |
| `init`    | `init` is not allowed to be changed during initialization  |

Note, `start, init` influence the state selection.

If an ODE state has neither a `start` nor an `init` value defined, a warning message is printed
(start/init value missing).

If a solved variable has an `init` value, an information message is printed in case `log=true`
(init value has no effect).

# Sketch of the Algorithm

Note, a statement *error xxx* means that an error occurs,
because it is not possible to transform the equations to an ODE.

## Generate equation sets

From BLT blocks that are not further differentiated (= highest order BLT blocks),
generate the equation sets as described by the dummy derivative method of
[Mattsson and Söderlind (1992)](https://ieeexplore.ieee.org/document/274429).
See also sections 4.5 and 4.6 of the paper
[Otter and Elmqvist (2017)](https://modelica.org/events/modelica2017/proceedings/html/submissions/ecp17132565_OtterElmqvist.pdf).

## Define initial ODE states

Initially, all differentiated variables are defined to be ODE states.
Whenever an initial ODE state is explicitly solved for in the course of the
algorithm, the variable is no longer an ODE state.

## Analyse equation sets

For every equation set on every differentiation level perform the following actions:

- If the equation set consists of *one* equation in *one* unknown, the equation
  is solved for this unknown (*error, if this is not possible*).

- If the equation set consists of *N linear* equations in *N* unknowns solve this
  equation system (*error, if no linear system*):\\
  The equations are first teared to reduce the number of iteration variables and
  afterwards the teared equation system is solved with a special iterator loop that
  solves a linear equation system with an LU decomposition with column pivoting
  (for details see [`Modia.LinearEquationsIteration!`](@ref)).\\

- If the equation set consists of *N linear* equations in *M* unknowns (*M > N*) perform
  the following actions (*error, if no linear system*).\\
  Note, due to the structure of BLT, this equation set cannot be
  on highest derivative level; due to the Pantelides algorithm,
  M < N is not possible).

  1. Solve the equation system explicitly via tearing for N unknowns
     (*error, if this is not possible*). Use start/init attributes of the M unknowns
     to influence tearing (e.g. it is tried to eliminate first variables that have neither start not init values).

  2. The N solved variables have been initially ODE states and are now defined to be
     no ODE states (so called dummy states, that is algebraic variables).

A BLT transformation of **all equations** is made under the assumption that
the selected ODE states are known. This phase is needed, because the dummy-derivative
ordering at the beginning does not necessarily provide already the right ordering.
In this phase the information from the previous phase is used (that is the already determined
tearing variables for systems of equations are utilized).
"""
function getSortedAndSolvedAST(Goriginal,     # Typically ::Vector{Vector{Int}}
                               G,
                               BLT,   # Typically ::Vector{Vector{Int}}
                               assign::Vector{Int},
                               A::Vector{Int},
                               B::Vector{Int},
                               stateSelectionFunctions::StateSelectionFunctions;
                               log::Bool        = false,
                               logDetails::Bool = false,
                               logStates::Bool  = false,
                               modelName        = "???",
                               unitless::Bool   = false,
                               defaultParameterAndStartValues = nothing)

    if log
        println("\n=== getSortedAndSolvedAST(...) started for $modelName.")
    end

    # Define error/warning function
    showMessage2(message; severity=0, details="",variables=Int[],equations=Int[]) =
    stateSelectionFunctions.showMessage(message; from="getSortedAndSolvedAST for model $modelName",
                                            severity=severity, details=details, variables=variables, equations=equations)

    # Construct data structure eq that holds all needed information for the state selection
    # (this includes constructing the equation/constraint sets and initialization of tearing)
    eq = EquationGraph(G, BLT, assign, A, B, stateSelectionFunctions, log, logDetails, showMessage2)
    if !eq.success
        @goto ERROR_RETURN
    end
    eq.equationInfo.defaultParameterAndStartValues = defaultParameterAndStartValues

    # Inspect every equation/constraint set in sequence from lowest-order to highest-order derivative
    for j in eachindex(eq.eConstraintsVec)
        # Equation/Constraint set and its unknowns that are inspected
        eConstraints = eq.eConstraintsVec[j]
        vConstraints = eq.vConstraintsVec[j]

        # Analyze all equation sets eConstraints[i] from lowest-order to highest-order derivatives
        for i in eachindex(eConstraints)
            if log
                println("\n... Equation set $j.$i ..............................")
                println("Equations: ")
                printEquations(eq, eConstraints[i])
                println("Unknown variables: ")
                printVariables(eq, vConstraints[i])
            end

            solveLinearEquation = true
            if length(eConstraints[i]) == 1 && length(vConstraints[i]) == 1
                # One equation in one unknown
                if i < length(eConstraints)
                    # Not on highest derivative level. Unknown is a dummy state.
                    v = vConstraints[i][1]
                    @assert(!eq.vActive[v])
                    eq.vActive[v] = true
                    if log
                        println("One equation in one unknown variable: Unknown is a dummy state.")
                    end
                else
                    if log
                        println("One equation in one unknown variable. Solve the equation:")
                    end
                end

                # Solve equation for unknown v and add solved equation to AST
                try
                    addSolvedEquations!(eq, eConstraints[i], vConstraints[i])
                    solveLinearEquation = false
                catch
                    if log
                        println("Not possible to solve the equation directly. Try so solve it as linear equation:")
                    end
                    solveLinearEquation = true
                    eq.fullAssignRev[eConstraints[i][1]] = -1
                end
            end

            if solveLinearEquation
                # N equations in M unknowns (M >= N)
                @assert(length(vConstraints[i]) >= length(eConstraints[i]))
                if log
                    N = length(eConstraints[i])
                    M = length(vConstraints[i])
                    println("$N equation(s) in $M unknown variable(s). Tear the system of equations:")
                end

                # Tear systems of equations
                if i == 1
                    # On the lowest derivative level, initialize tearing with empty vectors
                    empty!(eq.der_eSolvedLower)
                    empty!(eq.der_vSolvedLower)
                    empty!(eq.der_vTearLower)

                    # Variables/equations to be inspected by tearing
                    eq.vc = vConstraints[i]
                    eq.ec = eConstraints[i]

                else
                    # On derivative level i > 1, initialize tearing with the differentiated eSolved, vSolved, vTear
                    # from the previous level i-1.
                    # (eConstraints[i] is the derivative of eConstraints[i-1] + potentially additional equations)
                    eq.der_eSolvedLower = derivativeOfVlower(eq.eSolved, B)
                    eq.der_vSolvedLower = derivativeOfVlower(eq.vSolved, A)
                    eq.der_vTearLower   = derivativeOfVlower(eq.vTear  , A)

                    # Variables/equations to be inspected by tearing (remove the parts that are known from level i-1)
                    eq.vc = setdiff(vConstraints[i], eq.der_vSolvedLower)
                    eq.ec = setdiff(eConstraints[i], eq.der_eSolvedLower)
                end

                if log && i < length(eConstraints)
                    vWithStart = Int[]
                    for v in vConstraints[i]
                        (v_startOrInit,v_fixed) = eq.fc.var_startInitFixed(v)
                        if !(v_startOrInit isa Nothing)
                            push!(vWithStart, v)
                        end
                    end
                    println("    Unknowns with start or init values: ", getNames(eq, vWithStart))
                end

                tearEquationsWithCandidates!(eq)
                if log
                    println("    Tearing  variables: ", getNames(eq, eq.vTear))
                    if length(eq.eResiduals) > 0
                        println("    Residual equations: ");
                        printEquations(eq, eq.eResiduals; indent=8)
                    end
                end

                # Select dummy states
                if i < length(eConstraints)
                    # Not on highest derivative level. Solved variables are dummy states
                    for v in eq.vSolved
                        @assert(!eq.vActive[v])
                        eq.vActive[v] = true
                    end

                    if length(eConstraints[i]) == length(vConstraints[i])
                        # N equations in N unknowns - tearing variables are dummy states
                        for v in eq.vTear
                            @assert(!eq.vActive[v])
                            eq.vActive[v] = true
                        end
                        if log
                            println("    All unknowns are dummy states.")
                        end

                    else
                        if log
                            println("    All solved unknowns are dummy states.")
                        end
                    end
                elseif log
                    println("    All unknowns are solved.")
                end

                # Further actions depending on type of equations
                if length(eq.eResiduals) == 0
                    # There are no residual equations, add solved equations to AST
                    addSolvedEquations!(eq, eq.eSolved, eq.vSolved)

                elseif length(eConstraints[i]) == length(vConstraints[i])
                    # N equations in N unknowns with residual equations

                    # Check that equation system is linear in the unknowns
                    (isLinear, hasConstantCoefficients) = isLinearEquationSystem!(eq, eConstraints[i], vConstraints[i])
                    if !isLinear
                        if i == length(eConstraints)
                            # On highest derivative level:
                            # Assume that the equation system is linear, if at least one of the unknowns is a derivative
                            linearAssumption = true  # temporarily

                            #linearAssumption = false
                            #for v in vConstraints[i]
                            #    if eq.Arev[v] > 0
                            #        linearAssumption = true
                            #    end
                            #end

                            if linearAssumption
                                isLinear = true
                                hasConstantCoefficients = false
                                #showMessage2("It is heuristically assumed that equation system is linear (although isLinearEquation returned isLinear=false).";
                                #             severity  = WARNING,
                                #             variables = vConstraints[i],
                                #             equations = eConstraints[i])
                            end
                        end

                        if !isLinear
                            showMessage2("Cannot transform to ODE, because equation system is not linear.";
                                        severity  = ERROR,
                                        variables = vConstraints[i],
                                        equations = eConstraints[i])
                            @goto ERROR_RETURN
                        end
                    end

                    if log
                        println("Teared equation system is linear. Solve system with hasConstantCoefficients = $hasConstantCoefficients.")
                    end

                    # Generate AST to build-up and solve linear equation system of the teared equations
                    addLinearEquations!(eq, hasConstantCoefficients, unitless)

                else
                    # N equations in M unknowns with M > N and residual equations
                    @assert(length(vConstraints[i]) > length(eConstraints[i]))

                    # Assertion, that this case is not possible on highest derivative level
                    @assert(i < length(eConstraints))

                    showMessage2("Cannot transform to ODE because constraint equations cannot be explicitly solved by tearing",
                                 severity  = ERROR,
                                 variables = vConstraints[i],
                                 equations = eConstraints[i])
                    @goto ERROR_RETURN
                end
            end # if
        end  # for
    end # for

    # Determine ODE states
    ODE_states = Int[]
    for (ve,veActive) in enumerate(eq.vActive)
        if !veActive
            push!(ODE_states, ve)
        end
    end

    # Store state information in equationInfo
    x_info            = eq.equationInfo.x_info
    x_without_start   = Int[]
    x_vec             = Int[]
    x_vec_startOrInit = Any[]
    x_vec_fixed       = Bool[]
    for v in ODE_states
        (v_startOrInit, v_fixed) = eq.fc.var_startInitFixed(v)
        if isnothing(v_startOrInit)
            push!(x_without_start, v)
        end
        v_name = eq.fc.var_name(v)

        if isFixedLengthStartOrInit(v_startOrInit, v_name)
            v_der_x            = eq.A[v]
            v_unit             = unitless ? "" : eq.fc.var_unit(v)
            v_julia_name       = eq.fc.var_julia_name(v)
            v_der_x_name       = eq.fc.var_name(v_der_x)
            v_der_x_julia_name = eq.fc.var_julia_name(v_der_x)
            v_nominal          = eq.fc.var_nominal(v)
            v_unbounded        = eq.fc.var_unbounded(v)
            v_stateCategory    = Modia.XD

            push!(x_info, Modia.StateElementInfo(
                            v_name, v_julia_name,
                            v_der_x_name, v_der_x_julia_name,
                            v_stateCategory, v_unit, v_startOrInit, v_fixed,
                            v_nominal, v_unbounded))
        else
            # length(v_startOrInit) may change after compilation
            push!(x_vec            , v)
            push!(x_vec_startOrInit, v_startOrInit)
            push!(x_vec_fixed      , v_fixed)
        end
    end
    eq.equationInfo.nxFixedLength = length(x_info)
    for (i,v) in enumerate(x_vec)
        v_startOrInit      = x_vec_startOrInit[i]
        v_fixed            = x_vec_fixed[i]
        v_name             = eq.fc.var_name(v)
        v_der_x            = eq.A[v]
        v_unit             = unitless ? "" : eq.fc.var_unit(v)
        v_julia_name       = eq.fc.var_julia_name(v)
        v_der_x_name       = eq.fc.var_name(v_der_x)
        v_der_x_julia_name = eq.fc.var_julia_name(v_der_x)
        v_nominal          = eq.fc.var_nominal(v)
        v_unbounded        = eq.fc.var_unbounded(v)
        v_stateCategory    = Modia.XD

        push!(x_info, Modia.StateElementInfo(
                        v_name, v_julia_name,
                        v_der_x_name, v_der_x_julia_name,
                        v_stateCategory, v_unit, v_startOrInit, v_fixed,
                        v_nominal, v_unbounded))
    end

    # Handle systems with only algebraic variables, by introducing a dummy
    # differential equation der_x[1] = -x[1].
    if length(ODE_states) == 0
        if log
            println("Model has only algebraic variables.\n",
                    "Added a dummy differential equation der(_dummy_x) = -_dummy_x, _dummy_x(t0) = 0")
        end
        push!(eq.equationInfo.x_info, Modia.StateElementInfo(
              "_dummy_x", :(), "der(_dummy_x)", :(), XD, "", 0.0, true, NaN, false))
    end

    # Finalize equationInfo
    initEquationInfo!(eq.equationInfo)

    # Print ODE states
    if logStates
        println("\nSelected ODE states: ")
        x_table = Modia.get_x_table(x_info)
        show(stdout, x_table; allrows=true, allcols=true, summary=false, eltypes=false)
        print("\n\n")
    end

    # Sort equations and generate sorted AST (eq.AST)
    sortEquations!(eq, Goriginal)

    # Check that all ODE states have a start value
    if length(x_without_start) > 0
        showMessage2("Init/start values missing in the model for some ODE states.";
                     severity  = WARNING,
                     variables = x_without_start)
    end

    # Print info, if there are iteration variables without a start value
    if length(eq.vTearWithoutStart) > 0 && !unitless
        showMessage2("The following variables are iteration variables but have no start/init values defined.",
                     details = "If units are used in the model, start/init values with correct units should be defined\n" *
                               "to avoid unit errors during compilation.",
                     severity  = INFORMATION,
                     variables = eq.vTearWithoutStart)
    end

    # Print warning, if there are variables with fixed=true that are
    # explicitly solved for and log=true
    if log && length(eq.equationInfo.vSolvedWithFixedTrue) > 0
        showMessage2("The following variables have an 'init' initialization and are explicitly solved for.",
                     details = "Therefore, the 'init' values have no effect, but must exactly match the values,\n"*
                               "computed during initialization. Otherwise this gives a run-time error.\n"*
                               "It is adviced to use 'start' initialization or remove initialization for these variables." ,
                     severity  = WARNING,
                     variables = eq.vSolvedWithFixedTrue)
    end

    if log
        println("\n=== getSortedAndSolvedAST(...) terminated for $modelName.\n")
    end
    return (true, eq.AST, eq.equationInfo)

    @label ERROR_RETURN
    return (false, eq.AST, eq.equationInfo)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "libsvm"
version = v"3.25.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/cjlin1/libsvm/archive/v$(version.major)$(version.minor).tar.gz",
                  "1f587ec0df6fd422dfe50f942f8836ac179b0723b768fe9d2fabdfd1601a0963"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libsvm*/
for f in ${WORKSPACE}/srcdir/patches/*.patch; do
    atomic_patch -p1 ${f}
done
make
make lib
install -Dvm 0755 "libsvm.${dlext}" "${libdir}/libsvm.${dlext}"
install -Dvm 0755 "svm-train${exeext}" "${bindir}/svm-train${exeext}"
install -Dvm 0755 "svm-predict${exeext}" "${bindir}/svm-predict${exeext}"
install -Dvm 0755 "svm-scale${exeext}" "${bindir}/svm-scale${exeext}"
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libsvm", :libsvm),
    ExecutableProduct("svm-scale", :svm_scale),
    ExecutableProduct("svm-train", :svm_train),
    ExecutableProduct("svm-predict", :svm_predict)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    # For OpenMP we use libomp from `LLVMOpenMP_jll` where we use LLVM as compiler (BSD
    # systems), and libgomp from `CompilerSupportLibraries_jll` everywhere else.
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"); platforms=filter(!Sys.isbsd, platforms)),
    Dependency(PackageSpec(name="LLVMOpenMP_jll", uuid="1d63c593-3942-5779-bab2-d838dc0a180e"); platforms=filter(Sys.isbsd, platforms)),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

ext: jl

max_stars_repo_name: cmhamel/Yggdrasil

max_issues_repo_name: cmhamel/Yggdrasil

max_forks_repo_name: cmhamel/Yggdrasil


text:
module Solution
  export simpleMultiplication
  simpleMultiplication(number) = number * ( 8 + (number % 2) )
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
module PyFenicsAD

import ChainRulesCore
using PyCall

const fecr = PyNULL()
const to_numpy = PyNULL()
const from_numpy = PyNULL()
const evaluate_primal = PyNULL()
const evaluate_pullback = PyNULL()
const evaluate_pushforward = PyNULL()

function __init__()
    copy!(fecr, pyimport("fecr"))
    copy!(to_numpy, fecr.to_numpy)
    copy!(from_numpy, fecr.from_numpy)
    copy!(evaluate_primal, fecr.evaluate_primal)
    copy!(evaluate_pullback, fecr.evaluate_pullback)
    copy!(evaluate_pushforward, fecr.evaluate_pushforward)
end

export fecr, to_numpy, from_numpy, evaluate_primal, evaluate_pullback, evaluate_pushforward

"""
A convenience macro that creates a julia function and registers forward or reverse rules for the provided FEniCS/Firedrake `fem_function`.
"""
macro register_fem_function(julia_function_name, fem_templates, fem_function)
    name_str = "$julia_function_name"
    quote
        function $(esc(julia_function_name))(inputs...)
            return evaluate_primal($(esc(fem_function)), $(esc(fem_templates)), inputs...)[1]
        end

        function ChainRulesCore.rrule(::typeof($(esc(julia_function_name))), inputs...)
            pyout = pycall(
                evaluate_primal,
                PyObject,
                $(esc(fem_function)),
                $(esc(fem_templates)),
                inputs...,
            )
            numpy_output, fem_output, fem_inputs, tape =
                [get(pyout, PyObject, i) for i = 0:3]

            # Only single-output functions are supported for now
            function fem_pullback(g)
                vjp_out = evaluate_pullback(fem_output, fem_inputs, tape, g)
                return (ChainRulesCore.NO_FIELDS, vjp_out...)
            end
            return get(pyout, 0), fem_pullback
        end

        function ChainRulesCore.frule(
            (_, Δinputs),
            ::typeof($(esc(julia_function_name))),
            inputs...,
        )
            pyout = pycall(
                evaluate_primal,
                PyObject,
                $(esc(fem_function)),
                $(esc(fem_templates)),
                inputs...,
            )
            numpy_output, fem_output, fem_inputs, tape =
                [get(pyout, PyObject, i) for i = 0:3]

            ∂numpy_output = evaluate_pushforward(fem_output, fem_inputs, tape, Δinputs)
            return get(pyout, 0), ∂numpy_output
        end

        println("Registered FEniCS/Firedrake function with name ", $name_str, "!")
    end
end

export @register_fem_function

end # module


text:
module DataGraphs
using Dictionaries
using Graphs
using MultiDimDictionaries
using NamedGraphs
using SimpleTraits

using MultiDimDictionaries: tuple_convert, SliceIndex, ElementIndex

#
# imports
#

import Base:
  get, getindex, setindex!, convert, show, isassigned, eltype, copy, hvncat, hcat, vcat
import Graphs:
  adjacency_matrix,
  add_edge!,
  add_vertex!,
  bfs_parents,
  bfs_tree,
  dfs_parents,
  dfs_tree,
  edges,
  edgetype,
  has_edge,
  has_vertex,
  induced_subgraph,
  is_connected,
  is_cyclic,
  is_directed,
  is_strongly_connected,
  is_weakly_connected,
  ne,
  neighbors,
  nv,
  rem_edge!,
  rem_vertex!,
  vertices
import MultiDimDictionaries: IndexType
import NamedGraphs: disjoint_union, ⊔, to_vertex

# Dictionaries.jl patch
# TODO: delete once fixed in Dictionaries.jl
convert(::Type{Dictionary{I,T}}, dict::Dictionary{I,T}) where {I,T} = dict

# General functions
not_implemented() = error("Not implemented")

# Returns just the edges of a directed graph,
# but both edge directions of an undirected graph.
@traitfn function all_edges(g::::IsDirected)
  return edges(g)
end

@traitfn function all_edges(g::::(!IsDirected))
  e = edges(g)
  return Iterators.flatten(zip(e, reverse.(e)))
end

include("abstractdatagraph.jl")
include("datagraph.jl")
include("nameddimdatagraph.jl")

#
# exports
#

export DataGraph,
  NamedDimDataGraph,
  AbstractNamedDimDataGraph,
  AbstractDataGraph,
  map_vertex_data,
  map_edge_data,
  map_data,
  disjoint_union,
  ⊔

end # module DataGraphs


text:
"""
It iterates (approximately) equilateral triangles!
"""
struct EquilateralIterator{T}
    dists       ::Matrix{T}
    sortededges ::Vector{Tuple{Int, Int}}
    tol         ::T
end

Base.show(io::IO, ::EquilateralIterator{T}) where {T} =
    print(io, "EquilateralIterator{$T}(...)")

Base.IteratorSize(::EquilateralIterator) = Base.SizeUnknown()
Base.IteratorEltype(::EquilateralIterator) = Base.HasEltype()
Base.eltype(::EquilateralIterator{T}) where {T} = Tuple{Tuple{Int, Int, Int}, T}
Base.eltype(::Type{EquilateralIterator{T}}) where {T} = Tuple{Tuple{Int, Int, Int}, T}

"""
    equilaterals(dists, tol)

Iterator that returns approximately (up to `tol`) equilateral triangles in
distance matrix `dists`, sorted by longest side length.
"""
equilaterals(dists::AbstractMatrix{T}, tol) where {T} =
    EquilateralIterator{T}(dists, sortededges(dists), T(tol))

"""
    sortededges(dists)

Return edges `(i, j)` ordered by the values of `dists[i, j]`.
"""
sortededges(dists) =
    sort!([(i, j) for i in 2:size(dists, 1) for j in 1:i-1],
          lt = (e1, e2) -> isless(dists[e1...], dists[e2...]))

function Base.iterate(ei::EquilateralIterator, st = (1, 1))
    dists = ei.dists
    n = size(dists, 1)
    sortededges = ei.sortededges
    tol = ei.tol

    eindex, kstart = st
    if kstart > n
        kstart = 1
        eindex += 1
    end

    @inbounds while eindex ≤ length(sortededges)
        i, j = sortededges[eindex]
        a = dists[i, j]

        for k in kstart:n
            (k == i || k == j) && continue
            b = dists[k, i]
            c = dists[k, j]
            # This makes sure each triangle is returned only once when some of
            # the sides have the same length.
            a > b || k > i || continue
            a > c || k > j || continue

            if 0 ≤ a - b ≤ tol && 0 ≤ a - c ≤ tol && abs(b - c) ≤ tol
                return ((i, j, k), a), (eindex, k + 1)
            end
        end
        kstart = 1
        eindex += 1
    end
    nothing
end

max_stars_repo_name: mtsch/IntrinsicPersistence

max_issues_repo_name: mtsch/IntrinsicPersistence

max_forks_repo_name: mtsch/IntrinsicPersistence


text:
using Plots, Test

@testset "Plotly" begin
    @testset "Basic" begin
        @test plotly() == Plots.PlotlyBackend()
        @test backend() == Plots.PlotlyBackend()

        p = plot(rand(10))
        @test isa(p, Plots.Plot) == true
        @test_nowarn Plots.plotly_series(plot())
    end

    @testset "Contours" begin
        x = (-2π):0.1:(2π)
        y = (-π):0.1:π
        z = cos.(y) .* sin.(x')

        @testset "Contour numbers" begin
            @testset "Default" begin
                @test Plots.plotly_series(contour(x, y, z))[1][:ncontours] ==
                      Plots._series_defaults[:levels] + 2
            end
            @testset "Specified number" begin
                @test Plots.plotly_series(contour(x, y, z, levels = 10))[1][:ncontours] ==
                      12
            end
        end

        @testset "Contour values" begin
            @testset "Range" begin
                levels = -1:0.5:1
                p = contour(x, y, z, levels = levels)
                @test p[1][1].plotattributes[:levels] == levels
                @test Plots.plotly_series(p)[1][:contours][:start] == first(levels)
                @test Plots.plotly_series(p)[1][:contours][:end] == last(levels)
                @test Plots.plotly_series(p)[1][:contours][:size] == step(levels)
            end

            @testset "Set of contours" begin
                levels = [-1, -0.25, 0, 0.25, 1]
                levels_range =
                    range(first(levels), stop = last(levels), length = length(levels))
                p = contour(x, y, z, levels = levels)
                @test p[1][1].plotattributes[:levels] == levels
                series_dict = @test_logs (
                    :warn,
                    "setting arbitrary contour levels with Plotly backend " *
                    "is not supported; use a range to set equally-spaced contours or an " *
                    "integer to set the approximate number of contours with the keyword " *
                    "`levels`. Setting levels to -1.0:0.5:1.0",
                ) Plots.plotly_series(p)
                @test series_dict[1][:contours][:start] == first(levels_range)
                @test series_dict[1][:contours][:end] == last(levels_range)
                @test series_dict[1][:contours][:size] == step(levels_range)
            end
        end
    end
end


text:
@testset "Partial derivatives of a solution w.r.t. to initial conditions and parameters" begin

    test_cases = []
    DType = Union{fmpq_mpoly, Generic.Frac{fmpq_mpoly}}

    R, (x, a) = PolynomialRing(QQ, ["x", "a"])
    push!(
        test_cases,
        Dict(
          "ODE" => ODE{fmpq_mpoly}(Dict{fmpq_mpoly, DType}(x => x + a), Array{fmpq_mpoly, 1}()),
          "ic" => Dict(x => QQ(rand(1:10))),
          "param_vals" => Dict(a => QQ(rand(1:10))),
          "inputs" => Dict{fmpq_mpoly, Array{fmpq, 1}}(),
          "prec" => 20
        )
    )

    R, (x, a) = PolynomialRing(QQ, ["x", "a"])
    push!(
        test_cases,
        Dict(
          "ODE" => ODE{fmpq_mpoly}(Dict{fmpq_mpoly, DType}(x => x^2 + a),Array{fmpq_mpoly, 1}()),
          "ic" => Dict(x => QQ(rand(1:10))),
          "param_vals" => Dict(a => QQ(rand(1:10))),
          "inputs" => Dict{fmpq_mpoly, Array{fmpq, 1}}(),
          "prec" => 20
        )
    )

    R, (x, y, a, b) = PolynomialRing(QQ, ["x", "y", "a", "b"])
    push!(
        test_cases,
        Dict(
          "ODE" => ODE{fmpq_mpoly}(Dict{fmpq_mpoly, DType}(x => x^2 + 2 * x * y - 3 * a * y, y =>  x^2 + a * b - b^2 + 4 * b * x), Array{fmpq_mpoly, 1}()),
          "ic" => Dict(x => QQ(rand(1:10)), y => QQ(rand(1:10))),
          "param_vals" => Dict(a => QQ(rand(1:10)), b => QQ(rand(1:10))),
          "inputs" => Dict{fmpq_mpoly, Array{fmpq, 1}}(),
          "prec" => 8
        )
    )

    R, (x, a, u) = PolynomialRing(QQ, ["x", "a", "u"])
    push!(
        test_cases,
        Dict(
          "ODE" => ODE{fmpq_mpoly}(Dict{fmpq_mpoly, DType}(x => u + a), [u]),
          "ic" => Dict(x => QQ(rand(1:10))),
          "param_vals" => Dict(a => QQ(rand(1:10))),
          "inputs" => Dict(u => [QQ(rand(-3:3)) for i in 1:20]),
          "prec" => 20
        )
    )

    F = GF(2^31 - 1)
    P = Generic.MPoly{Nemo.gfp_elem}
    DType = Union{P, Generic.Frac{P}}

    varnames = vcat(
        ["x_$i" for i in 1:3],
        ["p_$i" for i in 1:3],
        ["u_$i" for i in 1:2],
    )
    R, vars = PolynomialRing(F, varnames)
    push!(
        test_cases,
        Dict(
          "ODE" => ODE{P}(Dict{P, DType}(vars[i] => rand_poly(1, vars) for i in 1:3), vars[7:8]),
          "ic" => Dict(vars[i] => F(rand(1:50)) for i in 1:3),
          "param_vals" => Dict(vars[i + 3] => F(rand(1:50)) for i in 1:3),
          "inputs" => Dict(u => [F(rand(-30:30)) for i in 1:6] for u in vars[7:end]),
          "prec" => 6
        )
    )

    varnames = vcat(
        ["x_$i" for i in 1:3],
        ["p_$i" for i in 1:3],
        ["u_$i" for i in 1:2],
    )
    R, vars = PolynomialRing(F, varnames)
    push!(
        test_cases,
        Dict(
          "ODE" => ODE{P}(Dict{P, DType}(vars[i] => rand_poly(2, vars) for i in 1:3), vars[7:8]),
          "ic" => Dict(vars[i] => F(rand(1:50)) for i in 1:3),
          "param_vals" => Dict(vars[i + 3] => F(rand(1:50)) for i in 1:3),
          "inputs" => Dict(u => [F(rand(-30:30)) for i in 1:6] for u in vars[7:end]),
          "prec" => 6
        )
    )

    varnames = vcat(
        ["x_$i" for i in 1:2],
        ["p_$i" for i in 1:2],
        "u",
    )
    R, vars = PolynomialRing(F, varnames)
    push!(
        test_cases,
        Dict(
          "ODE" => ODE{P}(Dict{P, DType}(vars[i] => rand_poly(1, vars) // (vars[1] + vars[3]) for i in 1:2), [vars[end]]),
          "ic" => Dict(vars[i] => F(rand(1:50)) for i in 1:2),
          "param_vals" => Dict(vars[i + 2] => F(rand(1:50)) for i in 1:2),
          "inputs" => Dict(vars[end] => [F(rand(-30:30)) for i in 1:4]),
          "prec" => 4
        )
    )

    for case in test_cases
        ode, prec = case["ODE"], case["prec"]
        @time sol1 = differentiate_solution(ode, case["param_vals"], case["ic"], case["inputs"], prec)[2]
        sol2 = diff_sol_Lie_derivatives(ode, ode.x_vars, case["param_vals"], case["ic"], case["inputs"], prec)
        for i in 1:length(ode.x_vars)
            x = ode.x_vars[i]
            for v in vcat(ode.x_vars, ode.parameters)
                @test sol2[i][v] == [base_ring(ode.poly_ring)(coeff(sol1[(x, v)], i) * factorial(i)) for i in 0:(prec - 1)]
            end
        end
    end

end

max_stars_repo_name: pogudingleb/ExperimentsBound

max_issues_repo_name: pogudingleb/ExperimentsBound

max_forks_repo_name: pogudingleb/ExperimentsBound


text:
using SIMD

@inline mylog(n) = if n == 1 0 else 1 + mylog(n>>1) end

@inline concat(a::Vec{N,T}, b::Vec{M,T}) where {N,M,T} =  Vec((NTuple{N,T}(a)..., NTuple{M,T}(b)...))


text:
using Polynomials: fit
using Statistics: quantile!

function removeMSC(aout,ain,NpY::Integer)
    #Start loop through all other variables
    tmsc, tnmsc = zeros(Union{Float64,Missing},NpY),zeros(Int,NpY)
    fillmsc(1,tmsc,tnmsc,ain,NpY)
    subtractMSC(tmsc,ain,aout,NpY)
    nothing
end

"""
    removeMSC(c)

Removes the mean annual cycle from each time series of a data cube.

**Input Axis** `Time`axis

**Output Axis** `Time`axis
"""
function removeMSC(c;kwargs...)
    NpY = getNpY(c)
    mapCube(
        removeMSC,
        c,
        NpY;
        indims  = InDims("Time"),
        outdims = OutDims("Time"),
        kwargs...
    )
end

"""
    gapFillMSC(c;complete_msc=false)

Fills missing values of each time series in a cube with the mean annual cycle.
If `complete_msc` is set to `true`, the MSC will be filled with a polynomial
in case it still contains missing values due to systematic gaps.

**Input Axis** `Time`axis

**Output Axis** `Time`axis
"""
function gapFillMSC(c;kwargs...)
  NpY=getNpY(c)
  mapCube(gapFillMSC,c,NpY;indims=InDims("Time"),outdims=OutDims("Time"),kwargs...)
end

function gapFillMSC(aout::AbstractVector,ain::AbstractVector,NpY::Integer;complete_msc=false)
  tmsc, tnmsc = zeros(Union{Float64,Missing},NpY),zeros(Int,NpY)
  fillmsc(1,tmsc,tnmsc,ain,NpY)
  if complete_msc
    fill_msc_poly!(tmsc)
  end
  replaceMisswithMSC(tmsc,ain,aout,NpY)
end

function fill_msc_poly!(tmsc)
  mscrep = [tmsc;tmsc;tmsc]
  n = length(tmsc)
  a = gapfillpoly!(mscrep, max_gap = n÷2, nbefore_after = max(3,n÷30))
  tmsc .= view(a,(n+1):(n+n))
end

gapfillpoly(x;max_gap=30,nbefore_after=10, polyorder = 2) =
  mapslices(gapfillpoly!,x,dims="Time",max_gap=max_gap, nbefore_after=nbefore_after, polyorder=polyorder)

"""
    fillgaps_poly(x;max_gap=30,nbefore_after=10, polyorder = 2)

Function for polnomial gap filling. Whenever a gap smaller than `max_gap` is found
the algorithm uses `nbefore_after` time steps before and after the gap to fit
a polynomial of order `polyorder`. The missing alues are then replaced by the
fitted polynomial.
"""
function gapfillpoly!(x;max_gap=30,nbefore_after=10, polyorder = 2)
    x = replace(i->(!ismissing(i) && isfinite(i)) ? i : missing,x)
    a = copy(x)
    workx = Float64[]
    worky = Float64[]
    # find first nonmissing value
    idxstart = findfirst(!ismissing,a)
    idxstart === nothing && return a
    while true
        #find next missing value
        gapstart = findnext(ismissing,a,idxstart)
        gapstart === nothing && break
        gapstop = findnext(!ismissing,a,gapstart)
        gapstop === nothing && break
        if gapstop-gapstart < max_gap
            idxfirst = max(1,gapstart - nbefore_after)
            idxlast  = min(length(a), gapstop + nbefore_after - 1)
            idxr = idxfirst:idxlast
            for (ii,idx) in enumerate(idxr)
                if !ismissing(x[idx])
                    push!(workx,idxr[ii])
                    push!(worky,x[idx])
                end
            end
            if length(workx)>polyorder
                p = fit(workx,worky,polyorder)
                for idx in gapstart:(gapstop-1)
                    a[idx] = p(idx)
                end
            end
        end
        idxstart = gapstop
        empty!(workx)
        empty!(worky)
    end
    a
end

"""
    getMSC(c)

Returns the mean annual cycle from each time series.

**Input Axis** `Time`axis

**Output Axis** `MSC`axis

"""
function getMSC(c;kwargs...)
  N = getNpY(c)
  outdims = OutDims(RangeAxis("MSC",DateTime(1900):Day(ceil(Int,366/N)):DateTime(1900,12,31,23,59,59)),
  outtype = mscouttype(eltype(c)))
  indims = InDims("Time")
  mapCube(getMSC,c,getNpY(c);indims=indims,outdims=outdims,kwargs...)
end

mscouttype(T) = Base.nonmissingtype(T)
mscouttype(::Type{<:Union{Missing,Integer}}) = Float64

function getMSC(aout::AbstractVector,ain::AbstractVector,NpY;imscstart::Int=1)
    nmsc = zeros(Int,NpY)
    fillmsc(imscstart,aout,nmsc,ain,NpY)
end

"Subtracts given msc from input vector"
function subtractMSC(msc::AbstractVector,xin::AbstractVector,xout,NpY)
    imsc=1
    ltime=length(xin)
    for i in 1:ltime
        xout[i] = xin[i]-msc[imsc]
        imsc =imsc==NpY ? 1 : imsc+1 # Increase msc time step counter
    end
end

"Replaces missing values with mean seasonal cycle"
function replaceMisswithMSC(msc::AbstractVector,xin::AbstractArray,xout::AbstractArray,NpY::Integer)
  imsc=1
  for i in eachindex(xin)
    if ismissing(xin[i]) && !ismissing(msc[imsc])
      xout[i]=msc[imsc]
    else
      xout[i]=xin[i]
    end
    imsc= imsc==NpY ? 1 : imsc+1 # Increase msc time step counter
  end
end

"""
    getMedMSC(c)

Returns the median annual cycle from each time series.

**Input Axis** `Time`axis

**Output Axis** `MSC`axis
"""
function getMedSC(c;kwargs...)
  N = getNpY(c)
  outdims = OutDims(RangeAxis("MSC",DateTime(1900):Day(ceil(Int,366/N)):DateTime(1900,12,31,23,59,59)),
  outtype = mscouttype(eltype(c)))
  indims = InDims("Time")
  mapCube(getMedSC,c;indims=indims,outdims=outdims,kwargs...)
end

function getMedSC(aout::AbstractVector{Union{T,Missing}},ain::AbstractVector) where T
    #Reshape the cube to squeeze unimportant variables
    NpY=length(aout)
    yvec=T[]
    q=[convert(T,0.5)]
    for doy=1:NpY
        empty!(yvec)
        for i=doy:NpY:length(ain)
            ismissing(ain[i]) || push!(yvec,ain[i])
        end
        aout[doy] = isempty(yvec) ? missing : quantile!(yvec,q)[1]
    end
    aout
end

"Calculates the mean seasonal cycle of a vector"
function fillmsc(imscstart::Integer,msc::AbstractVector,nmsc::AbstractVector,xin::AbstractVector,NpY)
  imsc=imscstart
  fill!(msc, 0)
  fill!(nmsc,0)
  for v in xin
    if !ismissing(v)
      msc[imsc]  += v
      nmsc[imsc] += 1
    end
    imsc=imsc==NpY ? 1 : imsc+1 # Increase msc time step counter
  end
  for i in 1:NpY msc[i] = nmsc[i] > 0 ? msc[i]/nmsc[i] : missing end # Get MSC by dividing by number of points
end


text:
using SparseJacobian
using Test

@testset "SparseJacobian.jl" begin
    # Write your tests here.
end