text:
"""
Module for symbolic transformation of models.

* Developer: Hilding Elmqvist, Mogram AB
* First version: July-September 2016
* License: MIT (expat)

Reference:
For matrix formulas, see:
Petersen, K.B, Pedersen M.S.: The Matrix Cookbook
https://www.math.uwaterloo.ca/~hwolkowi/matrixcookbook.pdf

"""
module SymbolicTransform

#using Debug
using ..Modia
using ..Instantiation
import ..Instantiation: GetField, This, Der, Symbolic, time_global, simulationModel_symbol
using Base.Meta: quot, isexpr
using DataStructures
using Unitful
@static if !(VERSION < v"0.7.0-DEV.2005")
    import LinearAlgebra
end
@static if VERSION < v"0.7.0-DEV.2005"
    notFound = 0
else
    notFound = nothing
end

using ..Synchronous
using ..BLTandPantelidesUtilities
#using ..Utilities
import ModiaMath
using ..ModiaLogging

export initSynchronousCounters, substituteClocked, solve, solveEquations, solveSortedEquations, differentiate, differentiateEquations, differentiateSortedEquations, residue_der, newDifferentiateEquations
export setDerAsFunction, derAsFunction, setTimeInvariants, setRealStates

const logSolve = false
const logDifferentiate = true
const logDifferentiateVariable = false
const logSolved = false
const noUnits = true
const pushResiduals = false
const checkUnits = false

#=
@show logSolve
@show logSolved
=#

global derAsFunction = true
setDerAsFunction(v) = (global derAsFunction = v)

global nClocks
global nsamples
global nprevious
global npositive

function initSynchronousCounters()
    global nClocks = 0
    global nsamples = 0
    global nprevious = 0
    global npositive = 0
end

substituteClocked(previousVariables, ex) = ex

function substituteClocked(previousVariables, ex::Expr)
    global npositive
    if isexpr(ex, :quote)
        ex
    elseif isexpr(ex, :call) && ex.args[1] == Synchronous.Clock
        global nClocks += 1
        Expr(ex.head, ex.args[1], ex.args[2], simulationModel_symbol, nClocks)
    elseif isexpr(ex, :call) && ex.args[1] == Synchronous.sample
        global nsamples += 1
        Expr(ex.head, ex.args[1], ex.args[2], substituteClocked(previousVariables, ex.args[3]), simulationModel_symbol, nsamples)
    elseif isexpr(ex, :call) && ex.args[1] == Synchronous.previous
        global nprevious += 1
        push!(previousVariables, ex.args[2])
        Expr(ex.head, ex.args[1], ex.args[2], substituteClocked(previousVariables, ex.args[3]), simulationModel_symbol, nprevious)
    elseif isexpr(ex, :call) && ex.args[1] == Synchronous.positive
        global npositive += 1
        Expr(ex.head, ModiaMath.ModiaToModiaMath.positive!, ex.args[2], simulationModel_symbol, npositive)
    else
        Expr(ex.head, [substituteClocked(previousVariables, arg) for arg in ex.args]...)
    end
end

zero = Float64(0.0)
one = Float64(1.0)

function add(e1, e2)
    if e1 == zero
        return e2
    elseif e2 == zero
        return e1
    else
        Expr(:call, +, e1, e2)
    end
end

function sub(e1, e2)
    if typeof(e1) in [Float64, Int64] && typeof(e2) in [Float64, Int64]
        e1 - e2
    elseif e1 == zero
        return Expr(:call, -, e2)
    elseif e2 == zero
        return e1
    else
        Expr(:call, -, e1, e2)
    end
end

function mult(e1, e2)
    if e1 == zero || e2 == zero
        return zero
    elseif e1 == one
        return e2
    elseif e2 == one
        return e1
    else
        Expr(:call, *, e1, e2)
    end
end

function power(e1, e2)
    if e2 == one
        return e1
    else
        Expr(:call, ^, e1, e2)
    end
end

# operators = [call, +, :+, -, :-, *, :*, (==), ^, :^, :(=), \, /, :/, transpose]
operators = [:call, +, :+, -, :-, *, :*, (==), ^, :^, :(=), \, /, :/, transpose] # v 0.6

checkincidence(s::Function) = s in operators ? nothing : s
checkincidence(e::Instantiation.GetField) = e.name
# checkincidence(e::Instantiation.Der) = e.base
# checkincidence(e::Expr) = map(checkincidence, e.args)
checkincidence(e::Number) = nothing
checkincidence(e::LineNumberNode) = nothing
checkincidence(x) = x

# From: https://rosettacode.org/wiki/Flatten_a_list#Julia
@static if VERSION < v"0.7.0-DEV.2005"
    flat(A) = mapreduce(x -> isa(x, Array) ? flat(x) : x, vcat, [], A)
else
    flat(A) = mapreduce(x -> isa(x, Array) ? flat(x) : x, vcat, A, init=[])
end

function getincidence(e)
    # @show e typeof(e)
    if isa(e, Function)
        incidence = [e in operators ? nothing : e]
    elseif typeof(e) <: Unitful.Quantity
        # dump(e)
        incidence = nothing
    elseif typeof(e) == Instantiation.GetField
        incidence = [e, e.name]
    elseif typeof(e) == Instantiation.Der
        incidence = e # getincidence(e.base)
    elseif typeof(e) == Expr
        if e.head == :call && e.args[1] == :der
            incidence = e
        elseif !derAsFunction && e.head == :call && e.args[1] == Synchronous.previous
            incidence = map(getincidence, e.args[3:end])
        else
            incidence = map(getincidence, e.args[2:end])
        end
    elseif e == time_global || e == :time
        incidence = nothing
    elseif typeof(e) == Symbol
        # if e in [:kg, :m, :s, :N]
            # incidence = nothing
        # else
        incidence = e
        # end
    elseif typeof(e) in [Number, Float64, Int64, LineNumberNode]
        incidence = nothing
    # elseif typeof(e) in [block]
            # incidence = checkincidence(e.args[2])
    elseif e != nothing
        incidence = map(checkincidence, e)
    else
        incidence = nothing
    end

    if isa(incidence, Array)
      ### This allocates a lot of memory and should be rewritten.
        incidence = flat(incidence)
    else
        incidence = [incidence]
    end
    incidence = filter(e -> e != nothing, incidence)
    # println("getincidence: ", incidence)
    return incidence
end

find_incidence!(incidence::Array{Any,1}, ex) = nothing
find_incidence!(incidence::Array{Any,1}, der::Der) = (push!(incidence, der); nothing)
find_incidence!(incidence::Array{Any,1}, get::GetField) = (push!(incidence, get.name); nothing)

function find_incidence!(incidence::Array{Any,1}, ex::Expr)
    if !isexpr(ex, :quote)
        if ex.head == :call
            if ex.args[1] in [(Synchronous.previous), :(Synchronous.previous)]
                for arg in ex.args[3:end]    # 3:end ?
                    find_incidence!(incidence, arg)
                end
            elseif ex.args[1] == :der
                println(":der incidence")
                push!(incidence, ex)
            end
        else
            for arg in ex.args
                find_incidence!(incidence, arg)
            end
        end
    end
    nothing
end

#=
inc = []
find_incidence!(inc, LHS)
if testIncidence(inc, x)
=#

function testIncidence(e, x)
    inc = getincidence(e)
    return findfirst(isequal(x), inc) != notFound
end

# --------------------------------------------------------------

function solve(eq::Expr, x)
    #=
    Solve equation eq for x

    Rewriting rules:
    ex: expression containing x
    e = ex  ->  ex = e
    -ex = e  ->  ex = -e

    ex + e1 = e2  ->  ex = e2 - e1
    e1 + ex = e2  ->  ex = e2 - e1

    ex - e1 = e2  ->  ex = e2 + e1
    e1 - ex = e2  ->  ex = e1 - e2

    ex * e1 = e2  ->  ex = e2 / e1
    e1 * ex = e2  ->  ex = e1 \ e2

    ex / e1 = e2  ->  ex = e2 * e1
    e1 / ex = e2  ->  e1 = e2 * ex  ->  ex = e2 \ e1

    ex \ e1 = e2  ->  e1 = ex * e2  ->  ex = e1 / e2
    e1 \ ex = e2  ->  ex = e1 * e2

    transpose(ex) = e  ->  ex = transpose(e)
    =#

    #  @assert eq.head in [:(=), :(:=)] "solve only solves equations"
    if !(eq.head in [:(=), :(:=)])
        return (eq, false)
    end

    if logSolve
        println("\nSOLVE: ", string(x), " from:   ", prettyPrint(eq))
    end

    LHS = eq.args[1]
    RHS = eq.args[2]

    if testIncidence(LHS, x) && testIncidence(RHS, x)
        println("Multiple occurances of unknown: ", prettyPrint(eq))
        return eq, false
    elseif RHS == x || testIncidence(RHS, x)
        # Unknown x in RHS
        if logSolve println("  right incidence, swap") end
        # Swap LHS and RHS
        eq, solved = solve(Expr(:(=), RHS, LHS), x)
        return eq, solved
    elseif LHS == x || testIncidence(LHS, x)
        # Unknown x in LHS
        if LHS == x || typeof(LHS) == Symbol || typeof(LHS) == Instantiation.GetField || typeof(LHS) == Instantiation.Der
            if logSolve println("  Already solved: ", prettyPrint(eq)) end
            return eq, true
        elseif LHS.head == :block && (isa(LHS.args[1], LineNumberNode) || isexpr(LHS.args[1], :line))
            eq = Expr(:(=), LHS.args[2], RHS)
            eq, solved = solve(eq, x)
            return eq, solved
        end

        op = LHS.args[1]
        arguments = LHS.args[2:end]
        er = RHS

        if size(arguments, 1) == 1 && op in [-, :-, :transpose]
            # Unary operators - and transpose
            #  - ex = er  ->  ex = -er
            #  transpose(ex) = er  ->  ex = transpose(er)
            ex = LHS.args[2]
            eq = Expr(:(=), ex, Expr(:call, op, er))

            if ex != x
                eq, solved = solve(eq, x)
                return eq, solved
            else
                return eq, true
                dump(ex)
                dump(x)
                error("Should not happen")
            end

        elseif op in [*, :*, +, :+, -, :-, /, :/]   # ^, :^]
            e1 = arguments[1]
            if size(arguments, 1) == 2
                e2 = arguments[2]
            else
                # e2 = Expr(:call, op, arguments[2:end])  # Not correct
                e2 = Expr(:call) # , op, arguments[2])
                e2.args = [op; arguments[2:end]]  # Fix
            end

            if testIncidence(e1, x) && testIncidence(e2, x)
                println("Multiple occurances of unknown: ", prettyPrint(LHS))
                return eq, false
            elseif testIncidence(e1, x)
                if logSolve println("  left expression incidence") end
                ex = e1
                if op in [*, :*]
                    # ex * e2 * ... = er  ->  ex = er / (e2 * ...)
                    eq = Expr(:(=), ex, Expr(:call, :/, er, e2))
                elseif op in [+, :+]
                    # ex + e2 + ... = er  ->  ex = er - (e2 + ...)
                    eq = Expr(:(=), ex, sub(er, e2))
                elseif op in [-, :-]
                    # ex - e2 - ... = er  ->  ex = er + e2 + ...
                    if size(arguments, 1) > 2
                        e2.args[1] = +
                    end
                    eq = Expr(:(=), ex, add(er, e2))
                elseif op in [/, :/]
                    # ex / e1 / ... = er  ->  ex = er * e1 * ...
                    if size(arguments, 1) > 2
                        e2.args[1] = *
                    end
                    eq = Expr(:(=), ex, mult(er, e2))
                elseif op in [^, :^]
                    # ex ^ 2 = er  ->  ex = sqrt(er)  # Generalize, could be negative
                    # eq = Expr(:(=), ex, Expr(:call, :-, Expr(:call, :sqrt, er)))
                    return eq, false
                end
            else
                if logSolve println("  right expression incidence") end
                ex = e2
                if op in [*, :*]
                    # e1 * ex * ... = er  ->  ex * ... = e1 \ er
                    eq = Expr(:(=), ex, Expr(:call, :\, e1, er))
                elseif op in [+, :+]
                    # e1 + ex + ... = er  ->  ex + ... = er - e1
                    eq = Expr(:(=), ex, sub(er, e1))
                elseif op in [-, :-]
                    # e1 - ex - ... = er  ->  ex + ... = e1 - er
                    if size(arguments, 1) > 2
                        ex.args[1] = +
                    end
                    eq = Expr(:(=), ex, sub(e1, er))
                elseif op in [/, :/]
                    # e1 / ex / ... = er  ->  ex * ... = e1 / er
                    if size(arguments, 1) > 2
                        ex.args[1] = *
                    end
                    eq = Expr(:(=), ex, Expr(:call, /,  e1, er))
                end
            end

            if ex != x
                eq, solved = solve(eq, x)
                return eq, solved
            else
                if logSolve println("  Solved: ", prettyPrint(eq)) end
                return eq, true
            end
        else
            if logSolve
                println("Equations involving this operator are not solved: ", op)
                println("  Not solved: ", prettyPrint(eq))
            end
            # dump(eq)
            return eq, false
        end

        if logSolve println("  Never solved: ", prettyPrint(eq)) end
        error("Should not happen")
        return eq, false
    else
        # println("Should not happen:")
        # dump(eq.args[1])
        # dump(eq.args[2])
        # dump(x)
        return eq, false
        error("Should not happen")  # Check up!!!
    end

    if logSolve println("  Not solved: ", prettyPrint(eq)) end
    error("Should not happen")
    return eq, false
end

