text:
module TensegrityEquilibria

export stableEquilibria,
       start_demo,
       animateTensegrity

import LinearAlgebra: norm, nullspace, zeros, eigvals
import HomotopyContinuation: @var, Expression, evaluate, solve, differentiate, System, InterpretedSystem, randn, target_parameters!, real_solutions, nparameters, jacobian!, monodromy_solve, Variable, solutions, ParameterHomotopy, solver_startsolutions, parameters, variables
import GLMakie: scatter!, Node, @lift, limits!, linesegments!, record, Point3f0, on, Point2f0, Point, FRect, Scene, cam3d!, xlims!, ylims!, zlims!
import GLMakie.AbstractPlotting: layoutscene, labelslidergrid!, Box, Label, LScene, MakieLayout, set_close_to!, Axis
import HomotopyOpt

function energyfunc(currentvalue, unknownvalues, vertices, parameters, targetparameters, unknownCables)
    Q=0.0
    p = verticesToMatrix(currentvalue, vertices, unknownvalues)
    for cable in unknownCables
        i,j,r,e = Int64(cable[1]), Int64(cable[2]), Float64(cable[3]), Float64(cable[4]) # unpack
        cablelength = norm(p[i,:] - p[j,:])
        Q += 1/2 * e * (cablelength - r)^2
    end
    return(Q)
end

function verticesToMatrix(currentvalue, vertices, unknownvalues)
    output = zeros(typeof(currentvalue[1]),length(vertices), length(vertices[1]))
    for i in 1:size(output)[1], j in 1:size(output)[2]
        if typeof(vertices[i][j])!=Float64 && typeof(vertices[i][j])!=Int64
            for k in 1:length(unknownvalues)
                if unknownvalues[k]==vertices[i][j]
                    output[i,j] = currentvalue[k]
                end
            end
        else
            output[i,j] = vertices[i][j]
        end
    end
    return(output)
end
#=
@input vertices::[[p_11,p_12,...], ..., [p_m1,p_m2,...]], unknownBars::[[i,j,l_ij],...] non-redundant list, unknowncables::[[i,j,r_ij,e_ij]] non-redundant list,
listOfInternalVariables::[a,...], listOfControlParams::[a,...], targetParams::[a,...], knownBars::[[i,j],...] for plotting, knownCables::[[i,j],...] for plotting.
If an animation is wished, use the optional argument timestamps::[[q_11,q_12,...],...] to provide a path (given by points) in the space listOfControlParams.
Calculate configurations in equilibrium for the given tensegrity framework
=#
function stableEquilibria(vertices::Array, unknownBars::Array, unknownCables::Array, listOfInternalVariables::Array{Variable,1}, listOfControlParams::Array, targetParams::Array, knownBars::Array, knownCables::Array, timestamps=[]; thresholdForIteration = 30, thresholdForEDStep=6)
    assertCorrectInput(vertices, unknownBars, unknownCables, listOfInternalVariables, listOfControlParams, targetParams, knownBars, knownCables)
    equilibriaarray = []
    if length(listOfInternalVariables)>=thresholdForIteration
        B = []
        for index in 1:length(unknownBars)
            l_ij=unknownBars[index][3]; p_i=vertices[Int64(unknownBars[index][1])]; p_j=vertices[Int64(unknownBars[index][2])];
            bar = !isempty(listOfControlParams) ? Expression(evaluate(sum((p_i-p_j).^2)-l_ij^2, listOfControlParams=>targetParams)) : Expression(sum((p_i-p_j).^2)-l_ij^2)
            push!(B, bar)
        end
        B = Vector{Expression}(B)
        ConV = HomotopyOpt.ConstraintVariety(listOfInternalVariables, B, length(listOfInternalVariables), length(listOfInternalVariables)-length(B), 2^(length(vertices[1]))*6)
        Q = x->energyfunc(x, listOfInternalVariables, vertices, listOfControlParams, targetParams, unknownCables)
        for q in ConV.samples
        	resultminimum = HomotopyOpt.findminima(q, 1e-4, ConV, Q; maxseconds=1000, whichstep="gaussnewtonstep", initialstepsize=0.3);
        	if !any( point -> norm(resultminimum.computedpoints[end]-point)<0.01, equilibriaarray ) && resultminimum.lastpointisminimum
        	    push!(equilibriaarray,resultminimum.computedpoints[end])
        	end
        end
        realSol = plotStaticFramework(vertices, equilibriaarray, unknownBars, knownBars, unknownCables, knownCables, listOfInternalVariables, listOfControlParams, targetParams)
        return(realSol, listOfInternalVariables)
    elseif length(listOfInternalVariables)>=thresholdForEDStep
        B = []
        for index in 1:length(unknownBars)
            l_ij=unknownBars[index][3]; p_i=vertices[Int64(unknownBars[index][1])]; p_j=vertices[Int64(unknownBars[index][2])];
            bar = !isempty(listOfControlParams) ? Expression(evaluate(sum((p_i-p_j).^2)-l_ij^2, listOfControlParams=>targetParams)) : Expression(sum((p_i-p_j).^2)-l_ij^2)
            push!(B, bar)
        end
        B = Vector{Expression}(B)
        ConV = HomotopyOpt.ConstraintVariety(listOfInternalVariables, B, length(listOfInternalVariables), length(listOfInternalVariables)-length(B), 2^(length(vertices[1]))*6)
        Q = x->energyfunc(x, listOfInternalVariables, vertices, listOfControlParams, targetParams, unknownCables)
        for q in ConV.samples
        	resultminimum = HomotopyOpt.findminima(q, 1e-4, ConV, Q; maxseconds=1000, whichstep="EDStep", initialstepsize=0.3);
        	if !any( point -> norm(resultminimum.computedpoints[end]-point)<0.01, equilibriaarray ) && resultminimum.lastpointisminimum
        	    push!(equilibriaarray,resultminimum.computedpoints[end])
        	end
        end
        realSol = plotStaticFramework(vertices, equilibriaarray, unknownBars, knownBars, unknownCables, knownCables, listOfInternalVariables, listOfControlParams, targetParams)
        return(realSol, listOfInternalVariables)
    end
    @var delta[1:length(unknownCables)]
    B, C, G = [], [], []
    for index in 1:length(unknownBars)
        l_ij=unknownBars[index][3]; p_i=vertices[Int64(unknownBars[index][1])]; p_j=vertices[Int64(unknownBars[index][2])];
        push!(B, sum((p_i-p_j).^2)-l_ij^2)
    end
    for index in 1:length(unknownCables)
        r_ij=unknownCables[index][3]; p_i=vertices[Int64(unknownCables[index][1])]; p_j=vertices[Int64(unknownCables[index][2])]; e_ij=unknownCables[index][4];
        push!(C, sum((p_i-p_j).^2)-delta[index]^2)
    end
    append!(G,B); append!(G,C);
    e = map(t->t[4], unknownCables); r = map(t->t[3], unknownCables);
    Q = sum(e.*(r-delta).^2/2)

    @var lambda[1:length(G)]
    L = Q + lambda'*G
    ∇L = differentiate(L, [listOfInternalVariables; delta; lambda])
    catastropheWitness = catastrophePoints([entry for entry in vertices], listOfInternalVariables, listOfControlParams, targetParams, unknownCables, unknownBars)
    if(!isempty(listOfInternalVariables))
        F = isempty(listOfControlParams) ? System(∇L; variables = [listOfInternalVariables; delta; lambda]) : System(∇L; variables = [listOfInternalVariables; delta; lambda], parameters = listOfControlParams)
        params₀ = randn(ComplexF64, nparameters(F))
        #TODO compute a local solution via gradient descent
        res₀ = isempty(listOfControlParams) ? solve(F) : solve(F; target_parameters = params₀)

        S₀ = solutions(res₀)
        H = ParameterHomotopy(F, params₀, params₀)
        solver, _ = solver_startsolutions(H, S₀)
        if(isempty(timestamps))
            # produces a non-animate, interactive plots
            realSol = plotStaticFramework(vertices, unknownBars, knownBars, unknownCables, knownCables, solver, S₀, listOfInternalVariables, listOfControlParams, targetParams, delta, lambda, L, G, catastropheWitness)
            return(realSol, listOfInternalVariables)
        else
            # produces an animation
            realSol = animateSolvedFramework(vertices, unknownBars, knownBars, unknownCables, knownCables, solver, S₀, listOfInternalVariables, listOfControlParams, targetParams, delta, lambda, L, G, catastropheWitness, timestamps)
            return(realSol, listOfInternalVariables)
        end
    else
        throw(error("Internal variables need to be provided!"))
    end
    return(nothing)
end

