text:
include("utils.jl")

@testset "MultinomialMatrices manifold" begin
    M = MultinomialMatrices(3, 2)
    @testset "MultinomialMatrices manifold Basics" begin
        @test ProbabilitySimplex(2)^3 === MultinomialMatrices(3, 3)
        @test ProbabilitySimplex(2)^(3,) === PowerManifold(ProbabilitySimplex(2), 3)
        @test ^(ProbabilitySimplex(2), 2) === MultinomialMatrices(3, 2)
        @test typeof(^(ProbabilitySimplex(2), 2)) == MultinomialMatrices{3,2,2}
        @test repr(M) == "MultinomialMatrices(3,2)"
        @test representation_size(M) == (3, 2)
        @test manifold_dimension(M) == 4
        @test Manifolds.get_iterator(ProbabilitySimplex(3)^4) === Base.OneTo(4)
        p = [2, 0, 0]
        p2 = [p p]
        @test !is_manifold_point(M, p)
        @test_throws DomainError is_manifold_point(M, p, true)
        @test !is_manifold_point(M, p2)
        @test_throws CompositeManifoldError is_manifold_point(M, p2, true)
        @test !is_tangent_vector(M, p2, 0.0)
        @test_throws DomainError is_tangent_vector(M, p2, [-1.0, 0.0, 0.0], true)
        @test !is_tangent_vector(M, p2, [-1.0, 0.0, 0.0])
        @test_throws DomainError is_tangent_vector(M, p, [-1.0, 0.0, 0.0], true)
        @test injectivity_radius(M) ≈ 0
        x = [0.5 0.4 0.1; 0.5 0.4 0.1]'
        y = [0.6 0.3 0.1; 0.4 0.5 0.1]'
        z = [0.3 0.6 0.1; 0.6 0.3 0.1]'
        test_manifold(
            M,
            [x, y, z],
            test_injectivity_radius = false,
            test_forward_diff = false,
            test_reverse_diff = false,
            test_vector_spaces = true,
            test_project_tangent = false,
            test_musical_isomorphisms = true,
            test_default_vector_transport = false,
            is_tangent_atol_multiplier = 5.0,
        )
    end
end


