text:
module CartPole

using MDPs
using DelimitedFiles
using PGFPlots

export state_space, action_space, transition, reward, simulate,
       viz_policy, viz_trajectory
export THETAMIN, THETAMAX, THETADOTMIN, THETADOTMAX, XMIN, XMAX, XDOTMIN, XDOTMAX, T

# mdp
const MTHETA = 10
const MTHETADOT = 10
const MX = 10
const MXDOT = 10
const N = 1000

# state
const THETAMIN = float(-pi / 2)
const THETAMAX = float(pi / 2)
const THETADOTMIN = -3.0
const THETADOTMAX = 3.0
const XMIN = -2.4
const XMAX = 2.4
const XDOTMIN = -3.5
const XDOTMAX = 3.5

# action
const AMIN = -10
const AMAX = 10

# transition
const GRAVITY = 9.8
const MASSCART = 1.0
const MASSPOLE = 0.1
const LENGTH = 0.5
# const FORCE_MAG = 10
const DT = 0.1  # time intervals

# reward
const K = 15

# test
const T = 200

# visualization
const EPS = 1e-4

function state_space(mtheta::Int64=MTHETA, mthetadot::Int64=MTHETADOT, mx::Int64=MX, mxdot::Int64=MXDOT)
    theta = range(THETAMIN, stop=THETAMAX, length=mtheta)
    thetadot = range(THETADOTMIN, stop=THETADOTMAX, length=mthetadot)
    x = range(XMIN, stop=XMAX, length=mx)
    xdot = range(XDOTMIN, stop=XDOTMAX, length=mxdot)
    return RectangleGrid(theta, thetadot, x, xdot)
end

function action_space(n::Int64=N)
    return range(AMIN, stop=AMAX, length=n)
end

function wrap_around(angle::Float64)
    while angle < -pi / 2
        angle = -pi / 2
    end
    while angle > pi / 2
        angle = pi / 2
    end
    return angle
end

function clip(val::Float64, minval::Float64, maxval::Float64)
    return min(maxval, max(minval, val))
end

function transition(state::Vector{Float64}, action::Float64)
    temp = (action + MASSPOLE * LENGTH * state[2] * state[2] * sin(state[1])) / (MASSPOLE + MASSCART)
    theta_acc = (GRAVITY * sin(state[1]) - cos(state[1]) * temp) / (LENGTH * (4/3 - MASSPOLE * cos(state[1]) * cos(state[1]) / (MASSPOLE + MASSCART)))
    x_acc = temp - MASSPOLE * LENGTH * theta_acc * cos(state[1]) / (MASSPOLE + MASSCART)

    snext = zeros(length(state))
    snext[1] = wrap_around(state[1] + DT * state[2])
    snext[2] = clip(state[2] + DT * theta_acc, THETADOTMIN, THETADOTMAX)
    snext[3] = clip(state[3] + DT * state[4], XMIN, XMAX)
    snext[4] = clip(state[4] + DT * x_acc, XDOTMIN, XDOTMAX)
    return snext
end

function reward(state::Vector{Float64}, action::Float64)
    return (cos(K * state[1]))^4
end

function simulate(mdp::MDP, policy::Policy, state::Vector{Float64})
    trajectory = zeros(T, dimensions(mdp.S))
    actions = zeros(T)
    for t = 1:T
        trajectory[t, :] = state
        if state[1] == 0.0
            trajectory = trajectory[1:t, :]
            actions = actions[1:t - 1]
            break
        end
        action = policy!(policy, get_belief(mdp, state))
        state = mdp.transition(state, action)
        actions[t] = action
    end
    return trajectory, actions
end

function viz_policy(mdp::MDP, policy::Policy, title::String="Policy heatmap", saveplot::Bool=false, name::String="")
    function getmap(p::Float64, v::Float64)
        return policy!(policy, get_belief(mdp, [p, v]))
    end

    p = Axis([Plots.Image(getmap, (THETAMIN + EPS, THETAMAX - EPS),
                      (THETADOTMIN + EPS, THETADOTMAX - EPS),
                      xbins = 250, ybins = 250,
                      colormap = ColorMaps.Named("Blues"))],
          width="14cm", height="14cm",
          xlabel="angle", ylabel="angular speed",
          title=title)
    p
    if saveplot
      PGFPlots.save(name, p)
    end
end

function rad2deg_vec(trajectory::Vector{Float64})
    l = length(trajectory)
    out = zeros(l)
    for i = 1:l
        out[i] = rad2deg(trajectory[i])
    end
    return out
end

function viz_trajectory(trajectory::Matrix{Float64}, actions::Vector{Float64},
                        title1::String="angular position over time", title2::String="control input over time",
                        saveplot::Bool=false, name::String="")
    g = GroupPlot(2, 1, groupStyle="horizontal sep=1.5cm")
    push!(g, Axis([Plots.Linear(rad2deg_vec(trajectory[:, 1]))],
                  width="12cm", height="12cm",
                  xlabel="time", ylabel="angle (deg)",
                  title=title1))
    push!(g, Axis([Plots.Linear(actions)],
                  width="12cm", height="12cm",
                  xlabel="time", ylabel="input",
                  title=title2))
    g
    if saveplot
      PGFPlots.save(name, g)
    end
end

end

max_stars_repo_name: tsamtsam/SV-RL

max_issues_repo_name: tsamtsam/SV-RL

max_forks_repo_name: tsamtsam/SV-RL


text:
#=
Entry point of Modia

* Developer: Hilding Elmqvist, Mogram AB
* First version: December 2020
* License: MIT (expat)

=#

function assert(condition, message)
    if ! condition
        println("ASSERTION: ", message)
    end
end

function printArray(array, heading; order=1:length(array), log=false, number=true, indent=false, spacing=" ")
    if length(array) > 0 && log
        println(heading)
        n = 0
        for v in array
            n += 1
            if number
                print(lpad(order[n], 4), ": ")
            end
            if indent
                print(("|"*spacing)^(n-1))
            end
            println(replace(replace(string(v), "\n" => " "), "  " => " "))
        end
    end
end

function performConsistencyCheck(G, Avar, vActive, parameters, unknowns, states, equations, log=false)
    nUnknowns = length(unknowns) - length(states)
    nEquations = length(equations)

    if nEquations > 0 && nUnknowns != nEquations
        printstyled("The number of unknowns ($nUnknowns) and equations ($nEquations) is not equal!\n", bold=true, color=:red)
        printArray(parameters, "Parameters:", log=true)
        printArray(states, "Potential states:", log=true)
        printArray(setdiff(unknowns,states), "Unknowns:", log=true)
        printArray(equations, "Equations:", log=true)
    end

    ok = nUnknowns == nEquations

    # Check consistency
    EG::Array{Array{Int,1},1} = [G; buildExtendedSystem(Avar)]
    assignEG = matching(EG, length(Avar))
#        @assert all(assignEG .> 0) "  The DAE is structurally singular\n  $assignEG."
    assigned = [a > 0 for a in assignEG]
    if ! (all(assigned)) || ! ok
        if nUnknowns > nEquations
            missingEquations = nUnknowns - nEquations
            printstyled("\nThe DAE has no unique solution, because $missingEquations equation(s) missing!\n", bold=true, color=:red)
        elseif nEquations > nUnknowns
            tooManyEquations = nEquations - nUnknowns
            printstyled("\nThe DAE has no unique solution, because $tooManyEquations equation(s) too many\n", bold=true, color=:red)
        else
            printstyled("\nThe DAE has no unique solution, because it is structurally inconsistent\n" *
                        "  (the number of equations is identical to the number of unknowns)\n", bold=true, color=:red)
        end

        singular = true
        hequ = []
        # Create dummy equations for printing: integrate(x, der(x)) = 0
        for i in 1:length(Avar)
            a = Avar[i]
            if a > 0
                n = unknowns[i]
                der = unknowns[a]
                push!(hequ, "integrate($n, $der) = 0")
            end
        end
        equationsEG = vcat(equations, hequ)
        println("\nAssigned equations after consistency check:")
        printAssignedEquations(equationsEG, unknowns, 1:length(EG), assignEG, Avar, fill(0, length(EG)))
        unassigned = []
        for i in eachindex(assignEG)
            if assignEG[i] == 0
                push!(unassigned, unknowns[i])
            end
        end
        printArray(unassigned, "Not assigned variables after consistency check:", log=true)
        if nUnknowns > nEquations
            missingEquations = nUnknowns - nEquations
            error("Aborting, because $missingEquations equation(s) missing!")
        elseif nEquations > nUnknowns
            tooManyEquations = nEquations - nUnknowns
            error("Aborting, because $tooManyEquations equation(s) too many!")
        else
            error("Aborting, because DAE is structurally inconsistent")
        end
    elseif log
        println("  The DAE is structurally nonsingular.")
    end
end

function assignAndBLT(model, modelName, modelModule, equations, unknowns, modelStructure, Avar, G, states, logDetails, log=false, logTiming=false)
    parameters = modelStructure.parameters
    vActive = [a == 0 for a in Avar]

    if drawIncidence
        showIncidence(equations, unknowns, G, vActive, [], [], [], [[e] for e in 1:length(G)])
    end

    if false # logAssignmentBeforePantelides
        @show G length(Avar) vActive
        assign = matching(G, length(Avar), vActive)
        @show assign

        println("\nAssigned equations:")
        Bequ = fill(0, length(G))
        printAssignedEquations(equations, unknowns, 1:length(G), assign, Avar, Bequ)

        printUnassigned(equations, unknowns, assign, Avar, Bequ, vActive)
    end

    if logTiming
        println("\nCheck consistency of equations by matching extended equation set")
        @timeit to "performConsistencyCheck" performConsistencyCheck(G, Avar, vActive, parameters, unknowns, states, equations, log)
    else
        performConsistencyCheck(G, Avar, vActive, parameters, unknowns, states, equations, log)
    end

    if logTiming
        println("\nIndex reduction (Pantelides)")
        @timeit to "pantelides!" assign, Avar, Bequ = pantelides!(G, length(Avar), Avar)
    else
        assign, Avar, Bequ = pantelides!(G, length(Avar), Avar)
    end
    if logDetails
        @show assign Avar Bequ
    end

    moreVariables = length(Avar) > length(unknowns)
    unknowns = vcat(unknowns, fill(:x, length(Avar) - length(unknowns)))
    for i in 1:length(Avar)
        if Avar[i] > 0
            unknowns[Avar[i]] = derivative(unknowns[i], keys(parameters))
        end
    end
    vActive = [a == 0 for a in Avar]

    moreEquations = length(Bequ) > length(equations)
    equations = vcat(equations, fill(:(a=b), length(Bequ) - length(equations)))
    for i in 1:length(Bequ)
        if Bequ[i] > 0
            if log
                println("Differentiating: ", equations[i])
            end
            equations[Bequ[i]] = simplify(derivative(equations[i], keys(parameters)))
            if log
                println("  Differentiated equation: ", equations[Bequ[i]])
            end
        end
    end

    if moreVariables && log
        printArray(unknowns, "Unknowns after index reduction:", log=logDetails)
    end
    if moreEquations && log
        printArray(equations, "Equations after index reduction:", log=logDetails)
    end

    HG = Array{Array{Int64,1},1}()
    originalEquationIndex = []
    highestDerivativeEquations = []
    for i in 1:length(G)
        if Bequ[i] == 0
            push!(HG, G[i])
            push!(highestDerivativeEquations, equations[i])
            push!(originalEquationIndex, i)
        end
    end

    if logTiming
        println("Assign")
        @timeit to "matching" assign = matching(HG, length(Avar), vActive)
    else
        assign = matching(HG, length(Avar), vActive)
    end
    if logTiming
        println("BLT")
        @timeit to "BLT" bltComponents = BLT(HG, assign)
    else
       bltComponents = BLT(HG, assign)
    end
    if logDetails
        @show HG
        @show bltComponents
    end

    if log
        println("\nSorted highest derivative equations:")
        printSortedEquations(highestDerivativeEquations, unknowns, bltComponents, assign, Avar, Bequ)
    end

    if experimentalTranslation
        return specialTranslation(model, modelName, modelModule, highestDerivativeEquations, unknowns, modelStructure, states, bltComponents, assign, Avar, Bequ)
    end

    if drawIncidence
        assignedVar, unAssigned = invertAssign(assign)
        showIncidence(equations, unknowns, HG, vActive, [], assign, assignedVar, bltComponents, sortVariables=true)
    end

    blt = Vector{Int}[]  # TODO: Change type
    for c in bltComponents
        push!(blt, originalEquationIndex[c])
    end
    for i in eachindex(assign)
        if assign[i] != 0
            assign[i] = originalEquationIndex[assign[i]]
        end
    end
    return (unknowns, equations, G, Avar, Bequ, assign, blt, parameters)
