text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    mcismol, mcismolsize,
    mcesmol, mcesmolsize

"""
    mcismol(mol1::UndirectedGraph, mol2::UndirectedGraph; kwargs...
        ) -> Tuple{Dict{Int,Int},Symbol}

Compute maximum common induced substructure (MCIS) of mol1 and mol2.

## Keyword arguments

- connected(Bool): if true, apply connected MCS constraint.
- topological(Bool): if true, apply topological constraint.
- diameter(Int): distance cutoff for topological constraint.
- tolerance(Int): distance mismatch tolerance for topological constraint.
- timeout(Int): abort calculation and return suboptimal results if the execution
time has reached the given value (default=60, in seconds).
- targetsize(Int): abort calculation and return suboptimal result so far if the
given mcs size achieved.

# References

1. Kawabata, T. (2011). Build-Up Algorithm for Atomic Correspondence between
Chemical Structures. Journal of Chemical Information and Modeling, 51(8),
1775–1787. https://doi.org/10.1021/ci2001023
1. https://www.jstage.jst.go.jp/article/ciqs/2017/0/2017_P4/_article/-char/en
"""
function mcismol(mol1::UndirectedGraph, mol2::UndirectedGraph; kwargs...)
    afunc = atommatch(mol1, mol2)
    bfunc = bondmatch(mol1, mol2)
    return findmcis(
        mol1, mol2, nodematcher=afunc, edgematcher=bfunc; kwargs...)
end

mcismolsize(mol1, mol2; kwargs...) = length(mcismol(mol1, mol2; kwargs...)[1])

"""
    mcesmol(mol1::UndirectedGraph, mol2::UndirectedGraph; kwargs...
        ) -> Tuple{Dict{Int,Int},Symbol}

Compute maximum common edge induced substructure (MCES) of mol1 and mol2.
"""
function mcesmol(mol1::UndirectedGraph, mol2::UndirectedGraph; kwargs...)
    afunc = atommatch(mol1, mol2)
    bfunc = bondmatch(mol1, mol2)
    return findmces(
        mol1, mol2, nodematcher=afunc, edgematcher=bfunc; kwargs...)
end

mcesmolsize(mol1, mol2; kwargs...) = length(mcesmol(mol1, mol2; kwargs...)[1])


text:
function SatelliteAntenna(arg0::JString, arg1::JString, arg2::Map, arg3::SatelliteSystem, arg4::jint, arg5::SatelliteType, arg6::jint, arg7::JString, arg8::AbsoluteDate, arg9::AbsoluteDate)
    return SatelliteAntenna((JString, JString, Map, SatelliteSystem, jint, SatelliteType, jint, JString, AbsoluteDate, AbsoluteDate), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function get_cospar_id(obj::SatelliteAntenna)
    return jcall(obj, "getCosparID", JString, ())
end

function get_prn_number(obj::SatelliteAntenna)
    return jcall(obj, "getPrnNumber", jint, ())
end

function get_satellite_code(obj::SatelliteAntenna)
    return jcall(obj, "getSatelliteCode", jint, ())
end

function get_satellite_system(obj::SatelliteAntenna)
    return jcall(obj, "getSatelliteSystem", SatelliteSystem, ())
end

function get_satellite_type(obj::SatelliteAntenna)
    return jcall(obj, "getSatelliteType", SatelliteType, ())
end

function get_valid_from(obj::SatelliteAntenna)
    return jcall(obj, "getValidFrom", AbsoluteDate, ())
end

function get_valid_until(obj::SatelliteAntenna)
    return jcall(obj, "getValidUntil", AbsoluteDate, ())
end

ext: jl


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module TestConstraintGeomeanToRelentr

using Test

using MathOptInterface
const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

include("../utilities.jl")

function test_GeoMeantoRelEntr()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.GeoMeantoRelEntr{Float64}(mock)
    MOI.Test.test_basic_VectorOfVariables_GeometricMeanCone(
        bridged_mock,
        config,
    )
    MOI.empty!(bridged_mock)
    MOI.Test.test_basic_VectorAffineFunction_GeometricMeanCone(
        bridged_mock,
        config,
    )
    MOI.empty!(bridged_mock)
    return MOI.Test.test_basic_VectorQuadraticFunction_GeometricMeanCone(
        bridged_mock,
        config,
    )
end

function test_conic_GeometricMeanCone_VectorAffineFunction()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.GeoMeantoRelEntr{Float64}(mock)
    var_primal = [1, 1, 1, 1, 0]
    relentr_dual = [1, 1, 1, 1, -1, -1, -1] / 3
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            var_primal,
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [-inv(3)],
            (MOI.VectorOfVariables, MOI.Nonnegatives) => [[1]],
            (MOI.VectorAffineFunction{Float64}, MOI.RelativeEntropyCone) =>
                [relentr_dual],
        )

    MOI.Test.test_conic_GeometricMeanCone_VectorAffineFunction(
        bridged_mock,
        config,
    )
    var_names = ["t", "x", "y", "z"]
    MOI.set(
        bridged_mock,
        MOI.VariableName(),
        MOI.get(bridged_mock, MOI.ListOfVariableIndices()),
        var_names,
    )
    nonneg = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{MOI.VectorOfVariables,MOI.Nonnegatives}(),
    )
    relentr = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.RelativeEntropyCone,
        }(),
    )
    aux = MOI.get(mock, MOI.ListOfVariableIndices())[end]
    MOI.set(mock, MOI.VariableName(), aux, "aux")
    @test length(nonneg) == 1
    MOI.set(mock, MOI.ConstraintName(), nonneg[1], "nonneg")
    @test length(relentr) == 1
    MOI.set(mock, MOI.ConstraintName(), relentr[1], "relentr")
    less = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{Float64},
            MOI.LessThan{Float64},
        }(),
    )
    @test length(less) == 1
    MOI.set(mock, MOI.ConstraintName(), less[1], "less")

    s = """
    variables: t, x, y, z, aux
    less: x + y + z in MathOptInterface.LessThan(3.0)
    nonneg: [aux] in MathOptInterface.Nonnegatives(1)
    relentr: [0.0, x, y, z, t + aux, t + aux, t + aux] in MathOptInterface.RelativeEntropyCone(7)
    maxobjective: t
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        mock,
        model,
        vcat(var_names, "aux"),
        ["less", "nonneg", "relentr"],
    )
    geomean = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.GeometricMeanCone,
        }(),
    )
    @test length(geomean) == 1
    MOI.set(bridged_mock, MOI.ConstraintName(), geomean[1], "geomean")
    less = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{Float64},
            MOI.LessThan{Float64},
        }(),
    )
    @test length(less) == 1
    MOI.set(bridged_mock, MOI.ConstraintName(), less[1], "less")

    s = """
    variables: t, x, y, z
    less: x + y + z in MathOptInterface.LessThan(3.0)
    geomean: [1.0t, x, y, z] in MathOptInterface.GeometricMeanCone(4)
    maxobjective: t
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        bridged_mock,
        model,
        var_names,
        ["less", "geomean"],
    )
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.GeometricMeanCone,
            }(),
        ),
    )
    for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        value =
            (attr isa MOI.ConstraintPrimalStart) ? ones(4) :
            vcat(-1, fill(inv(3), 3))
        MOI.set(bridged_mock, attr, ci, value)
        @test MOI.get(bridged_mock, attr, ci) ≈ value
        if attr isa MOI.ConstraintPrimalStart
            @test MOI.get(mock, MOI.VariablePrimalStart(), aux) == 0
            @test MOI.get(mock, attr, nonneg[1]) == [0]
            @test MOI.get(mock, attr, relentr[1]) == [0, 1, 1, 1, 1, 1, 1]
        else
            @test MOI.get(mock, attr, nonneg[1]) == [1]
            @test MOI.get(mock, attr, relentr[1]) ≈ relentr_dual
        end
    end
    _test_delete_bridge(
        bridged_mock,
        ci,
        4,
        (
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}, 1),
            (MOI.VectorOfVariables, MOI.Nonnegatives, 0),
            (MOI.VectorAffineFunction{Float64}, MOI.RelativeEntropyCone, 0),
        ),
    )
    return
end

function test_conic_GeometricMeanCone_VectorAffineFunction_2()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.GeoMeantoRelEntr{Float64}(mock)
    var_primal = vcat(ones(10), 0)
    relentr_dual = vcat(fill(inv(9), 10), fill(-inv(9), 9))
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            var_primal,
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) =>
                fill(-inv(9), 9),
            (MOI.VectorOfVariables, MOI.Nonnegatives) => [[1]],
            (MOI.VectorAffineFunction{Float64}, MOI.RelativeEntropyCone) =>
                [relentr_dual],
        )

    MOI.Test.test_conic_GeometricMeanCone_VectorAffineFunction_2(
        bridged_mock,
        config,
    )
    var_names = vcat("t", ["x$(i)" for i in 1:9])
    MOI.set(
        bridged_mock,
        MOI.VariableName(),
        MOI.get(bridged_mock, MOI.ListOfVariableIndices()),
        var_names,
    )
    nonneg = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{MOI.VectorOfVariables,MOI.Nonnegatives}(),
    )
    relentr = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.RelativeEntropyCone,
        }(),
    )
    aux = MOI.get(mock, MOI.ListOfVariableIndices())[end]
    MOI.set(mock, MOI.VariableName(), aux, "aux")
    @test length(nonneg) == 1
    MOI.set(mock, MOI.ConstraintName(), nonneg[1], "nonneg")
    @test length(relentr) == 1
    MOI.set(mock, MOI.ConstraintName(), relentr[1], "relentr")
    equalto = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{Float64},
            MOI.EqualTo{Float64},
        }(),
    )
    @test length(equalto) == 9
    equalto_names = ["equalto$(i)" for i in 1:9]
    MOI.set.(bridged_mock, MOI.ConstraintName(), equalto, equalto_names)

    s = """
    variables: t, x1, x2, x3, x4, x5, x6, x7, x8, x9, aux
    equalto1: 1.0x1 in MathOptInterface.EqualTo(1.0)
    equalto2: 1.0x2 in MathOptInterface.EqualTo(1.0)
    equalto3: 1.0x3 in MathOptInterface.EqualTo(1.0)
    equalto4: 1.0x4 in MathOptInterface.EqualTo(1.0)
    equalto5: 1.0x5 in MathOptInterface.EqualTo(1.0)
    equalto6: 1.0x6 in MathOptInterface.EqualTo(1.0)
    equalto7: 1.0x7 in MathOptInterface.EqualTo(1.0)
    equalto8: 1.0x8 in MathOptInterface.EqualTo(1.0)
    equalto9: 1.0x9 in MathOptInterface.EqualTo(1.0)
    nonneg: [aux] in MathOptInterface.Nonnegatives(1)
    relentr: [0.0, x1, x2, x3, x4, x5, x6, x7, x8, x9, t + aux, t + aux, t + aux, t + aux, t + aux, t + aux, t + aux, t + aux, t + aux] in MathOptInterface.RelativeEntropyCone(19)
    maxobjective: t
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        mock,
        model,
        vcat(var_names, "aux"),
        vcat(equalto_names, "nonneg", "relentr"),
    )
    geomean = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.GeometricMeanCone,
        }(),
    )
    @test length(geomean) == 1
    MOI.set(bridged_mock, MOI.ConstraintName(), geomean[1], "geomean")
    equalto = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{Float64},
            MOI.EqualTo{Float64},
        }(),
    )
    @test length(equalto) == 9
    equalto_names = ["equalto$(i)" for i in 1:9]
    MOI.set.(bridged_mock, MOI.ConstraintName(), equalto, equalto_names)

    s = """
    variables: t, x1, x2, x3, x4, x5, x6, x7, x8, x9
    equalto1: 1.0x1 in MathOptInterface.EqualTo(1.0)
    equalto2: 1.0x2 in MathOptInterface.EqualTo(1.0)
    equalto3: 1.0x3 in MathOptInterface.EqualTo(1.0)
    equalto4: 1.0x4 in MathOptInterface.EqualTo(1.0)
    equalto5: 1.0x5 in MathOptInterface.EqualTo(1.0)
    equalto6: 1.0x6 in MathOptInterface.EqualTo(1.0)
    equalto7: 1.0x7 in MathOptInterface.EqualTo(1.0)
    equalto8: 1.0x8 in MathOptInterface.EqualTo(1.0)
    equalto9: 1.0x9 in MathOptInterface.EqualTo(1.0)
    geomean: [1.0t, x1, x2, x3, x4, x5, x6, x7, x8, x9] in MathOptInterface.GeometricMeanCone(10)
    maxobjective: t
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        bridged_mock,
        model,
        var_names,
        vcat(equalto_names, "geomean"),
    )
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.GeometricMeanCone,
            }(),
        ),
    )
    for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        value =
            (attr isa MOI.ConstraintPrimalStart) ? ones(10) :
            vcat(-1, fill(inv(9), 9))
        MOI.set(bridged_mock, attr, ci, value)
        @test MOI.get(bridged_mock, attr, ci) ≈ value
        if attr isa MOI.ConstraintPrimalStart
            @test MOI.get(mock, MOI.VariablePrimalStart(), aux) == 0
            @test MOI.get(mock, attr, nonneg[1]) == [0]
            @test MOI.get(mock, attr, relentr[1]) == vcat(0, ones(18))
        else
            @test MOI.get(mock, attr, nonneg[1]) == [1]
            @test MOI.get(mock, attr, relentr[1]) ≈ relentr_dual
        end
    end
    _test_delete_bridge(
        bridged_mock,
        ci,
        10,
        (
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}, 9),
            (MOI.VectorOfVariables, MOI.Nonnegatives, 0),
            (MOI.VectorAffineFunction{Float64}, MOI.RelativeEntropyCone, 0),
        ),
    )
    return