#= Check if the input configuration (a critical point) is a local minimum of the energy function or not.
=#
function isLocalMinimum(listOfInternalVariables, listOfControlParams, delta, lambda, L, G)
    ∇L = differentiate(L, [listOfInternalVariables; delta])
    if(isempty(listOfControlParams))
        HL = InterpretedSystem(System(∇L, [listOfInternalVariables; delta], lambda))
        dG = InterpretedSystem(System(G, [listOfInternalVariables; delta], lambda))
    else
        HL = InterpretedSystem(System(∇L, [listOfInternalVariables; delta], [lambda; listOfControlParams]))
        dG = InterpretedSystem(System(G, [listOfInternalVariables; delta], [lambda; listOfControlParams]))
    end

    (s, params) -> begin
      v = s[1:length(delta)+length(listOfInternalVariables)]
      q = [s[1+length(delta)+length(listOfInternalVariables):end]; params]
      W = zeros(length(delta)+length(listOfInternalVariables), length(delta)+length(listOfInternalVariables))
      jacobian!(W, HL, v, q)
      J = zeros(length(delta)+length(listOfInternalVariables)-1, length(delta)+length(listOfInternalVariables))
      jacobian!(J, dG, v, q)
      V = nullspace(J)
      all(e -> e ≥ 1e-14, eigvals(V' * W * V))
    end
end

#= assert that the inputs of stableEquilibria are as expected
=#
function assertCorrectInput(vertices, unknownBars, unknownCables, listOfInternalVariables, listOfControlParams, targetParams, knownBars, knownCables)
    @assert(length(vertices)>1)
    D = length(vertices[1])
    foreach(vx->@assert(length(vx)==D), vertices)
    foreach(bar->@assert(length(bar)==3), unknownBars)
    foreach(cable->@assert(length(cable)==4), unknownCables)
    foreach(bar->@assert(length(bar)==2), knownBars)
    foreach(cable->@assert(length(cable)==2), knownCables)
    foreach(var->@assert(typeof(var)==Variable), listOfInternalVariables)
    foreach(par->@assert(typeof(par)==Variable), listOfControlParams)
    @assert(length(listOfControlParams)==length(targetParams))
end

#= This function extracts the position of the parameters in the vertices and replaces them
 with their respective computed value. It returns the array of vertices in the form of Point2f0/Point3f0
 and the array of all other minimizers of the energy function Q, that were not yet used, in shadowPoints.
 WARNING: This cannot deal with Polynomial-type entries of the vertices yet. Only pure variables are allowed.=#
function arrangeArray(array, listOfInternalVariables, realSol, listOfControlParams, targetParams)
    #TODO add polynomial substitution for coordinates on a curve (polynomials)
    subsArray=[[] for sol in realSol]
    for i in 1:length(array)
        helper=[];
        for j in 1:length(array[i])
            # We check, if the current array element is a number or a variable
            if(typeof(array[i][j])!=Float64 && typeof(array[i][j])!=Int64 && typeof(array[i][j])!=ComplexF64)
                index=findfirst(en->en==array[i][j], listOfInternalVariables)
                if(index==nothing || typeof(index) == Nothing)
                    index=findfirst(en->en==array[i][j], listOfControlParams)
                    if(index==nothing || typeof(index) == Nothing)
                        throw(error("The variable in the array is neither a control nor an internal parameter."))
                    end
                    push!(helper,targetParams[index])
                else
                    push!(helper,[realSol[i][index] for i in 1:length(realSol)])
                end
            else
                push!(helper,array[i][j])
            end
        end
        # All possible configurations are added to an array. Later, one is chosen to be the vertex to be displayed and the others are chosen to be shadow points.
        for i in 1:length(realSol)
            if length(helper)==2
                push!(subsArray[i], Point2f0(typeof(helper[1])==Array{Float64,1} ? helper[1][i] : helper[1], typeof(helper[2])==Array{Float64,1} ? helper[2][i] : helper[2]))
            else
                push!(subsArray[i], Point3f0(typeof(helper[1])==Array{Float64,1} ? helper[1][i] : helper[1], typeof(helper[2])==Array{Float64,1} ? helper[2][i] : helper[2], typeof(helper[3])==Array{Float64,1} ? helper[3][i] : helper[3]))
            end
        end
    end
    return(subsArray)
end

#= Creates a static parametric plot of the vertices (red) with corresponding bars (black) and cables (blue) for larget systems where HC.jl fails.
The other minima are displayed in grey.
The method allows for swapping between the different local minima of the energy function with the 'n' key.=#
function plotStaticFramework(vertices, solutions, unknownBars, knownBars, unknownCables, knownCables, listOfInternalVariables, listOfControlParams, targetParams)
    bars=vcat(unknownBars,knownBars); cables=vcat(unknownCables,knownCables)
    # Make the variable params interactive.
    params=Node(targetParams)

    scene, layout = layoutscene(resolution = (1400, 850))
    ax = layout[1:4, 1] = length(vertices[1])==3 ? LScene(scene, width=750, height=750, camera = cam3d!, raw = false) : MakieLayout.Axis(scene,width=750,height=750)

    # Initialization of the scene's layout.
    layout[1:4, 2] = Box(scene, color = :white, strokecolor = :transparent, width=50)
    layout[1, 3] = Box(scene, color = :white, strokecolor = :transparent, height=200)
    layout[2, 3] = Box(scene, color = (:white, 0.1), strokecolor = :red, height=80, width=400)
    layout[2, 3] = Label(scene, "Press the 'n' Key to iterate through\nthe different vertex configurations", textsize = 20, halign=:center, valign=:center, color=:teal)
    layout[3, 3] = Box(scene, color = (:white, 0.1), strokecolor = :red, height=80, width=400)
    layout[3, 3] = Label(scene, "Due to the given framework's size,\n only some equilibria are displayed.", textsize = 20, halign=:center, valign=:center, color=:teal)
    layout[4, 3] = Box(scene, color = :white, strokecolor = :transparent, height=200)

    solutions = arrangeArray(vertices, listOfInternalVariables, solutions, listOfControlParams, [])

    # Make the index of the current configuration in the space of all possible configurations `allVertices` interactively choosable. Change it during runtime by pressing 'n'
    currentIndex = Node(1)
    on(scene.events.unicode_input) do input
        'n' in input ? (currentIndex[]-1==0 ? currentIndex[]=length(solutions) : currentIndex[]=currentIndex[]-1) : nothing
    end
    fixedVertices=@lift((solutions)[$(currentIndex) > length(solutions) ? 1 : $(currentIndex)])
    # Plot all given edges and vertices. They are modified upon registering change in `allVertices`, which is equivalent to change in params.
    foreach(bar->linesegments!(ax, @lift([($fixedVertices)[Int64(bar[1])], ($fixedVertices)[Int64(bar[2])]]) ; linewidth = length(vertices[1])==2 ? 4.0 : 5.0, color=:black), bars)
    foreach(cable->linesegments!(ax, @lift([($fixedVertices)[Int64(cable[1])], ($fixedVertices)[Int64(cable[2])]]) ; color=:blue), cables)
    foreach(solution -> scatter!(ax, [p for p in solution]; color=:goldenrod3, marker = :diamond, alpha=0.1, markersize = length(vertices[1])==2 ? 12 : 75), solutions)
    scatter!(ax, @lift([p for p in $fixedVertices]); color=:red, markersize = length(vertices[1])==2 ? 12 : 75)
    if(length(vertices[1])==2)
        xlimiter = Node([Inf,-Inf]); xlimiter = @lift(computeMinMax($fixedVertices, solutions, $xlimiter, 1));
        ylimiter = Node([Inf,-Inf]); ylimiter = @lift(computeMinMax($fixedVertices, solutions, $ylimiter, 2));
        @lift(limits!(ax, FRect((($xlimiter)[1]-0.5,($ylimiter)[1]-0.5), (($xlimiter)[2]-($xlimiter)[1]+1.0,($ylimiter)[2]-($ylimiter)[1]+1.0))))
    end
    display(scene)
    return(fixedVertices[])
end

#= Creates a dynamic parametric plot of the vertices (red) with corresponding bars (black) and cables (blue).
The shadow vertices are plotted in grey. If there are no parameters given, this plot is static.
It allows for swapping between the different local minima of the energy function with the 'n' key.=#
function plotStaticFramework(vertices, unknownBars, knownBars, unknownCables, knownCables, solver, S₀, listOfInternalVariables, listOfControlParams, targetParams, delta, lambda, L, G, catastrophePoints)
    bars=vcat(unknownBars,knownBars); cables=vcat(unknownCables,knownCables)
    # Make the variable params interactive.
    params=Node(targetParams)

    scene, layout = layoutscene(resolution = (1400, 850))
    ax = layout[1:4, 1] = length(vertices[1])==3 ? LScene(scene, width=750, height=750, camera = cam3d!, raw = false) : MakieLayout.Axis(scene,width=750,height=750)

    # Initialization of the scene's layout.
    layout[1:4, 2] = Box(scene, color = :white, strokecolor = :transparent, width=50)
    layout[1, 3] = Box(scene, color = :white, strokecolor = :transparent, height=50)
    layout[2, 3] = Box(scene, color = (:white, 0.1), strokecolor = :red, height=80, width=400)
    layout[2, 3] = Label(scene, "Press the 'n' Key to iterate through\nthe different vertex configurations", textsize = 20, halign=:center, valign=:center, color=:teal)
    layout[3, 3] = Box(scene, color = (:white, 0.1), strokecolor = :transparent, height=100)
    layout[3, 3] = Label(scene, "Control Parameters:", textsize = 20, halign=:left, valign=:bottom)

    # Adding sliders for interactively choosing the parameters during runtime.
    if length(listOfControlParams)>0
        lsgrid = labelslidergrid!(
            scene,
            [string(string(para),": ") for para in listOfControlParams],
            [targetParams[i]-2:0.05:targetParams[i]+2 for i in 1:length(targetParams)];
            formats = x -> "$(x==(round(x, digits = 1)) ? string(string(x),"0") : string(x))",
            width = 400,
            valign=:top
        )
        layout[4,3]=lsgrid.layout
        for i in 1:length(listOfControlParams)
            set_close_to!(lsgrid.sliders[i], targetParams[i])
            params = @lift begin
                param = ($params)
                param[i]=$(lsgrid.sliders[i].value)
                return(param)
            end
        end
    end

    # Upon change of params, recalculate the possible vertex configurations of this system
    allVertices=@lift begin
        target_parameters!(solver, $params)
        res₀ = solve(solver, S₀, threading = true)
        realSol = (map(t->t[1:length(listOfInternalVariables)+length(delta)], filter(en->isempty(filter(x->x<0, en[length(listOfInternalVariables)+1:length(listOfInternalVariables)+length(delta)])) &&
            isLocalMinimum(listOfInternalVariables, listOfControlParams, delta, lambda, L, G)(real.(en), $params), real_solutions(res₀))))
        arrangeArray(vertices, listOfInternalVariables, [sol[1:length(listOfInternalVariables)] for sol in realSol], listOfControlParams, $params)
    end

    # Make the index of the current configuration in the space of all possible configurations `allVertices` interactively choosable. Change it during runtime by pressing 'n'
    currentIndex = Node(1)
    on(scene.events.unicode_input) do input
        'n' in input ? (currentIndex[]-1==0 ? currentIndex[]=length(allVertices[]) : currentIndex[]=currentIndex[]-1) : nothing
    end
    fixedVertices=@lift(($allVertices)[$(currentIndex) > length($allVertices) ? 1 : $(currentIndex)])
    shadowPoints=@lift(vcat(($allVertices)[1:end]...))

    # Plot all given edges and vertices. They are modified upon registering change in `allVertices`, which is equivalent to change in params.
    foreach(bar->linesegments!(ax, @lift([($fixedVertices)[Int64(bar[1])], ($fixedVertices)[Int64(bar[2])]]) ; linewidth = length(vertices[1])==2 ? 4.0 : 5.0, color=:black), bars)
    foreach(cable->linesegments!(ax, @lift([($fixedVertices)[Int64(cable[1])], ($fixedVertices)[Int64(cable[2])]]) ; color=:blue), cables)
    scatter!(ax, @lift([f for f in $shadowPoints]); color=:goldenrod3, marker = :diamond, alpha=0.1, markersize = length(vertices[1])==2 ? 12 : 75)
    scatter!(ax, @lift([f for f in $fixedVertices]); color=:red, markersize = length(vertices[1])==2 ? 12 : 75)
    #mesh!(ax, catastrophePoints, color=:lightgrey, alpha=0.15)
    !isempty(catastrophePoints) ? scatter!(ax, catastrophePoints; alpha=0.1, markersize=0.5, color=:lightgrey) : nothing
    if(length(vertices[1])==2)
        xlimiter = Node([Inf,-Inf]); xlimiter = @lift(computeMinMax($fixedVertices, $shadowPoints, $xlimiter, 1));
        ylimiter = Node([Inf,-Inf]); ylimiter = @lift(computeMinMax($fixedVertices, $shadowPoints, $ylimiter, 2));
        @lift(limits!(ax, FRect((($xlimiter)[1]-0.5,($ylimiter)[1]-0.5), (($xlimiter)[2]-($xlimiter)[1]+1.0,($ylimiter)[2]-($ylimiter)[1]+1.0))))
    end
    display(scene)
    return(fixedVertices[])
end

#= This method omputes the minimal and maximal element of the position xyz a nested array (so all values of the form
array[:][index]). Afterwards it checks if the previous minimum limiter[1] or the maximum limter[2] is beaten. If so,
it returns the new optima.=#
function computeMinMax(fixedVertices,shadowPoints,limiter,xyz)
    for i in 1:length(fixedVertices)
        fixedVertices[i][xyz] > limiter[2] ? limiter[2]=fixedVertices[i][xyz] : nothing
        fixedVertices[i][xyz] < limiter[1] ? limiter[1]=fixedVertices[i][xyz] : nothing
    end
    for i in 1:length(shadowPoints)
        shadowPoints[i][xyz] > limiter[2] ? limiter[2]=shadowPoints[i][xyz] : nothing
        shadowPoints[i][xyz] < limiter[1] ? limiter[1]=shadowPoints[i][xyz] : nothing
    end
    return(limiter)
end

#= The method catastrophePoints computes a witness set of the catastrophe discriminant using monodromy loops
 If the computation fails for some reason, it returns nothing, so the plot can still be displayed.=#
function catastrophePoints(vertices, internalVariables, controlParameters, targetParams, unknownCables, unknownBars)
    cp = let
        #TODO find control node
        @var delta[1:length(unknownCables)]
        B, C, G = [], [], []
        for index in 1:length(unknownBars)
            l_ij=unknownBars[index][3]; p_i=vertices[Int64(unknownBars[index][1])]; p_j=vertices[Int64(unknownBars[index][2])];
            push!(B, sum((p_i-p_j).^2)-l_ij^2)
        end
        for index in 1:length(unknownCables)
            r_ij=unknownCables[index][3]; p_i=vertices[Int64(unknownCables[index][1])]; p_j=vertices[Int64(unknownCables[index][2])]; e_ij=unknownCables[index][4];
            push!(C, sum((p_i-p_j).^2)-delta[index]^2)
        end
        append!(G,B); append!(G,C);
        e = map(t->t[4], unknownCables); r = map(t->t[3], unknownCables);
        Q = sum(e.*(r-delta).^2/2)

        @var lambda[1:length(G)]
        L = Q + lambda'*G
        ∇L = differentiate(L, [internalVariables; delta; lambda])
        @var v[1:length(∇L)]
        α = randn(length(∇L))
        J_L_v = [differentiate(∇L, [internalVariables; delta; lambda]) * v; α'v - 1]
        @var a[1:length(controlParameters),1:length(controlParameters)-1] b[1:length(controlParameters)-1]
        L1 = a' * controlParameters + b
        try
            P = System(
              [∇L; J_L_v; L1];
              variables = [controlParameters; internalVariables; delta; lambda; v],
              parameters = [collect(Iterators.flatten(a)); b]
            )
            res = monodromy_solve(P)
            #res=solve(P, target_parameters=randn(ComplexF64, nparameters(P)))
            rand_lin_space = let
                () -> randn(nparameters(P))
            end
            N = 500
            alg_catastrophe_points = solve(
                P,
                solutions(res),
                start_parameters = parameters(res),#startParams,
                target_parameters = [rand_lin_space() for i = 1:N],
                transform_result = (r, p) -> real_solutions(r),
                flatten = true
            )
            filter!(p -> all(k->k ≥ 0, p[length(controlParameters)+length(internalVariables)+1:length(controlParameters)+length(internalVariables)+length(delta)]), alg_catastrophe_points);
            cp = map(p -> length(vertices[1])==2 ? Point2f0(p[1], p[2]) : Point3f0(p[1], p[2], p[3]), alg_catastrophe_points)
            return(cp)
        catch e
            display(e)
        end
        return([])
    end
    return(cp)
end

#= This method allows for the animation of a tensegrity framework. Given a path in parameter space defined by point samples, it animates the movement
 of the energy functions' minima along this path.=#
function animateSolvedFramework(vertices, unknownBars, knownBars, unknownCables, knownCables, solver, S₀, listOfInternalVariables, listOfControlParams, targetParams, delta, lambda, L, G, cp, timestamps)
    bars=vcat(unknownBars,knownBars); cables=vcat(unknownCables,knownCables)
    params=Node(targetParams)
    scene, layout = layoutscene(resolution = (850, 850))
    ax = layout[1:length(listOfControlParams)+2, 1] = length(vertices[1])==3 ? LScene(scene, width=750, height=750, camera = cam3d!, raw = false) : Axis(scene,width=750,height=750)

    allVertices=@lift begin
        target_parameters!(solver, $params)
        res₀ = solve(solver, S₀, threading = true)
        realSol = (map(t->t[1:length(listOfInternalVariables)], filter(en->isempty(filter(x->x<0, en[length(listOfInternalVariables)+1:length(listOfInternalVariables)+length(delta)])) &&
            isLocalMinimum(listOfInternalVariables, listOfControlParams, delta, lambda, L, G)(real.(en), $params), real_solutions(res₀))))
        arrangeArray(vertices, listOfInternalVariables, realSol, listOfControlParams, $params)
    end
    currentIndex=Node(1)
    fixedVertices=@lift(($allVertices)[$(currentIndex) > length($allVertices) ? 1 : $(currentIndex)])
    shadowPoints=@lift(vcat(($allVertices)[1:end]...))

    foreach(bar->linesegments!(ax, @lift([($fixedVertices)[Int64(bar[1])], ($fixedVertices)[Int64(bar[2])]]) ; linewidth = length(vertices[1])==2 ? 4.0 : 5.0, color=:black), bars)
    foreach(cable->linesegments!(ax, @lift([($fixedVertices)[Int64(cable[1])], ($fixedVertices)[Int64(cable[2])]]) ; color=:blue), cables)
    scatter!(ax, @lift([f for f in $shadowPoints]); color=:goldenrod3, marker = :diamond, alpha=0.1, markersize = length(vertices[1])==2 ? 12 : 75)
    scatter!(ax, @lift([f for f in $fixedVertices]); color=:red, markersize = length(vertices[1])==2 ? 12 : 75)
    scatter!(ax, cp; alpha=0.1, markersize=0.5, color=:teal)
    if(length(vertices[1])==2)
        xlimiter = Node([Inf,-Inf]); xlimiter = @lift(computeMinMax($fixedVertices, $shadowPoints, $xlimiter, 1));
        ylimiter = Node([Inf,-Inf]); ylimiter = @lift(computeMinMax($fixedVertices, $shadowPoints, $ylimiter, 2));
        @lift(limits!(ax, FRect((($xlimiter)[1]-0.5,($ylimiter)[1]-0.5), (($xlimiter)[2]-($xlimiter)[1]+1.0,($ylimiter)[2]-($ylimiter)[1]+1.0))))
    end

    # Animate the framework by giving params different values and thus triggering change in the framework.
    record(scene, "time_animation.gif", timestamps; framerate = 15) do t
        params[] = t
    end
    display(scene)
    return(fixedVertices[])
end

#= transform a list of vertices into a list of points.
=#
function to_point(listOfVertices)
    if(length(listOfVertices[1])==2)
        return([Point2f0([vx[1],vx[2]]) for vx in listOfVertices])
    else
        return([Point3f0([vx[1],vx[2],vx[3]]) for vx in listOfVertices])
    end
end

#= Given a list of lists of vertices (on a parameter path) as a list of Point2f0 or Point3f0 and a list of edges ([[i,j], ...]). fr is the preferred framerate
=#
function animateTensegrity(listOfVertexSets::Array, listOfEdges::Array, frames::Int)
    currentConfiguration = Node(listOfVertexSets[1])
    scene=Scene()
    xmin, xmax = minimum([minimum([p[1] for p in vertices]) for vertices in listOfVertexSets]), maximum([maximum([p[1] for p in vertices]) for vertices in listOfVertexSets])
    ymin, ymax = minimum([minimum([p[2] for p in vertices]) for vertices in listOfVertexSets]), maximum([maximum([p[2] for p in vertices]) for vertices in listOfVertexSets])

    scatter!(scene, @lift([f for f in to_point($currentConfiguration)]))
    xlims!(scene, (xmin-2e-1, xmax+2e-1))
    ylims!(scene, (ymin-2e-1, ymax+2e-1))
    for edge in listOfEdges
        linesegments!(scene, @lift([to_point($currentConfiguration)[edge[1]], to_point($currentConfiguration)[edge[2]]]))
    end

    timestamps = [i for i in 1:length(listOfVertexSets)]
    record(scene, "animatedTensegrity.gif", timestamps; framerate = frames) do i
        currentConfiguration[] = listOfVertexSets[i]
    end
    display(scene)
end

end


text:
"""
$(TYPEDEF)

An expression representing the application of a function to symbolic arguments.

# Fields
$(FIELDS)

# Examples

Operations can be built by application of most built-in mathematical functions
to other [`Expression`](@ref) instances:

```jldoctest
julia> using ModelingToolkit

julia> @variables x y;

julia> op1 = sin(x)
sin(x())

julia> typeof(op1.op)
typeof(sin)

julia> op1.args
1-element Array{Expression,1}:
 x()

julia> op2 = x + y
x() + y()

julia> typeof(op2.op)
typeof(+)

julia> op2.args
2-element Array{Expression,1}:
 x()
 y()
```
"""
struct Operation <: Expression
    """The function to be applied."""
    op::Function
    """The arguments the function is applied to."""
    args::Vector{Expression}
end

Base.isequal(x::Operation,y::Operation) =
    x.op == y.op && length(x.args) == length(y.args) && all(isequal.(x.args,y.args))
Base.isequal(::Operation, ::Number   ) = false
Base.isequal(::Number   , ::Operation) = false
Base.isequal(::Operation, ::Variable ) = false
Base.isequal(::Variable , ::Operation) = false
Base.isequal(::Operation, ::Constant ) = false
Base.isequal(::Constant , ::Operation) = false

Base.show(io::IO, O::Operation) = print(io, convert(Expr, O))

# For inv
Base.convert(::Type{Operation}, x::Number) = Operation(identity, Expression[Constant(x)])
Base.convert(::Type{Operation}, x::Operation) = x
Base.convert(::Type{Operation}, x::Expression) = Operation(identity, Expression[x])
Operation(x) = convert(Operation, x)

#convert to Expr
Base.Expr(op::Operation) = simplified_expr(op)
Base.convert(::Type{Expr},x::Operation) = Expr(x)

ext: jl


text:
using Test
using DiffEqBase
using Sundials

@testset "Test Reactors" begin

phaseDict = readinput("../src/testing/liquid_phase.rms")
spcs = phaseDict["phase"]["Species"]; #mechanism dictionaries index:  phaseDict[phasename]["Species" or "Reactions"]
rxns = phaseDict["phase"]["Reactions"];
solv = phaseDict["Solvents"][1];
liq = IdealDiluteSolution(spcs,rxns,solv;name="phase",diffusionlimited=true) #Define the phase (how species thermodynamic and kinetic properties calculated)

@testset "Test liquid phase Constant T Constant V reactor simulation" begin
#Constant T and V Ideal Dilute Liquid
initialconds = Dict(["T"=>450.0,"V"=>1.0e-6*1e6,"octane"=>6.154e-3*1e6,"oxygen"=>4.953e-6*1e6]) #Set simulation Initial Temp and Pressure
domain,y0,p = ConstantTVDomain(phase=liq,initialconds=initialconds,constantspecies=["oxygen"]) #Define the domain (encodes how system thermodynamic properties calculated)
react = Reactor(domain,y0,(0.0,140000.01);p=p) #Create the reactor object

sol = solve(react.ode,CVODE_BDF(),abstol=1e-20,reltol=1e-8); #solve the ode associated with the reactor

spcnames = getfield.(liq.species,:name)
octaneind = findfirst(isequal("octane"),spcnames)
y = sol(32977.61568)
@test y[octaneind]/sum(y) ≈ 0.461599061 rtol=3e-2 #from RMG simulator I believe the slight difference is due to better calculation of diffusion limits in RMS

#analytic jacobian vs. ForwardDiff jacobian
t=32977.61568;
y=sol(t)
ja=jacobiany(y,p,t,domain,[],nothing);
j=jacobianyforwarddiff(y,p,t,domain,[],nothing);
@test all((abs.(ja.-j) .> 1e-4.*abs.(j).+1e-16).==false)

jpa=jacobianp(y,p,t,domain,[],nothing);
jp=jacobianpforwarddiff(y,p,t,domain,[],nothing);
@test all((abs.(jpa.-jp) .> 1e-4.*abs.(jp).+1e-16).==false)
end;

@testset "Test liquid phase Parametrized T Constant V reactor jacobian" begin
#Parametrized T constant V Ideal Dilute Liquid
initialconds = Dict(["ts"=>[0.,600.,1200.],"T"=>[450.0,490.,500.],"V"=>1.0e-6*1e6,"octane"=>6.154e-3*1e6,"oxygen"=>4.953e-6*1e6])
domain,y0,p = ParametrizedTConstantVDomain(phase=liq,initialconds=initialconds) #Define the domain (encodes how system thermodynamic properties calculated)

react = Reactor(domain,y0,(0.0,140000.01),p=p) #Create the reactor object
sol = solve(react.ode,CVODE_BDF(),abstol=1e-20,reltol=1e-8); #solve the ode associated with the reactor

#analytic jacobian vs. ForwardDiff jacobian
t=32977.61568;
y=sol(t)
ja=jacobiany(y,p,t,domain,[],nothing);
j=jacobianyforwarddiff(y,p,t,domain,[],nothing);
@test all((abs.(ja.-j) .> 1e-4.*abs.(j).+1e-16).==false)

jpa=jacobianp(y,p,t,domain,[],nothing);
jp=jacobianpforwarddiff(y,p,t,domain,[],nothing);
@test all((abs.(jpa.-jp) .> 1e-4.*abs.(jp).+1e-16).==false)
end;

@testset "Test ConstantTAPhi Reactor Simulation" begin
    phaseDict = readinput("../src/testing/ORR.rms")
    spcs = phaseDict["phase"]["Species"]; #mechanism dictionaries index:  phaseDict[phasename]["Species" or "Reactions"]
    rxns = phaseDict["phase"]["Reactions"];
    AdivV = 1.0
    is = IdealSurface(spcs,rxns,2.486e-5;name="surf");
    initialconds = Dict(["T"=>298.0,"A"=>1.0,"Phi"=>-0.5,"A*"=>0.5,"B*"=>0.5,"OHA"=>0.0,"O2"=>0.25*AdivV,"H+"=>1.81e-6*AdivV,"H2O2aq"=>0.0*AdivV,"H2O"=>5.55e4*AdivV])
    domain,y0,p =ConstantTAPhiDomain(phase=is,initialconds=initialconds;sensitivity=false,constantspecies=["H+","O2","H2O2aq","H2O"]);
    react = Reactor(domain,y0,(0.0,1e-5),p=p);
    sol = solve(react.ode,CVODE_BDF(),abstol=1e-20,reltol=1e-6);
    sim = Simulation(sol,domain);
    @test concentrations(sim,"OA",1e-6) ≈ 0.00013773350978007822 rtol=1e-5
    @test concentrations(sim,"A*",1e-6) ≈ 1.3559599927418817e-5 rtol=1e-5
    @test concentrations(sim,"OHA",1e-6) ≈ 0.4998487068452803 rtol=1e-5
end

#Use superminimal example to test
phaseDict = readinput("../src/testing/superminimal.rms") #load mechanism dictionary
spcs = phaseDict["phase"]["Species"]; #mechanism dictionaries index:  phaseDict[phasename]["Species" or "Reactions"]
rxns = phaseDict["phase"]["Reactions"];
ig = IdealGas(spcs,rxns,name="phase")

#Constant T and P Ideal Gas
@testset "Test constant T and P reactor simulation" begin
 #Define the phase (how species thermodynamic and kinetic properties calculated)
initialconds = Dict(["T"=>1000.0,"P"=>1e5,"H2"=>0.67,"O2"=>0.33]) #Set simulation Initial Temp and Pressure
domain,y0,p = ConstantTPDomain(phase=ig,initialconds=initialconds) #Define the domain (encodes how system thermodynamic properties calculated)

react = Reactor(domain,y0,(0.0,150.11094);p=p) #Create the reactor object
sol = solve(react.ode,CVODE_BDF(),abstol=1e-20,reltol=1e-12); #solve the ode associated with the reactor
sim = Simulation(sol,domain);

spcnames = getfield.(ig.species,:name)
h2ind = findfirst(isequal("H2"),spcnames)
o2ind = findfirst(isequal("O2"),spcnames)
h2oind = findfirst(isequal("H2O"),spcnames)
y = sol(20.44002454)
N = sim.N(20.44002454)
@test y[h2ind]/N ≈ 0.412883111 rtol=1e-4 #from RMG simulator
@test y[o2ind]/N ≈ 0.200419093 rtol=1e-4
@test y[h2oind]/N ≈ 0.386618602 rtol=1e-4

#analytic jacobian vs. ForwardDiff jacobian
t=20.44002454;
y=sol(t)
ja=jacobiany(y,p,t,domain,[],nothing);
j = jacobianyforwarddiff(y,p,t,domain,[],nothing);
@test all((abs.(ja.-j) .> 1e-4.*abs.(j).+1e-16).==false)

jpa=jacobianp(y,p,t,domain,[],nothing);
jp=jacobianpforwarddiff(y,p,t,domain,[],nothing);
@test all((abs.(jpa.-jp) .> 1e-4.*abs.(jp).+1e-16).==false)

#sensitivities
dps = getadjointsensitivities(sim,"H2",CVODE_BDF(linear_solver=:GMRES);sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)),abstol=1e-16,reltol=1e-6)
react2 = Reactor(domain,y0,(0.0,150.11094);p=p,forwardsensitivities=true)
sol2 = solve(react2.ode,CVODE_BDF(linear_solver=:GMRES),abstol=1e-21,reltol=1e-7); #solve the ode associated with the reactor
sim2 = Simulation(sol2,domain)