"""
    solveEquations(equations, variables, indices, assign, A, B)

Solve assigned equations.

* `equations`: infix string for original equations
* `variables`: infix string for original variables
* `indices`: indices for the equations to be printed
* `assign`: assign[j] contains the E-node to which V-node j is assigned or 0 if V-node j not assigned
* `A`: A[j] = if V[k] = der(V[j]) then k else 0
* `B`: B[i] = if E[l] = der(E[l]) then l else 0

"""
function solveEquations(equations, variables, indices, assign, A, B)
    (orgIndexVar, derOrderVar) = invertDer(A)
    (orgIndexEqu, derOrderEqu) = invertDer(B)
    (assignedVar,) = invertAssign(assign)

    for i in indices
        j = assignedVar[i]
        if j > 0
            if derOrderVar[j] == 1
                print("der(")
            elseif derOrderVar[j] > 1
                print("der", derOrderVar[j], "(")
            end

            print(variables[orgIndexVar[j]])

            if derOrderVar[j] > 0
                print(")")
            end
        end
        print(": ")

        if derOrderEqu[i] == 1
            print("DER() ")
        elseif derOrderEqu[i] > 1
            print("DER", derOrderEqu[i], "() ")
        end

        println(equations[orgIndexEqu[i]])
        println(solve(equations[i], variables[j]))
    end
end

"""
    solveSortedEquations(equations, variables, components, assign, A, B)

Solve sorted equations.

* `equations`: infix string for original equations
* `variables`: infix string for original variables
* `components`: cell array of components. Each component is a list of indices to E-nodes
* `assign`: assign[j] contains the E-node to which V-node j is assigned or 0 if V-node j not assigned
* `A`: A[j] = if V[k] = der(V[j]) then k else 0
* `B`: B[i] = if E[l] = der(E[l]) then l else 0

"""
function solveSortedEquations(equations, variables, components, assign, A, B)
    println("[assigned variable]: [differentiation] equation")
    println("Strongly connected components are enclosed in []")

    for c in components
        if length(c) > 1
            println("[")
        end

        solveEquations(equations, variables, c, assign, A, B)

        if length(c) > 1
            println("]")
        end
    end
end

# ----------------------------------------------------------------------------------------

# Derivatives of functions
divide_der_1(x, y) = 1 / y
divide_der_2(x, y) = x / y^2

power_der_1(x, y) = y * x^(y - 1)
power_der_2(x, y) =  x^y * log(x)

exp_der(x) = exp(x)
log_der(x) = 1 / x
sin_der(x) = cos(x)
cos_der(x) = -sin(x)
tan_der(x) = 1 / cos(x)^2
arcsin_der(x) = 1 / sqrt(1 - x^2)
arccos_der(x) = -1 / sqrt(1 - x^2)
arctan_der(x) = 1 / (1 + x^2)
sqrt_der(x) = 0.5 / sqrt(x)
sqrt_der_der(x) = -0.25 / (x * sqrt(x))

global timeInvariants = []
setTimeInvariants(tinv) = (global timeInvariants; timeInvariants = tinv)

global realStates = []
setRealStates(st) = (global realStates; realStates = st)

global dummyDerivatives

differentiate(e::Irrational) = zero

function differentiate(e)
    if logDifferentiate
        # print("\nDIFFERENTIATE: "); println(prettyPrint(e))
        # @show typeof(e)
    end

    if typeof(e) in [Float64, Int64, String, Bool]
        # diff = if ! derAsFunction; :($zero*$e) else zero*e end
        diff = if !derAsFunction; :($zero) else zero * e end
        if !noUnits
            diff = diff / SIUnits.Second
        end
    elseif typeof(e) <: Unitful.Quantity ||  typeof(e) <: Unitful.Unitlike
        diff = zero * e / Unitful.s
    elseif typeof(e) in [Array{Float64}, Array{Int64}]
        # diff = :($zero*$e) # zero*e
        diff = :($zero) # zero*e
    elseif e == :time
        diff = one
    elseif e in timeInvariants
        # println("TIME INVARIANT" )
        diff = zero
    elseif typeof(e) == Symbol
        diff = Expr(:call, :der, e)
    elseif typeof(e) == Instantiation.GetField
        # @show e.name
        if e.name in timeInvariants
            # println("TIME INVARIANT" )
            # diff = :($zero*$e)
            diff = :($zero)
        else
            if logDifferentiateVariable
                      print("DIFFERENTIATE VARIABLE:    "); println(e)
            end
            # println("\nGetField:")
            # @show e
            # dump(e)
            i = findfirst(isequal(e), realStates)

            if i != notFound
                diff = Der(e) # / Unitful.s
                # @show realStates i e diff
            else
                # diff = GetField(This(), Symbol("DER("*string(e.name)*")"))
                # dummyDerivatives[Symbol("DER("*string(e.name)*")")] = nothing
                diff = GetField(This(), Symbol("der_" * string(e.name))) # / Unitful.s # With units
                dummyDerivatives[Symbol("der_" * string(e.name))] = nothing

                if logDifferentiateVariable
                    println("  Dummy derivative: ", diff)
                end
            end

            # diff = Der(e)
            if logDifferentiateVariable
                print("  Differentiated variable: "); println(diff)
            end
        end

    elseif typeof(e) == Instantiation.Der
        if logDifferentiateVariable
                print("DIFFERENTIATE DERIVATIVE:    "); println(e)
        end

#        @show e e.base e.base.name realStates
        state = GetField(This(), Symbol("der_" * string(e.base.name)))
#        @show state
        i = findfirst(isequal(state), realStates)
#        @show i
        if i == notFound
            diff = GetField(This(), Symbol("der_der_" * string(e.base.name)))
#            @show diff
            if !noUnits
                diff = diff / Unitful.s
            end
            dummyDerivatives[Symbol("der_der_" * string(e.base.name))] = nothing
        else
            diff = Der(state)
            dummyDerivatives[Symbol("der_" * string(e.base.name))] = nothing
        end

        # @show e diff
        if logDifferentiateVariable
            print("  Differentiated derivative: "); println(diff)
        end

    elseif e.head in [:(=), :(:=)] # Equation
        diff = Expr(e.head, differentiate(e.args[1]), differentiate(e.args[2]))
        if logDifferentiate
            logModia("DIFFERENTIATE:    "); loglnModia(prettyPrint(e))
            logModia("  Differentiated: "); loglnModia(prettyPrint(diff)); loglnModia()
        end
    elseif e.head == :block && (isa(e.args[1], LineNumberNode) || isexpr(e.args[1], :line))
        diff = differentiate(e.args[2])
    elseif e.head in [:vect, Symbol("vect"), :row, Symbol("row"), :vcat, Symbol("vcat")]
        diff = Expr(e.head)
        for i in 1:length(e.args)
            push!(diff.args, differentiate(e.args[i]))
        end
        # @show e
        # @show diff
    elseif e.head in [:ref]
        diff = copy(e)
        diff.args[1] = differentiate(e.args[1])  # Probably not complete.
    elseif e.head in [:comprehension]
        diff = copy(e)
        diff.args[2].args[2] = differentiate(e.args[2].args[2])
    elseif e.head in [:if]
        diff = Expr(e.head, e.args[1], differentiate(e.args[2]), differentiate(e.args[3]))
    elseif e.head == :call
        op = e.args[1]
        arguments = e.args[2:end]
        if length(arguments) == 1 && op in [:der, :transpose, transpose, Symbol("'"), :residue, Modia.BasicStructuralTransform.residue, :skew, Modia.BasicStructuralTransform.skew] # skewCoords
            diffarg = differentiate(arguments[1])
            # @show e
            # @show diffarg
            if diffarg == zero
                diff = zero
            else
                diff = Expr(:call, op, diffarg)
            end
            # @show diff
        elseif (VERSION < v"0.7.0-DEV.2005" && op in [:zeros, zeros, :eye, eye] ) ||
              ( ! (VERSION < v"0.7.0-DEV.2005") && op in [:zeros, zeros, :eye, LinearAlgebra.I, :identity, identity] )
            diff = zero # !!!!
            # diff = Expr(:call, *, e, zero)  # = e*zero
        elseif op in [+, :+, -, :-]
            # der(e1 + e2 + e3) = der(e1) + der(e2) + der(e3)
            # der(e1 - e2 - e3) = der(e1) - der(e2) - der(e3)
            # diff = Expr(:call, op, [differentiate(arg) for arg in arguments])  # Not correct
            # diff = Expr(:call, op)
            # diff.args = [op; [differentiate(arg) for arg in arguments]]  # Fix
            diff = Expr(:call, op)
            firstNonZero = false

            for i in 1:length(arguments)
                term = differentiate(arguments[i])
                if term != zero
                    if i == 1
                        firstNonZero = true
                    end
                    push!(diff.args, term)
                end
            end

            if length(diff.args) == 1
                diff = zero
            elseif op in [+, :+] && length(diff.args) == 2
                # The cases:
                # der(x+5) -> der(x)
                diff = diff.args[2]
            elseif op in [-, :-] && length(arguments) > 1 && length(diff.args) == 2 && firstNonZero
                # The cases:
                # der(x-5) -> der(x)
                # der(5-x) -> -der(x)
                # der(-x) -> -der(x)
                diff = diff.args[2]
            end
        elseif op in [*, :*]
            # der(e1 * e2 * e3) = der(e1)*e2*e3 + e1*der(e2)*e3 + e1*e2*der(e3)
            arguments = e.args[2:end]
            diff = Expr(:call, +)
            for i in 1:length(arguments)
                terms = []
                for j in 1:length(arguments)
                    term = if i == j; differentiate(arguments[j]) else arguments[j] end
                    if term == zero
                        terms = []
                        break
                    else
                        push!(terms, term)
                    end
                end

                if terms != []
                    product = Expr(:call, *)
                    for t in terms
                        push!(product.args, t)
                    end

                    if length(diff.args) > 0
                        push!(diff.args, product)
                    end
                end
            end

            if length(diff.args) == 1 # Only + operator
                diff = zero
            elseif length(diff.args) == 2 # Only + zero, remove +
                diff = diff.args[2]
            end
        elseif op in [/, :/]
            # der(e1/e2) = der(e1)/e2 + e1/e2^2*der(e2)
            e1 = arguments[1]
            e2 = arguments[2]
            d1 = differentiate(e1)
            d2 = differentiate(e2)

            if d1 == zero
                diff = :($e1 / $e2^2 * $d2)
            elseif d2 == zero
                diff = :($d1 / $e2)
            else
                diff = :($d1 / $e2 + $e1 / $e2^2 * $d2)
            end
        elseif op in [^, :^] && typeof(arguments[2]) in [Int64]
            # der(e1 ^ n) = n*e1^(n-1)*der(e1)
            e1 = arguments[1]
            n = arguments[2]
            s = sub(n, 1)
            p = power(e1, s)
            d = differentiate(e1)
            diff = :($n * $p * $d)
        elseif op in [^, :^]
            # der(e1^e2) = e2*e1^(e2-1)*der(e1) + e1^e2*log(e1)*der(e2)
            e1 = arguments[1]
            e2 = arguments[2]
            d1 = differentiate(e1)
            d2 = differentiate(e2)
            if d1 == zero
                diff = :($e1^$e2 * log($e1) * $d2)
            elseif d2 == zero
                diff = :($e2 * $e1^($e2 - 1) * $d1)
            else
                diff = :($e2 * $e1^($e2 - 1) * $d1 + $e1^$e2 * log($e1) * $d2)
            end
        elseif op in [log, :log]
            # der(log(x)) = 1/x*der(x)
            e1 = arguments[1]
            d = differentiate(e1)
            diff = :(1 / $e1 * $d)
        elseif op in [exp, :exp]
            # der(exp(x)) = exp(x)*der(x)
            e1 = arguments[1]
            d = differentiate(e1)
            diff = :(exp($e1) * $d)
        elseif op in [cos, :cos]
            # der(cos(e1)) = -sin(e1)*der(e1)
            e1 = arguments[1]
            d = differentiate(e1)
            diff = :(-sin($e1) * $d)
        elseif op in [sin, :sin]
            # der(sin(e1)) = cos(e1)*der(e1)
            e1 = arguments[1]
            d = differentiate(e1)
            diff = :(cos($e1) * $d)
        elseif op in [tan, :tan]
            # der(tan(x)) = 1/cos(x)^2*der(x)
            e1 = arguments[1]
            d = differentiate(e1)
            diff = :(1 / cos($e1)^2 * $d)
        elseif op in [asin, :asin]
            # der(arcsin(x)) = 1/sqrt(1-x^2)*der(x)
            e1 = arguments[1]
            d = differentiate(e1)
            diff = :(1 / sqrt(1 - $e1^2) * $d)
        elseif op in [acos, :acos]
            # der(arccos(x)) = -1/sqrt(1-e1^2)*der(x)
            e1 = arguments[1]
            d = differentiate(e1)
            diff = :(-1 / sqrt(1 - $e1^2) * $d)
        elseif op in [atan, :atan]
            # der(arctan(x)) = 1/(1+x^2)*der(x)
            e1 = arguments[1]
            d = differentiate(e1)
            diff = :(1 / (1 + $e1^2) * $d)
        else
            # der(f(e1, e2, e3)) = f_der_1(e1, e2, e3)*der(e1) + f_der_2(e1, e2, e3)*der(e2) + f_der_3(e1, e2, e3)*der(e3)
            arguments = e.args[2:end]
            diff = Expr(:call, +)
            if op in [/, :/]
                op = :divide
            elseif op in [^, :^]
                op = :power
            end

            for i in 1:length(arguments)
                arg_der = differentiate(arguments[i])
                if arg_der != zero
                    path = split(string(op), ".")
                    if length(path) >= 2
                        mod = path[end-1]
                    else
                        mod = ""
                    end
                    func = path[end]
                    f_der = Symbol(string(func * "_der_", i))

                    if mod == ""
                        println("Derivative function ", string(f_der), " not found.")
                    elseif ! (f_der in names(getfield(Main, Symbol(mod))))
                        error("Derivative function ", string(f_der), " not found.")
                    else
                        f_der = getfield(getfield(Main, Symbol(mod)), f_der)
                    end
                    f_der_i = Expr(:call, f_der)
                    for a in arguments
                        push!(f_der_i.args, a)
                    end

                    term = mult(f_der_i, arg_der)
                    push!(diff.args, term)
                end
            end

            if length(diff.args) == 1
                push!(diff.args, zero)
            elseif length(diff.args) == 2
                diff = diff.args[2]
            end
        end
    else
        println("Error: ", prettyPrint(e))
        println("Not able to differentiate the above expression yet. ")
        dump(e)
        diff = nothing
    end
    # println("diff=", diff)
    return diff