end

function test_conic_GeometricMeanCone_VectorAffineFunction_3()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.GeoMeantoRelEntr{Float64}(mock)
    var_primal = [2, 2, 0]
    relentr_dual = [2, 2, -2]
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            var_primal,
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [-2],
            (MOI.VectorOfVariables, MOI.Nonnegatives) => [[2]],
            (MOI.VectorAffineFunction{Float64}, MOI.RelativeEntropyCone) =>
                [relentr_dual],
        )

    MOI.Test.test_conic_GeometricMeanCone_VectorAffineFunction_3(
        bridged_mock,
        config,
    )
    var_names = ["t", "x"]
    MOI.set(
        bridged_mock,
        MOI.VariableName(),
        MOI.get(bridged_mock, MOI.ListOfVariableIndices()),
        var_names,
    )
    nonneg = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{MOI.VectorOfVariables,MOI.Nonnegatives}(),
    )
    relentr = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.RelativeEntropyCone,
        }(),
    )
    aux = MOI.get(mock, MOI.ListOfVariableIndices())[end]
    MOI.set(mock, MOI.VariableName(), aux, "aux")
    @test length(nonneg) == 1
    MOI.set(mock, MOI.ConstraintName(), nonneg[1], "nonneg")
    @test length(relentr) == 1
    MOI.set(mock, MOI.ConstraintName(), relentr[1], "relentr")
    less = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{Float64},
            MOI.LessThan{Float64},
        }(),
    )
    @test length(less) == 1
    MOI.set(mock, MOI.ConstraintName(), less[1], "less")
    s = """
    variables: t, x, aux
    less: 1.0x in MathOptInterface.LessThan(2.0)
    nonneg: [aux] in MathOptInterface.Nonnegatives(1)
    relentr: [0.0, x, t + aux] in MathOptInterface.RelativeEntropyCone(3)
    maxobjective: 2.0t
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        mock,
        model,
        vcat(var_names, "aux"),
        ["less", "nonneg", "relentr"],
    )
    geomean = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.GeometricMeanCone,
        }(),
    )
    @test length(geomean) == 1
    MOI.set(bridged_mock, MOI.ConstraintName(), geomean[1], "geomean")
    less = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{Float64},
            MOI.LessThan{Float64},
        }(),
    )
    @test length(less) == 1
    MOI.set(bridged_mock, MOI.ConstraintName(), less[1], "less")

    s = """
    variables: t, x
    less: 1.0x in MathOptInterface.LessThan(2.0)
    geomean: [1.0t, x] in MathOptInterface.GeometricMeanCone(2)
    maxobjective: 2.0t
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        bridged_mock,
        model,
        var_names,
        ["less", "geomean"],
    )
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.GeometricMeanCone,
            }(),
        ),
    )
    for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        value = (attr isa MOI.ConstraintPrimalStart) ? [2, 2] : [-2, 2]
        MOI.set(bridged_mock, attr, ci, value)
        @test MOI.get(bridged_mock, attr, ci) ≈ value
        if attr isa MOI.ConstraintPrimalStart
            @test MOI.get(mock, MOI.VariablePrimalStart(), aux) == 0
            @test MOI.get(mock, attr, nonneg[1]) == [0]
            @test MOI.get(mock, attr, relentr[1]) == [0, 2, 2]
        else
            @test MOI.get(mock, attr, nonneg[1]) == [2]
            @test MOI.get(mock, attr, relentr[1]) ≈ relentr_dual
        end
    end
    _test_delete_bridge(
        bridged_mock,
        ci,
        2,
        (
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}, 1),
            (MOI.VectorOfVariables, MOI.Nonnegatives, 0),
            (MOI.VectorAffineFunction{Float64}, MOI.RelativeEntropyCone, 0),
        ),
    )
    return
end

end  # module

TestConstraintGeomeanToRelentr.runtests()

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
# ---
# title: Phase Quantization
# id: demo_quantizephase
# cover: assets/quantizephase.gif
# author: Peter Kovesi
# date: 2018-10-26
# ---

# Phase values in an image are important.  However, despite this, phase can be quantized
# very heavily with little perceptual loss.  It can be quantized to a few as four levels, or
# even three.  Quantizing to two levels still gives an image that can be interpreted.

using TestImages
using Images
using ImagePhaseCongruency

img = Float64.(restrict(testimage("mandril_gray")))

results = map((8, 4, 3, 2)) do nlevels
    out = quantizephase(img, nlevels)
    clamp01!(Gray.(out))
end
mosaic(results; nrow=1)

# save cover image #src
isdir("assets") || mkdir("assets") #src
save(joinpath("assets", "quantizephase.gif"), Images.gif([results...]); fps=1) #src

max_stars_repo_name: peterkovesi/ImagePhaseCongruency

max_issues_repo_name: peterkovesi/ImagePhaseCongruency

max_forks_repo_name: peterkovesi/ImagePhaseCongruency


text:
@testset "199.binary-tree-right-side-view.jl" begin
    root = TreeNode(
        1,
        TreeNode(2,
            nothing,
            TreeNode(5)
            ),
        TreeNode(3,
            nothing,
            TreeNode(4)
            ),
        )
    @test right_side_view(root) == [1, 3, 4]
end


text:
include("modelrunner.jl")
using Hwloc
Hwloc.num_physical_cores()

### https://github.com/crstnbr/JuliaWorkshop19/blob/master/3_Three/1_parallel_computing.ipynb

function batch_runner()
    batches = 10

end

using Distributed

using Distributed, BenchmarkTools; rmprocs(workers()); addprocs(4); nworkers()

# serial version - count heads in a series of coin tosses
function add_serial(n)
    c = 0
    for i = 1:n
        c += rand(Bool)
    end
    c
end

@btime add_serial(200_000_000);

function add_distributed(n)
    c = @distributed (+) for i in 1:n
        Int(rand(Bool))
    end
    c
end

@btime add_distributed(200_000_000);

function add_distributed(n)
    c = @distributed (+) for i in 1:n
        x = Int(rand(Bool))
        println(x);
        x
    end
    c
end

add_distributed(8);


text:
module DFMS
using DataFrames
using LsqFit
using PyPlot
using HDF5
using Spice

include("typedefs.jl")
include("io.jl")
include("spice_funcs.jl")
include("fitting.jl")
include("plot_funcs.jl")
include("data_transform.jl")

export crop,
       clean,
       clean_h5,
       findPeaks,
       get_dfms_data,
       getGainFactor,
       get_gain_factor,
       getPixelGain,
       get_pixel_gain,
       get_eph,
       get_ephemeris_spice,
       get_outliers,
       load_h5_data,
       load_h5_data!,
       baseline_model,
       moving_avg,
       parseDataFileBothRows,
       parseDataFile,
       pds2h5,
       peakFit,
       polyFit,
       remove_outliers,
       remove_outliers!,
       show_me,
       show_time_series,
       show_single_spec,
       single_spec,
       time_series,
       time_series_autospec,
       update!,
       findSimilar,
       doCO2Correction

c_adc = 2.5 / (2.^12 - 1)
c_leda = 4.22 * 10.0^-12
q = 1.60217657 * 10.0^-19
const global t_integration = 19.8
const global C = c_adc * c_leda / q

const global _rowA = 1
const global _rowB = 2

#=
function findSimilar(t1, n1, t2, n2)
  const dtMax = DateTime(2000,1,1,0,35) - DateTime(2000,1,1,0)
  t_sync = DateTime[]
  n1_sync = Float64[]
  n2_sync = Float64[]

  for i=1:length(t2)
    for j=1:length(t1)
      if ((t2[i] > t1[j]) & ((t2[i] - t1[j]) < dtMax))
        push!(n1_sync, n1[j])
        push!(n2_sync, n2[i])
        push!(t_sync, t2[i])
        break
      end
    end
  end
  return t_sync, n1_sync, n2_sync

end

function findSimilar(df1::DataFrame, df2::DataFrame)
  const dtMax = DateTime(2000,1,1,0,35) - DateTime(2000,1,1,0)
  t_sync = DateTime[]
  n1_sync = Float64[]
  n2_sync = Float64[]
  dxPeaks_sync = Float64[]
  peakCenter1_sync = Float64[]
  peakCenter2_sync = Float64[]

  df1[:date] = DateTime[DateTime(tStr) for tStr in df1[:date]]
  df2[:date] = DateTime[DateTime(tStr) for tStr in df2[:date]]

  t1 = convert(Array, df1[:date])
  t2 = convert(Array, df2[:date])
  n1 = convert(Array, df1[:peakArea])
  n2 = convert(Array, df2[:peakArea])
  peakCenter1 = convert(Array, df1[:peakCenter])
  peakCenter2 = convert(Array, df2[:peakCenter])

  for i=1:length(t2)
    for j=1:length(t1)
      if ((t2[i] > t1[j]) & ((t2[i] - t1[j]) < dtMax))
        push!(n1_sync, n1[j])
        push!(n2_sync, n2[i])
        push!(t_sync, t2[i])
        push!(dxPeaks_sync, (peakCenter1[j] - peakCenter2[i]))
        push!(peakCenter1_sync, peakCenter1[j])
        push!(peakCenter2_sync, peakCenter2[i])
        break
      end
    end
  end
  return t_sync, n1_sync, n2_sync, dxPeaks_sync, peakCenter1_sync, peakCenter2_sync

end
=#

end


text:
#unum.jl - a julia implementation of the unum
# this file is the module definition file and also contains
# includes for all of the components which make it work

#for now, only compatible with 64-bit architectures.
@assert(sizeof(Int) == 8, "currently only compatible with 64-bit architectures")

module Unums

#engage the unpms options engine.
include("./options/options.jl")
#various tools to help coding.
include("./tools.jl")

################################################################################
#TYPE DEFINITION FILES
#type definitions for int64 array.
include("./int64op/i64o-typedefs.jl")
#type definition of unum.
include("./unum/unum-typedefs.jl")
#type definition of ubound
include("./ubound/ubound-typedefs.jl")
#type definition of gnum.
include("./gnum/gnum-typedefs.jl")

################################################################################
#IMPLEMENTATION FILES
#implementation of int64 and int64 array utility code.
include("./int64op/int64ops.jl")
#implementation of unums.
include("./unum/unum.jl")
#ubound-related code
include("./ubound/ubound.jl")
#utype - the type wrapper that is better than Union{Unum, Ubound}
include("./utype/utype.jl")
#gnums-related code
#include("./gnum/gnum.jl")

#we may need to patch julia in certain places to make this work smoothly.
include("./julia-patches.jl")

#utility files
include("./ubox/ubox.jl")
include("./ubox/uslice.jl")
#include("unum-bitwalk.jl")
#include("unum-promote.jl")
#include("unum-expwalk.jl")
#include("unum_solver.jl")

end #module

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
function Relativity(arg0::jdouble)
    return Relativity((jdouble,), arg0)
end

function acceleration(obj::Relativity, arg0::FieldSpacecraftState, arg1::Vector{RealFieldElement})
    return jcall(obj, "acceleration", FieldVector3D, (FieldSpacecraftState, Vector{RealFieldElement}), arg0, arg1)
end

function acceleration(obj::Relativity, arg0::SpacecraftState, arg1::Vector{jdouble})
    return jcall(obj, "acceleration", Vector3D, (SpacecraftState, Vector{jdouble}), arg0, arg1)
end

function depends_on_position_only(obj::Relativity)
    return jcall(obj, "dependsOnPositionOnly", jboolean, ())
end

function get_events_detectors(obj::Relativity)
    return jcall(obj, "getEventsDetectors", Stream, ())
end

function get_field_events_detectors(obj::Relativity, arg0::Field)
    return jcall(obj, "getFieldEventsDetectors", Stream, (Field,), arg0)
end

function get_parameters_drivers(obj::Relativity)
    return jcall(obj, "getParametersDrivers", Vector{ParameterDriver}, ())
end


text:
export convFFTKernel, getConvFFTKernel
## For the functions nImgIn, nImgOut, nFeatIn, nFeatOut, nTheta, getOp, initTheta : see AbstractConvKernel.jl
## All convKernel types are assumed to have fields nImage and sK
mutable struct convFFTKernel{T} <: AbstractConvKernel{T}
    nImg :: Array{Int,1}
    sK   :: Array{Int,1}
	Kp   :: Array{T}
	I    :: Array{Int}
end

function getConvFFTKernel(TYPE::Type,nImg,sK)
	return convFFTKernel{TYPE}(nImg,sK,(TYPE)[],(Int64)[])
end