x,dp = extract_local_sensitivities(sol2,150.11094);
ind = findfirst(isequal("H2"),sim2.names)
dpvs = [v[ind] for v in dp]
dpvs[length(domain.phase.species)+1:end] .*= domain.p[length(domain.phase.species)+1:end]
dpvs ./= sol2(150.11094)[ind]
rerr = (dpvs .- dps')./dpvs
rerr = [isinf(x) ? 0.0 : x for x in rerr]
@test all((abs.(rerr) .> 1e-1).==false)
end;

#Constant T and P Ideal Gas
@testset "Test constant T and P reactor with interfaces simulation" begin
    #Define the phase (how species thermodynamic and kinetic properties calculated)
   initialconds = Dict(["T"=>1000.0,"P"=>1e5,"H2"=>0.67,"O2"=>0.33]) #Set simulation Initial Temp and Pressure
   domain,y0,p = ConstantTPDomain(phase=ig,initialconds=initialconds) #Define the domain (encodes how system thermodynamic properties calculated)

   interfaces = [Inlet(domain,Dict{String,Float64}("H2"=>0.67,"O2"=>0.33,"T"=>1000.0,"P"=>1e5),x->0.001),
                Outlet(domain,x->0.001)]

   react = Reactor(domain,y0,(0.0,150.11094),interfaces;p=p) #Create the reactor object
   sol = solve(react.ode,CVODE_BDF(),abstol=1e-20,reltol=1e-12); #solve the ode associated with the reactor
   sim = Simulation(sol,domain,interfaces);

   #analytic jacobian vs. ForwardDiff jacobian
   t = 20.44002454;
   y = sol(t)
   ja = jacobiany(y,p,t,domain,interfaces,nothing);
   j = jacobianyforwarddiff(y,p,t,domain,interfaces,nothing);
   @test all((abs.(ja.-j) .> 1e-4.*abs.(j).+1e-16).==false)

   jpa = jacobianp(y,p,t,domain,interfaces,nothing);
   jp = jacobianpforwarddiff(y,p,t,domain,interfaces,nothing);
   @test all((abs.(jpa.-jp) .> 1e-4.*abs.(jp).+1e-16).==false)

   #sensitivities
   dps = getadjointsensitivities(sim,"H2",CVODE_BDF(linear_solver=:GMRES);sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)),abstol=1e-16,reltol=1e-6)
   react2 = Reactor(domain,y0,(0.0,150.11094),interfaces;p=p,forwardsensitivities=true)
   sol2 = solve(react2.ode,CVODE_BDF(linear_solver=:GMRES),abstol=1e-21,reltol=1e-7); #solve the ode associated with the reactor
   sim2 = Simulation(sol2,domain,interfaces)

   x,dp = extract_local_sensitivities(sol2,150.11094);
   ind = findfirst(isequal("H2"),sim2.names)
   dpvs = [v[ind] for v in dp]
   dpvs[length(domain.phase.species)+1:end] .*= domain.p[length(domain.phase.species)+1:end]
   dpvs ./= sol2(150.11094)[ind]
   rerr = (dpvs .- dps')./dpvs
   rerr = [isinf(x) ? 0.0 : x for x in rerr]
   @test all((abs.(rerr) .> 1e-1).==false)
   end;

#Constant V adiabatic Ideal Gas
#uses superminimal.yml mechanism
@testset "Constant volume adiabatic reactor simulation" begin
initialconds = Dict(["T"=>1000.0,"P"=>10.0e5,"H2"=>0.67,"O2"=>0.33]) #Set simulation Initial Temp and Pressure
domain,y0,p = ConstantVDomain(phase=ig,initialconds=initialconds) #Define the domain (encodes how system thermodynamic properties calculated)

react = Reactor(domain,y0,(0.0,0.101),p=p) #Create the reactor object
sol = solve(react.ode,CVODE_BDF(),abstol=1e-20,reltol=1e-12); #solve the ode associated with the reactor

ts = exp.(range(log(1e-15),length=10000,stop=log(0.1)))
IDT = ts[argmax(diff([sol(t)[end] for t in ts]))] #Ignition Delay Time based on argmax(dTdt(t))

@test IDT ≈ 0.038384723436228063 rtol=1e-5 #from Cantera simulation

#analytic jacobian vs. ForwardDiff jacobian
t=0.01;
y = sol(t);
ja=jacobiany(y,p,t,domain,[],nothing);
j = jacobianyforwarddiff(y,p,t,domain,[],nothing);
@test all((abs.(ja.-j) .> 1e-4.*abs.(j).+1e-16).==false)

jpa=jacobianp(y,p,t,domain,[],nothing);
jp=jacobianpforwarddiff(y,p,t,domain,[],nothing);
@test all((abs.(jpa.-jp) .> 1e-4.*abs.(jp).+1e-16).==false)

#sensitivities
react = Reactor(domain,y0,(0.0,0.02),p=p) #Create the reactor object
sol = solve(react.ode,CVODE_BDF(),abstol=1e-20,reltol=1e-12); #solve the ode associated with the reactor
sim = Simulation(sol,domain)
dps = getadjointsensitivities(sim,"H2",CVODE_BDF();sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(false)),abstol=1e-16,reltol=1e-6)
react2 = Reactor(domain,y0,(0.0,0.02);p=p,forwardsensitivities=true)
sol2 = solve(react2.ode,CVODE_BDF(),abstol=1e-16,reltol=1e-6); #solve the ode associated with the reactor
sim2 = Simulation(sol2,domain)