end

function differentiateEquations(equations, variables, indices, assign, A, B, residualIndex)
    (orgIndexVar, derOrderVar) = invertDer(A)
    (orgIndexEqu, derOrderEqu) = invertDer(B)
    (assignedVar,) = invertAssign(assign)
    nSolved = 0
    component = []
    algebraic = []
    for i in indices
        if logSolve
            println("\nNEXT EQUATION")
        end

        j = assignedVar[i]
        assigned = j > 0
        if assigned
            if logSolved
                if derOrderVar[j] == 1
                    print("der(")
                elseif derOrderVar[j] > 1
                    print("der", derOrderVar[j], "(")
                end

                print(variables[orgIndexVar[j]])
                if derOrderVar[j] > 0
                    print(")")
                end
            end
        end

        if logSolved
            print(": ")

            if derOrderEqu[i] == 1
                print("DER() ")
            elseif derOrderEqu[i] > 1
                print("DER", derOrderEqu[i], "() ")
            end
            println(prettyPrint(equations[orgIndexEqu[i]]))
        end

        diff = equations[orgIndexEqu[i]]
        for d in 1:derOrderEqu[i]
            diff = differentiate(diff)
        end

        if logSolved
            println(prettyPrint(diff))
        end

        if assigned
            solveFor = variables[orgIndexVar[j]]
            solveForDummyDer = solveFor
            stateVar = GetField(This(), solveFor)
            # @show stateVar
            for d in 1:derOrderVar[j]
                # @show d
                if derAsFunction
                    solveFor = Expr(:call, :der, solveFor)
                    solveForDummyDer = Symbol("der_" * string(solveForDummyDer))
                else
                    # @show solveFor
                    # solveFor = Der(GetField(This(), solveFor))
                    # @show solveFor
                    # k = findfirst(realStates, stateVar)
                    # @show k
                    # if k == 0
                    if !(stateVar in realStates) ### Patch
                        # @show solveForDummyDer
                        # @show realStates
                        if d == 1
                            # solveForDummyDer = GetField(This(), Symbol("der_"*string(solveForDummyDer)))
                            solveForDummyDer = GetField(This(), Symbol("der_" * string(solveFor)))
                        else
                            solveForDummyDer = GetField(This(), Symbol("der_der_" * string(solveFor)))
                        end

                        if logDifferentiateVariable
                            println("DUMMY DERIVATIVE: ", solveForDummyDer)
                        end
                    else
                        solveForDummyDer = if !Modia.BasicStructuralTransform.useKinsol; nothing else Der(GetField(This(), solveFor)) end ### nothing # solveFor
                    end
                end
            end
            solveFor = solveForDummyDer
        else
            solveFor = nothing
        end

        if assigned && solveFor != nothing && length(indices) == 1
            (solution, solved) = solve(diff, solveFor)
        else
            solution = diff
            # solved = diff.head == :(:=)  # false  # To allow torn equations
            solved = false
        end

        if (assigned && solved && length(indices) == 1)  ||  !derAsFunction

    else
            push!(algebraic, solveFor)
            residualIndex += 1
            if noUnits && pushResiduals
                push!(component, (:(show = @show length(RESIDUALS)), true))
                push!(component, (:(show = @show $(diff.args[1]) - $(diff.args[2])), true))
                solution = :(dummy = push!(RESIDUALS, $(diff.args[1]) - $(diff.args[2])...))
            elseif noUnits
                solution = :(RESIDUALS[$residualIndex] = $(diff.args[1]) - $(diff.args[2]))
            else
                solution = :(RESIDUALS[$residualIndex] = float($(diff.args[1]) - $(diff.args[2])))
            end
            solved = false
        end

        push!(component, (solution, solved))
        if solved
            nSolved += 1
        end

        if logSolved
            if !solved
                print("NOT SOLVED: ")
                println()
                @show solveFor assigned solution
            end
            println(prettyPrint(solution))
        end
    end
    return component, algebraic, residualIndex
end

function differentiateSortedEquations(equations, variables, components, assign, A, B)
    if logSolved
        println("[assigned variable]: [differentiation] equation")
        println("Strongly connected components are enclosed in []")
    end

    global dummyDerivatives = OrderedDict()
    solvedComponents = []
    algebraic = []
    residualIndex = 0

    for c in components
        if logSolved && length(c) > 1
            println("[")
        end

        (solvedComponent, newAlgebraic, residualIndex) = differentiateEquations(equations, variables, c, assign, A, B, residualIndex)
        push!(solvedComponents, solvedComponent)

        for a in newAlgebraic
            push!(algebraic, a)
        end

        if logSolved && length(c) > 1
            println("]")
        end
    end

    # @show algebraic
    # @show solvedComponents
    # println("End of symbolic processing.")
    return solvedComponents, algebraic, residualIndex
end

function calcLinearCoeffcients(e)

end

function getLinearCoeffcients(e)
#=

=#
    linearCoefficients = Dict()
    vars = []
    BasicStructuralTransform.findIncidence!(vars, :(5x + 5 * (2y + x * 8 + z / 3 = 0)))
    linearCoefficients[:1] = 0

    for i in vars
        linearCoefficients[i] = 0
    end
    @show linearCoefficients
end

function newDifferentiateEquations(equations, variables, A, B, ESorted, ESolved)
    global dummyDerivatives = OrderedDict()
    println("newDifferentiateEquations")
    (orgIndexVar, derOrderVar) = invertDer(A)
    (orgIndexEqu, derOrderEqu) = invertDer(B)
    nSolved = 0
    component = []
    algebraic = []

    for k in 1:length(ESorted)
        i = ESorted[k]
        if logSolve
            println("\nNEXT EQUATION: $i")
        end
        j = ESolved[k]
        # @show ESolved[k]
        assigned = j > 0
        if assigned
            if logSolved
                if derOrderVar[j] == 1
                    print("der(")
                elseif derOrderVar[j] > 1
                    print("der", derOrderVar[j], "(")
                end

                print(variables[orgIndexVar[j]])

                if derOrderVar[j] > 0
                    print(")")
                end
            end
        end

        if logSolved
            print(": ")

            if derOrderEqu[i] == 1
                print("DER() ")
            elseif derOrderEqu[i] > 1
                print("DER", derOrderEqu[i], "() ")
            end
            println(prettyPrint(equations[orgIndexEqu[i]]))
        end

        diff = equations[orgIndexEqu[i]]

        for d in 1:derOrderEqu[i]
            diff = differentiate(diff)
        end

        if logSolved
            println(prettyPrint(diff))
        end

        # @show assigned
        if assigned
            solveForOrg = variables[orgIndexVar[j]]
            solveFor = variables[orgIndexVar[j]]

            for d in 1:derOrderVar[j]
                # @show d
                # @show solveFor
                # dump(solveFor)

                if d == 1
                    # solveForDummyDer = GetField(This(), Symbol("der_"*string(solveForDummyDer)))
                    if false ###
                        solveFor = GetField(This(), Symbol("der_" * string(solveFor)))
                    else
                        solveFor = Der(GetField(This(), solveFor))
                    end
                elseif d == 2
                    solveFor = GetField(This(), Symbol("der_der_" * string(solveForOrg)))
                else
                    solveFor = GetField(This(), Symbol("der_der_" * string(solveFor)))
                end
            end

        else
            solveFor = nothing
        end

        # @show solveFor
        if assigned && solveFor != nothing
            # @show solveFor
            (solution, solved) = solve(diff, solveFor)
            # @show solved
            if !solved
                # dump(solveFor)
                # dump(diff)
            end
        else
            solution = diff
            # solved = diff.head == :(:=)  # false  # To allow torn equations
            solved = false
        end

        if (assigned && solved) # ||  ! derAsFunction

        else
            push!(algebraic, solveFor)
            residualIndex = -j

            if noUnits && pushResiduals
                push!(component, (:(show = @show length(_r)), true))
                push!(component, (:(show = @show $(diff.args[1]) - $(diff.args[2])), true))
                solution = :(dummy = push!(_r, $(diff.args[1]) - $(diff.args[2])...))
            elseif noUnits
                solution = :(_r[$residualIndex] = $(diff.args[1]) - $(diff.args[2]))
            else
                solution = :(_r[$residualIndex] = float($(diff.args[1]) - $(diff.args[2])))
            end
            solved = false
        end

        push!(component, (solution, solved))

        if solved
            nSolved += 1
        end

        if logSolved
            if !solved
                print("NOT SOLVED: ")
            end
            println(prettyPrint(solution))
        end
    end

    return component
end

end


text:
struct BinnedCalibrator
    marginal_grid::Vector{Float64}
    marginal_h::Float64
    Q::Vector{Float64}
    Qo::Float64
end

function BinnedCalibrator(marginal_grid, Q, Qo)
    marginal_h = marginal_grid[2] - marginal_grid[1]
    BinnedCalibrator(marginal_grid, marginal_h, Q, Qo)
end

function BinnedCalibrator(marginal_grid, Q)
    BinnedCalibrator(marginal_grid, Q, zero(Float64))
end

function (c::BinnedCalibrator)(x)
    marginal_grid = c.marginal_grid
    marginal_h = c.marginal_h
    Q = c.Q
    # TODO: Some checks that we don't violate stuff..
    # Vectorized version of this..?
    idx = searchsortedlast(marginal_grid .- marginal_h/2, x)
    idx = max(1, min(idx, length(marginal_grid)))
    Q[idx] + c.Qo
end

mutable struct MinimaxCalibrator <: LinearEstimator
    Q::BinnedCalibrator
    max_bias::Float64
    sd::Float64 # not to be trusted though (i.e. recalculate)
    f1::Vector{Float64}
    f2::Vector{Float64}
    π1::Vector{Float64}
    π2::Vector{Float64}
    ds::MixingNormalConvolutionProblem
    f::BinnedMarginalDensity #have separate f for neighborhood and calibration?
    m::Int64
    target::LinearInferenceTarget
    ε_reg::Float64
    C::Float64
    δ::Float64
end

pretty_label(ma::MinimaxCalibrator) = "Minimax"

function (c::MinimaxCalibrator)(x)
    (c.Q)(x)
end

function MinimaxCalibrator(ds::MixingNormalConvolutionProblem,
                  f::BinnedMarginalDensity, m,
                  target = LFSRNumerator(2.0);
                  C=0.1,
                  max_iter=300,
                  ε = isinf(C) ? 0.001 : C,
                  tol=1e-5,
                  bias_check=false,
                  solver = GurobiSolver(OutputFlag=0),
                  max_smoother=false,
                  max_smoother_κ=0.1)

    n_priors = length(ds.priors)

    n_marginal_grid = length(ds.marginal_grid)
    h_marginal_grid = ds.marginal_grid[2] - ds.marginal_grid[1]

    ε_trunc = ε*h_marginal_grid
    f_marginal = f.marginal
    if max_smoother
        f_marginal_reg = smooth_poly_max.(f.marginal,
                                          ε_trunc,
                                          max_smoother_κ);
    else
        f_marginal_reg = max.(f.marginal,ε_trunc)
    end

    # Let us get linear functional
    #post_stats = posterior_stats(ds, target)
    #L = [z[1] for z in post_stats]
    L = posterior_stats(ds, target)

    A = ds.marginal_map;

    jm = Model(solver=solver)
    @variable(jm, π1[1:n_priors] >= 0)
    @variable(jm, π2[1:n_priors] >= 0)
    @variable(jm, t)
    @variable(jm, f1[1:n_marginal_grid])
    @variable(jm, f2[1:n_marginal_grid])

    @constraint(jm, sum(π1)== 1 )
    @constraint(jm, sum(π2)== 1 )
    @constraint(jm, A*π1 .== f1)
    @constraint(jm, A*π2 .== f2)

    if (C < Inf)
        #@constraint(jm, f1 .- f_marginal .<= C*f_marginal_reg)
        #@constraint(jm, f1 .- f_marginal .>= -C*f_marginal_reg)
        #@constraint(jm, f2 .- f_marginal .<= C*f_marginal_reg)
        #@constraint(jm, f2 .- f_marginal .>= -C*f_marginal_reg)
        @constraint(jm, f1 .- f_marginal .<= C*h_marginal_grid)
        @constraint(jm, f1 .- f_marginal .>= -C*h_marginal_grid)
        @constraint(jm, f2 .- f_marginal .<= C*h_marginal_grid)
        @constraint(jm, f2 .- f_marginal .>= -C*h_marginal_grid)
    end

    @objective(jm, Max, dot(L,π1-π2))
    status = solve(jm)

    t_curr = getobjectivevalue(jm)
    K_squared_curr = norm(sqrt(m)*(getvalue(f1).-getvalue(f2))./sqrt.(f_marginal_reg), 2)^2

    curr_obj = t_curr^2/(K_squared_curr+4)

    t_min = 0.0;
    t_max = t_curr;
    t_mid = (t_min + t_curr)/2;

    @constraint(jm, bias, dot(L,π1-π2) == t_mid)
    # update objective as well
    @objective(jm, Min, sum((f1-f2).^2 ./f_marginal_reg))
    status = solve(jm)

    for j=1:max_iter
        deriv = (2*t_mid*(4+getobjectivevalue(jm)*m)- getdual(bias)*t_mid^2*m)/(4+getobjectivevalue(jm)*m)^2

        if (deriv >= 0)
            t_min = t_mid
        else
            t_max = t_mid
        end
        t_mid = (t_min + t_max)/2

        JuMP.setRHS( bias, t_mid)
        status = solve(jm)

        curr_obj = t_mid^2/(4+getobjectivevalue(jm)*m)

        (t_max - t_min < tol) && break # actually could easily calibrate at beginning
    end

    fg1 = getvalue(f1)
    fg2 = getvalue(f2)

    K_squared = getobjectivevalue(jm)*m
    δ = sqrt(getobjectivevalue(jm))

    tg = t_mid

    π1 = getvalue(π1)
    π2 = getvalue(π2)

    L_g1 = dot(L, π1)
    L_g2 = dot(L, π2)
    L_g0 = (L_g1 + L_g2)/2

    Q = vec(m/(4+K_squared).*(fg1.-fg2)./f_marginal_reg .* tg)
    # ok the other bias term which I thought was going to 0...
    Qo = - 0.5*tg*sum( (fg1.^2 .- fg2.^2)./f_marginal_reg)
    Qo = Qo/(4 + K_squared)*m

    Qo = L_g0 + Qo

    max_bias = 2*tg/(4+K_squared)
    sd = tg*sqrt(K_squared)/(4+K_squared)
    #(max_bias, sd, max_bias^2+sd^2, next_obj)

    Q_c = BinnedCalibrator(ds.marginal_grid, Q, Qo)

    # Let us recalculate the max bias to be on the safe side?

    ma = MinimaxCalibrator(Q_c,
                      max_bias,
                      sd,
                      fg1,
                      fg2,
                      π1,
                      π2,
                      ds,
                      f,
                      m,
                      target,
                      ε,
                      C,
                      δ )

    if bias_check
        opt1 = check_bias(ma, solver=solver)
        opt2 = check_bias(ma; maximization= false, solver=solver)
        max_bias = max(abs(opt1), abs(opt2))
        ma.max_bias = max_bias
    end

    ma