function getKp(this::convFFTKernel{T}) where {T<:Number}
	# setup the padded convolution kernel and get the indices of non-zeros
	if isempty(this.Kp)
		theta = reshape(T.(collect(1:prod(this.sK))),this.sK[1],this.sK[2],prod(this.sK[3:4]))
		Kp    = zeros(T,this.nImg[1],this.nImg[2],size(theta,3))
		Kp[1:this.sK[1],1:this.sK[2],:] = theta;
		center = (this.sK[1:2]+1)./2
		Kp    = circshift(Kp,1-center);
		I = find(Kp)
		this.Kp = zero(T)*Kp;
		idp     = sortperm(Kp[I])
		this.I  = I[idp];
	end
	return this.Kp,this.I
end
function getK1(this::convFFTKernel{T},theta::Array{T}) where {T<:Number}
	# get first columns of convolultion matrices
	Kp,I = getKp(this)
	for k=1:length(I)
		Kp[I[k]]  = theta[k]
	end
	return Kp
end

# methods for A*x
function multRed!(Zkh::Array{Complex{T},2},S::Array{Complex{T},4},Yh::Array{Complex{T},3},k::Int) where {T<:Number}
	# compute  Zkh[i1,i2,k] = S[i1,i2,:,k]'*Yh[i1,i2,:]
	Zkh[:]=Complex128(0.0)
	for i3=1:size(Yh,3)
		for i2=1:size(Zkh,2)
			for i1=1:size(Zkh,1)
				@inbounds Zkh[i1,i2] += S[i1,i2,i3,k].*Yh[i1,i2,i3]
			end
		end
	end
end

function Amv!(this::convFFTKernel{T},Z::AbstractArray{T,3},S::Array{Complex{T}},Y::AbstractArray{T,3},Yh::Array{Complex{T},3},Zkh::Array{Complex{T},2}) where {T<:Number}
	# 2D convolution for a single example.
	Yh[:]=Y; ifft2!(Yh)
    for k=1:this.sK[4]
		multRed!(Zkh,S,Yh,k)
        Z[:,:,k]  = real(fft2!(Zkh))
    end
	return Z
end
function Amv(this::convFFTKernel{T},theta::Array{T},Y::Array{T}) where {T<:Number}
		nex   = div(length(Y),prod(nImgIn(this)))
		Z = zeros(T,tuple([nImgOut(this); nex]...))
	    Amv!(this,Z,theta,Y)
		return Z
end
function Amv!(this::convFFTKernel{T},Z::Array{T},theta::Array{T},Y::Array{T}) where {T<:Number}
    nex   = div(length(Y),prod(nImgIn(this)))
	Y     = reshape(Y,tuple([nImgIn(this);nex]...))
	# pre-allocation for temps
	Ykh = zeros(Complex{T},this.nImg[1],this.nImg[2],this.sK[3])
	Zk  = zeros(T, this.nImg[1],this.nImg[2], this.sK[4])
    Zik = zeros(Complex{T}, this.nImg[1],this.nImg[2])

	# get kernel
	S = reshape( fft2(getK1(this,theta)), this.nImg[1],this.nImg[2], this.sK[3], this.sK[4])

    # compute convolution
    for k=1:nex
		Amv!(this,view(Z,:,:,:,k),S,view(Y,:,:,:,k),Ykh,Zik)
	end
	Z = reshape(Z,:,nex)
	return Z
end

# methods for A'*x
function multRedT!(Ykh::Array{Complex{T},2},S::Array{Complex{T},4},Zh::Array{Complex{T},3},j::Int) where {T<:Number}
	# compute  Ykh[i1,i2,k] = S[i1,i2,j,:]'*Zh[i1,i2,:]
	Ykh[:]=Complex{T}(0.0)
	for i3=1:size(Zh,3)
		for i2=1:size(Ykh,2)
			for i1=1:size(Ykh,1)
				@inbounds Ykh[i1,i2] += S[i1,i2,j,i3].*Zh[i1,i2,i3]
			end
		end
	end
	return Ykh
end

function ATmv!(this::convFFTKernel{T},Y::AbstractArray{T,3},S::Array{Complex{T}},Z::AbstractArray{T,3},Zh::Array{Complex{T},3},Ykh::Array{Complex{T},2}) where {T<:Number}
	# 2D convolution for a single example.
	Zh[:]=Z; Zh =fft2!(Zh)
    for j=1:this.sK[3]
		multRedT!(Ykh,S,Zh,j)
        Y[:,:,j]  = real(ifft2!(Ykh))
    end
end
function ATmv(this::convFFTKernel{T},theta::Array{T},Z::Array{T}) where {T<:Number}
		nex = div(length(Z),prod(nImgOut(this)))
		Y   = zeros(T,tuple([nImgIn(this); nex]...))
	    ATmv!(this,Y,theta,Z)
		return Y
end

function ATmv!(this::convFFTKernel{T},Y::Array{T},theta::Array{T},Z::Array{T}) where {T<:Number}
    nex   = div(length(Y),prod(nImgIn(this)))
	Z     = reshape(Z,tuple([nImgOut(this);nex]...))

	# pre-allocation for temps
	Zkh = zeros(Complex{T},this.nImg[1],this.nImg[2],this.sK[4])
	Yk  = zeros(T, this.nImg[1],this.nImg[2], this.sK[3])
    Yik = zeros(Complex{T}, this.nImg[1],this.nImg[2])

	# get kernel
	S = reshape( fft2(getK1(this,theta)), this.nImg[1],this.nImg[2], this.sK[3], this.sK[4])

    # compute convolution
    for k=1:nex
		ATmv!(this,view(Y,:,:,:,k),S,view(Z,:,:,:,k),Zkh,Yik)
	end
	Y = reshape(Y,:,nex)
	return Y
end

function Jthetamv(this::convFFTKernel{T},dtheta::Array{T},dummy::Array{T},Y::Array{T},temp=nothing) where {T<:Number}
    return Amv(this,dtheta,Y)
end

function JthetaTmv(this::convFFTKernel{T},Z::Array{T},dummy::Array{T},Y::Array{T}) where {T<:Number}
	nex    =  div(length(Y),nFeatIn(this))

	Y    = reshape(Y,this.nImg[1],this.nImg[2], this.sK[3],nex)
	Z    = reshape(Z,this.nImg[1],this.nImg[2], this.sK[4],nex)

	# temps
	Yh   = zeros(Complex{T},this.nImg[1],this.nImg[2],this.sK[3])
	Zh   = zeros(Complex{T},this.nImg[1],this.nImg[2],this.sK[4])
	tt   = zeros(Complex{T},this.nImg[1],this.nImg[2],this.sK[3],this.sK[4])

	for k=1:nex
		Yh[:] = Y[:,:,:,k]
		Zh[:] = Z[:,:,:,k]

		ifft2!(Yh)
		fft2!(Zh)

		for i4=1:this.sK[4]
			for i3=1:this.sK[3]
				for i2=1:this.nImg[2]
					for i1=1:this.nImg[1]
						@inbounds tt[i1,i2,i3,i4] += Yh[i1,i2,i3].*Zh[i1,i2,i4]
					end
				end
			end
		end
	end
	tt = real(fft2(reshape(tt,this.nImg[1],this.nImg[2],:)))
	return tt[this.I]
end


text:
isdefined(Base, :__precompile__) && __precompile__(false)

module SMC
    using BenchmarkTools, DataFrames, Dates, Distributed, Distributions
    using FileIO, HDF5, JLD2, LinearAlgebra, Random, Test
    using ModelConstructors, SparseArrays

    using Roots: fzero, ConvergenceFailed
    using StatsBase: sample, Weights

    import Base.<, Base.isempty, Base.min, Base.max
    import Calculus, ModelConstructors
    import SparseArrays.SparseMatrixCSC

    export
        compute_parameter_covariance, get_estimation_output_files,
        compute_moments, find_density_bands, mutation, resample, smc,
        mvnormal_mixture_draw, nearest_spd, marginal_data_density,
        initial_draw!, Cloud, get_cloud

    const VERBOSITY   = Dict(:none => 0, :low => 1, :high => 2)
    const DATE_FORMAT = "yymmdd"

    include("particle.jl")
    include("initialization.jl")
    include("helpers.jl")
    include("util.jl")
    include("mutation.jl")
    include("resample.jl")
    include("smc_main.jl")
end


text:
module Dionysus
using Flux, FileIO, Pages, Mux

# Build a RNN

rnn = Flux.RNNCell(10, 5)

x = rand(10) # dummy data
h = rand(5)  # initial hidden state

# h, y = rnn2(h, x)

m = Flux.Recur(rnn, h)

seq = [BigFloat(π^i) for i = 1:200]
m.(seq)
seq2 = [BigFloat(e^i) for i = 1:200]

y = m(x)
end # module
msg = hex2bytes("Hello, World!")

max_stars_repo_name: jempyre/Dionysus

max_issues_repo_name: jempyre/Dionysus

max_forks_repo_name: jempyre/Dionysus


text:
using Rainflow
using Base.Test

srand(1234)
signal = randn(100)
ext, t = sort_peaks(signal)
result = count_cycles(ext,t)
bins = sum_cycles(result,13,1)
bins1 = sum_cycles(result, collect(linspace(0, 100, 13+1)), collect(linspace(0, 100, 1+1)))
bins2 = sum_cycles(result, linspace(0, 100, 13+1), linspace(0, 100, 1+1))
@test bins==bins1==bins2
@test squeeze(bins,2) == [4.0, 6.0, 7.5, 3.0, 2.5, 5.0, 0.0, 2.5, 1.0, 0.0, 0.5, 0.0, 1.0]


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ServiceSasParameters <: SwaggerModel
    canonicalizedResource::Any # spec type: Union{ Nothing, String } # spec name: canonicalizedResource
    signedResource::Any # spec type: Union{ Nothing, String } # spec name: signedResource
    signedPermission::Any # spec type: Union{ Nothing, String } # spec name: signedPermission
    signedIp::Any # spec type: Union{ Nothing, String } # spec name: signedIp
    signedProtocol::Any # spec type: Union{ Nothing, String } # spec name: signedProtocol
    signedStart::Any # spec type: Union{ Nothing, DateTime } # spec name: signedStart
    signedExpiry::Any # spec type: Union{ Nothing, DateTime } # spec name: signedExpiry
    signedIdentifier::Any # spec type: Union{ Nothing, String } # spec name: signedIdentifier
    startPk::Any # spec type: Union{ Nothing, String } # spec name: startPk
    endPk::Any # spec type: Union{ Nothing, String } # spec name: endPk
    startRk::Any # spec type: Union{ Nothing, String } # spec name: startRk
    endRk::Any # spec type: Union{ Nothing, String } # spec name: endRk
    keyToSign::Any # spec type: Union{ Nothing, String } # spec name: keyToSign
    rscc::Any # spec type: Union{ Nothing, String } # spec name: rscc
    rscd::Any # spec type: Union{ Nothing, String } # spec name: rscd
    rsce::Any # spec type: Union{ Nothing, String } # spec name: rsce
    rscl::Any # spec type: Union{ Nothing, String } # spec name: rscl
    rsct::Any # spec type: Union{ Nothing, String } # spec name: rsct

    function ServiceSasParameters(;canonicalizedResource=nothing, signedResource=nothing, signedPermission=nothing, signedIp=nothing, signedProtocol=nothing, signedStart=nothing, signedExpiry=nothing, signedIdentifier=nothing, startPk=nothing, endPk=nothing, startRk=nothing, endRk=nothing, keyToSign=nothing, rscc=nothing, rscd=nothing, rsce=nothing, rscl=nothing, rsct=nothing)
        o = new()
        validate_property(ServiceSasParameters, Symbol("canonicalizedResource"), canonicalizedResource)
        setfield!(o, Symbol("canonicalizedResource"), canonicalizedResource)
        validate_property(ServiceSasParameters, Symbol("signedResource"), signedResource)
        setfield!(o, Symbol("signedResource"), signedResource)
        validate_property(ServiceSasParameters, Symbol("signedPermission"), signedPermission)
        setfield!(o, Symbol("signedPermission"), signedPermission)
        validate_property(ServiceSasParameters, Symbol("signedIp"), signedIp)
        setfield!(o, Symbol("signedIp"), signedIp)
        validate_property(ServiceSasParameters, Symbol("signedProtocol"), signedProtocol)
        setfield!(o, Symbol("signedProtocol"), signedProtocol)
        validate_property(ServiceSasParameters, Symbol("signedStart"), signedStart)
        setfield!(o, Symbol("signedStart"), signedStart)
        validate_property(ServiceSasParameters, Symbol("signedExpiry"), signedExpiry)
        setfield!(o, Symbol("signedExpiry"), signedExpiry)
        validate_property(ServiceSasParameters, Symbol("signedIdentifier"), signedIdentifier)
        setfield!(o, Symbol("signedIdentifier"), signedIdentifier)
        validate_property(ServiceSasParameters, Symbol("startPk"), startPk)
        setfield!(o, Symbol("startPk"), startPk)
        validate_property(ServiceSasParameters, Symbol("endPk"), endPk)
        setfield!(o, Symbol("endPk"), endPk)
        validate_property(ServiceSasParameters, Symbol("startRk"), startRk)
        setfield!(o, Symbol("startRk"), startRk)
        validate_property(ServiceSasParameters, Symbol("endRk"), endRk)
        setfield!(o, Symbol("endRk"), endRk)
        validate_property(ServiceSasParameters, Symbol("keyToSign"), keyToSign)
        setfield!(o, Symbol("keyToSign"), keyToSign)
        validate_property(ServiceSasParameters, Symbol("rscc"), rscc)
        setfield!(o, Symbol("rscc"), rscc)
        validate_property(ServiceSasParameters, Symbol("rscd"), rscd)
        setfield!(o, Symbol("rscd"), rscd)
        validate_property(ServiceSasParameters, Symbol("rsce"), rsce)
        setfield!(o, Symbol("rsce"), rsce)
        validate_property(ServiceSasParameters, Symbol("rscl"), rscl)
        setfield!(o, Symbol("rscl"), rscl)
        validate_property(ServiceSasParameters, Symbol("rsct"), rsct)
        setfield!(o, Symbol("rsct"), rsct)
        o
    end