x,dp = extract_local_sensitivities(sol2,0.02);
ind = findfirst(isequal("H2"),sim2.names)
dpvs = [v[ind] for v in dp]
dpvs[length(domain.phase.species)+1:end] .*= domain.p[length(domain.phase.species)+1:end]
dpvs ./= sol2(0.02)[ind]
rerr = (dpvs .- dps')./dpvs
rerr = [isinf(x) ? 0.0 : x for x in rerr]
@test all((abs.(rerr) .> 2e-1).==false)
end;

#Constant P adiabatic Ideal Gas
#uses ethane.rms mechanism
@testset "Constant pressure adiabatic reactor simulation" begin

phaseDict = readinput("../src/testing/ethane.rms")
spcs = phaseDict["phase"]["Species"]
rxns = phaseDict["phase"]["Reactions"]
ig = IdealGas(spcs,rxns,name="phase")

initialconds = Dict(["T"=>1000.0,"P"=>2.0e5,"ethane"=>1.0,"Ar"=>1.0,"O2"=>3.5]) #Set simulation Initial Temp and Pressure
domain,y0,p = ConstantPDomain(phase=ig,initialconds=initialconds) #Define the domain (encodes how system thermodynamic properties calculated)

react = Reactor(domain,y0,(0.0,1.0);p=p) #Create the reactor object
sol = solve(react.ode,CVODE_BDF(),abstol=1e-16,reltol=1e-6); #solve the ode associated with the reactor

ts = exp.(range(log(1e-15),length=10000,stop=log(0.2)))
IDT = ts[argmax(diff([sol(t)[end] for t in ts]))] #Ignition Delay Time based on argmax(dTdt(t))

@test IDT ≈ 0.07324954954380769 rtol=1e-5

#analytic jacobian vs. ForwardDiff jacobian
t=0.01;
y=sol(t)
ja=jacobiany(y,p,t,domain,[],nothing);
j = jacobianyforwarddiff(y,p,t,domain,[],nothing);
@test all((abs.(ja.-j) .> 1e-4.*abs.(j).+1e-16).==false)

jpa=jacobianp(y,p,t,domain,[],nothing);
jp=jacobianpforwarddiff(y,p,t,domain,[],nothing);
@test all((abs.(jpa.-jp) .> 1e-4.*abs.(jp).+1e-16).==false)

end;

#Parametrized T and P Ideal Gas
#uses superminimal.yml mechanism
@testset "Parametrized T and P reactor jacobian" begin
initialconds = Dict(["T"=>[1000.0,1400.0,1500.0],"ts"=>[0.,100.,200.],"P"=>[1.0e5,1.8e5,2.0e5],"H2"=>0.67,"O2"=>0.33]) #Set simulation Initial Temp/Pressure and Volume (function/array)
domain,y0,p = ParametrizedTPDomain(phase=ig,initialconds=initialconds) #Define the domain (encodes how system thermodynamic properties calculated)

react = Reactor(domain,y0,(0.0,0.101),p=p) #Create the reactor object
sol = solve(react.ode,CVODE_BDF(),abstol=1e-20,reltol=1e-12); #solve the ode associated with the reactor

#analytic jacobian vs. ForwardDiff jacobian
t=0.1;
y = sol(t);
ja=jacobiany(y,p,t,domain,[],nothing);
j = jacobianyforwarddiff(y,p,t,domain,[],nothing);
@test all((abs.(ja.-j) .> 1e-4.*abs.(j).+1e-16).==false)

jpa=jacobianp(y,p,t,domain,[],nothing);
jp=jacobianpforwarddiff(y,p,t,domain,[],nothing);
@test all((abs.(jpa.-jp) .> 1e-4.*abs.(jp).+1e-16).==false)
end;

#Parametrized V adiabatic Ideal Gas
#uses superminimal.yml mechanism
@testset "Parametrized volume adiabatic reactor jacobian" begin
initialconds = Dict(["T"=>1000.0,"ts"=>[0.,100.,200.],"V"=>[0.01,0.05,0.1],"H2"=>0.67,"O2"=>0.33]) #Set simulation Initial Temp/Pressure and Volume (function/array)
domain,y0,p = ParametrizedVDomain(phase=ig,initialconds=initialconds) #Define the domain (encodes how system thermodynamic properties calculated)

react = Reactor(domain,y0,(0.0,0.101),p=p) #Create the reactor object
sol = solve(react.ode,CVODE_BDF(),abstol=1e-20,reltol=1e-12); #solve the ode associated with the reactor

#analytic jacobian vs. ForwardDiff jacobian
t=0.1;
y = sol(t);
ja=jacobiany(y,p,t,domain,[],nothing);
j = jacobianyforwarddiff(y,p,t,domain,[],nothing);
@test all((abs.(ja.-j) .> 1e-4.*abs.(j).+1e-16).==false)

jpa=jacobianp(y,p,t,domain,[],nothing);
jp=jacobianpforwarddiff(y,p,t,domain,[],nothing);
@test all((abs.(jpa.-jp) .> 1e-4.*abs.(jp).+1e-16).==false)
end;

#Parametrized P adiabatic Ideal Gas
#uses ethane.yml mechanism
@testset "Parametrized pressure adiabatic reactor jacobian" begin
initialconds = Dict(["T"=>1000.0,"ts"=>[0.,100.,200.],"P"=>[2.0e5,3.0e5,5.0e5],"H2"=>0.67,"O2"=>0.33])
domain,y0,p = ParametrizedPDomain(phase=ig,initialconds=initialconds) #Define the domain (encodes how system thermodynamic properties calculated)

react = Reactor(domain,y0,(0.0,0.101),p=p) #Create the reactor object
sol = solve(react.ode,CVODE_BDF(),abstol=1e-20,reltol=1e-12); #solve the ode associated with the reactor

#analytic jacobian vs. ForwardDiff jacobian
t=0.1;
y = sol(t);
ja=jacobiany(y,p,t,domain,[],nothing);
j = jacobianyforwarddiff(y,p,t,domain,[],nothing);
@test all((abs.(ja.-j) .> 1e-4.*abs.(j).+1e-16).==false)

jpa=jacobianp(y,p,t,domain,[],nothing);
jp=jacobianpforwarddiff(y,p,t,domain,[],nothing);
@test all((abs.(jpa.-jp) .> 1e-4.*abs.(jp).+1e-16).==false)
end;

@testset "Multi-domain ConstantV sensitivity analysis" begin
    phaseDict = readinput("../src/testing/superminimal.rms")
    spcs = phaseDict["phase"]["Species"]
    rxns = phaseDict["phase"]["Reactions"]
    ig = IdealGas(spcs,rxns,name="phase")

    initialcondsV = Dict(["T"=>1000.0,"P"=>10.0e5,"H2"=>0.67,"O2"=>0.33])
    domainV,y0V,pV = ConstantVDomain(phase=ig,initialconds=initialcondsV) #Define the domain (encodes how system thermodynamic properties calculated)

    reactV = Reactor(domainV,y0V,(0.0,0.037);p=pV) #Create the reactor object
    solV = solve(reactV.ode,CVODE_BDF(),abstol=1e-16,reltol=1e-6); #solve the ode associated with the reactor
    simV = Simulation(solV,domainV)

    initialcondsV1 = Dict(["T"=>1000.0,"P"=>10.0e5,"H2"=>0.67,"O2"=>0.33])
    domainV1,y0V1,pV1 = ConstantVDomain(phase=ig,initialconds=initialcondsV1) #Define the domain (encodes how system thermodynamic properties calculated)
    initialcondsV2 = Dict(["T"=>1000.0,"P"=>10.0e5,"H2"=>0.67,"O2"=>0.33])
    domainV2,y0V2,pV2 = ConstantVDomain(phase=ig,initialconds=initialcondsV2) #Define the domain (encodes how system thermodynamic properties calculated)

    react,y0,p = Reactor((domainV1,domainV2),(y0V1,y0V2),(0.0,0.037),[],(pV1,pV2));
    sol = solve(react.ode,CVODE_BDF(),abstol=1e-16,reltol=1e-6);
    sysim = SystemSimulation(sol,(domainV1,domainV2),[],p);

    t = 0.03
    @test sol(t)[1:length(spcs)] ≈ solV(t)[1:end-2] rtol=1e-5
    @test sol(t)[length(spcs)+1:end-4] ≈ solV(t)[1:end-2] rtol=1e-5

    dpsV = getadjointsensitivities(simV,"H2",CVODE_BDF();sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(false)),abstol=1e-16,reltol=1e-6)
    dps = getadjointsensitivities(sysim,sysim.sims[1],"H2",CVODE_BDF();sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(false)),abstol=1e-16,reltol=1e-6)
    @test dpsV ≈ dps rtol=1e-4
