text:
# A Simple Computational Algebra System

# scale
Base.:(-)(x::AbstractBlock{N}) where {N} = Scale(Val(-1), x)
Base.:(-)(x::Scale{Val{-1}}) = content(x)
Base.:(-)(x::Scale{Val{S}}) where S = Scale(Val(-S), content(x))
Base.:(-)(x::Scale) = Scale(-x.alpha, content(x))
Base.:(+)(x::AbstractBlock) = x

Base.:(*)(x::AbstractBlock, α::Number) = α * x

# NOTE: ±,±im should be identical
Base.:(*)(α::Val{S}, x::AbstractBlock) where S = Scale(α, x)

function Base.:(*)(α::T, x::AbstractBlock) where T <: Number
    return α ==  one(T) ? x                 :
    α == -one(T) ? Scale(Val(-1), x)   :
    α ==      im ? Scale(Val(im), x)        :
    α ==     -im ? Scale(Val(-im), x)       :
    Scale(α, x)
end

Base.:(*)(α::T, x::Scale) where {T <: Number} = α == one(T) ? x : Scale(x.alpha * α, content(x))
Base.:(*)(α::T, x::Scale{Val{S}}) where {T <: Number, S} = α * S * content(x)

Base.:(*)(α::Val{S}, x::Scale) where S = (S * x.alpha) * content(x)
Base.:(*)(α::Val{S1}, x::Scale{Val{S2}}) where {S1, S2} = (S1 * S2) * content(x)

Base.:(*)(x::Scale, y::Scale) = (x.alpha * y.alpha) * (content(x) * content(y))
Base.:(*)(x::Scale{Val{S1}}, y::Scale{Val{S2}}) where {S1, S2} = (S1 * S2) * (content(x) * content(y))
Base.:(*)(x::Scale, y::Scale{Val{S}}) where S = (x.alpha * S) * (content(x) * content(y))
Base.:(*)(x::Scale{Val{S}}, y::Scale) where S = (S * y.alpha) * (content(x) * content(y))
Base.:(*)(x::Scale, y::AbstractBlock) = x.alpha * chain(y, content(x))
Base.:(*)(y::AbstractBlock, x::Scale) = x.alpha * chain(content(x), y)

Base.:(+)(xs::AbstractBlock...) = Sum(xs...)
Base.:(*)(xs::AbstractBlock...) = chain(Iterators.reverse(xs)...)
Base.:(/)(A::AbstractBlock, x::Number) = (1/x)*A
# reduce
Base.sum(a::AbstractBlock{N}, blocks::AbstractBlock{N}...) where N = Sum(a, blocks...)
Base.prod(a::AbstractBlock{N}, blocks::AbstractBlock{N}...) where N = chain(Iterators.reverse(blocks)..., a)

Base.:(-)(lhs::AbstractBlock, rhs::AbstractBlock) = Sum(lhs, -rhs)
Base.:(^)(x::AbstractBlock, n::Int) = chain((copy(x) for k in 1:n)...)

"""
circuit optimisation
"""
module Optimise
using SimpleTraits
using YaoBlocks, YaoBlocks.ConstGate

export is_pauli, IsPauliGroup

"""
    IsPauliGroup{X}

Trait to check if `X` is an element of Pauli group.
"""
@traitdef IsPauliGroup{X}
IsPauliGroup(x) = IsPauliGroup{typeof(x)}()

"""
    is_pauli(x)

Check if `x` is an element of pauli group.

!!! note
    this function is just a binding of `SimpleTraits.istrait`, it will not work
    if the type is not registered as a trait with `@traitimpl`.
"""
is_pauli(x::T) where T = SimpleTraits.istrait(IsPauliGroup{T})
is_pauli(xs...) = all(is_pauli, xs)

for G in [:I2, :X, :Y, :Z]
    ImG = Symbol(:Im, G)
    nImG = Symbol(:nIm, G)
    nG = Symbol(:n, G)

    @eval const $ImG = im * $G
    @eval const $nImG = -im * $G
    @eval const $nG = -$G

    @eval @traitimpl IsPauliGroup{typeof($G)}
    @eval @traitimpl IsPauliGroup{typeof($ImG)}
    @eval @traitimpl IsPauliGroup{typeof($nImG)}
    @eval @traitimpl IsPauliGroup{typeof($nG)}
end

export merge_pauli
merge_pauli(x) = x
merge_pauli(x::AbstractBlock, y::AbstractBlock) = x * y

function merge_pauli(ex::ChainBlock{1})
    L = length(ex)
    new_ex = chain(1)

    # find all contiguous pauli and merge them
    # note we need to iterate in inverse order
    iterm = L
    while iterm > 0
        if iterm > 1 && is_pauli(ex[iterm], ex[iterm-1])
            pushfirst!(new_ex, merge_pauli(ex[iterm], ex[iterm-1]))
            iterm = iterm - 2
        else
            # search next
            pushfirst!(new_ex, ex[iterm])
            iterm -= 1
        end
    end
    return new_ex
end

merge_pauli(::XGate, ::YGate) = ImZ
merge_pauli(::XGate, ::ZGate) = -ImY
merge_pauli(::YGate, ::XGate) = -ImZ
merge_pauli(::YGate, ::ZGate) = ImX
merge_pauli(::ZGate, ::XGate) = ImY
merge_pauli(::ZGate, ::YGate) = nImX

for G in [:X, :Y, :Z]
    GT = Symbol(G, :Gate)

    @eval merge_pauli(::I2Gate, x::$GT) = x
    @eval merge_pauli(x::$GT, ::I2Gate) = x
    @eval merge_pauli(::$GT, ::$GT) = I2
end

merge_pauli(::I2Gate, ::I2Gate) = I2

export eliminate_nested
eliminate_nested(ex::AbstractBlock) = ex

# TODO: eliminate nested expr e.g chain(X, chain(X, Y))
function eliminate_nested(ex::T) where {T <: Union{ChainBlock, Sum}}
    _flatten(x) = (x, )
    _flatten(x::T) = subblocks(x)

    isone(length(ex)) && return first(subblocks(ex))
    return chsubblocks(ex, Iterators.flatten(map(_flatten, subblocks(ex))))
end

# temporary utils
_unscale(x::AbstractBlock) = x
_unscale(x::Scale) = content(x)
merge_alpha(alpha, x::AbstractBlock) = alpha
merge_alpha(alpha, x::Scale) = alpha * x.alpha
merge_alpha(alpha, x::Scale{Val{S}}) where S = alpha * S

# since we don't have T in blocks, this is a workaround
# to get correct identity in type stable term
merge_alpha(::Nothing, x::Scale) = x.alpha
merge_alpha(::Nothing, x::Scale{Val{S}}) where S = S
merge_alpha(::Nothing, x::AbstractBlock) = nothing

merge_scale(ex::AbstractBlock) = ex

# a simple function to find one for Val and Number
_one(x) = one(x)
_one(::Type{Val{S}}) where S = one(S)
_one(::Val{S}) where S = one(S)

export merge_scale

function merge_scale(ex::Union{Scale{S, N}, ChainBlock{N}}) where {S, N}
    alpha = nothing
    for each in subblocks(ex)
        alpha = merge_alpha(alpha, each)
    end
    ex = chsubblocks(ex, map(_unscale, subblocks(ex)))
    if alpha === nothing
        return ex
    else
        return alpha * ex
    end
end

export combine_similar

combine_similar(ex::AbstractBlock) = ex

combine_alpha(alpha, x) = alpha
combine_alpha(alpha, x::AbstractBlock) = alpha + 1
combine_alpha(alpha, x::Scale) = alpha + x.alpha
combine_alpha(alpha, x::Scale{Val{S}}) where S = alpha + S

function combine_similar(ex::Sum{N}) where N
    table = zeros(Bool, length(ex))
    list = []; p = 1
    while p <= length(ex)
        if table[p] == true
            # checked term, skip
            p += 1
        else
            # check similar term
            term = ex[p]
            table[p] = true # mark it in the table
            alpha = 1
            for (k, each) in enumerate(ex)
                if table[k] == true # checked term, skip
                    continue
                else
                    # check if unscaled term is the same
                    # merge them if they are
                    if _unscale(term) == _unscale(each)
                        alpha = combine_alpha(alpha, each)
                        # mark checked term in the table
                        table[k] = true
                    end
                end
            end

            # eliminate zeros
            if alpha != 0
                alpha = imag(alpha) == 0 ? real(alpha) : alpha
                alpha = isinteger(alpha) ? Integer(alpha) : alpha
                push!(list, alpha * term)
            end
        end
    end

    if isempty(list)
        return Sum{N}(())
    else
        return Sum(list...)
    end
end

export simplify

const __default_simplification_rules__ = Function[
    merge_pauli,
    eliminate_nested,
    merge_scale,
    combine_similar]

# Inspired by MasonPotter/Symbolics.jl
"""
    simplify(block[; rules=__default_simplification_rules__])

Simplify a block tree accroding to given rules, default to use
[`__default_simplification_rules__`](@ref).
"""
function simplify(ex::AbstractBlock; rules=__default_simplification_rules__)
    out1 = simplify_pass(rules, ex)
    out2 = simplify_pass(rules, out1)
    counter = 1
    while (out1 isa AbstractBlock) && (out2 isa AbstractBlock) && (out2 != out1)
        out1 = simplify_pass(rules, out2)
        out2 = simplify_pass(rules, out1)
        counter += 1
        if counter > 1000
            @warn "possible infinite loop in simplification rules. Breaking"
            return out2
        end
    end
    return out2
end

function simplify_pass(rules, ex)
    ex = chsubblocks(ex, map(x->simplify_pass(rules, x), subblocks(ex)))

    for rule in rules
        ex = rule(ex)
    end
    return ex
end

end


text:
using Documenter
using PCL

makedocs(
    modules = [PCL],
    clean   = false,
    format   = Documenter.Formats.HTML,
    sitename = "PCL.jl",
    pages = Any[
        "Home" => "index.md",
        "Installation" => "installation.md",
        "Getting started" => "getting-started.md",
        "Modules" => Any[
            "Common" => "modules/common.md"
            "Features" => "modules/features.md"
            "Filters" => "modules/filters.md"
            "IO" => "modules/io.md"
            "KDTree" => "modules/kdtree.md"
            "KeyPoints" => "modules/keypoints.md"
            "Octree" => "modules/octree.md"
            "Recognition" => "modules/recognition.md"
            "Registration" => "modules/registration.md"
            "SampleConsensus" => "modules/sampleconsensus.md"
            "Search" => "modules/search.md"
            "Segmentation" => "modules/segmentation.md"
            "Surface" => "modules/surface.md"
            "Tracking" => "modules/tracking.md"
            "Visualization" => "modules/visualization.md"
            "LibPCL" => "modules/libpcl.md"
            ]
        ],
)

#deploydocs(
#    target = "build",
#    deps = nothing,
#    make = nothing,
#    repo = "github.com/r9y9/WORLD.jl.git",
#)

ext: jl


text:
using Vulnerability
using Test

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


text:
module GraphEvolve

using Random

export

    # graphs
    AbstractGraph,
    AbstractNetwork,
    AbstractLattice,
    Observables,
    Network,
    Lattice2D,
    Lattice3D,

    # edge_methods
    choose_edge,
    add_edge!,
    nearest_neighbors,
    cart_to_lin,

    # cluster_methods
    get_cluster,
    get_largest_cluster_size,
    get_avg_cluster_size,
    get_largest_clusters,
    update_clusters!,
    update_cluster_sizes!,
    update_cluster_ids!,
    merge_clusters!,
    update_observables!,

    # evolution_methods
    erdos_renyi!,
    bohman_frieze!,
    product_rule!,
    stochastic_edge_acceptance!,

    # analysis_methods
    compute_delta!,
    finalize_observables!

# includes
include("./AbstractGraphs.jl")
include("./edge_methods.jl")
include("./cluster_methods.jl")
include("./evolution_methods.jl")
include("./analysis_methods.jl")

end # module


text:
# # Trade-off curves
using Random
Random.seed!(1)
m = 25;
n = 10;
A = randn(m, n);
b = randn(m, 1);

#-

using Convex, SCS, LinearAlgebra

gammas = exp10.(range(-4, stop = 2, length = 100));