end # type ServiceSasParameters

const _property_map_ServiceSasParameters = Dict{Symbol,Symbol}(Symbol("canonicalizedResource")=>Symbol("canonicalizedResource"), Symbol("signedResource")=>Symbol("signedResource"), Symbol("signedPermission")=>Symbol("signedPermission"), Symbol("signedIp")=>Symbol("signedIp"), Symbol("signedProtocol")=>Symbol("signedProtocol"), Symbol("signedStart")=>Symbol("signedStart"), Symbol("signedExpiry")=>Symbol("signedExpiry"), Symbol("signedIdentifier")=>Symbol("signedIdentifier"), Symbol("startPk")=>Symbol("startPk"), Symbol("endPk")=>Symbol("endPk"), Symbol("startRk")=>Symbol("startRk"), Symbol("endRk")=>Symbol("endRk"), Symbol("keyToSign")=>Symbol("keyToSign"), Symbol("rscc")=>Symbol("rscc"), Symbol("rscd")=>Symbol("rscd"), Symbol("rsce")=>Symbol("rsce"), Symbol("rscl")=>Symbol("rscl"), Symbol("rsct")=>Symbol("rsct"))
const _property_types_ServiceSasParameters = Dict{Symbol,String}(Symbol("canonicalizedResource")=>"String", Symbol("signedResource")=>"String", Symbol("signedPermission")=>"String", Symbol("signedIp")=>"String", Symbol("signedProtocol")=>"String", Symbol("signedStart")=>"DateTime", Symbol("signedExpiry")=>"DateTime", Symbol("signedIdentifier")=>"String", Symbol("startPk")=>"String", Symbol("endPk")=>"String", Symbol("startRk")=>"String", Symbol("endRk")=>"String", Symbol("keyToSign")=>"String", Symbol("rscc")=>"String", Symbol("rscd")=>"String", Symbol("rsce")=>"String", Symbol("rscl")=>"String", Symbol("rsct")=>"String")
Base.propertynames(::Type{ ServiceSasParameters }) = collect(keys(_property_map_ServiceSasParameters))
Swagger.property_type(::Type{ ServiceSasParameters }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ServiceSasParameters[name]))}
Swagger.field_name(::Type{ ServiceSasParameters }, property_name::Symbol) =  _property_map_ServiceSasParameters[property_name]

const _allowed_ServiceSasParameters_signedResource = ["b", "c", "f", "s"]

const _allowed_ServiceSasParameters_signedPermission = ["r", "d", "w", "l", "a", "c", "u", "p"]

const _allowed_ServiceSasParameters_signedProtocol = ["https,http", "https"]

function check_required(o::ServiceSasParameters)
    (getproperty(o, Symbol("canonicalizedResource")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ServiceSasParameters }, name::Symbol, val)
    if name === Symbol("signedResource")
        Swagger.validate_param(name, "ServiceSasParameters", :enum, val, _allowed_ServiceSasParameters_signedResource)
    end
    if name === Symbol("signedPermission")
        Swagger.validate_param(name, "ServiceSasParameters", :enum, val, _allowed_ServiceSasParameters_signedPermission)
    end
    if name === Symbol("signedProtocol")
        Swagger.validate_param(name, "ServiceSasParameters", :enum, val, _allowed_ServiceSasParameters_signedProtocol)
    end
    if name === Symbol("signedIdentifier")
        Swagger.validate_param(name, "ServiceSasParameters", :maxLength, val, 64)
    end
end


text:
module VectorSpaces

using ..Combinatorics: AbstractCombinatorics
using ...Prerequisites: rawtype
using ...Prerequisites.TypeTraits: efficientoperations, forder, corder, subtoind, indtosub

import ...Interfaces: dimension, ⊕, rank, dims, inds

export VectorSpace
export HasTable, TableSorted, IsMultiIndexable, MultiIndexOrderStyle
export SimpleVectorSpace
export OrderedIndices, SimpleIndices, TabledIndices
export NamedVectorSpace

"""
    VectorSpace{B} <: AbstractVector{B}

Abstract vector space.
"""
abstract type VectorSpace{B} <: AbstractVector{B} end
Base.:(==)(vs1::VectorSpace, vs2::VectorSpace) = ==(efficientoperations, vs1, vs2)
Base.isequal(vs1::VectorSpace, vs2::VectorSpace) = isequal(efficientoperations, vs1, vs2)
Base.size(vs::VectorSpace) = (dimension(vs),)
Base.iterate(vs::VectorSpace, state::Integer=1) = (state > length(vs)) ? nothing : (vs[state], state+1)

"""
    HasTable(B::Bool)
    HasTable(::Type{<:VectorSpace})

Trait of whether a subtype of `VectorSpace` has the attribute `:table`.

Only two instances are allowed, the first of which is the default for a subtype:
* `HasTable(false)`: indication of not having the attribute `:table`
* `HasTable(true)`: indication of having the attribute `:table`
"""
struct HasTable{B} end
HasTable(B::Bool) = HasTable{B}()
HasTable(::Type{<:VectorSpace}) = HasTable{false}()

"""
    TableSorted(B::Bool)
    TableSorted(::Type{<:VectorSpace})

Trait of whether the attribute `:table` of a subtype of `VectorSpace` is sorted.

Only two instances are allowed, the first of which is the default for a subtype:
* `TableSorted(false)`: indication of unsorted attribute `:table`
* `TableSorted(true)`: indication of sorted attribute `:table`
"""
struct TableSorted{B} end
TableSorted(B::Bool) = TableSorted{B}()
TableSorted(::Type{<:VectorSpace}) = TableSorted{false}()

"""
    IsMultiIndexable(B::Bool)
    IsMultiIndexable(::Type{<:VectorSpace})

Trait of whether the bases of a subtype of `VectorSpace` can be represented by multiindices (Cartesian indices).

Only two instances are allowed, the first of which is the default for a subtype:
* `IsMultiIndexable(false)`: indication of irrepresentability by Cartesian indices
* `IsMultiIndexable(true)`: indication of representability by Cartesian indices
"""
struct IsMultiIndexable{B} end
IsMultiIndexable(B::Bool) = IsMultiIndexable{B}()
IsMultiIndexable(::Type{<:VectorSpace}) = IsMultiIndexable{false}()

"""
    MultiIndexOrderStyle(M::Char)
    MultiIndexOrderStyle(::Type{<:VectorSpace})

Trait of the order style of the Cartesian-index representation of the bases of a subtype of `VectorSpace`.

Only two instances are allowed, the first of which is the default for a subtype:
* `MultiIndexOrderStyle(`C`)`: indication of C/C++ order style
* `MultiIndexOrderStyle(`F`)`: indication of Fortran order style
"""
struct MultiIndexOrderStyle{M} end
MultiIndexOrderStyle(M::Char) = (M ∈ ('C', 'F')) ? MultiIndexOrderStyle{M}() : error("MultiIndexOrderStyle error: wrong input M($M).")
MultiIndexOrderStyle(::Type{<:VectorSpace}) = MultiIndexOrderStyle{'C'}()

"""
    dimension(vs::VectorSpace) -> Int

The dimension of a vector space.
"""
dimension(vs::VectorSpace) = _dimension_(vs, vs|>typeof|>HasTable)
_dimension_(vs::VectorSpace, ::HasTable{true}) = length(vs.table)
_dimension_(vs::VectorSpace, ::HasTable{false}) = _dimension_(vs, vs|>typeof|>IsMultiIndexable)
_dimension_(vs::VectorSpace, ::IsMultiIndexable{true}) = prod(dims(vs))

"""
    getindex(vs::VectorSpace, i::Int)

Get the ith basis of a vector space by the `[]` operator.
"""
Base.getindex(vs::VectorSpace, i::Int) = _getindex_(vs, i, vs|>typeof|>HasTable, vs|>typeof|>IsMultiIndexable)
_getindex_(vs::VectorSpace, i::Int, ::HasTable{true}, ::IsMultiIndexable{false}) = vs.table[i]
_getindex_(vs::VectorSpace, i::Int, ::HasTable{true}, ::IsMultiIndexable{true}) = rawtype(eltype(vs))(vs.table[i], vs)
_getindex_(vs::VectorSpace, i::Int, ::HasTable{false}, ::IsMultiIndexable{true}) = _getindex_(vs, i, vs|>typeof|>MultiIndexOrderStyle)
_getindex_(vs::VectorSpace, i::Int, ::MultiIndexOrderStyle{'F'}) = rawtype(eltype(vs))(indtosub(dims(vs), i, forder), vs)
_getindex_(vs::VectorSpace, i::Int, ::MultiIndexOrderStyle{'C'}) = rawtype(eltype(vs))(indtosub(dims(vs), i, corder), vs)

"""
    searchsortedfirst(vs::VectorSpace{B}, basis::B) where B -> Int
    searchsortedfirst(vs, basis) -> Int

Search the index of a basis in a vector space.
"""
Base.searchsortedfirst(vs::VectorSpace{B}, basis::B) where {B} = _search_(vs, basis, vs|>typeof|>HasTable)
_search_(vs::VectorSpace{B}, basis::B, ::HasTable{true}) where {B} = _search_(vs, basis, vs|>typeof|>TableSorted, vs|>typeof|>IsMultiIndexable)
_search_(vs::VectorSpace{B}, basis::B, ::TableSorted{true}, ::IsMultiIndexable{false}) where {B} = searchsortedfirst(vs.table, basis)
_search_(vs::VectorSpace{B}, basis::B, ::TableSorted{false}, ::IsMultiIndexable{false}) where {B} = findfirst(isequal(basis), vs.table)
_search_(vs::VectorSpace{B}, basis::B, ::TableSorted{true}, ::IsMultiIndexable{true}) where {B} = searchsortedfirst(vs.table, inds(basis, vs))
_search_(vs::VectorSpace{B}, basis::B, ::TableSorted{false}, ::IsMultiIndexable{true}) where {B} = findfirst(isequal(inds(basis, vs)), vs.table)
_search_(vs::VectorSpace{B}, basis::B, ::HasTable{false}) where {B} = _search_(vs, basis, vs|>typeof|>IsMultiIndexable)
_search_(vs::VectorSpace{B}, basis::B, ::IsMultiIndexable{true}) where {B} = _search_(vs, basis, vs|>typeof|>MultiIndexOrderStyle)
function _search_(vs::VectorSpace{B}, basis::B, ::MultiIndexOrderStyle{'F'}) where B
    (basis ∈ vs) ? subtoind(dims(vs), inds(basis, vs), forder) : error("searchsortedfirst error: input basis($basis) not found.")
end
function _search_(vs::VectorSpace{B}, basis::B, ::MultiIndexOrderStyle{'C'}) where B
    (basis ∈ vs) ? subtoind(dims(vs), inds(basis, vs), corder) : error("searchsortedfirst error: input basis($basis) not found.")
end
function _search_(vs::VectorSpace{B}, basis::B, ::IsMultiIndexable{false}) where B
    for i = 1:dimension(vs)
        (vs[i] == basis) && return i
    end
end
function Base.searchsortedfirst(vs, basis)
    lo, hi = 0, length(vs)+1
    @inbounds while lo < hi-1
        m = (lo+hi) >>> 1
        if vs[m] < basis
            lo = m
        else
            hi = m
        end
    end
    return hi
end

"""
    findfirst(basis::B, vs::VectorSpace{B}) where B -> Int
    findfirst(bases, vs::VectorSpace) -> NTuple{length(bases), Int}

Get the index of a basis or the indices of a couple of bases in a vector space.
"""
Base.findfirst(basis::B, vs::VectorSpace{B}) where {B} = _find_(basis, vs, vs|>typeof|> HasTable)
Base.findfirst(bases, vs::VectorSpace) = NTuple{length(bases), Int}(findfirst(basis, vs)::Int for basis in bases)
_find_(basis::B, vs::VectorSpace{B}, ::HasTable{true}) where {B} = _find_(basis, vs, vs|>typeof|>TableSorted)
_find_(basis::B, vs::VectorSpace{B}, ::TableSorted{true}) where {B} = _findwithmissing_(basis, vs)
_find_(basis::B, vs::VectorSpace{B}, ::TableSorted{false}) where {B} = _findwithoutmissing_(basis, vs)
_find_(basis::B, vs::VectorSpace{B}, ::HasTable{false}) where {B} = _find_(basis, vs, vs|>typeof|>IsMultiIndexable)
_find_(basis::B, vs::VectorSpace{B}, ::IsMultiIndexable{true}) where {B} = _findwithoutmissing_(basis, vs)
_find_(basis::B, vs::VectorSpace{B}, ::IsMultiIndexable{false}) where {B} = _findwithmissing_(basis, vs)
function _findwithmissing_(basis::B, vs::VectorSpace{B}) where B
    index = searchsortedfirst(vs, basis)
    ((0 < index <= dimension(vs)) && (vs[index] == basis)) ? index : nothing