end;

@testset "Multi-domain ConstantV and ConstantTP simulation" begin
    phaseDict = readinput("../src/testing/superminimal.rms")
    spcs = phaseDict["phase"]["Species"]
    rxns = phaseDict["phase"]["Reactions"]
    ig = IdealGas(spcs,rxns,name="phase")

    initialcondsTP = Dict(["T"=>1000.0,"P"=>10.0e5,"H2"=>0.67,"O2"=>0.33])
    domainTP,y0TP,pTP = ConstantTPDomain(phase=ig,initialconds=initialcondsTP) #Define the domain (encodes how system thermodynamic properties calculated)

    reactTP = Reactor(domainTP,y0TP,(0.0,0.04);p=pTP) #Create the reactor object
    solTP = solve(reactTP.ode,CVODE_BDF(),abstol=1e-16,reltol=1e-6); #solve the ode associated with the reactor

    initialcondsV = Dict(["T"=>1000.0,"P"=>10.0e5,"H2"=>0.67,"O2"=>0.33])
    domainV,y0V,pV = ConstantVDomain(phase=ig,initialconds=initialcondsV) #Define the domain (encodes how system thermodynamic properties calculated)

    reactV = Reactor(domainV,y0V,(0.0,0.04);p=pV) #Create the reactor object
    solV = solve(reactV.ode,CVODE_BDF(),abstol=1e-16,reltol=1e-6); #solve the ode associated with the reactor

    initialcondsTP = Dict(["T"=>1000.0,"P"=>10.0e5,"H2"=>0.67,"O2"=>0.33])
    domainTP,y0TP,pTP = ConstantTPDomain(phase=ig,initialconds=initialcondsTP) #Define the domain (encodes how system thermodynamic properties calculated)
    initialcondsV = Dict(["T"=>1000.0,"P"=>10.0e5,"H2"=>0.67,"O2"=>0.33])
    domainV,y0V,pV = ConstantVDomain(phase=ig,initialconds=initialcondsV) #Define the domain (encodes how system thermodynamic properties calculated)

    react,y0,p = Reactor((domainTP,domainV),(y0TP,y0V),(0.0,0.04),[],(pTP,pV));
    sol = solve(react.ode,CVODE_BDF(),abstol=1e-16,reltol=1e-6);

    t = 0.03
    @test sol(t)[1:length(spcs)] ≈ solTP(t)[1:end-1] rtol=1e-5
    print(length(sol(t)[length(spcs)+1:end-3]))
    print(length(solV(t)[1:end-2]))
    @test sol(t)[length(spcs)+1:end-3] ≈ solV(t)[1:end-2] rtol=1e-5
end;

@testset "Multi-domain Gas-Surface ConstantTP and ConstantTAPhi Simulation" begin
    phaseDict = readinput("../src/testing/ch4o2cat.rms")
    gasspcs = phaseDict["gas"]["Species"];
    gasrxns = phaseDict["gas"]["Reactions"];
    surfacespcs = phaseDict["surface"]["Species"]
    surfacerxns = phaseDict["surface"]["Reactions"]
    interfacerxns = phaseDict[Set(["gas","surface"])]["Reactions"];

    ig = IdealGas(gasspcs,gasrxns;name="gas");
    cat = IdealSurface(surfacespcs,surfacerxns,2.486e-5;name="surface");

    initialconds = Dict(["T"=>800.0,"P"=>1.0e5,"O2"=>0.2,"N2"=>0.7,"CH4"=>0.1]);
    domaingas,y0gas,pgas = ConstantTPDomain(phase=ig,initialconds=initialconds,);

    V = 8.314*800.0/1.0e5
    A = 1.0e5*V
    initialconds = Dict(["T"=>800.0,"A"=>A,"vacantX"=>cat.sitedensity*A]);
    domaincat,y0cat,pcat = ConstantTAPhiDomain(phase=cat,initialconds=initialconds,);

    inter,pinter = ReactiveInternalInterfaceConstantTPhi(domaingas,domaincat,interfacerxns,800.0,A);

    react,y0,p = Reactor((domaingas,domaincat),(y0gas,y0cat),(0.0,0.1),(inter,),(pgas,pcat,pinter));

    sol = solve(react.ode,CVODE_BDF(),abstol=1e-20,reltol=1e-6);

    ssys = SystemSimulation(sol,(domaingas,domaincat,),(inter,),p);

    @test concentrations(ssys,"OX",0.5e-5) ≈ 8.033191655902819e-6 rtol=1e-5
    @test molefractions(ssys.sims[1],"H2O",0.5e-5) ≈ 0.10899527627867926 rtol=1e-5
end;

@testset "Multi-domain Gas-Surface ConstantV and ConstantTAPhi Simulation" begin
    phaseDict = readinput("../src/testing/ch4o2cat.rms")
    gasspcs = phaseDict["gas"]["Species"];
    gasrxns = phaseDict["gas"]["Reactions"];
    surfacespcs = phaseDict["surface"]["Species"]
    surfacerxns = phaseDict["surface"]["Reactions"]
    interfacerxns = phaseDict[Set(["gas","surface"])]["Reactions"];

    ig = IdealGas(gasspcs,gasrxns;name="gas");
    cat = IdealSurface(surfacespcs,surfacerxns,2.486e-5;name="surface");

    initialconds = Dict(["T"=>800.0,"P"=>1.0e5,"O2"=>0.2,"N2"=>0.7,"CH4"=>0.1]);
    domaingas,y0gas,pgas = ConstantVDomain(phase=ig,initialconds=initialconds,);

    V = 8.314*800.0/1.0e5
    A = 1.0e5*V
    initialconds = Dict(["T"=>800.0,"A"=>A,"vacantX"=>cat.sitedensity*A]);
    domaincat,y0cat,pcat = ConstantTAPhiDomain(phase=cat,initialconds=initialconds,);

    inter,pinter = ReactiveInternalInterface(domaingas,domaincat,interfacerxns,A);

    react,y0,p = Reactor((domaingas,domaincat),(y0gas,y0cat),(0.0,0.1),(inter,),(pgas,pcat,pinter));

    sol = solve(react.ode,CVODE_BDF(),abstol=1e-20,reltol=1e-6);

    ssys = SystemSimulation(sol,(domaingas,domaincat,),(inter,),p);

    @test concentrations(ssys,"OX",0.5e-5) ≈ 1.9165723392283484e-5 rtol=1e-5
    @test molefractions(ssys.sims[1],"H2O",1e-3) ≈ 0.12732345278036702 rtol=1e-5
end;
end;


text:
command_pool = CommandPool(device, 0)
buffer = first(unwrap(allocate_command_buffers(device, CommandBufferAllocateInfo(command_pool, VK_COMMAND_BUFFER_LEVEL_PRIMARY, 1))))


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

function apply_smoothing!(mesh::UnstructuredMesh2D, alpha, alpha_tmp, dg, cache)
  # Diffuse alpha values by setting each alpha to at least 50% of neighboring elements' alpha
  # Copy alpha values such that smoothing is indpedenent of the element access order
  alpha_tmp .= alpha

  # Loop over interfaces
  for interface in eachinterface(dg, cache)
    # Get neighboring element ids
    left  = cache.interfaces.element_ids[1, interface]
    right = cache.interfaces.element_ids[2, interface]

    # Apply smoothing
    alpha[left]  = max(alpha_tmp[left],  0.5 * alpha_tmp[right], alpha[left])
    alpha[right] = max(alpha_tmp[right], 0.5 * alpha_tmp[left],  alpha[right])
  end
end

end # @muladd


text:
function con(q)
    # joint constraint function
    # c(q_0) should be all zeros
    c_ = zeros(eltype(q), 23)

    rb = q[1:3]
    Qb = q[4:7]
    r0 = q[8:10]
    Q0 = q[11:14]
    r1 = q[15:17]
    Q1 = q[18:21]
    r2 = q[22:24]
    Q2 = q[25:28]
    r3 = q[29:31]
    Q3 = q[32:35]

    pb = rb + rotate(Qb, l_cb) # position vector from world frame to *JOINTS* 0 and 2

    c_[1:3] = pb - r0 - rotate(Q0, -l_c0)
    c_[4:5] = [0 1 0 0; 0 0 0 1]*L(Qb)'*Q0  # y axis rotation constraint: set x & z = 0
    c_[6:8] = r0 + rotate(Q0,  l0-l_c0) - r1 - rotate(Q1, -l_c1)
    c_[9:10] = [0 1 0 0; 0 0 0 1]*L(Q0)'*Q1
    c_[11:13] = pb - r2 - rotate(Q2, -l_c2)
    c_[14:15] = [0 1 0 0; 0 0 0 1]*L(Qb)'*Q2
    c_[16:18] = r2 + rotate(Q2, l2 - l_c2) - r3 - rotate(Q3, -l_c3)
    c_[19:20] = [0 1 0 0; 0 0 0 1]*L(Q2)'*Q3
    c_[21:23] = r1 + rotate(Q1, l1 - l_c1) - r3 - rotate(Q3, lc-l_c3)
    return c_
end

function Dc(q)
    ForwardDiff.jacobian(dq->con(dq),q)*Ḡ(q)
end

function ϕ(q)
    ϕ = zeros(eltype(q), 1)
    # signed distance
    r3 = q[29:31]
    Q3 = q[32:35]
    rf = r3 + rotate(Q3, lee-l_c3) # position of foot
    ϕ[1] = rf[3] - 0.025  # subtract radius of foot
    return ϕ
end

function Dϕ(q)
    # dϕ = reshape(ForwardDiff.jacobian(dq->ϕ(dq),q), 1, 35)*Ḡ(q)  # 30x1
    dϕ = ForwardDiff.jacobian(dq->ϕ(dq),q)*Ḡ(q)  # 30x1
    return dϕ
end