end

function setAvar(unknowns)
    variablesIndices = OrderedDict(key => k for (k, key) in enumerate(unknowns))
    Avar = fill(0, length(unknowns))
    states = []
    derivatives = []
    for i in 1:length(unknowns)
        v = unknowns[i]
        if isexpr(v, :call)
            push!(derivatives, v)
            x = v.args[2]
            ix = variablesIndices[x]
            Avar[ix] = i
            push!(states, unknowns[ix])
        end
    end
    return Avar, states, derivatives
end

mutable struct ModelStructure
    parameters::OrderedDict{Any,Any}
    mappedParameters::OrderedDict{Symbol,Any}
    init::OrderedDict{Any,Any}
    start::OrderedDict{Any,Any}
    variables::OrderedDict{Any,Any}
    flows::OrderedDict{Any,Any}
    inputs::OrderedDict{Any,Any}
    outputs::OrderedDict{Any,Any}
#    components
#    extends
    equations::Array{Expr,1}
    hideResults::OrderedSet{Any}   # Do not store these variables in the result data structure
end

ModelStructure() = ModelStructure(OrderedDict(), OrderedDict{Symbol,Any}(), OrderedDict(), OrderedDict(), OrderedDict(), OrderedDict(), OrderedDict(), OrderedDict(), Expr[], OrderedSet())

function printDict(label, d)
    if length(d) > 0
        print(label)
        for (k,v) in d
            print(k)
            if v != nothing
                print(" => ", v)
            end
            print(", ")
        end
        println()
    end
end

function printModelStructure(m, label=""; log=false)
    if log
        println("ModelStructure ", label)
        printDict("  parameters:  ", m.parameters)
        printDict("  mappedParameters:  ", m.mappedParameters)
        printDict("  init:        ", m.init)
        printDict("  start:       ", m.start)
        printDict("  variables:   ", m.variables)
        printDict("  flows:       ", m.flows)
        printDict("  inputs:      ", m.inputs)
        printDict("  outputs:     ", m.outputs)
        println("  equations: ")
        for e in m.equations
            println("    ", e)
        end
    end
end

prependDict(dict, prefix) = OrderedDict([prepend(k, prefix) => prepend(v, prefix) for (k,v) in dict])

function mergeModelStructures(parent::ModelStructure, child::ModelStructure, prefix)
    if length(child.mappedParameters) > 0
        merge!(parent.parameters, child.parameters)
        parent.mappedParameters[prefix] = child.mappedParameters

        merge!(parent.init, child.init)
        parent.mappedParameters[prefix] = child.mappedParameters

        merge!(parent.start, child.start)
        parent.mappedParameters[prefix] = child.mappedParameters
    end

    merge!(parent.variables, child.variables)
    merge!(parent.flows, child.flows)
    union!(parent.hideResults, child.hideResults)

    push!(parent.equations, prepend(child.equations, prefix)...)
end

function replaceLinearIntegerEquations(Gint, eInt, GcInt, unknowns, equations)
    for i = 1:length(Gint)
        e = Gint[i]
        if length(e) > 0
            # Construct equation
            rhs = 0
            for j = 1:length(e)
                v  = e[j]
                vc = GcInt[i][j]
                rhs = add(rhs, mult(vc, unknowns[v]))
            end
            equ = :(0 = $rhs)
            equations[eInt[i]] = equ
        else
            equations[eInt[i]] = :(0 = 0)
        end
    end
end

function performAliasReduction(unknowns, equations, Avar, logDetails, log)

    @timeit to "enumerate(unknowns)" variablesIndices = OrderedDict(key => k for (k, key) in enumerate(unknowns))

    linearVariables = []
    nonlinearVariables = []
    linearEquations = Vector{Int}()
    nonlinearEquations = []
    G = Vector{Int}[] # Array{Array{Int64,1},1}()
    Gint = Array{Array{Int64,1},1}()
    GcInt = Array{Array{Int64,1},1}()

    @timeit to "build graph information" for i in 1:length(equations)
        e = equations[i]
        nameIncidence, coefficients, rest, linear = getCoefficients(e)
        incidence = [variablesIndices[n] for n in nameIncidence if n in keys(variablesIndices)]
        unique!(incidence)
        push!(G, incidence)

        if linear && all([c == 1 || c == -1 for c in coefficients]) && rest == 0 && all(n in keys(variablesIndices) for n in nameIncidence)
            push!(linearVariables, nameIncidence...)
            push!(linearEquations, i)
            push!(Gint, incidence)
            push!(GcInt, coefficients)
        else
            push!(nonlinearVariables, nameIncidence...)
            push!(nonlinearEquations, e)
        end
    end

    if logDetails
        @show G
        @show Avar
    end

    @timeit to "unique!(linearVariables)" unique!(linearVariables)
    @timeit to "unique!(nonlinearVariables)" unique!(nonlinearVariables)
    @timeit to "setdiff" vSolveInLinearEquations = setdiff(linearVariables, nonlinearVariables)
    if logDetails
        printArray(linearVariables, "linearVariables:", number=false)
        printArray(equations[linearEquations], "linearEquations:", number=false)
#        printArray(nonlinearVariables, "nonlinearVariables:", number=false)
#        printArray(nonlinearEquations, "nonlinearEquations:", number=false)
        @show linearEquations vSolveInLinearEquations Gint GcInt
    end

    GCopy = deepcopy(G)
    AvarCopy = deepcopy(Avar)
    @timeit to "simplifyLinearIntegerEquations!" (vEliminated, vProperty, nvArbitrary, redundantEquations) = simplifyLinearIntegerEquations!(GCopy, linearEquations, GcInt, AvarCopy)

    if logDetails
        @show vEliminated vProperty nvArbitrary redundantEquations
    end
    if log
        printSimplifiedLinearIntegerEquations(GCopy, linearEquations, GcInt, vEliminated, vProperty, nvArbitrary, redundantEquations, v->string(unknowns[v]), printTest=false)
    end

    substitutions = OrderedDict()
    @timeit to "build substitutions" if length(vEliminated) - nvArbitrary > 0
        for v in vEliminated[nvArbitrary+1:end]
            if isZero(vProperty, v)
                substitutions[unknowns[v]] = 0.0
            elseif isAlias(vProperty, v)
                substitutions[unknowns[v]] = unknowns[alias(vProperty,v)]
            else
                substitutions[unknowns[v]] = :(- $(unknowns[negAlias(vProperty,v)]))
            end
        end
    end

    @timeit to "replaceLinearIntegerEquations" replaceLinearIntegerEquations(GCopy[linearEquations], linearEquations, GcInt, unknowns, equations)
#    printArray(equations, "new equations:", log=log)

    reducedEquations = []
    @timeit to "substitute" for ei in 1:length(equations)
        e = equations[ei]
        if e != :(0 = 0)
            push!(reducedEquations, substitute(substitutions, e))
        end
    end

    reducedUnknowns = []
    for vi in 1:length(unknowns)
        if isNotEliminated(vProperty, vi)
            push!(reducedUnknowns, unknowns[vi])
        end
    end

    reducedVariablesIndices = OrderedDict(key => k for (k, key) in enumerate(reducedUnknowns))
    reducedG = Vector{Int}[] # Array{Array{Int64,1},1}()
    @timeit to "build reducedG" for i in 1:length(reducedEquations)
        e = reducedEquations[i]
        nameIncidence = Incidence[]
        findIncidence!(e, nameIncidence)
        incidence = [reducedVariablesIndices[n] for n in nameIncidence if n in keys(reducedVariablesIndices)]
        unique!(incidence)
        push!(reducedG, incidence)
    end

    reducedAvar, reducedStates, reducedDerivatives = setAvar(reducedUnknowns)
    return reducedEquations, reducedUnknowns, reducedAvar, reducedG, reducedStates, vEliminated, vProperty
end

function stateSelectionAndCodeGeneration(modStructure, Gexplicit, name, modelModule, buildDict, FloatType, TimeType, init, start, inputs, outputs, vEliminated, vProperty, unknownsWithEliminated, mappedParameters, hideResults;
    unitless=false, logStateSelection=false, logCode=false, logExecution=false, logCalculations=false, logTiming=false, evaluateParameters=false, saveCodeOnFile="")
    (unknowns, equations, G, Avar, Bequ, assign, blt, parameters) = modStructure
    Goriginal = deepcopy(G)
    function getSolvedEquationAST(e, v)
        (solution, solved) = solveEquation(equations[e], unknowns[v])
        unknown = solution.args[1]
        expr = solution.args[2]
###        solution = :(_variables[$(v+1)] = $unknown = $expr)
        equ = string(equations[e])
        @assert solved "Equation not solved for $(unknowns[v]): $(equations[e])"
#        sol = string(solution)
#        solution = prepend(makeDerivativeVar(solution, components), :m)
        solution = substituteForEvents(solution)
        if false
    #        solution = :(begin println("Executing: ", $sol); $solution end)
            if sol != equ
                solution = :(try $solution; catch e; println("\n\nFailure executing: ", $sol); println("Original equation: ", $equ, "\n\n"); rethrow(e) end)
            else
                solution = :(try $solution; catch e; println("\n\nFailure executing: ", $sol, "\n\n"); rethrow(e) end)
            end
    #        solution = :(try $solution; catch e; println("Failure executing: ", $sol); println(sprint(showerror,e)); rethrow(e) end)
    #        solution = :(try $solution; catch e; println("Failure executing: ", $sol); printstyled(stderr,"ERROR: ", bold=true, color=:red);
    #        printstyled(stderr,sprint(showerror,e), color=:light_red); println(stderr); end)
        end
        solution = makeDerVar(solution, parameters, inputs, FloatType, evaluateParameters)
        if logCalculations
            var = string(unknowns[v])
            solutionString = string(solution)
            return :(println("Calculating: ", $solutionString); $solution; println("  Result: ", $var, " = ", upreferred.($(solution.args[1]))))
        else
            return solution
        end
    end

    function getResidualEquationAST(e, residualName)
        eq = equations[e] # prepend(makeDerivativeVar(equations[e], components), :m)
        lhs = eq.args[1]
        rhs = eq.args[2]
        if isexpr(rhs, :call) && rhs.args[1] == :_DUPLICATEEQUATION
            return nothing
        end
        if isexpr(lhs, :tuple) && all(a == 0 for a in lhs.args) || lhs == :(0)
            eq_rhs = makeDerVar(:($rhs), parameters, inputs, FloatType, evaluateParameters)
            if unitless
                eqs = eq_rhs
            else
                eqs = :(Modia.Unitful.ustrip.($eq_rhs))
            end
        #elseif isexpr(lhs, :tuple) && isexpr(rhs, :call) && unitless
        #    eq_rhs = makeDerVar(:($rhs), parameters, inputs, evaluateParameters)
        #    eq_lhs = makeDerVar(:($lhs), parameters, inputs, evaluateParameters)
        #    eqs =  :( ($eq_rhs .-= $eq_lhs) )
        else
            eq_rhs = makeDerVar(:($rhs), parameters, inputs, FloatType, evaluateParameters)
            eq_lhs = makeDerVar(:($lhs), parameters, inputs, FloatType, evaluateParameters)
            if unitless
                eqs = :( $eq_rhs .- $eq_lhs )
            else
                eqs = :( Modia.Unitful.ustrip.($eq_rhs) .- Modia.Unitful.ustrip.($eq_lhs))
            end
        end
        residual = :(Modia.appendVariable!(_leq_mode.residuals, $eqs))
        residString = string(eqs)
        if logCalculations
             return :(println("Calculating residual: ", $residString); $residualName = $eqs; println("  Residual: ", $residualName) )