x_values = zeros(n, length(gammas));
x = Variable(n);
for i in 1:length(gammas)
    cost = sumsquares(A * x - b) + gammas[i] * norm(x, 1)
    problem = minimize(cost, [norm(x, Inf) <= 1])
    solve!(problem, SCS.Optimizer; silent_solver = true)
    x_values[:, i] = evaluate(x)
end

#-

# Plot the regularization path.

using Plots
plot(
    title = "Entries of x vs lambda",
    xaxis = :log,
    xlabel = "lambda",
    ylabel = "x",
)
for i in 1:n
    plot!(gammas, x_values[i, :], label = "x$i")
end
plot!()


text:
let
    standard_duck_volume = 0.0014 # m^3
    @test isapprox(rubber_ducks_in_earth(0.0014), 3.5769429547897145e21)
end

let
    @test RUBBER_DUCK_PER_SEC == 200
end

let
    @test add_one(1) == 2
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -13.266290625856113 <= q <= 11.094199160852137)

                     add_NL_constraint(m, :(log(1 + exp(-0.2977454197977769 + 0.7722469146312885*log(1 + exp(0.9100824262656433 + 0.0704650714228956*log(1 + exp(-0.7110372864748489 + -0.3189000958940382*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.039301419137926175*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + -0.19790129907631693*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))) + 0.8966302258455006*log(1 + exp(0.7299768289984225 + 0.7045346620386517*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.8068588624409387*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + 0.8888276792068996*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))) + 0.7890296052704904*log(1 + exp(-0.5099367345236407 + -0.3810722083004521*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.4850891682061631*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + -0.158929561309892*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))) + 0.5815923034237733*log(1 + exp(-0.15598180224250147 + 0.7627892034632455*log(1 + exp(-0.7110372864748489 + -0.3189000958940382*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.039301419137926175*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + -0.19790129907631693*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))) + 0.6806349200893864*log(1 + exp(0.7299768289984225 + 0.7045346620386517*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.8068588624409387*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + 0.8888276792068996*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))) + -0.163059274878961*log(1 + exp(-0.5099367345236407 + -0.3810722083004521*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.4850891682061631*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + -0.158929561309892*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))) + -0.731486602090937*log(1 + exp(-0.9577813343303001 + 0.660774759467007*log(1 + exp(-0.7110372864748489 + -0.3189000958940382*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.039301419137926175*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + -0.19790129907631693*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))) + -0.16269598559087406*log(1 + exp(0.7299768289984225 + 0.7045346620386517*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.8068588624409387*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + 0.8888276792068996*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))) + 0.9308521240523335*log(1 + exp(-0.5099367345236407 + -0.3810722083004521*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.4850891682061631*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + -0.158929561309892*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))))) + log(1 + exp(0.23088458854123672 + 0.9025271253252471*log(1 + exp(0.9100824262656433 + 0.0704650714228956*log(1 + exp(-0.7110372864748489 + -0.3189000958940382*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.039301419137926175*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + -0.19790129907631693*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))) + 0.8966302258455006*log(1 + exp(0.7299768289984225 + 0.7045346620386517*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.8068588624409387*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + 0.8888276792068996*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))) + 0.7890296052704904*log(1 + exp(-0.5099367345236407 + -0.3810722083004521*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.4850891682061631*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + -0.158929561309892*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))) + 0.9612479731086991*log(1 + exp(-0.15598180224250147 + 0.7627892034632455*log(1 + exp(-0.7110372864748489 + -0.3189000958940382*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.039301419137926175*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + -0.19790129907631693*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))) + 0.6806349200893864*log(1 + exp(0.7299768289984225 + 0.7045346620386517*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.8068588624409387*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + 0.8888276792068996*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))) + -0.163059274878961*log(1 + exp(-0.5099367345236407 + -0.3810722083004521*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.4850891682061631*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + -0.158929561309892*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))) + 0.9633061422344626*log(1 + exp(-0.9577813343303001 + 0.660774759467007*log(1 + exp(-0.7110372864748489 + -0.3189000958940382*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.039301419137926175*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + -0.19790129907631693*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))) + -0.16269598559087406*log(1 + exp(0.7299768289984225 + 0.7045346620386517*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.8068588624409387*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + 0.8888276792068996*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))) + 0.9308521240523335*log(1 + exp(-0.5099367345236407 + -0.3810722083004521*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.4850891682061631*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + -0.158929561309892*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))))) + log(1 + exp(-0.8299198997441475 + 0.7709117851978275*log(1 + exp(0.9100824262656433 + 0.0704650714228956*log(1 + exp(-0.7110372864748489 + -0.3189000958940382*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.039301419137926175*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + -0.19790129907631693*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))) + 0.8966302258455006*log(1 + exp(0.7299768289984225 + 0.7045346620386517*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.8068588624409387*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + 0.8888276792068996*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))) + 0.7890296052704904*log(1 + exp(-0.5099367345236407 + -0.3810722083004521*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.4850891682061631*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + -0.158929561309892*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))) + 0.5338511667423695*log(1 + exp(-0.15598180224250147 + 0.7627892034632455*log(1 + exp(-0.7110372864748489 + -0.3189000958940382*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.039301419137926175*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + -0.19790129907631693*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))) + 0.6806349200893864*log(1 + exp(0.7299768289984225 + 0.7045346620386517*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.8068588624409387*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + 0.8888276792068996*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))) + -0.163059274878961*log(1 + exp(-0.5099367345236407 + -0.3810722083004521*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.4850891682061631*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + -0.158929561309892*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))) + 0.3695524647159538*log(1 + exp(-0.9577813343303001 + 0.660774759467007*log(1 + exp(-0.7110372864748489 + -0.3189000958940382*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.039301419137926175*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + -0.19790129907631693*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))) + -0.16269598559087406*log(1 + exp(0.7299768289984225 + 0.7045346620386517*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.8068588624409387*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + 0.8888276792068996*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))) + 0.9308521240523335*log(1 + exp(-0.5099367345236407 + -0.3810722083004521*log(1 + exp(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4]))) + 0.4850891682061631*log(1 + exp(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4]))) + -0.158929561309892*log(1 + exp(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module FrostHelperCustomFeatherModule

using ..Ahorn, Maple

@mapdef Entity "FrostHelper/CustomFeather" FrostHelperCustomFeatherPlacement(x::Integer, y::Integer, shielded::Bool=false, singleUse::Bool=false, respawnTime::Number=3.0, flyColor::String="ffd65c", spriteColor="White", flyTime::Number=2.0, maxSpeed::Number=190.0, lowSpeed::Number=140.0, neutralSpeed::Number=91.0, spritePath::String="objects/flyFeather/")

const colors = sort(collect(keys(Ahorn.XNAColors.colors)))

const placements = Ahorn.PlacementDict(
	"Custom Feather (Frost Helper)" => Ahorn.EntityPlacement(
		FrostHelperCustomFeatherPlacement
	)
)

Ahorn.editingOptions(entity::FrostHelperCustomFeatherPlacement) = Dict{String, Any}(
    "flyColor" => colors,
	"spriteColor" => colors
)

sprite = "objects/flyFeather/idle00.png"

function Ahorn.selection(entity::FrostHelperCustomFeatherPlacement)
    x, y = Ahorn.position(entity)
    return Ahorn.getSpriteRectangle(sprite, x, y)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::FrostHelperCustomFeatherPlacement, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, 0)

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
using SymmetricFormats, Test, LinearAlgebra

A = collect(reshape(1:9.0,3,3))

@testset "setindex!" begin
    APL = SymmetricPacked(A, :L)
    APL[1,1]=3
    @test APL[1,1] == 3
    @test_throws ArgumentError APL[1,3]=3
end

@testset "mul!" begin
    for uplo in [:U, :L]
        y = Float64[1,2,3]
        AP = SymmetricPacked(A, uplo)
        x = Float64[4,5,6]
        α = 0.5
        β = 0.5

        y_julia = α*AP*x + β*y
        y_blas = copy(y);  mul!(y_blas, AP, x, α, β)
        @test isapprox(y_julia, y_blas)
    end
end

@testset "packedsize" begin
    @test packedsize(A) == 6
    AP = SymmetricPacked(A)
    @test packedsize(AP) == 6
end

@testset "vector constructor" begin
    AP = SymmetricPacked(A)
    BP = SymmetricPacked(AP.tri)
    @test AP == BP
end

VERSION<v"1.8.0-DEV.1049" && include("blas.jl")


text:
using LinearAlgebra

ex=Ten([1,2,3],:i)*Ten([90,80,70],:i)
@test simplify(ex)==460

ex=Ten(:A,:i)*Ten(:A,:j)
@test simplify(ex&@equs(A=[1,2,3], j=i))==14

a=Alt([:i,:j])
@test a.x[1,2]==1
@test a.x[2,1]==-1
a=Alt([:i,:j,:k])
@test a.x[1,2,3]==1
@test a.x[3,2,1]==-1
a=Alt([:i,:j,:k,:l])
@test a.x[1,2,3,4]==1
@test a.x[4,2,3,1]==-1

ex=Alt([:i,:j,:k])*Ten(:a,:j)*Ten(:b,:k);a=[5,123,-12];b=[90,80,70];ex=ex&Equation(:a,a)&Equation(:b,b)
@test ex==Ten(Any[9570.0,-1430.0,-10670.0],Any[:i])
c=cross(a,b)
@test c[1]==ex&@equ(i=1)
@test c[2]==ex&@equ(i=2)
@test c[3]==ex&@equ(i=3)
ex=Alt([1,2,:k])*Ten([1,2,3],2)*Ten([90,80,70],:k);sex=simplify(ex)
@test sex==140
ex=Alt([1,3,:k])*Ten([1,2,3],3)*Ten([90,80,70],:k);sex=simplify(ex)
@test sex==-240
ex=Alt([1,:j,:k])*Ten([1,2,3],:j)*Ten([90,80,70],:k);sex=simplify(ex)
@test sex==-100
r=Alt([:i,:j,:k])*Ten([:a1,:a2,:a3],:j)*Ten([:b1,:b2,:b3],:k)&@equ i=1
@test r==simplify(:a2*:b3+-1.0*:a3*:b2)

ex=:c+Ten([:a1,:a2,:a3],:i)*Ten([:b1,:b2,:b3],:i)+:c;nex=simplify(ex)
@test nex==simplify(2*:c+:a1*:b1+:a2*:b2+:a3*:b3)

@test simplify(Ten(zeros(3,3).+Diagonal([1,1,1]),[:i,:i]))==3
@test simplify(Ten(zeros(4,4).+Diagonal([1,1,1,1]),[:i,:i]))==4
@test simplify(Ten(zeros(5,5).+Diagonal([1,1,1,1,1]),[:i,:i]))==5
r=Ten(:A,[:i,:i])&@equ A=[:a 0;0 :b]
@test r==:a+:b

io=IOBuffer()
print(io,simplify(Ten([1,2],:j)*Ten([3 2;1 -1],[:i,:j])))
str=String(take!(io))
#@test str=="Any[3,1](i) + 2 [2,-1](i)"
io=IOBuffer()
print(io,simplify(Ten([1,2],:j)*Ten([3 2;1 -1],[:j,:i])))
str=String(take!(io))
#@test str=="Any[3 2](i) + 2 [1 (-1)](i)"

@test duplicates([1,2,3,4,2])==[2,5]
@test duplicates([1,2,3,4,2,2])==[5,6]
@test duplicates([1,2,3],[0,3,0])==[3,2]
@test duplicates([1,3,3],[3,3,3])==[3,3]
@test Equations.arrduplicates([1,3,5],[12,241,13,3],[3,1,1])==([1,2],[2,4])
@test Equations.arrduplicates([1,3,55],[12,241,13],[33,1324,5,6,1])==([1,3],[1,5])

ex=Ten(:A,[:j,:i,:i])*Ten(:B,:j);r=ex&@equs(A=ones(3,3,3), B=[1,2,3])
@test r==18

eye3=zeros(3,3).+Diagonal([1,1,1])
ex=Ten(:A,[:i,:i]);eq=@equ A=$eye3;r=ex&eq
@test r==3

r=simplify(Ten([1,0],:i)+Ten([0,1],:i))
@test r==Ten([1,1],:i)
r=simplify(Ten([:c,:d],:i)+Ten([:a,:b],:i))
@test r==Ten(Any[:c+:a,:d+:b],Any[:i])
r=simplify(Alt([:i,:j,:k])*Ten([:a1,:a2,:a3],:j)*Ten([:b1,:b2,:b3],:k))

tt=Term[Factor[3,Ten([:a,:b],:i)],Factor[5,Ten([:c,:d],:i)]];stt=Equations.sumlify(tt)
@test stt[1][1].x[1]==simplify(Any[3*:a+5*:c,3*:b+5*:d][1])
@test stt[1][1].x[2]==simplify(Any[3*:a+5*:c,3*:b+5*:d][2])
#stt==simplify(Term[Factor[Equations.Ten(Any[3*:a+5*:c,3*:b+5*:d],Any[:i])]])

ex=Equations.sumlify(Equations.untensify!(sumconv(Alt([:i,:j])*Ten([10,100],:j)).terms))
@test sum(ex[1][1].x)==90
#ex==Term[Factor[Equations.Ten(Any[100.0,-10.0],Any[:i]),10]]

ex=Alt([:i,:j,:k])*Ten([1,0,0],:j)*Ten([0,0,1],:k);ex=sumconv(ex);ex=sumconv(ex);tt=Equations.untensify!(ex.terms);ttt=Equations.sumlify(tt)
@test ttt[1][1]&(@equ i=2)==-1

a=rand(Int,3).%9;b=rand(Int,3).%9
ab=[Equation(:a,a),Equation(:b,b)]
ex=Ten(:a,:i)*Ten(:b,:i)
@test ex&ab==dot(a,b)
ex=Alt([:i,:j,:k])*Ten(:a,:j)*Ten(:b,:k)
c=cross(a,b)
x=(ex&ab).x
@test c==x
ex=Alt([:i,:j,:k])*Ten(:a,:j)*Ten(:b,:k)*Ten(:a,:m)*Ten(:b,:m)
r=cross(a,b)*dot(a,b)
x=(ex&ab).x
@test r==x
ex=Alt([:i,:j,:k])*Ten(:a,:j)*Ten(:b,:k)*Ten(:a,:m)*Ten(:b,:m)+Ten(:b,:i)
r=cross(a,b)*dot(a,b)+b
x=(ex&ab).x
@test r==x
c=rand(Int,3).%9;d=rand(Int,3).%9
abcd=[Equation(:a,a),Equation(:b,b),Equation(:c,c),Equation(:d,d)]
ex=Alt([:i,:j,:k])*Ten(:a,:j)*Ten(:b,:k)*Alt([:m,:n,:o])*Ten(:c,:n)*Ten(:d,:o)
r=dot(cross(a,b),cross(c,d))
x=ex&abcd&Equation(:m,:i)
@test r==x

B=[:b1 :b2;:b3 :b4];ex=Alt([:i,:j])*Alt([:k,:l])*Ten(B,[:i,:k])*Ten(B,[:j,:l]);r=simplify(ex)
@test r==simplify(2*:b1*:b4+-2.0*:b2*:b3)


text:
cd("/Users/adrienpajon/Documents/all.nosync/git_repo/pulse_analysis/2020-11-12 xp pulse all seuil/Julia_analysis")

using Plots

using StatsPlots

using JLD

# using Pyplot
pyplot()

include("../functions/analysisPULSE.jl")

####################################################################
gel=load("results/PULSE_analysis_gel.jld")
pleine=load("results/PULSE_analysis_pleine.jld")
pneu=load("results/PULSE_analysis_pneu.jld")

filenameXsensor=[gel["filenameXsensor"];pleine["filenameXsensor"];pneu["filenameXsensor"]]
Energysignal_Xsensor_seat=merge(gel["Energysignal_Xsensor_seat"],pleine["Energysignal_Xsensor_seat"],pneu["Energysignal_Xsensor_seat"])
rms_Xsensor_seat=merge(gel["rms_Xsensor_seat"],pleine["rms_Xsensor_seat"],pneu["rms_Xsensor_seat"])
vdv_Xsensor_seat=merge(gel["vdv_Xsensor_seat"],pleine["vdv_Xsensor_seat"],pneu["vdv_Xsensor_seat"])
crest_factor_Xsensor_seat=merge(gel["crest_factor_Xsensor_seat"],pleine["crest_factor_Xsensor_seat"],pneu["crest_factor_Xsensor_seat"])
std_Xsensor_seat=merge(gel["std_Xsensor_seat"],pleine["std_Xsensor_seat"],pneu["std_Xsensor_seat"])
min_Xsensor_seat=merge(gel["min_Xsensor_seat"],pleine["min_Xsensor_seat"],pneu["min_Xsensor_seat"])
max_Xsensor_seat=merge(gel["max_Xsensor_seat"],pleine["max_Xsensor_seat"],pneu["max_Xsensor_seat"])
range_Xsensor_seat=merge(gel["range_Xsensor_seat"],pleine["range_Xsensor_seat"],pneu["range_Xsensor_seat"])

splitName,
plot_order,
perturbLabel,
perturbLabel_unique,
perturbIndices,
wheelType,
wheelLabel_unique,
wheelIndices=splitFileName(filenameXsensor)

bool_perturb_choice="4cm"#"3cm"
bool_perturb=perturbLabel.==bool_perturb_choice

filenameXsensor_=filenameXsensor[bool_perturb]

splitName,
plot_order,
perturbLabel,
perturbLabel_unique,
perturbIndices,
wheelType,
wheelLabel_unique,
wheelIndices=splitFileName(filenameXsensor_,0,0.05)

include("../functions/analysisPULSE.jl")
plotAnalysisPULSE_All(
        Energysignal_Xsensor_seat,
        rms_Xsensor_seat,
        vdv_Xsensor_seat,
        crest_factor_Xsensor_seat,
        filenameXsensor_,
        "Xsensor",
        wheelLabel_unique,
        plot_order,
        wheelIndices,
        [:blue :red :green :orange :yellow],
        [:rect :rect :rect :rect :circle],
        string("\n","perturbation : ",bool_perturb_choice))

include("../functions/analysisPULSE.jl")
plotAnalysisPULSE_All_2(
            std_Xsensor_seat,
            min_Xsensor_seat,
            max_Xsensor_seat,
            range_Xsensor_seat,
            filenameXsensor_,
            "Xsensor",
            ["x";"y"],
            wheelLabel_unique,
            plot_order,
            wheelIndices,
            [:blue :red :green :orange :yellow],
            [:rect :rect :rect :rect :circle],
            string("\n","perturbation : ",bool_perturb_choice))

println("PULSE Analysis Finished")

# using JLD
# save("results/PULSE_analysis_gel.jld",
#         "filenameXsensor_seat",filenameXsensor_seat,
#         "sensorXsensNames",sensorXsensNames,
#         "step_max_perturb_Xsensor_seat",step_max_perturb_Xsensor_seat,
#         "Energysignal_Xsensor_seat",Energysignal_Xsensor_seat,
#         "rms_Xsensor_seat",rms_Xsensor_seat,
#         "vdv_Xsensor_seat",vdv_Xsensor_seat,
#         "crest_factor_Xsensor_seat",crest_factor_Xsensor_seat,
#         "Energysignal_Xsensor_seat_polar",Energysignal_Xsensor_seat_polar,
#         "rms_Xsensor_seat_polar",rms_Xsensor_seat_polar,
#         "vdv_Xsensor_seat_polar",vdv_Xsensor_seat_polar)


text:
"""
    Flux.testmode!(lh::LearnedHeuristic, mode = true)

Make it possible to change the mode of the LearnedHeuristic: training or testing. This makes sure
you stop updating the weights or the approximator once the training is done. It is used at the beginning
of the `train!` function to make sure it's training and changed automatically at the end of it but a user can
manually change the mode again if he wants.
"""
Flux.testmode!(lh::LearnedHeuristic, mode = true) = Flux.testmode!(lh.agent, mode)

"""
    update_with_cpmodel!(lh::LearnedHeuristic{SR, R, A}, model::CPModel)

This function initializes the fields of a LearnedHeuristic which are useful to do reinforcement learning
and which depend on the CPModel considered. It is called at the beginning of the `search!` function (in the
InitializingPhase).
"""
function update_with_cpmodel!(lh::LearnedHeuristic{SR, R, A}, model::CPModel) where {
    SR <: AbstractStateRepresentation,
    R <: AbstractReward,
    A <: ActionOutput
}

    # construct the action_space
    valuesOfVariables = sort(branchable_values(model))

    lh.action_space = valuesOfVariables
    # state rep construction
    lh.current_state = SR(model; action_space=lh.action_space)

    # create and initialize the reward
    lh.reward = R(model)

    lh.search_metrics = SearchMetrics(model)

    lh
end

"""
    sync!(lh::LearnedHeuristic, model::CPModel, x::AbstractIntVar)

Synchronize the environment part of the LearnedHeuristic with the CPModel.
"""
function sync_state!(lh::LearnedHeuristic, model::CPModel, x::AbstractIntVar)
    update_representation!(lh.current_state, model, x)
    nothing
end

"""
    get_observation!(lh::LearnedHeuristic, model::CPModel, x::AbstractIntVar, done = false)

This function retrieve all the elements that are useful for doing reinforcement learning.
- The reward, which has been incremented through the set_reward!(PHASE, ...) functions.
- The terminal boolean (is the episode done or not)
- The current state, which is the Array representation of a StateRepresentation of the CPModel at this stage of the solving.
This state is what will then be given to the agent to make him proposed an action (a value to assign.)
- The agent uses indexes of branchable values instead of values directly in order to facilitate the interface understanding
- The legal_actions & legal_actions_mask used to make sure the agent won't propose a value which isn't
in the domain of the variable we're branching on.

The result is given as a namedtuple for convenience with ReinforcementLearning.jl interface.
"""
function get_observation!(lh::LearnedHeuristic, model::CPModel, x::AbstractIntVar, done = false)
    #æet action_space_index. CAUTION : the action space is never updated.
    action_space_index=collect(1:size(lh.action_space)[1])

    # get legal_actions_mask
    legal_actions_mask = [value in x.domain for value in lh.action_space]

    # compute legal actions
    legal_actions = lh.action_space[legal_actions_mask]

    #get reward
    reward = lh.reward.value

    # Initialize reward for the next state: not compulsory with DefaultReward, but maybe useful in case the user forgets it
    model.statistics.AccumulatedRewardBeforeReset += lh.reward.value
    lh.reward.value = 0

    # synchronize state:
    sync_state!(lh, model, x)
    state = trajectoryState(lh.current_state)

    if !wears_mask(lh)
        return unmaskedCPEnv(reward, done, state, action_space_index)
    end
    # return the observation as a named tuple (useful for interface understanding)
    return CPEnv(reward, done, state, action_space_index, legal_actions, legal_actions_mask)
end

"""
    set_metrics!(PHASE::T, lh::LearnedHeuristic, model::CPModel, symbol::Union{Nothing, Symbol}, x::Union{Nothing, AbstractIntVar}) where T <: LearningPhase

Call set_metrics!(::SearchMetrics, ...) on env.search_metrics to simplify synthax.
Could also add it to basicheuristic !
"""
function set_metrics!(PHASE::T, lh::LearnedHeuristic, model::CPModel, symbol::Union{Nothing, Symbol}, x::Union{Nothing, AbstractIntVar}) where T <: LearningPhase
    set_metrics!(PHASE, lh.search_metrics, model, symbol, x::Union{Nothing, AbstractIntVar})
end

function wears_mask(valueSelection::LearnedHeuristic)
    if (hasfield(typeof(valueSelection.agent.policy.learner.approximator),:actor))
        wears_mask(valueSelection.agent.policy.learner.approximator.actor)      #A2C
    else
        wears_mask(valueSelection.agent.policy.learner.approximator.model)              #DQN
    end
end

"""
    from_order_to_id(state::AbstractArray, value_order::Int64)

Return the ids of the valid indexes from the Array representation of the AbstractStateRepresentation. Used to be able to work with
ActionOutput of variable size (VariableOutput).
"""
function from_order_to_id(state::AbstractTrajectoryState, value_order::Int64, SR::Type{<:AbstractStateRepresentation})
    @assert !isnothing(state.possibleValuesIdx)
    return state.possibleValuesIdx[value_order]
end

"""
    action_to_value(vs::LearnedHeuristic{SR, R, VariableOutput}, action::Int64, state::AbstractArray, model::CPModel)

Mapping action taken to corresponding value when handling VariableOutput type of ActionOutput.
"""
function action_to_value(vs::LearnedHeuristic{SR, R, VariableOutput}, action::Int64, state::AbstractTrajectoryState, model::CPModel) where {SR <: DefaultStateRepresentation, R}
    value_id = from_order_to_id(state, action, SR)
    cp_vertex = cpVertexFromIndex(vs.current_state.cplayergraph, value_id)
    @assert isa(cp_vertex, ValueVertex)
    return cp_vertex.value
end

function action_to_value(vs::LearnedHeuristic{SR, R, VariableOutput}, action::Int64, state::AbstractTrajectoryState, model::CPModel) where {SR <: TsptwStateRepresentation, R}
    return from_order_to_id(state, action, SR)
end

"""
    action_to_value(vs::LearnedHeuristic{SR, R, FixedOutput}, action::Int64, state::AbstractArray, model::CPModel)

Mapping index of Q-value vector to value in the action space when using a FixedOutput.
"""
function action_to_value(vs::LearnedHeuristic{SR, R, FixedOutput}, action::Int64, state::AbstractTrajectoryState, model::CPModel) where {SR <: DefaultStateRepresentation, R}
    return vs.action_space[action]
end

"""
    function branchable_values(cpmodel::CPModel)

Return an array of all possible values taken by the variables we can branch on.
"""
function branchable_values(cpmodel::CPModel)
    setOfValues = Set{Int}()
    for (k, x) in branchable_variables(cpmodel)
        for value in x.domain
            push!(setOfValues, value)
        end
    end
    return collect(setOfValues)
end


text:
# =================================================================
# Quadrotor model
# See https://easychair.org/publications/paper/gjfh
# =================================================================

using Reachability: Options, SafeStatesProperty
using MathematicalSystems, LazySets
using DynamicPolynomials, SemialgebraicSets

# ==============================
# Load model
# ==============================

"""
    quadrotor(; [T], [X0], [variables], [controller_inputs])

Construct the Quadrotor model.

### Input

- `T`                 -- (optional, default: `5.0`) time horizon
- `X0`                -- (optional, default: position is uncertain in all directions
                          within `[-0.4, 0.4]m` and velocity is uncertain in all directions
                          within `[-0.4, 0.4]m/s`) set of initial states
- `variables`         -- (optional, default: `PolyVar` variables) the set of polynomal
                          variables that are used in the equations
- `controller_inputs` -- (optional, default: `(1.0, 0.0, 0.0)`) tuple with
                         the controller inputs `u₁`, `u₂` and `u₃` which
                         correspond to the desired values for height, roll and
                         pitch respectively

### Output

The tuple `(𝑃, 𝑂)` where `𝑃` is an initial-value problem and `𝑂` are the options.

### Notes

 Physical variables:

-  x₁ : interitial (north) position
-  x₂ : intertial (east) position
-  x₃ : altitude
-  x₄ : longitudinal velocity
-  x₅ : lateral velocity
-  x₆ : vertical velocity
-  x₇ : roll angle
-  x₈ : pitch angle
-  x₉ : yaw angle
-  x₁₀ : roll rate
-  x₁₁ : pitch rate
-  x₁₂ : yaw rate

### Specification

The task is to change the height from `0[m]` to `1[m]` within `5[s]`. A goal region
`[0.98, 1.02]` of the height `x₃` has to be reached within `5[s]` and the height
has to stay below 1.4 for all times. After `1[s]` the height should stay above
`0.9[m]`.
"""
function quadrotor(; T=5.0,
                     X0=Hyperrectangle(zeros(12), [fill(0.4, 6); fill(0.0, 6)]),
                     variables=(@polyvar x[1:12]),
                     controller_inputs=(1.0, 0.0, 0.0))

    # parameters of the model
    g = 9.81           # gravity constant in m/s^2
    R = 0.1            # radius of center mass in m
    l = 0.5            # distance of motors to center mass in m
    Mrotor = 0.1       # motor mass in kg
    M = 1.             # center mass in kg
    m = M + 4*Mrotor   # total mass in kg

    # moments of inertia
    Jx = 2/5*M*R^2 + 2*l^2*Mrotor
    Jy = Jx
    Jz = 2/5*M*R^2 + 4*l^2*Mrotor

    𝑂 = Options()

    # unrwap the variables and the inputs
    x₁, x₂, x₃, x₄, x₅, x₆, x₇, x₈, x₉, x₁₀, x₁₁, x₁₂ = variables[1] # or variables
    u₁, u₂, u₃ = controller_inputs

    𝑂[:variables] = variables
    𝑂[:vars] = [1:12;]

    # equations of the controllers
    F = m*g - 10*(x₃ - u₁) + 3*x₆  # height control
    τϕ = -(x₇-u₂) - x₁₀            # roll control
    τθ = -(x₈ - u₃) - x₁₁          # pitch control
    τψ = 0                         # heading is uncontrolled

    # differential equations for the quadrotor
    ẋ₁ = cos(x₈)*cos(x₉)*x₄ + (sin(x₇)*sin(x₈)*cos(x₉) - cos(x₇)*sin(x₉))*x₅
         + (cos(x₇)*sin(x₈)*cos(x₉) + sin(x₇)*sin(x₉))*x₆
    ẋ₂ = cos(x₈)*sin(x₉)*x₄ + (sin(x₇)*sin(x₈)*sin(x₉) + cos(x₇)*cos(x₉))*x₅
         + (cos(x₇)*sin(x₈)*sin(x₉) - sin(x₇)*cos(x₉))*x₆
    ẋ₃ = sin(x₈)*x₄ - sin(x₇)*cos(x₈)*x₅ - cos(x₇)*cos(x₈)*x₆
    ẋ₄ = x₁₂*x₅ - x₁₁*x₆ - g*sin(x₈)
    ẋ₅ = x₁₀*x₆ - x₁₂*x₄ + g*cos(x₈)*sin(x₇)
    ẋ₆ = x₁₁*x₄ - x₁₀*x₅ + g*cos(x₈)*cos(x₇) - F/m
    ẋ₇ = x₁₀ + sin(x₇)*tan(x₈)*x₁₁ + cos(x₇)*tan(x₈)*x₁₂
    ẋ₈ = cos(x₇)*x₁₁ - sin(x₇)*x₁₂
    ẋ₉ = (sin(x₇)/cos(x₈))*x₁₁ + (cos(x₇)/cos(x₈))*x₁₂
    ẋ₁₀ = (Jy - Jz)/Jx * x₁₁ * x₁₂ + τϕ/Jx
    ẋ₁₁ = (Jz - Jx)/Jy * x₁₀ * x₁₂ + τθ/Jy
    ẋ₁₂ = (Jx - Jy)/Jz * x₁₀ * x₁₁ + τψ/Jz

    𝐹 = PolynomialContinuousSystem(f)

    # instantiate the IVP
    𝑃 = InitialValueProblem(𝐹, X0)

    # time horizon
    𝑂[:T] = T

    # variables to plot
    𝑂[:plot_vars] = [0, 3]

    # safety property
    # TODO: add specification
    #𝑂[:property] = LinearConstraintProperty([0, 0, 0, -1., 0, 0, 0], -unsafe_bound)
    # @set x₄ ≥ 0.01, vars=(x₁, x₂, x₃, x₄, x₅, x₆, x₇)
    return (𝑃, 𝑂)
end

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
using InvariantMeasures

@testset "Chebyshev assembler" begin

N = 16
B = InvariantMeasures.Chebyshev(N)

D = mod1_dynamic(x->2*x)
L(ϕ, x) = (ϕ(x/2)+ϕ(x/2+0.5))/2

M = assemble(B, D)
using LinearAlgebra

all_true = true
for i in 1:N
   w = mid.(L.(B[i], B.p))
   z = InvariantMeasures.chebtransform(w)
   all_true = all_true && norm(z-M[:, i], Inf)< 10^-13
end

@test all_true

end


text:
#Copyright Nathan Parrish 2018
#University of Colorado, Boulder
#Colorado Center for Astrodynamics Research
#Celestial and Spaceflight Mechanics Laboratory
#
#Function to do direct multiple shooting in CRTBP
#
#This version:
#   - Direct method.
#   - Can optimize endpoints and time of flight.
#   - Works great for "short" transfers, but not for "long" transfers. Uses
#     linear endpoint constraints w/ quadratic endpoint penalty function
#     parameterized by β input term.
#   - Allows impulsive maneuvers at endpoints.
#   - Fixed step integrator, with mesh refinement
#
#Inputs:
#   - X_all     = [6-7 x n_nodes] Initial guess array of states. Each column is
#                 the state at time t[i]. [Nondimensional units of CRTBP]
#   - u_all     = [3 x n_nodes] Initial guess array of controls. [Newtons]
#   - τ1        = Initial guess of position on the initial orbit. In range [0-1]
#   - τ2        = Initial guess of position on the final orbit. In range [0-1]
#   - t_TU      = Guess vector of times [Nondimensional time units of CRTBP]
#   - dV1       = Initial guess of initial impulsive maneuver
#   - dV2       = Initial guess of final impulsive maneuver
#   - MU        = μ for CRTBP
#   - DU        = Distance unit for CRTBP [km]
#   - TU        = Time unit for CRTBP [sec]
#   - n_nodes   = Number of nodes
#   - nsteps    = Number of intergration steps per node
#   - mass      = Initial spacecraft mass [kg]
#   - Isp       = Specific impulse [sec]
#   - X0_times  = Vector of times on which X0_states are defined (usually
#                 linspace(0, 1, n_nodes))
#   - X0_states = [6 x n_nodes] States given to define initial orbit.
#   - Xf_times  = Vector of times on which Xf_states are defined (usually
#                 linspace(0, 1, n_nodes))
#   - Xf_states = [6 x n_nodes] States given to define final orbit.
#   - plot_yn   = Boolean, true to plot solution at every iteration.
#   - flagEnd   = Boolean, true to allow the endpoints and time of flight to
#                 vary.
#   - β         = Float to determine magnitude of quadratic endpoint cost.
#   - allowImpulsive = Boolean, true to allow impulsive maneuvers at endpoints.
#   - maxIter   = Max number of iterations.
#
#Outputs: (X_all, u_all, τ1, τ2, t_TU, dV1, dV2, defect)
#   - X_all     = [6-7 x n_nodes] Solution array of states.
#   - u_all     = [3 x n_nodes] Solution array of controls.
#   - τ1        = Solution position on the initial orbit. In range [0-1]
#   - τ2        = Solution position on the final orbit. In range [0-1]
#   - t_TU      = Solution vector of times [Nondimensional time units of CRTBP]
#   - dV1       = Initial guess of initial impulsive maneuver
#   - dV2       = Initial guess of final impulsive maneuver
#   - defect    = Differential defects at solution.

function multiShoot_CRTBP_direct(X_all, u_all, τ1, τ2, t_TU, dV1, dV2, MU, DU, TU, n_nodes,
    nsteps, mass, Isp, X0_times, X0_states, Xf_times, Xf_states, plot_yn,
    flagEnd, β, allowImpulsive, maxIter)

    global day

    ## DEFINE SUB-FUNCTIONS

    function defectCalc(X_all, u_all, t_TU, nstate, n_nodes, nsteps, Isp, odefun)
        #Calculate defect constraint violations

        #find mid-point times:
        t_mid_TU = t_TU[1:end-1] + diff(t_TU)/2

        #pre-allocate defect:
        defect1 = zeros(nstate, n_nodes-1)

        #pre-allocate integration errors:
        errors = zeros(n_nodes-1)
        for i = 1:n_nodes-1
            #propagate node i forward from time i to time between (i+1) and (i)
            #propagate node i+1 backward from time i+1 to time between (i+1) and (i)

            ########### propagate forward:
            x0 = copy(X_all[:,i])
            control = copy(u_all[:,i])
            tspan = linspace(t_TU[i], t_mid_TU[i], nsteps)
            time_direction = 1.0 #forward
            (state_for, maxErr_for) = ode7_8(odefun, tspan, x0, MU, DU, TU, Isp, control, time_direction)

            ########### propagate backward:
            #can use the same tspan because we just need the length of time to be correct
            x0 = copy(X_all[:,i+1])
            #reverse velocity for going backward:
            x0[4:6] = -x0[4:6]
            control = copy(u_all[:,i+1]) #should this be negative? No.
            time_direction = -1.0 #backward
            (state_back, maxErr_back) = ode7_8(odefun, tspan, x0, MU, DU, TU, Isp, control, time_direction)
            stateF_back = state_back[:,end]
            #reverse velocity to be forward again:
            stateF_back[4:6] = -stateF_back[4:6]

            ########## defect between forward and backward results:
            defect1[:,i] = state_for[:,end] - stateF_back

            #save largest error:
            errors[i] = max(maxErr_for, maxErr_back)
        end

        #outputs:
        (defect1, errors)
    end

    function jacobianCalc(X_all, u_all, t_TU, defect, nstate, n_nodes, nsteps, Isp, odefun, pert)
        #Calculate Jacobian of defect constraints with respect to all states and
        #controls

        #number of variables that affect the defect between 2 nodes:
        nvar = 2*(nstate+3)

        # --- JACOBIAN (found first as a smaller, dense matrix) ---
        #The i'th defect vector is affected by X_all[:,i:i+1] and by control[:,i:i+1]

        #size of Jacobian is [# defects,    # variables]
        Jac_temp = zeros(nstate*(n_nodes-1), nvar)
        for i = 1:n_nodes-1

            XU_nom = vcat(X_all[:, i:i+1][:], u_all[:, i:i+1][:])
            for j = 1:nvar
                #create modified arrays
                XU_mod = copy(XU_nom)
                XU_mod[j] = XU_mod[j] + pert

                #reshape
                X_temp = reshape(XU_mod[1:(nstate*2)], nstate, 2)
                u_temp = reshape(XU_mod[(nstate*2)+1 : end], 3, 2)

                #compute single defect between the two affected nodes:
                (defect_temp, errors) = defectCalc(X_temp, u_temp, t_TU[i:i+1], nstate, 2, nsteps, Isp, odefun)

                #index of the rows which correspond to this defect
                idx = (i-1)*nstate+1 : i*nstate
                Jac_temp[idx,j] = (defect_temp - defect[:,i]) / pert;

            end
        end

        #re-create full Jacobian (band diagonal):
        Jac_full = zeros(nstate*(n_nodes-1), n_nodes*(nstate+3));
        for ind = 1:size(Jac_full,1) #loop through rows
            #Which node's defect are we on?
            idx = Int32(ceil(ind/(nstate)));

            #offset (in # of columns) for the state-related ones
            idx2 = (idx-1)*(nstate);

            #offset (in # of columns) for the control-related ones
            idx3 = nstate*n_nodes + (idx-1)*(3);

            #set the appropriate columns of Jac_full from Jac_temp:
            #state-related terms:
            Jac_full[ind, (idx2+1):(idx2+2*nstate)] = Jac_temp[ind,1:2*nstate];
            #control-related terms:
            Jac_full[ind, (idx3+1):(idx3+2*3)] = Jac_temp[ind,(2*nstate+1):end];
        end

        #outputs:
        Jac_full
    end

    function endpointPartials(X_inner, u_all, τ1, τ2, t_TU, dV1, dV2, Jac_XU,
        defect, X0_times, X0_states, Xf_times, Xf_states)
        #Computes partial derivatives of the defect constraints with respect to
        #the endpoint parameters (τ1 and τ2) and time of flight.

        (state_0, state_f) = interpEndStates(τ1, τ2, X0_times, X0_states,
            Xf_times, Xf_states, MU)

        ############ partial of all defects wrt tf
        pert = 1e-5
        t0_TU = copy(t_TU[1])
        tf_TU = copy(t_TU[end])
        tf_TU_mod = t_TU[end] + pert

        tau = (t_TU - t0_TU) / (tf_TU - t0_TU) * 2 - 1 #transform t -> tau
        t_TU_mod = t0_TU + (tau+1)/2*(tf_TU_mod-t0_TU) #transform tau -> t

        (defect_mod, errors) = defectCalc(X_all, u_all, t_TU_mod, nstate,
            n_nodes, nsteps, Isp, odefun)

        ddefect_dt = (defect_mod[:] - defect[:]) / pert

        ############ partial of defect1 wrt dV1, and defectN wrt dV2
        pert = 1e-5 # (DU/TU)
        d_defect1_dV1 = zeros(nstate, 3)
        d_defectN_dV2 = zeros(nstate, 3)
        for ind = 1:3

            #apply impulsive delta-v at beginning and end
            state_0_mod = state_0
            state_0_mod[ind+3] += pert

            state_f_mod = state_f
            state_f_mod[ind+3] += pert

            ############ defect1 wrt dV1
            X_temp = hcat(state_0_mod, X_inner[:,1])
            u_temp = u_all[:, 1:2]

            (defect_temp, errors) = defectCalc(X_temp, u_temp, t_TU[1:2], nstate, 2, nsteps, Isp, odefun)

            d_defect1_dV1[:,ind] = (defect_temp - defect[:,1]) / pert

            ############ defectN wrt dV2
            X_temp = hcat(X_inner[:,end], state_f_mod)
            u_temp = u_all[:, end-1 : end]

            (defect_temp, errors) = defectCalc(X_temp, u_temp, t_TU[end-1 : end], nstate, 2, nsteps, Isp, odefun)

            d_defectN_dV2[:,ind] = (defect_temp - defect[:,end]) / pert
        end

        ############ Combine together with the rest of the Jacobian:
        n1 = (n_nodes-2)*nstate + n_nodes*3

        Jac_full = hcat(Jac_XU, zeros(nstate*(n_nodes-1), 9) )

        #partial of defect1 wrt τ1
        Jac_full[1:nstate, n1 + 1] = d_defect1_dp1

        #partial of defectN wrt τ2
        Jac_full[end-(nstate-1):end, n1 + 2] = d_defectN_dp2

        #partial of all defects wrt tf
        Jac_full[:, n1 + 3] = ddefect_dt

        #partial of defect1 wrt dV1
        Jac_full[1:nstate, (n1+4):(n1+6)] = d_defect1_dV1

        #partial of defectN wrt dV2
        Jac_full[end-(nstate-1):end, (n1+7):(n1+9)] = d_defectN_dV2

        ############ Outputs:
        Jac_full
    end

    function optimizeTraj(X_all, u_all, τ1, τ2, tf_TU, tau, dV1, dV2, defect, Jac_full,
        nstate, n_nodes, X0_times, X0_states, Xf_times, Xf_states, flagEnd, β)
        #Optimizes 1 iteration of the QP/SOCP problem.

        # global day

        #JuMP linear/quadratic constrained solver
        #
        #See link for more solvers:
        # http://www.juliaopt.org/JuMP.jl/0.18/installation.html#getting-solvers
        #Any solver capable of solving SOCP or NLP problems should work. The
        #following options have been used successfully. Gurobi is very fast, but
        #has a commercial license (free for academics). Ipopt is usable, but
        #has suffered frequently from bugs when updating Julia.
        #
        # m = Model(solver=GurobiSolver(BarConvTol=1e-10, OutputFlag=0))
        m = Model(solver=IpoptSolver(print_level=0))

        #State update variables:
        @variable(m,    X_jump[1:length(X_all)] )
        #constrain initial mass:
        if nstate == 7
            @constraint(m, (X_jump[7] + X_all[7, 1]) == mass)
        end

        #control update variables:
        #limits on u_jump have to be determined by u_all
        @variable(m, u_jump[1:length(u_all)])

        #endpoint parameter variables:
        if flagEnd
            d = 0.1 #unfreeze endpoints
        else
            d = 0.0 #freeze endpoints
        end
        @variable(m, (- d) <= p1_jump <= (d))
        @variable(m, (- d) <= p2_jump <= (d))

        #time of flight variable:
        if flagEnd
            d = 1 * day / TU #unfreeze endpoints
        else
            d = 0.0 #freeze endpoints
        end
        d = 0.
        @variable(m, -d <= tf_jump <= d)
        @constraint(m, (tf_TU + tf_jump) <= 40 * day/TU) #absolute max bound
        @constraint(m, (tf_TU + tf_jump) >= 0.0 * day/TU) #absolute min bound

        #impulsive maneuvers at endpoints:
        @variable(m, dV1_jump[1:3])
        @variable(m, dV2_jump[1:3])
        if !allowImpulsive
            @constraint(m, dV1_jump .== 0.) #disallow impulsive maneuvers
            @constraint(m, dV2_jump .== 0.) #
        end

        #Thrust magnitude constraints:
        #
        #Note: Thrust magnitude constraint is disabled because it was found
        #more effective to constrain thrust with the indirect method (used after
        #converging on the direct solution).
        #
        #Since u_jump is just the *update* to control, the value we want to actually
        #limit to be within constraints is (u_jump + u_all).
        #Create a vector of symbolic quadratic constraints:
        # umag = (u_jump[1:3:end] + u_all[1:3:end]).^2 +
        #     (u_jump[2:3:end] + u_all[2:3:end]).^2 +
        #     (u_jump[3:3:end] + u_all[3:3:end]).^2
        # @constraint(m, umag .<= thrustLimit^2 )

        #Need 'dt' so that unequal time steps gives the right result. If we
        #don't use 'dt' like this, then high thrust for a long segment is
        #weighted equally with high thrust for a short segment.
        t_TU_fixed = t0_TU + (tau+1)/2*(tf_TU-t0_TU) #transform tau -> t
        dt = diff(t_TU_fixed)
        dt_temp = vcat(dt/2, dt[end]/2) + vcat(0, dt[1:end-1]/2, 0)
        dt_repeated = ones(3,1) * dt_temp' #same size as control

        #replace this by actual mass vector (repeated like dt)
        if nstate == 7
            mass_repeated = repmat(X_all[7,:], 3, 1)
        else
            mass_repeated = 1000.0 * ones(size(u_all))
        end

        #Dynamics constraints:
        #constrain bringing the defects to zero (according to linearization)
        @constraint(m, -Jac_full*vcat(X_jump, u_jump, tf_jump) .== defect[:])

        #Constraints for endpoints:

        #Finite difference derivatives of endpoints wrt parameters
        pert = 0.05 #needs to be relatively large because the endpoint orbits are only specified at 100 points (so, every 0.01 tau)
        (state_0, state_f) = interpEndStates(τ1, τ2, X0_times, X0_states, Xf_times, Xf_states, MU)
        (state_0_mod1, state_f_mod1) = interpEndStates(τ1+pert, τ2+pert, X0_times, X0_states, Xf_times, Xf_states, MU)
        (state_0_mod2, state_f_mod2) = interpEndStates(τ1-pert, τ2-pert, X0_times, X0_states, Xf_times, Xf_states, MU)
        dstate0_dp1 = (state_0_mod1 - state_0_mod2) / (2*pert)
        dstatef_dp2 = (state_f_mod1 - state_f_mod2) / (2*pert)
        #2nd derivatives:
        ddstate0_dp1 = (state_0_mod1 - 2*state_0 + state_0_mod2) / (pert^2)
        ddstatef_dp2 = (state_f_mod1 - 2*state_f + state_f_mod2) / (pert^2)

        if flagEnd
            #linear:
            linearEnd1 = (state_0 + dstate0_dp1*(p1_jump))
            linearEnd2 = (state_f + dstatef_dp2*(p2_jump))
            @constraint(m, (X_all[1:6, 1]   + X_jump[1:6]                     + vcat(zeros(3), dV1_jump+dV1) ) - linearEnd1 .== 0 )
            @constraint(m, (X_all[1:6, end] + X_jump[(end-nstate+1):end][1:6] + vcat(zeros(3), dV2_jump+dV2) ) - linearEnd2 .== 0 )

            #quadratic:
            quadEnd1 = (state_0 + dstate0_dp1*(p1_jump) + ddstate0_dp1 / 2 * (p1_jump)^2 )
            quadEnd2 = (state_f + dstatef_dp2*(p2_jump) + ddstatef_dp2 / 2 * (p2_jump)^2 )
            c1_quad = (X_all[1:6,1]   + X_jump[1:6])                     - (state_0 + dstate0_dp1*(p1_jump) + ddstate0_dp1 / 2 * (p1_jump)^2 )
            c2_quad = (X_all[1:6,end] + X_jump[(end-nstate+1):end][1:6]) - (state_f + dstatef_dp2*(p2_jump) + ddstatef_dp2 / 2 * (p2_jump)^2 )

            #Difference from linear:
            quadCost = norm(ddstate0_dp1) * 1/2 * (p1_jump)^2 +
                       norm(ddstatef_dp2) * 1/2 * (p2_jump)^2

        else #hard fix endpoints
            quadCost = 0.
            (state_0, state_f) = interpEndStates(τ1, τ2, X0_times, X0_states, Xf_times, Xf_states, MU)

            @constraint(m, (X_all[1:6, 1]   + X_jump[1:6]                     + vcat(zeros(3), dV1_jump+dV1) ) - state_0 .== 0 )
            @constraint(m, (X_all[1:6, end] + X_jump[(end-nstate+1):end][1:6] + vcat(zeros(3), dV2_jump+dV2) ) - state_f .== 0 )
        end

        costEnd = sum( ( (dV1 + dV1_jump) * DU/TU ).^2) + sum( ( (dV2 + dV2_jump) * DU/TU ).^2)

        #quadratic cost, with quadratic endpoint constraint violation cost added:
        cost = sum( ((u_all[:] + u_jump).^2) .* dt_repeated[:] ) + β * quadCost + costEnd

        @objective(m, Min, cost)

        #Run QP solver:
        status = JuMP.solve(m)

        x_update = reshape(getvalue(X_jump), nstate, n_nodes)
        u_update = reshape(getvalue(u_jump), 3, n_nodes)

        p1_update = getvalue(p1_jump)
        p2_update = getvalue(p2_jump)
        tf_update = getvalue(tf_jump)
        dV1_update = getvalue(dV1_jump)
        dV2_update = getvalue(dV2_jump)

        tf_TU = tf_TU + getvalue(tf_jump)

        cost = getvalue(cost)

        #Outputs:
        (x_update, u_update, p1_update, p2_update, tf_update, dV1_update, dV2_update, status, cost, dstate0_dp1, dstatef_dp2, ddstate0_dp1, ddstatef_dp2)
    end

    function lineSearch(X_all, x_update, u_all, u_update, t_TU, nstate, n_nodes, nsteps, Isp, odefun)
        #Very simple line search function. Seeks to minimize the sum of squared
        #defect constraint violations, while forcing some step to be taken.

        alpha = 1.0

        alpha_all = linspace(0.1, 1, 10)

        er = zeros(size(alpha_all))

        for ind = 1:length(alpha_all)
            alpha = alpha_all[ind]

            X_all_trial = X_all + x_update*alpha
            u_all_trial = u_all + u_update *alpha

            #check defects
            (defect, errors) = defectCalc(X_all_trial, u_all_trial, t_TU, nstate, n_nodes, nsteps, Isp, odefun)

            er[ind] = sum(defect[:].^2)
        end

        #Output:
        alpha = alpha_all[ er .== minimum(er) ]
        alpha = alpha[1]
    end

    function interpEndStates(τ1, τ2, X0_times, X0_states, Xf_times, Xf_states, MU)
        #Interpolate the initial and final states from the given set of states

        #check bounds on endpoint parameters (wrap to 0,1)
        while τ1 > 1
            τ1 -= 1
        end
        while τ1 < 0
            τ1 += 1
        end
        while τ2 > 1
            τ2 -= 1
        end
        while τ2 < 0
            τ2 += 1
        end

        #build interpolation object:
        X0_itp = Interpolations.scale(interpolate(X0_states', (BSpline(Cubic(Natural())), NoInterp()), OnGrid()), X0_times, 1:6)
        Xf_itp = Interpolations.scale(interpolate(Xf_states', (BSpline(Cubic(Natural())), NoInterp()), OnGrid()), Xf_times, 1:6)

        #endpoints from fixed parameters given:
        state_0 = [X0_itp[τ1, 1], X0_itp[τ1, 2], X0_itp[τ1, 3], X0_itp[τ1, 4], X0_itp[τ1, 5], X0_itp[τ1, 6]]
        state_f = [Xf_itp[τ2, 1], Xf_itp[τ2, 2], Xf_itp[τ2, 3], Xf_itp[τ2, 4], Xf_itp[τ2, 5], Xf_itp[τ2, 6]]

        #outputs:
        (state_0, state_f)
    end

################################################################################
    ## Done defining sub-functions. Now, run the algorithm.

    #number of state elements (should be 6 or 7)
    nstate = size(X_all, 1)

    #derivatives function to use
    odefun = CRTBP_prop_EP_deriv

    pert = 1e-8;

    #find mid-point times:
    t_mid_TU = t_TU[1:end-1] + diff(t_TU)/2

    t0_TU = copy(t_TU[1])
    tf_TU = copy(t_TU[end])
    tau = (t_TU - t0_TU) / (tf_TU - t0_TU) * 2 - 1; #transform t -> tau

    (state_0, state_f) = interpEndStates(τ1, τ2, X0_times, X0_states, Xf_times, Xf_states, MU)

    ############ First nominal run:
    (defect, errors) = defectCalc(X_all, u_all, t_TU, nstate, n_nodes, nsteps, Isp, odefun)

    #iterate!
    iterCount = 0;
    er = 1.0 #initialize error
    while er > 1e-6
        iterCount += 1
        if iterCount > maxIter
            warn("Reached max iteration count at ", iterCount," iterations")
            break
        end

        ############ Compute Jacobian for dynamics constraints
        Jac_full = jacobianCalc(X_all, u_all, t_TU, defect, nstate, n_nodes, nsteps, Isp, odefun, pert)

        #Partial wrt tf
        pert_tf = 1e-3
        t0_TU = copy(t_TU[1])
        tf_TU = copy(t_TU[end])
        tf_TU_mod1 = tf_TU + pert_tf
        tf_TU_mod2 = tf_TU - pert_tf
        t_TU_mod1 = t0_TU + (tau+1)/2*(tf_TU_mod1-t0_TU) #transform tau -> t
        t_TU_mod2 = t0_TU + (tau+1)/2*(tf_TU_mod2-t0_TU) #transform tau -> t
        (defect_mod1, errors) = defectCalc(X_all, u_all, t_TU_mod1, nstate, n_nodes, nsteps, Isp, odefun)
        (defect_mod2, errors) = defectCalc(X_all, u_all, t_TU_mod2, nstate, n_nodes, nsteps, Isp, odefun)

        ddefect_dt = (defect_mod1 - defect_mod2) / (2*pert_tf)

        Jac_full = hcat(Jac_full, ddefect_dt[:])

        ############ OPTIMIZE QP SUBPROBLEM

        #every other iteration, allow the endpoints to vary
        flagEnd_temp = copy(flagEnd)
        if mod(iterCount,2) == 0
            flagEnd_temp = false * flagEnd
        end

        (x_update, u_update, p1_update, p2_update, tf_update, dV1_update, dV2_update,
            status, cost, dstate0_dp1, dstatef_dp2, ddstate0_dp1, ddstatef_dp2) =
            optimizeTraj(X_all, u_all, τ1, τ2, tf_TU, tau, dV1, dV2, defect, Jac_full,
            nstate, n_nodes, X0_times, X0_states, Xf_times, Xf_states, flagEnd_temp, β)

        #Plot endpoint constraint lines:
        #check the linear/quadratic approximations of the initial and final orbits:
        (state_0, state_f) = interpEndStates(τ1, τ2, X0_times, X0_states, Xf_times, Xf_states, MU)

        num_test = 30
        x1_1 = zeros(6, num_test)
        x2_1 = zeros(6, num_test)
        x1_2 = zeros(6, num_test)
        x2_2 = zeros(6, num_test)
        temp = linspace(-1, 1, num_test)
        for ind = 1:num_test
            #linear:
            x1_1[:,ind] = (state_0 + dstate0_dp1*(p1_update)*temp[ind])
            x2_1[:,ind] = (state_f + dstatef_dp2*(p2_update)*temp[ind])

            #quadratic:
            x1_2[:,ind] = (state_0 + dstate0_dp1*(p1_update)*temp[ind] + ddstate0_dp1/2*( (p1_update)*temp[ind] )^2)
            x2_2[:,ind] = (state_f + dstatef_dp2*(p2_update)*temp[ind] + ddstatef_dp2/2*( (p2_update)*temp[ind] )^2)
        end

        ############ LINE SEARCH
        alpha = 1.0
        #Call line search:
        if iterCount > 10
            alpha = lineSearch(X_all, x_update, u_all, u_update, t_TU, nstate, n_nodes, nsteps, Isp, odefun)
        end

        X_all = X_all + x_update * alpha
        u_all = u_all + u_update * alpha
        τ1 = τ1 + p1_update * alpha
        τ2 = τ2 + p2_update * alpha
        tf_TU = tf_TU + tf_update * alpha
        dV1 = dV1 + dV1_update * alpha
        dV2 = dV2 + dV2_update * alpha

        ############ PLOT UPDATE
        if plot_yn
            #Can plot either with Julia (slow) or with MATLAB (less slow)

            #Plot with Julia:
            plotTrajPlotly_direct(X_all, zeros(size(u_all)), X0_states, Xf_states, 0.2)
            gui()
        end

        t_TU = t0_TU + (tau+1)/2*(tf_TU-t0_TU) #transform tau -> t

        ############ CHECK UPDATE
        (defect, errors) = defectCalc(X_all, u_all, t_TU, nstate, n_nodes, nsteps, Isp, odefun)

        #print progress report
        er = norm(defect[:], Inf)
        @printf "Iter %d. Max defect = %.2e. Cost = %.5f. tf = %.2f days. α = %.3f.\n" iterCount er cost tf_TU*TU/day alpha
    end

    #outputs:
    (X_all, u_all, τ1, τ2, t_TU, dV1, dV2, defect)
end

function meshRefine_direct(X_all, u_all, t_TU, nstate, n_nodes, nsteps, Isp, odefun)
    ############ MESH REFINEMENT

    (defect, errors) = defectCalc(X_all, u_all, t_TU, nstate, n_nodes, nsteps, Isp, odefun)

    #These can be adjusted up or down depending on the problem.
    tol_min = 1e-20
    tol_max = 1e-18

    println("Starting with ",n_nodes," nodes.")
    n_nodes_old = n_nodes

    # First, remove extra points
    while minimum(errors) < tol_min
        #find the index of the minimum & maximum errors
        idx_min = find(errors .== minimum(errors))[1]

        #ensure that we dont remove the endpoints with this.
        if idx_min == 1
            idx_min = 2
        end

        #Remove the node because it's just making things messy:
        X_all = hcat(X_all[:, 1:(idx_min-1)], X_all[:, (idx_min+1):end])
        u_all = hcat(u_all[:, 1:(idx_min-1)], u_all[:, (idx_min+1):end])
        t_TU = vcat(t_TU[1:(idx_min-1)], t_TU[(idx_min+1):end])
        n_nodes = n_nodes - 1;

        #Update the defects & error estimates:
        (defect, errors) = defectCalc(X_all, u_all, t_TU, nstate, n_nodes, nsteps, Isp, odefun)
    end

    if ~(n_nodes == n_nodes_old)
        println("Removed nodes, now n_nodes = ", n_nodes)
    end

    n_nodes_old2 = n_nodes

    #Add more points where needed:
    while maximum(errors) > tol_max

        #find the index of the minimum & maximum errors
        idx_max = find(errors .== maximum(errors))[1]

        #find mid-point times:
        t_mid_TU = t_TU[1:end-1] + diff(t_TU)/2

        #new time to add:
        t_new = t_TU[idx_max] + diff( t_TU[idx_max:idx_max+1] )[1]/2

        #Propagate to the midpoint:
        x0 = X_all[:,idx_max]
        control = u_all[:,idx_max]
        tspan = [t_TU[idx_max],    t_new ]
        time_direction = 1.0 #forward
        state_new = ode7(odefun, tspan, x0, MU, DU, TU, Isp, control, time_direction)

        #new control is the average of the control on either side:
        u_new = (u_all[:, idx_max] + u_all[:, idx_max+1]) / 2

        #add a node because it's not good enough
        X_all = hcat(X_all[:, 1:idx_max], state_new[:,2], X_all[:,(idx_max+1:end)])
        u_all = hcat(u_all[:, 1:idx_max], u_new, u_all[:,(idx_max+1:end)])
        t_TU    = vcat(t_TU[1:idx_max], t_new, t_TU[idx_max+1:end])
        n_nodes = n_nodes + 1;

        #Update the defects & error estimates:
        (defect, errors) = defectCalc(X_all, u_all, t_TU, nstate, n_nodes, nsteps, Isp, odefun)
    end

    if !(n_nodes == n_nodes_old2)
        println("Added nodes, now n_nodes = ", n_nodes)
    end
    if n_nodes == n_nodes_old
        println("Did not need to add or remove nodes.")
    else
        println("Refined the mesh. Now have ",n_nodes," nodes.")
    end

    #outputs:
    (X_all, u_all, t_TU, n_nodes)
end

function plotTrajPlotly_direct(X_all, u_all, X0_states, Xf_states, u_scale)
    #Plot trajectory with PlotlyJS package

    Plots.plot(X0_states[1,:], X0_states[2,:], X0_states[3,:], w=3,
        xlabel="X (DU)", ylabel="Y (DU)", zlabel="Z (DU)",
        label="X₀")
    Plots.plot!(Xf_states[1,:], Xf_states[2,:], Xf_states[3,:], w=2, label="Xf")
    Plots.plot!([1-MU], [0.], [0.], m=(0,2,:black) , label="Moon")
    Plots.plot!(X_all[1,:],X_all[2,:],X_all[3,:], w=3, color=:black, label="Transfer trajectory")

    arrows_start = X_all[1:3,:]
    arrows_end = X_all[1:3,:] + u_all*u_scale

    #remove zero thrust vectors:
    arrows_start = arrows_start[:, norm_many(u_all) .!== 0.]
    arrows_end   = arrows_end[:, norm_many(u_all) .!== 0.]

    arrows_x = vcat(arrows_start[1,:]', arrows_end[1,:]')
    arrows_y = vcat(arrows_start[2,:]', arrows_end[2,:]')
    arrows_z = vcat(arrows_start[3,:]', arrows_end[3,:]')

    plot!([1-MU,], [0,], [0.2,], marker=1, color=:white)

    Plots.plot!(arrows_x, arrows_y, arrows_z, color=:red, show=true)

    #create data points for a sphere surface:
    (x,y,z) = sphere(32)

    #scale and move the sphere:
    r_moon = 1737. #km
    x = (x * r_moon/DU) + (1-MU)
    y = y * r_moon/DU
    z = z * r_moon/DU

    #plot the sphere as a surface:
    Plots.surface!(x,y,z, c=:blues)
end

max_stars_repo_name: travelingspaceman/LowThrustOpt

max_issues_repo_name: travelingspaceman/LowThrustOpt

max_forks_repo_name: travelingspaceman/LowThrustOpt


text:
export rotate, rotate!

"""
    rotate!(corrs,azi,baz)

Rotate cross-correlations from E,N,Z to R,T,Z coordinate system.
# Arguments
- `corrs::Array{CorrData,1}`: Array of EE, EN, EZ, NE, NN, NZ, ZE, ZN, ZZ CorrData
          for a single station-pair.
- `azi::Real`: Azimuth (in radians) from station 1 to station 2
- `baz::Real`: Back azimuth (in radians) from station 1 to station 2
"""
function rotate!(corrs::Array{CorrData,1},azi::Real,baz::Real)
    T = eltype(corrs[1].corr)
    azi = T(azi)
    baz = T(baz)
    N = length(corrs)
    # get time intersect of all components
    ts = [c.t for c in corrs]
    inter = intersect(ts)[1]

    # subset so all corrs match in time
    for ii = 1:N
        ind = findall(x -> x ∈ inter, corrs[ii].t)
        corrs[ii].t = corrs[ii].t[ind]
        corrs[ii].corr = corrs[ii].corr[:,ind]
    end

    comp_dict = Dict()
    for ii = 1:N
        comp_dict[corrs[ii].comp] = ii
    end

    # get EE EN NN NE comps for rotation to TT RR TR RT
    # see Lin et al., 2008: Surface wave tomography of the western
    # United States from ambient seismic noise, GJI
    big_rot = [[-cos(azi)*cos(baz) cos(azi)*sin(baz) -sin(azi)*sin(baz) sin(azi)*cos(baz)];
               [-sin(azi)*sin(baz) -sin(azi)*cos(baz) -cos(azi)*cos(baz) -cos(azi)*sin(baz)];
               [-cos(azi)*sin(baz) -cos(azi)*cos(baz) sin(azi)*cos(baz) sin(azi)*sin(baz)];
               [-sin(azi)*cos(baz) sin(azi)*sin(baz) cos(azi)*sin(baz) -cos(azi)*cos(baz)]]
    to_rot = Array{eltype(corrs[1].corr)}(undef,4,size(corrs[1].corr)[1],length(inter))
    to_rot[1,:,:] = corrs[comp_dict["EE"]].corr
    to_rot[2,:,:] = corrs[comp_dict["EN"]].corr
    to_rot[3,:,:] = corrs[comp_dict["NN"]].corr
    to_rot[4,:,:] = corrs[comp_dict["NE"]].corr
    rotated = rotation_kernel(big_rot,to_rot)

    # update correlations from EE EN NN NE to TT RR TR RT
    update_corr!(corrs[comp_dict["EE"]],rotated[1,:,:],"TT")
    update_corr!(corrs[comp_dict["EN"]],rotated[2,:,:],"RR")
    update_corr!(corrs[comp_dict["NN"]],rotated[3,:,:],"TR")
    update_corr!(corrs[comp_dict["NE"]],rotated[3,:,:],"RT")

    # rotate from EZ NZ to RZ TZ
    azi_rot = [[cos(azi) -sin(azi)];
               [sin(azi) cos(azi)]]
    to_rot = Array{eltype(corrs[1].corr)}(undef,2,size(corrs[1].corr)[1],length(inter))
    to_rot[1,:,:] = corrs[comp_dict["NZ"]].corr
    to_rot[2,:,:] = corrs[comp_dict["EZ"]].corr
    rotated = rotation_kernel(azi_rot,to_rot)
    update_corr!(corrs[comp_dict["NZ"]],rotated[1,:,:],"RZ")
    update_corr!(corrs[comp_dict["EZ"]],rotated[2,:,:],"TZ")

    # rotate from ZE ZN to ZR ZT
    baz_rot = [[-cos(baz) -sin(baz)];
               [sin(baz) -cos(baz)]]
    to_rot = Array{eltype(corrs[1].corr)}(undef,2,size(corrs[1].corr)[1],length(inter))
    to_rot[1,:,:] = corrs[comp_dict["ZN"]].corr
    to_rot[2,:,:] = corrs[comp_dict["ZE"]].corr
    rotated = rotation_kernel(baz_rot,to_rot)
    update_corr!(corrs[comp_dict["ZN"]],rotated[1,:,:],"ZR")
    update_corr!(corrs[comp_dict["ZE"]],rotated[2,:,:],"ZT")
    # update ZZ rotation
    corrs[comp_dict["ZZ"]].rotated = true

    return nothing
end
rotate(corrs::Array{CorrData,1},azi::Real,baz::Real) =
                      (U = deepcopy(corrs);rotate!(U,azi,baz);return U)

"""
  rotation_kernel(rot_matrix, to_rotate)

Rotation function for 2D correlation functions.

# Arguments
- `rot_matrix::Array`: Rotation matrix (either 2x2 or 4x4).
- `to_rotate::Array`: Array to be rotated column-wise.
"""
function rotation_kernel(rot_matrix::AbstractArray, to_rotate::AbstractArray)
    rotated = similar(to_rotate)
    for ii = 1:size(to_rotate,2)
        for jj = 1:size(to_rotate,3)
            rotated[:,ii,jj] .= rot_matrix * to_rotate[:,ii,jj]
        end
    end
    return rotated
end

"""
  update_corr!(C,comp)

Change component, name and corr of rotated correlation `C` to component `comp`.

# Arguments
- `C::CorrData`: CorrData object.
- `new_corr::Array`: Rotated correlation matrix.
- `comp::String`: New component after rotation, e.g. "RR"
"""
function update_corr!(C::CorrData,new_corr::AbstractArray,comp::String)
    # change comp of correlation
    C.comp = comp

    # change name of correlation
    newname = split(C.name,'.')
    newname[4] = newname[4][1:end-1] * comp[1]
    newname[8] = newname[8][1:end-1] * comp[2]
    C.name = join(newname,'.')

    # change rotated flag
    C.rotated = true

    # update correlation
    C.corr = new_corr
    return C
end


text:
#=
Emitters are defined by Pixels and SpatialLayouts. A emitter has a spectrum, and an optical power distribution over the hemisphere. These are
intrinsic physical properties of the emitter.
=#
#Ray distribution types
using Distributions

"""
Each AbstractSpectrum type defines a spectrumSample function which returns a uniformly sampled point from the spectrum of the light source and the power at that wavelength.
"""
abstract type AbstractSpectrum{T<:Real} end
""" Every ray direction distribution must implement the functions direction(a::Source)"""
abstract type AbstractRayDirectionDistribution{T<:Real} end
abstract type AbstractAngularPowerDistribution{T<:Real} end
abstract type AbstractRayOriginDistribution{T<:Real} end

"""
Create an emitter by choosing a type for the spectral distribution, ray origin distribution, ray direction distribution, and angular light power distribution.
For example, to create an emitter with uniform spectrum, rays originating on a regular grid, ray directions distributed around a vector v up to a maximum angle θmax, with Lambertian power distribution call the Source construction with these types and arguments:

a = Source{UniformSpectrum,GridOrigin, ConeDistribution, Lambertian}(UniformSpectrum(),GridOrigin(),ConeDistribution(θmax),Lambertian())
"""
struct Source{T<:Real,S<:AbstractSpectrum{T},O<:AbstractRayOriginDistribution{T},D<:AbstractRayDirectionDistribution{T},P<:AbstractAngularPowerDistribution{T}}
    spectrum::S
    rayorigin::O
    raydirection::D
    power::P

    Source(spectrum::S, origins::O, directions::D, power::P, ::Type{T} = Float64) where {S<:AbstractSpectrum,O<:AbstractRayOriginDistribution,D<:AbstractRayDirectionDistribution,P<:AbstractAngularPowerDistribution,T<:Real} = new{T,S,O,D,P}(spectrum, origins, directions, power)
end
export Source

const UNIFORMSHORT = 0.450 #um
const UNIFORMLONG = 0.680 #um

"""
flat spectrum from 450nm to 680nm
"""
struct UniformSpectrum{T} <: AbstractSpectrum{T} end
export UniformSpectrum
"""returns a tuple (power,wavelength) """
spectrumsample(::Source{T,UniformSpectrum{T},O,D,P}) where {O,D,P,T} = (T(1), rand(Distributions.Uniform(UNIFORMSHORT, UNIFORMLONG)))

struct DeltaFunctionSpectrum{T<:Real} <: AbstractSpectrum{T}
    λ::T
end
export DeltaFunctionSpectrum
spectrumsample(a::Source{T,DeltaFunctionSpectrum{T},O,D,P}) where {O,D,P,T} = (T(1), a.spectrum.λ)

"""
Use measured spectrum to compute emitter power. Create spectrum by reading CSV files.

Evaluate spectrum at arbitrary wavelength with [`spectrumpower`](@ref)
"""
struct MeasuredSpectrum{T<:Real} <: AbstractSpectrum{T}
    lowwavelength::Int
    highwavelength::Int
    wavelengthstep::Int
    powersamples::Vector{T}

    function MeasuredSpectrum(samples::DataFrame)
        colnames = names(samples)
        @assert "Wavelength" in colnames
        @assert "Power" in colnames
        wavelengths = samples[:Wavelength]
        @assert eltype(wavelengths) <: Integer

        power::Vector{T} where {T<:Real} = samples[:Power] # no missing values allowed and must be real numbers
        maxpower = maximum(power)
        T = eltype(power)
        p = sortperm(wavelengths)
        wavelengths = wavelengths[p]
        power = power[p] ./ maxpower #make sure power values have the same sorted permutation as wavelengths normalized with maximum power equal to 1
        λmin = wavelengths[p[1]]
        λmax = wavelengths[p[end]]
        step = wavelengths[2] - wavelengths[1]

        for i in 3:length(wavelengths)
            @assert wavelengths[i] - wavelengths[i - 1] == step  # no missing values allowed and step size must be consistent
        end

        return new{T}(λmin, λmax, step, power)
    end
end
export MeasuredSpectrum

"""expects wavelength in nm not um"""
function spectrumpower(spectrum::MeasuredSpectrum{T}, λ::T)::Union{Nothing,T} where {T<:Real}
    λ = λ #convert from um to nm
    if λ < spectrum.lowwavelength || λ > spectrum.highwavelength
        return nothing
    end

    lowindex = floor(Int, (λ - spectrum.lowwavelength)) ÷ spectrum.wavelengthstep + 1

    if lowindex == length(spectrum.powersamples)
        return T(spectrum.powersamples[end])
    else
        highindex = lowindex + 1
        α = mod(λ, spectrum.wavelengthstep) / spectrum.wavelengthstep

        return T((1 - α) * spectrum.powersamples[lowindex] + α * spectrum.powersamples[highindex])
    end
end

# Every spectrum type must define a function spectrumsample which randomly selects a wavelength and returns the power at that wavelength. For delta function spectra the same wavelength will be returned every time.

function spectrumsample(a::Source{T,MeasuredSpectrum{T},O,D,P}) where {O,D,P,T}
    spectrum = a.spectrum::MeasuredSpectrum{T}
    λ = rand(Distributions.Uniform{T}(T(spectrum.lowwavelength), T(spectrum.highwavelength)))
    power = spectrumpower(spectrum, λ)
    if power === nothing        #added this condition because the compiler was allocating if just returned values directly.
        return (nothing, nothing)
    else
        return (power, λ / T(1000)) #convert from nm to um
    end
end

"""stores the state of ray generation. Keeps track of ray number so you can have a PointOrigin which emits multiple rays. This can be used in the future to do phase summations on the image plane, since only point sources are perfectly spatially coherent."""
struct RayState
    originnumber::Int
    directionnumber::Int

    RayState() = new(1, 1)
    RayState(originnum::Int, directionnum::Int) = new(originnum, directionnum)
end

struct PointOrigin{T<:Real} <: AbstractRayOriginDistribution{T}
    position::SVector{3,T}
end
export PointOrigin

numsamples(a::PointOrigin{T}) where {T} = 1

origin(a::Source{T,S,PointOrigin{T},D,P}, ::Int) where {S,D,P,T} = a.rayorigin.position

struct UniformRandomOrigin{T} <: AbstractRayOriginDistribution{T}
    width::T
    height::T
    numsamples::Int
end
export UniformRandomOrigin

numsamples(a::UniformRandomOrigin{T}) where {T} = a.numsamples

function origin(a::Source{T,S,UniformRandomOrigin{T},D,P}, ::Int) where {S,D,P,T}
    origin = a.origin
    x = rand(Distributions.Uniform(-origin.width / 2.0, origin.width / 2.0))
    y = rand(Distributions.Uniform(-origin.height / 2.0, origin.height / 2.0))
    z = T(0)
end

struct GridOrigin{T} <: AbstractRayOriginDistribution{T}
    width::T
    height::T
    usamples::Int
    vsamples::Int
    ustep::T
    vstep::T

    GridOrigin(width::T, height::T, usamples::Int, vsamples::Int) where {T<:Real} = new{T}(width, height, usamples, vsamples, width / (usamples - 1), height / (vsamples - 1))
end
export GridOrigin

numsamples(a::GridOrigin{T}) where {T<:Real} = a.usamples * a.vsamples

function origin(a::Source{T,S,GridOrigin{T},D,P}, originnumber::Int) where {S,D,P,T}
    @assert originnumber <= numsamples(a.rayorigin)
    origin = a.rayorigin

    u = originnumber ÷ origin.vsamples
    v = mod(originnumber, origin.usamples)
    return SVector{3,T}(T(-origin.width / 2.0 + origin.ustep * u), T(-origin.height / 2.0 + origin.vstep * v), T(0))
end

struct HexapolarOrigin{T} <: AbstractRayOriginDistribution{T}
    position::SVector{3,T}
    direction::SVector{3,T}
    uvec::SVector{3,T}
    vvec::SVector{3,T}
    halfsizeu::T
    halfsizev::T
    nrings::Int
    function HexapolarOrigin{T}(nrings::Int, halfsizeu::T, halfsizev::T; position::SVector{3,T} = SVector{3,T}(0.0, 0.0, 0.0), direction::SVector{3,T} = SVector{3,T}(0.0, 0.0, -1.0), rotationvec::SVector{3,T} = SVector{3,T}(0.0, 1.0, 0.0)) where {T<:Real}
        new{T}(position, getuvvecs(direction, rotationvec)..., halfsizeu, halfsizev, nrings)
    end
end

struct ConstantRayDirection{T} <: AbstractRayDirectionDistribution{T}
    direction::SVector{3,T}
end
export ConstantRayDirection

numsamples(a::ConstantRayDirection{T}) where {T} = 1

direction(a::Source{T,S,O,ConstantRayDirection{T},P}, directionnumber::Int) where {S,O,P,T} = a.raydirection.direction

struct ConeDistribution{T} <: AbstractRayDirectionDistribution{T}
    θmax::T
    numsamples::Int
end
export ConeDistribution

numsamples(a::ConeDistribution{T}) where {T} = a.numsamples

"""
Generates a unit vector pointing somewhere within the cone with half angle `θmax` around `direction` which is the normal to the emitter surface.
"""
function direction(a::Source{T,S,O,ConeDistribution{T},P}, directionnumber::Int) where {S,O,P,T<:Real}
    direction = SVector{3,T}(T(0), T(0), T(1))
    θmax = a.raydirection.θmax
    uvec = SVector{3,T}(T(1), T(0), T(0))
    vvec = SVector{3,T}(T(0), T(1), T(0))

    ϕ = rand(T) * 2π
    θ = NaNsafeacos(one(T) + rand(T) * (cos(θmax) - 1))
    return normalize(sin(θ) * (cos(ϕ) * uvec + sin(ϕ) * vvec) + cos(θ) * direction)
end

Base.length(a::Source{T,S,O,D,P}) where {T,S,O,D,P} = numsamples(a.rayorigin) * numsamples(a.raydirection)

"""raysample is the function that can couple origin and direction generation if necessary. The default function couples them in a simple way but more complex coupling should be possible. For each origin numsamples(a.raydirection) direction samples are taken with identical origin. Then the origin number is incremented. This repeats till all rays have been generated. The origin and direction functions receive an integer indicating the origin or direction number so regular patterns such as rectangular and hexapolar grids can be generated properly."""
function raysample(a::Source{T,S,O,D,P}, r::RayState) where {T,S,O,D,P}
    newray = Ray(origin(a, r.originnumber), direction(a, r.directionnumber))
    dirnum = r.directionnumber + 1
    if dirnum > numsamples(a.raydirection)
        return (newray, RayState(r.originnumber + 1, 1))  #reset counter for direction number and move on to next origin
    else
        return (newray, RayState(r.originnumber, dirnum))
    end

end

struct Lambertian{T} <: AbstractAngularPowerDistribution{T} end
export Lambertian

directionpower(::Source{T,S,O,D,Lambertian{T}}, ::SVector{3,T}) where {S,O,D,T} = T(1)

function Base.iterate(a::Source)::Tuple{OpticalRay,RayState}
    state = RayState()
    generateray(a, state)
end

function Base.iterate(a::Source, b::RayState)::Union{Nothing,Tuple{OpticalRay,RayState}}
    if b.originnumber > numsamples(a.rayorigin)
        return nothing
    else
        return generateray(a, b)
    end
end

# """generates an optical ray with origin and spectrum distribution dictated by OriginDistribution and Spectrum type. Ray is generated in the emitter local coordinate frame."""
function generateray(a::Source, state::RayState, sourcenumber = 0)
    spectralpower, λ = spectrumsample(a)
    ray, newstate = raysample(a, state)
    optray = OpticalRay(ray, spectralpower * directionpower(a, direction(ray)), λ, sourcenum = sourcenumber)
    return (optray, newstate)
end


text:
using ZeroIndex

function count_change(mat::Array{Int64, 2},
                      denom::Int64,
                      row::Int64,
                      col::Int64)::Int64
    prev_col = col - denom
    prev_row = row - 1
    denom_in = 0
    denom_out = 0
    if prev_row >= 0
        denom_out = aget(mat, prev_row, col)
    end
    if prev_col >= 0
        denom_in = aget(mat, row, prev_col)
    end
    return denom_in + denom_out
end

function main()
    denoms = [1, 2, 5, 10, 20, 50, 100, 200]
    sought = 200
    ndenoms = length(denoms)
    nways = Array{Int64, 2}(ndenoms, sought + 1)
    for row = 0:(ndenoms-1)
        aset!(nways, row, 0, 1)
    end
    for col in 1:sought
        for row in 0:(ndenoms-1)
            denom = vget(denoms, row)
            val = count_change(nways, denom, row, col)
            aset!(nways, row, col, val)
        end
    end
    @printf("number of ways is = %d\n", aget(nways, ndenoms-1, sought))
end

main()


text:
@testset "Test Sampling" begin
    @info "Test: Sampling Points on Ellipse"

    @testset "Samples points satisfy ellipse equation" begin
        𝛏 = [10, 5, 25, 25, π/4]
        geo_to_alg = GeometricToAlgebraic()
        𝛉 = geo_to_alg(𝛏)
        α₁ = 0
        α₂ = 2*π
        N = 25
        observations = sample_points_on_ellipse(𝛏..., N, α₁, α₂)
        𝐫 = GuaranteedEllipseFitting.vector_valued_objective(observations, 𝛉)
        @test isapprox(dot(𝐫, 𝐫), 0.0, atol= 1e-10)
    end

end


text:
# Experimental functions for linear extensions

export linear_extension, linear_extension_empty!, all_linear_extensions, random_linear_extension

"""
`linear_extension(P)` returns a linear extension of the poset `P`.
"""
function linear_extension(P::SimplePoset{T}) where T
    PP = deepcopy(P)
    result = linear_extension_empty!(PP)
    return result
end

"""
`linear_extension_empty!(P)` returns a linear extension of the poset `P`,
but also empties the poset `P`. This is a destructive operation.
"""
function linear_extension_empty!(PP::SimplePoset{T}) where T
    result = T[]
    while card(PP)>0
        M = minimals(PP)
        append!(result, M)
        for x in M
            delete!(PP,x)
        end
    end
    return result
end

"""
`random_linear_extension(P)` generates a linear extension of `P`
at random. This is slower than `linear_extension(P)` if one just
wants some linear extension.

Note that linear extensions are not generated uniformly
at random by this function. Rather we proceed recursively: A
minimal element `x` of `P` is chosen uniformly at random from among
all minimal elements, and that becomes the first (lowest) element in
the linear extension. Then `x` is deleted and the same process is
repeated until all elements have been processed.
"""
function random_linear_extension(P::SimplePoset{T}) where T
    result = T[]
    PP = deepcopy(P)
    while card(PP) > 0
        x = rand(minimals(PP))  # choose a min el't at random
        push!(result,x)
        delete!(PP,x)
    end
    return result
end

_LX_table = Dict{SimplePoset, Set}()
export clear_LX_table

"""
`clear_LX_table()` releases cached results computed by
`all_linear_extensions`.
"""
function clear_LX_table()
    global _LX_table = Dict{SimplePoset, Set}()
    nothing
end

"""
`all_linear_extensions(P)` returns the `Set` of all linear extensions
of `P`. This can take a very long time and eat up a lot of memory
(which can be freed using `clear_LX_table`).
"""
function all_linear_extensions(P::SimplePoset)::Set
    # see if we already have an answer
    global _LX_table
    if haskey(_LX_table,P)
        return _LX_table[P]
    end

    T = eltype(P)
    result::Set = Set{Array{T,1}}()
    if card(P) == 0
        return result
    end
    if card(P) == 1
        L = elements(P)
        push!(result, L)
        return result
    end

    M = maximals(P)
    for x in M
        PP = deepcopy(P)
        delete!(PP,x)
        PP_exts = all_linear_extensions(PP)
        for L in PP_exts
            LL = deepcopy(L)
            append!(LL,[x])
            push!(result, LL)
        end
    end
    _LX_table[P] = result   # save this answer
    return result
end


text:
module PatternFolds

# Reexporr
using Reexport

# Import
@reexport using Intervals
using Lazy
using Random

# exports
export IntervalsFold
export IVectorFold
export VectorFold

export check_pattern
export fold
export folds
export gap
export length # TODO: rename span
export make_vector_fold
export pattern
export unfold

# includes
include("intervals.jl")
include("common.jl")
include("immutable_vector.jl")
include("vector.jl")

end