end
_findwithoutmissing_(basis::B, vs::VectorSpace{B}) where {B} = searchsortedfirst(vs, basis)

"""
    in(basis::B, vs::VectorSpace{B}) where B -> Bool

Judge whether a basis belongs to a vector space.
"""
Base.in(basis::B, vs::VectorSpace{B}) where {B} = _in_(basis, vs, vs|>typeof|>HasTable, vs|>typeof|>IsMultiIndexable)
_in_(basis::B, vs::VectorSpace{B}, ::HasTable, ::IsMultiIndexable) where {B} = isa(findfirst(basis, vs), Integer)
_in_(basis::B, vs::VectorSpace{B}, ::HasTable{false}, ::IsMultiIndexable{true}) where {B} = _in_(basis, vs, Val(vs|>typeof|>rank))
@generated function _in_(basis::B, vs::VectorSpace{B}, ::Val{N}) where {B, N}
    (N == 0) && return :(true)
    prepare = :(INDS = inds(basis, vs); DIMS = dims(vs))
    expr = :(0 < INDS[1] <= DIMS[1])
    for i = 2:N
        expr = Expr(:&&, expr, :(0 < INDS[$i] <= DIMS[$i]))
    end
    return Expr(:block, prepare, expr)
end

"""
    SimpleVectorSpace{S}(table::NTuple{N, B}) where {S, B, N}
    SimpleVectorSpace{S}(table...) where S

Simplest vector space, whose bases are stored in the attribute `:table` as an ntuple.

The `:table` attribute can be sorted or unsorted, which is determined by the type parameter `S`, with `'T'` for sorted and `'F'` for unsorted.
"""
struct SimpleVectorSpace{S, B, N} <: VectorSpace{B}
    table::NTuple{N, B}
    function SimpleVectorSpace{S}(table::NTuple{N, B}) where {S, B, N}
        @assert S ∈ ('T', 'F') "SimpleVectorSpace error: wrong type parameter S($S)."
        new{S, B, N}(table)
    end
end
SimpleVectorSpace{S}(table...) where {S} = SimpleVectorSpace{S}(table)
HasTable(::Type{<:SimpleVectorSpace}) = HasTable(true)
TableSorted(::Type{<:SimpleVectorSpace{'T'}}) = TableSorted(true)
TableSorted(::Type{<:SimpleVectorSpace{'F'}}) = TableSorted(false)

"""
    ⊕(basis1::B, basis2::B) where B -> SimpleVectorSpace{'F', B, 2}
    ⊕(basis::B, vs::SimpleVectorSpace{<:Any, B}) where B -> SimpleVectorSpace{'F', B}
    ⊕(vs::SimpleVectorSpace{<:Any, B}, basis::B) where B -> SimpleVectorSpace{'F', B}
    ⊕(vs1::SimpleVectorSpace{<:Any, B}, vs2::SimpleVectorSpace{<:Any, B}) where B -> SimpleVectorSpace{'F', B}

Get the direct sum between bases or simple vector spaces.
"""
⊕(basis1::B, basis2::B) where {B} = SimpleVectorSpace{'F'}(basis1, basis2)
⊕(basis::B, vs::SimpleVectorSpace{<:Any, B}) where {B} = SimpleVectorSpace{'F'}(basis, vs.table...)
⊕(vs::SimpleVectorSpace{<:Any, B}, basis::B) where {B} = SimpleVectorSpace{'F'}(vs.table..., basis)
⊕(vs1::SimpleVectorSpace{<:Any, B}, vs2::SimpleVectorSpace{<:Any, B}) where {B} = SimpleVectorSpace{'F'}(vs1.table..., vs2.table...)

"""
    OrderedIndices{N} <: VectorSpace{NTuple{N, Int}}

The simplest abstract class of multiindexable vector spaces, whose bases are just tuples of integers.

This class of vector spaces must have the following attribute:
`dims::NTuple{N, Int}`: the dimesnions of the Cartesian indices along all axes
"""
abstract type OrderedIndices{N} <: VectorSpace{NTuple{N, Int}} end
IsMultiIndexable(::Type{<:OrderedIndices}) = IsMultiIndexable(true)
dims(vs::OrderedIndices) = vs.dims
inds(basis::NTuple{N, Int}, ::OrderedIndices{N}) where {N} = basis
rank(::Type{<:OrderedIndices{N}}) where {N} = N
Tuple(index::NTuple{N, Int}, ::OrderedIndices{N}) where {N} = index

"""
    SimpleIndices{M}(dims::NTuple{N, Int}) where {M, N}

Simple ordered Cartesian indices.

!!! note
    1) It can be C/C++ ordered or Fortran ordered depending on the first type parameter `M`, with `'C'` for the former and `'F'` the latter.
    2) For its bases (Cartesian indices), there is no restriction except that they should be in the proper range defined by its `dims`.
"""
struct SimpleIndices{M, N} <: OrderedIndices{N}
    dims::NTuple{N, Int}
    function SimpleIndices{M}(dims::NTuple{N, Int}) where {M, N}
        @assert M ∈ ('C', 'F') "SimpleIndices error: wrong type parameter M($M)."
        new{M, N}(dims)
    end
end
SimpleIndices{M}(dims::Int...) where {M} = SimpleIndices{M}(dims)
SimpleIndices{M, N}(dims::NTuple{N, Int}) where {M, N} = SimpleIndices{M}(dims)
MultiIndexOrderStyle(::Type{<:SimpleIndices{M}}) where {M} = MultiIndexOrderStyle(M)

"""
    TabledIndices{S}(dims::NTuple{N, Int}, table::Vector{NTuple{N, Int}}) where {S, N}
    TabledIndices{N}(::Type{M}, len::Int) where {N, M<:AbstractCombinatorics}

Tabled ordered Cartesian indices.

Compared to [`SimpleIndices`](@ref), the bases of this kind of vector spaces are stored in the attribute `:table`, which must be a vector of tuple of integers. The `:table` attribute can be sorted or unsorted, which is determined by the type parameter `S`, with `'T'` for sorted and `'F'` for unsorted. This type suits the situations when the Cartesian indices are restricted by extra conditions except that propoesed by the attribute `:dims`.
"""
struct TabledIndices{S, N} <: OrderedIndices{N}
    dims::NTuple{N, Int}
    table::Vector{NTuple{N, Int}}
    TabledIndices{S}(dims::NTuple{N, Int}, table::Vector{NTuple{N, Int}}) where {S, N} = new{S, N}(dims, table)
end
TabledIndices{N}(::Type{M}, len::Int) where {N, M<:AbstractCombinatorics} = TabledIndices{'T'}(ntuple(i->len, N), M{N}(1:len)|>collect)
HasTable(::Type{<:TabledIndices}) = HasTable(true)
TableSorted(::Type{<:TabledIndices{'T'}}) = TableSorted(true)
TableSorted(::Type{<:TabledIndices{'F'}}) = TableSorted(false)

"""
    NamedVectorSpace{M, NS, BS<:Tuple, VS<:Tuple{Vararg{AbstractVector}}} <: VectorSpace{NamedTuple{NS, BS}}

Abstract named vector space.

This is a wrapper of multiindexable vector spaces, each of whose indexable dimensions is associated with a name.

It has four type parameters:
* `M`: mode of the named vector space. It specifies how the indexable dimensions are combined to form the bases of the named vector space, and must take one of the following values:
  - `:zip`: elements from each indexable dimensions are zipped together to form the bases,
  - `:product`: elements from each indexable dimensions are direct producted together to form the bases.
For the `:zip` mode, all the indexable dimensions should have the same number of elements, and the number of formed bases is equal to this number; for the `:product` mode, there are no restriction on the numbers of the indexable dimensions, and the number of the final bases is equal to their product.
* `NS::Tuple{Vararg{Symbol}}`: the names of the indexable dimensions
* `BS<:Tuple`: the eltypes of the indexable dimensions
* `VS<:Tuple{Vararg{AbstractVector}}`: the contents of the indexable dimensions

The concrete types must have the following attribute:
* `:contents::VS`: storage of the contents of the indexable dimensions

By default, a named vector space uses C order for the indexable dimensions when the mode is `:product`. You can change it to F order for your own subtypes by defining the [`MultiIndexOrderStyle`](@ref) trait.
"""
abstract type NamedVectorSpace{M, NS, BS<:Tuple, VS<:Tuple{Vararg{AbstractVector}}} <: VectorSpace{NamedTuple{NS, BS}} end
IsMultiIndexable(::Type{<:NamedVectorSpace}) = IsMultiIndexable(true)
MultiIndexOrderStyle(::Type{<:NamedVectorSpace}) = MultiIndexOrderStyle('C')
rank(::Type{<:NamedVectorSpace{:zip}}) = 1
rank(::Type{<:NamedVectorSpace{:product, NS}}) where NS = length(NS)
dims(nvs::NamedVectorSpace{:zip}) = (length(getfield(nvs, :contents)[1]),)
@generated dims(nvs::NamedVectorSpace{:product}) = Expr(:tuple, [:(length(getfield(nvs, :contents)[$i])) for i = 1:rank(nvs)]...)
function inds(basis::NamedTuple{NS}, nvs::NamedVectorSpace{:zip, NS}) where NS
    index = findfirst(isequal(basis[1]), getfield(nvs, :contents)[1])
    @assert isa(index, Integer) "inds error: input basis out of range."
    for i = 2:length(NS)
        @assert getfield(nvs, :contents)[i][index] == basis[i] "inds error: input basis out of range."
    end
    return (index,)
end
@generated function inds(basis::NamedTuple{NS, BS}, nvs::NamedVectorSpace{:product, NS, BS}) where {NS, BS<:Tuple}
    return Expr(:tuple, [:(findfirst(isequal(basis[$i]), getfield(nvs, :contents)[$i])) for i = 1:length(NS)]...)
end
@generated function NamedTuple(index::Tuple{Int}, nvs::NamedVectorSpace{:zip, NS}) where NS
    return Expr(:tuple, [:($name = getfield(nvs, :contents)[$i][index[1]]) for (i, name) in enumerate(NS)]...)
end
@generated function NamedTuple(index::NTuple{N, Int}, nvs::NamedVectorSpace{:product, NS}) where {N, NS}
    @assert N == length(NS) "NamedTuple error: dismatched input index and rank of named vector space."
    return Expr(:tuple, [:($name = getfield(nvs, :contents)[$i][index[$i]]) for (i, name) in enumerate(NS)]...)
end

"""
    keys(nvs::NamedVectorSpace) -> Tuple{Vararg{Symbol}}
    keys(::Type{<:NamedVectorSpace{M, NS} where M}) where NS -> Tuple{Vararg{Symbol}}

Get the names of a named vector space.
"""
Base.keys(nvs::NamedVectorSpace) = keys(typeof(nvs))
Base.keys(::Type{<:NamedVectorSpace{M, NS} where M}) where {NS} = NS

"""
    values(nvs::NamedVectorSpace) -> Tuple{Vararg{AbstractVector}}

Get the contents of a named vector space.
"""
Base.values(nvs::NamedVectorSpace) = getfield(nvs, :contents)

"""
    pairs(nvs::NamedVectorSpace) -> Base.Iterators.Pairs

Get the name-value pairs of a named vector space.
"""
Base.pairs(nvs::NamedVectorSpace) = Base.Generator(=>, nvs|>keys, nvs|>values)

"""
    eltype(nvs::NamedVectorSpace, i::Int)
    eltype(::Type{<:NamedVectorSpace{M, NS, BS} where {M, NS}}, i::Int) where BS

Get the eltype of the ith indexable dimensions of a named vector space.
"""
Base.eltype(nvs::NamedVectorSpace, i::Int) = eltype(nvs|>typeof, i)
Base.eltype(::Type{<:NamedVectorSpace{M, NS, BS} where {M, NS}}, i::Int) where BS = fieldtype(BS, i)

end # module


text:
__precompile__()
module ArviZ

using Base: @__doc__
using Requires
using PyCall
using PyPlot
using Pandas
using NamedTupleTools

import Base: convert, propertynames, getproperty, hash, show, summary, +
import Base.Docs: getdoc
import PyCall: PyObject

# Exports

## Plots
export plot_autocorr,
       plot_compare,
       plot_density,
       plot_dist,
       plot_elpd,
       plot_energy,
       plot_ess,
       plot_forest,
       plot_hpd,
       plot_joint,
       plot_kde,
       plot_khat,
       plot_loo_pit,
       plot_mcse,
       plot_pair,
       plot_parallel,
       plot_posterior,
       plot_ppc,
       plot_rank,
       plot_trace,
       plot_violin

## Stats
export compare, hpd, loo, loo_pit, psislw, r2_score, waic

## Diagnostics
export bfmi, geweke, ess, rhat, mcse

## Stats utils
export autocov, autocorr, make_ufunc, wrap_xarray_ufunc

## Data
export InferenceData,
       convert_to_inference_data,
       load_arviz_data,
       to_netcdf,
       from_netcdf,
       from_dict,
       from_cmdstan,
       concat,
       concat!