#            return makeDerVar(:(dump($(makeDerVar(eq.args[2]))); dump($(makeDerVar(eq.args[1]))); $residual; println($residualName, " = ", upreferred.(($(eq.args[2]) - $(eq.args[1]))))))
        else
            return residual
        end
    end

    function isLinearEquation(e_original, v_original)
        equ = equations[e_original]
        var = unknowns[v_original]
        return isLinear(equ, var)
    end

    function isSolvableEquation(e_original, v_original)
        equ = equations[e_original]
        var = unknowns[v_original]
        #(solution, solved) = solveEquation(equ, var)
        #return solved
        (rest, factor, linear) = ModiaBase.linearFactor(equ, var)
        return linear && (typeof(factor) <: Number && factor != 0)
    end

    hasParticles(value) = typeof(value) <: MonteCarloMeasurements.StaticParticles ||
                          typeof(value) <: MonteCarloMeasurements.Particles

    invAvar = ModiaBase.revertAssociation(Avar)

    function var_unit(v)
        int_v = invAvar[v]
        if int_v > 0
            # v is a derivative variable
            var = unknowns[int_v]
        else
            # v is a non-differentiated variable
            var = unknowns[v]
        end

        if var in keys(init)
            value = eval(init[var])
        elseif var in keys(start)
            value = eval(start[var])
        else
            value = 0.0
        end
        if int_v > 0
            value = value / u"s"
        end
        # if length(value) == 1
        if ! (typeof(value) <: AbstractArray)
            un = unit(value)
        else
            un = unit.(value)
            @assert all([un[i] == un[1] for i in 2:length(un)]) "The unit of all elements of state vector must be equal: $var::$(value)"
            un = un[1]
        end
        return unitAsString(un)
    end

    function var_startInitFixed(v_original)
        v = unknowns[v_original]
        value::Any  = nothing
        fixed::Bool = false
        if v in keys(init)
            value = eval(init[v])
            fixed = true
        elseif v in keys(start)
            value = eval(start[v])
        end
        return (value, fixed)
    end

    function eq_equation(e)
        return replace(replace(string(equations[e]), "\n" => " "), "  " => " ")
    end

    # ----------------------------------------------------------------------------
    solvedAST = []
    Gsolvable = copy(G)
    juliaVariables  = [Symbol(u) for u in unknowns]
    stringVariables = [String(Symbol(v)) for v in unknowns]
    allEquations    = equations
    #    @show stringVariables equations G blt assign Avar Bequ Gsolvable

    stateSelectionFunctions = StateSelectionFunctions(
        var_name               = v -> stringVariables[v],
        var_julia_name         = v -> juliaVariables[v],
        var_unit               = var_unit,
        var_startInitFixed     = var_startInitFixed,
        var_nominal            = v_original -> NaN,
        var_unbounded          = v_original -> false,
        equation               = eq_equation,
        isSolvableEquation     = isSolvableEquation,
        isLinearEquation       = isLinearEquation,
        getSolvedEquationAST   = getSolvedEquationAST,
        getResidualEquationAST = getResidualEquationAST,
        showMessage            = (message;severity=0,from="???",details="",variables=Int[],equations=Int[]) ->
                                    Modia.showMessage(message, severity, from, details,
                                                      stringVariables[variables],
                                                      string.(allEquations[equations]))
        )

    if length(equations) == 0
        return nothing
    end

    if logTiming
        println("Get sorted and solved AST")
        @timeit to "getSortedAndSolvedAST" (success,AST,equationInfo) = getSortedAndSolvedAST(
            Goriginal, Gexplicit, Array{Array{Int64,1},1}(blt), assign, Avar, Bequ, stateSelectionFunctions;
            unitless, log = logStateSelection, modelName = name)
    else
        (success,AST,equationInfo) = getSortedAndSolvedAST(
            Goriginal, Gexplicit, Array{Array{Int64,1},1}(blt), assign, Avar, Bequ, stateSelectionFunctions;
            unitless, log = logStateSelection, modelName = name)
    end

    if ! success
        error("Aborting")
    end

    if useNewCodeGeneration
        AST, newFunctions = prepareGenerateGetDerivatives(AST, modelModule, functionSize, juliaVariables)
    else
        newFunctions = Vector{Expr}()
    end

#    @show AST
#    @show equationInfo

#=
    selectedStates = [v.x_name_julia for v in equationInfo.x_info]

    startValues = []
    for v in selectedStates
        v = Meta.parse(string(v))
        if v in keys(init)
            value = eval(init[v])
#            value = [0.0 0.0]  # TODO: Fix
#            @show value
            if hasParticles(value)
                push!(startValues, value)
            else
                push!(startValues, value...)
            end
        elseif v in keys(start)
            value = eval(start[v])
            if hasParticles(value)
                push!(startValues, value)
            else
                push!(startValues, value...)
            end
        else
            push!(startValues, 0.0)
        end
    end
    if logCode
        println("startValues = ", startValues)
    end
=#

    vSolvedWithInit = equationInfo.vSolvedWithFixedTrue
    vSolvedWithInitValuesAndUnit = OrderedDict{String,Any}()
    for name in vSolvedWithInit
        nameAsExpr = Meta.parse(name)
        if haskey(init, nameAsExpr)
            vSolvedWithInitValuesAndUnit[name] = eval( init[nameAsExpr] )
        else
            @warn "Internal issue of Modia: $name is assumed to have an init-value, but it is not found."
        end
    end

    # Generate code
    nCrossingFunctions, nAfter, nClocks, nSamples, previousVars, preVars, holdVars = getEventCounters()
    previousVars = Symbol.(previousVars)
    preVars = Symbol.(preVars)
    holdVars = Symbol.(holdVars)

    # Variables added to result
    extraResults = vcat(:time, setdiff([Symbol(u) for u in unknowns],
                                       [Symbol(h) for h in hideResults],
                                        Symbol[Symbol(xi_info.x_name_julia)     for xi_info in equationInfo.x_info],
                                        Symbol[Symbol(xi_info.der_x_name_julia) for xi_info in equationInfo.x_info]))

    if true # logTiming
#        println("Generate code")
        if useNewCodeGeneration
            @timeit to "generate_getDerivativesNew!" code = generate_getDerivativesNew!(AST, newFunctions, modelModule, equationInfo, [:(_p)], vcat(:time, [Symbol(u) for u in unknowns]), previousVars, preVars, holdVars, :getDerivatives, hasUnits = !unitless)
        else
            @timeit to "generate_getDerivatives!" code = generate_getDerivatives!(FloatType, TimeType, AST, equationInfo, [:(_p)], extraResults, previousVars, preVars, holdVars, :getDerivatives, hasUnits = !unitless)
        end
    else
#        code = generate_getDerivatives!(AST, equationInfo, Symbol.(keys(parameters)), extraResults, :getDerivatives, hasUnits = !unitless)
        code = generate_getDerivativesNew!(AST, newFunctions, modelModule, equationInfo, [:(_p)], extraResults, previousVars, preVars, holdVars, :getDerivatives, hasUnits = !unitless)
    end
    if logCode
        #@show mappedParameters
        showCodeWithoutComments(code)
    end

    # Compile code

#    generatedFunction = @RuntimeGeneratedFunction(modelModule, code)
    #getDerivatives = Core.eval(modelModule, code)

    if saveCodeOnFile != ""
        println("  Save generated code on file ", joinpath(pwd(), saveCodeOnFile))
        open(saveCodeOnFile, "w") do io
            print(io, replace(sprint(show,code), r"( *#= .*=#\n)|(#= .*=#)" => ""))
        end
    end

    if logTiming
        println("eval code")
        @time @timeit to "eval(code)" Core.eval(modelModule, code)
    else
        Core.eval(modelModule, code)
    end
    getDerivatives = modelModule.getDerivatives

    # If generatedFunction is not packed inside a function, DifferentialEquations.jl crashes
#    getDerivatives(derx,x,m,time) = generatedFunction(derx, x, m, time)

#    convertedStartValues = convert(Vector{FloatType}, [ustrip(v) for v in startValues])  # ustrip.(value) does not work for MonteCarloMeasurements
#    @show mappedParameters
     x_startValues = initialStateVector(equationInfo, FloatType)

#    println("Build SimulationModel")

    model = @timeit to "build SimulationModel" SimulationModel{FloatType,TimeType}(modelModule, name, buildDict, getDerivatives, equationInfo, x_startValues, previousVars, preVars, holdVars,
                                         mappedParameters, extraResults;
                                         vSolvedWithInitValuesAndUnit, vEliminated, vProperty,
                                         var_name = (v)->string(unknownsWithEliminated[v]),
                                         nz=nCrossingFunctions, nAfter=nAfter,  unitless=unitless)

    # Execute code
    if logExecution
        println("\nExecute getDerivatives")
    #    @show startValues
        #derx = deepcopy(x_startValues)     # deepcopy(convertedStartValues) # To get the same type as for x (deepcopy is needed for MonteCarloMeasurements)
        println("First executions of getDerivatives")
        @timeit to "execute getDerivatives" try
            #@time Base.invokelatest(getDerivatives, derx, x_startValues, model, convert(TimeType, 0.0))
            #@time Base.invokelatest(getDerivatives, derx, x_startValues, model, convert(TimeType, 0.0))
            @time init!(model)  # getDerivatives is called
            #@time invokelatest_getDerivatives_without_der_x!(x_startValues, model, convert(TimeType, 0.0))
            #@time invokelatest_getDerivatives_without_der_x!(x_startValues, model, convert(TimeType, 0.0))
#           @show derx
        catch e
            error("Failed: ", e)
            return nothing
        end
    end
    return model
end

function duplicateMultiReturningEquations!(equations)
    duplicatedEquations = []
    for e in equations
        if isexpr(e, :(=)) && isexpr(e.args[1], :tuple) || typeof(e.args[1]) <: NTuple
            lhs = e.args[1]
            if typeof(lhs) <: NTuple
                nargs = length(lhs)
            else
                nargs = length(lhs.args)
            end

            if e.args[2].args[1] != :implicitDependency
                func = :_DUPLICATEEQUATION
            else
                func = :_DUPLICATIMPLICITDEPENDENCY
            end

            nameIncidence = Incidence[]
            findIncidence!(e, nameIncidence, false)
            bandWidth = length(nameIncidence)-nargs+2
            #@show length(nameIncidence) bandWidth
            for i in 1:(nargs-1)
                if false
                    rhs = Expr(:call, :_DUPLICATEEQUATION, nameIncidence...)
                else
                    indices =[]
                    for j in 1:bandWidth
                        push!(indices, mod(j+i-1,length(nameIncidence))+1)
                    end
                    rhs = Expr(:call, :_DUPLICATEEQUATION, nameIncidence[indices]...)
                end
                newE = :(0 = $rhs)
                push!(duplicatedEquations, newE)
            end
        end
    end
    append!(equations, duplicatedEquations)
end

appendSymbol(path::Nothing, name::Symbol) = name
appendSymbol(path         , name::Symbol) = :( $path.$name )