end

function check_bias(Q::BinnedCalibrator,
                  ds::MixingNormalConvolutionProblem,
                  f::BinnedMarginalDensity, m,
                  target = LFSRNumerator(2.0);
                  C=Inf, maximization=true,
                  solver=GurobiSolver(OutputFlag=0))

    n_priors = length(ds.priors)

    f_marginal = f.marginal

    n_marginal_grid = length(ds.marginal_grid)

    h_marginal_grid = ds.marginal_grid[2] - ds.marginal_grid[1]

    A = ds.marginal_map;
    post_stats = posterior_stats(ds, target)
    L = [z[1] for z in post_stats]

    jm = Model(solver=solver)

    @variable(jm, π3[1:n_priors] >= 0)
    @variable(jm, f3[1:n_marginal_grid])

    @constraint(jm, sum(π3)== 1 )
    @constraint(jm, A*π3 .== f3)

    if (C < Inf)
        @constraint(jm, f3 .- f_marginal .<= C*h_marginal_grid)
        @constraint(jm, f3 .- f_marginal .>= -C*h_marginal_grid)
    end

    if maximization
        @objective(jm, Max, Q.Qo + dot(Q.Q,f3)-dot(L,π3))
    else
        @objective(jm, Min, Q.Qo + dot(Q.Q,f3)-dot(L,π3))
    end

    status = solve(jm)

    t_curr = getobjectivevalue(jm)
end

function check_bias(ma::MinimaxCalibrator; maximization=true)
    ds = ma.ds
    f = ma.f
    C = ma.C
    ε = ma.ε_reg
    m = ma.m
    target = ma.target

    Q = ma.Q

    check_bias(Q, ds, f, m, target; C=C, maximization=maximization)
end

ext: jl


text:
################################################################################
# Options
################################################################################

@with_kw mutable struct Options{T}
    # Options
	"Gauss-newton convergence criterion: tolerance."
	θ::T=1e-2

	"Initialization function of the primal dual vector."
	f_init::Function=rand

	"Initialization amplitude of the primal dual vector."
	amplitude_init::T=1e-8

	"Shift of the trajectory for the initial guess (useful for MPC, ususally, shift=1 in this case)."
	shift::Int=2^10

	# Regularization
	"Regularization of the residual and residual Jacobian."
	regularize::Bool=true

	"Current Jacobian regularization for each primal-dual variables."
	reg::Regularizer{T}=Regularizer()

	"Initial Jacobian regularization."
	reg_0::T=1e-3

	# Line search
	"Initial step size."
	α_0::T=1.0

	"Line search increase step."
	α_increase::T=1.2

	"Line search decrease step."
	α_decrease::T=0.5

	"Expected residual improvement."
	β::T=0.01

	"Number of line search iterations."
	ls_iter::Int=25

	"Minimum progress allowed before breaking out of inner loop."
	Δ_min::T=1e-9

	# Augmented Lagrangian Penalty
	"Initial augmented Lagrangian penalty."
	ρ_0::T=1.0

	"Fixed augmented Lagrangian penalty on the residual used for the line search trials."
	ρ_trial::T=1.0

	"Penalty increase step."
	ρ_increase::T=10.0

	"Maximum augmented Lagrangian penalty."
	ρ_max::T=1e7

	"Maximum Lagrange multiplier."
	λ_max::T=1e7

	"Control Dual ascent step size for all players."
	α_dual::T=1e0

	"State Dual ascent step size for each player."
	αx_dual::Vector{T}=ones(10)

	"Active set tolerance."
	active_set_tolerance::T=1e-4

	# Constraint satisfaction criterion
	"Dynamics constraint satisfaction criterion."
	ϵ_dyn::T=1e-3

	"State constraint satisfaction criterion."
	ϵ_sta::T=1e-3

	"Control constraint satisfaction criterion."
	ϵ_con::T=1e-3

	"Optimality constraint satisfaction criterion."
	ϵ_opt::T=1e-3

	# Augmented Lagrangian iterations.
	"Outer loop iterations."
	outer_iter::Int=7

	"Inner loop iterations."
	inner_iter::Int=20

	# Problem Scaling
	"Objective function scaling."
	γ::T=1e0

	# MPC
	"Number of time steps simulated with MPC"
	mpc_horizon::Int=20

	"Rate of the upsampling used for simulaion and feedback control in the MPC"
	upsampling::Int=2

	# Printing
	"Displaying the inner step results."
	inner_print::Bool=true

	"Displaying the outer step results."
	outer_print::Bool=true

	"Solver random seed."
	seed::Int=100

	"Reseting the duals."
	dual_reset::Bool=true
end

################################################################################
# Iterative Best Response Solver Options
################################################################################

@with_kw mutable struct IBROptions{T}
    # Options
	"Number of iterative best response iterations."
	ibr_iter::Int=100

	"Ordering of players, in the iterative best response scheme."
	ordering::Vector{Int}=Vector(1:100)

	"Minimum progress allowed before breaking out of the iterative best response loop."
	Δ_min::T=1e-9

	"Live plotting of the trajectories obtained through the iterative best response scheme."
	live_plotting::Bool=false
end


text:
# Solves subsonic airfoil
#   ARGS[1] = path to input file

# load PDESolver and solve the PDE
using PDESolver
using EulerEquationMod

const dirname = "./work"
# create temporary directory to do the solve in
if MPI.Comm_rank(MPI.COMM_WORLD) == 1
  if isdir(dirname)
    error("$dirname already exists, please delete it before proceeding")
  end
  mkdir(dirname)
end
MPI.Barrier(MPI.COMM_WORLD)
cd(dirname)

solvePDE(joinpath("../", ARGS[1]))


text:
export AbstractPreconditioner
export AbstractPreconditionerCache
export AbstractPreconditionerType
export PCinitialize!, PCapply!
export PCside, PCright, PCleft

"""
    AbstractPreconditioner

This is an abstract type representing an abstract preconditioner.

The available concrete implementations are:

- [`Solvent.Identity`](@ref)
- [`Solvent.Jacobi`](@ref)
"""
abstract type AbstractPreconditioner end

abstract type AbstractPreconditionerCache end
abstract type AbstractPreconditionerType end

abstract type PCside end
struct PCright <: PCside end
struct PCleft <: PCside end

"""
    function Preconditioner(
        pc_type::AbstractPreconditionerType,
        linearoperator!,
        Q,
    )

Construct preconditioner for `linearoperator!`.
"""
mutable struct Preconditioner{
    pcType <: AbstractPreconditionerType,
    OP,
    pccType <: AbstractPreconditionerCache,
    stype <: PCside,
} <: AbstractPreconditioner
    pc_type::pcType
    linop!::OP
    cache::pccType
    pc_side::stype

    function Preconditioner(
        pc_type::AbstractPreconditionerType,
        linearoperator!,
        Q,
    )
        pccache = cache(pc_type, linearoperator!, Q)
        OP = typeof(linearoperator!)
        pcType = typeof(pc_type)
        pccType = typeof(pccache)
        stype = typeof(pc_type.pc_side)

        return new{pcType, OP, pccType, stype}(
            pc_type,
            linearoperator!,
            pccache,
            pc_type.pc_side,
        )
    end
end

"""
    PCinitialize!(
        pc::AbstractPreconditioner,
        Q,
        Qrhs,
        args...,
    )

Perform precomputations for `pc`.
"""
function PCinitialize!(
    pc::AbstractPreconditioner,
    Q,
    Qrhs,
    args...,
)
    PCinitialize!(pc.pc_type, pc, Q, Qrhs, args...)
end

"""
    PCapply!(
        pc::AbstractPreconditioner,
        W,
        Q,
        args...,
    )

Apply preconditioner `pc` to `Q` and store the result in `W`.
"""
function PCapply!(
    pc::AbstractPreconditioner,
    W,
    Q,
    args...,
)
    PCapply!(pc.pc_type, pc, W, Q, args...)
end

include("identity.jl")
include("jacobi.jl")

max_stars_repo_name: CliMA/SystemSolvers

max_issues_repo_name: CliMA/SystemSolvers


text:
"""
  linearInterp()

linear interpolates the point x0 which lies in (x1,x2) onto the domain (s1,s2)
"""
function linearInterp(x0,x1,x2,s1,s2)
  A = [x1 1; x2 1]; rhs=[s1,s2]
  mbtemp = A\rhs
  m = mbtemp[1]
  b = mbtemp[2]

  s0 = m*x0 + b

  return s0
end

function distance(a,b)
    return sqrt((a[1] - b[1])^2 + (a[2] - b[2])^2)
end

function onSegment(p,q,r)
  val = abs(distance(p, q) + distance(q, r) - distance(p, r))
  if val < 1e-12
    return true
  else
    return false
  end
end

"""
  orientation(p,q,r) -> Int

Computes orientation of 3 (ordered) points p,q,r. where each is a 2-array
2 -> counterclockwise (CCW)
1 -> clockwise        (CW)
0 -> collinear

Algorithm is explained in https://www.geeksforgeeks.org/orientation-3-ordered-points/
"""
function orientation(p,q,r)
    val = (q[2] - p[2]) * (r[1] - q[1]) -
              (q[1] - p[1]) * (r[2] - q[2]);

    if (val == 0) return 0 end  # colinear

    return (val > 0) ? 1 : 2 # clock or counterclock wise
end

"""
  doIntersect(p1,q1,p2,q2) -> bool

Checks whether two lines, defined by (p1,q1) and (p2,q2) where pi,qi are 2-arrays.

Algorithm explained in https://www.geeksforgeeks.org/check-if-two-given-line-segments-intersect/.

Note we have not implemented the special case for overlapping line segments, as that will never be the case in finite element applications
"""
function doIntersect(p1,q1,p2,q2)
  # Find the four orientations needed for general case
  o1 = orientation(p1, q1, p2)
  o2 = orientation(p1, q1, q2)
  o3 = orientation(p2, q2, p1)
  o4 = orientation(p2, q2, q1)

  # General case
  if (o1 != o2 && o3 != o4)
    return true
  else
    return false
  end
end

"""
  isInside(polygon, n, p)

Returns true if the point p lies inside the polygon with n vertices
"""
function isInside(polygon, n, p)
  # There must be at least 3 vertices in polygon
  if (n < 3)  return false end

  # Create a point for line segment from p to infinite
  extreme = [100000.0, p[2]]

  # Count intersections of the above line with sides of polygon
  count = 0
  for i=1:n
    next = mod(i,n)+1
    # Check if the line segment from 'p' to 'extreme' intersects
    # with the line segment from 'polygon[i]' to 'polygon[next]'
    if (doIntersect(polygon[i], polygon[next], p, extreme))
      # If the point 'p' is colinear with line segment 'i-next',
      # then check if it lies on segment. If it lies, return true,
      # otherwise false
      if (orientation(polygon[i], p, polygon[next]) == 0)
         return onSegment(polygon[i], p, polygon[next])
      end
      count += 1
    end
  end

  # Return true if count is odd, false otherwise
  if isodd(count)
    return true
  else
    return false
  end
end