## Utils
export interactive_backend

## rcParams
export rc_context

const arviz = PyNULL()

function __init__()
    copy!(arviz, pyimport_conda("arviz", "arviz", "conda-forge"))
    pyimport_conda("xarray", "xarray", "conda-forge")
    pyimport_conda("dask", "dask", "conda-forge")

    pytype_mapping(arviz.InferenceData, InferenceData)

    @require MCMCChains = "c7f686f2-ff18-58e9-bc7b-31028e88f75d" begin
        @require DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" include("mcmcchains.jl")
    end
end

include("utils.jl")
include("dataset.jl")
include("data.jl")
include("diagnostics.jl")
include("plots.jl")
include("stats.jl")
include("stats_utils.jl")

end # module


text:
module SubClonalSelection

# package code goes here
using CancerSeqSim
using DataFrames
using ApproxBayes
using Distributions
using Distances
using GLM
using Colors
using StatsBase
using Plots
using ProgressMeter
using CSV
using Statistics
using DelimitedFiles
using Printf
using Random

import Base.show, ApproxBayes.show

export
  # types
  Results,

  #functions
  fitABCmodels,
  plothistogram,
  plotmodelposterior,
  plotparameterposterior,
  saveallplots,
  timefunc,
  timefuncrand,
  saveresults,
  show

include("types.jl")
include("util.jl")
include("plotting.jl")
include("abcalgorithm.jl")
include("runabc.jl")

end # module


text:
using LinearAlgebra

function lineq(;m, v, invert=false)
    s = m \ v
    if ! invert
        return (; bigarrs=(; s, ))
    else
        invm=inv(m)
        detm = det(m)
        return (; detm, bigarrs=(; s, invm))
    end
end

function eigval(;m)
    s = ComplexF64.(eigen(m).values) # result is otherwise type unstable (float or complex), and LabVIEW can't handle it
    return (; bigarrs=(; s, ))
end

fns = (; lineq, eigval,)


text:
module _GridEn
export GridEn
using Plots
using Plots.PlotMeasures
using StatsBase: Histogram, fit
    """
        GDE, GDR, _ = GridEn(Sig)

    Returns the gridded distribution entropy (`GDE`) and the gridded
    distribution rate (`GDR`) estimated from the data sequence (`Sig`) using
    the default  parameters:
    grid coarse-grain = 3, time delay = 1, logarithm = base 2

        GDE, GDR, PIx, GIx, SIx, AIx = GridEn(Sig)

    In addition to GDE and GDR, GridEn returns the following indices
    estimated for the data sequence (`Sig`) using the default  parameters:
    [`PIx`]   - Percentage of points below the line of identity (LI)
    [`GIx`]   - Proportion of point distances above the LI
    [`SIx`]   - Ratio of phase angles (w.r.t. LI) of the points above the LI
    [`AIx`]   - Ratio of the cumulative area of sectors of points above the LI

        GDE, GDR, ..., = GridEn(Sig::AbstractArray{T,1} where T<:Real; m::Int=3, tau::Int=1, Logx::Real=exp(1), Plotx::Bool=false)

    Returns the gridded distribution entropy (`GDE`) estimated from the data
    sequence (`Sig`) using the specified 'keyword' arguments:

    # Arguments:
    `m`     - Grid coarse-grain (m x m sectors), an integer > 1 \n
    `tau`   - Time Delay, a positive integer    \n
    `Logx`  - Logarithm base, a positive scalar \n
    `Plotx` - When Plotx == true, returns gridded Poicaré plot and a bivariate
              histogram of the grid point distribution (default: false)  \n

    # See also `PhasEn`, `CoSiEn`, `SlopEn`, `BubbEn`, `MSEn`

    # References:
        [1] Chang Yan, et al.,
                "Novel gridded descriptors of poincaré plot for analyzing
                heartbeat interval time-series."
                Computers in biology and medicine
                109 (2019): 280-289.

        [2] Chang Yan, et al.
                "Area asymmetry of heart rate variability signal."
                Biomedical engineering online
                16.1 (2017): 1-14.

        [3] Alberto Porta, et al.,
                "Temporal asymmetries of short-term heart period variability
                are linked to autonomic regulation."
                American Journal of Physiology-Regulatory, Integrative and
                Comparative Physiology
                295.2 (2008): R550-R557.

        [4] C.K. Karmakar, A.H. Khandoker and M. Palaniswami,
                "Phase asymmetry of heart rate variability signal."
                Physiological measurement
                36.2 (2015): 303.

    """
    function GridEn(Sig::AbstractArray{T,1} where T<:Real; m::Int=3, tau::Int=1,
        Logx::Real=exp(1), Plotx::Bool=false)

    N = size(Sig,1)
    (N > 10) ? nothing : error("Sig:   must be a numeric vector")
    (m > 1) ? nothing :  error("m:     must be an integer > 1")
    (tau >0) ? nothing : error("tau:   must be an integer > 0")
    (Logx>0) ? nothing : error("Logx:  must be a positive number > 0")

    Sig_n = (Sig.-minimum(Sig))/(maximum(Sig)-minimum(Sig))
    Temp = hcat(Sig_n[1:end-tau], Sig_n[1+tau:end])
    Edges = collect(0:1/m:1); Edges[end] = 1.1
    N = fit(Histogram,(Temp[:,1],Temp[:,2]),(Edges,Edges)).weights
    Pj = reverse(transpose(N),dims=1)/size(Temp,1); Ppi = Pj[Pj.>0];
    if round(sum(Ppi),digits=4) != 1
        @warn("Potential error of estimated probabilities: P = $(sum(Ppi))")
    end
    GDE = -sum(Ppi.*(log.(Ppi)/log(Logx)));
    GDR = sum(N.!=0)/(m*m);

    T2   = atand.(Temp[:,2]./Temp[:,1])
    Dup  = sum(abs.(diff(Temp[T2.>45,:],dims=2)))
    Dtot = sum(abs.(diff(Temp[T2.!=45,:],dims=2)))
    Sup  = sum((T2[T2.>45].-45))
    Stot = sum(abs.(T2[T2.!=45].-45))
    Aup  = sum(abs.((T2[T2.>45].-45).*sqrt.(sum(Temp[T2.>45,:].^2,dims=2))))
    Atot = sum(abs.((T2[T2.!=45].-45).*sqrt.(sum(Temp[T2.!=45,:].^2,dims=2))))

    PIx = 100*sum(T2.<45)/sum(T2.!=45)
    GIx = 100*Dup/Dtot
    SIx = 100*Sup/Stot
    AIx = 100*Aup/Atot

    if Plotx
        ntrvl = range(0,1,length=m+1)
        p1 = plot(Sig_n[1:end-tau],Sig_n[tau+1:end],seriestype =:scatter,
            c=:magenta, markerstrokecolor=:magenta, markersize = 2)
        plot!(hcat(ntrvl,ntrvl)',hcat(zeros(m+1),ones(m+1))',c=:cyan,lw=3)
        plot!(hcat(zeros(m+1),ones(m+1))',hcat(ntrvl,ntrvl)',c=:cyan,lw=3)
        plot!([0, 1],[0, 1],c=:black, legend=false, aspect_ratio=:equal,
            xlim = (0, 1), ylim = (0, 1),
            xticks=[0, 1], yticks=[0, 1],  xlabel= "X ₙ" , ylabel="X ₙ₊ₜ")

        p2 = histogram2d(Temp[:,1],Temp[:,2], nbins=ntrvl,
            xticks=[0, 1], yticks=[0, 1], xlabel= "X ₙ" , ylabel="X ₙ₊ₜ",
            xlim = (0, 1), ylim = (0, 1), aspect_ratio=:equal,
            colorbar_entry=false, seriescolor=:cool,show_empty_bins=true)

        xx = plot(p1,p2,layout=(1,2),legend=false, framestyle=:grid, margin=7mm)
        display(xx)
    end

    return GDE, GDR, PIx, GIx, SIx, AIx
    end

end

"""
Copyright 2021 Matthew W. Flood, EntropyHub

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

For Terms of Use see https://github.com/MattWillFlood/EntropyHub

"""


text:
abstract type AbstractSqpModel end

mutable struct Model{T,TD} <: AbstractSqpModel
    n::Int  # Num vars
    m::Int  # Num cons
    x::TD  # Starting and final solution
    x_L::TD # Variables Lower Bound
    x_U::TD # Variables Upper Bound
    g::TD  # Final constraint values
    g_L::TD # Constraints Lower Bound
    g_U::TD # Constraints Upper Bound
    j_str::Array{Tuple{Int,Int}}
    h_str::Array{Tuple{Int,Int}}
    mult_g::TD # lagrange multipliers on constraints
    mult_x_L::TD # lagrange multipliers on lower bounds
    mult_x_U::TD # lagrange multipliers on upper bounds
    obj_val::T  # Final objective
    status::Int  # Final status

    # Callbacks
    eval_f::Function
    eval_g::Function
    eval_grad_f::Function
    eval_jac_g::Function
    eval_h::Union{Function,Nothing}

    num_linear_constraints::Int # number of linear constraints

    intermediate  # Can be nothing

    # For MathProgBase
    sense::Symbol

    parameters::Parameters
    statistics::Dict{String,Any}   # collects parameters of all iterations inside the algorithm if StatisticsFlag > 0

    Model(
        n::Int,
        m::Int,
        x_L::TD,
        x_U::TD,
        g_L::TD,
        g_U::TD,
        j_str::Array{Tuple{Int,Int}},
        h_str::Array{Tuple{Int,Int}},
        eval_f::Function,
        eval_g::Function,
        eval_grad_f::Function,
        eval_jac_g::Function,
        eval_h::Union{Function,Nothing},
        num_linear_constraints::Int,
        parameters::Parameters
    ) where {T, TD<:AbstractArray{T}} = new{T,TD}(
        n, m,
        zeros(n), x_L, x_U,
        zeros(m), g_L, g_U,
        j_str, h_str,
        zeros(m), zeros(n), zeros(n),
        0.0,
        -5,
        eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h,
        num_linear_constraints,
        nothing, :Min,
        parameters,
        Dict{String,Any}()
    )
end

function optimize!(model::Model)
    if isnothing(model.parameters.external_optimizer)
    	model.status = -12;
        @error "`external_optimizer` parameter must be set for subproblem solutions."
    else
        if model.parameters.algorithm == "SLP-LS"
            model.eval_h = nothing
            slp = SlpLS(model)
            run!(slp)
        elseif model.parameters.algorithm == "SLP-TR"
            model.eval_h = nothing
            slp = SlpTR(model)
            run!(slp)
        elseif model.parameters.algorithm == "SQP-TR"
            sqp = SqpTR(model)
            run!(sqp)
        else
            @warn "$(model.parameters.algorithm) is not defined"
        end
    end
    return nothing
end

function add_statistic(model::AbstractSqpModel, name::String, value)
    if model.parameters.StatisticsFlag == 0
        return
    end
    model.statistics[name] = value
end

function add_statistics(model::AbstractSqpModel, name::String, value::T) where T
    if model.parameters.StatisticsFlag == 0
        return
    end
    if !haskey(model.statistics, name)
        model.statistics[name] = Array{T,1}()
    end
    push!(model.statistics[name], value)
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    FunctionPartitioner(func)

A method for partitioning spatial objects with a given function `func`.
"""
struct FunctionPartitioner <: AbstractFunctionPartitioner
  func::Function
end

(p::FunctionPartitioner)(i, j) = p.func(i, j)


text:
"""
    cdad!(rbm, data)