"""
    modifiedModel = buildSubModels!(model, modelModule, FloatType, TimeType, buildDict::OrderedDict)

Traverse `model` and for every `<subModel>` that is a `Model(..)` and has a key-value pair
`:_buildFunction = <buildFunction>` and optionally `:_buildOption=<buildOption>`, call

```
buildCode = <buildFunction>(<subModel>, modelModule, FloatType::Type, TimeType::Type,
                            buildDict::OrderedDict{String,Any},
                            modelPath::Union{Expr,Symbol,Nothing},
                            buildOption = <buildOption>)
```

The`buildCode` is merged to the corresponding `<subModel>` in the calling environment.
The arguments of `<buildFunction>`are:

- `subModel`: The returned `buildCode` is merged to `submodel`
- `FloatType`, `TimeType`: Types used when instantiating `SimulationModel{FloatType,TimeType}`
- `modelPath`: Path upto `<subModel>`, such as: `:( a.b.c )`.
- `buildDict`: Dictionary, that will be stored in the corresponding SimulationModel instance and
               that allows to store information about the build-process,
               typically with key `string(modelPath)` (if modelPath==Nothing, key="" is used).
- `buildOption`: Option used for the generation of `buildCode`.

Note, keys `_buildFunction` and `_buildOption` are deleted from the corresponding `<subModel>`.
"""
function buildSubModels!(model::AbstractDict, modelModule, FloatType::Type, TimeType::Type,
                         buildDict::OrderedDict{String,Any}; path::Union{Expr,Symbol,Nothing}=nothing)
    if haskey(model, :_buildFunction)
        buildFunction = model[:_buildFunction]
        delete!(model, :_buildFunction)
        quotedPath = Meta.quot(path)
        if haskey(model, :_buildOption)
            buildOption = model[:_buildOption]
            delete!(model, :_buildOption)
            buildCode = Core.eval(modelModule, :($buildFunction($model, $FloatType, $TimeType, $buildDict, $quotedPath, buildOption=$buildOption)) )
        else
            buildCode = Core.eval(modelModule, :($buildFunction($model, $FloatType, $TimeType, $buildDict, $quotedPath)))
        end
        return  model | buildCode
    end

    for (key,value) in model
        if typeof(value) <: OrderedDict && haskey(value, :_class) && value[:_class] == :Model
            model[key] = buildSubModels!(value, modelModule, FloatType, TimeType, buildDict; path=appendSymbol(path,key))
        end
    end
    return model
end

"""
    modelInstance = @instantiateModel(model; FloatType = Float64, aliasReduction=true, unitless=false,
        evaluateParameters=false, saveCodeOnFile="", log=false, logModel=false, logDetails=false, logStateSelection=false,
        logCode=false,logExecution=logExecution, logCalculations=logCalculations, logTiming=false)

Instantiates a model, i.e. performs structural and symbolic transformations and generates a function for calculation of derivatives suitable for simulation.

* `model`: model (declarations and equations)
* `FloatType`: Variable type for floating point numbers, for example: Float64, Measurements{Float64}, StaticParticles{Float64,100}, Particles{Float64,2000}
* `aliasReduction`: Perform alias elimination and remove singularities
* `unitless`: Remove units (useful while debugging models and needed for MonteCarloMeasurements)
* `evaluateParameters`: Use evaluated parameters in the generated code.
* `saveCodeOnFile`: If non-empty string, save generated code in file with name `saveCodeOnFile`.
* `log`: Log the different phases of translation
* `logModel`: Log the variables and equations of the model
* `logDetails`: Log internal data during the different phases of translation
* `logStateSelection`: Log details during state selection
* `logCode`: Log the generated code
* `logExecution`: Log the execution of the generated code (useful for timing compilation)
* `logCalculations`: Log the calculations of the generated code (useful for finding unit bugs)
* `logTiming`: Log timing of different phases
* `return modelInstance prepared for simulation`
"""
macro instantiateModel(model, kwargs...)
    modelName = string(model)
    source = string(__source__.file)*":"*string(__source__.line)

    code = :( instantiateModel($model; modelName=$modelName, modelModule=@__MODULE__, source=$source, $(kwargs...) ) )
    return esc(code)
end

"""
See documentation of macro [`@instantiateModel`]
"""
function instantiateModel(model; modelName="", modelModule=nothing, source=nothing, FloatType = Float64, aliasReduction=true, unitless=false,
    log=false, logModel=false, logDetails=false, logStateSelection=false, logCode=false,
    logExecution=logExecution, logCalculations=logCalculations, logTiming=false, evaluateParameters=false, saveCodeOnFile="")
    if isMonteCarloMeasurements(FloatType) && !unitless
        unitless=true
        printstyled("  @instantiateModel(...,unitless=true, ..) set automatically, because\n  FloatType=MonteCarloMeasurements often fails if units are involved.\n", color=:red)
    end

    #try
    #    model = JSONModel.cloneModel(model, expressionsAsStrings=false)
        println("\nInstantiating model $modelName\n  in module: $modelModule\n  in file: $source")
        resetEventCounters()
        global to = TimerOutput()

        modelStructure = ModelStructure()

        if isexpr(model, :quote)
            model = eval(model) # model defined with macro
        end

        if typeof(model) <: NamedTuple || typeof(model) <: Dict || typeof(model) <: OrderedDict
            # Traverse model and execute functions _buildFunction(..), to include code into sub-models
            buildDict = OrderedDict{String,Any}()
            TimeType = if FloatType <: Measurements.Measurement ||
                          FloatType <: MonteCarloMeasurements.AbstractParticles;
                          baseType(FloatType) else FloatType end  # baseType(..) is defined in CodeGeneration.jl
            model = buildSubModels!(model, modelModule, FloatType, TimeType, buildDict)

            if logModel
                @showModel(model)
            end
            if false
                println("drawModel")
                jsonDiagram = drawModel(modelName, model)
                println()
#                println("jsonModel")
#                JSON.print(jsonDiagram, 2)
            end

            if logTiming
                println("Flatten")
                @timeit to "flatten" flattenModelTuple!(model, modelStructure, modelName, to; unitless, log)
            else
                flattenModelTuple!(model, modelStructure, modelName, to; unitless, log)
            end

            if experimentalTranslation
                interface = buildInterface(model, modelStructure)
            end

            pars = OrderedDict{Symbol, Any}([(Symbol(k),v) for (k,v) in modelStructure.parameters])
            if length(modelStructure.equations) > 0
                flatModel = Model() | pars | Model(equations = [removeBlock(e) for e in modelStructure.equations])
            else
                flatModel = Model() | pars
            end
            printModelStructure(modelStructure, log=false)
            name = modelName
        else
            @show model
            error("Invalid model format")
        end

        duplicateMultiReturningEquations!(modelStructure.equations)

        allVariables = Incidence[]
        if logTiming
            println("Find incidence")
            @timeit to "findIncidence!" findIncidence!(modelStructure.equations, allVariables)
        else
            findIncidence!(modelStructure.equations, allVariables)
        end
        unique!(allVariables)

        if ! experimentalTranslation
            unknowns = setdiff(allVariables, keys(modelStructure.parameters), keys(modelStructure.inputs), [:time, :instantiatedModel, :_leq_mode, :_x])
        else
            unknowns = setdiff(allVariables, keys(modelStructure.parameters), [:time, :instantiatedModel, :_leq_mode, :_x])
        end
        Avar, states, derivatives = setAvar(unknowns)
        vActive = [a == 0 for a in Avar]

        if logStatistics || log
            println("Number of states:    ", length(states))
            if length(unknowns)-length(states) != length(modelStructure.equations)
                println("Number of unknowns:  ", length(unknowns)-length(states))
            end
            println("Number of equations: ", length(modelStructure.equations))
        end
        printArray(modelStructure.parameters, "Parameters:", log=log)
        printArray(modelStructure.inputs, "Inputs:", log=log)
        printArray(modelStructure.outputs, "Outputs:", log=log)
        printArray(states, "Potential states:", log=log)
        printArray(setdiff(unknowns, states), "Unknowns:", log=log)
        printArray(modelStructure.equations, "Equations:", log=log)
        if logDetails
            @show modelStructure.parameters modelStructure.mappedParameters modelStructure.init modelStructure.start modelStructure.variables modelStructure.flows modelStructure.inputs modelStructure.outputs
        end

        unknownsWithEliminated = unknowns
        if aliasReduction
            if log || logTiming
                println("Perform alias elimination and remove singularities")
            end
#            @show unknowns
#            @show keys(modelStructure.inputs) keys(modelStructure.outputs)
###            nonInputAndOutputs = setdiff(unknowns, keys(modelStructure.inputs), keys(modelStructure.outputs))
            if logTiming
                @timeit to "performAliasReduction" equations, unknowns, Avar, G, states, vEliminated, vProperty = performAliasReduction(unknowns, modelStructure.equations, Avar, logDetails, log)
                println("Number of reduced unknowns:  ", length(unknowns)-length(states))
                println("Number of reduced equations: ", length(equations))
            else
                equations, unknowns, Avar, G, states, vEliminated, vProperty = performAliasReduction(unknowns, modelStructure.equations, Avar, logDetails, log)
            end
            printArray(states, "States:", log=log)
            printArray(setdiff(unknowns, states), "Unknowns after alias reduction:", log=log)
            printArray(equations, "Equations after alias reduction:", log=log)
        else
            equations = modelStructure.equations
            G = Vector{Int}[] # Array{Array{Int64,1},1}()
            variablesIndices = OrderedDict(key => k for (k, key) in enumerate(unknowns))

            println("Build incidence matrix")
            unknownsSet = Set(unknowns)
            @timeit to "build incidence matrix" for i in 1:length(equations)
                e = equations[i]
                nameIncidence, coefficients, rest, linear = getCoefficients(e)

                incidence = [] # [variablesIndices[n] for n in nameIncidence if n in unknownsSet]
                for n in nameIncidence
                    if n in keys(variablesIndices)
                        push!(incidence, variablesIndices[n])
                    end
                end
                unique!(incidence)
                push!(G, incidence)
            end
            vEliminated = Int[]
            vProperty   = Int[]
        end

        modStructure = assignAndBLT(model, modelName, modelModule, equations, unknowns, modelStructure, Avar, G, states, logDetails, log, logTiming)
        (unknowns, equations, G, Avar, Bequ, assign, blt, parameters) = modStructure

        Gexplicit = Vector{Int}[]
        variablesIndices = OrderedDict(key => k for (k, key) in enumerate(unknowns))
        @timeit to "build explicit incidence matrix" for i in 1:length(equations)
            e = equations[i]
            if isexpr(e.args[2], :call) && e.args[2].args[1] in [:implicitDependency, :_DUPLICATIMPLICITDEPENDENCY]
                if e.args[2].args[1] == :_DUPLICATIMPLICITDEPENDENCY
                    e.args[2] = e.args[2].args[2]
                    e.args[2].args[1] = :_DUPLICATEEQUATION
                else
                    e.args[2] = e.args[2].args[2]
                end

                nameIncidence = Incidence[]
                findIncidence!(e, nameIncidence)

                incidence = [] # [variablesIndices[n] for n in nameIncidence if n in unknownsSet]
                for n in nameIncidence
                    if n in keys(variablesIndices)
                        push!(incidence, variablesIndices[n])
                    end
                end
                unique!(incidence)
                push!(Gexplicit, incidence)
                equations[i] = e
            else
                push!(Gexplicit, G[i])
            end
        end

        if false
            println("Explicit equations:")
            for e in equations
                println(e)
            end

            @show G Gexplicit
        end

        if ! experimentalTranslation
            inst = stateSelectionAndCodeGeneration(modStructure, Gexplicit, name, modelModule, buildDict, FloatType, TimeType, modelStructure.init, modelStructure.start, modelStructure.inputs, modelStructure.outputs,
                vEliminated, vProperty, unknownsWithEliminated, modelStructure.mappedParameters, modelStructure.hideResults;
                unitless, logStateSelection, logCode, logExecution, logCalculations, logTiming, evaluateParameters, saveCodeOnFile)
        else
            interface[:equations] = modStructure[:equations]
            return interface
        end

        if logTiming
            show(to, compact=true)
            println()
        end

        inst #, flatModel