"""
  reverseQuadMap(quad,point) -> (s,t)

Computes the (s,t) computational domain point for an arbitrary point (x,y) in the physical element quad

For use with meshTransform
"""
function reverseQuadMap(quad,point)
  PRINT = false   # for testing

  x0 = point[1];   y0 = point[2]
  x1 = quad[1][1]; y1 = quad[1][2]
  x2 = quad[2][1]; y2 = quad[2][2]
  x3 = quad[3][1]; y3 = quad[3][2]
  x4 = quad[4][1]; y4 = quad[4][2]

  if PRINT
    println("(x0,y0) = (",x0,",",y0,")")
    println()
    println("(x1,y1) = (",x1,",",y1,")")
    println("(x2,y2) = (",x2,",",y2,")")
    println("(x3,y3) = (",x3,",",y3,")")
    println("(x4,y4) = (",x4,",",y4,")")
    println()
  end

  # define ci
  c1 =  x1 + x2 + x3 + x4 - 4*x0
  c5 =  y1 + y2 + y3 + y4 - 4*y0

  c2 = -x1 + x2 + x3 - x4
  c6 = -y1 + y2 + y3 - y4

  c3 = -x1 - x2 + x3 + x4
  c7 = -y1 - y2 + y3 + y4

  c4 =  x1 - x2 + x3 - x4
  c8 =  y1 - y2 + y3 - y4

  if PRINT
    println("c1 = ",c1)
    println("c2 = ",c2)
    println("c3 = ",c3)
    println("c4 = ",c4)
    println("c5 = ",c5)
    println("c6 = ",c6)
    println("c7 = ",c7)
    println("c8 = ",c8)
    println()
  end

  s = 0.0
  t = 0.0

  # case 1
  if abs(c4) < 1e-12 && abs(c8) < 1e-12
    #println("case 1")
    A = [c2 c3; c6 c7]; b = [-c1;-c5]

    x = A\b
    s = x[1]; t = x[2]

  # case 2A -- NOT CHECKED YET
  elseif abs(c4) < 1e-12
    #println("case 2A")
    if PRINT println("Case 2A!") end

    # c2 != 0
    if abs(c2) > 1e-12
      A = -c3*c8
      B = c2*c7 - c3*c6
      C = c2*c5 - c1*c6
      det = B^2 - 4*A*C

      println("A = ",A)
      println("B = ",B)
      println("C = ",C)

      if abs(A) < 1e-12
        t = -C/B

      else
        if PRINT println("det = ",det) end
        t1 = 0.5*(-B + sqrt(det))/A
        t2 = 0.5*(-B - sqrt(det))/A
        println("t1 = ",t1)
        println("t2 = ",t2)
        println()

        t = t1
      end

      s = -(c1+c3*t)/c2

    # c2==0
    else
      A = c3*c8
      B = c3*c6 - c2*c7 + c1*c8
      C = c1*c6 - c2*c5
      det = B^2 - 4*A*C

      println("A = ",A)
      println("B = ",B)
      println("C = ",C)

      if PRINT println("det = ",det) end

      t1 = 0.5*(-B + sqrt(det))/A
      t2 = 0.5*(-B - sqrt(det))/A
      println("t1 = ",t1)
      println("t2 = ",t2)
      println()

      s1 = -(c5+c7*t1)/(c6+c8*t1)
      s2 = -(c5+c7*t2)/(c6+c8*t2)

      println((c5+c7*t2))

      println("s1 = ",s1)
      println("s2 = ",s2)

      if (s1>1.0 || s1<-1.0)
        t = t2
        s = s2
      else
        t = t1
        s = s1
      end
    end
  # case 2B ---- NOT CHECKED YET
  elseif abs(c8) < 1e-12
    #println("case 2B")
    if PRINT println("Case 2B!") end
    # c6 != 0
    if abs(c6) > 1e-12
      A = -c7*c4
      B = c6*c3 - c7*c2
      C = c6*c1 - c5*c2
      det = B^2 - 4*A*C

      println("A = ",A)
      println("B = ",B)
      println("C = ",C)

      if abs(A) < 1e-12
        t = -C/B

      else
        if PRINT println("det = ",det) end
        t1 = 0.5*(-B + sqrt(det))/A
        t2 = 0.5*(-B - sqrt(det))/A
        println("t1 = ",t1)
        println("t2 = ",t2)
        println()

        t = t1
      end

      s = -(c1+c3*t)/c2

    # c6==0
    else
      A = c7*c4
      B = c7*c2 - c6*c3 + c5*c4
      C = c5*c2 - c6*c1
      det = B^2 - 4*A*C

      println("A = ",A)
      println("B = ",B)
      println("C = ",C)

      if PRINT println("det = ",det) end

      t1 = 0.5*(-B + sqrt(det))/A
      t2 = 0.5*(-B - sqrt(det))/A
      println("t1 = ",t1)
      println("t2 = ",t2)
      println()

      s1 = -(c1+c3*t1)/(c2+c4*t1)
      s2 = -(c1+c3*t2)/(c2+c4*t2)

      println((c1+c3*t2))

      println("s1 = ",s1)
      println("s2 = ",s2)

      if (s1>1.0 || s1<-1.0)
        t = t2
        s = s2
      else
        t = t1
        s = s1
      end
    end
  # case 3 ----- NOT CHECKED YET
  else
    #println("case 3")
    if PRINT println("Case 3!") end
    # define A,B,C
    A = c4*c7 - c3*c8                   # t^2
    B = c4*c5 - c3*c6 + c2*c7 - c1*c8   # t
    C = c2*c5 - c1*c6                   # 1

    # solve for solution t
    if abs(A)<1e-12
      t = -C/B
    else
      det = sqrt(B^2 - 4*A*C)
      t1 = 0.5*(-B + det)/A
      t2 = 0.5*(-B - det)/A

      if t1 > 1.0 || t1 < -1.0
        t = t2
      else
        t = t1
      end
    end

    # solve for solution s
    s = -(c1+c3*t)/(c2+c4*t)
  end

  # return s,t
  if PRINT println(s," ",t) end
  return s,t
end

"""
  isInsideDomain(xyA,cmA,point) -> bool

Transforms finite element solutions from domain A to point xyB, returns solution Uinterp at point xyB

Can handle 4 point (linear) and 9 point (quadratic) quadrilateral elements
"""
function isInsideDomain(xyA,cmA,point)
  nptsA   = size(xyA,1)
  nelA    = size(cmA,1)
  npe     = size(cmA,2); npe == 4 ? 4 : 8 # adjust for only 8 bdry pts
  order   = (npe==4 ? 1 : 2)
  Uinterp = 0.0
  foundInside = false
  foundVal = false

  # check over all nodes in Q1 for Q2 match
  #nodeTOL = 1e0
  #for ti=1:nptsA
  #    if sqrt((xyA[ti,1]-point[1])^2 + (xyA[ti,2]-point[2])^2) < nodeTOL
  #      println("found!")
  #      Uinterp = UA[ti]
  #      foundVal = true
  #      break
  #    end
  #end

  # loop over elements in A to find which one contains above point
  if foundVal==false
    for el=1:nelA
      if foundInside == false
        c = cmA[el,1:4]
        p1 = [xyA[c[1],1],xyA[c[1],2]]
        p2 = [xyA[c[2],1],xyA[c[2],2]]
        p3 = [xyA[c[3],1],xyA[c[3],2]]
        p4 = [xyA[c[4],1],xyA[c[4],2]]
        quad  = [p1,p2,p3,p4]

        # check whether point is in element
        if isInside(quad,4,point)
          return true
        # check whether point is on side of element
        else
          if onSegment(p1,point,p2)
            return true
          elseif onSegment(p2,point,p3)
            return true
          elseif onSegment(p3,point,p4)
            return true
          elseif onSegment(p4,point,p1)
            return true
          end
        end
      end
    end
  end

  return false
end

"""
  pointTransform(xyA,cmA,UA,point) -> Uinterp

Transforms finite element solutions from domain A to point xyB, returns solution Uinterp at point xyB

Can handle 4 point (linear) and 9 point (quadratic) quadrilateral elements
"""
function pointTransform(xyA,cmA,UA,point)
  nptsA   = size(xyA,1)
  nelA    = size(cmA,1)
  npe     = size(cmA,2); npe == 4 ? 4 : 8 # adjust for only 8 bdry pts
  order   = (npe==4 ? 1 : 2)
  Uinterp = 0.0
  foundInside = false
  foundVal = false

  #println("transformed point")

  # loop over elements in A to find which one contains above point
  if foundVal==false
    for el=1:nelA
      if foundInside == false
        c = cmA[el,1:4]
        p1 = [xyA[c[1],1],xyA[c[1],2]]
        p2 = [xyA[c[2],1],xyA[c[2],2]]
        p3 = [xyA[c[3],1],xyA[c[3],2]]
        p4 = [xyA[c[4],1],xyA[c[4],2]]
        quad  = [p1,p2,p3,p4]

        # check whether point is in element
        if isInside(quad,4,point)
          s,t = reverseQuadMap(quad,point)

          phi,_,_ = shape2D(s,t,1)

          x = shapeEval([p1[1],p2[1],p3[1],p4[1]],phi)
          y = shapeEval([p1[2],p2[2],p3[2],p4[2]],phi)

          maxerr = maximum([abs(x-point[1]),abs(y-point[2])])
          if maxerr > 1e-2
            #error("(s,t) error is too large, = ",maxerr)
          end
          #println("good transform")
          Uinterp = shapeEval(UA[c],phi)
          foundInside = true
        # check whether point is on side of element
        else
          if onSegment(p1,point,p2)
            s = 0.0; t=-1.0

            phi,_,_ = shape2D(s,t,1)

            x = shapeEval([p1[1],p2[1],p3[1],p4[1]],phi)
            y = shapeEval([p1[2],p2[2],p3[2],p4[2]],phi)

            maxerr = maximum([abs(x-point[1]),abs(y-point[2])])
            if maxerr > 1e-2
              #error("(s,t) error is too large, = ",maxerr)
            end
            #println("good transform")
            Uinterp = shapeEval(UA[c],phi)
            foundInside = true
          elseif onSegment(p2,point,p3)
            s = 1.0; t=0.0

            phi,_,_ = shape2D(s,t,1)

            x = shapeEval([p1[1],p2[1],p3[1],p4[1]],phi)
            y = shapeEval([p1[2],p2[2],p3[2],p4[2]],phi)

            maxerr = maximum([abs(x-point[1]),abs(y-point[2])])
            if maxerr > 1e-2
              #error("(s,t) error is too large, = ",maxerr)
            end
            #println("good transform")
            Uinterp = shapeEval(UA[c],phi)
            foundInside = true
          elseif onSegment(p3,point,p4)
            s = 0.0; t=1.0

            phi,_,_ = shape2D(s,t,1)

            x = shapeEval([p1[1],p2[1],p3[1],p4[1]],phi)
            y = shapeEval([p1[2],p2[2],p3[2],p4[2]],phi)

            maxerr = maximum([abs(x-point[1]),abs(y-point[2])])
            if maxerr > 1e-2
              #error("(s,t) error is too large, = ",maxerr)
            end
            #println("good transform")
            Uinterp = shapeEval(UA[c],phi)
            foundInside = true
          elseif onSegment(p4,point,p1)
            s = -1.0; t=0.0

            phi,_,_ = shape2D(s,t,1)

            x = shapeEval([p1[1],p2[1],p3[1],p4[1]],phi)
            y = shapeEval([p1[2],p2[2],p3[2],p4[2]],phi)

            maxerr = maximum([abs(x-point[1]),abs(y-point[2])])
            if maxerr > 1e-2
              #error("(s,t) error is too large, = ",maxerr)
            end
            #println("good transform")
            Uinterp = shapeEval(UA[c],phi)
            foundInside = true
          end
        end
      end
    end
  end

  if !(foundInside || foundVal)
    errorstr = string("point (",point[1],",",point[2],") not inside mesh")
    error(errorstr)
    return 0.0
  end

  return Uinterp
end
# point transform for mesh
function pointTransform(mesh::T,UA,point) where T<:AbstractMesh
  # transform meshA to xy and cm
  xyA = NodesToArray(mesh.xy)
  cmA = ElementsToArray(mesh.cm)

  Uinterp = pointTransform(xyA,cmA,UA,point)
  return Uinterp
end

"""
  meshTransform(xyA,cmA,UA,xyB) -> UB

Transforms finite element solutions from domain A to domain B, returns solutions on UB

Can handle 4 point (linear) and 9 point (quadratic) quadrilateral elements
"""
function meshTransform(xyA,cmA,UA,xyB)
  nptsB = size(xyB,1)
  UB    = zeros(nptsB)

  for i=1:nptsB
    point = [xyB[i,1],xyB[i,2]]
    UB[i] = pointTransform(xyA,cmA,UA,point)
  end

  return UB
end

# transform mesh using mesh objects and points
# xyB is in Nx2 array where xyB[:,1] = x-vals, xyB[:,2] = y-vals
function meshTransform(meshA::T,UA,xyB) where T<:AbstractMesh
  # transform meshA to xy and cm
  xyA = NodesToArray(meshA.xy)
  cmA = ElementsToArray(meshA.cm)

  # call meshTransform
  UB = meshTransform(xyA,cmA,UA,xyB)

  # return interpolated solutions
  return UB
end

# transform mesh using mesh objects
function meshTransform(meshA::T1,UA,meshB::T2) where
                  {T1<:AbstractMesh, T2<:AbstractMesh}
  # transform meshA to xy and cm
  xyA = NodesToArray(meshA.xy)
  cmA = ElementsToArray(meshA.cm)

  # transform meshB to xy
  xyB = NodesToArray(meshB.xy)

  # call meshTransform
  UB = meshTransform(xyA,cmA,UA,xyB)

  # return interpolated solutions
  return UB
end

"""
  car2pol(x,y) -> r,theta

Transforms cartesian coordinates into polar
"""
function car2pol(x,y)
  r = sqrt(x^2 + y^2)
  theta = atan2(y,x)

  return r,theta
end

"""
  pol2car(r,theta) -> x,y

Transforms polar coordinates into cartesian
"""
function pol2car(r,theta)
    x = r*cos(theta)
    y = r*sin(theta)

  return x,y
end


text:
using DistributedFactorGraphs
using Test

struct TestInferenceVariable1 <: InferenceVariable end
struct TestInferenceVariable2 <: InferenceVariable end

