text:
"""
    spectral_mixture_kernel(
        h::Kernel=SqExponentialKernel(),
        αs::AbstractVector{<:Real},
        γs::AbstractMatrix{<:Real},
        ωs::AbstractMatrix{<:Real},
    )

where αs are the weights of dimension (A, ), γs is the covariance matrix of
dimension (D, A) and ωs are the mean vectors and is of dimension (D, A).
Here, D is input dimension and A is the number of spectral components.

`h` is the kernel, which defaults to [`SqExponentialKernel`](@ref) if not specified.

Generalised Spectral Mixture kernel function. This family of functions is  dense
in the family of stationary real-valued kernels with respect to the pointwise convergence.[1]

```math
   κ(x, y) = αs' (h(-(γs' * t)^2) .* cos(π * ωs' * t), t = x - y
```

# References:
    [1] Generalized Spectral Kernels, by Yves-Laurent Kom Samo and Stephen J. Roberts
    [2] SM: Gaussian Process Kernels for Pattern Discovery and Extrapolation,
            ICML, 2013, by Andrew Gordon Wilson and Ryan Prescott Adams,
    [3] Covariance kernels for fast automatic pattern discovery and extrapolation
        with Gaussian processes, Andrew Gordon Wilson, PhD Thesis, January 2014.
        http://www.cs.cmu.edu/~andrewgw/andrewgwthesis.pdf
    [4] http://www.cs.cmu.edu/~andrewgw/pattern/.

"""
function spectral_mixture_kernel(
    h::Kernel,
    αs::AbstractVector{<:Real},
    γs::AbstractMatrix{<:Real},
    ωs::AbstractMatrix{<:Real},
)
    if !(size(αs, 1) == size(γs, 2) == size(ωs, 2))
        throw(DimensionMismatch("The dimensions of αs, γs, ans ωs do not match"))
    end
    if size(γs) != size(ωs)
        throw(DimensionMismatch("The dimensions of γs ans ωs do not match"))
    end

    return sum(zip(αs, eachcol(γs), eachcol(ωs))) do (α, γ, ω)
        a = TransformedKernel(h, LinearTransform(γ'))
        b = TransformedKernel(CosineKernel(), LinearTransform(ω'))
        return α * a * b
    end
end

function spectral_mixture_kernel(
    αs::AbstractVector{<:Real},
    γs::AbstractMatrix{<:Real},
    ωs::AbstractMatrix{<:Real}
)
    return spectral_mixture_kernel(SqExponentialKernel(), αs, γs, ωs)
end

"""
    spectral_mixture_product_kernel(
        h::Kernel=SqExponentialKernel(),
        αs::AbstractMatrix{<:Real},
        γs::AbstractMatrix{<:Real},
        ωs::AbstractMatrix{<:Real},
    )

where αs are the weights of dimension (D, A), γs is the covariance matrix of
dimension (D, A) and ωs are the mean vectors and is of dimension (D, A).
Here, D is input dimension and A is the number of spectral components.

Spectral Mixture Product Kernel. With enough components A, the SMP kernel
can model any product kernel to arbitrary precision, and is flexible even
with a small number of components [1]

`h` is the kernel, which defaults to [`SqExponentialKernel`](@ref) if not specified.

```math
   κ(x, y) = Πᵢ₌₁ᴷ Σ(αsᵢᵀ .* (h(-(γsᵢᵀ * tᵢ)²) .* cos(ωsᵢᵀ * tᵢ))), tᵢ = xᵢ - yᵢ
```

# References:
    [1] GPatt: Fast Multidimensional Pattern Extrapolation with GPs,
        arXiv 1310.5288, 2013, by Andrew Gordon Wilson, Elad Gilboa,
        Arye Nehorai and John P. Cunningham
"""
function spectral_mixture_product_kernel(
    h::Kernel,
    αs::AbstractMatrix{<:Real},
    γs::AbstractMatrix{<:Real},
    ωs::AbstractMatrix{<:Real},
)
    if !(size(αs) == size(γs) == size(ωs))
        throw(DimensionMismatch("The dimensions of αs, γs, ans ωs do not match"))
    end
    return TensorProduct(spectral_mixture_kernel(h, α, reshape(γ, 1, :), reshape(ω, 1, :))
        for (α, γ, ω) in zip(eachrow(αs), eachrow(γs), eachrow(ωs)))
end

function spectral_mixture_product_kernel(
    αs::AbstractMatrix{<:Real},
    γs::AbstractMatrix{<:Real},
    ωs::AbstractMatrix{<:Real}
)
    return spectral_mixture_product_kernel(SqExponentialKernel(), αs, γs, ωs)
end


text:
using Test
using JSExpr
using JSExpr: @crawl, deparse, JSTerminal, JSAST

# myjs = @js foo = nothing
# @test myjs.s == "foo = null"

@testset "special symbols" begin
    @testset "nothing" begin
        jsast = @crawl(foo = nothing)
        @test jsast.head == :(=)
        @test length(jsast.args) == 2
        @test jsast.args[1] == JSTerminal(:foo)
        @test jsast.args[2] == JSTerminal(:null)

        jsstr = deparse(jsast)
        @test jsstr == js"foo = null"
    end

    @testset "booleans" begin
        jsast = @crawl(foo = true)
        @test jsast.head == :(=)
        @test length(jsast.args) == 2
        @test jsast.args[1] == JSTerminal(:foo)
        @test jsast.args[2] == JSTerminal(Symbol(:true))

        jsstr = deparse(jsast)
        @test jsstr == js"foo = true"

        jsast = @crawl(foo = false)
        @test jsast.head == :(=)
        @test length(jsast.args) == 2
        @test jsast.args[1] == JSTerminal(:foo)
        @test jsast.args[2] == JSTerminal(Symbol(:false))

        jsstr = deparse(jsast)
        @test jsstr == js"foo = false"
    end
end

@testset "number literals" begin
    jsast = @crawl(foo = 123)
    @test jsast.head == :(=)
    @test length(jsast.args) == 2
    @test jsast.args[1] == JSTerminal(:foo)
    @test jsast.args[2] == JSTerminal(123)

    jsstr = deparse(jsast)
    @test jsstr == js"foo = 123"
end

@testset "string literals" begin
    jsast = @crawl(foo = "foo")
    @test jsast.head == :(=)
    @test length(jsast.args) == 2
    @test jsast.args[1] == JSTerminal(:foo)
    @test jsast.args[2] == JSTerminal(js"\"foo\"")

    jsstr = deparse(jsast)
    @test jsstr == js"foo = \"foo\""

    @test string(@js """foo is "fantastic"!""") == "\"foo is \\\"fantastic\\\"!\""
end

ext: jl


text:
@testset "fsm_stream_reassembler_dup.jl" begin

    @testset begin
        rs = StreamReassembler(65000)
        push_substring!(rs, "abcd", 0)
        @test rs.cur_idx == 4
        @test read!(rs.output) == "abcd"
        @test !eof(rs.output)

        push_substring!(rs, "abcd", 0)
        @test rs.cur_idx == 4
        @test read!(rs.output) == ""
        @test !eof(rs.output)
    end
    @testset begin
        rs = StreamReassembler(65000)
        push_substring!(rs, "abcd", 0)
        @test rs.cur_idx == 4
        @test read!(rs.output) == "abcd"
        @test !eof(rs.output)

        push_substring!(rs, "abcd", 4)
        @test rs.cur_idx == 8
        @test read!(rs.output) == "abcd"
        @test !eof(rs.output)

        push_substring!(rs, "abcd", 0)
        @test rs.cur_idx == 8
        @test read!(rs.output) == ""
        @test !eof(rs.output)

        push_substring!(rs, "abcd", 4)
        @test rs.cur_idx == 8
        @test read!(rs.output) == ""
        @test !eof(rs.output)
    end

    @testset begin
        rs = StreamReassembler(65000)
        push_substring!(rs, "abcdefgh", 0)
        @test rs.cur_idx == 8
        @test read!(rs.output) == "abcdefgh"
        @test !eof(rs.output)
        data = "abcdefgh"
        for _ in 1:1000
            st = rand(0:8)
            ed = rand(st:8)
            push_substring!(rs, @view(data[st+1:ed]), st)
            @test rs.cur_idx == 8
            @test read!(rs.output) == ""
            @test !eof(rs.output)
        end
    end

    @testset begin
        rs = StreamReassembler(65000)
        push_substring!(rs, "abcd", 0)
        @test rs.cur_idx == 4
        @test read!(rs.output) == "abcd"
        @test !eof(rs.output)

        push_substring!(rs, "abcdef", 0)
        @test rs.cur_idx == 6
        @test read!(rs.output) == "ef"
        @test !eof(rs.output)
    end

end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
############################################################################################
#                              OBJECTIVES                                                  #
############################################################################################

abstract type AbstractObjective end
Base.length(obj::AbstractObjective) = length(obj.cost)
state_dim(obj::AbstractObjective) = throw(ErrorException("state_dim not implemented"))
control_dim(obj::AbstractObjective) = throw(ErrorException("control_dim not implemented"))
get_J(obj::AbstractObjective) = throw(ErrorException("get_J not implemented"))

"""$(TYPEDEF)
Objective: stores stage cost(s) and terminal cost functions

Constructors:
```julia
Objective(cost, N)
Objective(cost, cost_term, N)
Objective(costs::Vector{<:CostFunction}, cost_term)
Objective(costs::Vector{<:CostFunction})
```
"""
struct Objective{C} <: AbstractObjective
    cost::Vector{C}
    J::Vector{Float64}
    const_grad::BitVector
    const_hess::BitVector
    function Objective(cost::Vector{C}) where C <: CostFunction
        N = length(cost)
        J = zeros(N)
        grad = zeros(Bool,N)
        hess = zeros(Bool,N)
        new{C}(cost, J, grad, hess)
    end
end

state_dim(obj::Objective) = state_dim(obj.cost[1])
control_dim(obj::Objective) = control_dim(obj.cost[1])
Base.size(obj::Objective) = (state_dim(obj), control_dim(obj))
@inline ExpansionCache(obj::Objective) = ExpansionCache(obj[1])

"""
    is_quadratic(obj::Objective)

Only valid for a cost expansion, i.e. an objective containing the 2nd order expansion of
    another objective. Determines if the original objective is a quadratic function, or
    in other words, if the hessian of the objective is constant.

For example, if the original cost function is an augmented Lagrangian cost function, the
    result will return true only if all constraints are linear.
"""
is_quadratic(obj::Objective) = all(obj.const_hess)

# Constructors
function Objective(cost::CostFunction,N::Int)
    Objective([cost for k = 1:N])
end

function Objective(cost::CostFunction, cost_terminal::CostFunction, N::Int)
    stage, term = promote(cost, cost_terminal)
    Objective([k < N ? stage : term for k = 1:N])
end

function Objective(cost::Vector{<:CostFunction},cost_terminal::CostFunction)
    N = length(cost) + 1
    Objective([cost...,cost_terminal])
end

# Methods
"Get the vector of costs at each knot point. `sum(get_J(obj))` is equal to the cost"
get_J(obj::Objective) = obj.J

Base.copy(obj::Objective) = Objective(copy.(obj.cost))

Base.getindex(obj::Objective,i::Int) = obj.cost[i]

@inline Base.firstindex(obj::Objective) = firstindex(obj.cost)
@inline Base.lastindex(obj::Objective) = lastindex(obj.cost)
@inline Base.iterate(obj::Objective, start=1) = Base.iterate(obj.cost, start)
@inline Base.eltype(obj::Objective) = eltype(obj.cost)
@inline Base.length(obj::Objective) = length(obj.cost)
Base.IteratorSize(obj::Objective) = Base.HasLength()
Base.eachindex(obj::Objective) = Base.OneTo(length(obj))

Base.show(io::IO, obj::Objective{C}) where C = print(io,"Objective")

"""
    CostExpansion{n,m,T}

A vector of `Expansion`s, combined with some bit vectors for storing whether the
terms are constant or not.

# Constructors
    CostExpansion{T}(n,m,N)
    CostExpansion(n,m,T)  # defaults to Float64
    CostExpansion(E::CostExpansion, model)

where the last constructor is for allocating the expansion for the "raw" expansion,
    before accounting for the Lie group structure. For a standard (Euclidiean) state
    vector, it will create an alias to the original expansion. For a `LieGroupModel`,
    it will allocate new storage.
"""
struct CostExpansion{n,m,T} <: AbstractArray{Expansion{n,m,T},1}
    data::Vector{Expansion{n,m,T}}
    const_hess::BitVector
    const_grad::BitVector
    function CostExpansion{T}(n::Int, m::Int, N::Int) where T
        data = [Expansion{T}(n,m) for k = 1:N]
        const_hess = BitVector(undef, N)
        const_grad = BitVector(undef, N)
        new{n,m,T}(data, const_hess, const_grad)
    end
end
@inline CostExpansion(n,m,N) = CostExpansion{Float64}(n,m,N)
function CostExpansion(E::CostExpansion, model::AbstractModel)
    # Create QuadraticObjective linked to error cost expansion
    @assert RobotDynamics.state_diff_size(model) == size(model)[1]
    return E
end

function CostExpansion(E::CostExpansion{n,m,T}, model::LieGroupModel) where {n,m,T}
    # Create an expansion for the full state dimension
    @assert length(E[1].q) == RobotDynamics.state_diff_size(model)
    n0 = state_dim(model)
    return CostExpansion{T}(n0,m,length(E))
end
Base.size(E::CostExpansion) = size(E.data)
Base.getindex(E::CostExpansion, i::Int) = Base.getindex(E.data, i)

# ############################################################################################
# #                            Quadratic Objectives (Expansions)
# ############################################################################################
# const QuadraticObjective{n,m,T} = Objective{QuadraticCost{n,m,T,SizedMatrix{n,n,T,2,Matrix{T}},SizedMatrix{m,m,T,2,Matrix{T}}}}
# const QuadraticExpansion{n,m,T} = Objective{<:QuadraticCostFunction{n,m,T}}
# const DiagonalCostFunction{n,m,T} = Union{DiagonalCost{n,m,T},QuadraticCost{n,m,T,<:Diagonal,<:Diagonal}}

# function QuadraticObjective(n::Int, m::Int, N::Int, isequal::Bool=false)
#     Objective([QuadraticCost{Float64}(n,m, terminal=(k==N) && !isequal) for k = 1:N])
# end

# function QuadraticObjective(obj::QuadraticObjective, model::AbstractModel)
#     # Create QuadraticObjective linked to error cost expansion
#     @assert RobotDynamics.state_diff_size(model) == size(model)[1]
#     return obj
# end

# function QuadraticObjective(obj::QuadraticObjective, model::LieGroupModel)
#     # Create QuadraticObjective partially linked to error cost expansion
#     @assert length(obj[1].q) == RobotDynamics.state_diff_size(model)
#     n,m = size(model)
#     costfuns = map(obj.cost) do costfun
#         Q = SizedMatrix{n,n}(zeros(n,n))
#         R = costfun.R
#         H = SizedMatrix{m,n}(zeros(m,n))
#         q = @MVector zeros(n)
#         r = costfun.r
#         c = costfun.c
#         QuadraticCost(Q,R,H,q,r,c, checks=false, terminal=costfun.terminal)
#     end
#     Objective(costfuns)
# end

# Convenience constructors
@doc raw"""
    LQRObjective(Q, R, Qf, xf, N)

Create an objective of the form
`` (x_N - x_f)^T Q_f (x_N - x_f) + \sum_{k=0}^{N-1} (x_k-x_f)^T Q (x_k-x_f) + u_k^T R u_k``

Where `eltype(obj) <: DiagonalCost` if `Q`, `R`, and `Qf` are
    `Union{Diagonal{<:Any,<:StaticVector}}, <:StaticVector}`
"""
function LQRObjective(Q::AbstractArray, R::AbstractArray, Qf::AbstractArray,
        xf::AbstractVector, N::Int; checks=true, uf=@SVector zeros(size(R,1)))
    @assert size(Q,1) == length(xf)
    @assert size(Qf,1) == length(xf)
    @assert size(R,1) == length(uf)
    n = size(Q,1)
    m = size(R,1)
    H = SizedMatrix{m,n}(zeros(m,n))
    q = -Q*xf
    r = -R*uf
    c = 0.5*xf'*Q*xf + 0.5*uf'R*uf
    qf = -Qf*xf
    cf = 0.5*xf'*Qf*xf

    ℓ = QuadraticCost(Q, R, H, q, r, c, checks=checks)
    ℓN = QuadraticCost(Qf, R, H, q, r, cf, checks=false, terminal=true)

    Objective(ℓ, ℓN, N)
end

function LQRObjective(
        Q::Union{<:Diagonal, <:AbstractVector},
        R::Union{<:Diagonal, <:AbstractVector},
        Qf::Union{<:Diagonal, <:AbstractVector},
        xf::AbstractVector, N::Int;
        uf=(@SVector zeros(size(R,1))),
        checks=true)
    n,m = size(Q,1), size(R,1)
    @assert size(Q,1) == length(xf)
    @assert size(Qf,1) == length(xf)
    @assert size(R,1) == length(uf)
    Q,R,Qf = Diagonal(Q), Diagonal(R), Diagonal(Qf)
    q = -Q*xf
    r = -R*uf
    c = 0.5*xf'*Q*xf + 0.5*uf'R*uf
    qf = -Qf*xf
    cf = 0.5*xf'*Qf*xf

    ℓ = DiagonalCost(Q, R, q, r, c, checks=checks, terminal=false)

    ℓN = DiagonalCost(Qf, R, qf, r, cf, checks=false, terminal=true)

    Objective(ℓ, ℓN, N)
end

"""
    TrackingObjective(Q, R, Z; [Qf])

Generate a quadratic objective that tracks the reference trajectory specified by `Z`.
"""
function TrackingObjective(Q,R,Z::AbstractTrajectory; Qf=Q)
    costs = map(Z) do z
        LQRCost(Q, R, state(z), control(z))
    end
    costs[end] = LQRCost(Qf, R, state(Z[end]))
    Objective(costs)
end

"""
    update_trajectory!(obj, Z, [start=1])

For use with a tracking-style trajectory (see [`TrackingObjective`](@ref)).
Update the costs to track the new trajectory `Z`. The `start` parameter specifies the
index of reference trajectory that should be used as the starting point of the reference
tracked by the objective. This is useful when a single, long time-horizon trajectory is given
but the optimization only tracks a portion of the reference at each solve (e.g. MPC).
"""
function update_trajectory!(obj::Objective{<:QuadraticCostFunction}, Z::AbstractTrajectory, start=1)
    inds = (start-1) .+ (1:length(obj))
    for (i,k) in enumerate(inds)
        set_LQR_goal!(obj[i], state(Z[k]), control(Z[k]))
    end
end


text:
#. Abstract type for coefficients with a special type in which we perform the addmul operation
# Coefficients are stored in T, while addmul is done in Tbuf

"""
Expect:
    - Base.isone(::T) :: Bool
    - Base.iszero(::T) :: Bool
    - add(::NmodLikeΓ{T}, ::T, ::T) :: T
    - mul(::NmodLikeΓ{T}, ::T, ::T) :: T
    - inv(::NmodLikeΓ{T}, ::T) :: T
    - submul(::NmodLikeΓ{T, Tbuf}, ::Tbuf, ::T, ::T) :: Tbuf
    - normal(::NmodLikeΓ{T, Tbuf}, ::Tbuf) :: Tbuf
    - inflate(::NmodLikeΓ{T, Tbuf}, ::T) :: Tbuf
    - deflate(::NmodLikeΓ{T, Tbuf}, ::Tbuf) :: T
"""
abstract type NmodLikeΓ{T, Tbuf}<:Context{T} end

buftype(::NmodLikeΓ{T, Tbuf}) where {T, Tbuf} = Tbuf

inflate(::NmodLikeΓ{T, Tbuf}, x :: T) where {T, Tbuf} = Tbuf(x)
deflate(::NmodLikeΓ{T, Tbuf}, x :: Tbuf) where {T, Tbuf} = x%T

#. Z/nZ

#.. Abstract type

abstract type NmodΓ{T<:Unsigned, Tbuf<:Unsigned}<:NmodLikeΓ{T,Tbuf} end

(ctx::NmodΓ)(a::AA.ResElem) = ctx(a.data)
rand_elem(ctx::NmodΓ{T}) where T = rand(one(T):ctx.char - one(T))
abstractalgebra(ctx :: NmodΓ) = Nemo.NmodRing(UInt64(ctx.char))

#.. 32-bit modulus and 64-bit buffer

mutable struct Nmod32Γ<:NmodΓ{UInt32, UInt64}
    char::UInt32
    maxshift::UInt64

    function Nmod32Γ(char)
        uchar = UInt64(char)
        nbits = 64 - leading_zeros(uchar)
        @assert 2*nbits <= 64
        if 2*char > typemax(UInt32)
            error("Addition won't work properly in this characteristic. Choose a higher bitcount for the field element type.")
        end
        new(uchar, uchar << leading_zeros(uchar))
    end
end

#... arithmetic operations
# Must implement inv, mul, addmul and normal

function add(ctx::Nmod32Γ, a::UInt32, b::UInt32)
    c0 = a+b
    # c1 = c0 - ctx.char
    # d = max(c0, c1)
    #@assert d%ctx.char == (a%ctx.char + b%ctx.char)%ctx.char
    return c0 % ctx.char
end

@inline function addmul(ctx::Nmod32Γ, a::UInt64, b::UInt32, c::UInt32)
    z0 = a + UInt64(b)*UInt64(c)
    z1 = z0 - ctx.maxshift
    z = (z0 < a) ? z1 : z0
    return z
end

@inline function submul(ctx::Nmod32Γ, a::UInt64, b::UInt32, c::UInt32)
    z0 = a - UInt64(b)*UInt64(c)
    z1 = z0 + ctx.maxshift
    z = (z0 > a) ? z1 : z0
    # if (z%ctx.char + (UInt64(b)%ctx.char)*(UInt64(c)%ctx.char)%ctx.char)%ctx.char != a%ctx.char
    #    error((ctx, a, b, c))
    # end
    return z
end

function opp(ctx::Nmod32Γ, a::UInt32)
    @assert ctx.char > a
    return ctx.char - a
end

function inv(ctx::Nmod32Γ, a::UInt32)
    invmod(UInt64(a), UInt64(ctx.char)) % UInt32
end

function mul(ctx::Nmod32Γ, a, b)
    UInt32((UInt64(a)*UInt64(b) % ctx.char))
end

function normal(ctx::Nmod32Γ, a::UInt64)
    a % ctx.char
end

#.. 32-bit modulus with 128-bits buffer

mutable struct Nmod32xΓ<:NmodΓ{UInt32, UInt128}
    char::UInt32
    twotothe64::UInt32        # 2^64 mod char

    function Nmod32xΓ(char)
        uchar = UInt32(char)
        p = (UInt64(2^32) % uchar)^2 % uchar % UInt32

        new(char, p)
    end
end

#... arithmetic operations
# Must implement inv, mul, submul and normal

function add(ctx::Nmod32xΓ, a::UInt32, b::UInt32)
    c0 = a+b
    c1 = c0 - ctx.char
    d = max(c0, c1)
    #@assert d%ctx.char == (a%ctx.char + b%ctx.char)%ctx.char
    return d
end

@inline function addmul(ctx::Nmod32xΓ, a::UInt128, b::UInt32, c::UInt32)
    z0 = a + UInt128(UInt64(b)*UInt64(c))
    return z0
end

@inline function submul(ctx::Nmod32xΓ, a::UInt128, b::UInt32, c::UInt32)
    z0 = a - UInt128(UInt64(b)*UInt64(c))
    return z0
end

function inv(ctx::Nmod32xΓ, a::UInt32)
    UInt32(invmod(UInt64(a), UInt64(ctx.char)))
end

function mul(ctx::Nmod32xΓ, a::UInt32, b::UInt32)
    UInt32((UInt64(a)*UInt64(b) % ctx.char) % UInt32)
end

function normal(ctx::Nmod32xΓ, a::UInt128)
    return a % ctx.char
    # hi = (a >> 64) % UInt64
    # lo = a % UInt64
    # add(ctx,
    #     ctx.twotothe64 * hi % ctx.char % UInt32,
    #     lo % ctx.char % UInt32)
end


text:
module SignalTemps

using CUDA, Distributed
using LinearAlgebra
using DataFrames,Query
using StatsBase, Distributions, Random
using UMAP
using LightGraphs,MetaGraphs

export getalpha, getbeta, getkappa, gettemps, safediv, getup, getupavg, groupavg, getgroups
export filterdf,filteravg,filtersum,umapexpr,sublabels,fplot,filterexpressed
export lrtManual

include("temps.jl") # calculate interactions

include("data.jl") # load in the data

include("lrtManual.jl") # manual lrt ident, uses OrthoDB ortholog (eg danio) search from PCquery given primary (eg human) pathway genes
end


text:
module DP

using JuMP
using Gurobi

include("queries.jl")

""" Returns the closest nondecreasing vector to data """
function isotonic(data::Vector; nonneg=false)
    model = Model(Gurobi.Optimizer)
    set_silent(model)
    n = length(data)
    @variable(model, x[1:n])
    if nonneg
        @constraint(model, x .>= 0)
    end
    @constraint(model, x[2:end] .>= x[1:end-1])
    @objective(model, Min, sum((x - data).^2))
    optimize!(model)
    value.(x)
end

""" Solves a least square problem to fit the queries
assumes data is some array and datashape is its shape
the info variable is a vector of tuples of query, noisy query answer, and weight
to use in least squares optimization.
"""

function fit(datashape, info::Array{Tuple{Q, A, F}}; nonneg=false) where {Q <: Query, A <: Array, F <: Number}
    model = Model(Gurobi.Optimizer)
    set_silent(model)
    numvars = reduce(*, datashape)
    @variable(model, x[1:numvars])
    mainvars = reshape(x, datashape)
    obj = 0
    for (query, noisyanswers, weight) in info
        numanswers = reduce(*, size(noisyanswers))
        tmpvars = @variable(model, [1:numanswers])
        newvars = reshape(tmpvars, size(noisyanswers))
        obj += sum((newvars - noisyanswers).^2)*weight
        @constraint(model, answer(query,mainvars) .== newvars)
    end
    @objective(model, Min, obj)
    optimize!(model)
    reshape(value.(x), datashape)
end

function mwem(queries, T, data, epsilon)
    datashape = size(data)
    synthdata = zeros(datashape...)
    sensitivity = maximum([sens(q) for q in queries])
    info = Nothing
    for i in 1:T
        epsilon_round = epsilon/T
        epsilon_select = epsilon_round/2
        epsilon_measure = epsilon_round/2
        errors = [sum(abs, answer(q, data) - answer(q, synthdata)) for q in queries]
        worst_query_id = noisy_max(errors, sensitivity=sensitivity, epsilon=epsilon_select)
        worst_query = queries[worst_query_id]
        worst_query_answer = protect(worst_query, data, epsilon_measure)
        if info == Nothing
            info = [(worst_query, worst_query_answer, 1.0)]
        else
            info = cat(info, (worst_query, worst_query_answer, 1.0), dims=1)
        end
        synthdata = fit(datashape, info, nonneg=true)
    end
    synthdata
end

end # end module

max_stars_repo_name: dkifer/dpcourse

max_issues_repo_name: dkifer/dpcourse

max_forks_repo_name: dkifer/dpcourse


text:
# FileDataset

struct FileDataset
    tree::FileTree
    nodes::Vector{FileTrees.File}
end

function FileDataset(args...; kwargs...)
    tree = FileTree(args...; kwargs...)
    return FileDataset(tree, FileTrees.files(tree))
end

Base.show(io::IO, data::FileDataset) = print(
    io,
    "FileDataset(\"", data.tree.name, "\", ", nobs(data), " observations)")

LearnBase.nobs(ds::FileDataset) = length(ds.nodes)
LearnBase.getobs(ds::FileDataset, idx::Int) = path(ds.nodes[idx])

# File utilities

"""
    loadfile(file)

Load a file from disk into the appropriate format.
"""
function loadfile(file::String)
    if isimagefile(file)
        # faster image loading
        return FileIO.load(file, view = true)
    else
        return FileIO.load(file)
    end
end

loadfile(file::AbstractPath) = loadfile(string(file))
loadfile(file::FileTrees.File) = loadfile(path(file))

isimagefile(file::AbstractPath) = isimagefile(string(file))
isimagefile(file::File) = isimagefile(file.name)
isimagefile(file::String) = occursin(IMAGEFILE_REGEX, lowercase(file))
const IMAGEFILE_REGEX = r"\.(gif|jpe?g|tiff?|png|webp|bmp)$"

## TODO: TableDataset


text:
using GalacticOptim
using Test

@testset "GalacticOptim.jl" begin
    include("rosenbrock.jl")
    include("ADtests.jl")
end


text:
"""
Translated to Python and modified based on:
https://www.mathworks.com/matlabcentral/fileexchange/34208-uniform-distribution-over-a-convex-polytope
Copyright (c) 2011-2012 Tim J. Benham, School of Mathematics and Physics, University of Queensland.

CHEBYCENTER Compute Chebyshev center of polytope Ax <= b.
The Chebyshev center of a polytope is the center of the largest
hypersphere enclosed by the polytope.
"""
function chebyshev_center(A, b, opt=Ipopt.Optimizer)
    m, p = size(A)

    an = sqrt.(sum(A.^2, dims=2))

    A_r = zeros((m, p+1))
    A_r[:, 1:p] = A
    A_r[:, p+1] = an

    f = zeros((p+1, 1))
    f[p+1] = -1

	model = JuMP.Model(opt)
	JuMP.set_silent(model)
	JuMP.@variable(model, 0 <= x[1:p+1])
	JuMP.@constraint(model, A_r * x .<= b)
	JuMP.@objective(model, Max, x[p+1])  # Was f'*x, but not supported by JuMP?
	#TASK: Postive or negative for objective?
    JuMP.optimize!(model)
	result = JuMP.value.(x)

	center = result[1:p]'
    radius = result[p+1]

    return center, radius
end

"""
# NOT WORKING
Finds the maximum volume inscribed ellipsoids in the given polytope Ax<b
"""
function maximum_volume_inscribed_ellipsoid(A, b, opt=GLPK.Optimizer)
    m, p = size(A)

	model = JuMP.Model(opt)
	JuMP.set_silent(model)
	JuMP.@variable(model, 0 <= d[1:p])
	JuMP.@variable(model, 0 <= B[1:p,1:p])
	for i in 1:m
		#JuMP.@constraint(model, norm(B * A[i,:]) + dot(A[i,:], d) .<= b[i])
		JuMP.@constraint(model, [b[i] - dot(A[i,:], d) ; B * A[i,:]] in SecondOrderCone())
	end
	JuMP.@objective(model, Max, logdet(B))
    JuMP.optimize!(model)

	center = value.(d)
    return center, radius
end

"""
Translated to Julia and modified based on:
https://www.mathworks.com/matlabcentral/fileexchange/34208-uniform-distribution-over-a-convex-polytope
Copyright (c) 2011-2012 Tim J. Benham, School of Mathematics and Physics, University of Queensland.

CPRND Draw from the uniform distribution over a convex polytope.
	X = cprnd(N,A,b) Returns an N-by-P matrix of random vectors drawn
	from the uniform distribution over the interior of the polytope
	defined by Ax <= b. A is a M-by-P matrix of constraint equation
	coefficients. b is a M-by-1 vector of constraint equation
	constants

...
- 'method::String': 'gibbs': Gibbs sampler, 'hitandrun': Hit and Run, 'achr': Adaptive Centering Hit-and-Run
- 'iso::Integer': Isotropic Transformation (0: no xfrm, 1: Xfrm during runup, 2: Xfrm throughout sampling
- 'runup::Integer': # of initial iterations of the algorithm in the untransformed space for gibbs or hitandrun
- 'discard::Integer': # of initial samples (post run_up) to discard. Randomly selects n_samples from the
		discard + n_samples samples.
...
"""
function sample_polytope_cprnd(A, b, x0, n_samples=1; seed=nothing, method="achr", runup=nothing, iso=nothing, discard=nothing)
    m, p = size(A)

    # Setting up default options
    if isnothing(iso)
        if method != "achr"
            iso = 2
        else
            iso = 0
		end
	end

    if isnothing(runup)
		if method == "achr"
			runup = 10 * (p+1)
		elseif iso != 0
			runup = 10 * p * (p + 1)
		else
			runup = 0
		end
	end

    if isnothing(discard)
        if method == "achr"
            discard = 1 * (p + 1)
        else
            discard = runup
		end
	end

    # Initializations
    rng = Random.MersenneTwister(seed)

    rnd_unit_vector(rng, p) = normalize(randn(rng, (p,1)))

    X = zeros((p, n_samples + runup + discard))

    n = 0  # Number of samples generated
    x = x0

    # Containers for mean, covariance, and isotropic transform tracking
    M_inc = zeros((p, 1))  # Incremental mean
    S_inc = zeros((p, p))  # Incremental sum of products

    S = I(p)
    T = I(p)
    W = A

    while n < (n_samples + runup + discard)
        y = x

        if iso > 0
            if n == runup | (iso == 2 & n > runup)
                T = cholesky(S).L
				W = A * T
			end
            y = T \ y
		end

        if method == "gibbs"
            for i in 1:p
                e = circshift(Matrix(I, p, 1), i - 1)
                z = W * e  # Selects row i?
                c = (b - W * y) ./ z
                t_min = maximum(c[z .< 0])
                t_max = minimum(c[z .> 0])

                y = y + (t_min + (t_max - t_min) * Random.rand(rng)) * e
			end
        elseif method == "hitandrun"
			u = rnd_unit_vector(rng, p)
            z = W * u
            c = (b - W * y) ./ z
            t_min = maximum(c[z .< 0])
            t_max = minimum(c[z .> 0])

            y = y + (t_min + (t_max - t_min) * Random.rand(rng)) * u
        else  # achr
            if n < runup
                u = rnd_unit_vector(rng, p)
            else
                v = X[:, rand(rng, 1:n)]
                u = normalize(v - M_inc)
			end

            z = A * u
            c = (b - W * y) ./ z
            t_min = maximum(c[z .< 0])
            t_max = minimum(c[z .> 0])

            y = y + (t_min + (t_max - t_min) .* Random.rand(rng)) .* u
		end

        if iso > 0
            x = T * y
        else
            x = y
		end
        X[:, n+1] = x
        n = n + 1

        delta_0 = x - M_inc  # Delta between new point & old mean
        M_inc = M_inc + delta_0 / n  # Sample mean
        delta_1 = x - M_inc  # Delta between new point & new mean

        if iso > 0
			if n > 1
				# Time Bottleneck here for iso
				S_inc = S_inc + (n - 1) / (n * n) * delta_0 * delta_0' + delta_1 * delta_1'

				S_0 = S
				S = S_inc / (n - 1)
			else
				S = Matrix(I,p,p)
			end
			S = (S + S') / 2  # Make array Hermitian Symmetric since computational error is likely
		end
	end
    X = X[:, runup + 1: n_samples + discard + runup]

	#TEST: Make discard pick random values from after runup to n_samples
	keep = sample(rng, 1:(n_samples + discard), n_samples, replace=false)
	X = X[:, keep]
    return X
end

"Used to call sample_uniform from create_samples"
function sample_uniform(A, b, x0, n_samples=1; nominal_load=nothing, dist_perc=0.2)
	isnothing(nominal_load) && error("nominal_load must not be nothing. Include nominal load in the
									  sampler_opts dictionary with the key ':nominal_load'")
	samples = sample_uniform(nominal_load, dist_perc=dist_perc)
	for _ in 1:n_samples
		s = sample_uniform(nominal_load, dist_perc=dist_perc)
		samples = vcat(samples, s)
	end
	return samples
end

"""
Samples from a uniform distribution of the given percentage
around the given base load, x_0 = (pd\\_0 + qd\\_0). Power factors
are held constant.

x_il = U[0.8 x_i0, 1.2 x_i0] for all i in load buses
"""
function sample_uniform(nominal_load; dist_perc=0.2)
	@assert length(nominal_load) % 2 == 0 "Nominal load profile must contain an even number of load values (active and reactive load pairs). Currently contains $(length(nominal_load))."
	num_loads = Int(length(nominal_load) / 2)
	sampled_dist = (1 .+ dist_perc .* (rand(num_loads) .* 2 .- 1))

	new_load = nominal_load .* [sampled_dist ; sampled_dist]
	return new_load
end

"Used to call sample_uniform_w_pf from create_samples"
function sample_uniform_w_pf(A, b, x0, n_samples=1; nominal_load=nothing, dist_perc=0.2, pf_range=(0.8,1.0))
	isnothing(nominal_load) && error("nominal_load must not be nothing. Include nominal load in the
									  sampler_opts dictionary with the key ':nominal_load'")

	samples = sample_uniform_w_pf(nominal_load, dist_perc=dist_perc, pf_range=pf_range)
	for _ in 1:n_samples
		s = sample_uniform_w_pf(nominal_load, dist_perc=dist_perc, pf_range=pf_range)
		samples = vcat(samples, s)
	end
	return samples
	return
end

"""
Samples from a uniform distribution of the given percentage
around the given base load real power, pd_0. Then sample a
power factor from the given pf_range to determine the reactive
load of the sample. By default the pf_range refers to lagging
power factors. If a value larger than 1.0 is given in the range,
it will be converted to a leading power factor of pf = (1.0 - x).

- p\\_il = U[(1-dist) * p_i0, (1+dist) * p\\_i0] for all i in load buses
- d\\_il = U[pf\\_range] for all i in load buses
- q\\_il = tan(d\\_il) * p\\_il
"""
function sample_uniform_w_pf(nominal_load; dist_perc=0.2, pf_range=(0.8,1.0))
	@assert length(nominal_load) % 2 == 0 "Nominal load profile must contain an even number of load values (active and reactive load pairs). Currently contains $(length(nominal_load))."
	num_loads = Int(length(nominal_load) / 2)
	sampled_pl = nominal_load[1:num_loads] .* (1 .+ dist_perc .* (rand(num_loads) .* 2 .- 1))
	sampled_pf = (pf_range[2] - pf_range[1]) .* rand(num_loads) .+ pf_range[1]
	if sampled_pf > 1.0
		sampled_pf = 1.0 - sampled_pf
	end
	sampled_ql = sampled_pl .* tan.(acos.(sampled_pf))

	new_load = [sampled_pl; sampled_ql]
	return new_load
end

"Return the unit vector of vector, v"
function normalize(v)
    norm_ = norm(v)
    if norm_ == 0
        norm_ = eps()
	end
    return v / norm_
end


text:
# Code taken from FluxML/Optimisers.jl

module Optimisers

include("apply.jl")
include("rules.jl")

end # module


text:
@testset "4.6.7 (d trig)^m (a+b (c csc)^n)^p" begin
    (a, b, c, d, x, ) = @variables a b c d x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*(a+b*csc(c+d*x)^m)^n=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(a+b*csc(c+d*x)^m)^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(a+b*csc(c+d*x)^2)^n=#

    #= ::Subsubsection::Closed:: =#
    #=n=#

    @test_int [(a + b*csc(c + d*x)^2)^4, x, 4, a^4*x - (b*(2*a + b)*(2*a^2 + 2*a*b + b^2)*cot(c + d*x))/d - (b^2*(6*a^2 + 8*a*b + 3*b^2)*cot(c + d*x)^3)/(3*d) - (b^3*(4*a + 3*b)*cot(c + d*x)^5)/(5*d) - (b^4*cot(c + d*x)^7)/(7*d)]
    @test_int [(a + b*csc(c + d*x)^2)^3, x, 4, a^3*x - (b*(3*a^2 + 3*a*b + b^2)*cot(c + d*x))/d - (b^2*(3*a + 2*b)*cot(c + d*x)^3)/(3*d) - (b^3*cot(c + d*x)^5)/(5*d)]
    @test_int [(a + b*csc(c + d*x)^2)^2, x, 4, a^2*x - (b*(2*a + b)*cot(c + d*x))/d - (b^2*cot(c + d*x)^3)/(3*d)]
    @test_int [(a + b*csc(c + d*x)^2)^1, x, 3, a*x - (b*cot(c + d*x))/d]
    @test_int [1/(a + b*csc(c + d*x)^2)^1, x, 3, x/a - (sqrt(b)*atan((sqrt(a + b)*tan(c + d*x))/sqrt(b)))/(a*sqrt(a + b)*d)]
    @test_int [1/(a + b*csc(c + d*x)^2)^2, x, 5, x/a^2 + (sqrt(b)*(3*a + 2*b)*atan((sqrt(b)*cot(c + d*x))/sqrt(a + b)))/(2*a^2*(a + b)^(3/2)*d) + (b*cot(c + d*x))/(2*a*(a + b)*d*(a + b + b*cot(c + d*x)^2))]
    @test_int [1/(a + b*csc(c + d*x)^2)^3, x, 6, x/a^3 + (sqrt(b)*(15*a^2 + 20*a*b + 8*b^2)*atan((sqrt(b)*cot(c + d*x))/sqrt(a + b)))/(8*a^3*(a + b)^(5/2)*d) + (b*cot(c + d*x))/(4*a*(a + b)*d*(a + b + b*cot(c + d*x)^2)^2) + (b*(7*a + 4*b)*cot(c + d*x))/(8*a^2*(a + b)^2*d*(a + b + b*cot(c + d*x)^2))]
    @test_int [1/(a + b*csc(c + d*x)^2)^4, x, 7, x/a^4 + (sqrt(b)*(35*a^3 + 70*a^2*b + 56*a*b^2 + 16*b^3)*atan((sqrt(b)*cot(c + d*x))/sqrt(a + b)))/(16*a^4*(a + b)^(7/2)*d) + (b*cot(c + d*x))/(6*a*(a + b)*d*(a + b + b*cot(c + d*x)^2)^3) + (b*(11*a + 6*b)*cot(c + d*x))/(24*a^2*(a + b)^2*d*(a + b + b*cot(c + d*x)^2)^2) + (b*(19*a^2 + 22*a*b + 8*b^2)*cot(c + d*x))/(16*a^3*(a + b)^3*d*(a + b + b*cot(c + d*x)^2))]

    #= ::Subsubsection::Closed:: =#
    #=n/2=#

    @test_int [(a + b*csc(c + d*x)^2)^(5/2), x, 8, -((a^(5/2)*atan((sqrt(a)*cot(c + d*x))/sqrt(a + b + b*cot(c + d*x)^2)))/d) - (sqrt(b)*(15*a^2 + 10*a*b + 3*b^2)*atanh((sqrt(b)*cot(c + d*x))/sqrt(a + b + b*cot(c + d*x)^2)))/(8*d) - (b*(7*a + 3*b)*cot(c + d*x)*sqrt(a + b + b*cot(c + d*x)^2))/(8*d) - (b*cot(c + d*x)*(a + b + b*cot(c + d*x)^2)^(3/2))/(4*d)]
    @test_int [(a + b*csc(c + d*x)^2)^(3/2), x, 7, -((a^(3/2)*atan((sqrt(a)*cot(c + d*x))/sqrt(a + b + b*cot(c + d*x)^2)))/d) - (sqrt(b)*(3*a + b)*atanh((sqrt(b)*cot(c + d*x))/sqrt(a + b + b*cot(c + d*x)^2)))/(2*d) - (b*cot(c + d*x)*sqrt(a + b + b*cot(c + d*x)^2))/(2*d)]
    @test_int [(a + b*csc(c + d*x)^2)^(1/2), x, 6, -((sqrt(a)*atan((sqrt(a)*cot(c + d*x))/sqrt(a + b + b*cot(c + d*x)^2)))/d) - (sqrt(b)*atanh((sqrt(b)*cot(c + d*x))/sqrt(a + b + b*cot(c + d*x)^2)))/d]
    @test_int [1/(a + b*csc(c + d*x)^2)^(1/2), x, 3, -(atan((sqrt(a)*cot(c + d*x))/sqrt(a + b*csc(c + d*x)^2))/(sqrt(a)*d)), -(atan((sqrt(a)*cot(c + d*x))/sqrt(a + b + b*cot(c + d*x)^2))/(sqrt(a)*d))]
    @test_int [1/(a + b*csc(c + d*x)^2)^(3/2), x, 4, -(atan((sqrt(a)*cot(c + d*x))/sqrt(a + b + b*cot(c + d*x)^2))/(a^(3/2)*d)) + (b*cot(c + d*x))/(a*(a + b)*d*sqrt(a + b + b*cot(c + d*x)^2))]
    @test_int [1/(a + b*csc(c + d*x)^2)^(5/2), x, 6, -(atan((sqrt(a)*cot(c + d*x))/sqrt(a + b + b*cot(c + d*x)^2))/(a^(5/2)*d)) + (b*cot(c + d*x))/(3*a*(a + b)*d*(a + b + b*cot(c + d*x)^2)^(3/2)) + (b*(5*a + 3*b)*cot(c + d*x))/(3*a^2*(a + b)^2*d*sqrt(a + b + b*cot(c + d*x)^2))]
    @test_int [1/(a + b*csc(c + d*x)^2)^(7/2), x, 7, -(atan((sqrt(a)*cot(c + d*x))/sqrt(a + b + b*cot(c + d*x)^2))/(a^(7/2)*d)) + (b*cot(c + d*x))/(5*a*(a + b)*d*(a + b + b*cot(c + d*x)^2)^(5/2)) + (b*(9*a + 5*b)*cot(c + d*x))/(15*a^2*(a + b)^2*d*(a + b + b*cot(c + d*x)^2)^(3/2)) + (b*(33*a^2 + 40*a*b + 15*b^2)*cot(c + d*x))/(15*a^3*(a + b)^3*d*sqrt(a + b + b*cot(c + d*x)^2))]

    @test_int [(1 + csc(x)^2)^(3/2), x, 6, -2*asinh(cot(x)/sqrt(2)) - atan(cot(x)/sqrt(2 + cot(x)^2)) - (1/2)*cot(x)*sqrt(2 + cot(x)^2)]
    @test_int [sqrt(1 + csc(x)^2), x, 5, -asinh(cot(x)/sqrt(2)) - atan(cot(x)/sqrt(2 + cot(x)^2))]
    @test_int [1/sqrt(1 + csc(x)^2), x, 3, -atan(cot(x)/sqrt(2 + cot(x)^2))]

    @test_int [(1 - csc(x)^2)^(3/2), x, 4, (1/2)*cot(x)*sqrt(-cot(x)^2) + sqrt(-cot(x)^2)*log(sin(x))*tan(x)]
    @test_int [sqrt(1 - csc(x)^2), x, 3, sqrt(-cot(x)^2)*log(sin(x))*tan(x)]
    @test_int [1/sqrt(1 - csc(x)^2), x, 3, -((cot(x)*log(cos(x)))/sqrt(-cot(x)^2))]

    @test_int [(-1 + csc(x)^2)^(3/2), x, 4, (-(1/2))*(cot(x)^2)^(3/2)*tan(x) - sqrt(cot(x)^2)*log(sin(x))*tan(x)]
    @test_int [sqrt(-1 + csc(x)^2), x, 3, sqrt(cot(x)^2)*log(sin(x))*tan(x)]
    @test_int [1/sqrt(-1 + csc(x)^2), x, 3, -((cot(x)*log(cos(x)))/sqrt(cot(x)^2))]

    @test_int [(-1 - csc(x)^2)^(3/2), x, 7, -2*atan(cot(x)/sqrt(-2 - cot(x)^2)) - atanh(cot(x)/sqrt(-2 - cot(x)^2)) + (1/2)*cot(x)*sqrt(-2 - cot(x)^2)]
    @test_int [sqrt(-1 - csc(x)^2), x, 6, atan(cot(x)/sqrt(-2 - cot(x)^2)) + atanh(cot(x)/sqrt(-2 - cot(x)^2))]
    @test_int [1/sqrt(-1 - csc(x)^2), x, 3, -atanh(cot(x)/sqrt(-2 - cot(x)^2))]

    #= ::Subsection:: =#
    #=Integrands*of*the*form*(a+b*csc(c+d*x)^3)^n=#

    #= ::Subsection:: =#
    #=Integrands*of*the*form*(a+b*csc(c+d*x)^4)^n=#
    end

max_stars_repo_path: test/4 Trig functions/4.6 Cosecant/4.6.7 (d trig)^m (a+b (c csc)^n)^p.jl


text:
using Plots
using LinearAlgebra
include("mess.jl")
using SparseArrays
using Arpack
#komvos1 komvos2 komvos3
xy = [ 0 0;
       0 1]; #orismos komvon san xyz
nodes = size(xy)[2]; #arithmos komvon

#komvos1_start komvos1_end
numel = [1 2

]; #orismos ton stixion os grammes me toys akraious komvous
N = size(numel)[1]; #arithmos stoixion
dofs = 3*nodes;
#orismos elastikotitas stixion
a = 20e-2;
E = 220e9.*ones(N,1); I = a^4 .*ones(N,1)./12; A = a^2 .*ones(N,1) #is square with axa
rho = 7.8e3;

#orismos ton desmeumenon komvon
BC = zeros(dofs,1)
#ignored = 1:3:dofs;#ignore turn of axis=ignore freedom
BC[1:3] .= 1; #BC[ignored].=1;
fixedofs = findall(x->x==1,BC[:])
freedofs = findall(x->x==0,BC[:])
#boundary conditions
U = zeros(dofs+3*N,1);
F = zeros(dofs+3*N,1); F[[4,5]] .= 1e3;
#euresi olon ton deikton pou tha prospelasoun ton K
alldofs = zeros(N,6);
#orismos olikou mitroou stivarotitas
global K = sparse([],[],Float64[],dofs+3*N,dofs+3*N);
global M = sparse([],[],Float64[],dofs+3*N,dofs+3*N);
ignored = []
#global K = zeros(dofs,dofs);
#start main loop to calculate stiffness matrix
for i=1:N
    i1 = numel[i,1]
    i2 = numel[i,2]
    #calclulate directions
    x1 = xy[1,i1]; x2 = xy[1,i2]; dx = x2-x1;
    y1 = xy[2,i1]; y2 = xy[2,i2]; dy = y2-y1;
    L = sqrt(dx^2+dy^2);
    nx = dx/L; ny = dy/L;
    #direction matrix
    c = nx; s = ny;
    Bbend = [-s c 0  0 0 0  0 0 0;
              0 0 1  0 0 0  0 0 0;
              0 0 0 -s c 0  0 0 0;
              0 0 0  0 0 1  0 0 0;
              0 0 0  0 0 0 -s c 0;
              0 0 0  0 0 0  0 0 1];
    Btens = [c s 0 0 0 0;
             0 0 c s 0 0;
             0 0 0 0 c s];
    kbend = E[i]*I[i]/L^3;

    #insert tensile
    ktens = E[i]*A[i]/L
    global Ketens = (ktens/3)*[ 7 -8  1;
                               -8 16 -8;
                                1 -8  7];
    Ketens = Btens'*Ketens*Btens
    global Metens = (rho*A[i]*L/15)*[2  1 -1;
                                     1  8  1;
                                    -1  1  2];
    Metens = Btens'*Metens*Btens;

    Kebend = (kbend/35/(12+13*L)^2).*[4*(25776+55560*L+30695*L^2) 2*L*(80592+173240*L+95025*L^2)   -512*(927+2040*L+1085*L^2) 512*L*(537+1175*L+665*L^2) -4*(56016+127800*L+75985*L^2) 2*L*(18384+42680*L+25585*L^2);
                     0  L^2*4*(118080+25400*L+13835*L^2) -256*L*(492+1040*L+555*L^2) 128*L^2*(354+770*L+435*L^2) -2*L*(17616+40120*L+23985*L^2) 2*L^2*(2832+6600*L+3995*L^2);
                     0  0     8192*(126+270*L+145*L^2)  1024*L*(12+40*L+25*L^2) -512*(1044+2280*L+1235*L^2) 256*L*(516+1120+605*L^2);
                     0  0     0   1024*L^2*(183+405*L+230*L^2) -512*L*(561+1255*L+715*L^2) 128*L^2*(378+850*L+485*L^2);
                     0  0     0     0   4*(189648+419640*L+234065*L^2) -2*L*(84432+186040*L+103025*L^2);
                     0  0     0     0  0   L^2*4*(12192+26680*L+14635*L^2)];
    Kebend = Symmetric(Kebend);
    Kebend = Bbend'*Kebend*Bbend;
    Mebend = rho*A[i]*(L/(12+13*L)^2)*[(337248+750276*L+417883*L^2)/15015   L*(57552+132764*L+76767*L^2)/45045  8*(18408+41602*L+23423*L^2)/30030     -2*L*(31848+60578*L+27535*L^2)/45045 (81168+181304*L+10061*L^2)/30030 -L*(48720+104128*L+54703*L^2)/45045;
                                       0   L^2*(4236+10167*L+6131*L^2)/45045 8*L*(5772+13573*L+7907*L^2)/45045 -2*L^2*(764+1381*L+556*L^2)/15015 L*(58704+139712*L+81263*L^2)/180180 -L^2*(1840+4208*L+2343*L^2)/60060;
                                       0     0   1024*(858+1846+993*L^2)/15015    128*L*(156+556*L+415*L^2)/45045 8*(15912+32706*L+16783*L^2)/15015 -8*L*(4524+9125*L+4587*L^2)/45045;
                                       0     0    0     128*L^2*(471+1029*L+568*L^2)/45045 2*L*(42456+98782*L+56745*L^2)/45045 -2*L^2*(1076+2515*L+1441*L^2)/15015;
                                       0     0    0    0     (316032+673868*L+359463*L^2)/15015 -L*(49440+103456*L+54197*L^2)/45045;
                                       0     0    0    0    0    L^2*(300+615*L+316*L^2)/4095];
    Mebend = Symmetric(Mebend);
    Mebend = Bbend'*Mebend*Bbend;

    tensiledofs = [1 2 4 5 7 8];
    Ke = Kebend .+0;
    Ke[tensiledofs,tensiledofs] = Ke[tensiledofs,tensiledofs][1,:,1,:]  + Ketens;
        #println(size(Ke),"\n")
    ending = dofs+3*N
    edofs = [3*i1-2 3*i1-1 ending-3*i+3 ending-3*i+2 ending-3*i+1 3*i1+3  3*i2-2 3*i2-1 3*i2]
    push!(ignored,ending-3*i+3, ending-3*i+2, ending-3*i+1)
    K[edofs,edofs] = K[edofs,edofs][1,:,1,:] + Ke;
    #println(K,"\n")
    Me = Mebend .+0;
    Me[tensiledofs,tensiledofs] = Me[tensiledofs,tensiledofs][1,:,1,:]  + Metens;
    M[edofs,edofs] = M[edofs,edofs][1,:,1,:] + Me;
end
Kfinal = K[vcat(freedofs,ignored),vcat(freedofs,ignored)]
Mfinal = M[vcat(freedofs,ignored),vcat(freedofs,ignored)]
rhs = F
rhs = rhs[vcat(freedofs,ignored)]
rhs = rhs - K[vcat(freedofs,ignored),fixedofs]*U[fixedofs]

Ufree = Kfinal\rhs
U[freedofs] = Ufree[freedofs]
#Ufr = reshape(U,(3,nodes))
show(U[:])
#println(varinfo())
f = inv(Matrix(Mfinal))*Matrix(Kfinal)
omega = sqrt.(abs.(eigvals(f)))/(2*pi)
s = eigvecs(f)
plot(omega)

#calculate

max_stars_repo_path: axnai julia/test.jl

max_stars_repo_name: billyon/first

max_issues_repo_path: axnai julia/test.jl

max_issues_repo_name: billyon/first

max_forks_repo_path: axnai julia/test.jl

max_forks_repo_name: billyon/first


text:
#---------------------------------------------------------------------#
#This function computes the Derivative Mapping for General 2D Quad Grids.
#Written by F.X. Giraldo on 5/2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function map_deriv(ψ,dψ,f,Np,Nq,DFloat)

    #Allocate arrays
    f_ξ=zeros(DFloat,Nq,Nq,Nq)
    f_η=zeros(DFloat,Nq,Nq,Nq)
    f_ζ=zeros(DFloat,Nq,Nq,Nq)

    #Loop through Integration Points
    for n=1:Nq, m=1:Nq, l=1:Nq
        sum_ξ=0
        sum_η=0
        sum_ζ=0
        for k=1:Np, j=1:Np, i=1:Np
            sum_ξ=sum_ξ + dψ[i,l]*ψ[j,m]*ψ[k,n]*f[i,j,k]
            sum_η=sum_η + ψ[i,l]*dψ[j,m]*ψ[k,n]*f[i,j,k]
            sum_ζ=sum_ζ + ψ[i,l]*ψ[j,m]*dψ[k,n]*f[i,j,k]
        end #i,j,k
        f_ξ[l,m,n]=sum_ξ
        f_η[l,m,n]=sum_η
        f_ζ[l,m,n]=sum_ζ
    end #l,m,n
    return (f_ξ,f_η,f_ζ)
end


text:
"""
    BeliefMCTSSolver(mcts_solver, updater)

The belief mcts solver solves POMDPs by modeling them as an MDP on the belief space. The `updater` is used to update the belief as part of the belief MDP generative model.

Example:

    using ParticleFilters
    using POMDPModels
    using MCTS

    pomdp = BabyPOMDP()
    updater = SIRParticleFilter(pomdp, 1000)

    solver = BeliefMCTSSolver(DPWSolver(), updater)
    planner = solve(solver, pomdp)

    simulate(HistoryRecorder(max_steps=10), pomdp, planner, updater)
"""
mutable struct BeliefMCTSSolver
    solver::AbstractMCTSSolver
    updater::Updater
end

function POMDPs.solve(sol::BeliefMCTSSolver, p::POMDP)
    bmdp = GenerativeBeliefMDP(p, sol.updater)
    return solve(sol.solver, bmdp)
end


text:
@testset "Transformer" begin
    t = Transformer(10, 3, 15, 20)
    x = randn(10, 7, 3)

    @test size(t(x)) == (10, 7, 3)
    @test size(t(x[:, :, 2])) == (10, 7)
end


text:
using Pkg
cd(@__DIR__)
cd("..")
Pkg.activate(".")

using Base: runtests

using Behavior
using Behavior.Gherkin
using TypeDBClient

g = TypeDBClient
client = g.CoreClient("localhost",1729)

rootpath = joinpath(@__DIR__,"behaviour")
featurepath = joinpath(@__DIR__,"behaviour/features")
stepspath = joinpath(@__DIR__,"behaviour")
configpath = joinpath(@__DIR__,"behaviour/config/ConfigEnvironment.jl")

p = ParseOptions(allow_any_step_order = true)

function run_tests(tag::String = "")
    runspec(rootpath; featurepath = featurepath, stepspath = stepspath,  parseoptions=p, execenvpath = configpath, tags=tag)
end

results = run_tests("not @ignore-typedb-core")
# results = run_tests("@failure")
# results = run_tests("@actual, not @ignore-typedb-core")

if !all(values(results))
    throw("TestSuite failed. Please review the results")
else
    @info "All tests passed. Well done! 😄 "
end


text:
module Trade

using JSON3, HTTP, Dates

export getTrades, getOpenTrades, getTrade, closeTrade, clientExtensions, setTradeOrders, cancelTradeOrders

"Detailed Trade struct from Oanda"
mutable struct trade
    averageClosePrice # The average closing price of the Trad
    clientExtensions # The client extensions of the Trade
    closingTransactionIDs # The IDs of the Transactions that have closed portions of this Trade
    closeTime # The date/time when the Trade was fully closed
    currentUnits # Current units of the trade (- is short + is long)
    dividend # The dividend paid for this Trade
    financing # The financing paid / collected for this trade
    id # The id of the trade
    initialUnits # Initial opening units of the trade (- is short + is long)
    initialMarginRequired # The margin required at the time the Trade was created
    instrument # Instrument of the trade
    marginUsed # Margin currently used by the Trade
    openTime # The time the trade was opened
    price # The price the trade is set at
    realizedPL # The profit / loss of the trade that has been incurred
    state # current state of the trade
    takeProfitOrder # Full representation of the Trade’s Take Profit Order
    stopLossOrder # Full representation of the Trade’s Stop Loss Order
    trailingStopLossOrder # Full representation of the Trade’s Trailing Stop Loss Order
    unrealizedPL # The profit / loss of the trade that hasnt been incurred
    #= Better alternative that requires complete implementation of Order.jl
    takeProfitOrder::takeProfitOrder
    stopLossOrder::stopLossOrder
    trailingStopLossOrder::trailingStopLossOrder
    =#
    trade() = new()
end

"Coerce a given 'trade' into its proper types (Used internally)"
function coerceTrade(trade::trade)
    RFC = Dates.DateFormat("yyyy-mm-ddTHH:MM:SS.sssssssssZ")

    isdefined(trade,:averageClosePrice) && (trade.averageClosePrice = parse(Float32,trade.averageClosePrice))
    isdefined(trade,:closeTime) && (trade.closeTime = DateTime(first(trade.closeTime,23),RFC))
    trade.currentUnits = parse(Float32, trade.currentUnits)
    trade.initialUnits = parse(Float32, trade.initialUnits)
    trade.initialMarginRequired = parse(Float32,trade.initialMarginRequired)
    trade.financing = parse(Float32, trade.financing)
    # ID is left as a string, makes more sense to me for usage
    isdefined(trade,:marginUsed) && (trade.marginUsed = parse(Float32, trade.marginUsed))
    trade.openTime = DateTime(first(trade.openTime,23),RFC)
    trade.price = parse(Float32, trade.price)
    trade.realizedPL = parse(Float32, trade.realizedPL)
    isdefined(trade,:unrealizedPL) && (trade.unrealizedPL = parse(Float32, trade.unrealizedPL))
    #= Requires complete implementation of Order.jl
    isdefined(trade,:takeProfitOrder) && (trade.takeProfitOrder = coerceTakeProfitOrder(trade.takeProfitOrder))
    isdefined(trade,:stopLossOrder) && (trade.stopLossOrder = coerceStopLossOrder(trade.stopLossOrder))
    isdefined(trade,:trailingStopLossOrder) && (trade.trailingStopLossOrder = coerceTrailingStopLossOrder(trade.trailingStopLossOrder))
    =#
    return trade
end

mutable struct trades
    trades::Vector{trade}
    lastTransactionID

    trades() = new()
end

# Declaring JSON3 struct types
JSON3.StructType(::Type{trades}) = JSON3.Mutable()
JSON3.StructType(::Type{trade}) = JSON3.Mutable()

# ------------------------------------------------------------------------------------
# /accounts/{accountID}/trades Endpoint
# /accounts/{accountID}/openTrades Endpoint
# ------------------------------------------------------------------------------------
"""
    getTrades(config::config, instrument::String, state::String="OPEN", count::Int=50; kwargs...)

Return an array of trade struct

# Arguments
- 'config::config': a valid struct with user configuration data
- 'state::String": the state to filter the requested Trades by (OPEN, CLOSED, CLOSE_WHEN_TRADEABLE, ALL)
- 'count::Int': max number of trades to return

# Keyword Arguments
- 'instrument::String": a valid instrument (e.g. "EUR_USD")
- 'ids::String' List of trades to retrieve as ID values separated by commas
- 'beforeID::String' The maximum trade ID to return

# Examples

    getTrades(userConfig,"CLOSED";instrument="EUR_USD")

"""
function getTrades(config, state::String="ALL", count::Int=50; kwargs...)

    r = HTTP.get(string("https://", config.hostname, "/v3/accounts/", config.account, "/trades"),
        ["Authorization" => string("Bearer ", config.token)];
        query = push!(Dict(), "state" => state, "count" => count, kwargs...))

    temp = JSON3.read(r.body,trades)
    # type coersions
    for t in temp.trades
        t = coerceTrade(t)
    end
    temp.trades # Not returning lastTransactionID. That info belongs to Transaction.jl
end

"""
    getOpenTrades(config::config)

Return an array of trade struct

# Arguments
- 'config::config': a valid struct with user configuracion data

# Examples

    getOpenTrades(userconfig)

"""
function getOpenTrades(config)

    r = HTTP.get(string("https://", config.hostname, "/v3/accounts/", config.account, "/openTrades"),
        ["Authorization" => string("Bearer ", config.token)])

    temp = JSON3.read(r.body,trades)

    # type coersions
    for t in temp.trades
        t = coerceTrade(t)
    end
    return temp.trades # Not returning lastTransactionID. That info belongs to Transaction.jl
end

# ------------------------------------------------------------------------------------
# /accounts/{accountID}/trades/{tradeSpecifier} Endpoint
# ------------------------------------------------------------------------------------
mutable struct singleTrade
    trade::trade
    lastTransactionID

    singleTrade() = new()
end

# Declaring JSON3 struct types
JSON3.StructType(::Type{singleTrade}) = JSON3.Mutable()

"""
   getTrade(config::config, tradeID::String)

Return a specific trade

# Arguments
- 'config::config': a valid struct with user configuracion data
- 'id::string': a valid trade ID

# Examples

    getTrades(userconfig,"66")

"""
function getTrade(config, tradeID::String)

    r = HTTP.get(string("https://", config.hostname, "/v3/accounts/", config.account,"/trades/",tradeID),
        ["Authorization" => string("Bearer ", config.token)])

    temp = JSON3.read(r.body,singleTrade)

    return coerceTrade(temp.trade) # Not returning lastTransactionID. That info belongs to Transaction.jl
end

# ------------------------------------------------------------------------------------
# //accounts/{accountID}/trades/{tradeSpecifier}/close Endpoint
# ------------------------------------------------------------------------------------

# close trade endpoint response struct
mutable struct closeUnitsResp
    orderCreateTransaction::Dict{String,Any}
    orderFillTransaction::Dict{String,Any}
    orderCancelTransaction::Dict{String,Any}
    relatedTransactionIDs::Vector{String}
    lastTransactionID

    closeUnitsResp() = new()
end

# close trade endpoint request struct
struct closeUnits
    units::String

    closeUnits(x::String) = new(x)

    function closeUnits(x::Real)
        str = string(x)
        new(str)
    end

end

# Declaring JSON3 struct types
JSON3.StructType(::Type{closeUnits}) = JSON3.Struct()
JSON3.StructType(::Type{closeUnitsResp}) = JSON3.Mutable()

"""
    closeTrade(config::config, tradeID::String, units::Union{Real,String}="ALL")

Return an array of trade struct

# Arguments
- 'config::config': a valid struct with user configuracion data
- 'tradeID::string': a valid trade ID
- 'units::Union{Number,String}': how much of the Trade to close in units or "ALL"

# Examples

    closeTrade(userconfig,"66","ALL")

"""
function closeTrade(config, tradeID::String, units::Union{Real,String}="ALL")

    r = HTTP.put(string("https://", config.hostname, "/v3/accounts/", config.account,"/trades/",tradeID,"/close"),
        ["Authorization" => string("Bearer ", config.token),"Content-Type" => "application/json"], JSON3.write(closeUnits(units)))

    return JSON3.read(r.body,closeUnitsResp)

end

# ------------------------------------------------------------------------------------
# //accounts/{accountID}/trades/{tradeSpecifier}/clientExtensions Endpoint
# ------------------------------------------------------------------------------------

# clientExtensions response struct
mutable struct clientExtensionsResp
    tradeClientExtensionsModifyTransaction::Dict{String,Any}
    relatedTransactionIDs::Vector{String}
    lastTransactionID

    clientExtensionsResp() = new()
end

# clientExtension request structs
struct extensions
    id::String
    tag::String
    comment::String
end

# For JSON parsing
struct clientExtensions
    clientExtensions::extensions
end

# Declaring JSON3 struct types
JSON3.StructType(::Type{clientExtensions}) = JSON3.Struct()
JSON3.StructType(::Type{extensions}) = JSON3.Struct()
JSON3.StructType(::Type{clientExtensionsResp}) = JSON3.Mutable()

"""
    clientExtensions(config::config, tradeID::String; clientID::String="", tag::String="", comment::String="")

Lets add user information to a specific Trade

# Arguments
- 'config::config': a valid struct with user configuracion data
- 'tradeID::string': a valid trade ID
- clientID, tag and comment: strings with the user information

# Example

    clientExtensions(userconfig,"66", clientID="007",tag="foo")

"""
function clientExtensions(config, tradeID::String; clientID::String="", tag::String="", comment::String="")

    data = clientExtensions(extensions(clientID, tag, comment))

    r = HTTP.put(string("https://", config.hostname, "/v3/accounts/", config.account,"/trades/",tradeID,"/clientExtensions"),
        ["Authorization" => string("Bearer ", config.token),"Content-Type" => "application/json"], JSON3.write(data))

    return JSON3.read(r.body, clientExtensionsResp)

end

# ------------------------------------------------------------------------------------
# //accounts/{accountID}/trades/{tradeSpecifier}/orders Endpoint
# ------------------------------------------------------------------------------------

# orders endpoint response struct
mutable struct tradeOrdersResponse
    takeProfitOrderCancelTransaction::Dict{String,Any}
    takeProfitOrderTransaction::Dict{String,Any}
    takeProfitOrderFillTransaction::Dict{String,Any}
    takeProfitOrderCreatedCancelTransaction::Dict{String,Any}
    stopLossOrderCancelTransaction::Dict{String,Any}
    stopLossOrderTransaction::Dict{String,Any}
    stopLossOrderCreatedCancelTransaction::Dict{String,Any}
    trailingStopLossOrderCancelTransaction::Dict{String,Any}
    trailingStopLossOrderTransaction::Dict{String,Any}
    relatedTransactionIDs::Vector{String}
    lastTransactionID

    tradeOrdersResponse() = new()
end

# orders endpoint request structs
mutable struct takeProfit
    price::Real
    timeInForce::String
    gtdTime::String
    # clientExtensions::extensions -> TODO
    takeProfit() = new()
end

mutable struct stopLoss
    price::Real
    distance::Real
    timeInForce::String
    gtdTime::String
    # clientExtensions::extensions -> TODO

    stopLoss() = new()
end

mutable struct trailingStopLoss
    distance::Real
    timeInForce::String
    gtdTime::String
    # clientExtensions::extensions -> TODO

    trailingStopLoss() = new()
end

mutable struct tradeOrders
    takeProfit::takeProfit
    stopLoss::stopLoss
    trailingStopLoss::trailingStopLoss

    tradeOrders()=new()
end

# Declaring JSON3 struct types and setting fields to ignore in JSON3.write if # undef
JSON3.StructType(::Type{tradeOrders}) = JSON3.Mutable()
JSON3.omitempties(::Type{tradeOrders})=(:takeProfit,:stopLoss,:trailingStopLoss)

JSON3.StructType(::Type{takeProfit}) = JSON3.Mutable()
JSON3.omitempties(::Type{takeProfit})=(:price,:timeInForce,:gtdTime)

JSON3.StructType(::Type{stopLoss}) = JSON3.Mutable()
JSON3.omitempties(::Type{stopLoss})=(:price,:distance,:timeInForce,:gtdTime)

JSON3.StructType(::Type{trailingStopLoss}) = JSON3.Mutable()
JSON3.omitempties(::Type{trailingStopLoss})=(:price,:timeInForce,:gtdTime)

JSON3.StructType(::Type{tradeOrdersResponse}) = JSON3.Mutable()

"""
    function setTradeOrders(config::config, tradeID::String; [TP::NamedTuple, SL::NamedTuple, tSL::NamedTuple ])

Create or modify the linked orders for a specific trade

# Arguments
- 'config::config': a valid struct with user configuracion data
- 'tradeID::string': a valid trade ID

# Additional Arguments
- 'TP::NamedTuple' Take Profit parameters
- 'SL::NamedTuple': Stop Loss parameters
- 'tSL::NamedTuple': Trailing Stop Loss parameters

At least one type of order parameters must be provided

# Valid order parameters
- 'price = :Real' :price to create o modify for the specific order. Valir for Stop Loss and Take Profit
- 'distance = :Real' :price distance to create o modify for the specific order. Valir for Stop Loss and Trailing Stip Loss
- 'TIF= String': time in force for the order. Valid options are: GTC, GTD, GFD, FOK, IOC. Defaults to GTC
- 'gtdTime = DateTime': time for GTD (Good unTill Date)

Price and distance are incompatible. Only one can be set for a given order.

# Example

    setTradeOrders(userconfig, "34"; TP=(price=109.5,), SL=(distance=10,TIF="FOK")) # Do not forget the comma for 1 element NamedTuples

"""
function setTradeOrders(config, tradeID::String; TP::NamedTuple=NamedTuple(),SL::NamedTuple=NamedTuple(),tSL::NamedTuple=NamedTuple())

    data = tradeOrders()

    if !isempty(TP)
        TPdetails = takeProfit()
        haskey(TP, :price) && (TPdetails.price = TP.price)
        haskey(TP, :timeInForce) && (TPdetails.timeInForce = TP.timeInForce)
        haskey(TP, :gtdTime) && (TPdetails.price = Dates.format(TP.gtdTime,"yyyy-mm-ddTHH:MM:SS.sss000000Z"))

        data.takeProfit = TPdetails
    end

    if !isempty(SL)
        SLdetails = stopLoss()
        haskey(SL, :price) && (SLdetails.price = SL.price)
        haskey(SL, :distance) && (SLdetails.distance = SL.distance)
        haskey(SL, :timeInForce) && (SLdetails.timeInForce = SL.timeInForce)
        haskey(SL, :gtdTime) && (SLdetails.price = Dates.format(SL.gtdTime,"yyyy-mm-ddTHH:MM:SS.sss000000Z"))

        data.stopLoss = SLdetails
    end

    if !isempty(tSL)
        tSLdetails = trailingStopLoss()
        haskey(tSL, :distance) && (tSLdetails.distance = tSL.distance)
        haskey(tSL, :timeInForce) && (tSLdetails.timeInForce = tSL.timeInForce)
        haskey(tSL, :gtdTime) && (tSLdetails.price = Dates.format(tSL.gtdTime,"yyyy-mm-ddTHH:MM:SS.sss000000Z"))

        data.trailingStopLoss = tSLdetails
    end

    r = HTTP.put(string("https://", config.hostname, "/v3/accounts/", config.account,"/trades/",tradeID,"/orders"),
        ["Authorization" => string("Bearer ", config.token),"Content-Type" => "application/json"], JSON3.write(data))

    return JSON3.read(r.body,tradeOrdersResponse)

end

# Definition of cancelTradeOrders structs
mutable struct nullTradeOrders
    takeProfit
    stopLoss
    trailingStopLoss

    nullTradeOrders()=new()
end

# Declaring JSON3 struct types and setting fields to ignore in JSON3.write if # undef
JSON3.StructType(::Type{nullTradeOrders}) = JSON3.Mutable()
JSON3.omitempties(::Type{nullTradeOrders})=(:takeProfit,:stopLoss,:trailingStopLoss)

"""
    function cancelTradeOrders(config::config, tradeID::String, orders2cancel::Vector{String})

Cancel linked orders of a specific trade

# Arguments
- 'config::config': a valid struct with user configuracion data
- 'tradeID::string': a valid trade ID
- 'orders2cancel::Vector{String}': list of orders to cancel.

order2cancel valid fields are "TP" for Take Profit, "SL" for Stop Loss and "tSL" for Trailing Stop Loss

# Example

    cancelTradeOrders(userconfig, "34", ["SL", "TP"])

end
"""
function cancelTradeOrders(config, tradeID::String, orders2cancel::Vector{String})

    data=nullTradeOrders()

    in("TP",orders2cancel) && (data.takeProfit=missing)
    in("SL",orders2cancel) && (data.stopLoss=missing)
    in("tSL",orders2cancel) && (data.trailingStopLoss=missing)

    r = HTTP.put(string("https://", config.hostname, "/v3/accounts/", config.account,"/trades/",tradeID,"/orders"),
    ["Authorization" => string("Bearer ", config.token),"Content-Type" => "application/json"], JSON3.write(data))

    return JSON3.read(r.body,tradeOrdersResponse)
end

end # Module

max_stars_repo_name: CalebDepatie/Julianda

max_forks_repo_name: CalebDepatie/Julianda


text:
# generate an intensity trace from a blinking fluorophore

using Revise
using SMLMSim
using Plots

## create a blinking fluorophore
γ=10000.0
q=[0 10
   1e-1 0]
f=SMLMSim.GenericFluor(γ,q)

## Simulate intensity trace

# setup sim
nframes = 1000
framerate = 10.0

## generate CTMC
starttime=0.0
endtime = (nframes) / framerate
state1=1
ctmc=SMLMSim.CTMC(q,endtime,state1)
Plots.plot(ctmc.transitiontimes,ctmc.states)

##  generate integrated photons
photons=SMLMSim.intensitytrace(f,nframes,framerate)
Plots.plot((1:nframes),photons)


text:
using DashBootstrapComponents, DashHtmlComponents

inputs = html_div([
    dbc_input(placeholder="A large input...", bs_size="lg", className="mb-3"),
    dbc_input(placeholder="A medium input...", bs_size="md", className="mb-3"),
    dbc_input(placeholder="A small input...", bs_size="sm"),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
# [https://github.com/sjkelly/MeshSlicer.jl](https://github.com/sjkelly/MeshSlicer.jl)
# ![](../img/sliced_cylinder.png)

module MeshSlicer

using ImmutableArrays
import Base.push!
import Base.show

export Bounds3, Bounds2, Face, PolygonMesh, LineSegment,
       update!, rotate!, rotate, MeshSlice, Polygon

type Bounds3{T<:Number}
    xmax::T
    ymax::T
    zmax::T
    xmin::T
    ymin::T
    zmin::T
end

type Bounds2{T<:Number}
    xmax::T
    ymax::T
    xmin::T
    ymin::T
end

type Face
    vertices::Array{Vector3{Float64}}
    normal::Vector3{Float64}
    zmin::Float64
    zmax::Float64
end

function Face(vertices::Array, normal::Vector3{Float64})
    zmax = max(max(vertices[1].e3, vertices[2].e3), vertices[3].e3)
    zmin = min(min(vertices[1].e3, vertices[2].e3), vertices[3].e3)
    return Face(vertices, normal, zmin, zmax)
end

type PolygonMesh
    bounds::Bounds3
    faces::Array{Face}
    patios::Array{Float64}
end

type LineSegment
    start::Vector2{Float64}
    finish::Vector2{Float64}
    normal::Vector3{Float64}
end

type Polygon
    bounds::Bounds2
    segments::Array{LineSegment}
end

type MeshSlice
    bounds::Bounds2
    polygons::Array{Polygon}
    layer::Float64
end

function show(io::IO, poly::Polygon)
    println("Polygon:")
    println("\tBounds:", poly.bounds)
    println("\tSegments:")
    for seg in poly.segments
        println("\t\tStart:", seg.start, " Finish:", seg.finish)
    end
end

# ##MeshSlice(*mesh::PolygonMesh, heights::Array{Float64}*)
#
# Create an array of MeshSlice at heights given by a
# monotonically increasing array of `Float64`.
#
# ----------------------------------------------------------------------------
# Parameters:
# -------------------- --------------------------------------------------
# `mesh`               The `PolygonMesh` to slice.
#
# `heights`            An array of heights at which to slice.
# ----------------------------------------------------------------------------
# ----------------------------------------------------------------------------
# Returns:
# -------------------- --------------------------------------------------
# `Array{MeshSlice}`   An `Array` of `MeshSlice` at the requested height.
# ----------------------------------------------------------------------------
function MeshSlice(mesh::PolygonMesh, heights::Array{Float64}; eps=0.00001, autoeps=true)
    slices = [LineSegment[] for i = 1:length(heights)]
    bounds = [Bounds2() for i = 1:length(heights)]

    for face in mesh.faces
        i = 1
        for height in heights
            if height > face.zmax
                break
            elseif face.zmin <= height
                seg = LineSegment(face, height)
                if seg != nothing
                    push!(slices[i], seg)
                    update!(bounds[i], seg)
                end
            end
            i = i + 1
        end
    end

    polys = MeshSlice[]

    for i = 1:length(heights)
        push!(polys, MeshSlice(bounds[i], Polygon(slices[i], eps=eps, autoeps=autoeps), heights[i]))
    end

    return polys
end

# ##Polygon()
#
# ----------------------------------------------------------------------------
# Returns:
# -------------------- --------------------------------------------------
# `Polygon`            An empty `Polygon`.
# ----------------------------------------------------------------------------
Polygon() = Polygon(Bounds2(), LineSegment[])

# ##Polygon(*lines::Array{LineSegment}, eps::Real*)
#
# Construct an `Array` of `Polygon` from an `Array` of `LineSegments`.
#
# ----------------------------------------------------------------------------
# Parameters:
# -------------------- --------------------------------------------------
# `lines`              An unorder `Array` of `LineSegment`.
#
# `eps`                The tolerance in which to count starting and finishing points as identical.
# ----------------------------------------------------------------------------
# ----------------------------------------------------------------------------
# Returns:
# -------------------- --------------------------------------------------
# `Array{Polygon}`     An `Array` of `Polygon`.
# ----------------------------------------------------------------------------
function Polygon(lines::Array{LineSegment}; eps=0.00001, autoeps=true)
    n = length(lines)
    if n == 0
        return [Polygon()]
    end
    polys = Polygon[]
    paired = [false for i = 1:n]
    start = 1
    seg = 1
    paired[seg] = true

    if autoeps
        for segment in lines
            eps = min(eps, norm(segment.start-segment.finish)/2)
        end
    end

    while true
        #Start new polygon with seg
        poly = Polygon()
        push!(poly, lines[seg])

        #Pair lines until we get to start point
        lastseg = seg
        while norm(lines[start].start - lines[seg].finish) >= eps
            lastseg = seg

            for i = 1:n
                if !paired[i]
                    if norm(lines[seg].finish - lines[i].start) <= eps
                        push!(poly, lines[i])
                        paired[i] = true
                        seg = i
                    end
                end
            end

            if (seg == start #We couldn't pair the segment
                || seg == lastseg) #The polygon can't be closed
                break
            end
        end

        if length(poly.segments) > 2
            closed = true
            if poly.segments[1].start != poly.segments[end].finish
                closed = false
            end
            for i = 1:length(poly.segments)-2
                if closed
                    break
                end
                for j = i+2:length(poly.segments)
                    if poly.segments[i].start == poly.segments[j].finish
                        poly.segments = poly.segments[i:j]
                        closed = true
                        break
                    end
                end
            end
            push!(polys,poly)
        end
        #start new polygon
        for i = 1:length(lines)
            if !paired[i] #Find next unpaired seg
                start = i
                paired[i] = true
                seg = start
                break
            elseif i == length(lines) #we have paired each segment
                return polys
            end
        end
    end
end

# ##PolygonMesh()
#
# ----------------------------------------------------------------------------
# Returns:
# -------------------- --------------------------------------------------
# `PolygonMesh`        An empty `PolygonMesh`.
# ----------------------------------------------------------------------------
PolygonMesh() = PolygonMesh(Bounds3(), Face[], Float64[])

function binarySTLvertex(file)
    return Vector3(float64(read(file, Float32)),
                   float64(read(file, Float32)),
                   float64(read(file, Float32)))
end

function binarySTLvertices(file)
    return Vector3[binarySTLvertex(file),
            binarySTLvertex(file),
            binarySTLvertex(file)]
end

# ##PolygonMesh(*path::String*)
#
# ----------------------------------------------------------------------------
# Parameters:
# -------------------- --------------------------------------------------
# `mesh`               The `PolygonMesh` to slice.
#
# `height`             The height at which to slice.
# ----------------------------------------------------------------------------
# ----------------------------------------------------------------------------
# Returns:
# -------------------- --------------------------------------------------
# `MeshSlice`          A MeshSlice at the requested height.
# ----------------------------------------------------------------------------
function PolygonMesh(path::String)

    file = open(path, "r")

    mesh = PolygonMesh()

    try
        #Discover file type and construct mesh
        #STL
        if endswith(path, ".stl")
            header = ascii(readbytes(file, 5))

            #ASCII STL
            #https://en.wikipedia.org/wiki/STL_%28file_format%29#ASCII_STL
            if lowercase(header) == "solid"
                while !eof(file)
                    line = split(lowercase(readline(file)))
                    if line[1] == "facet"
                        normal = Vector3(float64(line[3:5]))
                        readline(file) # Throw away outerloop
                        vertices = [Vector3(float64(split(readline(file))[2:4])) for i = 1:3]
                        readline(file) # throwout endloop
                        readline(file) # throwout endfacet
                        push!(mesh, Face(vertices,normal))
                    end
                end

            #Binary STL
            #https://en.wikipedia.org/wiki/STL_%28file_format%29#Binary_STL
            else
                readbytes(file, 75) # throw out header
                read(file, Uint32) # throwout triangle count
                while !eof(file)
                    normal = binarySTLvertex(file)
                    vertices = binarySTLvertices(file)
                    skip(file, 2) # throwout 16bit attribute
                    push!(mesh, Face(vertices,normal))
                end
            end
        end

    finally
        close(file)
    end

    return mesh
end

# ##rotate!(*mesh::PolygonMesh, angle::Float64, axis::Array{Float64}*, through::Array{Float64})
#
# Rotates a `PolygonMesh` around an arbitrary axis.
#
# ----------------------------------------------------------------------------
# Parameters:
# -------------------- --------------------------------------------------
# `mesh`               The `PolygonMesh` to rotate.
#
# `angle`              The rotation angle in radians.
#
# `axis`               The axis  to rotate around.
#
# `through`            The `[x, y, z]` point the `axis` passes through.
# ----------------------------------------------------------------------------
function rotate!(mesh::PolygonMesh, angle::Float64, axis::Array{Float64}, through::Array{Float64})
    axis = axis/norm(axis) # normalize axis
    a, b, c = through
    u, v, w = axis
    for face in mesh.faces
        face.vertices = [begin
                            x, y, z = face.vertices[i];
                            rotate(x, y, z, a, b, c, u, v, w, angle)
                         end for i = 1:3]
        face.zmax = max(max(face.vertices[1].e3, face.vertices[2].e3), face.vertices[3].e3)
        face.zmin = min(min(face.vertices[1].e3, face.vertices[2].e3), face.vertices[3].e3)
        update!(mesh.bounds, face)
    end
end

# ##rotate(*x, y, z, a, b, c, u, v, w, angle*)
#
# Create a MeshSlice from a PolygonMesh at a given height.
# See: [http://inside.mines.edu/~gmurray/ArbitraryAxisRotation/#x1-10011](http://inside.mines.edu/~gmurray/ArbitraryAxisRotation/#x1-10011)
#
# ----------------------------------------------------------------------------
# Parameters:
# -------------------- --------------------------------------------------
# `x, y, z`            The position to rotate.
#
# `u, v, w`            The axis to rotate around.
#
# `a, b, c`            The point to rotate through.
# ----------------------------------------------------------------------------
# ----------------------------------------------------------------------------
# Returns:
# -------------------- --------------------------------------------------
# `Vector3`            The rotation of `x, y, z`.
# ----------------------------------------------------------------------------
function rotate(x, y, z, a, b, c, u, v, w, angle)
    return Vector3((a*(v^2+w^2)-u*(b*v+c*w-u*x-v*y-w*z))*(1-cos(angle))+x*cos(angle)+(-c*v+b*w-w*y+v*z)*sin(angle),
            (b*(u^2+w^2)-v*(a*u+c*w-u*x-v*y-w*z))*(1-cos(angle))+y*cos(angle)+(c*u-a*w+w*x-u*z)*sin(angle),
            (c*(u^2+v^2)-w*(a*u+b*v-u*x-v*y-w*z))*(1-cos(angle))+z*cos(angle)+(-b*u+a*v-v*x+u*y)*sin(angle))

end

# ##push!(*mesh::PolygonMesh, f::Face*)
#
# Add a `Face` to a `PolygonMesh` object.
#
# ----------------------------------------------------------------------------
# Parameters:
# -------------------- --------------------------------------------------
# `mesh`               The `PolygonMesh` to mutate.
#
# `f`                  The `Face` to add to mesh.
# ----------------------------------------------------------------------------
function push!(mesh::PolygonMesh, f::Face)
    if abs(f.normal.e3) == 1 && !in(f.vertices[1].e3, mesh.patios)
        push!(mesh.patios,f.vertices[1].e3)
        sort!(mesh.patios) #Keep monotonic
    end
    push!(mesh.faces, f)
    update!(mesh.bounds, f)
end

function push!(poly::Polygon, seg::LineSegment)
    push!(poly.segments, seg)
    update!(poly.bounds, seg)
end

# ##LineSegment(*f::Face, z::Float64*)
#
# Get a `LineSegment` in the X-Y plane from a `Face` at a requested height.
#
# ----------------------------------------------------------------------------
# Parameters:
# -------------------- --------------------------------------------------
# `f`                  The `Face` to slice.
#
# `z`                  The height to slice the face.
# ----------------------------------------------------------------------------
# Returns:
# -------------------- --------------------------------------------------
# `LineSegment`        A `LineSegment` at the requested height in the slice plane.
# ----------------------------------------------------------------------------
function LineSegment(f::Face, z::Float64)
    p0 = f.vertices[1]
    p1 = f.vertices[2]
    p2 = f.vertices[3]

    if p0.e3 < z && p1.e3 >= z && p2.e3 >= z
        return LineSegment(p0, p2, p1, z, f.normal)
    elseif p0.e3 > z && p1.e3 < z && p2.e3 < z
        return LineSegment(p0, p1, p2, z, f.normal)
    elseif p1.e3 < z && p0.e3 >= z && p2.e3 >= z
        return LineSegment(p1, p0, p2, z, f.normal)
    elseif p1.e3 > z && p0.e3 < z && p2.e3 < z
        return LineSegment(p1, p2, p0, z, f.normal)
    elseif p2.e3 < z && p1.e3 >= z && p0.e3 >= z
        return LineSegment(p2, p1, p0, z, f.normal)
    elseif p2.e3 > z && p1.e3 < z && p0.e3 < z
        return LineSegment(p2, p0, p1, z, f.normal)
    end

end

# ##LineSegment(*p0::Vector3, p1::Vector3, p2::Vector3, z::Float64, normal::Vector3*)
#
# Get a `LineSegment` at `z` generated by the rays $\overline{p0p1}$ and $\overline{p0p2}$.
#
# ----------------------------------------------------------------------------
# Parameters:
# -------------------- --------------------------------------------------
# `p0`                 The base point.
#
# `p1`                 A point such that `z` is between `p0[3]` and `p1[3]`.
#
# `p2`                 A point such that `z` is between `p0[3]` and `p2[3]`.
#
# `z`                  The height to slice the rays $\overline{p0p1}$ and $\overline{p0p2}$.
#
# `normal`             The normal of the face, retained for insetting and offsetting.
# ----------------------------------------------------------------------------
# Returns:
# -------------------- --------------------------------------------------
# `LineSegment`        A `LineSegment` at the requested height in the slice plane.
# ----------------------------------------------------------------------------
function LineSegment(p0::Vector3, p1::Vector3, p2::Vector3, z::Float64, normal::Vector3)
    start = Vector2(p0.e1 + (p1.e1 - p0.e1) * (z - p0.e3) / (p1.e3 - p0.e3),
                    p0.e2 + (p1.e2 - p0.e2) * (z - p0.e3) / (p1.e3 - p0.e3))
    finish = Vector2(p0.e1 + (p2.e1 - p0.e1) * (z - p0.e3) / (p2.e3 - p0.e3),
                     p0.e2 + (p2.e2 - p0.e2) * (z - p0.e3) / (p2.e3 - p0.e3))
    return LineSegment(start, finish, normal);
end

# ##==(*a::LineSegment, b::LineSegment*)
#
# ----------------------------------------------------------------------------
# Returns:
# -------------------- --------------------------------------------------
# `Bool`               Whether two `LineSegment`s are equal.
# ----------------------------------------------------------------------------
function ==(a::LineSegment, b::LineSegment)
    return (a.start == b.start &&
            a.finish == b.finish)
end

# ##Bounds3()
#
# ----------------------------------------------------------------------------
# Returns:
# -------------------- --------------------------------------------------
# `Bound3`              An empty `Bounds3` of type `Float64`.
# ----------------------------------------------------------------------------
Bounds3() = Bounds3{Float64}(-Inf,-Inf,-Inf,Inf,Inf,Inf)
Bounds2() = Bounds2{Float64}(-Inf,-Inf,Inf,Inf)
# ##update!(*box::Bounds, face::Face*)
#
# ----------------------------------------------------------------------------
# Parameters:
# -------------------- --------------------------------------------------
# `box`                The `Bounds` to update.
#
# `face`               The `Face` to update the bounds against.
# ----------------------------------------------------------------------------
function update!(box::Bounds3, face::Face)
    #update the bounds against a face
    x1, y1, z1 = face.vertices[1]
    x2, y2, z2 = face.vertices[2]
    x3, y3, z3 = face.vertices[3]

    box.xmin = min(min(min(box.xmin, x1), x2), x3)
    box.ymin = min(min(min(box.ymin, y1), y2), y3)
    box.zmin = min(box.zmin, face.zmin)
    box.xmax = max(max(max(box.xmax, x1), x2), x3)
    box.ymax = max(max(max(box.ymax, y1), y2), y3)
    box.zmax = max(box.zmax, face.zmax)
end

function update!(box::Bounds2, line::LineSegment)
    #update the bounds against a line segement
    box.xmin = min(box.xmin, min(line.start[1], line.finish[1]))
    box.ymin = min(box.ymin, min(line.start[2], line.finish[2]))
    box.xmax = max(box.xmax, max(line.start[1], line.finish[1]))
    box.ymax = max(box.xmax, max(line.start[2], line.finish[2]))
end

function update!(b1::Bounds2, b2::Bounds2)
    b1.xmin = min(b1.xmin, b2.xmin)
    b1.ymin = min(b1.ymin, b2.ymin)
    b1.xmax = max(b1.xmax, b2.xmax)
    b1.ymax = max(b1.ymax, b2.ymax)
end

# ##==(*a::Bounds3, b::Bounds3*)
#
# ----------------------------------------------------------------------------
# Returns:
# -------------------- --------------------------------------------------
# `Bool`               Whether two `Bounds3` are equal.
# ----------------------------------------------------------------------------
function ==(a::Bounds3, b::Bounds3)
    return (a.xmax == b.xmax &&
            a.ymax == b.ymax &&
            a.zmax == b.zmax &&
            a.xmin == b.xmin &&
            a.ymin == b.ymin &&
            a.zmin == b.zmin)
end

function ==(a::Bounds2, b::Bounds2)
    return (a.xmax == b.xmax &&
            a.ymax == b.ymax &&
            a.xmin == b.xmin &&
            a.ymin == b.ymin)
end

end # module


text:
import Unity: UnityMesh, UnityText
using ColorTypes, GeometryTypes

#define vertices
verts = [Point3f0(0,0,0),Point3f0(1,0,0),Point3f0(1,1,0),Point3f0(0,1,0)]
#define indices for surface mesh
surf_inds = Int32[0,2,1,0,3,2]
#define indices for line mesh
line_inds = Int32[0,1,1,2,2,3,3,0]
#define indices for vertex mesh
vert_inds = Int32[0,1,2,3]
#define colors
colors = rand(RGBA{Float32},length(verts))
#colors = [RGBA{Float32}(1,1,1,.001) for vert in verts]

#define optons
options = ["surface_shader = flat"]
transparent = ["surface_shader = transparent"]
point_size = ["point_size = 0.5"]

#define text
unity_text = UnityText("Hello World!", Point3f0(1.5,2,0), Point3f0(0.1,0.1,.15), Point3f0(1,1,1), rand(RGBA{Float32}))

#create meshes
surface_mesh = UnityMesh("Mesh 1", verts, Int32[], Int32[], surf_inds, colors, options, UnityText[unity_text])
line_mesh = UnityMesh("Mesh 2", map(x->Point3f0(x[1]+3,x[2],x[3]),verts), Int32[], line_inds, Int32[], colors, String[])
point_mesh = UnityMesh("Mesh 3", map(x->Point3f0(x[1],x[2]+3,x[3]),verts), vert_inds, Int32[], Int32[], colors, point_size)
mesh = UnityMesh("Mesh 4", map(x->Point3f0(x[1]+3,x[2]+3,x[3]),verts), vert_inds, line_inds, surf_inds, colors, vcat(point_size,transparent))

#send meshes to unity
socket = connect(8052)
write(socket, surface_mesh)
write(socket, line_mesh)
write(socket, point_mesh)
write(socket, mesh)
close(socket)


text:
export FastStandardSyndrome, AbstractFastSyndrome

abstract type AbstractFastSyndrome <: SyndromeCircuit end
struct FastStandardSyndrome <: AbstractFastSyndrome end

function apply_fast_error!(noise_params::Nothing,
                           state::ChpState, zx_counts_out::Vector{Int},
                           kind::Symbol, qubits)
    # Apply no error when no noise model
end
function apply_fast_error!(noise_params::Dict{Symbol, Float64},
                           state::ChpState, zx_counts_out::Vector{Int},
                           kind::Symbol, qubits)
    p = noise_params[kind]
    p == 0 && return
    # TODO: rng
    z_count = x_count = 0
    for q in qubits
        r = rand(rng, Float64)
        r >= p && return
        # Apply X, Y, or Z with (1/3)p probability
        if r < (2/3)*p
            state.x[q] ⊻= true  # Z gate
            z_count += 1
        end
        if r >= (1/3)*p
            state.z[q] ⊻= true  # X gate
            x_count += 1
        end
    end
    zx_counts_out[2] += z_count  # Z errors cause X syndromes
    zx_counts_out[1] += x_count
    nothing
end

function apply_fast_cnot_error!(noise_params::Nothing,
                                state::ChpState,
                                zx_counts_out1::Vector{Int},
                                zx_counts_out2::Vector{Int},
                                kind::Symbol, q1::Int, q2::Int)
    # Apply no error when no noise model
end
function apply_fast_cnot_error!(noise_params::Dict{Symbol, Float64},
                                state::ChpState,
                                zx_counts_out1::Vector{Int},
                                zx_counts_out2::Vector{Int},
                                kind::Symbol, q1::Int, q2::Int)
    p = noise_params[kind]
    p == 0 && return
    # TODO: rng
    z_count1 = x_count1 = 0
    z_count2 = x_count2 = 0
    r = rand(rng, Float64)
    r >= p && return
    # Apply IX, IY, IZ, XI, XX, XY, ..., ZY, ZZ with (1/15)p probability
    c = rand(0:14)
    c1 = div(c, 4)
    c2 = mod(c, 4)
    if c1 < 2
        state.x[q1] ⊻= true  # Z gate
        z_count1 += 1
    end
    if 1 <= c1 < 3
        state.z[q1] ⊻= true  # X gate
        x_count1 += 1
    end
    if c2 < 2
        state.x[q2] ⊻= true  # Z gate
        z_count2 += 1
    end
    if 1 <= c2 < 3
        state.z[q2] ⊻= true  # X gate
        x_count2 += 1
    end
    zx_counts_out1[2] += z_count1  # Z errors cause X syndromes
    zx_counts_out1[1] += x_count1
    zx_counts_out2[2] += z_count2  # Z errors cause X syndromes
    zx_counts_out2[1] += x_count2
    nothing
end

function exec_syndrome_layer(
        noise_params::Union{Nothing, Dict{Symbol, Float64}},
        syndrome_circuit::AbstractFastSyndrome,
        run::CodeDistanceRun,
        layer_i::Int)
    ctx = run.ctx
    state = run.state
    # Apply errors
    p_data = layer_i == 1 ? :p_data_layer1 : :p_data
    for q in ctx.data_qubits
        apply_fast_error!(noise_params, state, run.zx_error_counts,
                          p_data, (q,))
    end
    # Run circuit
    for info in ctx.z_plaq_info
        anc = info.ancilla
        for dat in info.data
            state.z[anc] ⊻= state.z[dat]  # CNOT
        end
    end
    for info in ctx.x_plaq_info
        anc = info.ancilla
        for dat in info.data
            state.x[anc] ⊻= state.x[dat]  # CNOT
        end
    end
    # Apply errors
    for info in ctx.z_plaq_info
        anc = info.ancilla
        for (i, dat) in enumerate(info.data)
            apply_fast_cnot_error!(noise_params, state, run.zx_error_counts,
                                   run.zx_meas_error_counts,
                                   (i==1 ? :p_cnot1 : :p_cnot), dat, anc)
        end
    end
    for info in ctx.x_plaq_info
        anc = info.ancilla
        for (i, dat) in enumerate(info.data)
            apply_fast_cnot_error!(noise_params, state, run.zx_meas_error_counts,
                                   run.zx_error_counts,
                                   (i==1 ? :p_cnot1 : :p_cnot), anc, dat)
        end
    end
    for info in ctx.z_plaq_info
        apply_fast_error!(noise_params, state, run.zx_meas_error_counts,
                          :p_anc_z, (info.ancilla,))
    end
    for info in ctx.x_plaq_info
        apply_fast_error!(noise_params, state, run.zx_meas_error_counts,
                          :p_anc_x, (info.ancilla,))
    end
    # Measure
    for (i, info) in enumerate(ctx.z_plaq_info)
        meas = state.z[info.ancilla]  # Measure
        state.z[info.ancilla] = false  # Reset
        run.z_syndromes[i, layer_i] = run.z_prev[i] ⊻ meas
        run.z_prev[i] = meas
    end
    for (i, info) in enumerate(ctx.x_plaq_info)
        meas = state.x[info.ancilla]  # Measure
        state.x[info.ancilla] = false  # Reset
        run.x_syndromes[i, layer_i] = run.x_prev[i] ⊻ meas
        run.x_prev[i] = meas
    end
    nothing
end

function simulate_syndrome_run(syndrome_circuit::AbstractFastSyndrome,
                               run::CodeDistanceRun)
    ctx = run.ctx
    run.state.z[1:run.ctx.num_qubits] .= false  # Clean start
    run.state.x[1:run.ctx.num_qubits] .= false
    run.zx_error_counts .= 0
    run.zx_meas_error_counts .= 0
    exec_syndrome_layer(nothing, syndrome_circuit, run, 1)
    for i in 1:ctx.m_dist+1
        # Noise-free end layer
        noise = i == ctx.m_dist+1 ? nothing : run.sim_noise_params
        exec_syndrome_layer(noise, syndrome_circuit, run, i)
    end
    run.z_syndromes, run.x_syndromes
end

max_stars_repo_name: cduck/VLQ

max_issues_repo_name: cduck/VLQ

max_forks_repo_name: cduck/VLQ


text:
using StaticArrays: SHermitianCompact

export MetricTensor
export directioncosine, directionangle, distance, interplanar_spacing

struct MetricTensor{T} <: AbstractMatrix{T}
    data::SHermitianCompact{3,T,6}
end
MetricTensor(m::AbstractMatrix) = MetricTensor(SHermitianCompact{3}(m))
"""
    MetricTensor(𝐚::AbstractVector, 𝐛::AbstractVector, 𝐜::AbstractVector)

Generate a `MetricTensor` from the three basis vectors.
"""
function MetricTensor(𝐚::AbstractVector, 𝐛::AbstractVector, 𝐜::AbstractVector)
    vecs = (𝐚, 𝐛, 𝐜)
    return MetricTensor([dot(vᵢ, vⱼ) for vᵢ in vecs, vⱼ in vecs])
end
"""
    MetricTensor(lattice::Lattice)

Generate a `MetricTensor` from a `Lattice`.
"""
function MetricTensor(lattice::Lattice)
    data = lattice.data
    return MetricTensor(transpose(data) * data)
end
"""
    MetricTensor(a, b, c, α, β, γ)

Generate a `MetricTensor` from the six cell parameters.
"""
function MetricTensor(a, b, c, α, β, γ)
    g₁₂ = a * b * cosd(γ)
    g₁₃ = a * c * cosd(β)
    g₂₃ = b * c * cosd(α)
    return MetricTensor(SHermitianCompact(SVector(a^2, g₁₂, g₁₃, b^2, g₂₃, c^2)))
end
@functor MetricTensor

"""
    Lattice(g::MetricTensor)

Construct a `Lattice` from a `MetricTensor`.
"""
Lattice(g::MetricTensor) = Lattice(cellparameters(g))

"""
    cellparameters(g::MetricTensor)

Get the six cell parameters from a `MetricTensor`.
"""
function cellparameters(g::MetricTensor)
    a², b², c², ab, ac, bc = g[1, 1], g[2, 2], g[3, 3], g[1, 2], g[1, 3], g[2, 3]
    a, b, c = map(sqrt, (a², b², c²))
    γ, β, α = acosd(ab / (a * b)), acosd(ac / (a * c)), acosd(bc / (b * c))
    return a, b, c, α, β, γ
end

"""
    directioncosine(𝐚::AbstractVector, g::MetricTensor, 𝐛::AbstractVector)

Get the direction cosine of two vectors and a `MetricTensor`.
"""
directioncosine(𝐚::AbstractVector, g::MetricTensor, 𝐛::AbstractVector) =
    dot(𝐚, g, 𝐛) / (norm(𝐚, g) * norm(𝐛, g))

"""
    directionangle(𝐚::AbstractVector, g::MetricTensor, 𝐛::AbstractVector)

Get the direction angle of two vectors and a `MetricTensor`.
"""
directionangle(𝐚::AbstractVector, g::MetricTensor, 𝐛::AbstractVector) =
    acosd(directioncosine(𝐚, g, 𝐛))

"""
    distance(𝐚::AbstractVector, g::MetricTensor, 𝐛::AbstractVector)

Get the distance between two coordinates using a `MetricTensor`.
"""
distance(𝐚::AbstractVector, g::MetricTensor, 𝐛::AbstractVector) = norm(𝐛 - 𝐚, g)

"""
    interplanar_spacing(𝐚::AbstractVector, g::MetricTensor)

Get the interplanar spacing from a `MetricTensor`.
"""
interplanar_spacing(𝐚::AbstractVector, g::MetricTensor) = inv(norm(𝐚, g))

Base.size(::MetricTensor) = (3, 3)

Base.IndexStyle(::Type{<:MetricTensor}) = IndexLinear()

Base.getindex(g::MetricTensor, I::Vararg) = getindex(g.data, I...)

Base.inv(g::MetricTensor) = MetricTensor(inv(g.data))

dot(𝐚::AbstractVector, g::MetricTensor, 𝐛::AbstractVector) = 𝐚' * g.data * 𝐛
norm(𝐚::AbstractVector, g::MetricTensor) = sqrt(dot(𝐚, g, 𝐚))


text:
using UnPack, LightGraphs

"""
    LazyWeightMatrix{T,F1,F2} <: AbstractMatrix{T}

A structure that computes the weights for the graph `g`
as generated by `make_graph(M, N)` lazily, for use in shortest path algorithms.
"""
struct LazyWeightMatrix{T,F1,F2} <: AbstractMatrix{T}
    M::Int
    N::Int
    node_weight::F1
    edge_weight::F2
    nv::Int
end

function LazyWeightMatrix{T}(M, N, node_weight::F1, edge_weight::F2) where {F1,F2,T}
    nv = 2 * N * M
    return LazyWeightMatrix{T,F1,F2}(M, N, node_weight, edge_weight, nv)
end

Base.size(L::LazyWeightMatrix) = (L.nv, L.nv)

function Base.getindex(LWM::LazyWeightMatrix{T,F1,F2}, i, j) where {T,F1,F2}
    @unpack M, N, node_weight, edge_weight = LWM
    CIs = CartesianIndices((1:2, 1:M, 1:N))
    src_edge = Base.tail(Tuple(CIs[i]))
    dst_edge = Base.tail(Tuple(CIs[j]))
    if src_edge == dst_edge
        return T(node_weight(src_edge...))
    else
        return T(edge_weight(src_edge, dst_edge))
    end
    return T(Inf)
end

make_graph(M, N) = make_graph(M, N, Int)

"""
    make_graph(M, N, ::Type{IT}) where {IT <: Integer} -> SimpleDiGraph

Generates a `LightGraphs.SimpleDiGraph` corresponding to the problem.
The graph has `2*N*M` vertices: since LightGraphs does not support node weights, we split each vertex
into two, and mimic a node weight by having an edge weight between these two parts.
"""
function make_graph(M, N, ::Type{IT}) where {IT<:Integer}
    LIs = LinearIndices((1:2, 1:M, 1:N))
    IN = 1
    OUT = 2

    g = SimpleDiGraph(2 * N * M)

    # Add the "node weights" by connecting the `IN` vertex to the `OUT` vertex
    for t = 1:N, j = 1:M
        add_edge!(g, LIs[IN, j, t], LIs[OUT, j, t])
    end

    # Add the real edge weights, from each `OUT` to `IN`.
    for t = 1:N-1
        for j = 1:M, k = 1:M
            add_edge!(g, LIs[OUT, j, t], LIs[IN, k, t+1])
        end
    end

    return g
end

function get_path!(warp, M, N, τ, g, wts, algo)
    # Get the optimal path via `algo`
    a = algo(g, 1, wts)
    src = a.parents[end]

    if src == 0
        error("Could not find shortest path; problem seems to be infeasible.")
    end
    # Now we want to convert it to a time-warping path
    # by removing the doubled vertices and getting the associated
    # values of `τ`.
    CIs = CartesianIndices((1:2, 1:M, 1:N))
    dst = nv(g)
    dst_ij = Base.tail(Tuple(CIs[dst]))
    src_ij = Base.tail(Tuple(CIs[src]))
    warp[N] = τ[src_ij...]
    i = N - 1
    # We don't need to accumulate the cost
    # but we can for testing to check
    # cost = 0.0
    while true
        # cost = cost + wts[src, dst]
        if src_ij != dst_ij
            warp[i] = τ[src_ij...]
            i = i - 1
            src_ij == (1, 1) && break
        end
        dst = src
        dst_ij = src_ij
        src = a.parents[src]
        src_ij = Base.tail(Tuple(CIs[src]))
    end
    # @assert cost ≈ a.dists[end]
    return a.dists[end]
end

function DynamicAxisWarping.single_gdtw!(data::D, algo::F) where {D, F}
    @unpack warp, M, N, τ, node_weight, edge_weight = data
    g = make_graph(M, N)
    lwm = LazyWeightMatrix{Float64}(M, N, node_weight, edge_weight)
    get_path!(warp, M, N, τ, g, lwm, algo)
end


text:
module SeisNoise

using Dates, DSP, FFTW, Glob, JLD2, LinearAlgebra, SeisIO, SeisIO.Nodal
using Statistics, StatsBase, CUDA, Adapt, GPUArrays, RecipesBase

# check use of cuda
const use_cuda = Ref(false)
if !CUDA.functional()
  else
    use_cuda[] = true
end

# import types first
include("Types/NoiseData.jl")
include("Types/NodalData.jl")
include("Types/show.jl")

# import pre and post processing tools
include("distance.jl")
include("io.jl")
include("ArrayFuncs.jl")
include("tools.jl")
include("slicing.jl")
include("filter.jl")
include("phase_shift.jl")
include("stacking.jl")

# import  routines for doin' stuff
include("compute_fft.jl")
include("correlation.jl")
include("nodalcorrelation.jl")
include("rotation.jl")
include("Plotting/plotting.jl")

# out with the old
include("deprecated.jl")

end # module


text:
using IJulia
notebook()

max_stars_repo_name: pabloferz/dotfiles

max_issues_repo_name: pabloferz/dotfiles

max_forks_repo_name: pabloferz/dotfiles


text:
# This file is a part of UpROOT.jl, licensed under the MIT License (MIT).

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

isdir(testdatadir()) || mkdir(testdatadir())

for k in _testfile_keys
    download("https://github.com/scikit-hep/uproot/raw/3.5.1/tests/samples/$k.root", _testfilename(k))
end