#=
    catch e
        if isa(e, ErrorException)
            println()
            printstyled("Model error: ", bold=true, color=:red)
            printstyled(e.msg, "\n", bold=true, color=:red)
            printstyled("Aborting @instantiateModel($modelName,...) in $modelModule.\n", bold=true, color=:red)
            println()
#            Base.rethrow()
        else
            Base.rethrow()
        end
        if logTiming
            show(to, compact=true)
            println()
        end
    end
=#

end

ext: jl


text:
# Cauchy distribution

export Cauchy, HalfCauchy

@parameterized Cauchy(μ,σ) ≃ (1/π) * Lebesgue(ℝ)

@kwstruct Cauchy()
@kwstruct Cauchy(μ)
@kwstruct Cauchy(σ)
@kwstruct Cauchy(μ,σ)
@kwstruct Cauchy(ω)
@kwstruct Cauchy(μ,ω)

for N in AFFINEPARS
    @eval begin
        proxy(d::Cauchy{$N}) = affine(params(d), Cauchy())
        logdensity(d::Cauchy{$N}, x) = logdensity(proxy(d), x)
        basemeasure(d::Cauchy{$N}) = basemeasure(proxy(d))
    end
end

# @affinepars Cauchy

function logdensity(d::Cauchy{()} , x)
    return -log1p(x^2)
end

function density(d::Cauchy{()} , x)
    return inv(1 + x^2)
end

Base.rand(rng::AbstractRNG, T::Type, μ::Cauchy{()}) = randn(rng, T) / randn(rng, T)

≪(::Cauchy, ::Lebesgue{X}) where X <: Real = true

TV.as(::Cauchy) = asℝ

@half Cauchy

HalfCauchy(σ) = HalfCauchy(σ=σ)

distproxy(d::Cauchy{()}) = Dists.Cauchy()
distproxy(d::Cauchy{(:μ,)}) = Dists.Cauchy(d.μ, 1.0)
distproxy(d::Cauchy{(:σ,)}) = Dists.Cauchy(0.0, d.σ)
distproxy(d::Cauchy{(:μ,:σ)}) = Dists.Cauchy(d.μ, d.σ)
distproxy(d::Cauchy{(:ω,)}) = Dists.Cauchy(0.0, inv(d.ω))
distproxy(d::Cauchy{(:μ,:ω)}) = Dists.Cauchy(d.μ, inv(d.ω))


text:
using SubspacePreservingSparsification
using SparseArrays
using LinearAlgebra
M = [17.05 16.98 0.3 6.99 7; 16.98 0.2 7.1 6.9 0; 0.3 7.1 -12 0.01 17; 6.99 6.9 0.01 -11.97 0; 7 0 17 0 -0.1]

M_id = sparsity_pattern(M, 0.6, 2)

bin_sparse_matrix!(M, M_id, 200)

binned_minimization(M_id, zeros(5, 5), Matrix{Float64}(I, 5, 5), M)

sparsify(M, 0.6, 2, 200)

sparsify(M, 0.6, 2, 200, true)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
#this example is modified from https://numericalenvironmental.wordpress.com/2018/05/26/a-steady-state-variably-saturated-flow-model-in-vertical-cross-section-a-finite-difference-approach-using-julia/
import NLsolve
import NonlinearEquations
import PyPlot
import Random
import Zygote

include("utilities.jl")
include("inputdeck.jl")

@NonlinearEquations.equations exclude=(coords, dirichletnodes, neighbors, areasoverlengths) function reqnumenv(psi, Ks, neighbors, areasoverlengths, dirichletnodes, dirichletpsis, coords, alphas, Ns, Qs)
	NonlinearEquations.setnumequations(length(psi))
	for i = 1:length(psi)
		if i in dirichletnodes
			NonlinearEquations.addterm(i, psi[i] - dirichletpsis[i])
		else
			NonlinearEquations.addterm(i, Qs[i])
		end
	end
	for (i, (node_a, node_b)) in enumerate(neighbors)
		for (node1, node2) in [(node_a, node_b), (node_b, node_a)]
			if !(node1 in dirichletnodes) && !(node2 in dirichletnodes)
				relperm = hm(kr(psi[node1], alphas[i], Ns[i]), kr(psi[node2], alphas[i], Ns[i]))
				NonlinearEquations.addterm(node1, relperm * Ks[i] * (psi[node2] + coords[2, node2] - psi[node1] - coords[2, node1]) * areasoverlengths[i])
			elseif !(node1 in dirichletnodes) && node2 in dirichletnodes
				relperm = hm(kr(psi[node1], alphas[i], Ns[i]), kr(dirichletpsis[node2], alphas[i], Ns[i]))
				NonlinearEquations.addterm(node1, relperm * Ks[i] * (dirichletpsis[node2] + coords[2, node2] - psi[node1] - coords[2, node1]) * areasoverlengths[i])
			end
		end
	end
end

@NonlinearEquations.equations exclude=(coords, dirichletnodes, neighbors, areasoverlengths) function req(psi, Ks, neighbors, areasoverlengths, dirichletnodes, dirichletpsis, coords, alphas, Ns, Qs)
	NonlinearEquations.setnumequations(length(psi))
	for i = 1:length(psi)
		if i in dirichletnodes
			NonlinearEquations.addterm(i, psi[i] - dirichletpsis[i])
		else
			NonlinearEquations.addterm(i, Qs[i])
		end
	end
	for (i, (node_a, node_b)) in enumerate(neighbors)
		for (node1, node2) in [(node_a, node_b), (node_b, node_a)]
			if !(node1 in dirichletnodes) && !(node2 in dirichletnodes)
				NonlinearEquations.addterm(node1, hm(kr(psi[node1], alphas[i], Ns[i]), kr(psi[node2], alphas[i], Ns[i])) * Ks[i] * (psi[node2] + coords[2, node2] - psi[node1] - coords[2, node1]) * areasoverlengths[i])
			elseif !(node1 in dirichletnodes) && node2 in dirichletnodes
				NonlinearEquations.addterm(node1, hm(kr(psi[node1], alphas[i], Ns[i]), kr(dirichletpsis[node2], alphas[i], Ns[i])) * Ks[i] * (dirichletpsis[node2] + coords[2, node2] - psi[node1] - coords[2, node1]) * areasoverlengths[i])
			end
		end
	end
end

psi0 = -ones(size(coords, 2))
psi0[dirichletnodes] = zeros(length(dirichletnodes))

function solveforpsi(Ks, psi0=psi0; doplot=false, donewtonish=true)
	resnumenv(psi) = reqnumenv_residuals(psi, Ks, neighbors, areasoverlengths, dirichletnodes, dirichletpsis, coords, alphas, Ns, Qs)
	jacnumenv(psi) = reqnumenv_psi(psi, Ks, neighbors, areasoverlengths, dirichletnodes, dirichletpsis, coords, alphas, Ns, Qs)
	res(psi) = req_residuals(psi, Ks, neighbors, areasoverlengths, dirichletnodes, dirichletpsis, coords, alphas, Ns, Qs)
	jac(psi) = req_psi(psi, Ks, neighbors, areasoverlengths, dirichletnodes, dirichletpsis, coords, alphas, Ns, Qs)
	j!(J, psi) = req_psi!(J, psi, Ks, neighbors, areasoverlengths, dirichletnodes, dirichletpsis, coords, alphas, Ns, Qs)
	f!(residuals, psi) = copy!(residuals, res(psi))
	if donewtonish
		psi1 = NonlinearEquations.newtonish(resnumenv, jacnumenv, psi0; numiters=10, rate=0.05)
	else
		psi1 = psi0
	end
	psi1 = soln_dbe[:, end]
	#callback(psi1, res(psi1), jac(psi1), 0)
	df = NLsolve.OnceDifferentiable(f!, j!, psi1, res(psi0), jac(psi0))
	nls = NLsolve.nlsolve(df, psi1; show_trace=false, iterations=200, ftol=1e-12)
	if doplot
		callback(nls.zero, res(nls.zero), jac(nls.zero), 0)
	end
	return nls.zero
end

psi = solveforpsi(Ks; doplot=true)

#compute the gradient
Random.seed!(1)
obsnode = rand(1:length(psi))
g(psi, Ks) = psi[obsnode]
function g_psi(psi, p)
	retval = zeros(length(psi))
	retval[obsnode] = 1.0
	return retval
end
g_Ks(psi, Ks) = zeros(length(Ks))
f_psi(psi, Ks) = req_psi(psi, Ks, neighbors, areasoverlengths, dirichletnodes, dirichletpsis, coords, alphas, Ns, Qs)
f_Ks(psi, Ks) = req_Ks(psi, Ks, neighbors, areasoverlengths, dirichletnodes, dirichletpsis, coords, alphas, Ns, Qs)
grad = NonlinearEquations.gradient(psi, Ks, g_psi, g_Ks, f_psi, f_Ks)
plotgradient(grad, obsnode)

#=
fdgrad = similar(grad)
sortedgradindices = sort(1:length(grad); by=i->abs(grad[i]), rev=true)
@time for i = 1:length(fdgrad)
	global fdgrad
	global Ks
	global psi
	dk = 1e-8
	psi0 = psi
	Ks0 = copy(Ks)
	thisKs = copy(Ks)
	thisKs[i] += dk
	thispsi = solveforpsi(thisKs, psi0; donewtonish=false)
	fdgrad[i] = (thispsi[obsnode] - psi0[obsnode]) / dk
end
plotgradient(fdgrad, obsnode)
=#


text:
function longest_maximal_repeated_no_overlap_all(text, nodes, root, minlen = nothing)

    if isnothing(minlen)
        minlen = 2
    end

    @assert minlen >= 1

    lenlongest = nothing
    ids = []

    #strlongest = Dict{Any,Array{Integer,1}}()
    strlongest = Dict()

    tovisit = [nd for nd in nodes if !isnothing(nd.children) || length(nd.at) > 1]

    for nd in tovisit

        if nd == root
            continue
        end

        repnode = length(nd.at) > 1
        if repnode
            @assert False # should not happen now as all strings differ
            mylen = nd.pflen + self.lengt - nd.tail

        else
            if length(nd.children) > 0
                maxcplen = maximum(child.cplen for child in nd.children)
            else
                maxcplen = 0
            end
            mylen = nd.pflen + maxcplen
        end

        if mylen < minlen
            continue
        end

        if !isnothing(lenlongest) && mylen < lenlongest

            continue
        end

        loci = collect(nd.at)

        found = false
        for i in range(mylen, 1, step = -1)

            newkids = [child for child in nd.children if child.cplen == i]

            newloci = []
            for child in newkids
                append!(newloci, child.at)
            end

            append!(loci, newloci)
            sort!(loci)

            p0 = loci[1]
            p1 = loci[end]
            nonlap = p1 - p0 >= i
            if nonlap

                mylen = i
                found = true
                break
            end

        end

        if !found
            continue  # next node
        end

        if mylen < minlen
            continue
        end

        if isnothing(lenlongest) || mylen > lenlongest
            lenlongest = mylen
            ids = [nd]
            empty!(strlongest)
            key = text[nd.sorg:nd.sorg+lenlongest-1]

            strlongest[key] = collect(loci)

        elseif mylen == lenlongest
            push!(ids, nd)
            key = text[nd.sorg:nd.sorg+lenlongest-1]
            if haskey(strlongest, key)
                append!(strlongest[key], loci)
            else
                strlongest[key] = collect(loci)
            end
        else
            # nd is smaller
        end
    end

    return (lenlongest, strlongest)