@testset "Needs-a-home tests" begin
    dfg = LightDFG{NoSolverParams}(solverParams=NoSolverParams())

    # Build a graph
    v1 = DFGVariable(:a, TestInferenceVariable1())
    v2 = DFGVariable(:b, TestInferenceVariable1())
    f1 = DFGFactor{TestFunctorInferenceType1}(:f1)
    addVariable!(dfg, v1)
    addVariable!(dfg, v2)
    addFactor!(dfg, [v1, v2], f1)

    # Standard tests
    newdfg = buildSubgraphFromLabels!(dfg, [:b])
    @test ls(newdfg) == [:b]
    # Okay it looks like this function only accepts variables, is that right?
    @test_throws Exception newdfg = buildSubgraphFromLabels!(dfg, [:b, :f1])
    newdfg = buildSubgraphFromLabels!(dfg, [:b, :a])
    @test symdiff(ls(newdfg), [:b, :a]) == []
    @test lsf(newdfg) == [:f1]

    # Check solvable filter
    @test setSolvable!(getVariable(dfg, :a), 1) == 1
    @test setSolvable!(getVariable(dfg, :b), 0) == 0
    @test setSolvable!(getFactor(dfg, :f1), 1) == 1
    newdfg = buildSubgraphFromLabels!(dfg, [:b, :a]; solvable = 1)
    @test ls(newdfg) == [:a]
end


text:
# pyplot_multiaxis.jl
#
# Demonstrate use of multiple overlap    ping axes
#
# gizmaa (https://gist.github.com/gizmaa/7214002)
# Julia 0.6.0
# 19.12.13
# Last Edit: 20.07.17

using PyPlot

##################
#  Make Up Data  #
##################
n = 15 # Number of values

x = [1:n;]
y1 = 10rand(n,1) - 2
y2 = 1000rand(n,1)
y3 = 30rand(n,1) - 15

##########
#  Plot  #
##########
fig = figure("pyplot_multiaxis",figsize=(10,10))
p = plot(x,y1,linestyle="-",marker="o",label="First") # Plot a basic line
ax = gca()
title("Multi-axis Plot")

xlabel("X Axis")
font1 = Dict("color"=>"blue")
ylabel("Left Axis",fontdict=font1)
setp(ax[:get_yticklabels](),color="blue") # Y Axis font formatting

################
#  Other Axes  #
################
new_position = [0.06;0.06;0.77;0.91] # Position Method 2
ax[:set_position](new_position) # Position Method 2: Change the size and position of the axis
#fig[:subplots_adjust](right=0.85) # Position Method 1

ax2 = ax[:twinx]() # Create another axis on top of the current axis
font2 = Dict("color"=>"purple")
ylabel("Right Axis",fontdict=font2)
p = plot_date(x,y2,color="purple",linestyle="-",marker="o",label="Second") # Plot a basic line
ax2[:set_position](new_position) # Position Method 2
setp(ax2[:get_yticklabels](),color="purple") # Y Axis font formatting

ax3 = ax[:twinx]() # Create another axis on top of the current axis
ax3[:spines]["right"][:set_position](("axes",1.12)) # Offset the y-axis label from the axis itself so it doesn't overlap the second axis
font3 = Dict("color"=>"green")
ylabel("Far Right Axis",fontdict=font3)
p = plot_date(x,y3,color="green",linestyle="-",marker="o",label="Third") # Plot a basic line
ax3[:set_position](new_position) # Position Method 2
setp(ax3[:get_yticklabels](),color="green") # Y Axis font formatting

axis("tight")

# Enable just the right part of the frame
ax3[:set_frame_on](true) # Make the entire frame visible
ax3[:patch][:set_visible](false) # Make the patch (background) invisible so it doesn't cover up the other axes' plots
ax3[:spines]["top"][:set_visible](false) # Hide the top edge of the axis
ax3[:spines]["bottom"][:set_visible](false) # Hide the bottom edge of the axis

fig[:canvas][:draw]() # Update the figure


text:
#! /bin/bash
#=
filedir=`dirname $( realpath $0 )`
exec /usr/bin/env julia \
    --project=$filedir \
    -J ${filedir}/build/sys_toml.so \
    -e "include(popfirst!(ARGS))" \
    ${BASH_SOURCE[0]} $@
=#

using Base.Iterators

""" Gets the real path of a file that may or may not be a symlink """
resolving_links(path::String) = if islink(path) readlink(path) else path end

# If the TOML sysimage isn't being used, include TOML normally.
if !isdefined(Main, :TOML)
    using Pkg; Pkg.instantiate()
    using TOML
end

const SYSTEM76_KBD_BACKLIGHT_DIR = "/sys/class/leds/system76::kbd_backlight/"
const KBD_COLOR_FILES = map(
    f -> joinpath(SYSTEM76_KBD_BACKLIGHT_DIR::String, f),
    ["color_left", "color_center", "color_right"]
)
const COLOR_MACROS_FILE = joinpath(dirname(resolving_links(@__FILE__)), "ColorMacros.toml")

openwrite(file, contents) = open(f -> write(f, contents), file, "w")

expand_colormacro(config::Dict{AbstractString,Any}, m::String) =
    get(config["colors"], m, m)

expand_patternmacro(config::Dict{AbstractString,Any}, m::String) =
    map(
        p -> expand_colormacro(config, p),
        get(config["patterns"], m, [m])
    )

function main()
    # Parse the color macros from a file
    color_macros = TOML.parsefile(COLOR_MACROS_FILE::String)

    # Expand the color macros into corresponding hex codes
    expanded_args = map(
        t -> expand_patternmacro(color_macros, t),
        ARGS
    ) |> flatten

    # Take the first three expanded args
    # This cycles the args and takes 3, so given N args, the colors will be:
    # N  | L C R
    # -----------
    # 1  | 1 1 1
    # 2  | 1 2 1
    # 3+ | 1 2 3
    colors = take(expanded_args |> cycle, 3)

    # Finally, write the colors into their respective files
    foreach(
        kv -> openwrite(kv...),
        zip(KBD_COLOR_FILES::Vector{String}, colors)
    )
end

main()

max_stars_repo_name: periodicaidan/oryx-kbd-colors

max_issues_repo_name: periodicaidan/oryx-kbd-colors

max_forks_repo_name: periodicaidan/oryx-kbd-colors


text:
module CUDAint

using ..CUDA

import ..Flux: Flux
import Zygote
using Zygote: @adjoint
import NNlib, NNlibCUDA

include("cudnn.jl")

end


text:
using Oceananigans.MultiRegion
using Oceananigans.MultiRegion: reconstruct_global_grid, reconstruct_global_field, getname
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom, GridFittedBoundary

devices(::CPU, num) = nothing
devices(::GPU, num) = Tuple(0 for i in 1:num)

@testset "Testing multi region grids" begin
    for arch in archs

        region_num   = [2, 4, 5]
        partitioning = [XPartition]

        grids = [LatitudeLongitudeGrid(arch, size=(20, 20, 1), latitude=(-80, 80), longitude=collect(range(-180, 180, length=21)), z=(0, 1)),
                RectilinearGrid(arch, size=(20, 20, 1), x=(0, 1), y=collect(range(0, 1, length=21)), z=(0, 1))]

        immersed_boundaries = [GridFittedBottom((x, y)->0.5),
                               GridFittedBottom(arch_array(arch, [0.5 for i in 1:20, j in 1:20])),
                               GridFittedBoundary((x, y, z)->z>0.5),
                               GridFittedBoundary(arch_array(arch, [false for i in 1:20, j in 1:20, k in 1:1]))]

        for grid in grids, P in partitioning, regions in region_num
            @info "Testing multi region $(getname(grid)) on $regions $(P)s"
            mrg = MultiRegionGrid(grid, partition = P(regions), devices = devices(arch, regions))

            @test reconstruct_global_grid(mrg) == grid

            for FieldType in [CenterField, XFaceField, YFaceField]
                @info "Testing multi region $(FieldType) on $(getname(grid)) on $regions $(P)s"

                par_field = FieldType(mrg)
                ser_field = FieldType(grid)

                set!(ser_field, (x, y, z) -> x)
                @apply_regionally set!(par_field, (x, y, z) -> x)

                fill_halo_regions!(ser_field)
                fill_halo_regions!(par_field)

                rec_field = reconstruct_global_field(par_field)

                @test all(Array(rec_field.data.parent) .≈ Array(ser_field.data.parent))
            end

            for immersed_boundary in immersed_boundaries
                @info "Testing multi region immersed boundaries on $(getname(grid)) on $regions $(P)s"
                ibg = ImmersedBoundaryGrid(grid, immersed_boundary)
                mrg = MultiRegionGrid(ibg, partition = P(regions), devices = devices(arch, regions))

                @test on_architecture(arch, reconstruct_global_grid(mrg)) == ibg
            end
        end
    end
end


text:
#Generic abstract sampler for Parallel Tempering
abstract type PTSampler <: AbstractMCMC.AbstractSampler end
AS = AbstractMCMC.AbstractSampler

mutable struct DEO{T,S<:AS, } <: PTSampler
    β::Vector{T}
    κ_exp::S
    nexp::Int
end

mutable struct FullSweep{T, S<:AS} <: PTSampler
    β::Vector{T}
    κ_exp::S
    nexp::Int
end


text:
abstract RegressionModel

function Model(X::AbstractMatrix, y::AbstractVector, regression_type::Symbol, kernel::Symbol)
    if kernel == :linear
        Model(X, y, regression_type)
    else
        error("kernel=$(kernel) is not implemented yet")
    end
end

function Model(X::AbstractMatrix, y::AbstractVector, regression_type::Symbol, f::RegressionFunction=LinearRegressionFunction())
    if regression_type==:ordinal
        w0 = vec(mean(X,1))
        OrdinalModel(f, w0)
    elseif regression_type==:binomial
        w0 = vec(mean(X[y.==1,:],1).-mean(X[y.==-1,:]))
        BinomialModel(f, w0)
    elseif regression_type==:multinomial
        Y = unique(y)
        k = length(Y)
        n, m = size(X)
        w0 = zeros(m,k-1)
        w0_k = mean(X[y.==Y[k],:],1)
        for i in 1:k-1
            w0[:, i] = mean(X[y.==Y[i],:],1)-w0_k
        end
        MultinomialModel(f, vec(w0), k)
    else
        error("regression_type=$(regression_type) is not implemented yet")
    end
end

values(model::RegressionModel, X::AbstractMatrix, theta::AbstractVector=model.theta) = values(model.f, X, theta)

# Ordinal regression

type OrdinalModel <: RegressionModel
    f::RegressionFunction
    theta::Vector            # parameter vector of regression function
end

predict(model::OrdinalModel, X::AbstractMatrix) = values(model.f, X)

# binomial regression

type BinomialModel <: RegressionModel
    f::RegressionFunction
    theta::Vector            # parameter of regression function
end

predict(model::BinomialModel, X::AbstractMatrix) = sign(values(model, X))

# multinomial regression

type MultinomialModel <: RegressionModel
    f::RegressionFunction
    theta::Vector    # k-1 stacked parameter vectors of regression function
    k::Int           # number of classes
end
theta(model::MultinomialModel) = reshape(model.theta, length(model.theta) / (model.k-1), model, model.k-1)

function values(model::MultinomialModel, X::AbstractMatrix, theta::AbstractVector=model.theta)
    n = size(X, 1)
    v = zeros(n, model.k-1)
    m = length(theta)/(model.k-1)
    for k = 1:model.k-1
        v[:,k] = values(model.f, X, theta[(k-1)*m+1:k*m])
    end
    return v
end
predict(model::MultinomialModel, X::AbstractMatrix) = (v = values(model, X); n=size(X,1); [indmax([v[i,:] 0]) for i in 1:n])


text:
# Custom structure for holding data for sparse, variational Gaussian process
struct SVGP_data
    x::Matrix{Float64}   # x values for observed data
    y::Vector{Float64}   # y values for observed data
    n::Int64              # Number of observations
    p::Int64              # Number of input dimensions
end

# Constructor for building by feeding X and Y data.
#    This should be the most common constructor
SVGP_data(x::Matrix{Float64}, y::Vector{Float64}) = SVGP_data(x, y, size(x)[1], size(x)[2])

# Constructor for dummy data object
#    Required for building SVGP object with full data not accessible
SVGP_data(dummy::String) = SVGP_data(zeros(1,2), [0], 1, 2)

# Custom structure for holding parameters for sparse, variational Gaussian process
mutable struct SVGP_params
    const_mean::Vector{Float64}                           # constant mean term
    log_rho::Matrix{Float64}                              # log correlation length for GP
    log_kappa::Vector{Float64}                            # log GP variance
    log_sigma::Vector{Float64}                            # log GP error standard deviation
    inducing_L::LowerTriangular{Float64,Matrix{Float64}}  # Cholesky of variational covariance
    inducing_mean::Vector{Float64}                        # Inducing point mean parameters
    inducing_locs::Matrix{Float64}                        # Inducing point location parameters
    function SVGP_params(inp_data::SVGP_data, n_inducing::Int64)
        ni   = n_inducing
        np   = size(inp_data.x)[2]
        # Initialize the correlation lengths to 1/4 of the domain in each coordinate dimension
        #     This requires the data to exist, hence the dummy data structure for now
        # TODO: More elegent way to do this. Maybe just an argument but somewhat want something better
        lrho = reshape([log( (maximum(inp_data.x[:,ii]) - minimum(inp_data.x[:,ii]))/4. ) for ii in 1:np], (1,np))
        # Initialize the variance based on the variance of the data
        #     This requires the data to exist, hence the dummy data structure for now
        # TODO: More elegent way to do this.
        lkap = [log(var(inp_data.y))]
        lssq = [-4.0]
        icov = 0.1 * I + zeros(ni, ni);

        # Initialize inducing point locations to be a Latin hypercube across the domain
        x_range = [(minimum(inp_data.x[:,ii]), maximum(inp_data.x[:,ii])) for ii in 1:np]
        xi      = scaleLHC(LHCoptim(ni, np, 1000)[1], x_range);

        # Initialize the inducing point means using the nearest neighbor observations
        kdtree      = KDTree(transpose(inp_data.x))
        idxs, dists = knn(kdtree, transpose(xi), 1);
        init_mean   = [inp_data.y[ii[1]] for ii in idxs];

        c_mean = [mean(inp_data.y)]

        return new(c_mean, lrho, lkap, lssq, cholesky(icov).L, init_mean, xi)
    end