Trains the RBM on data using contrastive divergence.
Computes gradients with Zygote.
"""
function cdad!(rbm::RBM, data::AbstractArray;
    batchsize = 1,
    epochs = 1,
    optim = ADAM(),
    wts = nothing,
    steps::Int = 1,
    callback = empty_callback
)
    @assert size(data) == (size(visible(rbm))..., size(data)[end])
    @assert isnothing(wts) || _nobs(data) == _nobs(wts)
    for epoch in 1:epochs
        batches = minibatches(data, wts; batchsize = batchsize)
        for (batch_idx, (vd, wd)) in enumerate(batches)
            vm = sample_v_from_v(rbm, vd; steps = steps)
            ∂, = Zygote.gradient(rbm) do rbm
                contrastive_divergence(rbm, vd, vm; wd = wd, wm = wd)
            end
            update!(∂, rbm, optim)
            update!(rbm, ∂)
            callback(; rbm, optim, epoch, batch_idx, vm, vd, wd)
        end
    end
    return rbm
end


text:
@testset "Convex Hull" begin
    A = zeros(50, 30)
    A= convert(Array{Bool}, A)
    A[25,1]=1
    A[1,10]=1
    A[10,10]=1
    A[10,30]=1
    A[40,30]=1
    A[40,10]=1
    A[50,10]=1
    B = @inferred convexhull(A)
    C = CartesianIndex{}[]
    push!(C, CartesianIndex{}(25,1))
    push!(C, CartesianIndex{}(1,10))
    push!(C, CartesianIndex{}(10,30))
    push!(C, CartesianIndex{}(40,30))
    push!(C, CartesianIndex{}(50,10))
    @test typeof(B)==Array{CartesianIndex{2},1}
    @test sort(B)==sort(C)

    A = [0.0, 0.0, 1.0, 0.0, 0.0,
         0.0, 1.0, 1.0, 0.0, 0.0,
         1.0, 0.0, 0.0, 1.0, 1.0,
         0.0, 0.0, 0.0, 0.0, 0.0,
         0.0, 0.0, 1.0, 0.0, 0.0]
    A = reshape(A, 5, 5)
    A = convert(Array{Bool}, A)
    B = B = @inferred convexhull(A)
    C = CartesianIndex{}[]
    push!(C, CartesianIndex{}(1,3))
    push!(C, CartesianIndex{}(3,1))
    push!(C, CartesianIndex{}(3,5))
    push!(C, CartesianIndex{}(5,3))
    @test typeof(B)==Array{CartesianIndex{2},1}
    @test sort(B)==sort(C)
end


text:
using DashBootstrapComponents

buttons = html_div([
    dbc_button("Regular", color = "primary", className = "me-1"),
    dbc_button("Active", color = "primary", active = true, className = "me-1"),
    dbc_button("Disabled", color = "primary", disabled = true),
]);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
function prepare_heartrate(directory::String, annotations_filename::String)
    annotations_path = joinpath(directory, annotations_filename)
    eda_path = joinpath(directory, "EDA.csv")
    hr_path = joinpath(directory, "HR.csv")

    #data = CSV.File(annotations_path; header = ["label","start","stop","first","step","last"]) |> DataFrame
    annotations = CSV.File(annotations_path) |> DataFrame
    eda = CSV.File(eda_path; header = ["EDA"]) |> DataFrame
    hr = CSV.File(hr_path; header = ["HR"]) |> DataFrame

    eda_timestamp = eda[1,1]
    eda_frequency = eda[2,1]

    hr_timestamp = hr[1,1]
    hr_frequency = hr[2,1]

    # Difference between the start of heart rate measurement versus eda measurement
    Δ = Int(hr_timestamp - eda_timestamp)

    # Number of actual heartrate samples
    N = size(hr,1) - 2
    # Number of annotations
    L, _ = size(annotations)

    for l = 1:L
        eda_start_index = annotations[l,:start]
        eda_stop_index = annotations[l,:stop]

        eda_start_seconds = floor(Int, eda_start_index / 4)
        eda_stop_seconds = floor(Int, eda_stop_index / 4)

        hr_start_seconds = eda_start_seconds - Δ
        hr_stop_seconds = eda_stop_seconds - Δ

        hr_start_index = Int(hr_start_seconds * hr_frequency)
        hr_stop_index = Int(hr_stop_seconds * hr_frequency)

        indicator = [hr_start_index <= i <= hr_stop_index ? 1 : 0 for i = 1:N]
        key = annotations[l,:label]
        hr = hcat(hr, DataFrame([vcat(hr_timestamp, hr_frequency, indicator)], [Symbol(key)]))
    end
    hr
    CSV.write(joinpath(directory, "HR_annotated.csv"), hr)
    #hr

    #eda
    #annotations
    #hr

    #hr[3:end,:]
    #hr = hr[3:end]

    #df = load_dataframe(annotations_path, ["EDA"])
end

function prepare_empatica(in_directory::String, out_directory::String, names::NamedTuple, offsets::NamedTuple)
    eda_path = joinpath(in_directory, "EDA.csv")
    hr_path = joinpath(in_directory, "HR.csv")
    tags_path = joinpath(in_directory, "tags.csv")
    ibi_path = joinpath(in_directory, "IBI.csv")

    eda = CSV.File(eda_path; header = ["EDA"]) |> DataFrame
    hr = CSV.File(hr_path; header = ["HR"]) |> DataFrame
    tags = CSV.File(tags_path; header = ["TAGS"]) |> DataFrame
    ibi = CSV.File(ibi_path; header = ["Timestamps", "IBI"]) |> DataFrame

    eda_timestamp = eda[1,1]
    eda_frequency = eda[2,1]

    Δt = (1 / eda_frequency) * 1000
    time₀ = Δt
    time₁ = Δt * length(eda[3:end,1])

    hr_timestamp = hr[1,1]
    hr_frequency = hr[2,1]
    hr_data = hr[3:end,1]

    ibi_timestamp = ibi[1,1]
    ibi_data = ibi[2:end,:]

    #tags

    eda_data = eda[3:end,1]

    eda_processed = prepare_eda(eda_data, eda_timestamp, eda_frequency, tags, names, offsets)
    CSV.write(joinpath(out_directory, "EDA_annotated.csv"), eda_processed)

    ibi_processed = prepare_ibi(ibi_data, ibi_timestamp, tags, names, offsets)
    CSV.write(joinpath(out_directory, "IBI_annotated.csv"), ibi_processed)

    ibi_processed = prepare_heart_rate(hr_data, hr_timestamp, hr_frequency, tags, names, offsets)
    CSV.write(joinpath(out_directory, "HR_annotated.csv"), ibi_processed)

    eda_processed
    #time₀, time₁, timestamps
end

function extract_tags(tags::DataFrame)
    tag_data = tags.TAGS
    if length(tag_data) != 2
        error("There need to be exactly two entries in tags.csv")
    end
    return tag_data[1], tag_data[2]
end

function prepare_eda(eda, unix_start, eda_frequency, tags::DataFrame, names::NamedTuple, offsets::NamedTuple)
    tag₁, tag₂ = extract_tags(tags)
    seconds_to_tag₁ = tag₁ - unix_start
    seconds_to_tag₂ = tag₂ - unix_start

    #@show seconds_to_tag₁,  seconds_to_tag₂

    interval₁, interval₂, interval₃, interval₄ = determine_intervals(seconds_to_tag₁, seconds_to_tag₂, offsets)

    N = length(eda)
    #@show N
    # Timestep expressed in fractions of a second.
    Δt = (1 / eda_frequency)
    timestamps = range(Δt, step = Δt, length = N)

    indicator₁ = [interval₁.start <= i <= interval₁.stop ? 1.0 : 0.0 for i in timestamps]
    indicator₂ = [interval₂.start <= i <= interval₂.stop ? 1.0 : 0.0 for i in timestamps]
    indicator₃ = [interval₃.start <= i <= interval₃.stop ? 1.0 : 0.0 for i in timestamps]
    indicator₄ = [interval₄.start <= i <= interval₄.stop ? 1.0 : 0.0 for i in timestamps]

    response_type = Lowpass(1; fs = eda_frequency)
    design_method = Butterworth(2)
    eda_filtered = filt(digitalfilter(response_type, design_method), eda)

    # TODO add filtered EDA
    df = DataFrame(Raw = vcat(unix_start, eda_frequency, eda))
    df = hcat(df, DataFrame(Filtered = vcat(unix_start, eda_frequency, eda_filtered)))
    df = hcat(df, DataFrame([vcat(unix_start, eda_frequency, indicator₁)], [Symbol(names.name_1)]))
    df = hcat(df, DataFrame([vcat(unix_start, eda_frequency, indicator₂)], [Symbol(names.name_2)]))
    df = hcat(df, DataFrame([vcat(unix_start, eda_frequency, indicator₃)], [Symbol(names.name_3)]))
    df = hcat(df, DataFrame([vcat(unix_start, eda_frequency, indicator₄)], [Symbol(names.name_4)]))

    return df
end

function prepare_ibi(ibi_data, unix_start, tags::DataFrame, names::NamedTuple, offsets::NamedTuple)
    tag₁, tag₂ = extract_tags(tags)
    seconds_to_tag₁ = tag₁ - unix_start
    seconds_to_tag₂ = tag₂ - unix_start

    interval₁, interval₂, interval₃, interval₄ = determine_intervals(seconds_to_tag₁, seconds_to_tag₂, offsets)

    N = size(ibi_data,2)
    #@show N
    timestamps = ibi_data[:,1]
    interbeat_intervals = ibi_data[:,2]

    indicator₁ = [interval₁.start <= i <= interval₁.stop ? 1.0 : 0.0 for i in timestamps]
    indicator₂ = [interval₂.start <= i <= interval₂.stop ? 1.0 : 0.0 for i in timestamps]
    indicator₃ = [interval₃.start <= i <= interval₃.stop ? 1.0 : 0.0 for i in timestamps]
    indicator₄ = [interval₄.start <= i <= interval₄.stop ? 1.0 : 0.0 for i in timestamps]

    df = DataFrame(Timestamp = vcat(unix_start, timestamps))
    df = hcat(df, DataFrame(IBI = vcat("IBI", interbeat_intervals)))
    df = hcat(df, DataFrame([vcat(unix_start, indicator₁)], [Symbol(names.name_1)]))
    df = hcat(df, DataFrame([vcat(unix_start, indicator₂)], [Symbol(names.name_2)]))
    df = hcat(df, DataFrame([vcat(unix_start, indicator₃)], [Symbol(names.name_3)]))
    df = hcat(df, DataFrame([vcat(unix_start, indicator₄)], [Symbol(names.name_4)]))

    return df
end

function prepare_heart_rate(hr_data, unix_start, hr_frequency, tags::DataFrame, names::NamedTuple, offsets::NamedTuple)
    tag₁, tag₂ = extract_tags(tags)
    seconds_to_tag₁ = tag₁ - unix_start
    seconds_to_tag₂ = tag₂ - unix_start

    #@show seconds_to_tag₁,  seconds_to_tag₂

    interval₁, interval₂, interval₃, interval₄ = determine_intervals(seconds_to_tag₁, seconds_to_tag₂, offsets)

    N = length(hr_data)
    #@show N
    # Timestep expressed in fractions of a second.
    Δt = (1 / hr_frequency)
    timestamps = range(Δt, step = Δt, length = N)

    indicator₁ = [interval₁.start <= i <= interval₁.stop ? 1.0 : 0.0 for i in timestamps]
    indicator₂ = [interval₂.start <= i <= interval₂.stop ? 1.0 : 0.0 for i in timestamps]
    indicator₃ = [interval₃.start <= i <= interval₃.stop ? 1.0 : 0.0 for i in timestamps]
    indicator₄ = [interval₄.start <= i <= interval₄.stop ? 1.0 : 0.0 for i in timestamps]

    df = DataFrame(Raw = vcat(unix_start, hr_frequency, hr_data))
    df = hcat(df, DataFrame([vcat(unix_start, hr_frequency, indicator₁)], [Symbol(names.name_1)]))
    df = hcat(df, DataFrame([vcat(unix_start, hr_frequency, indicator₂)], [Symbol(names.name_2)]))
    df = hcat(df, DataFrame([vcat(unix_start, hr_frequency, indicator₃)], [Symbol(names.name_3)]))
    df = hcat(df, DataFrame([vcat(unix_start, hr_frequency, indicator₄)], [Symbol(names.name_4)]))

    return df
end

function determine_intervals(seconds_to_tag₁::Number, seconds_to_tag₂::Number, offsets::NamedTuple)
    interval₁ = (start = seconds_to_tag₁ - offsets.offset_1,
                 stop =  seconds_to_tag₁)

    interval₂ = (start = seconds_to_tag₁,
                 stop =  seconds_to_tag₁ + offsets.offset_2)

    interval₃ = (start = seconds_to_tag₂ - offsets.offset_3,
                 stop =  seconds_to_tag₂)

    interval₄ = (start = seconds_to_tag₂,
                 stop =  seconds_to_tag₂ + offsets.offset_4)

    return interval₁, interval₂, interval₃, interval₄
end

# offsets = (offset_1 = seconds_prior_marker_1,
#            offset_2 = seconds_after_marker_1,
#            offset_3 = seconds_prior_marker_2,
#            offset_4 = seconds_after_marker_2)


text:
module Viability

using Base
using NearestNeighbors

# type definitions
COORDINATE_TYPE = Float64
POINT_TYPE = Array{COORDINATE_TYPE,1}
POINTS_ARRAY_TYPE = Array{COORDINATE_TYPE,2}
STATE_TYPE = UInt8
STATE_ARRAY_TYPE = Array{STATE_TYPE,1}
STATE_ARRAY_ARGUMENT_TYPE = AbstractArray{STATE_TYPE,1}
SHARED_STATE_ARRAY_TYPE = SharedArray{STATE_TYPE,1}

@enum SP_COMPUTATION_TYPE sp_parallel sp_mparallel

# computational globals
MAX_SP_VIABILITY_ITERATION_STEPS = 100

RequiredArgument(name) = error("$name is a required keyword argument")
RequiredArgument(name, explanation) = error("$name is a required keyword argument. $explanation")

function map_over_points{T}(::Type{T}, f, points::POINTS_ARRAY_TYPE)
    @assert ndims(points) == 2
    npoints = size(points, 2)
    res = Array{T}(npoints)
    for i = 1:npoints
        res[i] = f(points[:, i])
    end
    return res
end

function create_normalized_2d_grid(
    n_per_dim::Unsigned;
    epsilon::COORDINATE_TYPE=0.0 # how much space should be left around the boundary
    )
    create_2d_grid(
        n_per_dim,
        COORDINATE_TYPE[0 1; 0 1],
        epsilon=epsilon
    )
end

function create_2d_grid(
    n_per_dim::Unsigned,
    bounds::Array{COORDINATE_TYPE,2};
    epsilon::COORDINATE_TYPE=0.0 # how much space should be left around the boundary
    )
    dim = 2
    @assert ndims(bounds) == 2
    @assert size(bounds, 1) == dim
    @assert size(bounds, 2) == 2
    @assert epsilon >= 0
    grid = POINTS_ARRAY_TYPE(dim, n_per_dim^dim)
    start_points = bounds[:,1] + epsilon
    steps = (bounds[:,2] - bounds[:,1] - 2*epsilon) / (n_per_dim - 1)
    @assert all(steps .> 0.0) "check you bounds and epsilon values"
    for i = 1:n_per_dim, j = 1:n_per_dim
        index  = i + (j-1)*n_per_dim
        grid[1, index] = start_points[1] + steps[1] * (i-1)
        grid[2, index] = start_points[2] + steps[2] * (j-1)
    end
    return grid
end

function point2bounds(
    point::POINT_TYPE,
    bounds::Array{COORDINATE_TYPE,2};
    epsilon::COORDINATE_TYPE=0.0
    )
    offsets = bounds[:,1] + epsilon
    diffs = bounds[:,2] - bounds[:,1] - 2*epsilon
    point.*diffs + offsets
end

function points2bounds(
    points::POINTS_ARRAY_TYPE,
    bounds::Array{COORDINATE_TYPE,2};
    epsilon::COORDINATE_TYPE=0.0
    )
    new_points = zeros(points)
    n = size(points,2)
    for i=1:n
        new_points[:, i] = point2bounds(points[:, i], bounds, epsilon=epsilon)
    end
    return new_points
end

function temporal_homogenization(f; epsilon::COORDINATE_TYPE=1e-6)
    function f_homogenized(x::POINT_TYPE)
        xdot = f(x)
        xdot / (norm(xdot) + epsilon)
    end
end

function spatial_normalization(
        f,
        bounds::Array{COORDINATE_TYPE,2};
        epsilon::COORDINATE_TYPE=0.0 # how much space should be left around the boundary
    )
    offsets = bounds[:,1] + epsilon
    diffs = bounds[:,2] - bounds[:,1] - 2*epsilon
    # println(offsets)
    # println(diffs)
    function f_normalized(y::POINT_TYPE)
        f(y.*diffs + offsets)./diffs
    end
end

function spatial_compactification(
        f,
        xi_mid::COORDINATE_TYPE,
        axis::Int
    )
    function f_compactified(y)
        x = copy(y)
        x[axis] = y[axis] * xi_mid / (1-y[axis])
        p = copy(y)
        p[axis] = (1-y[axis])^2 / xi_mid
        p*f(x)
    end
end

function spatial_compactification(
    f,
    x_mid::POINT_TYPE,
    axes::Array{Int,1}
    )
    @assert len(x_mid) == len(axes)
    for (xi_mid, axis) in zip(xs_mid, axes)
        f = spatial_compactification(f, xi_mid, axis)
    end
    f
end

function rhs2stepfct(rhs, delta_t)
    function stepfct(x)
        x + rhs(x)*delta_t
    end
    return stepfct
end

function get_indices_of_neighbors_in_ball(
    point::POINT_TYPE,
    points::POINTS_ARRAY_TYPE,
    epsilon::COORDINATE_TYPE # radius of the ball
    )
    @assert ndims(point) == 1
    dim = size(point, 1)
    @assert ndims(points) == 2
    @assert size(points, 1) == dim
    @assert epsilon > 0
    point_diffs = points .- point
    diffs = sum(point_diffs.^2, [1])
    neighbor_indices = find(diffs .< epsilon^2)
end

function intersects(a::Array, b::Array)
    for i=1:size(a, 1), j=1:size(b, 1)
        if a[i] == b[j]
            return true
        end
    end
    return false
end

function sp_viability_kernel(
        points::POINTS_ARRAY_TYPE,
        states::STATE_ARRAY_ARGUMENT_TYPE,
        step_functions::Array;
        good_states::STATE_ARRAY_ARGUMENT_TYPE=STATE_TYPE[1],
        bad_states::STATE_ARRAY_ARGUMENT_TYPE=STATE_TYPE[0],
        computation_mode::SP_COMPUTATION_TYPE=sp_mparallel,
        delta::COORDINATE_TYPE=RequiredArgument("delta"),
        r::COORDINATE_TYPE=0.0,
        leafsize=10
    )
    ret_val = 2

    # create a kdtree that is used for the find the surrounding indices during the computation
    kdtree = KDTree(points; leafsize=leafsize)
    # create a shared array copy of states for parallelization
    inp_states = states
    states = convert(SharedArray, states)
    if computation_mode == sp_parallel
        ret_val = _sp_viability_kernel_parallel(points, states, step_functions, kdtree,
            good_states=good_states,
            unsuccessful_states=bad_states,
            successful_states=good_states,
            work_states=good_states,
            delta=delta
        )
    elseif computation_mode == sp_mparallel
        if r == 0.0
            RequiredArgument(
                "r",
                "If unsure, choose r to be the (maximal) step length of the step_function.")
        end
        ret_val = _sp_viability_kernel_mparallel(points, states, step_functions, kdtree,
            good_states=good_states,
            unsuccessful_states=bad_states,
            successful_states=good_states,
            work_states=good_states,
            delta=delta,
            r=r
        )
    else
        throw(ArgumentError("unknown computation_mode=$computation_mode"))
    end
    # copy back from shared array
    inp_states[:] = states
    return ret_val
end

function sp_capture_basin(
        points::POINTS_ARRAY_TYPE,
        states::STATE_ARRAY_TYPE,
        step_functions::Array;
        target_states::STATE_ARRAY_TYPE=STATE_TYPE[1],
        work_states::STATE_ARRAY_TYPE=STATE_TYPE[0],
        reached_states::STATE_ARRAY_TYPE=STATE_TYPE[1],
        delta::COORDINATE_TYPE=RequiredArgument("delta")
    )
    _sp_viability_kernel_parallel(points, states, step_functions,
        good_states=target_states+reached_states,
        work_states=work_states,
        successful_states=reached_states,
        unsuccessful_states=work_states,
        delta=delta
    )
end

function _sp_viability_kernel_mparallel(
        points::POINTS_ARRAY_TYPE,
        states::STATE_ARRAY_ARGUMENT_TYPE,
        step_functions::Array,
        kdtree::NearestNeighbors.KDTree;
        good_states::STATE_ARRAY_ARGUMENT_TYPE=RequiredArgument("good_states"),
        work_states::STATE_ARRAY_ARGUMENT_TYPE=RequiredArgument("work_states"),
        successful_states::STATE_ARRAY_ARGUMENT_TYPE=RequiredArgument("successful_states"),
        unsuccessful_states::STATE_ARRAY_ARGUMENT_TYPE=RequiredArgument("unsuccessful_states"),
        delta::COORDINATE_TYPE=RequiredArgument("delta"),
        r::COORDINATE_TYPE=RequiredArgument("r")
    )
    # check the input for correctness and consistency
    @assert ndims(points) == 2
    @assert ndims(states) == 1
    dim = size(points,1)
    n = size(points,2)
    @assert size(states, 1) == n
    @assert size(successful_states) == size(work_states)
    @assert size(unsuccessful_states) == size(work_states)

    # compute r_min and r_max
    r_min = r - delta - 0.01
    # TODO: check whether r_min really makes sense
    r_max = r + delta + 0.01
    println("r_min $r_min")
    println("r_max $r_max")
    # create an array where the points are marked if they have to be checked
    # in the next run
    # initialize all to true, because in the beginning, all points have to be
    # visited at least once
    marked_points = SharedArray{Bool,1}(size(states))
    marked_points[:] = true
    retval = 1

    for iteration_step = 1:MAX_SP_VIABILITY_ITERATION_STEPS
        # get the indices of all points that need to be visited in this iteration
        run_indices = find(marked_points)
        if isempty(run_indices)
            println("yay, I am done")
            retval = 0
            break
        end
        marked_points[:] = false
        println("start iteration $iteration_step")
        @sync @parallel for i in run_indices
            # check whether the state of i is a work state and if yes find the index
            work_state_index = findfirst(work_states, states[i])
            if work_state_index != 0
                # yes, it is a work state

                old_state = states[i] # save it for the comparison later
                # evaluate which new state
                states[i] = _sp_get_new_state(
                    work_state_index, points[:, i], states, step_functions, kdtree,
                    good_states=good_states,
                    successful_states=successful_states,
                    unsuccessful_states=unsuccessful_states,
                    delta=delta
                )
                if old_state != states[i] # did something change?)
                    # mark all points with a distance between r_min and r_max,
                    # because these are the only ones that could reach points[:, i]
                    surrounding_indices = inrange(kdtree, points[:, i], r_max, false)
                    distances = map( index->norm(points[:, index] - points[:, i]), surrounding_indices )
                    for (new_index, new_distance) in zip(surrounding_indices, distances)
                        if !(marked_points[new_index]) && # it's not marked already
                                new_distance > r_min && # distance between r_min and r_max
                                states[new_index] in work_states # it's a work state
                            # mark point for revisiting if it can reach 'point'
                            marked_points[new_index] = true
                        end
                    end
                end
            end
        end
    end
    return retval
end

function _sp_viability_kernel_parallel(
        points::POINTS_ARRAY_TYPE,
        states::STATE_ARRAY_ARGUMENT_TYPE,
        step_functions::Array,
        kdtree::NearestNeighbors.KDTree;
        good_states::STATE_ARRAY_ARGUMENT_TYPE=RequiredArgument("good_states"),
        work_states::STATE_ARRAY_ARGUMENT_TYPE=RequiredArgument("work_states"),
        successful_states::STATE_ARRAY_ARGUMENT_TYPE=RequiredArgument("successful_states"),
        unsuccessful_states::STATE_ARRAY_ARGUMENT_TYPE=RequiredArgument("unsuccessful_states"),
        delta::COORDINATE_TYPE=RequiredArgument("delta")
    )
    # check the input for correctness and consistency
    @assert ndims(points) == 2
    @assert ndims(states) == 1
    dim = size(points,1)
    n = size(points,2)
    @assert size(states, 1) == n
    @assert size(successful_states) == size(work_states)
    @assert size(unsuccessful_states) == size(work_states)

    # changed = true
    retval = 1 # unsuccesful
    for it = 1:MAX_SP_VIABILITY_ITERATION_STEPS
        println("start iteration $it")
        changed =  @parallel (|) for i = 1:n
            # check whether the state of i is a work state and if yes find the index
            work_state_index = findfirst(work_states, states[i])
            if work_state_index != 0
                # yes, it is a work state

                old_state = states[i] # save it for the comparison later
                # evaluate which new state
                states[i] = _sp_get_new_state(
                    work_state_index, points[:, i], states, step_functions, kdtree,
                    good_states=good_states,
                    successful_states=successful_states,
                    unsuccessful_states=unsuccessful_states,
                    delta=delta
                )
                old_state != states[i] # did something change?
            else
                false # return that nothing changed
            end
        end
        if !changed # nothing changed? yay, I am done
            println("nothing changed, yay I am done")
            retval = 0 # succesfully finished
            break
        end
    end
    return retval
end

@inline function _sp_get_new_state(
    work_state_index::Int64,
    point::POINT_TYPE,
    states::AbstractArray{STATE_TYPE,1},
    step_functions::Array,
    kdtree::NearestNeighbors.KDTree;
    good_states::STATE_ARRAY_TYPE=RequiredArgument("good_states"),
    successful_states::STATE_ARRAY_TYPE=RequiredArgument("successful_states"),
    unsuccessful_states::STATE_ARRAY_TYPE=RequiredArgument("unsuccessful_states"),
    delta::COORDINATE_TYPE=RequiredArgument("delta")
    )
    for sf in step_functions
        next_point = sf(point)
        surrounding_indices = inrange(kdtree, next_point, delta, false)
        if isempty(surrounding_indices)
            # take the closest point if there was none found in the delta-neighborhood of next_point
            surrounding_indices = knn(kdtree, next_point, 1)[1]
        end
        if intersects(states[surrounding_indices], good_states)
            found_good_state = true
            return successful_states[work_state_index]
        end
    end
    # no good state found, so return the unsuccesful state
    unsuccessful_states[work_state_index]
end

STATE_IDENTIFICATION_VIABILITY_KERNEL = Dict{String,STATE_TYPE}(
    "good" => 1,
    "bad" => 0
)

STATE_IDENTIFICATION_TSM = Dict{String,STATE_TYPE}(
    "unknown" => 0,
    "shelter" => 1,
    "glade" => 2,
    "lake" => 3,
    "remaining sunny upstream" => 4,
    "dark upstream" => 5,
    "backwaters" => 6,
    "remaining sunny downstream" => 7,
    "dark downstream" => 8,
    "sunny eddies" => 9,
    "dark eddies" => 10,
    "sunny abyss" => 11,
    "dark abyss" => 12,
    "trench" => 13
)

function write_result_file(fname, model_info, points, states, delim="; ")
    # check the input for correctness and consistency
    @assert ndims(points) == 2
    @assert ndims(states) == 1
    n = size(points,2)
    @assert size(states, 1) == n
    # TODO: added tests for model_info
    open(fname, "w") do f
        write(f, "#file-version:0.1\n")
        write(f, "#$model_info\n")
        for i = 1:n
            point = points[:, i]
            state = states[i]
            write(f, "$point$delim$state\n")
        end
    end
end

end