end


text:
function calImplied∂TEMP∂t!(
    ocn::Ocean;
    cfgs...
)

    swflx    = ocn.in_flds.swflx
    nswflx   = ocn.in_flds.nswflx
    frwflx   = ocn.in_flds.frwflx
    qflx_T   = ocn.in_flds.qflx_T

    seaice_nudge_energy = ocn.seaice_nudge_energy

    TFLUX_DIV_implied  = ocn.TFLUX_DIV_implied
    qflx2atm           = ocn.qflx2atm
    TFLUX_bot          = ocn.TFLUX_bot
    TSAS_clim          = ocn.TSAS_clim
    dTEMPdt            = ocn.dTEMPdt

    @loop_hor ocn i j let
        TFLUX_DIV_implied[i, j] =  ( - ( nswflx[i, j] + swflx[i, j] ) + max(qflx2atm[i, j], 0.0) ) / ρc_sw + TSAS_clim[i, j] + TFLUX_bot[i, j] - dTEMPdt[i, j]
    end

    if cfgs[:do_qflx]
        @loop_hor ocn i j let
            TFLUX_DIV_implied[i, j] +=  qflx_T[i, j] / ρc_sw
        end
    end

    if cfgs[:do_seaice_nudging]
        @loop_hor ocn i j let
            TFLUX_DIV_implied[i, j] +=  seaice_nudge_energy[i, j] / ρc_sw
        end
    end

end

function calImplied∂SALT∂t!(
    ocn::Ocean;
    cfgs...
)

    SFLUX_DIV_implied  = ocn.SFLUX_DIV_implied
    SFLUX_top             = ocn.SFLUX_top
    SFLUX_bot             = ocn.SFLUX_bot
    SSAS_clim          = ocn.SSAS_clim
    dSALTdt            = ocn.dSALTdt

    qflx_S   = ocn.in_flds.qflx_S

    @loop_hor ocn i j let
        SFLUX_DIV_implied[i, j] = SSAS_clim[i, j] - SFLUX_top[i, j] + SFLUX_bot[i, j] - dSALTdt[i, j]
    end

    if cfgs[:do_qflx]
        @loop_hor ocn i j let
            SFLUX_DIV_implied[i, j] += qflx_S[i, j]
        end
    end

end

function calDirect∂TEMP∂t!(ocn::Ocean; Δt::Float64)

    @loop_hor ocn i j let

        tmp_TEMP = OC_getIntegratedTemperature(ocn, i, j; target_z = ocn.cols.zs[i, j][ocn.Nz[i, j]+1])
        ocn.TEMP[i, j], ocn.dTEMPdt[i, j] = tmp_TEMP, (tmp_TEMP - ocn.TEMP[i, j]) / Δt

    end

end

function calTEMP!(ocn::Ocean)

    @loop_hor ocn i j let
        ocn.TEMP[i, j] = OC_getIntegratedTemperature(ocn, i, j; target_z = ocn.cols.zs[i, j][ocn.Nz[i, j]+1])
    end
end

function calDirect∂SALT∂t!(ocn::Ocean; Δt::Float64)

    @loop_hor ocn i j let

        tmp_SALT = OC_getIntegratedSalinity(ocn, i, j; target_z = ocn.cols.zs[i, j][ocn.Nz[i, j]+1])
        ocn.SALT[i, j], ocn.dSALTdt[i, j] = tmp_SALT, (tmp_SALT - ocn.SALT[i, j]) / Δt

    end
i
end

function calSALT!(ocn::Ocean)

    @loop_hor ocn i j let

        ocn.SALT[i, j] = OC_getIntegratedSalinity(ocn, i, j; target_z = ocn.cols.zs[i, j][ocn.Nz[i, j]+1])

    end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
#=
def linear(a, b, x):
    return b + a*x
=#

linear(a, b, x) = b + a * x

#=
# a linear demand function is generated for every
# pair of coefficients in vectors a_vec and b_vec
def demand_hypotheses(a_vec, b_vec):
    for a, b in itertools.product(a_vec, b_vec):
        yield {
            'd': functools.partial(linear, a, b),
            'p_opt': -b/(2*a)
        }
=#

function demand_hypothesis(f, a, b)
	f1(x) = f(a, b, x)
	return DataFrame(
		:a => a,
		:b => b,
		:d => f1,
        :d_opt => f1(-b / (2a)),
		:p_opt => -b / (2a)
	)
end

function generate_demand_hypothesis(a_range, b_range)
    h_vec = DataFrame()
    for a in a_range
    	for b in b_range
    		df1 = demand_hypothesis(linear, a, b)
    		push!(h_vec, df1[1, :])
    	end
    end
    h_vec
end


text:
module NNLS

using LinearAlgebra
using Distributed

export nnls,
       nnls!,
       NNLSWorkspace,
       load!

"""
CONSTRUCTION AND/OR APPLICATION OF A SINGLE
HOUSEHOLDER TRANSFORMATION..     Q = I + U*(U**T)/B

The original version of this code was developed by
Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
1973 JUN 12, and published in the book
"SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
Revised FEB 1995 to accompany reprinting of the book by SIAM.
"""
function construct_householder!(u::AbstractVector{T}, up::T)::T where {T}
    m = length(u)
    if m <= 1
        return up
    end

    cl = maximum(abs, u)
    @assert cl > 0
    clinv = 1 / cl
    sm = zero(T)
    @inbounds for ui in u
        sm += (ui * clinv)^2
    end
    cl *= sqrt(sm)
    if u[1] > 0
        cl = -cl
    end
    result = u[1] - cl
    u[1] = cl

    return result
end

"""
CONSTRUCTION AND/OR APPLICATION OF A SINGLE
HOUSEHOLDER TRANSFORMATION..     Q = I + U*(U**T)/B

The original version of this code was developed by
Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
1973 JUN 12, and published in the book
"SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
Revised FEB 1995 to accompany reprinting of the book by SIAM.
"""
function apply_householder!(u::AbstractVector{T}, up::T, c::AbstractVector{T}) where {T}
    m = length(u)
    @inbounds if m > 1
        cl = abs(u[1])
        @assert cl > 0
        b = up * u[1]
        if b >= 0
            return
        end
        b = 1 / b

        sm = c[1] * up
        @simd for i in 2:m
            @inbounds sm += c[i] * u[i]
        end
        @inbounds if sm != 0
            sm *= b
            c[1] += sm * up
            @simd for i in 2:m
                @inbounds c[i] += sm * u[i]
            end
        end
    end
end

"""
   COMPUTE ORTHOGONAL ROTATION MATRIX..
The original version of this code was developed by
Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
1973 JUN 12, and published in the book
"SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
Revised FEB 1995 to accompany reprinting of the book by SIAM.

   COMPUTE.. MATRIX   (C, S) SO THAT (C, S)(A) = (SQRT(A**2+B**2))
                      (-S,C)         (-S,C)(B)   (   0          )
   COMPUTE SIG = SQRT(A**2+B**2)
      SIG IS COMPUTED LAST TO ALLOW FOR THE POSSIBILITY THAT
      SIG MAY BE IN THE SAME LOCATION AS A OR B .
"""
function orthogonal_rotmat(a::T, b::T)::Tuple{T, T, T} where {T}
    if abs(a) > abs(b)
        xr = b / a
        yr = sqrt(1 + xr^2)
        c = (1 / yr) * sign(a)
        s = c * xr
        sig = abs(a) * yr
    elseif b != 0
        xr = a / b
        yr = sqrt(1 + xr^2)
        s = (1 / yr) * sign(b)
        c = s * xr
        sig = abs(b) * yr
    else
        sig = zero(T)
        c = zero(T)
        s = one(T)
    end
    return c, s, sig
end

"""
The original version of this code was developed by
Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
1973 JUN 15, and published in the book
"SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
Revised FEB 1995 to accompany reprinting of the book by SIAM.
"""
function solve_triangular_system!(zz, A, idx, nsetp, jj)
    @inbounds for l in 1:nsetp
        ip = nsetp + 1 - l
        if (l != 1)
            @simd for ii in 1:ip
                @inbounds zz[ii] -= A[ii, jj] * zz[ip + 1]
                # @inbounds zz[ii] = muladd(-A[ii, jj], zz[ip + 1], zz[ii])
            end
        end
        jj = idx[ip]
        zz[ip] /= A[ip, jj]
    end
    return jj
end

mutable struct NNLSWorkspace{T, I <: Integer}
    QA::Matrix{T}
    Qb::Vector{T}
    x::Vector{T}
    w::Vector{T}
    zz::Vector{T}
    idx::Vector{I}
    rnorm::T
    mode::I
    nsetp::I
end

function NNLSWorkspace{T,I}(m, n) where {T, I<:Integer}
    NNLSWorkspace{T,I}(
        zeros(T, m, n), # A
        zeros(T, m),    # b
        zeros(T, n),    # x
        zeros(T, n),    # w
        zeros(T, m),    # zz
        zeros(I, n),    # idx
        zero(T),        # rnorm
        zero(I),        # mode
        zero(I),        # nsetp
    )
end

function Base.resize!(work::NNLSWorkspace{T}, m::Integer, n::Integer) where {T}
    work.QA = zeros(T, m, n)
    work.Qb = zeros(T, m)
    resize!(work.x, n)
    resize!(work.w, n)
    resize!(work.zz, m)
    resize!(work.idx, n)
end

function load!(work::NNLSWorkspace{T}, A::AbstractMatrix{T}, b::AbstractVector{T}) where {T}
    m, n = size(A)
    @assert size(b) == (m,)
    if size(work.QA, 1) != m || size(work.QA, 2) != n
        resize!(work, m, n)
    end
    work.QA .= A
    work.Qb .= b
    work
end

function NNLSWorkspace(
        m::Integer,
        n::Integer,
        eltype::Type{T} = Float64,
        indextype::Type{I} = Int
    ) where {T,I}
    NNLSWorkspace{T, I}(m, n)
end

function NNLSWorkspace(A::Matrix{T}, b::Vector{T}, indextype::Type{I} = Int) where {T,I}
    m, n = size(A)
    @assert size(b) == (m,)
    work = NNLSWorkspace{T, I}(m, n)
    load!(work, A, b)
    work
end

"""
Views in Julia still allocate some memory (since they need to keep
a reference to the original array). This type allocates no memory
and does no bounds checking. Use it with caution.
"""
struct UnsafeVectorView{T} <: AbstractVector{T}
    offset::Int
    len::Int
    ptr::Ptr{T}
end

UnsafeVectorView(parent::DenseArray{T}, start_ind::Integer, len::Integer) where {T} = UnsafeVectorView{T}(start_ind - 1, len, pointer(parent))
Base.size(v::UnsafeVectorView) = (v.len,)
Base.getindex(v::UnsafeVectorView, idx) = unsafe_load(v.ptr, idx + v.offset)
Base.setindex!(v::UnsafeVectorView, value, idx) = unsafe_store!(v.ptr, value, idx + v.offset)
Base.length(v::UnsafeVectorView) = v.len
@static if VERSION >= v"0.7-"
    Base.IndexStyle(::Type{V}) where {V <: UnsafeVectorView} = Base.IndexLinear()
else
    Base.IndexStyle{V <: UnsafeVectorView}(::Type{V}) = Base.IndexLinear()
end

"""
UnsafeVectorView only works for isbits types. For other types, we're already
allocating lots of memory elsewhere, so creating a new View is fine.

This function looks type-unstable, but the isbitstype(T) test can be evaluated
by the compiler, so the result is actually type-stable.
"""
function fastview(parent::DenseArray{T}, start_ind::Integer, len::Integer) where {T}
    if isbitstype(T)
        UnsafeVectorView(parent, start_ind, len)
    else
        @view(parent[start_ind:(start_ind + len - 1)])
    end