end

# Custom structure for holding sparse, variational Gaussian process object
struct Inference_obj
    n_inducing::Int64      # Number of inducing points
    n_functions::Int64     # Number of latent functions to model
    data::SVGP_data
    params::Vector{SVGP_params}
    likelihood::String     # Label for which likelihood to use. "gaussian" or "poisson" for now
    ll_function
end

# Constructor for custom likelihood
function Inference_obj(x::Matrix{Float64}, y::Vector{Float64}, ni::Int64, likelihood::String, ll_func, nf::Int64)
    Inference_obj(ni, nf, SVGP_data(x,y), [SVGP_params(SVGP_data(x,y), ni) for ii in 1:nf], likelihood, ll_func)
end

# Constructor from data for general likelihoods
#     Should often be the
function Inference_obj(x::Matrix{Float64}, y::Vector{Float64}, ni::Int64, likelihood::String)
    if likelihood == "gaussian"
        marg_like = gaussian_likelihood
        n_func    = 1
        param_arr = [SVGP_params(SVGP_data(x,y), ni)]
    elseif likelihood == "poisson"
        marg_like = poisson_likelihood
        n_func    = 1
        param_arr = [SVGP_params(SVGP_data(x,y), ni)]
        param_arr[1].inducing_mean = log.(param_arr[1].inducing_mean .+ 1.0);
        param_arr[1].const_mean    = log.(param_arr[1].const_mean);
    elseif likelihood == "gev"
        marg_like = gev_likelihood
        n_func    = 3
        param_arr = [SVGP_params(SVGP_data(x,y), ni) for ii in 1:3]
    elseif likelihood == "gumbel"
        marg_like = gumbel_likelihood
        n_func    = 2
        param_arr = [SVGP_params(SVGP_data(x,y), ni) for ii in 1:n_func]
    else
        error(string(likelihood," not implemented as a likelihood. It can be implemented by creating the custom log likelihood"))
    end

    Inference_obj(ni, n_func, SVGP_data(x,y), param_arr, likelihood, marg_like)
end

# Constructor from data defaulting to Gaussian likelihood
function Inference_obj(x::Matrix{Float64}, y::Vector{Float64}, ni::Int64)
    return Inference_obj(x, y, ni, "gaussian")
end

# Custom structure for holding sparse, variational Gaussian process object
struct SVGP_obj
    n_inducing::Int64      # Number of inducing points
    data::SVGP_data
    params::SVGP_params
    distrib::String        # Label for which likelihood to use. "gaussian" or "poisson" for now
    function SVGP_obj(n_inducing::Int64, data::SVGP_data, params::SVGP_params, distrib::String)
        return new(n_inducing, data, params, distrib)
    end
end

# Constructor from data defaulting to Gaussian likelihood
function SVGP_obj(x::Matrix{Float64}, y::Vector{Float64}, ni::Int64)
    return SVGP_obj(ni, SVGP_data(x,y), SVGP_params(SVGP_data(x,y), ni), "gaussian")
end

# Constructor from data for general likelihoods
#     Should often be the
function SVGP_obj(x::Matrix{Float64}, y::Vector{Float64}, ni::Int64, distrib::String)
    return SVGP_obj(ni, SVGP_data(x,y), SVGP_params(SVGP_data(x,y), ni), distrib)
end

# Hierarchical SVGP structure. Holds global and local SVGPs
#     Not for distributed calculations as all data is held within the object
struct HSVGP_obj
    ni_global::Int64
    global_obj::SVGP_obj
    local_svgps::Vector{SVGP_obj}
    n_parts::Int64
    function HSVGP_obj(
        x::Matrix{Float64},
        y::Vector{Float64},
        ni_local::Int64,
        ni_global::Int64,
        part_labels::Vector{Int64},
        distrib::String
    )
        n_parts     = maximum(part_labels)
        local_svgps = [SVGP_obj(x[part_labels .== ii, :], y[part_labels .== ii], ni_local, distrib) for ii in 1:n_parts]
        data        = SVGP_data(x,y)
        obj         = SVGP_obj(ni_global, SVGP_data("global"), SVGP_params(data, ni_global), "gaussian")
        return new(ni_global, obj, local_svgps, n_parts)
    end
end

# Constructor from data defaulting to Gaussian likelihood
function HSVGP_obj(x::Matrix{Float64}, y::Vector{Float64}, ni_local::Int64, ni_global::Int64, part_labels::Vector{Int64})
    return HSVGP_obj(x, y, ni_local, ni_global, part_labels, "gaussian")
end

# UNUSED SO FAR
# Possible struct for carrying optimization info to container inside the full structs
# mutable struct SVGP_opt
#     trace_elbo::Vector{Float64} # Optimization trace for fit
#     function SVGP_opt(x)
#         return new(x)
#     end
# end

max_stars_repo_name: ayanbiswas/PRISM

max_issues_repo_name: ayanbiswas/PRISM

max_forks_repo_name: ayanbiswas/PRISM


text:
module TopOpt88

include("FastConv.jl/src/FastConv.jl")

using LinearAlgebra, SparseArrays
using Plots, PlotUtils, Statistics
import .FastConv.fastconv

BLAS.set_num_threads(1)

abstract type Filter end
abstract type TopOptProblem end
abstract type TopOptResult end

include("utils.jl")

struct FilterBase <: Filter; end
struct FilterConv <: Filter; end

export TopOptProblemBase, TopOptProblemConv
export optimize
export line_load!, fix_dof!, symmetry_axis!
export use_elements!

# Define operator for Kronecker product
⊗(x, y) = kron(x, y)

"""
    fastconv_2(E::Array{T,N}, k::Array{T,N}) where {T,N}

Convolution that returns an array of the same size as the input.

# Arguments
- `E::Array{T`: TODO
- `N}`: TODO
- `k::Array{T`: TODO
- `N}`: TODO

"""
function fastconv_2(E::Array{T,N}, k::Array{T,N}) where {T,N}
    i, j = size(E)
    ki, kj = size(k)

    ret = fastconv(E, k)[ki ÷ 2 + 1:i + ki ÷ 2, kj ÷ 2 + 1:j + kj ÷ 2]

    return ret
end

struct FEProblem
    F::Array{Float64}
    U::Array{Float64}
    Kₑ::Array{Float64}
    iK::Array{Int}
    jK::Array{Int}
    edofMat::Array{Int}
    freedofs::Array{Integer}
    nx::Int
    ny::Int
end

mutable struct TopOptProblemBase <: TopOptProblem
    nx::Int
    ny::Int
    rmin::Float64
    volfrac::Float64
    p::Float64
    dofs::Int
    F::Array{Float64}
    passive::Array{Int,2}
    fixed_dofs::Array{Int}

    function TopOptProblemBase(nx, ny, rmin, volfrac, p)
        dofs::Int = 2 * (nx + 1) * (ny + 1)
        F::Array{Float64} = zeros(dofs)
        passive::Array{Int,2} = zeros(ny, nx)
        fixed_dofs::Array{Int} = Array([])
        new(nx, ny, rmin, volfrac, p, dofs, F, passive, fixed_dofs)
    end
end

mutable struct TopOptProblemConv <: TopOptProblem
    nx::Int
    ny::Int
    rmin::Float64
    volfrac::Float64
    p::Float64
    dofs::Int
    F::Array{Float64}
    passive::Array{Int,2}
    fixed_dofs::Array{Int}

    function TopOptProblemConv(nx, ny, rmin, volfrac, p)
        dofs::Int = 2 * (nx + 1) * (ny + 1)
        F::Array{Float64} = zeros(dofs)
        passive::Array{Int,2} = zeros(ny, nx)
        fixed_dofs::Array{Int} = Array([])
        new(nx, ny, rmin, volfrac, p, dofs, F, passive, fixed_dofs)
    end
end

mutable struct TopOptResultBase{T} <: TopOptResult
    prob::TopOptProblemBase
    iter::Int
    obj::Float64
    vol::Float64
    x::T
    x̃::T
    H::SparseMatrixCSC
    Hs::Array

    function TopOptResultBase(prob, x)
        nx = prob.nx
        ny = prob.ny
        iter = 0
        obj = maxintfloat(Float64)
        vol = 1.0
        x̃ = copy(x)
        H = sparse([], [], [], nx*ny, nx*ny)
        Hs = Array{Float64}([])
        new{typeof(x)}(prob, iter, obj, vol, x, x̃, H, Hs)
    end
end

mutable struct TopOptResultConv{T} <: TopOptResult
    prob::TopOptProblemConv
    iter::Int
    obj::Float64
    vol::Float64
    x::T
    x̃::T
    x̂::T
    Hs::Array
    h::Array

    function TopOptResultConv(prob, x)
        nx = prob.nx
        ny = prob.ny
        iter = 0
        obj = maxintfloat(Float64)
        vol = 1.0
        x̃ = copy(x)
        x̂ = copy(x)
        Hs = Array{Float64}([])
        kern_size = length(-ceil(prob.rmin) + 1:ceil(prob.rmin) - 1)
        h = zeros(kern_size, kern_size)
        new{typeof(x)}(prob, iter, obj, vol, x, x̃, x̂, Hs, h)
    end
end

TopOptResult(prob::TopOptProblemBase, x) = TopOptResultBase(prob, x)
TopOptResult(prob::TopOptProblemConv, x) = TopOptResultConv(prob, x)