function del(m, I, r1, r2, r3, Q1, Q2, Q3, grav)
    [(1/h)*m*(r2-r1) - (1/h)*m*(r3-r2) - m*grav*h;
     (2.0/h)*G(Q2)'*L(Q1)*H*I*H'*L(Q1)'*Q2 + (2.0/h)*G(Q2)'*T*R(Q3)'*H*I*H'*L(Q2)'*Q3]
end

function DEL(q_1,q_2,q_3,λ,F1,F2,grav,b,n)

    rb_1 = q_1[1:3]
    Qb_1 = q_1[4:7]
    r0_1 = q_1[8:10]
    Q0_1 = q_1[11:14]
    r1_1 = q_1[15:17]
    Q1_1 = q_1[18:21]
    r2_1 = q_1[22:24]
    Q2_1 = q_1[25:28]
    r3_1 = q_1[29:31]
    Q3_1 = q_1[32:35]

    rb_2 = q_2[1:3]
    Qb_2 = q_2[4:7]
    r0_2 = q_2[8:10]
    Q0_2 = q_2[11:14]
    r1_2 = q_2[15:17]
    Q1_2 = q_2[18:21]
    r2_2 = q_2[22:24]
    Q2_2 = q_2[25:28]
    r3_2 = q_2[29:31]
    Q3_2 = q_2[32:35]

    rb_3 = q_3[1:3]
    Qb_3 = q_3[4:7]
    r0_3 = q_3[8:10]
    Q0_3 = q_3[11:14]
    r1_3 = q_3[15:17]
    Q1_3 = q_3[18:21]
    r2_3 = q_3[22:24]
    Q2_3 = q_3[25:28]
    r3_3 = q_3[29:31]
    Q3_3 = q_3[32:35]

    del1 = [del(mb, Ib, rb_1, rb_2, rb_3, Qb_1, Qb_2, Qb_3, grav);
            del(m0, I0, r0_1, r0_2, r0_3, Q0_1, Q0_2, Q0_3, grav);
            del(m1, I1, r1_1, r1_2, r1_3, Q1_1, Q1_2, Q1_3, grav);
            del(m2, I2, r2_1, r2_2, r2_3, Q2_1, Q2_2, Q2_3, grav);
            del(m3, I3, r3_1, r3_2, r3_3, Q3_1, Q3_2, Q3_3, grav)]

    Fb = zeros(eltype(b), 30)
    Fb[25:26] = b  # should be last link's x and y translation terms

    return del1 + (h/2.0)*F1 + (h/2.0)*F2 +
        reshape(h*Dc(q_2)'*λ, 30) + reshape(h*Dϕ(q_2)'*n, 30) + h*Fb
end

#Objective and constraint functions for IPOPT

function objective(z)
    s = z[n_q+n_c+1:n_q+n_c+n_s]
    λj = z[n_q+1:n_q+n_c]
    α = 1e-3
    # regularizer (Tikhonov)
    return α*λj'*λj + sum(s) #Minimize slacks associated with complementarity conditions
    # return sum(s)
end

function constraint!(c,z)
    qn = z[1:n_q]
    λ = z[n_q+1:n_q+n_c]
    s = z[n_q+n_c+1:n_q+n_c+n_s]
    b = z[n_q+n_c+n_s+1:n_q+n_c+n_s+n_b]
    λf = z[n_q+n_c+n_s+n_b+1:n_q+n_c+n_s+n_b+n_λf]
    n = z[n_q+n_c+n_s+n_b+n_λf+1:end]
    # nonlinear DEL equation                (30 equality constraint)   c1
    # quaternion norm squared - 1 = 0       (5 equality constraints)   c2-c6
    # joint constraints                     (23 equality constraints)  c7
    # maximum dissipation (friction)        (2 equality constraints)   c8

    # signed distance of foot from ground   (1 inequality constraint)  c9
    # collision complementarity foot        (1 inequality constraint)  c10
    # friction cone                         (1 inequality constraint)  c11
    # friction relaxed complementarity      (1 inequality constraint)  c12

    μ = 0.5  # friction coefficient
    vm = (qn-qhist[:,k])/h

    # equality constraints
    c1 = DEL(qhist[:,k-1], qhist[:,k], qn, λ, F_prev, F, ghist[:, k], b, n)  # 30x1
    c2 = norm(qn[4:7])^2 - 1  # 1x1
    c3 = norm(qn[11:14])^2 - 1
    c4 = norm(qn[18:21])^2 - 1
    c5 = norm(qn[25:28])^2 - 1
    c6 = norm(qn[32:35])^2 - 1
    c7 = con(qn)
    # c8 = Jac(qn)[1:2, :]*vm + λf.*b/(smoothsqrt(b'*b)+1e-7)  # maximum dissipation
    c8  = [0.; 0.]
    # inequality constraints
    c9 = ϕ(qn)  # signed distance
    c10 = s[1] .- n.*ϕ(qn)  # relaxed complementarity (signed dist) 1x1
    c11 = μ*n .- smoothsqrt(b'*b) # friction cone
    c12 = s[2] .- λf.*(μ*n .- smoothsqrt(b'*b)) # relaxed complementarity (friction)
    c .= [c1; c2; c3; c4; c5; c6; c7; c8; c9; c10; c11; c12]

    return nothing
end

function primal_bounds(n)
    #Enforce simple bound constraints on the decision variables (e.g. positivity) here
    # x_l ≤ [q; λ; s; b; λf] ≤ x_u
    x_l = -Inf*ones(n)  # 60
    x_l[n_q+n_c+1:n_q+n_c+n_s] = zeros(n_s)  # slack variables
    x_l[n_q+n_c+n_s+n_b+1:n_q+n_c+n_s+n_b+n_λf] = zeros(n_λf)  # λf is an unsigned magnitude?
    x_l[n_q+n_c+n_s+n_b+n_λf+1:end] = zeros(n_n)  # normal force cannot be into ground

    x_u = Inf*ones(n)

    # x_l[1:2] = zeros(2)
    # x_u[1:2] = zeros(2)
    return x_l, x_u
end

function constraint_check(z, n_tol)
    qn = z[1:n_q]
    λ = z[n_q+1:n_q+n_c]
    s = z[n_q+n_c+1:n_q+n_c+n_s]
    b = z[n_q+n_c+n_s+1:n_q+n_c+n_s+n_b]
    λf = z[n_q+n_c+n_s+n_b+1:n_q+n_c+n_s+n_b+n_λf]
    n = z[n_q+n_c+n_s+n_b+n_λf+1:end]

    μ = 0.5  # friction coefficient
    vm = (qn-qhist[:,k])/h

    # equality constraints
    c1 = DEL(qhist[:,k-1], qhist[:,k], qn, λ, F_prev, F, ghist[:, k], b, n)  # 30x1
    c2 = norm(qn[4:7])^2 - 1  # 1x1
    c3 = norm(qn[11:14])^2 - 1
    c4 = norm(qn[18:21])^2 - 1
    c5 = norm(qn[25:28])^2 - 1
    c6 = norm(qn[32:35])^2 - 1
    c7 = con(qn)
    c8 = [0.0; 0.0]  # J(qn)[1:2, :]*vm + λf.*b/smoothsqrt(b'*b)  # maximum dissipation

    # inequality constraints
    c9 = ϕ(qn)  # signed distance
    c10 = s[1] .- n.*ϕ(qn)  # relaxed complementarity (signed dist) 1x1
    c11 = μ*n .- smoothsqrt(b'*b) # friction cone
    c12 = s[2] .- λf.*(μ*n .- smoothsqrt(b'*b)) # relaxed complementarity (friction)

    A = [c1; c2; c3; c4; c5; c6; c7; c8]
    B = [c9; c10; c11; c12]

    if !isapprox(A, zeros(size(A)[1]); atol=n_tol, rtol=0)  # 58
        e = 1
        #print("\n", A, "\n")
        print(findall(A .< -ones(size(A)[1])*n_tol), " is less than 0 \n")
        print(findall(A .> ones(size(A)[1])*n_tol), " is greater than 0 \n")
        print("\n Sim stopped due to ipopt infeasibility \n")
    elseif B < -ones(size(B)[1])*n_tol  #25
        e = 1
        #print("\n", B, "\n")
        print(findall(B .< -ones(size(B)[1])*n_tol))  # 25
        print("\n Sim stopped due to ipopt infeasibility \n")
    elseif (-angle_y(Q0, Q1)) < 18*pi/180  # constrain relative angle between links 0 and 1
        e = 1
        print("RoM lower limit surpassed")
    elseif 166*pi/180 < (-angle_y(Q0, Q1))
        e = 1
        print("RoM upper limit surpassed")
    else
        e = 0
    end

    return e
end

max_stars_repo_name: bbokser/hopper-sim

max_issues_repo_name: bbokser/hopper-sim

max_forks_repo_name: bbokser/hopper-sim


text:
using Test
using Random
using QuantileSummaries

@testset "String" begin
    # string
    qb = QuantileStat{String}()

    for i = 1:1000
        fit!(qb, randstring(10))
    end

    summary = value(qb)
    @test true #passed without erroring
end

@testset "Correctness" begin

    eps = 0.01
    qb = QuantileStat(eps)

    n = 10_000
    r = randn(n)
    fit!(qb, r)
    sort!(r)

    qs = summarize(qb)

    # Test if the right percentile is identified (within eps)
    for percentile = eps:eps/2:1

        valueforpercentile = r[ Int(floor(percentile * n)) ]
        predictedpercentile = QuantileSummaries.percentile(qs, valueforpercentile)
        @test abs(predictedpercentile-percentile) < eps

    end

    for i = 1:1000
        percentile = rand()
        idx = 1 + Int(floor(percentile * n))
        valueforpercentile = r[idx]
        predictedpercentile = QuantileSummaries.percentile(qs, valueforpercentile)
        @test abs(predictedpercentile-percentile) < eps
    end

    # check if computed value is within eps of expected
    # Percentiles & values are the same in this setup
    eps = 0.01
    qb = QuantileStat(eps)

    r = collect(0:0.0001:1)
    shuffle!(r)
    fit!(qb, r)
    qs = summarize(qb)

    for p in 0:0.001:1
        value = qvalue(qs, p)
        @test( abs(p-value) < eps )
    end
end

@testset "merge!" begin
    eps = 0.01
    qb1 = QuantileStat(eps)
    qb2 = QuantileStat(eps)

    n = 10_000
    r = randn(n)
    for value in r
        if rand() < 0.1
            fit!(qb1, value)
        else
            fit!(qb2, value)
        end
    end

    merge!(qb1, qb2)
    sort!(r)

    qs = summarize(qb1)
    for percentile = eps:eps/2:1
        valueforpercentile = r[ Int(floor(percentile * n)) ]
        predictedpercentile = QuantileSummaries.percentile(qs, valueforpercentile)
        @test abs(predictedpercentile-percentile) < eps
    end

end


text:
using PaddedMatrices, LinearAlgebra, Aqua
using Test

import InteractiveUtils
InteractiveUtils.versioninfo(stdout; verbose = true)

@show Threads.nthreads(), PaddedMatrices._nthreads() PaddedMatrices.VectorizationBase.REGISTER_COUNT
const START_TIME = time()

@inferred PaddedMatrices.matmul_params(Float32)
@inferred PaddedMatrices.matmul_params(Float64)
@inferred PaddedMatrices.matmul_params(Int16)
@inferred PaddedMatrices.matmul_params(Int32)
@inferred PaddedMatrices.matmul_params(Int64)

@time @testset "PaddedMatrices.jl" begin
    Aqua.test_all(PaddedMatrices)
    # @test isempty(detect_unbound_args(PaddedMatrices))
    @time include("misc.jl")
    @time include("matmul_tests.jl")
    @time include("broadcast_tests.jl")
end

const ELAPSED_MINUTES = (time() - START_TIME)/60
@test ELAPSED_MINUTES < 180


text:
module NeuroMetadata

using NeuroCore.AnatomicalAPI
using FieldProperties

export
    EncodingDirection,
    encoding_names,
    freqdim,
    freqdim!,
    phasedim,
    phasedim!,
    slice_start,
    slice_start!,
    slice_end,
    slice_end!,
    slicedim,
    slicedim!,
    slice_duration,
    slice_duration!,
    phase_encoding_direction,
    phase_encoding_direction!,
    slice_encoding_direction,
    slice_encoding_direction!,
    InstitutionInformation,
    HardwareMetadata,
    EncodingDirectionMetadata

"""
    ContrastIngredient

An enumerable type with the following possible values:
* `IODINE`
* `GADOLINIUM`
* `CARBON`
* `DIOXIDE`
* `BARIUM`
* `XENON`
* `UnkownContrast`
"""
@enum ContrastIngredient begin
    IODINE
    GADOLINIUM
    CARBON
    DIOXIDE
    BARIUM
    XENON
    UnkownContrast
end

ContrastIngredient(i::AbstractString) = ContrastIngredient(Symbol(i))

function ContrastIngredient(i::Symbol)
    if i === :IODINE
        return IODINE
    elseif i === :GADOLINIUM
        return GADOLINIUM
    elseif i === :CARBON
        return CARBON
    elseif i === :DIOXIDE
        return DIOXIDE
    elseif i === :BARIUM
        return BARIUM
    elseif i === :XENON
        return XENON
    else
        return UnkownContrast
    end
end

Base.String(i::ContrastIngredient) = String(Symbol(i))

include("mri.jl")
include("encoding_direction.jl")
include("magnetization_transfer.jl")
include("sequence.jl")
include("spatial_encoding.jl")
include("spoiling.jl")
include("time.jl")
include("electrophysiology.jl")
include("institution.jl")
include("hardware.jl")
include("task.jl")

end


text:
function plot_3d_view( x_3d::Array{Float64,2} ; xf = [] )
    if(isempty(xf)); xf = ones(size(x_3d)); end
    gr()
    p1 = scatter( x_3d[1,:] , x_3d[2 ,:] , x_3d[3,:] , zcolor = xf , ma = 0.2 , m=3, msw=0 , c=:lightrainbow , leg=:none)
    p2 = scatter( x_3d[1,:] , x_3d[2 ,:]  , zcolor = xf , ma = 0.2 , m=3, msw=0 , c=:lightrainbow , leg=:none)
    p3 = scatter( x_3d[1,:] , x_3d[3 ,:]  , zcolor = xf , ma = 0.2 , m=3, msw=0 , c=:lightrainbow , leg=:none )
    p4 = scatter( x_3d[2,:] , x_3d[3 ,:]  , zcolor = xf , ma = 0.2 , m=3, msw=0 , c=:lightrainbow , leg=:none )

    display(plot(p1,p2,p3,p4,layout=(2,2)))
end

# dist-to-spiral function
spiral_3d = reduce( hcat , [ [sin(0.5*z);cos(0.5*z);z] for z in linspace(-20.,20.,60)] )

li_spir   = SamplingHS3.create_line(spiral_3d)

#f_cost_spir_3d = (x::Vector{Float64}) -> maximum( [0.5 ; min( SamplingHS3.dist2line(li_spir,x[:,:] + [2.5;2.5;0] )[1] , SamplingHS3.dist2line(li_spir,x[:,:] - [2.5;2.5;0] )[1] ) ] )
f_cost_spir_3d = (x::Vector{Float64}) -> begin sdi = SamplingHS3.dist2line(li_spir,x[:,:])[1]; (sdi<0.5)? 0.5-(0.5-sdi)^4 : sdi end

# polytope
xlim = [-20.,20.]; ylim = [-20.,20.]; zlim = [-20.,20.];
bb_G =[eye(3);-eye(3)]
bb_h =[xlim[2];ylim[2];zlim[2];-xlim[1];-ylim[1];-zlim[1]]

using Plots
Plots.gr()
Plots.plot(spiral_3d[1,:],spiral_3d[2,:],spiral_3d[3,:])

X_init = SamplingHS3.unifsamp_box( [-20 20;-20 20;-20 20] , 400 )
ns_result_01 = SamplingHS3.nested_sampling( X_init[:,:] , f_cost_spir_3d , bb_G , bb_h , SamplingHS3.NSConfig(6,0.75,SamplingHS3.StopCriterion(max_time=30.0),SamplingHS3.MCMCConfig()) )

## Plot the animated sampling process (just because it's cool)
x_live_1 = SamplingHS3.reconstruct_live_samples(ns_result_01.ns_iter_data)
x_live  = x_live_1[1]
xf_live = x_live_1[2]
for zi in 1:length(xf_live)
    #xf_zi = 1:length(xf_live) )
    xf_zi = xf_live[zi]
    display( plot_3d_view( x_live[zi] , xf = xf_zi ) )
end


text:
using POLO: loss!, boost!, stepsize, smooth!, prox!, terminate, log
using LinearAlgebra
import POLO: getf, getx!, initialize!
import Base: copyto!

abstract type MultiThread <: ExecutionPolicy end

initialize!(::ProxGradient{<:MultiThread}) = nothing

mutable struct Consistent <: MultiThread
    mutex::Threads.Mutex
    k::Int
    fval::Float64
    x::Vector{Float64}
    g::Vector{Float64}

    function (::Type{Consistent})()
        return new(Threads.Mutex(),1,0.,Vector{Float64}(),Vector{Float64}())
    end
end

getk(consistent::Consistent) = consistent.k
getf(consistent::Consistent) = consistent.fval
setf!(consistent::Consistent, fval::AbstractFloat) = consistent.fval = fval
readx(consistent::Consistent) = consistent.x
readg(consistent::Consistent) = consistent.g
increment(consistent::Consistent) = consistent.k += 1

function initialize!(consistent::Consistent, x₀::AbstractVector)
    consistent.k = 1
    consistent.fval = 0.
    resize!(consistent.x,length(x₀))
    consistent.x .= x₀
    resize!(consistent.g,length(x₀))
    consistent.g .= zeros(length(x₀))
end

mutable struct Inconsistent <: MultiThread
    k::Threads.Atomic{Int}
    fval::Threads.Atomic{Float64}
    x::Vector{Threads.Atomic{Float64}}
    g::Vector{Threads.Atomic{Float64}}

    function (::Type{Inconsistent})()
        return new(Threads.Atomic{Int}(1),Threads.Atomic{Float64}(0.),Vector{Threads.Atomic{Float64}}(),Vector{Threads.Atomic{Float64}}())
    end
end

function copyto!(::IndexStyle, dest::AbstractArray{<:Threads.Atomic}, ::IndexStyle, src::AbstractArray)
    destinds, srcinds = LinearIndices(dest), LinearIndices(src)
    isempty(srcinds) || (checkbounds(Bool, destinds, first(srcinds)) && checkbounds(Bool, destinds, last(srcinds))) ||
        throw(BoundsError(dest, srcinds))
    @inbounds for i in srcinds
        dest[i][] = src[i]
    end
    return dest
end

getk(inconsistent::Inconsistent) = inconsistent.k[]
getf(inconsistent::Inconsistent) = inconsistent.fval[]
setf!(inconsistent::Inconsistent, fval::AbstractFloat) = inconsistent.fval[] = fval
readx(inconsistent::Inconsistent) = getindex.(inconsistent.x)
readg(inconsistent::Inconsistent) = getindex.(inconsistent.g)
increment(inconsistent::Inconsistent) = Threads.atomic_add!(inconsistent.k, 1)

function initialize!(inconsistent::Inconsistent, x₀::AbstractVector)
    inconsistent.k[] = 1.
    inconsistent.fval[] = 0.
    resize!(inconsistent.x,length(x₀))
    inconsistent.x .= Threads.Atomic{Float64}.(x₀)
    resize!(inconsistent.g,length(x₀))
    inconsistent.g .= Threads.Atomic{Float64}.(0.)
end

function (proxgrad::ProxGradient{<:MultiThread})(x₀::AbstractVector, loss::AbstractLoss, termination::AbstractTermination, logger::AbstractLogger)
    initialize!(proxgrad, x₀)

    Threads.@threads for i = 1:Threads.nthreads()
        kernel(proxgrad, Threads.threadid(), loss, termination, logger)
    end
end

getf(proxgrad::ProxGradient{<:MultiThread}) = getf(proxgrad.execution)
getx!(proxgrad::ProxGradient{<:MultiThread}) = nothing
getx(multithread::MultiThread) = readx(multithread)

function kernel(proxgrad::ProxGradient{<:MultiThread}, wid::Integer, loss::AbstractLoss, termination::AbstractTermination, logger::AbstractLogger)
    xlocal::Vector{Float64} = zeros(length(proxgrad.execution.x))
    glocal::Vector{Float64} = zeros(length(proxgrad.execution.g))

    read!(proxgrad, xlocal)
    flocal::Float64 = loss!(loss, xlocal, glocal)
    while true
        klocal::Int = getk(proxgrad.execution)
        if !iterate!(proxgrad, wid, klocal, flocal, glocal, termination, logger)
            return nothing
        end
        read!(proxgrad, xlocal)
        flocal = loss!(loss, xlocal, glocal)
    end
end

function iterate!(proxgrad::ProxGradient{<:MultiThread}, wid::Integer, klocal::Integer, flocal::Float64, glocal::AbstractVector, termination::AbstractTermination, logger::AbstractLogger, ::Val{false})
    k = getk(proxgrad.execution)
    x = proxgrad.execution.x
    g = proxgrad.execution.g
    if terminate(termination, k, flocal, readx(proxgrad.execution), glocal)
        setf!(proxgrad.execution, flocal)
        return false
    end
    boost!(POLO.boosting(proxgrad), wid, klocal, k, glocal, g)
    smooth!(POLO.smoothing(proxgrad), klocal, k, readx(proxgrad.execution), readg(proxgrad.execution), g)
    η = stepsize(POLO.stepsize(proxgrad), klocal, k, flocal, readx(proxgrad.execution), readg(proxgrad.execution))
    prox!(POLO.prox(proxgrad), η, readx(proxgrad.execution), readg(proxgrad.execution), x)
    increment(proxgrad.execution)
    return true
end
function iterate!(proxgrad::ProxGradient{Inconsistent}, wid::Integer, klocal::Integer, flocal::Float64, glocal::AbstractVector, termination::AbstractTermination, logger::AbstractLogger)
    res = iterate!(proxgrad, wid, klocal, flocal, glocal, termination, logger, Val{false}())
    if wid == 1
        log(logger, klocal, flocal, readx(proxgrad.execution), readg(proxgrad.execution))
    end
    return res
end

function iterate!(proxgrad::ProxGradient{<:MultiThread}, wid::Integer, klocal::Integer, flocal::Float64, glocal::AbstractVector, termination::AbstractTermination, logger::AbstractLogger, ::Val{true})
    lock(proxgrad.execution.mutex)
    res = iterate!(proxgrad, wid, klocal, flocal, glocal, termination, logger, Val{false}())
    log(logger, klocal, flocal, readx(proxgrad.execution), readg(proxgrad.execution))
    unlock(proxgrad.execution.mutex)
    return res
end
iterate!(proxgrad::ProxGradient{Consistent}, wid::Integer, klocal::Integer, flocal::Float64, glocal::AbstractVector, termination::AbstractTermination, logger::AbstractLogger) = iterate!(proxgrad, wid, klocal, flocal, glocal, termination, logger, Val{true}())

function read!(proxgrad::ProxGradient{<:MultiThread}, xlocal::AbstractVector, ::Val{false})
    xlocal .= readx(proxgrad.execution)
end
read!(proxgrad::ProxGradient{Inconsistent}, xlocal::AbstractVector) = read!(proxgrad, xlocal, Val{false}())

function read!(proxgrad::ProxGradient{<:MultiThread}, xlocal::AbstractVector, ::Val{true})
    lock(proxgrad.execution.mutex)
    read!(proxgrad, xlocal, Val{false}())
    unlock(proxgrad.execution.mutex)
end
read!(proxgrad::ProxGradient{Consistent}, xlocal::AbstractVector) = read!(proxgrad, xlocal, Val{true}())


text:
######################################################################
# binned_bearing.jl
#
# like bearing only, but observations are discretized into bins
######################################################################

export BinnedBearing
struct BinnedBearing <: Sensor
    noise_sigma::Float64		# noise std deviation (degrees)

    # for discretized measurements
    num_bins::Int
    bin_range::UnitRange{Int64}

    # TODO: technically must ensure noise_sigma is positive
    function BinnedBearing(noise_sigma::Real, num_bins::Int)
        new( float(noise_sigma), num_bins, 0:(num_bins-1))
    end
    BinnedBearing(noise_sigma::Real) = BinnedBearing(noise_sigma, 36)
    BinnedBearing() = BinnedBearing(10.0)
end

function observe(m::SearchDomain, s::BinnedBearing, p::Pose)
    truth = true_bearing(p, m.theta)
    noise = s.noise_sigma * randn()
    return obs2bin(mod(truth + noise, 360.0), s)
end

function O(bo::BinnedBearing, theta::LocTuple, xp::Pose, o::Int64)

    # Calculate true bearing, and find distance to bin edges
    ang_deg = true_bearing(xp, theta)
    rel_start, rel_end = rel_bin_edges(ang_deg, o, bo)

    # now look at probability
    sig = bo.noise_sigma
    p = my_cdf(sig, rel_end) - my_cdf(sig, rel_start)
    return p
end

# Find the relative offset
# TODO: must account for different discretizations
#function rel_bin_edges(bearing_deg, o::ObsBin, df::DF)
function rel_bin_edges(bearing_deg, o::ObsBin, bo::BinnedBearing)

    # calculate start, end degrees of bin
    start_deg, end_deg = bin2deg(o, bo)

    # compute relative distance to true bearing
    rel_start = fit_180(bearing_deg - start_deg)
    rel_end = fit_180(bearing_deg - end_deg)

    #rel_start = min(abs(rel_start), abs(rel_end))
    #rel_end = rel_start + 360.0 / df.num_bins

    # Make sure start is further left on number line
    if rel_end < rel_start
        temp = rel_start
        rel_start = rel_end
        rel_end = temp
    end

    # If we straddle the wrong point
    # Say df.num_bins = 10, and rel_start = -175, rel_end = 175
    # rel_end - rel_start would be 350 degrees, but this should be 10
    # so set rel_start to 175 and rel_end to 185
    # TODO: I'm pretty sure this doesn't work for df.num_bins = 2
    if (rel_end - rel_start) - 1e-3 > (360.0/bo.num_bins)
        rel_start = rel_end
        rel_end += 360.0/bo.num_bins
    end

    return rel_start, rel_end
end

# returns (start_deg, end_deg) integer tuple
function bin2deg(bin_deg::Int, bo::BinnedBearing)
    full_bin = 360.0 / bo.num_bins
    half_bin = full_bin / 2.0
    if bin_deg == 0
        start_val = -half_bin
        end_val = half_bin
    else
        start_val = full_bin * bin_deg - half_bin
        end_val  = full_bin * bin_deg + half_bin
    end
    return start_val, end_val
end

# 355 - 4.9999 = 0
# 5 - 14.9999 = 1
# 15 - 24.999 = 2
function obs2bin(o::Float64, s::BinnedBearing)
    full_bin = 360.0 / s.num_bins
    half_bin = full_bin / 2.0

    ob = round( Int, div((o + half_bin), full_bin) )
    if ob == s.num_bins
        ob = 0
    end
    return ob
end


text:
using Nemo
using AbstractAlgebra
using Distributions
using OffsetArrays

coefftype(ring::AbstractAlgebra.ResRing{T}) where T<:Union{nmod_poly,fmpz_mod_poly} = coefftype(ring.base_ring)
coefftype(ring::PolyRing) = ring.base_ring
Distributions.DiscreteUniform(R::Union{Nemo.NmodRing, AbstractAlgebra.Generic.ResRing}) = R
function (R::AbstractAlgebra.Generic.ResRing{T})(o::OffsetVector) where T<:Union{nmod_poly,fmpz_mod_poly}
    @assert first(axes(o)[1]) == 0
    R(base_ring(R)(collect(o)))
end
Base.oftype(x::AbstractAlgebra.Generic.Res, y::AbstractArray) = parent(x)(y)
function (R::NmodPolyRing)(o::OffsetVector)
    @assert first(axes(o)[1]) == 0
    R(collect(o))
end

struct CoeffView <: AbstractVector{Any}
    poly
    degree::Int
end
Base.getindex(p::CoeffView, i::Int) = coeff(p.poly, i)
Base.axes(p::CoeffView) = (0:p.degree,)
NTT.coeffs_primal(p::AbstractAlgebra.Generic.Res) = CoeffView(lift(p), degree(modulus(p)))


text:
using HackThatBase
using Base.Test

func(x,y) = x + y
args = lminfo(func, (Int, Float64))

cd(splitdir(Base.find_source_file("inference.jl"))[1])
@hack W inference
result = run_inference(W, args)
@test isa(result[2], CodeInfo)
@test result[3] == Float64

max_stars_repo_name: ihnorton/HackThatBase

max_issues_repo_name: ihnorton/HackThatBase

max_forks_repo_name: ihnorton/HackThatBase


text:
# write your own tests here
N = 10
ph = discrete_phs2(N, 0, 1);
num_freq = frequencies(ph)

@test abs(num_freq[1]- 3.51602) < 1e-3

N = 10;
K = 7;
a = 0;
b = 1;
ad = 0;
bd = 0.1;
phport = discrete_phs2_distports(N,K,a,b,ad,bd)
Bd = phport.Bd*ones(K,1)
@test abs(Bd[2] - 110.799) < 1e-2

max_stars_repo_name: flavioluiz/phjulia

max_issues_repo_name: flavioluiz/phjulia

max_forks_repo_name: flavioluiz/phjulia


text:
# transfer_operators.jl : tranform quantities from fine to coarse grid and back

# TransferKind
abstract type TransferKind end

struct Injection <: TransferKind end

struct FullWeighting <: TransferKind end

struct HalfWeighting <: TransferKind end

struct Cubic <: TransferKind end

Base.broadcastable(kind::TransferKind) = Ref(kind)

# TransferOperator
abstract type TransferOperator end

struct Restriction <: TransferOperator end

struct Interpolation <: TransferOperator end

# GridTransferOperator
struct GridTransferOperator{K<:TransferKind,O<:TransferOperator}
    kind::K
    op::O
end

# stencil definitions
stencil(op::GridTransferOperator{Injection,Restriction}) = @SVector [0,1.,0]

stencil(op::GridTransferOperator{FullWeighting,Restriction}) = @SVector [1/4,1/2,1/4]

stencil(op::GridTransferOperator{FullWeighting,Interpolation}) = @SVector [1/2,1,1/2]

stencil(op::GridTransferOperator{Cubic,Interpolation}) = @SVector [-1/16,0,9/16,1,9/16,0,-1/16]

# TODO: half weighting needs a 2d stencil...

# restriction
R₁(op::TransferKind,n::Int) = expand(GridTransferOperator(op,Restriction()),n) # 1d restriction operator
R(op::TransferKind,n::Int) = R₁(op,n)
R(op::TransferKind,n::Int...) = kron(R₁.(op,n)...)
restrict(u::Vector{T},op::TransferKind,n::Int...) where {T<:AbstractFloat} = R(op,n...)*u

# interpolation
P₁(op::TransferKind,n::Int) = expand(GridTransferOperator(op,Interpolation()),2n) # 1d interpolation operator
P(op::TransferKind,n::Int) = P₁(op,n)
P(op::TransferKind,n::Int...) = kron(P₁.(op,n)...)
interpolate(u::Vector{T},op::TransferKind,n::Int...) where {T<:AbstractFloat} = P(op,n...)*u

# prolongation = interpolation
prolongate(u,op) = interpolate(u,op)

# expand
function expand(op::GridTransferOperator{K,O} where {K,O},n)
    st = stencil(op)
    R = CartesianIndices((n-1,))
    I1, Iend = extrema(R)
    Is = Int[]
    Js = Int[]
    Vs = Float64[]
    @inbounds for I in R
        if iseven(I.I[1])
            is,js,vs = _stencil2mat(st,I,I1,Iend)
            push!(Is,is.>>1...)
            push!(Js,js...)
            push!(Vs,vs...)
        end
    end
    operator2matrix(op,Is,Js,Vs)
end

operator2matrix(op::GridTransferOperator{K,Interpolation} where {K},Is,Js,Vs) = sparse(Js,Is,Vs)
operator2matrix(op::GridTransferOperator{K,Restriction} where {K},Is,Js,Vs) = sparse(Is,Js,Vs)


text:
using RomeoApp
using Test

@testset "ROMEO function tests" begin

niread = RomeoApp.niread
savenii = RomeoApp.savenii

p = joinpath("data", "small")
phasefile_me = joinpath(p, "Phase.nii")
phasefile_me_nan = joinpath(p, "phase_with_nan.nii")
magfile_me = joinpath(p, "Mag.nii")
tmpdir = mktempdir()
phasefile_1eco = joinpath(tmpdir, "Phase.nii")
magfile_1eco = joinpath(tmpdir, "Mag.nii")
phasefile_1arreco = joinpath(tmpdir, "Phase.nii")
magfile_1arreco = joinpath(tmpdir, "Mag.nii")
maskfile = joinpath(tmpdir, "Mask.nii")
savenii(niread(magfile_me)[:,:,:,1] |> I -> I .> RomeoApp.MriResearchTools.median(I), maskfile)
savenii(niread(phasefile_me)[:,:,:,1], phasefile_1eco)
savenii(niread(magfile_me)[:,:,:,1], magfile_1eco)
savenii(niread(phasefile_me)[:,:,:,[1]], phasefile_1arreco)
savenii(niread(magfile_me)[:,:,:,[1]], magfile_1arreco)

phasefile_me_5D = joinpath(tmpdir, "phase_multi_channel.nii")
magfile_5D = joinpath(tmpdir, "mag_multi_channel.nii")
savenii(repeat(niread(phasefile_me),1,1,1,1,2), phasefile_me_5D)
savenii(repeat(niread(magfile_me),1,1,1,1,2), magfile_5D)

function test_romeo(args)
    folder = tempname()
    args = [args..., "-o", folder]
    try
        msg = unwrapping_main(args)
        @test msg == 0
        @test isfile(joinpath(folder, "unwrapped.nii"))
    catch e
        println(args)
        println(sprint(showerror, e, catch_backtrace()))
        @test "test failed" == "with error" # signal a failed test
    end
end

configurations_se(pf, mf) = vcat(configurations_se.([[pf], [pf, "-m", mf]])...)
configurations_se(pm) = [
    [pm...],
    [pm..., "-g"],
    [pm..., "-N"],
    [pm..., "-i"],
    [pm..., "-q"],
    [pm..., "-Q"],
    [pm..., "-m", magfile_me, "-u"],
    [pm..., "-w", "romeo"],
    [pm..., "-w", "bestpath"],
    [pm..., "-w", "1010"],
    [pm..., "-w", "101011"],
    [pm..., "--threshold", "4"],
    [pm..., "-s", "50"],
    [pm..., "-s", "50", "--merge-regions"],
    [pm..., "-s", "50", "--merge-regions", "--correct-regions"],
    [pm..., "--wrap-addition", "0.1"],
    [pm..., "-k", "robustmask"],
    [pm..., "-k", "nomask"],
    [pm..., "-k", "qualitymask"],
    [pm..., "-k", "qualitymask", "0.1"],
    [pm..., "-k", maskfile],
]
configurations_me(phasefile_me, magfile_me) = vcat(configurations_me.([[phasefile_me], [phasefile_me, "-m", magfile_me]])...)
configurations_me(pm) = [
    [pm..., "-e", "1:2", "-t", "[2,4]"], # giving two echo times for two echoes used out of three
    [pm..., "-e", "[1,3]", "-t", "[2,4,6]"], # giving three echo times for two echoes used out of three
    [pm..., "-e", "[1", "3]", "-t", "[2,4,6]"],
    [pm..., "-t", "[2,4,6]"],
    [pm..., "-t", "2:2:6"],
    [pm..., "-t", "[2.1,4.2,6.3]"],
    [pm..., "-t", "epi"], # shorthand for "ones(<num-echoes>)"
    [pm..., "-t", "epi", "5.3"], # shorthand for "5.3*ones(<num-echoes>)"
    [pm..., "-B", "-t", "[2,4,6]"],
    [pm..., "-B", "-t", "[2" ,"4", "6]"], # when written like [2 4 6] in command line
    [pm..., "--temporal-uncertain-unwrapping", "-t", "[2,4,6]"],
    [pm..., "--template", "1", "-t", "[2,4,6]"],
    [pm..., "--template", "3", "-t", "[2,4,6]"],
    [pm..., "--phase-offset-correction", "-t", "[2,4,6]"],
    [pm..., "--phase-offset-correction", "bipolar", "-t", "[2,4,6]"],
    [pm..., "--phase-offset-correction", "-t", "[2,4,6]", "--phase-offset-smoothing-sigma-mm", "[5,8,4]"],
    [pm..., "--phase-offset-correction", "-t", "[2,4,6]", "--write-phase-offsets"],
]

files = [(phasefile_1eco, magfile_1eco), (phasefile_1arreco, magfile_1arreco), (phasefile_1eco, magfile_1arreco), (phasefile_1arreco, magfile_1eco)]
for (pf, mf) in files, args in configurations_se(pf, mf)
    test_romeo(args)
end
for args in configurations_me(phasefile_me, magfile_me)
    test_romeo(args)
end
for args in configurations_se(["-p", phasefile_me, "-m", magfile_me, "-t", "[2,4,6]"])
    test_romeo(args)
end
for args in configurations_me(phasefile_me_5D, magfile_5D)[end-2:end] # test the last 3 configurations_me lines for coil combination
    test_romeo(args)
end
files_se = [(phasefile_1eco, magfile_1eco), (phasefile_1arreco, magfile_1arreco)]
for (pf, mf) in files_se
    b_args = ["-B", "-t", "3.06"]
    test_romeo(["-p", pf, b_args...])
    test_romeo(["-p", pf, "-m", mf, b_args...])
end

test_romeo([phasefile_me_nan, "-t", "[2,4]", "-k", "nomask"])

## Test error and warning messages
m = "multi-echo data is used, but no echo times are given. Please specify the echo times using the -t option."
@test_throws ErrorException(m) unwrapping_main(["-p", phasefile_me, "-o", tmpdir, "-v"])

m = "masking option '0.8' is undefined (Maybe '-k qualitymask 0.8' was meant?)"
@test_throws ErrorException(m) unwrapping_main(["-p", phasefile_1eco, "-o", tmpdir, "-v", "-k", "0.8"])

m = "masking option 'blub' is undefined"
@test_throws ErrorException(m) unwrapping_main(["-p", phasefile_1eco, "-o", tmpdir, "-v", "-k", "blub"])

m = "Phase offset determination requires all echo times!"
@test_throws ErrorException(m) unwrapping_main(["-p", phasefile_me_5D, "-o", tmpdir, "-v", "-t", "[1,2]", "-e", "[1,2]", "--phase-offset-correction"])

m = "echoes=[1,5]: specified echo out of range! Number of echoes is 3"
@test_throws ErrorException(m) unwrapping_main(["-p", phasefile_me, "-o", tmpdir, "-v", "-t", "[1,2,3]", "-e", "[1,5]"])

m = "echoes=[1,5} wrongly formatted!"
@test_throws ErrorException(m) unwrapping_main(["-p", phasefile_me, "-o", tmpdir, "-v", "-t", "[1,2,3]", "-e", "[1,5}"])

m = "Number of chosen echoes is 2 (3 in .nii data), but 5 TEs were specified!"
@test_throws ErrorException(m) unwrapping_main(["-p", phasefile_me, "-o", tmpdir, "-v", "-t", "[1,2,3,4,5]", "-e", "[1,2]"])

m = "size of magnitude and phase does not match!"
@test_throws ErrorException(m) unwrapping_main(["-p", phasefile_me, "-o", tmpdir, "-v", "-t", "[1,2,3]", "-m", magfile_1eco])

m = "robustmask was chosen but no magnitude is available. No mask is used!"
@test_logs (:warn, m) match_mode=:any unwrapping_main(["-p", phasefile_1eco, "-o", tmpdir])

@test_logs unwrapping_main(["-p", phasefile_1eco, "-o", tmpdir, "-m", magfile_1eco]) # test that no warning appears

## test no-rescale
readphase = RomeoApp.readphase
phasefile_me_uw = joinpath(tempname(), "unwrapped.nii")
phasefile_me_uw_wrong = joinpath(tempname(), "wrong_unwrapped.nii")
phasefile_me_uw_again = joinpath(tempname(), "again_unwrapped.nii")
unwrapping_main([phasefile_me, "-o", phasefile_me_uw, "-t", "[2,4,6]"])
unwrapping_main([phasefile_me_uw, "-o", phasefile_me_uw_wrong, "-t", "[2,4,6]"])
unwrapping_main([phasefile_me_uw, "-o", phasefile_me_uw_again, "-t", "[2,4,6]", "--no-rescale"])

@test readphase(phasefile_me_uw_again; rescale=false).raw == readphase(phasefile_me_uw; rescale=false).raw
@test readphase(phasefile_me_uw_wrong; rescale=false).raw != readphase(phasefile_me_uw; rescale=false).raw

## test mcpc3ds output files
testpath = joinpath(tmpdir, "test5d")
unwrapping_main([phasefile_me_5D, "-o", testpath, "-m", magfile_5D, "-t", "[2,4,6]"])
@test isfile(joinpath(testpath, "combined_mag.nii"))
@test isfile(joinpath(testpath, "combined_phase.nii"))

testpath = joinpath(tmpdir, "test4d")
unwrapping_main([phasefile_me, "-o", testpath, "-m", magfile_me, "-t", "[2,4,6]", "--phase-offset-correction"])
@test isfile(joinpath(testpath, "corrected_phase.nii"))

## test quality map
unwrapping_main([phasefile_me, "-m", magfile_me, "-o", tmpdir, "-t", "[2,4,6]", "-qQ"])
fns = joinpath.(tmpdir, ["quality.nii", ("quality_$i.nii" for i in 1:4)...])
for i in 1:length(fns), j in i+1:length(fns)
    @test niread(fns[i]).raw != niread(fns[j]).raw
end

end

## print version to verify
println()
unwrapping_main(["--version"])