end

@noinline function checkargs(work::NNLSWorkspace)
    m, n = size(work.QA)
    @assert size(work.Qb) == (m,)
    @assert size(work.x) == (n,)
    @assert size(work.w) == (n,)
    @assert size(work.zz) == (m,)
    @assert size(work.idx) == (n,)
end

function largest_positive_dual(
        w::AbstractVector{T},
        idx::AbstractVector{TI},
        range
    ) where {T,TI}
    wmax = zero(T)
    izmax = zero(TI)
    @inbounds for i in range
        j = idx[i]
        if w[j] > wmax
            wmax = w[j]
            izmax = i
        end
    end
    wmax, izmax
end

"""
Algorithm NNLS: NONNEGATIVE LEAST SQUARES

The original version of this code was developed by
Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
1973 JUN 15, and published in the book
"SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
Revised FEB 1995 to accompany reprinting of the book by SIAM.

GIVEN AN M BY N MATRIX, A, AND AN M-VECTOR, B,  COMPUTE AN
N-VECTOR, X, THAT SOLVES THE LEAST SQUARES PROBLEM
                 A * X = B  SUBJECT TO X .GE. 0
"""
function nnls!(
        work::NNLSWorkspace{T, TI},
        max_iter::Integer = 3*size(work.QA, 2)
    ) where {T, TI}

    checkargs(work)

    A = work.QA
    Ainds = LinearIndices(A)
    b = work.Qb
    x = work.x
    w = work.w
    zz = work.zz
    idx = work.idx
    factor = T(0.01)
    work.mode = 1

    m = convert(TI, size(A, 1))
    n = convert(TI, size(A, 2))

    iter = 0
    x .= 0
    idx .= 1:n

    iz2 = n
    iz1 = one(TI)
    iz = zero(TI)
    j = zero(TI)
    jj = zero(TI)
    nsetp = zero(TI)
    up = zero(T)

    terminated = false

    # ******  MAIN LOOP BEGINS HERE  ******
    @inbounds while true
        # println("jl main loop")
        # QUIT IF ALL COEFFICIENTS ARE ALREADY IN THE SOLUTION.
        # OR IF M COLS OF A HAVE BEEN TRIANGULARIZED.
        if (iz1 > iz2 || nsetp >= m)
            terminated = true
            break
        end

        # COMPUTE COMPONENTS OF THE DUAL (NEGATIVE GRADIENT) VECTOR W().
        @inbounds for i in iz1:iz2
            idxi = idx[i]
            sm = zero(T)
            @simd for l in (nsetp + 1):m
                @inbounds sm += A[l, idxi] * b[l]
            end
            w[idxi] = sm
        end

        @inbounds while true
            # FIND LARGEST POSITIVE W(J).
            wmax, izmax = largest_positive_dual(w, idx, iz1:iz2)

            # IF WMAX .LE. 0. GO TO TERMINATION.
            # THIS INDICATES SATISFACTION OF THE KUHN-TUCKER CONDITIONS.
            if wmax <= 0
                terminated = true
                break
            end

            iz = izmax
            j = idx[iz]

            # THE SIGN OF W(J) IS OK FOR J TO BE MOVED TO SET P.
            # BEGIN THE TRANSFORMATION AND CHECK NEW DIAGONAL ELEMENT TO AVOID
            # NEAR LINEAR DEPENDENCE.
            Asave = A[nsetp + 1, j]
            up = construct_householder!(
                 fastview(A, Ainds[nsetp + 1, j], m - nsetp), up)
            unorm = zero(T)
            @simd for l in 1:nsetp
                @inbounds unorm += A[l, j]^2
            end
            unorm = sqrt(unorm)

            if ((unorm + abs(A[nsetp + 1, j]) * factor) - unorm) > 0
                # COL J IS SUFFICIENTLY INDEPENDENT.  COPY B INTO ZZ, UPDATE ZZ
                # AND SOLVE FOR ZTEST ( = PROPOSED NEW VALUE FOR X(J) ).
                # println("copying b into zz")
                zz .= b
                apply_householder!(
                    fastview(A, Ainds[nsetp + 1, j], m - nsetp),
                    up,
                    fastview(zz, nsetp + 1, m - nsetp))
                ztest = zz[nsetp + 1] / A[nsetp + 1, j]

                # SEE IF ZTEST IS POSITIVE
                if ztest > 0
                    break
                end
            end

            # REJECT J AS A CANDIDATE TO BE MOVED FROM SET Z TO SET P.
            # RESTORE A(NPP1,J), SET W(J)=0., AND LOOP BACK TO TEST DUAL
            # COEFFS AGAIN.
            A[nsetp + 1, j] = Asave
            w[j] = 0
        end
        if terminated
            break
        end

        # THE INDEX  J=INDEX(IZ)  HAS BEEN SELECTED TO BE MOVED FROM
        # SET Z TO SET P.    UPDATE B,  UPDATE INDICES,  APPLY HOUSEHOLDER
        # TRANSFORMATIONS TO COLS IN NEW SET Z,  ZERO SUBDIAGONAL ELTS IN
        # COL J,  SET W(J)=0.
        b .= zz

        idx[iz] = idx[iz1]
        idx[iz1] = j
        iz1 += one(TI)
        nsetp += one(TI)

        if iz1 <= iz2
            @inbounds for jz in iz1:iz2
                jj = idx[jz]
                apply_householder!(
                    fastview(A, Ainds[nsetp, j], m - nsetp + 1),
                    up,
                    fastview(A, Ainds[nsetp, jj], m - nsetp + 1))
            end
        end

        if nsetp != m
            @inbounds for l in (nsetp + 1):m
                A[l, j] = 0
            end
        end

        w[j] = 0

        # SOLVE THE TRIANGULAR SYSTEM.
        # STORE THE SOLUTION TEMPORARILY IN ZZ().
        jj = solve_triangular_system!(zz, A, idx, nsetp, jj)

        # ******  SECONDARY LOOP BEGINS HERE ******
        #
        # ITERATION COUNTER.
        @inbounds while true
            iter += 1
            if iter > max_iter
                work.mode = 3
                terminated = true
                println("NNLS quitting on iteration count")
                break
            end

            # SEE IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE.
            # IF NOT COMPUTE ALPHA.
            alpha = convert(T, 2)
            @inbounds for ip in one(TI):nsetp
                l = idx[ip]
                if zz[ip] <= 0
                    t = -x[l] / (zz[ip] - x[l])
                    if alpha > t
                        alpha = t
                        jj = ip
                    end
                end
            end

            # IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE THEN ALPHA WILL
            # STILL = 2.    IF SO EXIT FROM SECONDARY LOOP TO MAIN LOOP.
            if alpha == 2
                break
            end

            # OTHERWISE USE ALPHA WHICH WILL BE BETWEEN 0 AND 1 TO
            # INTERPOLATE BETWEEN THE OLD X AND THE NEW ZZ.
            @inbounds for ip in one(TI):nsetp
                l = idx[ip]
                x[l] = x[l] + alpha * (zz[ip] - x[l])
            end

            # MODIFY A AND B AND THE INDEX ARRAYS TO MOVE COEFFICIENT I
            # FROM SET P TO SET Z.
            i = idx[jj]
            kk = 1
            @inbounds while true
                x[i] = 0

                if jj != nsetp
                    jj += one(TI)
                    @inbounds for j in jj:nsetp
                        ii = idx[j]
                        idx[j - 1] = ii
                        cc, ss, sig = orthogonal_rotmat(A[j - 1, ii], A[j, ii])
                        A[j - 1, ii] = sig
                        A[j, ii] = 0
                        @inbounds for l in one(TI):n
                            if l != ii
                                # Apply procedure G2 (CC,SS,A(J-1,L),A(J,L))
                                temp = A[j - 1, l]
                                A[j - 1, l] = cc * temp + ss * A[j, l]
                                A[j, l] = -ss * temp + cc * A[j, l]
                            end
                        end

                        # Apply procedure G2 (CC,SS,B(J-1),B(J))
                        temp = b[j - 1]
                        b[j - 1] = cc * temp + ss * b[j]
                        b[j] = -ss * temp + cc * b[j]
                    end
                end

                nsetp -= one(TI)
                iz1 -= one(TI)
                idx[iz1] = i

                # SEE IF THE REMAINING COEFFS IN SET P ARE FEASIBLE.  THEY SHOULD
                # BE BECAUSE OF THE WAY ALPHA WAS DETERMINED.
                # IF ANY ARE INFEASIBLE IT IS DUE TO ROUND-OFF ERROR.  ANY
                # THAT ARE NONPOSITIVE WILL BE SET TO ZERO
                # AND MOVED FROM SET P TO SET Z.
                allfeasible = true
                @inbounds for jj in one(TI):nsetp
                    i = idx[jj]
                    if x[i] <= 0
                        allfeasible = false
                        kk = jj
                        break
                    end
                end
                if allfeasible
                    break
                end
            end

            # COPY B( ) INTO ZZ( ).  THEN SOLVE AGAIN AND LOOP BACK.
            zz .= b
            jj = solve_triangular_system!(zz, A, idx, nsetp, kk)
        end
        if terminated
            break
        end
        # ******  END OF SECONDARY LOOP  ******

        @inbounds for i in 1:nsetp
            x[idx[i]] = zz[i]
        end
        # ALL NEW COEFFS ARE POSITIVE.  LOOP BACK TO BEGINNING.
    end

    # ******  END OF MAIN LOOP  ******
    # COME TO HERE FOR TERMINATION.
    # COMPUTE THE NORM OF THE FINAL RESIDUAL VECTOR.

    sm = zero(T)
    if nsetp < m
        @simd for i in (nsetp + 1):m
            @inbounds sm += b[i]^2
        end
    else
        w .= 0
    end
    work.rnorm = sqrt(sm)
    work.nsetp = nsetp
    return work.x
end

function nnls!(
        work::NNLSWorkspace{T},
        A::AbstractMatrix{T},
        b::AbstractVector{T},
        max_iter = 3*size(A, 2)
    ) where {T}
    load!(work, A, b)
    nnls!(work, max_iter)
    work.x
end

"""
x = nnls(A, b; ...)

Solves non-negative least-squares problem by the active set method
of Lawson & Hanson (1974).

Optional arguments:
    max_iter: maximum number of iterations (counts inner loop iterations)

References:
    Lawson, C.L. and R.J. Hanson, Solving Least-Squares Problems,
    Prentice-Hall, Chapter 23, p. 161, 1974.
"""
function nnls(
        A,
        b::AbstractVector{T};
        max_iter::Int = 3*size(A, 2)
    ) where {T}
    work = NNLSWorkspace(A, b)
    nnls!(work, max_iter)
    work.x
end

function nnls(
        A,
        B::AbstractMatrix{T};
        use_parallel = false,
        max_iter::Int = 3*size(A, 2)
    ) where {T}

    m, n = size(A)
    k = size(B, 2)

    if k > 1 && use_parallel && nprocs() > 1
        X = @distributed (hcat) for i = 1:k
            nnls(A, @view(B[:,i]); max_iter=max_iter)
        end
    else
        work = NNLSWorkspace(m, n, T)
        X = zeros(T, n, k)
        @inbounds for i = 1:k
            X[:, i] = nnls!(work, A, @view(B[:,i]), max_iter)
        end
        return X
    end
end

end # module