function optimize(prob::TopOptProblem; Δ=0.01, filter=1)
    # Material properties
    E₀ = 1.0
    E_min = 1e-9
    ν = 0.3

    nx::Int = prob.nx
    ny::Int = prob.ny
    rmin = prob.rmin
    #
    nodenrs::Array{Int} = reshape(1:(1 + nx) * (1 + ny), 1 + ny, 1 + nx);
    edofVec::Array{Int} = reshape(2 .* nodenrs[1:end - 1,1:end - 1] .+ 1, nx * ny, 1);

    # Array with indices relative to the first dof
    rel_ind = [0 1 (2 * ny .+ [2 3 0 1]) -2 -1]
    edofMat = repeat(edofVec, 1, 8) .+ repeat(rel_ind, nx * ny);

    iK::Array{Int} = reshape((edofMat ⊗ ones(8, 1))', 64 * nx * ny);
    jK::Array{Int} = reshape((edofMat ⊗ ones(1, 8))', 64 * nx * ny);

    dim = prob.dofs
    # Define loads and supports (half MBB-beam)
    alldofs = Array(1:dim)
    freedofs = setdiff(alldofs, prob.fixed_dofs)

    # Get element stiffness matrix
    Kₑ = element_K(E₀, ν)
    F = prob.F
    U = zeros(dim)
    # Define FE-problem
    fe_prob = FEProblem(prob.F, U, Kₑ, iK, jK, edofMat, freedofs, nx, ny)
    # Intatiate result object
    x = prob.volfrac .* ones(ny, nx)
    res = TopOptResult(prob, x)
    # Prepare filter
    set_filter!(res, prob, rmin)

    # Initial value for β
    β = 1.0

    # Initilaize sensitivities
    ∂c = similar(x)
    ∂V = zeros(ny, nx)

    loop = 0
    loopbeta = 0
    change = 1.0
    p = prob.p

    # Start iteration
    while change > Δ
        loop += 1
        loopbeta += 1
        # FE-analysis
        solve_fe!(fe_prob, res.x̃, p, E₀, E_min)
        # Objective function and sensitivity analysis
        cₑ = compliance(fe_prob, E₀, E_min)
        c = sum((E_min .+ res.x̃.^p .* (E₀ - E_min)) .* cₑ)
        # Update ∂c and ∂V
        sensitivities!(∂c, ∂V, fe_prob, res, β, p, E₀, E_min, cₑ, filter)
        # Optimality criteria update of design variables and physical densities
        # and filtering/modification of sensitivities
        x_new = optimal_crit(prob, res, β, ∂c, ∂V, filter)
        # Compute the hange
        change = maximum(abs.(x_new - res.x))
        # Update result
        res.x .= x_new
        # Update β parameter
        if filter == 3 && β < 512 && ((loopbeta ≥ 50) || change ≤ 0.01)
            β = 2 * β
            loopbeta = 0
            change = 1
            println("Parameter β increased to $β")
        end
        # Print results
        vol_x = mean(res.x̃)
        println("It. $loop  Obj.: $c  Vol.: $vol_x  ch.: $change")
        # Plot the result
        cmap = cgrad(:greys, rev=true)
        plot_1 = heatmap(res.x̃, c=cmap, aspect_ratio=:equal, yflip=true,
                        grid=false, axis=false, colorbar=false)
        display(plot_1)

    end
    # Mirror array to show the full structure
    x_full = [res.x̃[:, end:-1:1] res.x̃]
    cmap = cgrad(:greys, rev=true)
    plot_1 = heatmap(x_full, c=cmap, aspect_ratio=:equal, yflip=true,
                    grid=false, axis=false, colorbar=false)
    savefig(plot_1, "plot_88.pdf")
    display(plot_1)
    return res
end

function set_filter!(res::TopOptResult, prob::TopOptProblemBase, rmin)
    nx = prob.nx
    ny = prob.ny
    dim_1::Int = nx * ny * (2 * (ceil(rmin) - 1) + 1)^2
    iH::Array{Int} = ones(dim_1)
    jH::Array{Int} = ones(size(iH))
    sH = zeros(size(iH))
    k = 0
    for i1 in 1:nx, j1 in 1:ny
        e1 = (i1 - 1) * ny + j1
        r_1 = max(i1 - (ceil(rmin) - 1), 1)
        r_2 = min(i1 + (ceil(rmin) - 1), nx)
        s_1 = max(j1 - (ceil(rmin) - 1), 1)
        s_2 = min(j1 + (ceil(rmin) - 1), ny)
        for i2 in r_1:r_2, j2 in s_1:s_2
            e2 = (i2 - 1) * ny + j2
            k += 1
            iH[k] = e1
            jH[k] = e2
            sH[k] = max(0, rmin - sqrt((i1 - i2)^2 + (j1 - j2)^2))
        end
    end
    H = sparse(iH, jH, sH)
    res.H = H
    res.Hs = sum(H, dims=2)
    nothing
end

function set_filter!(res::TopOptResult, prob::TopOptProblemConv, rmin)
    nx = prob.nx
    ny = prob.ny

    kern_size = length(-ceil(rmin) + 1:ceil(rmin) - 1)
    dy = -ceil(rmin) + 1:ceil(rmin) - 1
    dx = -ceil(rmin) + 1:ceil(rmin) - 1

    dy = dy' .* ones(kern_size)
    dx = ones(kern_size)' .* dx

    res.h = max.(0, rmin .- sqrt.(dx.^2 + dy.^2))
    res.Hs = fastconv_2(ones(ny, nx), res.h)
    nothing
end

# Optimal criterion for the normal filter
function optimal_crit(prob::TopOptProblemBase, res::TopOptResult, β::Float64,
                      ∂c::Array{Float64}, ∂V::Array{Float64}, filter)
    nx = prob.nx
    ny = prob.ny
    λ₁::Float64 = 0.0
    λ₂::Float64 = 1e9
    move::Float64 = 0.2
    x_new = similar(res.x)
    tot_volfrac = porb.volfrac * nx * ny

    while (λ₂ - λ₁) / (λ₁ + λ₂) > 1e-3
        λₘ = 0.5 * (λ₂ + λ₁);
        xₑBₑ = @. res.x * sqrt(-∂c / ∂V / λₘ)
        x_new .= criteria.(res.x, move, xₑBₑ)
        if filter == 1
            res.x̃ = x_new
        elseif filter == 2
            res.x̃[:] = (res.H * view(x_new, :)) ./ res.Hs
        elseif filter == 3
            res.x̂[:] = (res.H * view(x_new, :)) ./ res.Hs
            res.x̃[:] = @. 1 - exp(-β * res.x̂) + res.x̂ * exp(-β)
        end

        # Consider passibe elements
        res.x̃[prob.passive .== 1] .= 0
        res.x̃[prob.passive .== 2] .= 1

        if sum(res.x̃) > tot_volfrac
            λ₁ = λₘ
        else
            λ₂ = λₘ
        end
    end

    return x_new
end

# Optimal criterion using convolution
function optimal_crit(prob::TopOptProblemConv, res::TopOptResult, β::Float64,
                      ∂c::Array{Float64}, ∂V::Array{Float64}, filter)
    nx = prob.nx
    ny = prob.ny
    λ₁::Float64 = 0.0
    λ₂::Float64 = 1e9
    move::Float64 = 0.2
    x_new = similar(res.x)
    tot_volfrac = prob.volfrac * nx * ny

    while (λ₂ - λ₁) / (λ₁ + λ₂) > 1e-3
        λₘ = 0.5 * (λ₂ + λ₁);
        # x_new .= max.(0, max.(x .- move, min.(1, min.(x .+ move, x .* sqrt.(-∂c ./ ∂V / λₘ)))))
        xₑBₑ = @. res.x * sqrt(-∂c / ∂V / λₘ)
        x_new .= criteria.(res.x, move, xₑBₑ)
        if filter == 1
            res.x̃ = x_new
        elseif filter == 2
            res.x̃[:] = fastconv_2(x_new, res.h) ./ res.Hs
        elseif filter == 3
            res.x̂[:] = fastconv_2(x_new, res.h) ./ res.Hs
            res.x̃[:] = @. 1 - exp(-β * res.x̂) + res.x̂ * exp(-β)
        end

        res.x̃[prob.passive .== 1] .= 0
        res.x̃[prob.passive .== 2] .= 1

        if sum(res.x̃) > tot_volfrac
            λ₁ = λₘ
        else
            λ₂ = λₘ
        end
    end

    return x_new
end

function criteria(x, move, xₑBₑ)
    if xₑBₑ ≤ max(0.0, x - move)
        x_new = max(0.0, x - move)
    elseif xₑBₑ ≥ min(1.0, x + move)
        x_new = min(1.0, x + move)
    else
        x_new = xₑBₑ
    end
    return x_new
end

function solve_fe!(fe::FEProblem, x̃::Array, p, E₀, E_min)
    sK = reshape(fe.Kₑ[:] * (E_min .+ x̃[:]'.^p .* (E₀ - E_min)), 64 * fe.nx * fe.ny)
    # Global stiffness matrix
    K = sparse(fe.iK, fe.jK, sK)

    freedofs = fe.freedofs
    F = fe.F
    U = fe.U
    # Solve system
    fe.U[freedofs] = Symmetric(K[freedofs, freedofs]) \ F[freedofs]
    nothing
end

function compliance(fe::FEProblem, E₀, E_min)
    nx = fe.nx
    ny = fe.ny
    edofMat = fe.edofMat

    cₑ = zeros(ny, nx)

    for ely in 1:ny, elx in 1:nx
        nₑ = ely + ny * (elx - 1)
        cₑ[ely, elx] = sum(fe.U[view(edofMat, nₑ, :)]' * fe.Kₑ .* fe.U[view(edofMat, nₑ, :)]')
    end

    return cₑ
end

function sensitivities!(∂c, ∂V, fe::FEProblem, res::TopOptResultBase, β, p, E₀, E_min, cₑ, filter)
    nx = fe.nx
    ny = fe.ny
    Hs = res.Hs
    H = res.H

    ∂c .= -p .* (E₀ - E_min) .* res.x̃.^(p - 1) .* cₑ
    ∂V .= ones(ny, nx);
    # Filtering/modification of sensitivities
    if filter == 1
        ∂c[:] = H * (res.x[:] .* ∂c[:]) ./ Hs ./ max.(1e-3, res.x)
    elseif filter == 2
        ∂c[:] = H * (∂c[:] ./ Hs)
        ∂V[:] = H * (∂V[:] ./ Hs)
    elseif filter == 3
        ∂x = @. β * exp(-β * x̄) + exp(-β)
        ∂c[:] = H * (∂c[:] .* ∂x[:] ./ Hs)
        ∂V[:] = H * (∂V[:] .* ∂x[:] ./ Hs)
    end

    return ∂c, ∂V
end

function sensitivities!(∂c, ∂V, fe::FEProblem, res::TopOptResultConv, β, p, E₀, E_min, cₑ, filter)

    nx = fe.nx
    ny = fe.ny
    Hs = res.Hs
    h = res.h

    ∂c .= -p .* (E₀ - E_min) .* res.x̃.^(p - 1) .* cₑ
    ∂V .= ones(ny, nx);
    # Filtering/modification of sensitivities
    if filter == 1
        ∂c[:] = fastconv_2(∂c .* res.x̃, res.h) ./ res.Hs ./ max.(1e-3, res.x̃)
    elseif filter == 2
        ∂c[:] = fastconv_2(∂c ./ res.Hs, res.h)
        ∂V[:] = fastconv_2(∂V ./ res.Hs, res.h)
    elseif filter == 3
        ∂x = @. β * exp(-β * res.x̂) + exp(-β)
        ∂c[:] = fastconv_2(∂c .* ∂x ./ res.Hs, res.h)
        ∂V[:] = fastconv_2(∂V .* ∂x ./ res.Hs, res.h)
    end

    return ∂c, ∂V
end

function element_K(E₀::Float64, ν::Float64)
    A₁₁ = [12  3 -6 -3;  3 12  3  0; -6  3 12 -3; -3  0 -3 12]
    A₁₂ = [-6 -3  0  3; -3 -6 -3 -6;  0 -3 -6  3;  3 -6  3 -6]
    B₁₁ = [-4  3 -2  9;  3 -4 -9  4; -2 -9 -4 -3;  9  4 -3 -4]
    B₁₂ = [ 2 -3  4 -9; -3  2  9 -2;  4  9  2  3; -9 -2  3  2]
    Kₑ = 1 / (1 - ν^2) / 24 .* ([A₁₁ A₁₂; A₁₂' A₁₁] + ν * [B₁₁ B₁₂; B₁₂' B₁₁]);
    return Kₑ
end

end

max_stars_repo_name: cristobaltapia/topopt-julia

max_issues_repo_name: cristobaltapia/topopt-julia

max_forks_repo_name: cristobaltapia/topopt-julia


text:
using DataFrames
#using DataFramesMeta

include(joinpath(@__DIR__, "lib.jl"))
block_devices = TableDirectories.TableRoot("/sys/dev/block/")
partition = TableDirectories.RowRoot(block_devices, "8:1")
file_value = TableDirectories.Cell(partition, "power")
sparse_row = TableDirectories.RowRoot(block_devices, "11:0")
sparse_column = TableDirectories.Cell(sparse_row, "partition")
#TableDirectories.cell_value(file_value)

columns_dir = joinpath(TableDirectories.path(partition), "")
all_column_names = TableDirectories.column_names(partition)
TableDirectories.dict(partition, all_column_names)
data_frame = TableDirectories.data_frame(block_devices, all_column_names)

top_regex = r"^[-0-9]+$"
bus_regex = r"usb([1-9]+)"
endpoint_regex = r"^[-0-9.]+:.+$"
live_usb_devices = TableDirectories.data_frame(TableDirectories.TableRoot("/sys/bus/usb/devices", top_regex))
endpoint_data_frame = TableDirectories.data_frame(TableDirectories.TableRoot("/sys/bus/usb/devices", endpoint_regex))
#@show TableDirectories.discriminant(endpoint_data_frame)
#@show live_usb_devices[!, TableDirectories.unique_column_subpaths(live_usb_devices)]
live_usb_devices.port = getindex.(splitpath.(live_usb_devices.port), 8)
live_usb_devices[!, ["devnum", "devpath", "dev", "port"]]
#TableDirectories.default_columns(live_usb_devices)
eth0_dataframe = TableDirectories.data_frame(TableDirectories.TableRoot("/sys/class/net", r".+"))


text:
using GeometryTypes
using SurfaceTopology
using LinearAlgebra

function quadraticform(vects,vnormal)

    Lx = [0 0 0; 0 0 -1; 0 1 0]
    Ly = [0 0 1; 0 0 0; -1 0 0]
    Lz = [0 -1 0; 1 0 0; 0 0 0]

    d = [0,0,1] + vnormal
    d /= norm(d)

    Ln = d[1]*Lx + d[2]*Ly + d[3]*Lz
    R = exp(pi*Ln)

    vects = copy(vects)
    for vj in 1:length(vects)
        vects[vj] = R*vects[vj]
    end

    ### Construction of the system
    A = Array{Float64}(undef,3,3)
    B = Array{Float64}(undef,3)

    vects_norm2 = Array{Float64}(undef,length(vects))
    for vj in 1:length(vects)
       vects_norm2[vj] = norm(vects[vj])^2
    end

    A[1,1] = sum((v[1]^4 for v in vects) ./ vects_norm2)
    A[1,2] = sum((v[1]^3*v[2] for v in vects) ./ vects_norm2)
    A[1,3] = sum((v[1]^2*v[2]^2 for v in vects) ./ vects_norm2)
    A[2,1] = A[1,2]
    A[2,2] = A[1,3]
    A[2,3] = sum( (v[2]^3*v[1] for v in vects) ./vects_norm2)
    A[3,1] = A[1,3]
    A[3,2] = A[2,3]
    A[3,3] = sum((v[2]^4 for v in vects) ./vects_norm2)

    B[1] = sum((v[3]*v[1]^2 for v in vects) ./vects_norm2)
    B[2] = sum((v[1]*v[2]*v[3] for v in vects) ./vects_norm2)
    B[3] = sum((v[2]^2*v[3] for v in vects) ./vects_norm2)

    C,D,E = A\B
    return C,D,E
end

function meancurvature(points,topology)
    curvatures = Array{Float64}(undef,length(points))
    for v in 1:length(points)

        s = Point(0,0,0)
        for (v1,v2) in EdgeRing(v,topology)
            s += cross(points[v2],points[v1])
        end
        normal = s ./ norm(s)

        vring = collect(VertexRing(v,topology))
        vects = [points[vi] - points[v] for vi in vring]

        C,D,E = quadraticform(vects,normal)

        A = [C D/2;D/2 E]
        k1,k2 = eigvals(-A)
        H = (k1 + k2)/2

        curvatures[v] = H
    end
    return curvatures
end

### Testing
t = ( 1 + sqrt( 5 ) ) / 2;

vertices = Point{3,Float64}[
    [ -1,  t,  0 ], [  1, t, 0 ], [ -1, -t,  0 ], [  1, -t,  0 ],
    [  0, -1,  t ], [  0, 1, t ], [  0, -1, -t ], [  0,  1, -t ],
    [  t,  0, -1 ], [  t, 0, 1 ], [ -t,  0, -1 ], [ -t,  0,  1 ]
] ./ sqrt(1 + t^2)

faces = Face{3,Int64}[
    [1, 12, 6], [1, 6, 2], [1, 2, 8], [1, 8, 11], [1, 11, 12], [2, 6, 10], [6, 12, 5],
    [12, 11, 3], [11, 8, 7], [8, 2, 9], [4, 10, 5], [4, 5, 3], [4, 3, 7], [4, 7, 9],
    [4, 9, 10], [5, 10, 6], [3, 5, 12], [7, 3, 11], [9, 7, 8], [10, 9, 2]
]

### As paraboloid grows faster than sphere the estimated curvature for a coarse mesh would be lower
curvatures = meancurvature(vertices,faces)