text:
function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{L, D},
    ::Type{S},
) where {
    L <: PSY.ControllableLoad,
    D <: AbstractControllablePowerLoadFormulation,
    S <: PM.AbstractPowerModel,
}
    devices = get_available_components(L, sys)

    if !validate_available_devices(L, devices)
        return
    end

    #Variables
    add_variables!(psi_container, ActivePowerVariable, devices)
    add_variables!(psi_container, ReactivePowerVariable, devices)

    #Constraints
    add_constraints!(
        psi_container,
        RangeConstraint,
        ActivePowerVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    add_constraints!(
        psi_container,
        RangeConstraint,
        ReactivePowerVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    feedforward!(psi_container, devices, model, get_feedforward(model))

    #Cost Function
    cost_function!(psi_container, devices, model, S)

    return
end

function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{L, D},
    ::Type{S},
) where {
    L <: PSY.ControllableLoad,
    D <: AbstractControllablePowerLoadFormulation,
    S <: PM.AbstractActivePowerModel,
}
    devices = get_available_components(L, sys)

    if !validate_available_devices(L, devices)
        return
    end

    #Variables
    add_variables!(psi_container, ActivePowerVariable, devices)

    #Constraints
    add_constraints!(
        psi_container,
        RangeConstraint,
        ActivePowerVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    feedforward!(psi_container, devices, model, get_feedforward(model))

    #Cost Function
    cost_function!(psi_container, devices, model, S)

    return
end

function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{L, InterruptiblePowerLoad},
    ::Type{S},
) where {L <: PSY.ControllableLoad, S <: PM.AbstractPowerModel}
    devices = get_available_components(L, sys)

    if !validate_available_devices(L, devices)
        return
    end

    #Variables
    add_variables!(psi_container, ActivePowerVariable, devices)
    add_variables!(psi_container, ReactivePowerVariable, devices)
    add_variables!(psi_container, OnVariable, devices)

    #Constraints
    add_constraints!(
        psi_container,
        RangeConstraint,
        ActivePowerVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    add_constraints!(
        psi_container,
        RangeConstraint,
        ReactivePowerVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    feedforward!(psi_container, devices, model, get_feedforward(model))

    #Cost Function
    cost_function!(psi_container, devices, model, S)

    return
end

function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{L, InterruptiblePowerLoad},
    ::Type{S},
) where {L <: PSY.ControllableLoad, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(L, sys)

    if !validate_available_devices(L, devices)
        return
    end

    #Variables
    add_variables!(psi_container, ActivePowerVariable, devices)
    add_variables!(psi_container, OnVariable, devices)

    #Constraints
    add_constraints!(
        psi_container,
        RangeConstraint,
        ActivePowerVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    feedforward!(psi_container, devices, model, get_feedforward(model))

    #Cost Function
    cost_function!(psi_container, devices, model, S)

    return
end

function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{L, StaticPowerLoad},
    ::Type{S},
) where {L <: PSY.ElectricLoad, S <: PM.AbstractPowerModel}
    devices = get_available_components(L, sys)

    if !validate_available_devices(L, devices)
        return
    end

    nodal_expression!(psi_container, devices, S)

    return
end

function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{L, D},
    ::Type{S},
) where {
    L <: PSY.StaticLoad,
    D <: AbstractControllablePowerLoadFormulation,
    S <: PM.AbstractPowerModel,
}
    if D != StaticPowerLoad
        @warn("The Formulation $(D) only applies to FormulationControllable Loads, \n Consider Changing the Device Formulation to StaticPowerLoad")
    end

    construct_device!(psi_container, sys, DeviceModel(L, StaticPowerLoad), S)
    return
end

ext: jl


text:
"""
# Percent Change

```julia
pctchange(ts::TS, periods::Int = 1)
```

Return the percentage change between successive row elements.
Default is the element in the next row. `periods` defines the number
of rows to be shifted over. The skipped rows are rendered as `missing`.

`pctchange` returns an error if column type does not have the method `/`.

# Examples
```jldoctest; setup = :(using TSx, DataFrames, Dates, Random, Statistics)
julia> using Random, Statistics;

julia> random(x) = rand(MersenneTwister(123), x);

julia> dates = collect(Date(2017,1,1):Day(1):Date(2017,1,10));

julia> ts = TS(random(length(dates)), dates)
julia> show(ts)
(10 x 1) TS with Date Index

 Index       x1
 Date        Float64
───────────────────────
 2017-01-01  0.768448
 2017-01-02  0.940515
 2017-01-03  0.673959
 2017-01-04  0.395453
 2017-01-05  0.313244
 2017-01-06  0.662555
 2017-01-07  0.586022
 2017-01-08  0.0521332
 2017-01-09  0.26864
 2017-01-10  0.108871

# Pctchange over successive rows
julia> pctchange(ts)
(10 x 1) TS with Date Index

 Index       x1
 Date        Float64?
────────────────────────────
 2017-01-01  missing
 2017-01-02        0.223915
 2017-01-03       -0.283415
 2017-01-04       -0.413238
 2017-01-05       -0.207886
 2017-01-06        1.11514
 2017-01-07       -0.115511
 2017-01-08       -0.911039
 2017-01-09        4.15295
 2017-01-10       -0.594733

# Pctchange over the third row
julia> pctchange(ts, 3)
(10 x 1) TS with Date Index

 Index       x1
 Date        Float64?
─────────────────────────────
 2017-01-01  missing
 2017-01-02  missing
 2017-01-03  missing
 2017-01-04       -0.485387
 2017-01-05       -0.666944
 2017-01-06       -0.0169207
 2017-01-07        0.4819
 2017-01-08       -0.83357
 2017-01-09       -0.59454
 2017-01-10       -0.814221

```
"""

# Pctchange
function pctchange(ts::TS, periods::Int = 1)
    if periods <= 0
        error("periods must be a positive int")
    end
    ddf = (ts.coredata[:, Not(:Index)] ./ TSx.lag(ts, periods).coredata[:, Not(:Index)]) .- 1
    insertcols!(ddf, 1, "Index" => ts.coredata[:, :Index])
    TS(ddf, :Index)
end


text:
function divide_kernel!(xs, ys, max_idx)
    index = threadIdx().x + (blockIdx().x - 1) * blockDim().x

    @inbounds if index <= max_idx
        xs[index] = xs[index] / ys[index]
    end
    return nothing
end

function divide_kernel!(xs, counts, max_idx, max_dims_idx, dims_size)
    index = threadIdx().x + (blockIdx().x - 1) * blockDim().x

    @inbounds if index <= max_idx
        j, k = divrem(index-1, max_dims_idx)
        dims_i = Tuple(CartesianIndices(dims_size)[k+1])
        CUDA.@atomic xs[dims_i..., j+1] = xs[dims_i..., j+1] / counts[j+1]
    end
    return nothing
end

function NNlib.reverse_indices(idx::AnyCuArray{<:Any,N}) where N
    max_dims = NNlib.maximum_dims(idx)
    T = CartesianIndex{N}
    rev = Array{Vector{T}}(undef, max_dims...)
    for i in eachindex(rev)
        rev[i] = T[]
    end
    NNlib.reverse_indices!(rev, idx)
    return map(cu, rev)
end


text:
struct Staggered_Dirac_operator{Dim,T,fermion} <: Dirac_operator{Dim}  where T <: AbstractGaugefields
    U::Array{T,1}
    boundarycondition::Vector{Int8}
    mass::Float64
    _temporary_fermi::Vector{fermion}
    eps_CG::Float64
    MaxCGstep::Int64
    verbose_level::Int8
    method_CG::String
    verbose_print::Verbose_print
    _temporary_fermion_forCG::Vector{fermion}
    #verbose::Union{Verbose_1,Verbose_2,Verbose_3}
end

include("./StaggeredFermion_4D_wing.jl")
include("./StaggeredFermion_2D_wing.jl")

function Staggered_Dirac_operator(U::Array{<: AbstractGaugefields{NC,Dim},1},x,parameters) where {NC,Dim}
    xtype = typeof(x)
    num = 6
    _temporary_fermi = Array{xtype,1}(undef,num)

    @assert haskey(parameters,"mass") "parameters should have the keyword mass"
    mass = parameters["mass"]
    if Dim == 4
        boundarycondition = check_parameters(parameters,"boundarycondition",[1,1,1,-1])
    elseif Dim == 2
        boundarycondition = check_parameters(parameters,"boundarycondition",[1,-1])
    else
        error("Dim should be 2 or 4!")
    end
    eps_CG = check_parameters(parameters,"eps",default_eps_CG)
    MaxCGstep = check_parameters(parameters,"MaxCGstep",default_MaxCGstep)

    verbose_level = check_parameters(parameters,"verbose_level",2)

    method_CG = check_parameters(parameters,"method_CG","bicg")

    for i=1:num
        _temporary_fermi[i] = similar(x)
    end

    numcg = 7
    _temporary_fermion_forCG= Array{xtype,1}(undef,numcg)
    for i=1:numcg
        _temporary_fermion_forCG[i] = similar(x)
    end

    verbose_print = Verbose_print(verbose_level)

    return Staggered_Dirac_operator{Dim,eltype(U),xtype }(U,boundarycondition,mass,_temporary_fermi,
        eps_CG,MaxCGstep,verbose_level,method_CG,verbose_print,_temporary_fermion_forCG)
end

function (D::Staggered_Dirac_operator{Dim,T,fermion })(U) where {Dim,T,fermion}
    return Staggered_Dirac_operator{Dim,T,fermion}(
        U,D.boundarycondition,D.mass,D._temporary_fermi,D.eps_CG,D.MaxCGstep,D.verbose_level,
        D.method_CG,D.verbose_print,D._temporary_fermion_forCG)
end

struct DdagD_Staggered_operator{Dim,T,fermion} <: DdagD_operator
    dirac::Staggered_Dirac_operator{Dim,T,fermion}

    function DdagD_Staggered_operator(U::Array{T,1},x,parameters) where  T <: AbstractGaugefields
        return new{Dim,T,fermion}(Staggered_Dirac_operator(U,x,parameters))
    end

    function DdagD_Staggered_operator(D::Staggered_Dirac_operator{Dim,T,fermion}) where  {Dim,T,fermion}
        return new{Dim,T,fermion}(D)
    end

end

struct Adjoint_Staggered_operator{T} <: Adjoint_Dirac_operator
    parent::T
end

function Base.adjoint(A::T) where T <: Staggered_Dirac_operator
    Adjoint_Staggered_operator{typeof(A)}(A)
end

function Initialize_StaggeredFermion(u::AbstractGaugefields{NC,Dim}) where {NC,Dim}
    _,_,NN... = size(u)
    return Initialize_StaggeredFermion(NC,NN...)
end

function Initialize_StaggeredFermion(NC,NN...)
    Dim = length(NN)
    if Dim == 4
        fermion = StaggeredFermion_4D_wing(NC,NN...)
    elseif Dim == 2
        fermion = StaggeredFermion_2D_wing(NC,NN...)
    else
        error("Dimension $Dim is not supported")
    end
    return fermion
end

function LinearAlgebra.mul!(y::T1,A::T2,x::T3) where {T1 <:AbstractFermionfields ,T2 <: Staggered_Dirac_operator, T3 <:AbstractFermionfields}

    @assert typeof(A._temporary_fermi[1]) == typeof(x) "type of A._temporary_fermi[1] $(typeof(A._temporary_fermi[1])) should be type of x: $(typeof(x))"
    temps = A._temporary_fermi
    temp = temps[4]
    Dx!(temp,A.U,x,[temps[1],temps[2],temps[3]],A.boundarycondition)
    clear_fermion!(y)
    add_fermion!(y,A.mass,x,1,temp)
    set_wing_fermion!(y,A.boundarycondition)

    #error("LinearAlgebra.mul!(y,A,x) is not implemented in type y:$(typeof(y)),A:$(typeof(A)) and x:$(typeof(x))")
end

function LinearAlgebra.mul!(y::T1,A::T2,x::T3) where {T1 <:AbstractFermionfields,T2 <: Adjoint_Staggered_operator, T3 <:  AbstractFermionfields}
    #error("LinearAlgebra.mul!(y,A,x) is not implemented in type y:$(typeof(y)),A:$(typeof(A)) and x:$(typeof(x))")

    temps = A.parent._temporary_fermi
    temp = temps[4]
    Dx!(temp,A.parent.U,x,[temps[1],temps[2],temps[3]],A.parent.boundarycondition)
    clear_fermion!(y)
    add_fermion!(y,A.parent.mass,x,-1,temp)
    set_wing_fermion!(y,A.parent.boundarycondition)
    #println(y[1,1,1,1,1,1])
    return
end

function LinearAlgebra.mul!(y::T1,A::T2,x::T3) where {T1 <:AbstractFermionfields,T2 <:  DdagD_Staggered_operator, T3 <:  AbstractFermionfields}
    @assert typeof(A.dirac._temporary_fermi[1]) == typeof(x) "type of A._temporary_fermi[1] $(typeof(A.dirac._temporary_fermi[1])) should be type of x: $(typeof(x))"
    temps = A.dirac._temporary_fermi
    temp = temps[5]
    temp2 = temps[6]
    Dx!(temp,A.dirac.U,x,[temps[1],temps[2],temps[3]],A.dirac.boundarycondition)
    Dx!(temp2,A.dirac.U,temp,[temps[1],temps[2],temps[3]],A.dirac.boundarycondition)

    clear_fermion!(y)
    add_fermion!(y,A.dirac.mass^2,x,-1,temp2)
    set_wing_fermion!(y,A.dirac.boundarycondition)

    #error("xout")
    return
end


text:
module CuArrays

using CUDAdrv, CUDAnative

using GPUArrays

export CuArray, CuVector, CuMatrix, CuVecOrMat, cu, cuzeros, cuones, cufill

import LinearAlgebra

using Adapt

const ext = joinpath(dirname(@__DIR__), "deps", "ext.jl")
isfile(ext) || error("CuArrays.jl has not been built, please run Pkg.build(\"CuArrays\").")
include(ext)
if !configured
    # default (non-functional) values for critical variables,
    # making it possible to _load_ the package at all times.
    const libcublas = nothing
    const libcusparse = nothing
    const libcusolver = nothing
    const libcufft = nothing
    const libcurand = nothing
    const libcudnn = nothing
end

include("memory.jl")
include("array.jl")
include("subarray.jl")
include("utils.jl")
include("indexing.jl")
include("broadcast.jl")
include("matmul.jl")
include("mapreduce.jl")

include("gpuarray_interface.jl")

# many libraries need to be initialized per-device (per-context, really, but we assume users
# of CuArrays and/or CUDAnative only use a single context), so keep track of the active one.
const active_context = Ref{CuContext}()

libcublas !== nothing   && include("blas/CUBLAS.jl")
libcusparse !== nothing && include("sparse/CUSPARSE.jl")
libcusolver !== nothing && include("solver/CUSOLVER.jl")
libcufft !== nothing    && include("fft/CUFFT.jl")
libcurand !== nothing   && include("rand/CURAND.jl")
libcudnn !== nothing    && include("dnn/CUDNN.jl")

include("nnlib.jl")

include("deprecated.jl")

function __init__()
    if !configured
        @warn("CuArrays.jl has not been successfully built, and will not work properly.")
        @warn("Please run Pkg.build(\"CuArrays\") and restart Julia.")
        return
    end

    function check_library(name, path)
        path === nothing && return
        if !ispath(path)
            error("$name library has changed. Please run Pkg.build(\"CuArrays\") and restart Julia.")
        end
    end
    check_library("CUBLAS", libcublas)
    check_library("CUSPARSE", libcusparse)
    check_library("CUSOLVER", libcusolver)
    check_library("CUFFT", libcufft)
    check_library("CURAND", libcurand)
    check_library("CUDNN", libcudnn)

    # update the active context when we switch devices
    callback = (::CuDevice, ctx::CuContext) -> begin
        active_context[] = ctx

        # wipe the active handles
        isdefined(CuArrays, :CUBLAS)   && (CUBLAS._handle[] = C_NULL)
        isdefined(CuArrays, :CUSOLVER) && (CUSOLVER._dense_handle[] = C_NULL)
        isdefined(CuArrays, :CURAND)   && (CURAND._generator[] = nothing)
        isdefined(CuArrays, :CUDNN)    && (CUDNN._handle[] = C_NULL)
    end
    push!(CUDAnative.device!_listeners, callback)

    # a device might be active already
    existing_ctx = CUDAdrv.CuCurrentContext()
    if existing_ctx !== nothing
        active_context[] = existing_ctx
    end

    __init_memory__()
end

end # module


text:
function FieldUnivariateDerivative1(arg0::FieldDerivativeStructure)
    return FieldUnivariateDerivative1((FieldDerivativeStructure,), arg0)
end

function FieldUnivariateDerivative1(arg0::RealFieldElement, arg1::RealFieldElement)
    return FieldUnivariateDerivative1((RealFieldElement, RealFieldElement), arg0, arg1)
end

function abs(obj::FieldUnivariateDerivative1)
    return jcall(obj, "abs", FieldUnivariateDerivative1, ())
end

function acos(obj::FieldUnivariateDerivative1)
    return jcall(obj, "acos", FieldUnivariateDerivative1, ())
end

function acosh(obj::FieldUnivariateDerivative1)
    return jcall(obj, "acosh", FieldUnivariateDerivative1, ())
end

function add(obj::FieldUnivariateDerivative1, arg0::FieldUnivariateDerivative1)
    return jcall(obj, "add", FieldUnivariateDerivative1, (FieldUnivariateDerivative1,), arg0)
end

function add(obj::FieldUnivariateDerivative1, arg0::RealFieldElement)
    return jcall(obj, "add", FieldUnivariateDerivative1, (RealFieldElement,), arg0)
end

function add(obj::FieldUnivariateDerivative1, arg0::jdouble)
    return jcall(obj, "add", FieldUnivariateDerivative1, (jdouble,), arg0)
end

function asin(obj::FieldUnivariateDerivative1)
    return jcall(obj, "asin", FieldUnivariateDerivative1, ())
end

function asinh(obj::FieldUnivariateDerivative1)
    return jcall(obj, "asinh", FieldUnivariateDerivative1, ())
end

function atan(obj::FieldUnivariateDerivative1)
    return jcall(obj, "atan", FieldUnivariateDerivative1, ())
end

function atan2(obj::FieldUnivariateDerivative1, arg0::FieldUnivariateDerivative1)
    return jcall(obj, "atan2", FieldUnivariateDerivative1, (FieldUnivariateDerivative1,), arg0)
end

function atanh(obj::FieldUnivariateDerivative1)
    return jcall(obj, "atanh", FieldUnivariateDerivative1, ())
end

function cbrt(obj::FieldUnivariateDerivative1)
    return jcall(obj, "cbrt", FieldUnivariateDerivative1, ())
end

function ceil(obj::FieldUnivariateDerivative1)
    return jcall(obj, "ceil", FieldUnivariateDerivative1, ())
end

function compose(obj::FieldUnivariateDerivative1, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(obj, "compose", FieldUnivariateDerivative1, (RealFieldElement, RealFieldElement), arg0, arg1)
end

function copy_sign(obj::FieldUnivariateDerivative1, arg0::FieldUnivariateDerivative1)
    return jcall(obj, "copySign", FieldUnivariateDerivative1, (FieldUnivariateDerivative1,), arg0)
end

function copy_sign(obj::FieldUnivariateDerivative1, arg0::RealFieldElement)
    return jcall(obj, "copySign", FieldUnivariateDerivative1, (RealFieldElement,), arg0)
end

function copy_sign(obj::FieldUnivariateDerivative1, arg0::jdouble)
    return jcall(obj, "copySign", FieldUnivariateDerivative1, (jdouble,), arg0)
end

function cos(obj::FieldUnivariateDerivative1)
    return jcall(obj, "cos", FieldUnivariateDerivative1, ())
end

function cosh(obj::FieldUnivariateDerivative1)
    return jcall(obj, "cosh", FieldUnivariateDerivative1, ())
end

function divide(obj::FieldUnivariateDerivative1, arg0::FieldUnivariateDerivative1)
    return jcall(obj, "divide", FieldUnivariateDerivative1, (FieldUnivariateDerivative1,), arg0)
end

function divide(obj::FieldUnivariateDerivative1, arg0::RealFieldElement)
    return jcall(obj, "divide", FieldUnivariateDerivative1, (RealFieldElement,), arg0)
end

function divide(obj::FieldUnivariateDerivative1, arg0::jdouble)
    return jcall(obj, "divide", FieldUnivariateDerivative1, (jdouble,), arg0)
end

function equals(obj::FieldUnivariateDerivative1, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function exp(obj::FieldUnivariateDerivative1)
    return jcall(obj, "exp", FieldUnivariateDerivative1, ())
end

function expm1(obj::FieldUnivariateDerivative1)
    return jcall(obj, "expm1", FieldUnivariateDerivative1, ())
end

function floor(obj::FieldUnivariateDerivative1)
    return jcall(obj, "floor", FieldUnivariateDerivative1, ())
end

function get_derivative(obj::FieldUnivariateDerivative1, arg0::jint)
    return jcall(obj, "getDerivative", RealFieldElement, (jint,), arg0)
end

function get_exponent(obj::FieldUnivariateDerivative1)
    return jcall(obj, "getExponent", jint, ())
end

function get_field(obj::FieldUnivariateDerivative1)
    return jcall(obj, "getField", FieldUnivariateDerivative1Field, ())
end

function get_first_derivative(obj::FieldUnivariateDerivative1)
    return jcall(obj, "getFirstDerivative", RealFieldElement, ())
end

function get_order(obj::FieldUnivariateDerivative1)
    return jcall(obj, "getOrder", jint, ())
end

function get_real(obj::FieldUnivariateDerivative1)
    return jcall(obj, "getReal", jdouble, ())
end

function get_value(obj::FieldUnivariateDerivative1)
    return jcall(obj, "getValue", RealFieldElement, ())
end

function get_value_field(obj::FieldUnivariateDerivative1)
    return jcall(obj, "getValueField", Field, ())
end

function hash_code(obj::FieldUnivariateDerivative1)
    return jcall(obj, "hashCode", jint, ())
end

function hypot(obj::FieldUnivariateDerivative1, arg0::FieldUnivariateDerivative1)
    return jcall(obj, "hypot", FieldUnivariateDerivative1, (FieldUnivariateDerivative1,), arg0)
end

function linear_combination(obj::FieldUnivariateDerivative1, arg0::FieldUnivariateDerivative1, arg1::FieldUnivariateDerivative1, arg2::FieldUnivariateDerivative1, arg3::FieldUnivariateDerivative1)
    return jcall(obj, "linearCombination", FieldUnivariateDerivative1, (FieldUnivariateDerivative1, FieldUnivariateDerivative1, FieldUnivariateDerivative1, FieldUnivariateDerivative1), arg0, arg1, arg2, arg3)
end

function linear_combination(obj::FieldUnivariateDerivative1, arg0::FieldUnivariateDerivative1, arg1::FieldUnivariateDerivative1, arg2::FieldUnivariateDerivative1, arg3::FieldUnivariateDerivative1, arg4::FieldUnivariateDerivative1, arg5::FieldUnivariateDerivative1)
    return jcall(obj, "linearCombination", FieldUnivariateDerivative1, (FieldUnivariateDerivative1, FieldUnivariateDerivative1, FieldUnivariateDerivative1, FieldUnivariateDerivative1, FieldUnivariateDerivative1, FieldUnivariateDerivative1), arg0, arg1, arg2, arg3, arg4, arg5)
end

function linear_combination(obj::FieldUnivariateDerivative1, arg0::FieldUnivariateDerivative1, arg1::FieldUnivariateDerivative1, arg2::FieldUnivariateDerivative1, arg3::FieldUnivariateDerivative1, arg4::FieldUnivariateDerivative1, arg5::FieldUnivariateDerivative1, arg6::FieldUnivariateDerivative1, arg7::FieldUnivariateDerivative1)
    return jcall(obj, "linearCombination", FieldUnivariateDerivative1, (FieldUnivariateDerivative1, FieldUnivariateDerivative1, FieldUnivariateDerivative1, FieldUnivariateDerivative1, FieldUnivariateDerivative1, FieldUnivariateDerivative1, FieldUnivariateDerivative1, FieldUnivariateDerivative1), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function linear_combination(obj::FieldUnivariateDerivative1, arg0::RealFieldElement, arg1::FieldUnivariateDerivative1, arg2::RealFieldElement, arg3::FieldUnivariateDerivative1, arg4::RealFieldElement, arg5::FieldUnivariateDerivative1)
    return jcall(obj, "linearCombination", FieldUnivariateDerivative1, (RealFieldElement, FieldUnivariateDerivative1, RealFieldElement, FieldUnivariateDerivative1, RealFieldElement, FieldUnivariateDerivative1), arg0, arg1, arg2, arg3, arg4, arg5)
end

function linear_combination(obj::FieldUnivariateDerivative1, arg0::Vector{FieldUnivariateDerivative1}, arg1::Vector{FieldUnivariateDerivative1})
    return jcall(obj, "linearCombination", FieldUnivariateDerivative1, (Vector{FieldUnivariateDerivative1}, Vector{FieldUnivariateDerivative1}), arg0, arg1)
end

function linear_combination(obj::FieldUnivariateDerivative1, arg0::Vector{RealFieldElement}, arg1::Vector{FieldUnivariateDerivative1})
    return jcall(obj, "linearCombination", FieldUnivariateDerivative1, (Vector{RealFieldElement}, Vector{FieldUnivariateDerivative1}), arg0, arg1)
end

function linear_combination(obj::FieldUnivariateDerivative1, arg0::Vector{jdouble}, arg1::Vector{FieldUnivariateDerivative1})
    return jcall(obj, "linearCombination", FieldUnivariateDerivative1, (Vector{jdouble}, Vector{FieldUnivariateDerivative1}), arg0, arg1)
end

function linear_combination(obj::FieldUnivariateDerivative1, arg0::jdouble, arg1::FieldUnivariateDerivative1, arg2::jdouble, arg3::FieldUnivariateDerivative1)
    return jcall(obj, "linearCombination", FieldUnivariateDerivative1, (jdouble, FieldUnivariateDerivative1, jdouble, FieldUnivariateDerivative1), arg0, arg1, arg2, arg3)
end

function linear_combination(obj::FieldUnivariateDerivative1, arg0::jdouble, arg1::FieldUnivariateDerivative1, arg2::jdouble, arg3::FieldUnivariateDerivative1, arg4::jdouble, arg5::FieldUnivariateDerivative1)
    return jcall(obj, "linearCombination", FieldUnivariateDerivative1, (jdouble, FieldUnivariateDerivative1, jdouble, FieldUnivariateDerivative1, jdouble, FieldUnivariateDerivative1), arg0, arg1, arg2, arg3, arg4, arg5)
end

function linear_combination(obj::FieldUnivariateDerivative1, arg0::jdouble, arg1::FieldUnivariateDerivative1, arg2::jdouble, arg3::FieldUnivariateDerivative1, arg4::jdouble, arg5::FieldUnivariateDerivative1, arg6::jdouble, arg7::FieldUnivariateDerivative1)
    return jcall(obj, "linearCombination", FieldUnivariateDerivative1, (jdouble, FieldUnivariateDerivative1, jdouble, FieldUnivariateDerivative1, jdouble, FieldUnivariateDerivative1, jdouble, FieldUnivariateDerivative1), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function log(obj::FieldUnivariateDerivative1)
    return jcall(obj, "log", FieldUnivariateDerivative1, ())
end

function log10(obj::FieldUnivariateDerivative1)
    return jcall(obj, "log10", FieldUnivariateDerivative1, ())
end

function log1p(obj::FieldUnivariateDerivative1)
    return jcall(obj, "log1p", FieldUnivariateDerivative1, ())
end

function multiply(obj::FieldUnivariateDerivative1, arg0::FieldUnivariateDerivative1)
    return jcall(obj, "multiply", FieldUnivariateDerivative1, (FieldUnivariateDerivative1,), arg0)
end

function multiply(obj::FieldUnivariateDerivative1, arg0::RealFieldElement)
    return jcall(obj, "multiply", FieldUnivariateDerivative1, (RealFieldElement,), arg0)
end

function multiply(obj::FieldUnivariateDerivative1, arg0::jdouble)
    return jcall(obj, "multiply", FieldUnivariateDerivative1, (jdouble,), arg0)
end

function multiply(obj::FieldUnivariateDerivative1, arg0::jint)
    return jcall(obj, "multiply", FieldUnivariateDerivative1, (jint,), arg0)
end

function negate(obj::FieldUnivariateDerivative1)
    return jcall(obj, "negate", FieldUnivariateDerivative1, ())
end

function new_instance(obj::FieldUnivariateDerivative1, arg0::jdouble)
    return jcall(obj, "newInstance", FieldUnivariateDerivative1, (jdouble,), arg0)
end

function pow(::Type{FieldUnivariateDerivative1}, arg0::jdouble, arg1::FieldUnivariateDerivative1)
    return jcall(FieldUnivariateDerivative1, "pow", FieldUnivariateDerivative1, (jdouble, FieldUnivariateDerivative1), arg0, arg1)
end

function pow(obj::FieldUnivariateDerivative1, arg0::FieldUnivariateDerivative1)
    return jcall(obj, "pow", FieldUnivariateDerivative1, (FieldUnivariateDerivative1,), arg0)
end

function pow(obj::FieldUnivariateDerivative1, arg0::jdouble)
    return jcall(obj, "pow", FieldUnivariateDerivative1, (jdouble,), arg0)
end

function pow(obj::FieldUnivariateDerivative1, arg0::jint)
    return jcall(obj, "pow", FieldUnivariateDerivative1, (jint,), arg0)
end

function reciprocal(obj::FieldUnivariateDerivative1)
    return jcall(obj, "reciprocal", FieldUnivariateDerivative1, ())
end

function remainder(obj::FieldUnivariateDerivative1, arg0::FieldUnivariateDerivative1)
    return jcall(obj, "remainder", FieldUnivariateDerivative1, (FieldUnivariateDerivative1,), arg0)
end

function remainder(obj::FieldUnivariateDerivative1, arg0::RealFieldElement)
    return jcall(obj, "remainder", FieldUnivariateDerivative1, (RealFieldElement,), arg0)
end

function remainder(obj::FieldUnivariateDerivative1, arg0::jdouble)
    return jcall(obj, "remainder", FieldUnivariateDerivative1, (jdouble,), arg0)
end

function rint(obj::FieldUnivariateDerivative1)
    return jcall(obj, "rint", FieldUnivariateDerivative1, ())
end

function root_n(obj::FieldUnivariateDerivative1, arg0::jint)
    return jcall(obj, "rootN", FieldUnivariateDerivative1, (jint,), arg0)
end

function scalb(obj::FieldUnivariateDerivative1, arg0::jint)
    return jcall(obj, "scalb", FieldUnivariateDerivative1, (jint,), arg0)
end

function signum(obj::FieldUnivariateDerivative1)
    return jcall(obj, "signum", FieldUnivariateDerivative1, ())
end

function sin(obj::FieldUnivariateDerivative1)
    return jcall(obj, "sin", FieldUnivariateDerivative1, ())
end

function sin_cos(obj::FieldUnivariateDerivative1)
    return jcall(obj, "sinCos", FieldSinCos, ())
end

function sinh(obj::FieldUnivariateDerivative1)
    return jcall(obj, "sinh", FieldUnivariateDerivative1, ())
end

function sqrt(obj::FieldUnivariateDerivative1)
    return jcall(obj, "sqrt", FieldUnivariateDerivative1, ())
end

function subtract(obj::FieldUnivariateDerivative1, arg0::FieldUnivariateDerivative1)
    return jcall(obj, "subtract", FieldUnivariateDerivative1, (FieldUnivariateDerivative1,), arg0)
end

function subtract(obj::FieldUnivariateDerivative1, arg0::RealFieldElement)
    return jcall(obj, "subtract", FieldUnivariateDerivative1, (RealFieldElement,), arg0)
end

function subtract(obj::FieldUnivariateDerivative1, arg0::jdouble)
    return jcall(obj, "subtract", FieldUnivariateDerivative1, (jdouble,), arg0)
end

function tan(obj::FieldUnivariateDerivative1)
    return jcall(obj, "tan", FieldUnivariateDerivative1, ())
end

function tanh(obj::FieldUnivariateDerivative1)
    return jcall(obj, "tanh", FieldUnivariateDerivative1, ())
end

function taylor(obj::FieldUnivariateDerivative1, arg0::RealFieldElement)
    return jcall(obj, "taylor", RealFieldElement, (RealFieldElement,), arg0)
end

function taylor(obj::FieldUnivariateDerivative1, arg0::jdouble)
    return jcall(obj, "taylor", RealFieldElement, (jdouble,), arg0)
end

function to_degrees(obj::FieldUnivariateDerivative1)
    return jcall(obj, "toDegrees", FieldUnivariateDerivative1, ())
end

function to_derivative_structure(obj::FieldUnivariateDerivative1)
    return jcall(obj, "toDerivativeStructure", FieldDerivativeStructure, ())
end

function to_radians(obj::FieldUnivariateDerivative1)
    return jcall(obj, "toRadians", FieldUnivariateDerivative1, ())
end


text:
using Dates: AbstractTime, Second, value

# Note: The IANA time zone database rounds offset precision to the nearest second
# See "America/New_York" notes in tzdata file "northamerica" for an example.
"""
    UTCOffset

A `UTCOffset` is an amount of time subtracted from or added to UTC to get the current
local time – whether it's standard time or daylight saving time.
"""
struct UTCOffset <: AbstractTime
    std::Second  # Standard time offset from UTC in seconds
    dst::Second  # Daylight saving time offset in seconds

    function UTCOffset(std_offset::Second, dst_offset::Second=Second(0))
        new(std_offset, dst_offset)
    end
end

function UTCOffset(std_offset::Integer, dst_offset::Integer=0)
    UTCOffset(Second(std_offset), Second(dst_offset))
end

Dates.value(offset::UTCOffset) = value(offset.std + offset.dst)

Base.:(+)(dt::DateTime, offset::UTCOffset) = dt + (offset.std + offset.dst)
Base.:(-)(dt::DateTime, offset::UTCOffset) = dt - (offset.std + offset.dst)
Base.:(-)(a::UTCOffset, b::UTCOffset) = UTCOffset(a.std - b.std, a.dst - b.dst)

# Determines if the given `UTCOffset` is an offset for daylight saving time.
isdst(offset::UTCOffset) = offset.dst != Second(0)

# Two `UTCOffset`s can be considered equal if the total offset is the same and they are
# both either offsets for standard time or daylight saving time.
function Base.isequal(x::UTCOffset, y::UTCOffset)
    x == y || value(x) == value(y) && isdst(x) == isdst(y)
end
Base.isless(x::UTCOffset, y::UTCOffset) = isless(value(x), value(y))

function offset_string(seconds::Second, iso8601::Bool=false)
    val = value(seconds)
    sig = val < 0 ? '-' : '+'
    hour, val = divrem(abs(val), 3600)
    minute, second  = divrem(val, 60)

    if !iso8601 && minute == 0 && second == 0
        return @sprintf("%c%01d", sig, hour)
    elseif second == 0
        return @sprintf("%c%02d:%02d", sig, hour, minute)
    else
        # Not in the ISO 8601 standard
        return @sprintf("%c%02d:%02d:%02d", sig, hour, minute, second)
    end
end
function offset_string(offset::UTCOffset, iso8601::Bool=false)
    offset_string(offset.std + offset.dst, iso8601)
end

Base.print(io::IO, o::UTCOffset) = print(io, offset_string(o, true))

function Base.show(io::IO, o::UTCOffset)
    if get(io, :compact, false)
        # Show DST as a separate offset since we want to distinguish between normal hourly
        # daylight saving time offsets and exotic DST offsets (e.g. midsummer time).
        print(io, "UTC", offset_string(o.std), "/", offset_string(o.dst))
    else
        # Fallback to calling the default show instead of reimplementing it.
        invoke(show, Tuple{IO, Any}, io, o)
    end
end

function Base.show(io::IO, ::MIME"text/plain", o::UTCOffset)
    show(IOContext(io, :compact => true), o)
end


text:
using
    MaxPlus, Test, Base, LinearAlgebra, SparseArrays

include("test_maxplus.jl")
#include("test_syslin.jl")
#include("test_flowshop.jl")


text:
# Banded Needleman-Wunsch Algorithm
# =================================
#
# Banded counterpart of the Needleman-Wunsch algorithm.
#
# The banded Needleman-Wunsch algorithm is similar to the normal
# Needleman-Wunsch algorithm but search space is limited to specific cells that
# form a lower and upper bounded band in the dynamic programming matrix.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

type BandedNeedlemanWunsch{T<:Union{Signed,AbstractFloat}}
    trace::Matrix{Trace}
    H::Vector{T}
    E::Vector{T}
    lower::Int
    upper::Int

    function BandedNeedlemanWunsch(
            m::Integer, n::Integer,
            lower::Integer, upper::Integer)
        lower = min(lower, m)
        upper = min(upper, n)
        width = lower + upper + 1
        trace = Matrix{Trace}(width, n + 1)
        H = Vector{T}(width)
        E = Vector{T}(width)
        return new(trace, H, E, lower, upper)
    end

    function BandedNeedlemanWunsch()
        return BandedNeedlemanWunsch{T}(0, 0, 0, 0)
    end
end

function ensureroom!(nw::BandedNeedlemanWunsch, m, n, lower, upper)
    lower = min(lower, m)
    upper = min(upper, n)
    width = lower + upper + 1
    if size(nw.trace, 1) < width || size(nw.trace, 2) < n + 1
        # TODO: resize the trace matrix if possible
        nw.trace = Matrix{Trace}(width, n + 1)
        nw.H = resize!(nw.H, width)
        nw.E = resize!(nw.E, width)
    end
    nw.lower = lower
    nw.upper = upper
    return nw
end

Base.max(x) = x

# Generate a code block that updates the current score and trace.
# For example, `@update e f g` generates:
#     # update the current score
#     e = E[i-j+U+1]
#     g = H[(i-1)-(j-1)+U+1] + submat[a[i],b_j]
#     h = max(e, f, g)
#     H[i-j1+U+1] = h
#     # update the trace
#     t = trace[i-j+U+1j+1] | ft
#     e == h && (t |= TRACE_DELETE)
#     f == h && (t |= TRACE_INSERT)
#     g == h && (t |= TRACE_MATCH)
#     trace[i-j+U+1,j+1] = t
macro update2(args...)
    @assert :g in args && length(args) ≤ 3

    ex = Expr(:block)

    # update the current score
    if :e in args
        push!(ex.args, :(e = E[i-j+U+1]))
    end
    push!(ex.args, quote
        g = H[(i-1)-(j-1)+U+1] + submat[a[i],b_j]
        h = max($(args...))
        H[i-j+U+1] = h
    end)

    # update the trace
    if :e in args && :f in args
        push!(ex.args, :(t = trace[i-j+U+1,j+1] | ft))
    elseif :e in args
        push!(ex.args, :(t = trace[i-j+U+1,j+1]))
    elseif :f in args
        push!(ex.args, :(t = ft))
    else
        push!(ex.args, :(t = TRACE_NONE))
    end
    if :e in args
        push!(ex.args, :(e == h && (t |= TRACE_DELETE)))
    end
    if :f in args
        push!(ex.args, :(f == h && (t |= TRACE_INSERT)))
    end
    push!(ex.args, quote
        g == h && (t |= TRACE_MATCH)
        trace[i-j+U+1,j+1] = t
    end)

    return esc(ex)
end

macro update2_nextcol(args...)
    if length(args) == 1
        gap_init, = args
        return esc(quote
            e = h + $gap_init
            E[i-(j+1)+U+1] = e
            et = TRACE_NONE
            trace[i-(j+1)+U+1,(j+1)+1] = et
        end)
    elseif length(args) == 2
        gap_init, gap_extend = args
        return esc(quote
            e′ = e + $gap_extend
            e = max(e′, h + $gap_init)
            E[i-(j+1)+U+1] = e
            et = TRACE_NONE
            e == e′ && (et |= TRACE_EXTDEL)
            trace[i-(j+1)+U+1,(j+1)+1] = et
        end)
    end
    @assert false
end

macro update2_nextrow(gap_init, gap_extend)
    esc(quote
        f′ = f + $gap_extend
        f = max(f′, h + $gap_init)
        ft = TRACE_NONE
        f == f′ && (ft |= TRACE_EXTINS)
    end)
end

islowerbound(i, j, L) = j == i - L

function run!{T}(
        nw::BandedNeedlemanWunsch{T},
        a, b,
        submat::AbstractSubstitutionMatrix{T},
        # a
        start_gap_open_a::T,
        start_gap_extend_a::T,
        middle_gap_open_a::T,
        middle_gap_extend_a::T,
        end_gap_open_a::T,
        end_gap_extend_a::T,
        # b
        start_gap_open_b::T,
        start_gap_extend_b::T,
        middle_gap_open_b::T,
        middle_gap_extend_b::T,
        end_gap_open_b::T,
        end_gap_extend_b::T)

    m = length(a)
    n = length(b)
    L = nw.lower
    U = nw.upper
    ensureroom!(nw, m, n, L, U)
    H = nw.H
    E = nw.E
    trace = nw.trace

    start_gap_init_a = start_gap_open_a + start_gap_extend_a
    middle_gap_init_a = middle_gap_open_a + middle_gap_extend_a
    end_gap_init_a = end_gap_open_a + end_gap_extend_a

    start_gap_init_b = start_gap_open_b + start_gap_extend_b
    middle_gap_init_b = middle_gap_open_b + middle_gap_extend_b
    end_gap_init_b = end_gap_open_b + end_gap_extend_b

    # In order to save the working space, the matrices are vertically sheared.
    # Namely, the coordinate is transformed as (i, j) → (i-j+U, j), where
    # (i, j) points to the subproblem of the optimal alignment between two
    # prefixes, a[1:i] and b[1:j], within the band.
    # The position (i, j) is accessed with:
    #   * trace[i-j+U+1,j+1]
    #   * H[i-j+U+1]
    #   * E[i-j+U+1]  (j ≥ 1)

    H[0-0+U+1] = T(0)
    trace[0-0+U+1,0+1] = TRACE_NONE
    if m == 0 && n == 0
        return H[0-0+U+1]
    elseif m == 0
        for j in 1:U
            trace[0-j+U+1,j+1] = TRACE_DELETE
        end
        return start_gap_open_a + start_gap_extend_a * n
    elseif n == 0
        for i in 1:L
            trace[i-0+U+1,0+1] = TRACE_INSERT
        end
        return start_gap_open_b + start_gap_extend_b * m
    else
        # initialize the first column of H, E, and trace
        @inbounds for i in 1:L
            H[i-0+U+1] = start_gap_open_b + start_gap_extend_b * i
            E[i-1+U+1] = H[i-0+U+1] + (i == m ? end_gap_init_a : middle_gap_init_a)
            trace[i-0+U+1,0+1] = TRACE_INSERT
            trace[i-1+U+1,1+1] = TRACE_NONE
        end

        # run dynamic programming column by column (except the last column)
        @inbounds for j in 1:n-1
            b_j = b[j]
            # vertical bounds along the j-th column
            lo = max(0, j - U)
            hi = min(m, j + L)

            # fill the first cell of the column in the band
            f, ft = let i = lo
                if lo == 0
                    h = H[i-j+U+1] = start_gap_open_a + start_gap_extend_a * j
                    trace[i-j+U+1,j+1] = TRACE_DELETE
                elseif islowerbound(i, j, L)
                    @update2 g
                else
                    @update2 e g
                end
                h + middle_gap_init_b, TRACE_NONE
            end

            for i in lo+1:hi-1
                @update2 e f g
                @update2_nextcol middle_gap_init_a middle_gap_extend_a
                @update2_nextrow middle_gap_init_b middle_gap_extend_b
            end

            # fill the last cell of the column
            if hi > lo
                let i = hi
                    et = TRACE_NONE
                    if islowerbound(i, j, L)
                        @update2 f g
                        @update2_nextcol (i == m ? end_gap_init_a : middle_gap_init_a)
                    else
                        @update2 e f g
                        @update2_nextcol (i == m ? end_gap_init_a : middle_gap_init_a) (i == m ? end_gap_extend_a : middle_gap_extend_a)
                    end
                end
            end
        end

        # fill the last column
        let j = n
            b_j = b[j]
            # vertical bounds along the j-th column
            lo = max(0, j - U)
            hi = min(m, j + L)

            # fill the first cell of the column
            f, ft = let i = lo
                if lo == 0
                    h = H[i-j+U+1] = start_gap_open_a + start_gap_extend_a * j
                    trace[i-j+U+1,j+1] = TRACE_DELETE
                elseif islowerbound(i, j, L)
                    @update2 g
                else
                    @update2 e g
                end
                h + end_gap_init_b, TRACE_NONE
            end

            @inbounds for i in lo+1:hi-1
                @update2 e f g
                @update2_nextcol end_gap_init_b end_gap_extend_b
            end

            if hi > lo
                let i = hi
                    if islowerbound(i, j, L)
                        @update2 f g
                    else
                        @update2 e f g
                    end
                end
            end
        end

        return H[m-n+U+1]
    end
end

function run!(nw::BandedNeedlemanWunsch, a, b, submat, gap_open, gap_extend)
    return run!(nw, a, b, submat,
        # a
        gap_open,
        gap_extend,
        gap_open,
        gap_extend,
        gap_open,
        gap_extend,
        # b
        gap_open,
        gap_extend,
        gap_open,
        gap_extend,
        gap_open,
        gap_extend
    )
end

function traceback(nw::BandedNeedlemanWunsch, a, b, endpos)
    U = nw.upper
    anchors = Vector{AlignmentAnchor}()
    i, j = endpos
    was_extins = false
    was_extdel = false
    @start_traceback
    while i ≥ 1 || j ≥ 1
        t = nw.trace[i-j+U+1,j+1]
        @assert !(was_extins && was_extdel)
        if was_extins
            was_extins = t & TRACE_EXTINS > 0
            @anchor OP_INSERT
        elseif was_extdel
            was_extdel = t & TRACE_EXTDEL > 0
            @anchor OP_DELETE
        elseif t & TRACE_INSERT > 0
            was_extins = t & TRACE_EXTINS > 0
            @anchor OP_INSERT
        elseif t & TRACE_DELETE > 0
            was_extdel = t & TRACE_EXTDEL > 0
            @anchor OP_DELETE
        elseif t & TRACE_MATCH > 0
            if a[i] == b[j]
                @anchor OP_SEQ_MATCH
            else
                @anchor OP_SEQ_MISMATCH
            end
        else
            error("failed to trace back")
        end
    end
    @finish_traceback
    return AlignedSequence(a, anchors)
end


text:
__precompile__()

module EnhancedGJK

using StaticArrays
import GeometryTypes
const gt = GeometryTypes
using CoordinateTransformations: Transformation,
                                 transform_deriv,
                                 IdentityTransformation
using LinearAlgebra
using Statistics: mean

export CollisionCache,
       gjk!,
       gjk,
       GJKResult,
       NeighborMesh,
       ReferenceDistance,
       closest_point_in_world,
       closest_point_in_body,
       separation_distance,
       simplex_penetration_distance

include("tagged_points.jl")
include("simplices.jl")
include("johnson_distance.jl")
include("neighbor_mesh.jl")
include("traits.jl")
include("gjk.jl")
include("reference_distance.jl")

end # module


text:
function apply_noise!(H::SparseMatrixCSC, Ad::SparseMatrixCSC, noisearray)

# counter telling which of the indipendent links we are considering
    counter = 1
    for i = 1 : Int(sqrt(length(H)))

        # this will be the value of the diagonal element [i,i],
        # which depends on the noise on each
        # link starting from the considered node i.
        diagNoisyValue = 0.

        # We perturb the laplacian matrix of the graph with the noise
        # We do this at low level by acting on the sparse matrix
        # representation

        # elements of the row vector indicating which rows
        # are occupied in column i
        for j = Ad.colptr[i] : (Ad.colptr[i + 1] - 1)

            # We insert noise only in the lower triangular part
            # of the Hamiltonian matrix
            if i < Ad.rowval[j]

                H[Ad.rowval[j], i] = noisearray[counter]

                 counter+=1

            elseif i > Ad.rowval[j]
                H[Ad.rowval[j], i] = H[i, Ad.rowval[j]]
            end

            diagNoisyValue += -1. * H[Ad.rowval[j],i]
        end

        # the value of the diagonal element of a column (node i)
        # must be = minus the sum of all the other values in the
        # column, because of the conservation of the overall
        # probability sum_j p(j|i)=1
        H[i,i] = diagNoisyValue
    end
end

max_stars_repo_name: matteoacrossi/QuantumSpatialSearch

max_issues_repo_name: matteoacrossi/QuantumSpatialSearch

max_forks_repo_name: matteoacrossi/QuantumSpatialSearch


text:
# Autogenerated wrapper script for cilkrts_jll for x86_64-apple-darwin
export libcilkrts

JLLWrappers.@generate_wrapper_header("cilkrts")
JLLWrappers.@declare_library_product(libcilkrts, "@rpath/libcilkrts.5.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcilkrts,
        "lib/libcilkrts.5.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
## Packages
using Roots, Calculus, Plots

## Function for finding root
function findZero(g, init_guess)

    xn1 = init_guess;
    tol = 10^(-8);
    distance = 1;
    iter = 1;

    while distance >= tol
        xn = xn1;
        xn1 = xn - g(xn) / derivative(g, xn);
        xn1 = max(xn1, 0.0001);

        distance = abs(xn1 - xn);
        iter = iter + 1;
    end

    return xn1
end

## Function for defining utility function
function u(c)
    if c <= 0
        u = -Inf;
    else
        u = c^(1/2);
    end
    return u
end

## Parameters
β = 0.8;
πs1 = 1/2;
πs2 = 1/2;
Ys1 = 1;
Ys2 = 1/2;
M = 20; # Number of grids

## Grid for v
vmax = (πs1 * u(Ys1) + πs2 * u(Ys2)) / (1 - β);
vgrid = range(0.0, stop = vmax, length = M);

## Iterate on Value function
Pi1 = zeros(M);
cs1_grid = range(0, stop = Ys1, length = M);
ws1_grid = vgrid;
ws2_grid = vgrid;

distance = 1;
tol = 10^(-8);
iter = 1;

policy_fun = Array{CartesianIndex, 1}(undef, M);
policy_cs1 = Array{Int64, 1}(undef, M);

while distance >= tol
    println("This is iteration $iter")
    Pi = copy(Pi1);

    for vIndex in 1:M
        U = Array{Float64, 3}(undef, M, M, M)
        v = vgrid[vIndex];
        for cs1Index in 1:M, ws1Index in 1:M, ws2Index in 1:M
            cs1 = cs1_grid[cs1Index];
            ws1 = ws1_grid[ws1Index];
            ws2 = ws2_grid[ws2Index];

            # Solve for cs2
            f(cs2) = v - (πs1 * (cs1^(1/2) + β * ws1)) - (πs2 * (cs2^(1/2) + β * ws2));
            cs2 = findZero(f, 1.0);

            U_temp = (πs1 * (u(Ys1 - cs1) + β * Pi[ws1Index])) +
                     (πs2 * (u(Ys2 - cs2) + β * Pi[ws2Index]));

            if U_temp <= 0.0
                U[cs1Index, ws1Index, ws2Index] = 0.0;
            else
                U[cs1Index, ws1Index, ws2Index] = U_temp;
            end
        end
        Pi1[vIndex], policy_fun[vIndex] = findmax(U[:,:,:]);
    end

    distance = maximum(abs.(Pi1 - Pi));
    iter = iter + 1;
end

## A Check
λ = range(0.0, stop = 0.999, length = 100);
u1 = zeros(100);
u2 = zeros(100);
for i in 1:100
    A = (λ[i] / (1 - λ[i]))^2;
    cs1 = Ys1 * (A / (1 + A));
    cs2 = Ys2 * (A / (1 + A));
    u1[i] = (πs1 * u(cs1) + πs2 * u(cs2)) / (1 - β);
    u2[i] = (πs1 * u(Ys1 - cs1) + πs2 * u(Ys2 - cs2)) / (1 - β);
end

plot(vgrid, Pi1, xlabel = "v", ylabel = "P", label = "VFI", legend =:bottomleft)
plot!(u1, u2, label = "Theoretical")


text:
# Test for weights with a small set of trials

using SubmodularMaximization
using PyPlot
using Base.Iterators
using Base.Threads
using JLD2
using Statistics
using HistogramFilters
using RosDataProcess

close("all")

data_folder = "./data"

reprocess = false

# Configuration
targets_per_robot = 1
variance_scaling_factor = 0.5
range_limit = Inf
grid_cells_per_robot = 10^2 / 8

experiment_name = string("minimal_weights_by_num_robots",
                         "_", targets_per_robot,
                         "_", variance_scaling_factor,
                         "_", range_limit,
                         "_", grid_cells_per_robot)

num_robots = 4:4:40

steps = 40
trial_steps = 20:steps

# Note: we will run trials in threads so the solvers do not have to be threaded
solver = x->solve_n_partitions(4, x, threaded=true)

#
# Code to run experiments
#

all_tests = num_robots

function trial_fun(num_robots)
  grid_size = round(Int64, sqrt(grid_cells_per_robot * num_robots))
  num_targets = round(Int64, targets_per_robot * num_robots)

  sensor = RangingSensor(variance_scaling_factor=variance_scaling_factor,
                         range_limit=range_limit)
  grid = Grid(grid_size, grid_size)

  configs = MultiRobotTargetTrackingConfigs(grid=grid, sensor=sensor)

  trial_data = target_tracking_experiment(steps=steps,
                                          num_robots=num_robots,
                                          num_targets=num_targets,
                                          configs=configs,
                                          solver=solver)

  weights = map(trial_data[trial_steps]) do step_data
        robot_states = step_data.robot_states
        target_filters = step_data.target_filters

        problem = MultiRobotTargetTrackingProblem(robot_states, target_filters,
                                                  configs)

        compute_weight_matrix(problem, threaded=true)
  end

  (trial_weights=weights, trial_data=trial_data, configs=configs)
end
function print_summary(num_robots)
  println("Num. robots: ", num_robots)
end

@time data = run_experiments(all_tests,
                             trial_fun=trial_fun,
                             print_summary=print_summary,
                             experiment_name=experiment_name,
                             data_folder=data_folder,
                             reprocess=reprocess,
                             threaded=false
                            )

normalized_weights = map(num_robots) do num_robots
  trial = data[num_robots]

  total_weights = map(total_weight, trial.trial_weights)
  objectives = map(x->x.objective, trial.trial_data[trial_steps])

  total_weights ./ objectives
end

weights_series = TimeSeries(num_robots, hcat(normalized_weights...)')
plot_trials(weights_series, mean=true, marker=".", markersize=25, linewidth=3)

title_string = string("Tgts. per rbt.:", targets_per_robot,
                      " Var. scale:", variance_scaling_factor,
                      " Range lim.:", range_limit,
                      " Cells per rbt.:", grid_cells_per_robot)
title(title_string)
ylabel("Cost bound over obj.")
xlabel("Number of robots")

save_fig("fig", experiment_name)

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
#=

In this file we simulate a simple rocket landing

Units in kg, m, s

=#

# RUN:
# ] activate TrajOptSOCPs

using TrajOptSOCPs
using LinearAlgebra, SparseArrays

# If you don't want plotting comment the include and pyplot statements.
# Then, sent runplots = false
include("../plotting/plotTrajectory.jl")
include("../plotting/plotConstraintViolation.jl")
include("../plotting/plotObjective.jl")
include("../plotting/batchPlots.jl")
pyplot()
runplots = true
saveplots = false

println("\n--------------------------------------------")
println("          Setting Up Rocket Landing ")
println("--------------------------------------------")

# Based on the Falcon 9
# 549,054 kg (Mass)
# 282 s (Specific Impulse)
# select "y" as the vertical
mass = 1
isp = 1
grav = [0; -9.81]
deltaTime = 1
rocket = rocket_simple(mass, isp, grav, deltaTime)

# in m
# The Karman Line (100 km)
const rocketStart = [2.0; 20.0; 0.0; -18.0]
const rocketEnd = [0.0; 0.0; 0.0; 0.0]#[-5.0; 0.0; 0.0; 0.0]

uHover = mass * grav

# Number of time steps to discretize the trajectory
const NSteps = 60
# Initialize the trajectory with a line
initTraj = initializeTraj(rocketStart, rocketEnd, uHover, uHover, NSteps)

# Use a Linear Quadratic Regulator as the cost function
const lqrQMat = 0.0001 * Diagonal(I, size(rocketStart, 1))
const lqrRMat = 0.0025 * Diagonal(I, Int64(size(rocketStart, 1) / 2))
costFun = makeLQR_TrajReferenced(lqrQMat, lqrRMat, NSteps, initTraj)

# Create the Dynamics Constraints
const ADyn, BDyn = rocketDynamicsFull(rocket, rocketStart, rocketEnd, NSteps)
dynConstraint = AL_AffineEquality(ADyn, BDyn)
lambdaInit = -1 * ones(size(BDyn))

# Create the Ground Constraints
groundConstraint = makeGroundConstraint(NSteps, size(grav, 1), size(grav, 1))
groundLambda = zeros(size(groundConstraint.b, 1))

# Create the Max Thrust Constraints
thrustMax = 20.0
maxThrustConstraint = makeMaxThrustConstraint(NSteps, size(grav, 1), thrustMax)
maxThrustLambda = zeros(size(maxThrustConstraint.indicatorList))

# Create the constraint manager

# cMRocket = constraintManager_Base([dynConstraint], [lambdaInit])
# cMRocket = constraintManager_Dynamics([], [], dynConstraint, lambdaInit)
# cMRocket = constraintManager_Dynamics([groundConstraint], [groundLambda],
#                                         dynConstraint, lambdaInit)
cMRocket = constraintManager_Dynamics([groundConstraint, maxThrustConstraint],
                                      [groundLambda, maxThrustLambda],
                                      dynConstraint, lambdaInit)

# Initialize the primal-dual vector
initTrajPD = [initTraj; lambdaInit]

# Initialize the Augmented Lagrangian penalty on the constraints
penaltyStart = 1.0

# Test that the evaluations work
println("\n--------------------------------------------")
println(" Testing constraints are inputted correctly ")
println("--------------------------------------------")
println("Starting constraint violation: ")
println([evalConstraints(cMRocket, initTraj, penaltyStart)])
println("Starting gradient of constraints: ")
println(size(evalGradConstraints(cMRocket, initTraj, penaltyStart)))
println("Starting hessian of constraints: ")
println(size(evalHessConstraints(cMRocket, initTraj)))

# Equiped with the constraint term and the objective term, I now build the
# Augmented Lagrangian
alRocket = augLag(costFun, cMRocket, penaltyStart)

# Test that the evaluations work
println("--------------------------------------------")
println("           Testing AL is Functional         ")
println("--------------------------------------------")
println("Evaluating augmented lagrangian: ")
println([evalAL(alRocket, initTrajPD)])
println("Evaluating gradient of augmented lagrangian: ")
println(size(evalGradAL(alRocket, initTrajPD)))
println("Evaluating hessian of augmented lagrangian: ")
println(size(evalHessAl(alRocket, initTrajPD)))

# Next we select resonable solver parameters
currSolveParams = solverParams(0.1, 0.5,
                                8, 10,
                                10^-4,
                                10, 10^6,
                                2.5, 2, 0.2, 0.2, 0.4)
solParamPrint(currSolveParams)
println()
println()

println("--------------------------------------------")
println("             Beginning Solve                ")
println("--------------------------------------------")

# Solve the Trajectory Optimization problem
useALMethod = true

if useALMethod
    # Use an augmented lagrangian method
    trajLambdaSolved, resArr = ALPrimalNewtonMain(initTrajPD, alRocket,
                                            currSolveParams)
    trajStatesAllPD = [parsePrimalDualVec(trajL, size(initTraj, 1))
                                                for trajL in trajLambdaSolved]
    trajStatesAll = [trajPDThis.primals for trajPDThis in trajStatesAllPD]
else
    # Solve only the objective and linear dynamics
    costQ = costFun.lqr_qp.QR_Full
    costP = - (initTraj' * costQ)'
    trajLambdaSolvedTrue = solveQP_AffineEq(costQ, costP, ADyn, BDyn)
    trajStateLastTrue = parsePrimalDualVec(trajLambdaSolvedTrue,
                                                        size(initTraj, 1))
end

# Blocked so that it can be run independently after the fact
if runplots

    hDual = heatmap(hcat(trajStatesAllPD[end].duals))
    yflip!()
    ylabel!("Column")
    title!("Dual Vector for Dynamics Constraints")
    display(hDual)

    # Get the parsed list of trajectories
    nDim = size(grav, 1)
    pltTraj, pltCV, pltCV2, pltObj, plts, pltv, pltu =
                                    batchPlot(trajStatesAll, cMRocket, nDim)
end

# Blocked so that it can be run independently after the fact
if runplots && saveplots
    header = "MaxThrustNew_" * string(currSolveParams.maxOuterIters) *
             "Outer_" * string(currSolveParams.maxNewtonSteps) * "Newton" *
             string(Int64(rocketStart[4])) * "Vel" * "_"
    saveBulk(pltTraj, pltCV, pltCV2, pltObj, plts, pltv, pltu, header)
end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
using Mocha

data_layer  = HDF5DataLayer(name="train-data", source="data/train.txt",
    batch_size=64, shuffle=true)

conv_layer = ConvolutionLayer(name="conv1", n_filter=20, kernel=(5,5),
    bottoms=[:data], tops=[:conv1])

pool_layer = PoolingLayer(name="pool1", kernel=(2,2), stride=(2,2),
    bottoms=[:conv1], tops=[:pool1])
conv2_layer = ConvolutionLayer(name="conv2", n_filter=50, kernel=(5,5),
    bottoms=[:pool1], tops=[:conv2])
pool2_layer = PoolingLayer(name="pool2", kernel=(2,2), stride=(2,2),
    bottoms=[:conv2], tops=[:pool2])

fc1_layer  = InnerProductLayer(name="ip1", output_dim=500,
    neuron=Neurons.ReLU(), bottoms=[:pool2], tops=[:ip1])
fc2_layer  = InnerProductLayer(name="ip2", output_dim=10,
    bottoms=[:ip1], tops=[:ip2])

loss_layer = SoftmaxLossLayer(name="loss", bottoms=[:ip2,:label])

backend = CPUBackend()
init(backend)

common_layers = [conv_layer, pool_layer, conv2_layer, pool2_layer,
    fc1_layer, fc2_layer]
net = Net("MNIST-train", backend, [data_layer, common_layers..., loss_layer])

exp_dir = "snapshots"
method = SGD()
params = make_solver_parameters(method, max_iter=10000, regu_coef=0.0005,
    mom_policy=MomPolicy.Fixed(0.9),
    lr_policy=LRPolicy.Inv(0.01, 0.0001, 0.75),
    load_from=exp_dir)
solver = Solver(method, params)

setup_coffee_lounge(solver, save_into="$exp_dir/statistics.hdf5", every_n_iter=1000)

add_coffee_break(solver, TrainingSummary(), every_n_iter=100)
add_coffee_break(solver, Snapshot(exp_dir), every_n_iter=5000)

data_layer_test = HDF5DataLayer(name="test-data", source="data/test.txt", batch_size=100)
acc_layer = AccuracyLayer(name="test-accuracy", bottoms=[:ip2, :label])
test_net = Net("MNIST-test", backend, [data_layer_test, common_layers..., acc_layer])

add_coffee_break(solver, ValidationPerformance(test_net), every_n_iter=1000)

solve(solver, net)

destroy(net)
destroy(test_net)
shutdown(backend)

max_stars_repo_path: Chapter 11/ch11-Deeplearning.jl

max_stars_repo_name: PacktPublishing/Julia-for-data-science

max_issues_repo_path: Chapter 11/ch11-Deeplearning.jl

max_issues_repo_name: PacktPublishing/Julia-for-data-science

max_forks_repo_path: Chapter 11/ch11-Deeplearning.jl

max_forks_repo_name: PacktPublishing/Julia-for-data-science


text:
"""
digivec = digitizing_3elesvec(vec, max, steps, minus)

3要素ベクトルを離散化する．かつ，最大値を上回る場合は最大値に抑える

# Argument
 - `vec` : 離散化対象の3要素ベクトル
 - `max` : ベクトルの要素最大値
 - `steps` : 離散化段階（正負方向への離散化であれば，正方向に何段階分割するか）
 - `minus` : 要素に負値を許容するか，許容するならtrue，しないなら非記入
# return
 - `digivec` : 離散化後のベクトル
"""

function digitizing_3elesvec(vec, max, steps, minus = false)
    resol = max / steps
    digivec = round.(vec / resol) * resol
    for i=1:3
        if norm(digivec[i]) > max
            digivec[i] = digivec[i] * max / norm(digivec[i])
        end

        if digivec[i] < 0. && minus == false
            digivec[i] = 0.
        end
    end

    return digivec
end

"""
Treq, M, n = cp_and_bdot(tarqua, qua_ad, omega_ad, kp, kr, magvec_scsf, boundangle)

クロスプロダクト則による所望トルクと地磁場のなす角が 90 ± n deg の時はクロスプロダクト則，他の時はBdot則による所望トルクに切り替えて制御を行う．

# Argument
 - `tarqua` : 目標姿勢クォータニオン
 - `qua_ad` : 現在姿勢クォータニオン
 - `omega_ad` : 現在の回転角速度
 - `kp` : クロスプロダクト則の比例ゲイン
 - `kr` : クロスプロダクト則の微分ゲイン
 - `magvec_scsf` : 地磁場ベクトル＠SCSF座標系
 - `boundangle` : クロスプロダクト則とB-dot則を切り替える角度（上記説明のnにあたる）
# Return
 - `Treq` : 所望トルク
 - `M` : 所望磁気モーメント
 - `n` : 制御方式（n=0でクロスプロダクト，n=1でBdot）

"""
function cp_and_bdot(tarqua, qua_ad, omega_ad, kp, kr, magvec_scsf, boundangle)
    n = 0
    Treq, M = cross_product(tarqua,qua_ad, kp, kr, omega_ad, magvec_scsf)
    φ = dot(Treq, magvec_scsf) / norm(Treq) / norm(magvec_scsf)
    if norm(φ) > cos(deg2rad(boundangle))
        if i==1
            M = B_dot(magvec_scsf, sat_ω[i, :], sat_ω[i, :])
        else
            M = B_dot(magvec_scsf, sat_ω[i, :], sat_ω[i-1, :])
        end

        Treq = cross(M, magvec_scsf)
        n = 1

    end

    return Treq, M, n
end

"""
vec1, θ = anglediffs(qua, targetqua, cam_origindir)

"""
function anglediffs(qua, targetqua, cam_origindir)
    qua1 = qua * cam_origindir / qua
    vec1 = [qua1.q1, qua1.q2, qua1.q3]
    qua2 = targetqua * cam_origindir / targetqua
    vec2 = [qua2.q1, qua2.q2, qua2.q3]
    cθ = dot(vec1, vec2) / norm(vec1) / norm(vec2)
    if cθ > 1.
        cθ = 1
    end
    θ = rad2deg(acos(cθ))

    return vec1, θ
end

"""
tarpos = cal_targetlocus(targetpos_ecef, x_ecef_log, x_geod_log, v_ecef_log, x_ecef_error, qua, i)
"""
function cal_targetlocus(targetpos_ecef, x_ecef_log, x_geod_log, v_ecef_log, x_ecef_error, qua, i)
    sat2tar_ecef = targetpos_ecef - (x_ecef_log[i, :] + x_ecef_error)
    sat2tar_seof = ecef_to_DCM(x_ecef_log[i, :] + x_ecef_error, v_ecef_log[i, :], true) * sat2tar_ecef
    sat2tar_seof = sat2tar_seof / norm(sat2tar_seof) # 衛星→撮影対象の単位方向ベクトル
    sat2tar_scsfqua = qua \ sat2tar_seof * qua
    sat2tar_scsf = [sat2tar_scsfqua.q1, sat2tar_scsfqua.q2, sat2tar_scsfqua.q3]
    sat2tar_scsf = sat2tar_scsf / norm(sat2tar_scsf)                  # カメラの単位方向ベクトル
    tarpos = [sat2tar_scsf[1] * x_geod_log[i, 3] / sat2tar_scsf[3], sat2tar_scsf[2] * x_geod_log[i, 3] / sat2tar_scsf[3]]

    return tarpos
end

max_stars_repo_name: yui-project/ACSim

max_issues_repo_name: yui-project/ACSim

max_forks_repo_name: yui-project/ACSim


text:
 #map(abspath,readdir()) # get the file location and file names
 # create log file
using DelimitedFiles
writedlm("log_output_sigma1.txt",  " ")
writedlm("log_output_sigma10.txt", " ")
include("GB6.jl")
import GLM: fit, predict, LinearModel
using Distributions,SpecialFunctions
using LinearAlgebra, MultivariateStats
using StatsBase, Random, Distributed
using DelimitedFiles, DataStructures
using CSV
using DataFrames

Random.seed!(2019)
rng = 2019

# 1) Simulate n (= 500 )× p ( = 1~1000) training data and test data, where σ² = 1 and 10 separately
# 2) Plot the broken line chart

gbdt = GB6.GBDT(;
  sampling_rate = 0.75,                # Sampling rate
  learning_rate = 0.1,                # Learning rate
  num_iterations = 100,               # NuGBer of iterations
  tree_options = Dict(
      :nsubfeatures => 0,
      :min_samples_leaf => 1,
      :max_depth => 5))
gbl = GB6.GBLearner(gbdt, :regression)

######## σ² = 1 #########
temp_storedim = size(collect(10:1000),1)
fitval = Array{Union{Missing, Nothing, Int, Float64}}(nothing, 500, temp_storedim)
rssfitted = []
predicted = Array{Union{Missing, Nothing, Int, Float64}}(nothing, 500, temp_storedim)
rmsefitted = []
Ytestval = Array{Union{Missing, Nothing, Int, Float64}}(nothing, 500, temp_storedim)
μ_testval = Array{Union{Missing, Nothing, Int, Float64}}(nothing, 500, temp_storedim)
rmsefitted2 = []
s
for pval in collect(10:1000)
    # reproducibility
    Random.seed!(2019)
    rng = 2019

    ### simulate the training and testing data
    X, Y = GB6.test_data(500; p = pval, sigma = 1.0)
    X_test, Y_test = GB6.test_data(500; p = pval)
    μ = GB6.fried(X)
    μ_test = GB6.fried(X_test)

    println("The last simulated data:
        \nX:$(X[500, pval]) \nY:$(Y[500])
        \nX_test: $(X_test[500, pval]) \nY_test: $(Y_test[500])
        \nμ: $(μ[500]) \nμ_test: $(μ_test[500])\n\n")

    μ_testval[:, pval - 9] = μ_test
    Ytestval[:, pval - 9 ] = Y_test

    # train the model
    mymodel=GB6.fit!(gbl, X, Y)
    # get the fitted value
    myfit  = GB6.predict!(gbl, X)
    # store the fitted value
    fitval[:, pval - 9] = myfit

    # get the residual sum of squares
    rssfit = sqrt(mean(Y .- myfit) .^2)
    push!(rssfitted, rssfit)
    # get prediction
    mypred = GB6.predict!(gbl, X_test)

    println("The last pred is $(mypred[500])\n\n ")

    predicted[:, pval - 9] = mypred

    # get mse and root of mean squared error
    rmse = sqrt(mean((Y_test .- mypred ) .^2))

    println("The rmse b.w. Y_test and pred for p = $(pval) is $(rmse) ")

    push!(rmsefitted, rmse)

    # get the rmse(μ_test, mypred)
    rmse2 = sqrt(mean((μ_test .- mypred).^2))

    println("The rmse b.w. μ_test and pred for p = $(pval)  is $(rmse2).")

    push!(rmsefitted2, rmse2)
end

rssfitted
rmsefitted
rmsefitted2
Ytestval
μ_testval

fitval = DataFrame(fitval)
CSV.write("fitval_sigma1.csv", fitval)
rssfitted = DataFrame(rssfitted = rssfitted)
CSV.write("fss_sigma1.csv", rssfitted)
rmsefitted = DataFrame(rmsefitted = rmsefitted)
CSV.write("rmse_sigma1.csv", rmsefitted)
predicted = DataFrame(predicted)
CSV.write("pred_sigma1.csv", predicted)
rmsefitted2 = DataFrame(rmsefitted2 = rmsefitted2)
CSV.write("rmse2_sigma1.csv", rmsefitted2)
Ytestval = DataFrame(Ytestval)
CSV.write("Ytest_sigma1.csv", Ytestval)
μ_testval = DataFrame(μ_testval)
CSV.write("mu_test_sigma1.csv", μ_testval)
##### not runned yet

######## σ² = 10 #########

fitval_sigma10 = Array{Union{Missing, Nothing, Int, Float64}}(nothing, 500, temp_storedim)
rssfitted_sigma10 = []
predicted_sigma10 = Array{Union{Missing, Nothing, Int, Float64}}(nothing, 500, temp_storedim)
rmsefitted_sigma10 = []
Ytestval_sigma10 = Array{Union{Missing, Nothing, Int, Float64}}(nothing, 500, temp_storedim)
μ_testval_sigma10 = Array{Union{Missing, Nothing, Int, Float64}}(nothing, 500, temp_storedim)
rmsefitted2_sigma10 = []
s
for pval in collect(10:1000)
    # reproducibility
    Random.seed!(2019)
    rng = 2019

    ### simulate the training and testing data
    X, Y = GB6.test_data(500; p = pval, sigma = sqrt(10.0))
    X_test, Y_test = GB6.test_data(500; p = pval, sigma = sqrt(10.0))
    μ = GB6.fried(X)
    μ_test = GB6.fried(X_test)

    println("Simulated data:
    \nX: $(X[500, pval])
    \nY: $(Y[500])
    \nX_test: $(X_test[500, pval])
    \nY_test: $(Y_test[500])
    \nμ: $(μ[500])
    \nμ_test: $(μ_test[500])")

    μ_testval_sigma10[:, pval - 9] = μ_test
    Ytestval_sigma10[:, pval - 9] = Y_test

    # train the model
    mymodel=GB6.fit!(gbl, X, Y)
    # get the fitted value
    myfit  = GB6.predict!(gbl, X)
    # store the fitted value
    fitval_sigma10[:, pval - 9] = myfit

    # get the residual sum of squares
    rssfit = sqrt(mean((Y .- myfit) .^2))
    push!(rssfitted_sigma10, rssfit)
    # get prediction
    mypred = GB6.predict!(gbl, X_test)

    println("The last pred is $(mypred[500])")

    predicted_sigma10[:, pval - 9] = mypred

    # get mse and root of mean squared error
    rmse = sqrt(mean((Y_test .- mypred ) .^2))
    println("The p = $(pval), and the rmse is $rmse.")
    push!(rmsefitted_sigma10, rmse)

    # get the rmse(Y_test, μ_test)
    rmse2 = sqrt(mean((mypred .- μ_test) .^2))
    println("The p = $(pval) dim, rmse2 is $rmse2.")
    push!(rmsefitted2_sigma10, rmse2)
end

fitval_sigma10
rssfitted_sigma10
predicted_sigma10
rmsefitted_sigma10
Ytestval_sigma10
μ_testval_sigma10
rmsefitted2_sigma10

#predicted_sigma10 = DataFrame(predicted_sigma10)
rmsefitted_sigma10 = DataFrame(rmsefitted_sigma10=rmsefitted_sigma10)

#CSV.write("fitval_sigma10.csv", fitval_sigma10)
#CSV.write("rssfitted_sigma10.csv", rssfitted_sigma10)
#CSV.write("predicted_sigma10.csv", predicted_sigma10)
#CSV.write("rmsefitted_sigma10.csv", rmsefitted_sigma10)

fitval_sigma10 = DataFrame(fitval_sigma10)
CSV.write("fitval_sigma10.csv", fitval)
rssfitted_sigma10 = DataFrame(rssfitted_sigma10 = rssfitted_sigma10)
CSV.write("fss_sigma10.csv", rssfitted_sigma10)

CSV.write("rmse_sigma10.csv", rmsefitted_sigma10)
predicted_sigma10 = DataFrame(predicted_sigma10)
CSV.write("pred_sigma10.csv", predicted_sigma10)
rmsefitted2_sigma10 = DataFrame(rmsefitted2_sigma10 = rmsefitted2_sigma10)
CSV.write("rmse2_sigma10.csv", rmsefitted2_sigma10)
Ytestval_sigma10 = DataFrame(Ytestval_sigma10)
CSV.write("Ytest_sigma10.csv", Ytestval_sigma10)
μ_testval_sigma10 = DataFrame(μ_testval_sigma10)
CSV.write("mu_test_sigma10.csv", μ_testval_sigma10)

# plot the rmsefitted for sigma = 1
using Gadfly

rmse_sigma1 = CSV.read("rmse_sigma1.csv")
rmse_sigma10 = CSV.read("rmse_sigma10.csv")

rmse2_sigma1 = CSV.read("rmse2_sigma1.csv")
rmse2_sigma10 = CSV.read("rmse2_sigma10.csv")

 # 991
i = collect(1:1:10)
ii = collect(20:10:100)
iii = collect(200:100:900)
ind = vcat(i, ii, iii)
rmse1_forplot = rmse_sigma1[ind,1]
rmse10_forplot = rmse_sigma10[ind,1]

p_1 = layer(x = ind, y = rmse1_forplot,
            Geom.point, Geom.line,
            Theme(default_color="red"))
p_2 = layer(x = ind, y = rmse10_forplot,
            Geom.point, Geom.line, Theme(default_color="blue"))
myplot1 = plot(p_1, p_2, Guide.xlabel("P"),
    Guide.ylabel("rmse"), Guide.title("Root Mean Squared Error"),
    Guide.manual_color_key("", ["σ² = 1", "σ² = 10"], ["red", "blue"]))

using Cairo, Fontconfig
draw(PNG("myplot1.png", 3inch, 3inch), myplot1)

rmse2_1_forplot = rmse2_sigma1[ind,1]
rmse2_10_forplot = rmse2_sigma10[ind,1]

pp_1 = layer(x = ind, y = rmse2_1_forplot,
            Geom.point, Geom.line,
            Theme(default_color="red"))
pp_2 = layer(x = ind, y = rmse2_10_forplot,
            Geom.point, Geom.line, Theme(default_color="blue"))
myplot2 = plot(pp_1, pp_2, Guide.xlabel("P"),
    Guide.ylabel("rmse"), Guide.title("Root Mean Squared Error"),
    Guide.manual_color_key("", ["σ² = 1", "σ² = 10"], ["red", "blue"]));

Gadfly.push_theme(style(line_width=1mm))

ppp_2 = layer(x = ind, y = rmse2_10_forplot,
            style(line_width=10cm),
            Geom.point, Geom.line, Theme(default_color="blue"))
pppp = plot(ppp_2, Guide.xlabel("P"),
    Guide.ylabel("rmse"), Guide.title("Root Mean Squared Error"),
    Guide.manual_color_key("", ["σ² = 1", "σ² = 10"], ["red", "blue"]));
draw(PNG("test.png", 15inch, 15inch), pppp)

Gadfly.pop_theme()

draw(PNG("myplot2.png", 15inch, 15inch), myplot2)

p = 10
x, y = GB6.test_data(500; p = 10, sigma = sqrt(10.0))
x_test, y_test = GB6.test_data(500; p = 10, sigma = sqrt(10.0))
mu = GB6.fried(x)
mu_test = GB6.fried(x_test)
x=DataFrame(x)
y=DataFrame(y=y)
x_test=DataFrame(x_test)
y_test=DataFrame(y_test=y_test)
mu = DataFrame(mu = mu)
mu_test = DataFrame(mu_test = mu_test)

## for comparing the data with bencnmark method, use CSV

CSV.write("_X.csv", x)
CSV.write("_Y.csv", y)
CSV.write("_X_test.csv", x_test)
CSV.write("_Y_test.csv", y_test)
CSV.write("_mu.csv", mu)
CSV.write("_mu_test.csv", mu_test)

for j in 1 : size(x, 2)
    x[:, j] = collect(skipmissing(x[:,j]))
end
x = convert(Matrix, x)
y = collect(skipmissing(y[:,1]))

for j in 1 : size(x_test, 2)
    x_test[:, j] = collect(skipmissing(x_test[:,j]))
end
x_test = convert(Matrix, x_test)
y_test = collect(skipmissing(y_test[:,1]))

mu_test = collect(skipmissing(mu_test[:,1]))

gbdt = GB6.GBDT(;
  sampling_rate = 0.75,                # Sampling rate
  learning_rate = 0.1,                # Learning rate
  num_iterations = 1,               # NuGBer of iterations
  tree_options = Dict(
      :nsubfeatures => 0,
      :min_samples_leaf => 1,
      :max_depth => 5))
gbl = GB6.GBLearner(gbdt, :regression)

hypers = GB6.Hypers()
hypers.σ = sqrt(1)
hypers.τ = 3 / 10 * var(y) / gbdt.num_iterations
mymodel=GB6.fit!(gbl, x, y, hypers)
# get the fitted value
myfit  = GB6.predict!(gbl, x)

# get the residual sum of squares
sqrt(mean(y .- myfit) .^2) # 0.0103

# get prediction
mypred = GB6.predict!(gbl, x_test)

# get mse and root of mean squared error
sqrt(mean((y_test .- mypred ) .^2)) #10.2

sqrt(mean((mypred .- mu_test ) .^2))

#### compare with R
# σ² = 10, σ = √10
xx = CSV.read("xx.csv")[:, 2:end]
yy = CSV.read("yy.csv")[:, 2]
xx_test = CSV.read("xx_test.csv")[:, 2:end]
yy_test = CSV.read("yy_test.csv")[:, 2]
μμ_test = CSV.read("mumu_test.csv")[:, 2]

for j in 1 : size(xx, 2)
    xx[:, j] = collect(skipmissing(xx[:,j]))
end
xx = convert(Matrix, xx)
yy = collect(skipmissing(yy[:,1]))

for j in 1 : size(xx_test, 2)
    xx_test[:, j] = collect(skipmissing(xx_test[:,j]))
end
xx_test = convert(Matrix, xx_test)
yy_test = collect(skipmissing(yy_test[:,1]))

μμ_test = collect(skipmissing(μμ_test[:,1]))

# set the same combination of parameters, though I did not use cross validation here
gbdt_new = GB6.GBDT(;
  sampling_rate = 0.8,                # Sampling rate
  learning_rate = 0.1,                # Learning rate
  num_iterations = 100,               # NuGBer of iterations
  tree_options = Dict(
      :nsubfeatures => 0,
      :min_samples_leaf => 5,
      :max_depth => 6))
gbl_new = GB6.GBLearner(gbdt_new, :regression)

mymodel2=GB6.fit!(gbl_new, xx, yy)
# get the fitted value
myfit2 = GB6.predict!(gbl_new, xx)

# get the residual sum of squares
sqrt(mean(yy .- myfit2) .^2) # 0.02446; 0.0045(σ = √10)

# get prediction
mypred2 = GB6.predict!(gbl_new, xx_test)

# get mse and root of mean squared error
mean((yy_test .- mypred2 ) .^2) # ; 14.1986(σ = √10)
sqrt(mean((yy_test .- mypred2 ) .^2)) #12.02; 3.7678(σ = √10)

sqrt(mean((mypred2 .- μμ_test ) .^2)) # 6.6 ; 1.93(σ = 1.9)

max_stars_repo_name: YinpuLi/SparsityBoosting

max_issues_repo_name: YinpuLi/SparsityBoosting

max_forks_repo_name: YinpuLi/SparsityBoosting


text:
using LinearAlgebra
function Vectorin(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,x5,y5,z5)

    v = zeros(Float64,6,3)
    CD = [x2-x1;y2-y1;z2-z1]
    CH = [x4-x1;y4-y1;z4-z1]
    na = LinearAlgebra.cross(CD,CH)
    BC = [x1-x3;y1-y3;z1-z3]
    j1 = dot(na,BC)
    if j1 > 0
        n1 = na
    else
        n1 = -na
    end
    union1 = n1/sqrt(n1[1]^2+n1[2]^2+n1[3]^2)
    v[1,1] = union1[1]
    v[1,2] = union1[2]
    v[1,3] = union1[3]
    DH = [x4-x2;y4-y2;z4-z2]
    DB = [x3-x2;y3-y2;z3-z2]
    nb = LinearAlgebra.cross(DH,DB)
    CD = [x2-x1;y2-y1;z2-z1]
    j2 = dot(nb,CD)
    if j2 > 0
        n2 = nb
    else
        n2 = -nb
    end
    union2 = n2/sqrt(n2[1]^2+n2[2]^2+n2[3]^2)
    v[2,1] = union2[1]
    v[2,2] = union2[2]
    v[2,3] = union2[3]
    HB = [x3-x4;y3-y4;z3-z4]
    HC = [x1-x4;y1-y4;z1-z4]
    nc = LinearAlgebra.cross(HB,HC)
    DH = [x4-x2;y4-y2;z4-z2]
    j3 = dot(nc,DH)
    if j3 > 0
        n3 = nc
    else
        n3 = -nc
    end
    union3 = n3/sqrt(n3[1]^2+n3[2]^2+n3[3]^2)
    v[3,1] = union3[1]
    v[3,2] = union3[2]
    v[3,3] = union3[3]
    CD = [x2-x1;y2-y1;z2-z1]
    CI = [x5-x1;y5-y1;z5-z1]
    nd = LinearAlgebra.cross(CD,CI)
    BC = [x1-x3;y1-y3;z1-z3]
    j4 = dot(nd,BC)
    if j4 > 0
        n4 = nd
    else
        n4 = -nd
    end
    union4 = n4/sqrt(n4[1]^2+n4[2]^2+n4[3]^2)
    v[4,1] = union4[1]
    v[4,2] = union4[2]
    v[4,3] = union4[3]
    BD = [x2-x3;y2-y3;z2-z3]
    BI = [x5-x3;y5-y3;z5-z3]
    ne = LinearAlgebra.cross(BD,BI)
    CB = [x3-x1;y3-y1;z3-z1]
    j5 = dot(ne,CB)
    if j5 > 0
        n5 = ne
    else
        n5 = -ne
    end
    union5 = n5/sqrt(n5[1]^2+n5[2]^2+n5[3]^2)
    v[5,1] = union5[1]
    v[5,2] = union5[2]
    v[5,3] = union5[3]
    IB = [x3-x5;y3-y5;z3-z5]
    IC = [x1-x5;y1-y5;z1-z5]
    nf = LinearAlgebra.cross(IB,IC)
    DI = [x5-x2;y5-y2;z5-z2]
    j6 = dot(nf,DI)
    if j6 > 0
        n6 = nf
    else
        n6 = -nf
    end
    union6 = n6/sqrt(n6[1]^2+n6[2]^2+n6[3]^2)
    v[6,1] = union6[1]
    v[6,2] = union6[2]
    v[6,3] = union6[3]

    return v
end

max_stars_repo_path: src/single-core juSFEM/vectorin.jl

max_stars_repo_name: Elliothuo/juSFEM

max_issues_repo_path: src/single-core juSFEM/vectorin.jl

max_issues_repo_name: Elliothuo/juSFEM

max_forks_repo_path: src/single-core juSFEM/vectorin.jl

max_forks_repo_name: Elliothuo/juSFEM


text:
kolmogorov_smirnov(cdf_1_evaluated, pm_1, cdf_2_evaluated, pm_2, grid) = maximum(abs.(cdf_1_evaluated-cdf_2_evaluated))
function L1(f1_evaluated, pm_1, f2_evaluated, pm_2, grid)
    d = f1_evaluated-f2_evaluated
    d = sum(abs.(d[1:end-1,:] + d[2:end,:])./2.0.*diff(grid))
    return log10(d)
end
function L2(f1_evaluated, pm_1, f2_evaluated, pm_2, grid)
    d = f1_evaluated-f2_evaluated
    d = sum(((d[1:end-1,:] + d[2:end,:]).^2)./2.0.*diff(grid))
    return log10(d)
end
function L1_cell_probs(f1_evaluated, pm_1, f2_evaluated, pm_2, grid)
    d = sum(abs.(f1_evaluated-f2_evaluated))
    d += sum(abs.(pm_1[:]-pm_2[:]))
    return log10(d)
end


text:
"""
    factorial_iterative(n)

Finds factorial of a number using Iterative method

# Example
```julia
factorial_iterative(5)      # returns 120
factorial_iterative(0.1)    # returns error
factorial_iterative(-1)     # returns error
```
# Reference
- factorial of a positive integer -- https://en.wikipedia.org/wiki/Factorial

Contributed By:- [Ashwani Rathee](https://github.com/ashwani-rathee)
"""
function factorial_iterative(n)
    if n != trunc(n) || n < 0
        throw(error("factorial_iterative() only accepts non-negative integral values"))
    end
    factorial::BigInt = 1
    map(i -> factorial *= i, 1:n)
    return factorial
end

"""
    factorial_recursive(n)

Finds factorial of anumber using recursive method

# Example
```julia
factorial_recursive(5)      # returns 120
factorial_recursive(0.1)    # returns error
factorial_recursive(-1)     # returns error
```
# Reference
- factorial of a positive integer -- https://en.wikipedia.org/wiki/Factorial

Contributed By:- [Ashwani Rathee](https://github.com/ashwani-rathee)
"""
function factorial_recursive(n)
    if n != trunc(n) || n < 0
        throw(error("factorial_iterative() only accepts non-negative integral values"))
    end
    if n == 0 || n == 1
        return 1
    else
        factorial::BigInt = n * factorial_recursive(n - 1)
        return factorial
    end
end

max_stars_repo_name: ashwani-rathee/Julia

max_issues_repo_name: ashwani-rathee/Julia

max_forks_repo_name: ashwani-rathee/Julia


text:
using Test

@testset "/src/ConvexPolytope.jl" begin

using BellScenario: ConvexPolytope

@testset "linear_interpolation()" begin

    @test isapprox(ConvexPolytope.linear_interpolation([1;0;0],[0;1;0],dw=0.2), map(
        x -> reshape(x,(3,1)),
        [[1.0;0;0],[0.8;0.2;0],[0.6;0.4;0],[0.4;0.6;0],[0.2;0.8;0],[0;1.0;0]]
    ))

    @test all(map(
        x -> ConvexPolytope.facet_contains([0 0 0 -1], x),
        ConvexPolytope.linear_interpolation([1;0;0],[0;1;0],dw=0.2)
    ))
end

@testset "facet_contains()" begin

    @test ConvexPolytope.facet_contains([-1 1 1 1],[1/3;1/3;1/3])
    @test !(ConvexPolytope.facet_contains([-1 1 1 1],[1;1;1]))

    @test ConvexPolytope.facet_contains([-1 1 0 -1 0 1 -1],[2/3;1/6;1/6;1/6;2/3;1/6])
end

@testset "violated_facets()" begin

    facets = [[0 -1 0 0],[0 0 -1 0],[0 0 0 -1],[-1 1 1 1]]
    @test ConvexPolytope.violated_facets(facets, [1;1;1]) == [[-1 1 1 1]]
    @test ConvexPolytope.violated_facets(facets, [1/4;1/4;1/4]) == []

    facets = [
        [0 -1 0 0 0 0 0],[0 0 -1 0 0 0 0],[0 0 0 -1 0 0 0],[0 0 0 0 -1 0 0],
        [0 0 0 0 0 -1 0],[0 0 0 0 0 0 -1],[-1 0 0 1 0 0 1],[-1 0 1 0 0 1 0],
        [-1 1 0 0 1 0 0],[-1 -1 0 1 -1 1 0],[-1 -1 1 0 -1 0 1],[-1 0 -1 1 1 -1 0],
        [-1 0 1 -1 1 0 -1], [-1 1 -1 0 0 -1 1], [-1 1 0 -1 0 1 -1]
    ]
    @test ConvexPolytope.violated_facets(facets, [1;0;0;0;1;0]) == [[-1 1 0 -1 0 1 -1]]
end

@testset "facet_violation()" begin

    @test ConvexPolytope.facet_violation([-1 1 0 -1 0 1 -1], [1;1;0;0;0;1;0]) == 1
    @test ConvexPolytope.facet_violation([-1 1 0 -1 0 1 -1], [1;0;0;0;1;0]) == 1

    @test_throws ArgumentError ConvexPolytope.facet_violation([-1 1 0 -1 0 1 -1],[1;1])

    @test ConvexPolytope.facet_violation([-1 1 1 1], [1;1;1]) == 2
end

@testset "facet_distance()" begin
    @test ConvexPolytope.facet_distance([0 -1 0 0],[1;1;1]) ≈ -1
    @test ConvexPolytope.facet_distance([-1 1 1 1], [1;1;1]) ≈ 2/sqrt(3)
end

@testset "contains()" begin

    vertices = [[0;0;0],[1;0;0],[0;1;0],[0;0;1]]
    facets = [[0 -1 0 0],[0 0 -1 0],[0 0 0 -1],[-1 1 1 1]]

    @test !(ConvexPolytope.facets_contain(facets, [1;1;1]))
    @test ConvexPolytope.facets_contain(facets,[1/4;1/4;1/4])
    @test ConvexPolytope.facets_contain(facets,[1/3;1/3;1/3])
    @test all(map(
        (v) -> ConvexPolytope.facets_contain(facets, v),
        vertices
    ))

end

@testset "collect_vertices()" begin
    facets = [[-1 1 0 -1 0 1 -1],[0 -1 0 0 0 0 0]]
    vertices = [
        [1;0;0;0;0;0;0],[1;1;1;1;0;0;0],[1;0;0;0;1;1;1],[1;1;1;0;0;0;0],
        [1;1;1;0;0;0;1],[1;1;0;1;0;0;0],[1;1;0;1;0;1;0],[1;0;1;1;0;0;0],
        [1;0;1;1;1;0;0]
    ]

    facet_dict = ConvexPolytope.collect_vertices(facets, vertices)

    @test length(facet_dict) == 2
    @test facet_dict[1]["facet"] == [-1 1 0 -1 0 1 -1]
    @test facet_dict[1]["vertices"] == [[1;1;1;0;0;0;0],[1;1;0;1;0;1;0]]
    @test facet_dict[2]["facet"] == [0 -1 0 0 0 0 0]
    @test facet_dict[2]["vertices"] == [[1;0;0;0;0;0;0],[1;0;0;0;1;1;1],[1;0;1;1;0;0;0],[1;0;1;1;1;0;0]]
end

@testset "collect_facets()" begin
    facets = [[-1 1 0 -1 0 1 -1],[0 -1 0 0 0 0 0]]
    vertices = [
        [1;0;0;0;0;0;0],[1;0;0;0;1;1;1],[1;1;1;0;0;0;1],[1;0;0;0;1;1;0]
    ]

    facet_dict = ConvexPolytope.collect_facets(vertices, facets)

    @test length(facet_dict) == 4
    @test facet_dict[1]["vertex"] == [1;0;0;0;0;0;0]
    @test facet_dict[1]["facets"] == [[0 -1 0 0 0 0 0]]
    @test facet_dict[2]["vertex"] == [1;0;0;0;1;1;1]
    @test facet_dict[2]["facets"] == [[0 -1 0 0 0 0 0]]
    @test facet_dict[3]["vertex"] == [1;1;1;0;0;0;1]
    @test facet_dict[3]["facets"] == []
    @test facet_dict[4]["vertex"] == [1;0;0;0;1;1;0]
    @test facet_dict[4]["facets"] == [[-1 1 0 -1 0 1 -1],[0 -1 0 0 0 0 0]]
end

@testset "contains()" begin
    facets = [[0 -1 -1 -1 0 0 0], [0 0 0 0 -1 -1 -1], [0 0 0 0 0 0 0]]

    @test ConvexPolytope.contains(facets, [1;0;0;0;0;0;0])
    @test ConvexPolytope.contains(facets, [1;1;1;1;0;0;0])
    @test ConvexPolytope.contains(facets, [1;0.5;0.5;0.5;0.5;0.5;0.5])
    @test !ConvexPolytope.contains(facets, [1;-1;0;0;0;0;0])
end

end


text:
using Test

using Spglib

function type2dict(dt)
    di = Dict{Symbol,Any}()
    for n in propertynames(dt)
        di[n] = getproperty(dt, n)
    end
    di
end

@testset "Construct a `Cell`" begin
    lattice = [
        [5.0759761474456697, 5.0759761474456697, 0],  # a
        [-2.8280307701821314, 2.8280307701821314, 0],  # b
        [0, 0, 8.57154746],  # c
    ]
    positions = [
        [0.0, 0.84688439, 0.1203133],
        [0.0, 0.65311561, 0.6203133],
        [0.0, 0.34688439, 0.3796867],
        [0.0, 0.15311561, 0.8796867],
        [0.5, 0.34688439, 0.1203133],
        [0.5, 0.15311561, 0.6203133],
        [0.5, 0.84688439, 0.3796867],
        [0.5, 0.65311561, 0.8796867],
    ]
    types = fill(35, length(positions))
    cell = Cell(lattice, positions, types)
    @test cell == Cell(
        [
            5.07597614744567 -2.8280307701821314 0.0
            5.07597614744567 2.8280307701821314 0.0
            0.0 0.0 8.57154746
        ],
        positions,
        types,
    )
    @test cell == Cell(
        lattice,
        [
            0.0 0.0 0.0 0.0 0.5 0.5 0.5 0.5
            0.84688439 0.65311561 0.34688439 0.15311561 0.34688439 0.15311561 0.84688439 0.65311561
            0.1203133 0.6203133 0.3796867 0.8796867 0.1203133 0.6203133 0.3796867 0.8796867
        ],
        types,
    )
    @test cell == Cell(
        [
            5.07597614744567 -2.8280307701821314 0.0
            5.07597614744567 2.8280307701821314 0.0
            0.0 0.0 8.57154746
        ],
        [
            0.0 0.0 0.0 0.0 0.5 0.5 0.5 0.5
            0.84688439 0.65311561 0.34688439 0.15311561 0.34688439 0.15311561 0.84688439 0.65311561
            0.1203133 0.6203133 0.3796867 0.8796867 0.1203133 0.6203133 0.3796867 0.8796867
        ],
        types,
    )
end

@testset "Test `get_spacegroup_type`" begin
    # Adapted from https://github.com/unkcpz/LibSymspg.jl/blob/53d2f6d/test/test_api.jl#L7-L12
    spacegroup_type = get_spacegroup_type(101)
    @test spacegroup_type.number == 15
    @test spacegroup_type.hall_symbol == "-I 2a"
    @test spacegroup_type.arithmetic_crystal_class_symbol == "2/mC"
    # These results are compared with Python's spglib results.
    @test type2dict(get_spacegroup_type(419)) == Dict(
        :number => 136,
        :international_short => "P4_2/mnm",
        :international_full => "P 4_2/m 2_1/n 2/m",
        :international => "P 4_2/m n m",
        :schoenflies => "D4h^14",
        :hall_symbol => "-P 4n 2n",
        :choice => "",
        :pointgroup_schoenflies => "D4h",
        :pointgroup_international => "4/mmm",
        :arithmetic_crystal_class_number => 36,
        :arithmetic_crystal_class_symbol => "4/mmmP",
    )
    @test type2dict(get_spacegroup_type(1)) == Dict(
        :number => 1,
        :international_short => "P1",
        :international_full => "P 1",
        :international => "P 1",
        :schoenflies => "C1^1",
        :hall_symbol => "P 1",
        :choice => "",
        :pointgroup_schoenflies => "C1",
        :pointgroup_international => "1",
        :arithmetic_crystal_class_number => 1,
        :arithmetic_crystal_class_symbol => "1P",
    )
    @test type2dict(get_spacegroup_type(525)) == Dict(
        :number => 227,
        :international_short => "Fd-3m",
        :international_full => "F 4_1/d -3 2/m",
        :international => "F d -3 m",
        :schoenflies => "Oh^7",
        :hall_symbol => "F 4d 2 3 -1d",
        :choice => "1",
        :pointgroup_schoenflies => "Oh",
        :pointgroup_international => "m-3m",
        :arithmetic_crystal_class_number => 72,
        :arithmetic_crystal_class_symbol => "m-3mF",
    )
    @test type2dict(get_spacegroup_type(485)) == Dict(
        :number => 191,
        :international_short => "P6/mmm",
        :international_full => "P 6/m 2/m 2/m",
        :international => "P 6/m m m",
        :schoenflies => "D6h^1",
        :hall_symbol => "-P 6 2",
        :choice => "",
        :pointgroup_schoenflies => "D6h",
        :pointgroup_international => "6/mmm",
        :arithmetic_crystal_class_number => 58,
        :arithmetic_crystal_class_symbol => "6/mmm",
    )
end

# @testset "Test MgB2 structure" begin
#     a = 3.07
#     c = 3.52
#     lattice = [
#         a 0 0
#         -a/2 a/2*sqrt(3) 0
#         0 0 c
#     ]
#     positions = [
#         0 0 0
#         1.0/3 2.0/3 0.5
#         2.0/3 1.0/3 0.5
#     ]
#     types = [12, 5, 5]
#     MgB2 = Cell(lattice, positions, types)
# end

# From https://github.com/unkcpz/LibSymspg.jl/blob/53d2f6d/test/test_api.jl#L34-L77
@testset "Get symmetry operations" begin
    @testset "Normal symmetry" begin
        lattice = [
            4.0 0.0 0.0
            0.0 4.0 0.0
            0.0 0.0 4.0
        ]
        positions = [
            0.0 0.5
            0.0 0.5
            0.0 0.5
        ]
        types = [1, 1]
        cell = Cell(lattice, positions, types)
        num_atom = length(types)
        max_size = num_atom * 48
        rotation = Array{Cint,3}(undef, 3, 3, max_size)
        translation = Array{Float64,2}(undef, 3, max_size)
        get_symmetry!(rotation, translation, cell, 1e-5)
        @test size(rotation) == (3, 3, 96)
        @test size(translation) == (3, 96)
        @test get_hall_number_from_symmetry(rotation, translation, max_size, 1e-5) == 529
    end
    # See https://github.com/spglib/spglib/blob/deb6695/python/test/test_collinear_spin.py#L18-L37
    @testset "Get symmetry with collinear spins" begin
        lattice = [
            4.0 0.0 0.0
            0.0 4.0 0.0
            0.0 0.0 4.0
        ]
        positions = [
            0.0 0.5
            0.0 0.5
            0.0 0.5
        ]
        types = [1, 1]
        @testset "Test ferromagnetism" begin
            magmoms = [1.0, 1.0]
            cell = Cell(lattice, positions, types, magmoms)
            rotation, translation, equivalent_atoms =
                get_symmetry_with_collinear_spin(cell, 1e-5)
            @test size(rotation) == (3, 3, 96)
            @test equivalent_atoms == [0, 0]
        end
        @testset "Test antiferromagnetism" begin
            magmoms = [1.0, -1.0]
            cell = Cell(lattice, positions, types, magmoms)
            rotation, translation, equivalent_atoms =
                get_symmetry_with_collinear_spin(cell, 1e-5)
            @test size(rotation) == (3, 3, 96)
            @test equivalent_atoms == [0, 0]
        end
        @testset "Test broken magmoms" begin
            magmoms = [1.0, 2.0]
            cell = Cell(lattice, positions, types, magmoms)
            rotation, translation, equivalent_atoms =
                get_symmetry_with_collinear_spin(cell, 1e-5)
            @test size(rotation) == (3, 3, 48)
            @test size(translation) == (3, 48)
            @test equivalent_atoms == [0, 1]
        end
    end

    @testset "Get multiplicity" begin
        lattice = [
            4.0 0.0 0.0
            0.0 4.0 0.0
            0.0 0.0 4.0
        ]
        positions = [
            0.0 0.5
            0.0 0.5
            0.0 0.5
        ]
        types = [1, 1]
        cell = Cell(lattice, positions, types)
        @test get_multiplicity(cell, 1e-5) == 96
    end
end


text:
#===============================================================================
# Promotion and conversion in Julia is no-automatic and uses multiple dispatch
#
# Author: Carl Cortright
# Date:12/31/2016
#
===============================================================================#
import Base.convert
import Base.promote_rule

# The convert function takes a type and object as arguments
convert(Float64, 20)

# Defining new convert functions is easy
convert(::Type{Bool}, x::Real) = x == 0 ? false : x==1 ? true : error("Must provide 1 or 0")
convert(Bool, 1)

# Promotion is desiged to convert to a "greater" type, mostly used on numerics
promote(1,2.5)

# Instead of permuting all promote declarations, we just create promote rules
promote_rule(::Type{Float32}, ::Type{Int16}) = Float64

max_issues_repo_name: carlcortright/Playground


text:
##################################################################################################
# Integration with UncertainData.jl (sampling from full supports of the furnishing distributions)
#################################################################################################

uvals_x = [UncertainValue(Normal, rand(Normal(0, 5)), abs(rand(Normal(0, 3)))) for i = 1:100]
uvals_y = [UncertainValue(Normal, rand(Normal(0, 5)), abs(rand(Normal(0, 3)))) for i = 1:100];

# UncertainValueDataset
UVX = UncertainValueDataset(uvals_x)
UVY = UncertainValueDataset(uvals_y)

# UncertainIndexDataset
UVX_idx = UncertainIndexDataset(uvals_x)
UVY_idx = UncertainIndexDataset(uvals_y)

# Real-valued vectors
x = resample.(uvals_x);
y = resample.(uvals_y);
binnings = [RectangularBinning(i) for i = 3:6]
vftest = VisitationFrequencyTest(binning = binnings, ηs = 1)

@test causality(x, y, vftest) isa Array{<:Real, 0}
@test causality(uvals_x, uvals_y, vftest) isa Array{<:Real, 0}
@test causality(x, uvals_y, vftest) isa Array{<:Real, 0}
@test causality(uvals_x, y, vftest) isa Array{<:Real, 0}
@test causality(UVX, UVY, vftest) isa Array{<:Real, 0}
@test causality(x, UVY, vftest) isa Array{<:Real, 0}
@test causality(UVX, y, vftest) isa Array{<:Real, 0}

vftest = VisitationFrequencyTest(binning = binnings, ηs = 1:5)

@test causality(x, y, vftest) isa Array{T, 1} where T
@test causality(uvals_x, uvals_y, vftest) isa Array{T, 1} where T
@test causality(x, uvals_y, vftest) isa Array{T, 1} where T
@test causality(uvals_x, y, vftest) isa Array{T, 1} where T
@test causality(UVX, UVY, vftest) isa Array{T, 1} where T
@test causality(x, UVY, vftest) isa Array{T, 1} where T
@test causality(UVX, y, vftest) isa Array{T, 1} where T

@test causality(x, y, vftest) |> length == 5
@test causality(uvals_x, uvals_y, vftest) |> length == 5
@test causality(x, uvals_y, vftest) |> length == 5
@test causality(uvals_x, y, vftest) |> length == 5
@test causality(UVX, UVY, vftest) |> length == 5
@test causality(x, UVY, vftest) |> length == 5
@test causality(UVX, y, vftest) |> length == 5

################################################################
# Integration with UncertainData.jl (with sampling constraints)
################################################################
onevar_constraints = ConstrainedResampling(TruncateStd(1))
twovar_constraints = ConstrainedResampling(TruncateStd(2), TruncateStd(1))

test_vf = VisitationFrequencyTest(binning = RectangularBinning(5), ηs = 1)
@test causality(x, y, test_vf, twovar_constraints) isa Array{<:Real, 0}
@test causality(uvals_x, uvals_y, test_vf, twovar_constraints) isa Array{<:Real, 0}
@test causality(x, uvals_y, test_vf, onevar_constraints) isa Array{<:Real, 0}
@test causality(uvals_x, y, test_vf, onevar_constraints) isa Array{<:Real, 0}
@test causality(UVX, UVY, test_vf, twovar_constraints) isa Array{<:Real, 0}
@test causality(uvals_x, UVY, test_vf, twovar_constraints) isa Array{<:Real, 0}
@test causality(UVX, uvals_y, test_vf, twovar_constraints) isa Array{<:Real, 0}
@test causality(x, UVY, test_vf, onevar_constraints) isa Array{<:Real, 0}
@test causality(UVX, y, test_vf, onevar_constraints) isa Array{<:Real, 0}

test_vf = VisitationFrequencyTest(binning = RectangularBinning(5), ηs = -3:3)
@test causality(x, y, test_vf, twovar_constraints) isa Array{<:Real, 1}
@test causality(uvals_x, uvals_y, test_vf, twovar_constraints) isa Array{<:Real, 1}
@test causality(x, uvals_y, test_vf, onevar_constraints) isa Array{<:Real, 1}
@test causality(uvals_x, y, test_vf, onevar_constraints) isa Array{<:Real, 1}
@test causality(UVX, UVY, test_vf, twovar_constraints) isa Array{<:Real, 1}
@test causality(uvals_x, UVY, test_vf, twovar_constraints) isa Array{<:Real, 1}
@test causality(UVX, uvals_y, test_vf, twovar_constraints) isa Array{<:Real, 1}
@test causality(x, UVY, test_vf, onevar_constraints) isa Array{<:Real, 1}
@test causality(UVX, y, test_vf, onevar_constraints) isa Array{<:Real, 1}


text:
module ISA

include("components/basicComps.jl")
export AMFMcomp

include("components/numericalComps.jl")
include("demod/demod.jl")
export derivApprox, AMFMdemod, AMFMcompN

include("models/basicModels.jl")
export AMFMmodel, fourierSeries

include("plots/isaPlots.jl")
export isaPlot3d

include("decomp/emd.jl")
export findLocalMaxima, SIFT, EMD, ℂSIFT, ℂEMD

include("TFDs/TFDs.jl")
export frameSignal, STFT

end # module

max_stars_repo_name: kricheso/ISA

max_issues_repo_name: kricheso/ISA

max_forks_repo_name: kricheso/ISA