text:
using Distributed
@everywhere begin
	import DPFEHM
	import GaussianRandomFields
	import Optim
	import PyPlot
	import Random
	import Zygote

	function inversemodel(sigma, lambda, mu, deltah, sqrtnumobs, alpha, beta, gamma, iterations)
		#sigma is the standard deviation of logK
		#lambda is the correlation length
		#mu is the mean of logK
		#deltah is the pressure drop from the left boundary to the right boundary
		#sqrtnumobs is the square root of the number of observation wells (which are distributed on a regular sqrtnumobs-by-sqrtnumobs grid)
		#alpha is the weight associated with head measurements
		#beta is the weight associated with the logK measurements
		#gamma is the weight associated with the darcy velocity measurements
		#iterations is the maximum number of iterations that the inverse analysis will perform

		function getobsnodes(coords, obslocs)
			obsnodes = Array{Int}(undef, length(obslocs))
			for i = 1:length(obslocs)
				obsnodes[i] = findmin(map(j->sum((obslocs[i] .- coords[:, j]) .^ 2), 1:size(coords, 2)))[2]
			end
			return obsnodes
		end

		mins = [0.0, 0.0]#meters
		maxs = [100.0, 100.0]#meters
		ns = [101, 101]
		dx = (maxs[1] - mins[1]) / (ns[1] - 1)
		dy = (maxs[2] - mins[2]) / (ns[2] - 1)
		num_eigenvectors = 200
		x0 = zeros(num_eigenvectors)

		eigenvalue_error = fill(NaN, num_eigenvectors) #RMS of eigenvalues
		logK_error = fill(NaN, iterations + 1) #RMS of logK
		obj_func_values = fill(NaN, iterations + 1) #RMS of Loss

		#set up the eigenvector parameterization of the geostatistical field
		cov = GaussianRandomFields.CovarianceFunction(2, GaussianRandomFields.Matern(lambda, 1; σ=sigma))
		x_pts = range(mins[1], maxs[1]; length=ns[1])
		y_pts = range(mins[2], maxs[2]; length=ns[2])
		grf = GaussianRandomFields.GaussianRandomField(cov, GaussianRandomFields.KarhunenLoeve(num_eigenvectors), x_pts, y_pts)
		logKs = GaussianRandomFields.sample(grf)
		parameterization = copy(grf.data.eigenfunc)
		sigmas = copy(grf.data.eigenval)

		#set up the mesh  and boundary conditions
		coords, neighbors, areasoverlengths, _ = DPFEHM.regulargrid2d(mins, maxs, ns, 1.0)
		Qs = zeros(size(coords, 2))
		boundaryhead(x, y) = deltah * (x - maxs[1]) / (mins[1] - maxs[1])
		dirichletnodes = Int[]
		dirichleths = zeros(size(coords, 2))
		for i = 1:size(coords, 2)
			if coords[1, i] == mins[1] || coords[1, i] == maxs[1] || coords[2, i] == mins[2] || coords[2, i] == maxs[2]
				push!(dirichletnodes, i)
				dirichleths[i] = boundaryhead(coords[1:2, i]...)
			end
		end

		logKs2Ks_neighbors(Ks) = exp.(0.5 * (Ks[map(p->p[1], neighbors)] .+ Ks[map(p->p[2], neighbors)]))#Zygote differentiates this efficiently but the definitions above are ineffecient with Zygote
		function solveforh(logKs, dirichleths)
			@assert length(logKs) == length(Qs)
			if maximum(logKs) - minimum(logKs) > 25
				return fill(NaN, length(Qs))#this is needed to prevent the solver from blowing up if the line search takes us somewhere crazy
			else
				Ks_neighbors = logKs2Ks_neighbors(logKs)
				return DPFEHM.groundwater_steadystate(Ks_neighbors, neighbors, areasoverlengths, dirichletnodes, dirichleths, Qs)
			end
		end
		solveforheigs(x) = solveforh(x2logKs(x), dirichleths)
		x2logKs(x) = reshape(parameterization * (sigmas .* x), ns...)

		#set up the observations
		obslocs_x = range(mins[1], maxs[1]; length=sqrtnumobs + 2)[2:end - 1]
		obslocs_y = range(mins[2], maxs[2]; length=sqrtnumobs + 2)[2:end - 1]
		obslocs = collect(Iterators.product(obslocs_x, obslocs_y))[:]
		x_true = randn(num_eigenvectors)
		logKs_true = x2logKs(x_true)
		h_true = solveforheigs(x_true)
		obsnodes = getobsnodes(coords, obslocs)
		h_obser = h_true[obsnodes]
		logK_obser = logKs_true[obsnodes]
		hmat = reshape(h_true, reverse(ns)...)
		Ks = exp.(logKs)
		vy = Ks .* vcat(fill(0.0, ns[1])', hmat[3:end, :] - hmat[1:end - 2, :], fill(0.0, ns[1])') / (2 * dy)
		vx = Ks .* hcat(fill(0.0, ns[2]), hmat[:, 3:end] - hmat[:, 1:end - 2], fill(0.0, ns[2])) / (2 * dx)
		vel_obser = map(i->sqrt(vx[i] ^ 2 + vy[i] ^ 2), obsnodes)

		#set up the loss function
		function f(logKs)
			h_obssigma = 1e-3
			logK_obssigma = 1e-1
			vel_obssigma = 1e-1
			h = solveforh(logKs, dirichleths)
			h_error = sum((h[obsnodes] - h_obser) .^ 2 ./ h_obssigma ^ 2)
			K_error = sum((logKs[obsnodes] - logK_obser) .^ 2 ./ logK_obssigma ^ 2)
			Ks = exp.(logKs)
			hmat = reshape(h, reverse(ns)...)
			vy = Ks .* vcat(fill(0.0, ns[1])', hmat[3:end, :] - hmat[1:end - 2, :], fill(0.0, ns[1])') / (2 * dy)
			vx = Ks .* hcat(fill(0.0, ns[2]), hmat[:, 3:end] - hmat[:, 1:end - 2], fill(0.0, ns[2])) / (2 * dx)
			vel_pred = map(i->sqrt(vx[i] ^ 2 + vy[i] ^ 2), obsnodes)
			vel_error = sum((vel_pred - vel_obser) .^ 2 ./ vel_obssigma ^ 2)
			error_sys = alpha * h_error + beta * K_error + gamma * vel_error
			return error_sys
		end
		feigs(x) = f(x2logKs(x)) + sum(x .^ 2)

		#solve the inverse problem
		options = Optim.Options(iterations=iterations, store_trace=true, show_trace=false, extended_trace=true, x_tol=1e-6, time_limit=60 * 60 * 2)
		opt = Optim.optimize(feigs, x->Zygote.gradient(feigs, x)[1], x0, Optim.LBFGS(), options; inplace=false)
		logK_error[1:length(opt.trace)] = map(j->(sqrt.((sum((x2logKs(opt.trace[j].metadata["x"])-x2logKs(x_true)).^2))/num_eigenvectors)),1:length(opt.trace))
		obj_func_values[1:length(opt.trace)] = map(j->(sqrt.(opt.trace[j].value)), 1:length(opt.trace))
		x_est = opt.minimizer
		eigenvalue_error .= x_est - x_true
		return eigenvalue_error, logK_error, obj_func_values
	end
end#end the distributed portion


text:
#include("")

```
 long-lived assets, tangible assets[1] or property, plant and equipment
 -[wikipedia](https://en.wikipedia.org/wiki/Fixed_asset)
 me: a fixed asset is used for:
  1. personal use as a use (to store now, & resell later) -
  -at a higher price: specific-case: house(Real-estate):Bull-market in general,
  2. busiess: an investment for day-to-day operations (i.e in work)
  #reason: to pay of its (paid )book price
 for (small) business: any thing that could be use
  E.g: Land, Building,Equipment(PC, printer, Phone,)

```
module mutable struct fixedAsset <: Asset
  ownershipDate:dateTime

  #fields
end

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
module quickPIV

using FFTW, FileIO, LIBTIFF, Statistics

export PIV

# Shortened type names to keep function signatures short ( when possible one-liners )
include("units.jl")

F      = Union{Float64,Float32} # considering 64 and 32 bit systems
I      = Union{Int64,Int32}
A      = Array
S      = String
C      = Complex
II     = Tuple{I,I}
III    = Tuple{I,I,I}
UnI    = Union{I,II,III}
SySySy = Tuple{Symbol,Symbol,Symbol}
DU     = DistanceUnits
TU     = TimeUnits
UNTS   = Tuple{DU,DU,DU,TU}

Base.convert( ::Type{Union{Int32,Int64}}, a::AbstractFloat ) = round( typeof(0), a )

struct ZNCC   end
struct FFT    end
struct NSQECC end

CORRTYPES = Union{ZNCC,FFT,NSQECC};

include("parameters.jl")
include("pivIO.jl")
include("utils.jl")

include("PIV/subpixel.jl")
include("PIV/signalToNoise.jl")
include("PIV/ccrZNCC.jl")
include("PIV/ccrFFT.jl")
include("PIV/ccrNSQECC.jl")

include("syntheticParticles.jl")

include("postprocessing.jl")

function PIV( data1::Array{T,N}, data2::Array{T,N}, params::PIVParameters;
              checkParams=true, units=nothing ) where {T<:AbstractFloat,N}

	checkParams && checkPIVParameters( size(data1), params )

	return _PIV( data1, data2, params )
end

""" Calling 2D PIV analyses """

function _PIV( img1::A{T,2}, img2::A{T,2}, p::PIVParameters; units=nothing
             ) where {T<:AbstractFloat}

	u, v, sn = PIV_2D( p.corr, img1, img2, p.args... );

	if units !== nothing
	    ratios = changeUnits.(  p.units[1:2], units[1:2] )./changeUnits( p.units, units[4]  )
	    for e in 1:length(u)
	        u[e] = u[e]*ratios[1]
	        v[e] = v[e]*ratios[2]
	    end
	end
	return u, v, sn;
end

""" Calling 3D PIV analyses """

function _PIV( vol1::A{T,3}, vol2::A{T,3}, p::PIVParameters; units=nothing
             ) where {T<:AbstractFloat}

	u, v, w, sn = PIV_3D( p.corr, vol1, vol2, p.args... );

	if units !== nothing
	    ratios = changeUnits.( p.units, units[1:3] )./changeUnits( p.units, units[4] )
	    for e in 1:length(u)
	        u[e] = u[e]*ratios[1]
	        v[e] = v[e]*ratios[2]
	        w[e] = w[e]*ratios[3]
	    end
	end

	return u, v, w, sn;
end

end

max_issues_repo_name: ErichZimmer/quickPIV

max_forks_repo_name: ErichZimmer/quickPIV


text:
# AS wrapper.

mutable struct AsNode <: AbstractSQLNode
    over::Union{SQLNode, Nothing}
    name::Symbol

    AsNode(;
           over = nothing,
           name::Union{Symbol, AbstractString}) =
        new(over, Symbol(name))
end

AsNode(name; over = nothing) =
    AsNode(over = over, name = name)

"""
    As(; over = nothing, name)
    As(name; over = nothing)
    name => over

In a scalar context, `As` specifies the name of the output column.  When
applied to a tabular node, `As` wraps the output of the node in a nested
record.

# Examples

```jldoctest
julia> person = SQLTable(:person, columns = [:person_id, :year_of_birth]);

julia> q = From(person) |>
           As(:p) |>
           Select(:birth_year => Get.p.year_of_birth);
```
"""
As(args...; kws...) =
    AsNode(args...; kws...) |> SQLNode

dissect(scr::Symbol, ::typeof(As), pats::Vector{Any}) =
    dissect(scr, AsNode, pats)

Base.convert(::Type{AbstractSQLNode}, p::Pair{<:Union{Symbol, AbstractString}}) =
    AsNode(name = first(p), over = convert(SQLNode, last(p)))

function PrettyPrinting.quoteof(n::AsNode, ctx::QuoteContext)
    ex = Expr(:call, nameof(As), quoteof(n.name))
    if n.over !== nothing
        ex = Expr(:call, :|>, quoteof(n.over, ctx), ex)
    end
    ex
end

label(n::AsNode) =
    n.name

rebase(n::AsNode, n′) =
    AsNode(over = rebase(n.over, n′), name = n.name)