text:
# ---
# title: 442. Find All Duplicates in an Array
# id: problem442
# author: Indigo
# date: 2021-06-21
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/find-all-duplicates-in-an-array/description/>
# hidden: true
# ---
#
# Given an array of integers, 1 ≤ a[i] ≤ _n_ ( _n_ = size of array), some
# elements appear **twice** and others appear **once**.
#
# Find all the elements that appear **twice** in this array.
#
# Could you do it without extra space and in O( _n_ ) runtime?
#
# **Example:**
#
#
#
#     Input:
#     [4,3,2,7,8,2,3,1]
#
#     Output:
#     [2,3]
#
#
#
## @lc code=start
using LeetCode

function find_duplicates_442(nums::Vector{Int})
    res = Int[]
    for num in nums
        absn = abs(num)
        nums[absn] < 0 ? push!(res, absn) : (nums[absn] = -nums[absn])
    end
    return res
end
## @lc code=end


text:
using Flux
using BSON
using LinearAlgebra
using Statistics
using Random
using Zygote: @adjoint
using Zygote: @nograd
using Zygote
using ProgressBars
using Distributions

include("ctc.jl")

Random.seed!(1)

const TRAINDIR = "train"

const EPOCHS = 100
const BATCH_SIZE = 1

losses = []

forward = LSTM(26, 100)
backward = LSTM(26, 100)
output = Dense(200, 62)

const NOISE = Normal(0, 0.6)

function m(x)
  h0f = collect(map(forward, x))
  h0b = reverse(map(backward, reverse(x)))
  h0 = map(i -> vcat(h0f[i], h0b[i]), 1:length(h0f)) |> collect
  o = collect(map(output, h0))
  return o
end

function loss(x, y)
  x = addNoise(x)
  Flux.reset!((forward, backward))
  yhat = m(x)
  yhat = reduce(hcat, yhat)
  l = ctc(yhat, y)
  addToGlobalLoss(l)
  return l
end

@nograd function addNoise(x)

  x = deepcopy(x)
  n = [rand(NOISE, 26) for i in 1:length(x)]
  for i in 1:length(x)
    x[i] .+= n[i]
  end

  return x
end

function readData(dataDir)
  # fnames = open(readlines, "shuffled_names.txt")
  fnames = readdir(dataDir)
  shuffle!(MersenneTwister(4), fnames)
  Xs = []
  Ys = []

  for fname in fnames
    BSON.@load joinpath(dataDir, fname) x y
    x = [Float32.(x[i,:]) for i in 1:size(x,1)]
    push!(Xs, x)
    push!(Ys, Array(y'))
  end

  m = mean(reduce(vcat, Xs))
  st = std(reduce(vcat, Xs))
  for (i, x) in enumerate(Xs)
    Xs[i] = [(xI .- m) ./ st for xI in x]
  end

  return (Xs, Ys)
end

function lev(s, t)
    m = length(s)
    n = length(t)
    d = Array{Int}(zeros(m+1, n+1))

    for i=2:(m+1)
        @inbounds d[i, 1] = i-1
    end

    for j=2:(n+1)
        @inbounds d[1, j] = j-1
    end

    for j=2:(n+1)
        for i=2:(m+1)
            @inbounds if s[i-1] == t[j-1]
                substitutionCost = 0
            else
                substitutionCost = 1
            end
            @inbounds d[i, j] = min(d[i-1, j] + 1, # Deletion
                            d[i, j-1] + 1, # Insertion
                            d[i-1, j-1] + substitutionCost) # Substitution
        end
    end

    @inbounds return d[m+1, n+1]
end

"""
    collapse(seq)

Collapses `seq` into a version that has symbols collapsed into one repetition and removes all instances
of the blank symbol.
"""
function collapse(seq)
  s = [x for x in seq if x != 62]
  if isempty(s) return s end
  s = [seq[1]]
  for ch in seq[2:end]
    if ch != s[end] && ch != 62
      push!(s, ch)
    end
  end
  return s
end

"""
    per(x, y)

Compute the phoneme error rate of the model for input `x` and target `y`. The phoneme error rate
is defined as the Levenshtein distance between the labeling produced by running `x` through
the model and the target labeling in `y`, all divided by the length of the target labeling
in `y`
"""
function per(x, y)
  Flux.reset!((forward, backward))
  yhat = m(x)
  yhat = reduce(hcat, yhat)
  yhat = mapslices(argmax, yhat, dims=1) |> vec |> collapse
  y = mapslices(argmax, y, dims=1) |> vec |> collapse
  return lev(yhat, y) / length(y)
end

function addToGlobalLoss(x)
	global losses
	push!(losses, x)
end

@adjoint function addToGlobalLoss(x)
	addToGlobalLoss(x)
	return nothing, () -> nothing
end

function main()
  println("Loading files")
  Xs, Ys = readData(TRAINDIR)
  data = collect(zip(Xs, Ys))

  Xs = [d[1] for d in data]
  Ys = [d[2] for d in data]

  println("Beginning training")

  data = zip(Xs, Ys) |> collect
  valData = data[1:186]
  data = data[187:end]

  opt = Momentum(1e-4)

  for i in 1:EPOCHS
    global losses
    losses = []
    println("Beginning epoch $i/$EPOCHS")
    Flux.train!(loss, Flux.params((forward, backward, output)), ProgressBar(data), opt)
    println("Calculating PER...")
    p = mean(map(x -> per(x...), valData))
    println("PER: $(p*100)")

    println("Mean loss: ", mean(losses))
  end
end

main()

ext: jl

max_stars_repo_name: maetshju/flux-ctc-grad

max_issues_repo_name: maetshju/flux-ctc-grad

max_forks_repo_name: maetshju/flux-ctc-grad


text:
using SMLMFrameConnection
using Documenter

DocMeta.setdocmeta!(SMLMFrameConnection, :DocTestSetup, :(using SMLMFrameConnection); recursive=true)

makedocs(;
    modules=[SMLMFrameConnection],
    authors="klidke@unm.edu",
    repo="https://github.com/JuliaSMLM/SMLMFrameConnection.jl/blob/{commit}{path}#{line}",
    sitename="SMLMFrameConnection.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaSMLM.github.io/SMLMFrameConnection.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaSMLM/SMLMFrameConnection.jl",
    devbranch = "main"
)


text:
"""
    __LOG_INDICATOR__()

Helper struct to specify a method as logging tool.
"""
struct __LOG_INDICATOR__
end

"""
# Notes
```julia
x = [(; d1 = "A1"),
     (; d1b=(; d2a="B1")),
     (; d1b=(; d2b="C1")),
     (; d1b=(; d2c=(; d3a="D1"))),
     (; d1b=(; d2c=(; d3b="E1"))),]
SimulationLogger.recursive_merge(x...)
julia> recursive_merge(x...)
(d1 = "A1", d1b = (d2a = "B1", d2b = "C1", d2c = (d3a = "D1", d3b = "E1")))
```
# References
https://discourse.julialang.org/t/multi-layer-dict-merge/27261/2?u=ihany
"""
function recursive_merge(x::NamedTuple...)
    merge!!(recursive_merge, x...)
end
function recursive_merge(x...)
    error("Invalid recursive_merge; did you assign values with the same key?
          For example,
          ```julia
          @nested_log :values x = 1
          @nested_log :values x = 2
          ```
          Otherwise, please report an issue.")
end

"""
    @Loggable(defun)

Generate two methods from a function definition (`defun`).

# Notes
Function definition with macro `@Loggable` will generate two methods with the same name (generic function).
For example,

```julia
@Loggable function dynamics!(dx, x, p, t; kwargs...)
    @log x
    dx = -x
end
```
will generate

```julia
function dynamics!(dx, x, p, t; kwargs...)
    @log x
    dx = -x
end
function dynamics!(dx, x, p, t, __log_indicator__::__LOG_INDICATOR__; kwargs...)
    __LOGGER__ = @isdefined(:__LOGGER__) ? __LOGGER__ : NamedTuple()  # if isdefined, nested logging will work
    @log x
    dx = -x
    return __LOGGER__
end
```
"""
macro Loggable(defun)
    _def = splitdef(defun)  # original
    def = deepcopy(_def)
    _body = _def[:body]
    push!(def[:args], :(__log_indicator::__LOG_INDICATOR__))
    def[:body] = quote
        __LOGGER__ = @isdefined($(:__LOGGER__)) ? __LOGGER__ : NamedTuple()
        $(_body.args...)  # remove the last line, return, to return __LOGGER__
        __LOGGER__  # return a NamedTuple
    end
    res = quote
        $(MacroTools.combinedef(_def))  # the original method
        $(MacroTools.combinedef(def))  # the modified method
        $(def[:name])  # return methods
    end
    esc(res)
end

"""
    @log(__LOGGER__, expr)

A macro which logs annotated data (based on `expr`) to a privileged data name (`__LOGGER__`).
For example,
```julia
@Loggable function dynamics!(dx, x, p, t; kwargs...)
    @log x
    dx = -x
end
```
will generate

```julia
function dynamics!(dx, x, p, t; kwargs...)
    @log x
    dx = -x
end
function dynamics!(dx, x, p, t, __log_indicator__::__LOG_INDICATOR__; kwargs...)
    __LOGGER__ = @isdefined(:__LOGGER__) ? __LOGGER__ : NamedTuple()  # if isdefined, nested logging will work
    @log x
    dx = -x
    return __LOGGER__
end
```
Then, if `@isdefined(__LOGGER__) == false`,
just evaluate given experssion `expr`.
Otherwise, @log(`expr`) will add a variable to `__LOGGER__` based on `expr`,
and also evaluate given experssion `expr`.
"""
macro log(__LOGGER__, expr)
    return if expr isa Symbol  # @log x
        quote
            local val = $(esc(expr))
            local var_name = $((expr,))[1]
            local logger_ = $(esc(__LOGGER__))
            logger_ = haskey(logger_, var_name) ? error("Already defined key: $(var_name)") : $setindex!!(logger_, val, var_name)
            $(esc(__LOGGER__)) = logger_
            nothing
        end
    elseif expr.head == :(=)
        if expr.args[1] isa Expr && expr.args[1].head == :tuple  # @log x, y = a, b
            quote
                local vals = $(esc(expr.args[2]))
                local var_names = $(esc(expr.args[1].args))
                local logger_ = $(esc(__LOGGER__))
                for (var_name, val) in zip(var_names, vals)
                    logger_ = haskey(logger_, var_name) ? error("Already defined key: $(var_name)") : $setindex!!(logger_, val, var_name)
                end
                $(esc(__LOGGER__)) = logger_
                $(esc(expr))
            end
        else  # @log x = a
            quote
                local val = $(esc(expr.args[2]))
                local var_name = $((expr.args[1],))[1]
                local logger_ = $(esc(__LOGGER__))
                logger_ = haskey(logger_, var_name) ? error("Already defined key: $(var_name)") : $setindex!!(logger_, val, var_name)
                $(esc(__LOGGER__)) = logger_
                $(esc(expr))
            end
        end
    elseif expr.args isa Array  # @log x, y
        quote
            local vals = $(esc(expr))
            local var_names = $(esc((expr.args)))
            local logger_ = $(esc(__LOGGER__))
            for (var_name, val) in zip(var_names, vals)
                logger_ = haskey(logger_, var_name) ? error("Already defined key: $(var_name)") : $setindex!!(logger_, val, var_name)
            end
            $(esc(__LOGGER__)) = logger_
            $(esc(expr))
        end
    else
        :(error("To log a variable, use either one of forms: `@log val` or `@log var_name = val`"))
    end
end

macro log(expr)
    esc(:(@isdefined($:__LOGGER__) ? @log($:__LOGGER__, $expr) : $expr))
end

"""
    @onlylog(expr)

A macro that activates given expression (`expr`) only when logging data.
Unlike `@log(expr)`,
this macro does not evaluate given experssion `expr`.
"""
macro onlylog(expr)
    esc(:(@isdefined($:__LOGGER__) ? @log($:__LOGGER__, $expr) : nothing))
end

"""
    @nested_log(symbol, expr)

A macro that enables us to log data in a nested sense.
# Examples
- Example 1
```julia
@nested_log :subsystem dynamics!(dx.sub, x.sub, p.sub, t)
```
will log data from `dynamics!(dx.sub, x.sub, p.sub, t)` as
`__LOGGER__.subsystem`.
- Example 2
```julia
@nested_log :subsystem state = x
```

# NOTICE
If you assign values with the two same keys, it will yield an error looks like:
```julia
ERROR: MethodError: no method matching recursive_merge(::Int64, ::Int64)
```
"""
macro nested_log(__LOGGER__, symbol, expr)
    if expr isa Symbol
        res = quote
            local logger_ = $(__LOGGER__)
            __TMP__ = NamedTuple()
            @log(__TMP__, $expr)
            if $symbol == nothing
                logger_ = $recursive_merge([logger_, __TMP__]...)
            else
                if haskey(logger_, $symbol)
                    logger_ = $setindex!!(logger_, $recursive_merge([logger_[$symbol], __TMP__]...), $symbol)
                else
                    logger_ = $setindex!!(logger_, __TMP__, $symbol)
                end
            end
            $(__LOGGER__) = logger_
            $expr
        end
        esc(res)
    elseif expr.head == :call
        _expr = copy(expr)
        push!(_expr.args, :(__LOG_INDICATOR__()))
        res = quote
            local logger_ = $(__LOGGER__)
            __TMP__ = NamedTuple()
            __TMP__ = $recursive_merge([__TMP__, $_expr]...)
            if $symbol == nothing
                logger_ = $recursive_merge([logger_, __TMP__]...)
            else
                if haskey(logger_, $symbol)
                    logger_ = $setindex!!(logger_, $recursive_merge([logger_[$symbol], __TMP__]...), $symbol)
                else
                    logger_ = $setindex!!(logger_, $_expr, $symbol)
                end
            end
            $(__LOGGER__) = logger_
            $expr
        end
        esc(res)
    elseif expr.head == :(=)  # @nested_log env_name x = a  or  @nested_log env_name x, y = a, b
        res = quote
            local logger_ = $(__LOGGER__)
            __TMP__ = NamedTuple()
            @log(__TMP__, $expr)
            if $symbol == nothing
                logger_ = $recursive_merge([logger_, __TMP__]...)  # 'cause it is not in-place.
            else
                if haskey(logger_, $symbol)
                    logger_ = $setindex!!(logger_, $recursive_merge([logger_[$symbol], __TMP__]...), $symbol)
                else
                    logger_ = $setindex!!(logger_, __TMP__, $symbol)
                end
            end
            $(__LOGGER__) = logger_
            $expr
        end
        esc(res)
    elseif expr.args isa Array  # @nested_log env_name a, b
        res = quote
            local logger_ = $(__LOGGER__)
            __TMP__ = NamedTuple()
            @log(__TMP__, $expr)
            if $symbol == nothing
                logger_ = $recursive_merge([logger_, __TMP__]...)  # 'cause it is not in-place.
            else
                if haskey(logger_, $symbol)
                    logger_ = $setindex!!(logger_, $recursive_merge([logger_[$symbol], __TMP__]...), $symbol)
                else
                    logger_ = $setindex!!(logger_, __TMP__, $symbol)
                end
            end
            $(__LOGGER__) = logger_
            $expr
        end
        esc(res)
    else
        error_invalid_expr_head()
    end
end

macro nested_log(symbol, expr)
    esc(:(@isdefined($:__LOGGER__) ? @nested_log($:__LOGGER__, $symbol, $expr) : $expr))
end

macro nested_log(expr)
    if expr isa Symbol  # @nested_log x
        res = quote
            @nested_log(esc($expr), $expr)  # __LOGGER__[:x] = x
        end
        esc(res)
    elseif expr.head == :call  # @nested_log my_func() where my_func() |> typeof :< NamedTuple
        _expr = copy(expr)
        push!(_expr.args, :(__LOG_INDICATOR__()))  # to distinguish whether call extended function with __LOGGER__ or not (in advance)
        res = quote
            if @isdefined($:__LOGGER__)
                __LOGGER__ = $recursive_merge([__LOGGER__, $_expr]...)
                # $expr
            # else
            #     $expr
            end
            $expr
        end
        esc(res)
    else
        error_invalid_expr_head()
    end
end

"""
    @nested_onlylog(symbol, expr)

A macro that activates given expression (`expr`) only when logging data.
Unlike `@nested_log(symbol, expr)`,
this macro does not evaluate given experssion `expr`.
"""
macro nested_onlylog(symbol, expr)
    esc(:(@isdefined($:__LOGGER__) ? @nested_log($:__LOGGER__, $symbol, $expr) : nothing))
end

"""
    @nested_onlylog(expr)

A macro that activates given expression (`expr`) only when logging data.
Unlike `@nested_log(expr)`,
this macro does not evaluate given experssion `expr`.
"""
macro nested_onlylog(expr)
    esc(:(@isdefined($:__LOGGER__) ? @nested_log($:__LOGGER__, nothing, $expr) : nothing))
end

# etc
function error_invalid_expr_head()
    error("Invalid expression head")
end


text:
function value(obj::ExtendedUnivariateDifferentiableFunction, arg0::Derivative)
    return jcall(obj, "value", Derivative, (Derivative,), arg0)
end

function value(obj::UnivariateDifferentiableFunction, arg0::DerivativeStructure)
    return jcall(obj, "value", DerivativeStructure, (DerivativeStructure,), arg0)
end

function value(obj::UnivariateFunction, arg0::jdouble)
    return jcall(obj, "value", jdouble, (jdouble,), arg0)
end


text:
# Tests formulations in prob/scopf
@testset "test scopf" begin

scopf_dc_cuts_soft_woc_objective = [14642.16, 26982.17]
@testset "scopf cuts dc soft 2, without cuts - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)
    network["gen_flow_cuts"] = []
    network["branch_flow_cuts"] = []

    result = run_scopf_cuts_soft(network, DCPPowerModel, lp_solver)

    @test isapprox(result["termination_status"], OPTIMAL)
    @test isapprox(result["objective"], scopf_dc_cuts_soft_woc_objective[i]; atol = 1e0)
end

@testset "scopf cuts dc soft 2, infeasible" begin
    network = deepcopy(network_infeasible)
    network["gen_flow_cuts"] = []
    network["branch_flow_cuts"] = []

    result = run_scopf_cuts_soft(network, DCPPowerModel, lp_solver)
    @test isapprox(result["termination_status"], INFEASIBLE)
end

scopf_dc_cuts_soft_wc_objective = [14642.16, 37403.76]
@testset "scopf cuts dc soft 2, with cuts - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)
    network["gen_flow_cuts"] = []
    network["branch_flow_cuts"] = []

    result = run_opf_cheap(network, DCPPowerModel, lp_solver)
    @test isapprox(result["termination_status"], OPTIMAL)

    update_active_power_data!(network, result["solution"])

    cuts = check_contingencies_branch_power(network, total_cut_limit=2, gen_flow_cuts=[], branch_flow_cuts=[])

    #println(length(cuts.gen_cuts) + length(cuts.branch_cuts))
    #cuts_found = sum(length(c.gen_cuts)+length(c.branch_cuts) for c in cuts)
    append!(network["gen_flow_cuts"], cuts.gen_cuts)
    append!(network["branch_flow_cuts"], cuts.branch_cuts)

    result = run_scopf_cuts_soft(network, DCPPowerModel, lp_solver)

    @test isapprox(result["termination_status"], OPTIMAL)
    @test isapprox(result["objective"], scopf_dc_cuts_soft_wc_objective[i]; atol = 1e0)
end

scopf_ac_cuts_soft_wc_objective = [14676.95, 37915.33]
@testset "scopf cuts ac soft 2, with cuts - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)
    network["gen_flow_cuts"] = []
    network["branch_flow_cuts"] = []

    result = run_opf_cheap(network, ACPPowerModel, nlp_solver)
    @test isapprox(result["termination_status"], LOCALLY_SOLVED)

    update_active_power_data!(network, result["solution"])

    cuts = check_contingencies_branch_power(network, total_cut_limit=2, gen_flow_cuts=[], branch_flow_cuts=[])

    #println(length(cuts.gen_cuts) + length(cuts.branch_cuts))
    #cuts_found = sum(length(c.gen_cuts)+length(c.branch_cuts) for c in cuts)
    append!(network["gen_flow_cuts"], cuts.gen_cuts)
    append!(network["branch_flow_cuts"], cuts.branch_cuts)

    result = run_scopf_cuts_soft(network, ACPPowerModel, nlp_solver)

    @test isapprox(result["termination_status"], LOCALLY_SOLVED)
    @test isapprox(result["objective"], scopf_ac_cuts_soft_wc_objective[i]; atol = 1e0)
end

scopf_dc_cuts_soft_woc_objective = [14642.16, 26982.17]
@testset "scopf cuts dc soft 2, without cuts - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)
    network["gen_flow_cuts"] = []
    network["branch_flow_cuts"] = []

    result = run_scopf_cuts_soft_bpv(network, DCPPowerModel, lp_solver)

    @test isapprox(result["termination_status"], OPTIMAL)
    @test isapprox(result["objective"], scopf_dc_cuts_soft_woc_objective[i]; atol = 1e0)
end

@testset "scopf cuts dc soft 2, infeasible" begin
    network = deepcopy(network_infeasible)
    network["gen_flow_cuts"] = []
    network["branch_flow_cuts"] = []

    result = run_scopf_cuts_soft_bpv(network, DCPPowerModel, lp_solver)
    @test isapprox(result["termination_status"], INFEASIBLE)
end

scopf_dc_cuts_soft_wc_objective = [14642.16, 30737.94]
@testset "scopf cuts dc soft 2, with cuts - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)
    network["gen_flow_cuts"] = []
    network["branch_flow_cuts"] = []

    result = run_opf_cheap(network, DCPPowerModel, lp_solver)
    @test isapprox(result["termination_status"], OPTIMAL)

    update_active_power_data!(network, result["solution"])

    cuts = check_contingencies_branch_power_bpv(network, total_cut_limit=2, gen_flow_cuts=[], branch_flow_cuts=[])

    #println(length(cuts.gen_cuts) + length(cuts.branch_cuts))
    #cuts_found = sum(length(c.gen_cuts)+length(c.branch_cuts) for c in cuts)
    append!(network["gen_flow_cuts"], cuts.gen_cuts)
    append!(network["branch_flow_cuts"], cuts.branch_cuts)

    result = run_scopf_cuts_soft_bpv(network, DCPPowerModel, lp_solver)

    @test isapprox(result["termination_status"], OPTIMAL)
    @test isapprox(result["objective"], scopf_dc_cuts_soft_wc_objective[i]; atol = 1e0)
end

scopf_ac_cuts_soft_wc_objective = [14676.95, 28346.36]
@testset "scopf cuts ac soft 2, with cuts - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)
    network["gen_flow_cuts"] = []
    network["branch_flow_cuts"] = []

    result = run_opf_cheap(network, ACPPowerModel, nlp_solver)
    @test isapprox(result["termination_status"], LOCALLY_SOLVED)

    update_active_power_data!(network, result["solution"])

    cuts = check_contingencies_branch_power_bpv(network, total_cut_limit=2, gen_flow_cuts=[], branch_flow_cuts=[])

    #println(length(cuts.gen_cuts) + length(cuts.branch_cuts))
    #cuts_found = sum(length(c.gen_cuts)+length(c.branch_cuts) for c in cuts)
    append!(network["gen_flow_cuts"], cuts.gen_cuts)
    append!(network["branch_flow_cuts"], cuts.branch_cuts)

    result = run_scopf_cuts_soft_bpv(network, ACPPowerModel, nlp_solver)

    @test isapprox(result["termination_status"], LOCALLY_SOLVED)
    @test isapprox(result["objective"], scopf_ac_cuts_soft_wc_objective[i]; atol = 1e0)
end

end # close test group


text:
#===============================================================================

    Abstract types reflecting conceptual inheritence structure that come
    in addition to those defined in StochasticProcessMCMC.jl

===============================================================================#
#NOTE possibly superfluous
"""
    DiffusionSpecificCallback <: Callback

Supertype of all callbacks specific to diffusion processes.
"""
abstract type DiffusionSpecificCallback <: eMCMC.Callback end

"""
    MCMCDiffusionImputation <: MCMCImputation

Supertype of all imputations related to diffusions.
"""
abstract type MCMCDiffusionImputation <: eMCMC.MCMCImputation end

"""
    BlockingType <: MCMCUpdateDecorator

Supertype of all blocking patterns.
"""
abstract type BlockingType <: eMCMC.MCMCUpdateDecorator end

"""
    GenericMCMCBackend

A flag that no specific backend is passed.
"""
struct DiffusionMCMCBackend <: eMCMC.MCMCBackend end


text:
#=
test_parser:
- Julia version: 1.1.0
- Author: bramb
- Date: 2019-02-13
=#
using Test

function has_no_errors(x)
    isempty(filter(x -> x.type == SimpleTAGMLTokenizer.ERROR, x))
end

@testset "tagml use cases" begin

@testset "use case 1: 2 layers, full overlap" begin
    using HyperGraphTools.SimpleTAGMLTokenizer
    using HyperGraphTools.CYK

    # use case 1: 2 layers, full overlap
    tagml1 = "[root>[tag|a>[other|b>one two three<other]<tag]<root]"
#     grammar1 = """
#     S -> ROOT
#     ROOT -> [root> TAG <root]
#     TAG -> [tag|a> OTHER <tag]
#     OTHER -> [other|b> TEXT <other]
#     TEXT -> "one two three"
#     """
#     chomsky_grammar1 = """
#     S -> ROOT
#     ROOT -> ROOT_START TAG_ROOT_END
#     TAG_ROOT_END -> TAG ROOT_END
#     ROOT_START -> [root>
#     ROOT_END -> <root]
#     TAG -> TAG_START OTHER_TAG_END
#     OTHER_TAG_END -> OTHER TAG_END
#     TAG_START -> [tag|a>
#     TAG_END -> <tag]
#     OTHER -> OTHER_START TEXT_OTHER_END
#     TEXT_OTHER_END -> TEXT OTHER_END
#     OTHER_START -> [other|b>
#     OTHER_END -> <other]
#     TEXT -> "one two three"
#     """
    tokens = tokenize(tagml1)
    @test has_no_errors(tokens)
    string_tokens = [t.content for t in tokens]
#     println("string_tokens=$string_tokens")
    grammar = Grammar()
    add_rule!(grammar, Rule("S",("ROOT_START","TAG_ROOT_END")))
    add_rule!(grammar, Rule("TAG_ROOT_END",("TAG","ROOT_END")))
    add_rule!(grammar, Rule("TAG", ("TAG_START", "OTHER_TAG_END")))
    add_rule!(grammar, Rule("OTHER_TAG_END", ("OTHER", "TAG_END")))
    add_rule!(grammar, Rule("OTHER", ("OTHER_START", "TEXT_OTHER_END")))
    add_rule!(grammar, Rule("TEXT_OTHER_END", ("TEXT", "OTHER_END")))
    add_rule!(grammar, Rule("ROOT_START", "[root>"))
    add_rule!(grammar, Rule("TAG_START", "[tag|a>"))
    add_rule!(grammar, Rule("OTHER_START", "[other|b>"))
    add_rule!(grammar, Rule("TEXT", "one two three"))
    add_rule!(grammar, Rule("OTHER_END", "<other]"))
    add_rule!(grammar, Rule("TAG_END", "<tag]"))
    add_rule!(grammar, Rule("ROOT_END", "<root]"))
#     println("grammar=$grammar")
    @test validate(string_tokens,grammar)
end

@testset "use case 2: 2 layers, partial overlap" begin
    using HyperGraphTools.SimpleTAGMLTokenizer
    using HyperGraphTools.CYK

    # use case 2: 2 layers, partial overlap
    tagml2 = "[root>[tag|a>one [other|b>two<tag] three<other]<root]"
#     grammar2 = """
#     S -> ROOT
#     ROOT -> [root> ROOTBODY <root]
#     ROOTBODY -> OPENTAG ROOTBODY | CLOSETAG ROOTBODY | TEXT ROOTBODY | ε
#     OPENTAG -> [tag|a> | [other|b>
#     CLOSETAG -> <tag] | <other]
#     TEXT -> "one " | "two" | " three"
#     """
    tokens = tokenize(tagml2)
    @test has_no_errors(tokens)

    string_tokens = [t.content for t in tokens]
    grammar = Grammar()
    add_rule!(grammar, Rule("S",("ROOTSTART","MIXED_ROOTEND")))
    add_rule!(grammar, Rule("MIXED_ROOTEND",("MIXED","ROOTEND")))
    add_rule!(grammar, Rule("MIXED", ("MIXED", "MIXED")))
    add_rule!(grammar, Rule("ROOTSTART", "[root>"))
    add_rule!(grammar, Rule("MIXED", "[tag|a>"))
    add_rule!(grammar, Rule("MIXED", "[other|b>"))
    add_rule!(grammar, Rule("MIXED", "one "))
    add_rule!(grammar, Rule("MIXED", "two"))
    add_rule!(grammar, Rule("MIXED", " three"))
    add_rule!(grammar, Rule("MIXED", "<other]"))
    add_rule!(grammar, Rule("MIXED", "<tag]"))
    add_rule!(grammar, Rule("ROOTEND", "<root]"))
#     println("grammar=$grammar")
    @test validate(string_tokens,grammar)
end

@testset "use case 3: 2 layers, no overlap" begin
    using HyperGraphTools.SimpleTAGMLTokenizer
    using HyperGraphTools.CYK
    # use case 3: 2 layers, no overlap
    tagml3 = "[root>[tag|a>one<tag] two [other|b>three<other]<root]"
#     grammar3 = """
#     S -> ROOT
#     ROOT -> [root> ROOTBODY <root]
#     ROOTBODY -> (OPENTAG | CLOSETAG | TEXT)+
#     OPENTAG -> [tag|a> | [other|b>
#     CLOSETAG -> <tag] | <other]
#     TEXT -> "one" | " two " | "three"
#     """
#     grammar_idlp3 = """
#     S -> {ROOT}
#     ROOT -> {ROOTOPEN,BODY,ROOTCLOSE}
#     ROOTOPEN < {BODY,ROOTCLOSE}
#     BODY < ROOTCLOSE
#     BODY -> {OPENTAG,CLOSETAG,TEXT,BODY}
#     OPENTAG -> {OPENTAG("tag"),OPENTAG("other")}
#     CLOSETAG -> {CLOSETAG("tag"),CLOSETAG("other")}
#     OPENTAG("tag") < CLOSETAG("tag")
#     OPENTAG("other") < CLOSETAG("other")
#     """
    tokens = tokenize(tagml3)
    @test has_no_errors(tokens)

    string_tokens = [t.content for t in tokens]
    grammar = Grammar()
    add_rule!(grammar, Rule("S",("ROOTSTART","MIXED_ROOTEND")))
    add_rule!(grammar, Rule("MIXED_ROOTEND",("MIXED","ROOTEND")))
    add_rule!(grammar, Rule("MIXED", ("MIXED", "MIXED")))
    add_rule!(grammar, Rule("ROOTSTART", "[root>"))
    add_rule!(grammar, Rule("MIXED", "[tag|a>"))
    add_rule!(grammar, Rule("MIXED", "[other|b>"))
    add_rule!(grammar, Rule("MIXED", "one"))
    add_rule!(grammar, Rule("MIXED", " two "))
    add_rule!(grammar, Rule("MIXED", "three"))
    add_rule!(grammar, Rule("MIXED", "<other]"))
    add_rule!(grammar, Rule("MIXED", "<tag]"))
    add_rule!(grammar, Rule("ROOTEND", "<root]"))
#     println("grammar=$grammar")
    @test validate(string_tokens,grammar)
end

@testset "use case 4: 2 layers, shared tags (a)" begin
    using HyperGraphTools.SimpleTAGMLTokenizer
    using HyperGraphTools.CYK

    # use case 4: 2 layers, shared tags
    tagml4a = "[root|+a,+b>[tag|a,b>[other|a,b>one two three<other]<tag]<root]"
    tokens = tokenize(tagml4a)
    @test has_no_errors(tokens)

    string_tokens = [t.content for t in tokens]
    grammar = Grammar()
    add_rule!(grammar, Rule("S",("ROOTSTART","MIXED_ROOTEND")))
    add_rule!(grammar, Rule("MIXED_ROOTEND",("MIXED","ROOTEND")))
    add_rule!(grammar, Rule("MIXED", ("MIXED", "MIXED")))
    add_rule!(grammar, Rule("ROOTSTART", "[root|+a,+b>"))
    add_rule!(grammar, Rule("MIXED", "[tag|a,b>"))
    add_rule!(grammar, Rule("MIXED", "[other|a,b>"))
    add_rule!(grammar, Rule("MIXED", "one two three"))
    add_rule!(grammar, Rule("MIXED", "<other]"))
    add_rule!(grammar, Rule("MIXED", "<tag]"))
    add_rule!(grammar, Rule("ROOTEND", "<root]"))
#     println("grammar=$grammar")
    @test validate(string_tokens,grammar)
end

@testset "use case 4: 2 layers, shared tags (b)" begin
    using HyperGraphTools.SimpleTAGMLTokenizer
    using HyperGraphTools.CYK

    # use case 4: 2 layers, shared tags
    tagml4b = "[root|+a,+b>[tag|a,b>[other|a,b>[w|a>one<w] [b|b>two<b] three<other]<tag]<root]"
    tokens = tokenize(tagml4b)
    @test has_no_errors(tokens)

    string_tokens = [t.content for t in tokens]
    grammar = Grammar()
    add_rule!(grammar, Rule("S",("ROOTSTART","MIXED_ROOTEND")))
    add_rule!(grammar, Rule("MIXED_ROOTEND",("MIXED","ROOTEND")))
    add_rule!(grammar, Rule("MIXED", ("MIXED", "MIXED")))
    add_rule!(grammar, Rule("ROOTSTART", "[root|+a,+b>"))
    add_rule!(grammar, Rule("MIXED", "[tag|a,b>"))
    add_rule!(grammar, Rule("MIXED", "[other|a,b>"))
    add_rule!(grammar, Rule("MIXED", "[w|a>"))
    add_rule!(grammar, Rule("MIXED", "one"))
    add_rule!(grammar, Rule("MIXED", "<w]"))
    add_rule!(grammar, Rule("MIXED", " "))
    add_rule!(grammar, Rule("MIXED", "[b|b>"))
    add_rule!(grammar, Rule("MIXED", "two"))
    add_rule!(grammar, Rule("MIXED", "<b]"))
    add_rule!(grammar, Rule("MIXED", " three"))
    add_rule!(grammar, Rule("MIXED", "<other]"))
    add_rule!(grammar, Rule("MIXED", "<tag]"))
    add_rule!(grammar, Rule("ROOTEND", "<root]"))
#     println("grammar=$grammar")
    @test validate(string_tokens,grammar)
end

end

max_stars_repo_name: brambg/HyperGraphTools

max_issues_repo_name: brambg/HyperGraphTools

max_forks_repo_name: brambg/HyperGraphTools


text:
using IntervalWavelets
import IntervalWavelets: inner, unit
using Base.Test

println("Testing reconstruction in wavelet basis...")

#=
- Scaling functions related to wavelets with p vanishing moments should reconstruct polynomials of degree p
- Unit vectors should reconstruct a single function
=#

R = 7
J = 4
S = R - J
DS = DaubSupport(0,1)

# ------------------------------------------------------------------------
# 1D

for p = 2:8
	for k = 0:2^J-1
		u = unit(2^J, k+1)
		xrecon, yrecon = weval(u, p, R, DaubSupport(0,1))

		y = IntervalScaling(p, k, J, R)
		#@show k, norm(yrecon - y)
		@test_approx_eq yrecon y
	end
end

# ------------------------------------------------------------------------
# 2D

for p = 2:8
	p = 2
	# The tests pass for for all l and k, but takes forever w/o inlining
	k = rand( [0:2^J-1;] )
	l = rand( [0:2^J-1;] )

	u1 = unit(2^J, k+1)
	u2 = unit(2^J, l+1)
	U = u1*u2'
	yrecon = weval(U, p, R, DaubSupport(0,1))

	y1 = IntervalScaling(p, k, J, R)
	y2 = IntervalScaling(p, l, J, R)
	Y = y1*y2'
	#@show k, l, norm(yrecon - Y)
	@test_approx_eq yrecon Y
end


text:
using DiskJockey.EnsembleSampler

# Rosenbrock function
@everywhere function lnprob(p)
    x, y = p

    return -((1.0 - x)^2 + 100 * (y - x^2)^2)
end

# A 2D Gaussian to check against
# @everywhere function lnprob(p)
#
#     a, b = p
#
#     if a < 0.0
#         return -Inf
#     end
#
#     mu = [1.0, -2.0]
#     Sigma = [[0.5^2,  0.00] [0.00, 0.3^2]]
#
#     R = p - mu # Residual vector
#     ld = logdet(Sigma)
#
#     lnp = -0.5 * ((R' * (Sigma\R))[1] + ld + 2 * log(2 * pi))
#     return lnp
#
# end

nwalkers = 40
ndim = 2

sampler = Sampler(nwalkers, ndim, lnprob)

# pos0 is the starting position, it needs to be a (ndim, nwalkers array)
pos0 = Array(Float64, ndim, nwalkers)
for i=1:nwalkers
    pos0[:,i] = randn(ndim)
end

function f(sampler::Sampler, outdir::AbstractString)
    println("calling function")
    try
        spawn(`plotly_walkers.py --name test`)
    catch
        println("Couldn't reach plotly server.")
    end
end

run_schedule(sampler, pos0, 10, 1, "", f)

# Now, let's try sampling a function with more parameters than walkers (just for testing purposes)
# A 4D Gaussian to check against
@everywhere function lnprob(p)

    mu = [1.0, -2.0, 0.0, 3.4]
    Sigma = [[0.5^2,  0.00, 0.0, 0.0] [0.00, 0.3^2, 0.0, 0.0] [0.0, 0.0, 0.7^2, 0.0] [0.0, 0.0, 0.0, 1.0^2] ]

    R = p - mu # Residual vector
    ld = logdet(Sigma)

    lnp = -0.5 * ((R' * (Sigma\R))[1] + ld + 2 * log(2 * pi))
    return lnp

end

nwalkers = 4
ndim = 4

sampler = Sampler(nwalkers, ndim, lnprob, true)

# pos0 is the starting position, it needs to be a (ndim, nwalkers array)
pos0 = Array(Float64, ndim, nwalkers)
for i=1:nwalkers
    pos0[:,i] = randn(ndim)
end

run_schedule(sampler, pos0, 1, 1, "", f)

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{SVF}) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{SAF{T}}) where {T<:Real} = true

"""
set and get function overloads
"""
MOI.get(model::Optimizer, ::MOI.ObjectiveSense) = model.inner.sense

function MOI.set(model::Optimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    model.inner.sense = sense
    return
end

function MOI.set(model::Optimizer, ::MOI.ObjectiveFunction, func::SVF)
    check_inbounds(model, func)
    model.inner.var_in_obj[func.variable.value] = true
    model.inner.objective =
        SingleVariableObjective(func, func.variable.value, [func.variable.value])
    return
end

function MOI.set(model::Optimizer, ::MOI.ObjectiveFunction, func::SAF{T}) where {T<:Real}
    check_inbounds(model, func)
    indices = [func.terms[i].variable_index.value for i in 1:length(func.terms)]
    coeffs = [func.terms[i].coefficient for i in 1:length(func.terms)]
    lc = LinearCombination(indices, coeffs)
    model.inner.var_in_obj[indices] .= true
    model.inner.objective = LinearCombinationObjective(func, lc, func.constant, indices)
    return
end


text:
using GLVisualize, GeometryTypes, Colors, GLAbstraction

if !isdefined(:runtests)
    window = glscreen()
    timesignal = loop(linspace(0f0, 1f0, 360))
end

description = """
Animating text particles color, scale and rotation.
"""

# Did we mention, that text is just a normal sprite particle system?
# This means, it supports all the attributes like the other particle systems.

# some text again
s = """The quick brown fox jumped over
some lazy text sample.
He wasn't really into numbers, but it's
really important to try out number rendering:
This number goes from 0 to π in no time!
And then back to 0 again... Wow!
This is real crazy stuff,
but it gets even more ludicrous:
∮ E⋅da = Q,  n → ∞, ∑ f(i) = ∏ g(i),
∀x∈ℝ: ⌈x⌉ = −⌊−x⌋, α ∧ ¬β = ¬(¬α ∨ β),
ℕ ⊆ ℕ₀ ⊂ ℤ ⊂ ℚ ⊂ ℝ ⊂ ℂ,
⊥ < a ≠ b ≡ c ≤ d ≪ ⊤ ⇒ (A ⇔ B),
2H₂ + O₂ ⇌ 2H₂O, R = 4.7 kΩ, ⌀
I can't even...
"""

# create a rotation from the time signal
rotation = map(timesignal) do t
    t2π = t*pi*2
    Vec3f0(cos(t2π),sin(t2π), 1)
end

# create some funcy scale change
scale = map(timesignal) do t
    circular = sin(t*pi)
    (1 + circular * 0.5) * 20f0
end
const len = length(s)

# per glyph color
color = map(timesignal) do t
    RGBA{Float32}[RGBA{Float32}(i/len,1, (sin(t)+1.)/2., 1) for i=1:len]
end

# per glyph stroke color
stroke_color = RGBA{Float32}[RGBA{Float32}(0,0,0, i/len) for i=1:len]

# _view and visualize it!
# you could also pass positions as a keyword argument or make
# the scale/rotation per glyph by supplying a Vector of them.
textvizz = visualize(s,
    model = translationmatrix(Vec3f0(0,600,0)), # move this up, since the text starts at 0 and goes down from there
    rotation = rotation,
    color = color,
    stroke_color = stroke_color,
    stroke_width = 1f0,
    relative_scale = scale
)

_view(textvizz, window)

if !isdefined(:runtests)
    renderloop(window)
end


text:
type JToc{T1<:AbstractString}
  curve_names::Vector{T1}
  multimesh_names::Vector{T1}
  multimeshadj_names::Vector{T1}
  multivar_names::Vector{T1}
  multimat_names::Vector{T1}
  multimatspecies_names::Vector{T1}
  csgmesh_names::Vector{T1}
  csgvar_names::Vector{T1}
  defvars_names::Vector{T1}
  qmesh_names::Vector{T1}
  qvar_names::Vector{T1}
  ucdmesh_names::Vector{T1}
  ucdvar_names::Vector{T1}
  ptmesh_names::Vector{T1}
  ptvar_names::Vector{T1}
  mat_names::Vector{T1}
  matspecies_names::Vector{T1}
  var_names::Vector{T1}
  obj_names::Vector{T1}
  dir_names::Vector{T1}
  array_names::Vector{T1}
  mrgtree_names::Vector{T1}
  groupelmap_names::Vector{T1}
  mrgvar_names::Vector{T1}

  namemapping::Dict{Symbol,Symbol}
  function JToc(dbfile::Ptr{DBfile})
    inttoc = new()
    const thesevecnames = @compat Dict(:ncurve => :curve_names,:nmultimesh=>:multimesh_names,:nmultimeshadj=>:multimeshadj_names,:nmultivar=>:multivar_names,:nmultimat=>:multimat_names,:nmultimatspecies=>:multimatspecies_names,:ncsgmesh=>:csgmesh_names,:ncsgvar=>:csgvar_names,:ndefvars=>:defvars_names, :nqmesh=>:qmesh_names,:nqvar=>:qvar_names,:nucdmesh=>:ucdmesh_names,:nucdvar=>:ucdvar_names,:nptmesh=>:ptmesh_names,:nptvar=>:ptvar_names,:nmat=>:mat_names, :nmatspecies=>:matspecies_names,:nvar=>:var_names,:nobj=>:obj_names,:ndir=>:dir_names,:narray=>:array_names,:nmrgtree=>:mrgtree_names,:ngroupelmap=>:groupelmap_names,:nmrgvar=>:mrgvar_names)

    const rawnames = unsafe_load(DBGetToc(dbfile))
    const thesenums = [keyname => Int32(rawnames.(keyname)) for keyname in keys(thesevecnames)]

    for keyname in keys(thesevecnames)
      if thesenums[keyname] == 0
        inttoc.(thesevecnames[keyname]) = Array(ASCIIString, 0,)
      else
        inttoc.(thesevecnames[keyname]) = map(bytestring,pointer_to_array(rawnames.(thesevecnames[keyname]), thesenums[keyname]))
      end
  end
    inttoc.namemapping = thesevecnames
  return inttoc
  end #function
end #type


text:
################################################################################
#
#             AbGrp.jl : Finitely generated abelian groups
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015, 2016 Tommy Hofmann
#
################################################################################

export AbGrp, AbGrpElem, parent, is_finite, is_infinite, rank, length,
       getindex, show, +, *, call

################################################################################
#
#  Parent and element types
#
################################################################################

# At the moment support only for diagonal relation matrices
# We hash the sorted diagonal

const AbGrpID = Dict{Array{fmpz, 1}, Nemo.Set}()

type AbGrp <: Nemo.Set
  diagonal::Array{fmpz, 1}
  inf_num::Int
  fin_num::Int

  function AbGrp(a::Array{fmpz, 1})
    sort!(a)
    #println(a, " ",hash(a))
    if haskey(AbGrpID, a)
      return AbGrpID[a]::AbGrp
    else
      z = new()
      z.diagonal = a
      z.inf_num = findfirst(x -> x != 0, a) - 1
      z.fin_num = length(a) - z.inf_num
      AbGrpID[a] = z
      return z
    end
  end

  function AbGrp(a::Array{Int, 1})
    return AbGrp(map(fmpz, a))
  end
end

type AbGrpElem
  parent::AbGrp
  coord::Array{fmpz, 1}

  function AbGrpElem(A::AbGrp, a::Array{fmpz, 1})
    z = new()
    z.parent = A
    z.coord = a
    return z
  end
end

################################################################################
#
#  Field access
#
################################################################################

doc"""
***
    parent(x::AbGrpElem) -> AbGrp

> Returns the parent of $x$.
"""
function parent(x::AbGrpElem)
  return x.parent
end

################################################################################
#
#  Predicates and basic invariants
#
################################################################################

doc"""
***
    is_finite(A::AbGrp) -> Bool

> Returns whether $A$ is finite.
"""
is_finite(A::AbGrp) = A.inf_num == 0 ? true : false

doc"""
***
    is_infinite(A::AbGrp) -> Bool

> Returns whether $A$ is infinite.
"""
is_infinite(A::AbGrp) = !is_finite(A)

doc"""
***
    rank(A::AbGrp) -> Int

> Returns the rank of $A$, that is, the dimension of the
> $\mathbf{Q}$-vectorspace $A \otimes_{\mathbf Z} \mathbf Q$.
"""
rank(A::AbGrp) = A.inf_num

doc"""
***
    length(A::AbGrp) -> Int

> ??
"""
length(A::AbGrp) = length(A.diagonal)

doc"""
***
    getindex(x::AbGrpElem, i::Int) -> fmpz

> Returns the $i$-th component of the element $x$.
"""
function getindex(x::AbGrpElem, i::Int)
  return x.coord[i]
end

doc"""
***
    order(A::AbGrpElem, i::Int) -> fmpz

> Returns the order of $A$. It is assumed that $A$ is not infinite.
"""
function order(A::AbGrp)
  is_infinite(A) && error("Group must be finite")
  return prod(A.diagonal)
end

################################################################################
#
#  I/O
#
################################################################################

function show(io::IO, A::AbGrp)
  print(io, "Abelian group\n")

  if A.inf_num > 0
    print(io, "Z")
  end

  for i in 1:A.inf_num - 1
    print(io, " x Z")
  end

  if A.fin_num > 0 && A.inf_num > 0
    print(io, " x ")
  end

  if A.fin_num > 0
    print(io, "Z/$(A.diagonal[A.inf_num + 1])")
  end

  for i in 2:A.fin_num
    print(io, " x Z/$(A.diagonal[A.inf_num + i])")
  end
end

function show(io::IO, a::AbGrpElem)
  print(io, "Element of\n$(a.parent)\n with components\n$(a.coord)")
end

################################################################################
#
#  Addition
#
################################################################################

doc"""
***
    +(x::AbGrpElem, y::AbGrpElem) -> AbGrpElem

> Computes $x + y$.
"""
function +(x::AbGrpElem, y::AbGrpElem)
  parent(x) != parent(y) && error("Parents must coincide")
  z = Array(fmpz, length(x.coord))
  for i in 1:length(z)
    z[i] = x.coord[i] + y.coord[i]
  end
  return parent(x)(z) # this gives the reduced element
end

################################################################################
#
#  Scalar multiplication
#
################################################################################

doc"""
***
    *(x::Integer, y::AbGrpElem) -> AbGrpElem

> Computes $x \cdot y$.
"""
function *(x::Integer, y::AbGrpElem)
  z = parent(x)()
  for i in 1:length(z)
    z.coord[i] = x*y.coord[i]
  end
  return z
end

*(x::AbGrpElem, y::Integer) = y*x

doc"""
***
    *(x::fmpz, y::AbGrpElem) -> AbGrpElem

> Computes $x \cdot y$.
"""
function *(x::fmpz, y::AbGrpElem)
  z = parent(x)()
  for i in 1:length(z)
    z.coord[i] = x*y.coord[i]
  end
  return z
end

*(x::AbGrpElem, y::fmpz) = y*x

################################################################################
#
#  Parent object overloading
#
################################################################################

doc"""
***
    call(A::AbGrp, x::Array{fmpz, 1}) -> AbGrpElem

> Given an array of elements of type `fmpz` of the same length as $A$, this
> function returns the element of $A$ with components `x`.
"""
function call(A::AbGrp, x::Array{fmpz, 1})
  length(A) != length(x) && error("Lengths do not coincide")
  y = deepcopy(x)
  for i in A.inf_num + 1:length(y)
    y[i] = mod(x[i], A.diagonal[i])
  end
  z = AbGrpElem(A, y)
  return z
end

doc"""
***
    call{T <: Integer}(A::AbGrp, x::Array{T, 1}) -> AbGrpElem

> Given an array of elements of type `Integer` of the same length as $A$, this
> function returns the element of $A$ with components `x`.
"""
function call{T <: Integer}(A::AbGrp, x::Array{T, 1})
  length(A) != length(x) && error("Lengths do not coincide")
  z = A(map(fmpz, x))
  return z
end

doc"""
***
    getindex(A::AbGrp, i::Int) -> AbGrpElem

> Returns the element of $A$ with components $(0,\dotsc,0,1,0,\dotsc,0)$,
> where the $1$ is at the $i$-th position.
"""
function getindex(A::AbGrp, i::Int)
  z = Array(fmpz, length(A.diagonal))
  for j in 1:length(A.diagonal)
    z[j] = fmpz()
  end
  z[i] = fmpz(1)
  return AbGrpElem(A, z)
end


text:
n = 10 # rand(1:10)
# no null vector for 1-by-1 chebspec.
if n == 1
    n = 2
end
A = matrixdepot("chebspec", n)
# A has null vector ones(n)
@test A*ones(n) ≈ zeros(n) atol = 1e-6
B = matrixdepot("chebspec", n, 1)
# B's eigenvalues have negative real parts.
v = real(eigvals(B))
for i = 1:n
    @test v[i] < 0
end
println("'chebspec' passed test...")


text:
using DiffEqBayesTuring
using Test

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


text:
# Extra variables of a triangle, that can be computed from the Triangle
mutable struct TriangleAux
    # Middle of the triangle
    midpoint::Vector{Float64}
    # Normal
    normal::Vector{Float64}
    # Tangent vectors to the edges
    tau1::Vector{Float64}
    tau2::Vector{Float64}
    tau3::Vector{Float64}
    # Normal vectors for the triangles
    nu1::Vector{Float64}
    nu2::Vector{Float64}
    nu3::Vector{Float64}

    # For quadrature
    gaussP::Matrix{Float64} # Gauss points
    gaussW::Vector{Float64} # Gauss weights
end

struct Triangle #  Note: unmutable, whereas aux is mutable
    # The vertex coordinates
    P1::Vector{Float64}
    P2::Vector{Float64}
    P3::Vector{Float64}
    # Area of the element
    area::Float64
    # Auxilary info that can be computed from P1,P2,P3
    aux::TriangleAux
end

# Constructor
function Triangle(P1,P2,P3;area=0)
    aux=TriangleAux([],[],[],[],[],[],[],[],zeros(3,3),[]);
    return Triangle(P1,P2,P3,area,aux);
end


text:
# let
    q = Quat(1.0,2.0,3.0,4.5)
    @test length(q) == 4
    @test imag(q) == VecE3(1.0,2.0,3.0)
    @test conj(q) == Quat(-1.0,-2.0,-3.0,4.5)
    @test vec(q) == [1.0,2.0,3.0,4.5]
    @test vec(copy(q)) == vec(q)
    @test convert(Quat, [1.0,2.0,3.0,4.5]) == q
    @test norm(q) == norm(vec(q))
    @test vec(normalize(q)) == vec(q) ./ norm(vec(q))

    q = Quat(-0.002, -0.756, 0.252, -0.604)
    sol = push!(vec(get_axis(q)), get_rotation_angle(q))
    @test isapprox(sol, [-0.00250956, -0.948614,  0.316205, mod2pi(-1.84447)], atol=1e-3) ||
          isapprox(sol, [ 0.00250956,  0.948614, -0.316205, mod2pi( 1.84447)], atol=1e-3)

    srand(0)
    for i in 1 : 5
        a = normalize(VecE3(rand(),rand(),rand()))
        b = normalize(VecE3(rand(),rand(),rand()))
        q = quat_for_a2b(a,b)
        @test norm(rot(q, a) - b) < 1e-8
    end

    @test isapprox(angledist(Quat(1,0,0,0), Quat(0,1,0,0)), π, atol=1e-6)
    @test norm(lerp(Quat(1,0,0,0), Quat(0,1,0,0), 0.5) - Quat(√2/2, √2/2, 0, 0)) < 1e-6

    srand(0)
    for i in 1 : 5
        q = rand(Quat)
        @test isapprox(norm(q), 1.0, atol=1e-8)
        @test norm(inv(q)*q - Quat(0,0,0,1)) < 1e-8
        @test norm(q*inv(q) - Quat(0,0,0,1)) < 1e-8
    end
# end


text:
puzzleinput = readlines("inputs/day18.txt")

# First Half

struct My1
    val::Int
end
# Define MY numbers that follow different rules

import Base.+, Base.-

+(x::My1, y::My1) = My1(x.val + y.val)
# For MY numbers, '+' acts like the regular one

-(x::My1, y::My1) = My1(x.val * y.val)
# For MY numbers, '*' has the same precedence of '-' so, to mantain
# conventional precedence rules, the '*' operator is replaced by '-'

function my1(s)
    s = replace(s, r"(\d)" => s"My1(\1)")
    replace(s, r"\*" => s"-")
end
# Replace in the expression, numbers with MY numbers
# and operators with MY operators

evaluate1(s) = eval(Meta.parse(my1(s))).val
# After reaplacement evaluation is done with MY rules

sum(evaluate1.(puzzleinput))

# Second Half

# As the firs half, only rules change: now for MY numbers, '+' has the
# precedence rule of '*' and viceversa

import Base.*

struct My2
    val::Int
end

+(x::My2, y::My2) = My2(x.val * y.val)
*(x::My2, y::My2) = My2(x.val + y.val)

function my2(s)
    s = replace(s, r"(\d)" => s"My2(\1)")
    s = replace(s, r"\*" => s".")
    s = replace(s, r"\+" => s"*")
    replace(s, r"\." => s"+")
end

evaluate2(s) = eval(Meta.parse(my2(s))).val

sum(evaluate2.(puzzleinput))


text:
module FirstRound
using JuMP
using ..RocketlandDefns
using Mosek
using MathOptInterface
using MosekTools
using StaticArrays
using Rotations
import ..Dynamics
const MOI=MathOptInterface

#=
<<<<<<< HEAD
# JuMP-based initial 3dof solver. Doesn't work as a good initalizaion, though, and written for an old JuMP versions
=======
>>>>>>> 6d18d726748568198fe41d32ebcf2420ce96670c
function solve_initial(prob::DescentProblem)
    #computed constants
    kf = prob.K
    N = kf
    a = prob.alpha
    tfs = prob.tf_guess
    dt = tfs/kf
    mu = [((kf - k)/kf) * prob.mwet + (k/kf)*prob.mdry for k=0:kf]
    gamm0 = prob.Tmin
    g0 = prob.g
    mdry = prob.mdry
    tmin = prob.Tmin
    tmax = prob.Tmax
    thetamax = prob.thetaMax
    gammags = prob.gammaGs
    m0 = prob.mwet
    r0 = prob.rIi
    rf = prob.rIf
    v0 = prob.vIi
    vf = prob.vIf
    n0 = [1 0 0]
    nf = [1 0 0]
    nsc = 20
    wkar = 100.0

    tggs = tand(prob.gammaGs)
    sqcm = (cosd(prob.thetaMax))
    delMax = cosd(prob.deltaMax)

    m = Model(with_optimizer(MosekOptimizer,MSK_IPAR_INFEAS_REPORT_AUTO=1))

    @variable(m, T[1:3,1:N+1])
    @variable(m, s[1:N+1])
    @variable(m, r[1:3,1:N+1])
    @variable(m, v[1:3,1:N+1])
    @variable(m, ma[1:N+1])
    @variable(m, ga[1:N+1])
    @variable(m, kaR[1:N+1])
    @variable(m, kaRr[1:N+1])
    @variable(m, ar[1:3,1:N+1])
    @variable(m, nkaR)

    #initial pos
    @constraint(m, r[i=1:3,1] .== r0[i])
    @constraint(m, v[i=1:3,1] .== v0[i])
    @constraint(m, ma[1] .== m0)
    @constraint(m, r[1:3,N+1] .== 0)
    @constraint(m, v[1:3,N+1] .== 0)
    @constraint(m, 0 == T[2,N+1])
    @constraint(m, 0 == T[3,N+1])

    @objective(m, Min, -ma[N+1] + wkar*nkaR)
    @constraint(m, [nkaR, kaRr[:]...] in MOI.SecondOrderCone(N+2))
    @constraint(m, kaRr[:] .== kaR[:])

    for i=1:N
        @constraint(m, ma[i+1] == ma[i]- a*(ga[i] + ga[i+1])*dt/2.0)
        ak = T[:,i]/mu[i] + ar[:,i] + [-g0,0,0]
        akp = T[:,i+1]/mu[i+1] + ar[:,i+1] + [-g0,0,0]
        @constraint(m, r[j=1:3,i+1] .== r[j,i] + v[j,i]*dt + 1/3 * (ak[j] + 1/2*akp[j])*dt^2)
        @constraint(m, v[j=1:3,i+1] .== v[j,i] + 1/2 * (ak[j] + akp[j])*dt)
    end

    for i=1:N+1
        @constraint(m, mdry <= ma[i])
        @constraint(m, [r[1,i]/tggs, r[2,i], r[3,i]] in MOI.SecondOrderCone(3))
        @constraint(m, [ga[i], T[:,i]...] in MOI.SecondOrderCone(4))
        @constraint(m, tmin <= ga[i])
        @constraint(m, ga[i] <= tmax)
        @constraint(m, ga[i]*sqcm <= T[1,i])
        @constraint(m, [kaR[i], ar[:,i]...] in MOI.SecondOrderCone(4))
    end
    JuMP.optimize(m)

    mass = JuMP.resultvalue.(ma)
    ri = JuMP.resultvalue.(r)
    vi = JuMP.resultvalue.(v)
    Ti = JuMP.resultvalue.(T)

    initial_points = Array{LinPoint,1}(N+1)
    for k=1:N+1
        mk = mass[k]
        rIk = ri[:,k]
        vIk = vi[:,k]

        rot = rotation_between([1,0,0], -Ti[:,k])
        qBIk = @SVector [rot.w, rot.x, rot.y, rot.z]
        state_init = vcat(mk,rIk,vIk,qBIk,(@SVector [0.0,0,0]))
        control_init = @SVector [norm(Ti[:,k]),0,0]
        initial_points[k] = LinPoint(state_init, control_init)
    end

    return initial_points,Dynamics.linearize_dynamics(initial_points, prob.tf_guess, 1.0/(N+1), ProbInfo(prob))
end
=#

function linear_points(problem::DescentProblem)
    K = problem.K
    initial_points = Array{LinPoint,1}(UndefInitializer(), K+1)
    for k=0:K
        mk = (K-k)/(K) * problem.mwet + (k/(K))*problem.mdry
        rIk = (K-k)/(K) * problem.rIi + (k/(K))*problem.rIf
        vIk = (K-k)/(K) * problem.vIi + (k/(K))*problem.vIf

        rot = rotation_between([1,0,0], -vIk)
        qBIk = [rot.w, rot.x, rot.y, rot.z]
        TBk = [mk*problem.g,0,0]
        state_init = vcat(mk,rIk,vIk,qBIk,([0.0,0,0]))
        control_init = [mk*problem.g,0,0]
        initial_points[k+1] = LinPoint(state_init, control_init)
    end
    return initial_points
end

function linear_initial(problem::DescentProblem, cache::Dynamics.IntegratorCache)
    initial_points = linear_points(problem)
    return initial_points, #Dynamics.linearize_dynamics(initial_points, problem.tf_guess, 1/(problem.K+1), ProbInfo(problem))
        Dynamics.linearize_dynamics(initial_points, problem.tf_guess, 1/(problem.K+1), cache)
end

end

max_stars_repo_name: BenChung/SuccessiveConvexification

max_issues_repo_name: BenChung/SuccessiveConvexification

max_forks_repo_name: BenChung/SuccessiveConvexification


text:
module adapters_mysql_placeholder_test

using Test # @test_throws
using Octo.Adapters.MySQL # from to_sql Schema.model PlaceHolder SELECT FROM WHERE
using Octo.Backends: UnsupportedError

Repo.debug_sql()

include("options.jl")
Repo.connect(;
    adapter = Octo.Adapters.MySQL,
    Options.arguments...
)

struct User
end
Schema.model(User, table_name="users")

Repo.execute([DROP TABLE IF EXISTS User])
Repo.execute(Raw("""CREATE TABLE IF NOT EXISTS users (
                     ID INT NOT NULL AUTO_INCREMENT,
                     name VARCHAR(255),
                     salary FLOAT(8),
                     PRIMARY KEY (ID)
                 )"""))

❔ = Octo.PlaceHolder
u = from(User)

userName = """ ' OR '1'='1 """
@test_throws UnsupportedError Repo.query([SELECT * FROM u WHERE u.name == ❔ AND u.salary > ❔], [userName, 2000])

Repo.disconnect()

end # module adapters_mysql_placeholder_test


text:
module operators_lazyproduct

export LazyProduct

import Base: ==, *, /, +, -
import ..operators

using ..bases, ..states, ..operators, ..operators_dense
using SparseArrays

"""
    LazyProduct(operators[, factor=1])
    LazyProduct(op1, op2...)

Lazy evaluation of products of operators.

The factors of the product are stored in the `operators` field. Additionally a
complex factor is stored in the `factor` field which allows for fast
multiplication with numbers.
"""
mutable struct LazyProduct{BL<:Basis,BR<:Basis,T<:Tuple{Vararg{AbstractOperator}}} <: AbstractOperator{BL,BR}
    basis_l::BL
    basis_r::BR
    factor::ComplexF64
    operators::T

    function LazyProduct{BL,BR,T}(operators::T, factor::Number=1) where {BL<:Basis,BR<:Basis,T<:Tuple{Vararg{AbstractOperator}}}
        for i = 2:length(operators)
            check_multiplicable(operators[i-1], operators[i])
        end
        new(operators[1].basis_l, operators[end].basis_r, factor, operators)
    end
end
function LazyProduct(operators::T, factor::Number=1) where {T<:Tuple{Vararg{AbstractOperator}}}
    BL = typeof(operators[1].basis_l)
    BR = typeof(operators[end].basis_r)
    LazyProduct{BL,BR,T}(operators, factor)
end
LazyProduct(operators::Vector{T}, factor::Number=1) where T<:AbstractOperator = LazyProduct((operators...,), factor)
LazyProduct(operators::AbstractOperator...) = LazyProduct((operators...,))
LazyProduct() = throw(ArgumentError("LazyProduct needs at least one operator!"))

Base.copy(x::T) where T<:LazyProduct = T(([copy(op) for op in x.operators]...,), x.factor)

operators.dense(op::LazyProduct) = op.factor*prod(dense.(op.operators))
operators.dense(op::LazyProduct{B1,B2,T}) where {B1<:Basis,B2<:Basis,T<:Tuple{AbstractOperator}} = op.factor*dense(op.operators[1])
SparseArrays.sparse(op::LazyProduct) = op.factor*prod(sparse.(op.operators))
SparseArrays.sparse(op::LazyProduct{B1,B2,T}) where {B1<:Basis,B2<:Basis,T<:Tuple{AbstractOperator}} = op.factor*sparse(op.operators[1])

==(x::LazyProduct, y::LazyProduct) = false
==(x::T, y::T) where T<:LazyProduct = (x.operators==y.operators && x.factor == y.factor)

# Arithmetic operations
-(a::T) where T<:LazyProduct = T(a.operators, -a.factor)

*(a::LazyProduct{B1,B2}, b::LazyProduct{B2,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis} = LazyProduct((a.operators..., b.operators...), a.factor*b.factor)
*(a::T, b::Number) where T<:LazyProduct = T(a.operators, a.factor*b)
*(a::Number, b::T) where T<:LazyProduct = T(b.operators, a*b.factor)

/(a::LazyProduct, b::Number) = LazyProduct(a.operators, a.factor/b)

operators.dagger(op::LazyProduct) = LazyProduct(dagger.(reverse(op.operators)), conj(op.factor))

operators.tr(op::LazyProduct) = throw(ArgumentError("Trace of LazyProduct is not defined. Try to convert to another operator type first with e.g. dense() or sparse()."))

operators.ptrace(op::LazyProduct, indices::Vector{Int}) = throw(ArgumentError("Partial trace of LazyProduct is not defined. Try to convert to another operator type first with e.g. dense() or sparse()."))

operators.permutesystems(op::LazyProduct, perm::Vector{Int}) = LazyProduct(([permutesystems(op_i, perm) for op_i in op.operators]...,), op.factor)

operators.identityoperator(::Type{LazyProduct}, b1::Basis, b2::Basis) = LazyProduct(identityoperator(b1, b2))

# Fast in-place multiplication
function operators.gemv!(alpha, a::LazyProduct{B1,B2}, b::Ket{B2}, beta, result::Ket{B1}) where {B1<:Basis,B2<:Basis}
    tmp1 = Ket(a.operators[end].basis_l)
    operators.gemv!(a.factor, a.operators[end], b, 0, tmp1)
    for i=length(a.operators)-1:-1:2
        tmp2 = Ket(a.operators[i].basis_l)
        operators.gemv!(1, a.operators[i], tmp1, 0, tmp2)
        tmp1 = tmp2
    end
    operators.gemv!(alpha, a.operators[1], tmp1, beta, result)
end

function operators.gemv!(alpha, a::Bra{B1}, b::LazyProduct{B1,B2}, beta, result::Bra{B2}) where {B1<:Basis,B2<:Basis}
    tmp1 = Bra(b.operators[1].basis_r)
    operators.gemv!(b.factor, a, b.operators[1], 0, tmp1)
    for i=2:length(b.operators)-1
        tmp2 = Bra(b.operators[i].basis_r)
        operators.gemv!(1, tmp1, b.operators[i], 0, tmp2)
        tmp1 = tmp2
    end
    operators.gemv!(alpha, tmp1, b.operators[end], beta, result)
end

end # module


text:
# # [Linear Time Continuous System with Controls](@id linear_continuous_controls)
#
# Now we will extend the [`previous example`](@ref linear_continuous) by adding some exegeneous control signals.
# As always, we will generate some data via `OrdinaryDiffEq.jl`

using DataDrivenDiffEq
using ModelingToolkit
using LinearAlgebra
using OrdinaryDiffEq
#md using Plots

A = [-0.9 0.2; 0.0 -0.2]
B = [0.0; 1.0]
u0 = [10.0; -10.0]
tspan = (0.0, 10.0)

f(u,p,t) = A*u .+ B .* sin(0.5*t)

sys = ODEProblem(f, u0, tspan)
sol = solve(sys, Tsit5(), saveat = 0.05);

# We will use the data provided by our problem, but add the control signal `U = sin(0.5*t)` to it.
X = Array(sol)
t = sol.t
control(u,p,t) = [sin(0.5*t)]
prob = ContinuousDataDrivenProblem(X, t, U = control)

# And plot the problems data.

#md plot(prob)

# Again, we will use `gDMD` to estimate the systems dynamics. Since we have a control signal
# defined in the problem, the algorithm will detect it automatically and use `gDMDc`:

res = solve(prob, DMDSVD(), digits = 1)
#md println(res)

# We see that the system has been recovered correctly, indicated by the small error and high AIC score of the result. We can confirm this by looking at the resulting [`Basis`](@ref)

#md
system = result(res)
#md println(system)

# And also plot the prediction of the recovered dynamics

#md plot(res)

# Again, we can have a look at the generator of the system, which is independent from the inputs.

generator(system)

# Sticking to the same procedure as earlier, we now use a linear sparse regression to solve the problem

@parameters t
@variables x[1:2](t) u[1:1](t)

basis = Basis([x; u], x, controls = u, independent_variable = t, name = :LinearBasis)
#md print(basis) #hide

# Note that we added a new variable `u[1](t)` as a control to both the equations and the basis constructor.
# Afterwards, we simply `solve` the already defined problem with our `Basis` and a `SparseOptimizer`

sparse_res = solve(prob, basis, STLSQ(1e-1))
#md println(sparse_res)

# Which holds the same equations
sparse_system = result(sparse_res)
#md println(sparse_system)

# Again, we can have a look at the result

#md plot(
#md     plot(prob), plot(sparse_res), layout = (1,2)
#md )

# Both results can be converted into an `ODESystem`. To include the control signal, we simply
# substitute the control variables in the corresponding equations.

subs_control = (u[1] => sin(0.5*t))

eqs = map(equations(sparse_system)) do eq
    eq.lhs ~ substitute(eq.rhs, subs_control)
end

@named sys = ODESystem(
    eqs,
    get_iv(sparse_system),
    states(sparse_system),
    parameters(sparse_system)
    );

# And simulated using `OrdinaryDiffEq.jl` using the (known) initial conditions and the parameter mapping of the estimation.

x0 = [x[1] => u0[1], x[2] => u0[2]]
ps = parameter_map(sparse_res)

ode_prob = ODEProblem(sys, x0, tspan, ps)
estimate = solve(ode_prob, Tsit5(), saveat = prob.t);

# And look at the result
#md plot(sol, color = :black)
#md plot!(estimate, color = :red, linestyle = :dash)

#md # ## [Copy-Pasteable Code](@id linear_continuous_copy_paste)
#md #
#md # ```julia
#md # @__CODE__
#md # ```

@test all(aic(sparse_res) .>= 1e3) #src
@test all(l2error(sparse_res) .<= 5e-1) #src
@test all(determination(sparse_res) .>= 0.97) #src
@test Array(sol) ≈ Array(estimate) rtol = 5e-2 #src


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

"""
    Ngon(p1, p2, ..., pN)

A N-gon is a polygon with `N` vertices `p1`, `p2`, ..., `pN`
oriented counter-clockwise (CCW). In this case the number of
vertices is fixed and known at compile time. Examples of N-gon
are `Triangle` (N=3), `Quadrangle` (N=4), `Pentagon` (N=5), etc.

### Notes

- Although the number of vertices `N` is known at compile time,
  we use abstract vectors to store the list of vertices. This
  design allows constructing N-gon from views of global vectors
  without expensive memory allocations.

- Type aliases are `Triangle`, `Quadrangle`, `Pentagon`, `Hexagon`,
  `Heptagon`, `Octagon`, `Nonagon`, `Decagon`.
"""
struct Ngon{N,Dim,T,V<:AbstractVector{Point{Dim,T}}} <: Polygon{Dim,T}
  vertices::V
end

Ngon{N}(vertices::AbstractVector{Point{Dim,T}}) where {N,Dim,T} =
  Ngon{N,Dim,T,typeof(vertices)}(vertices)

Ngon(vertices::AbstractVector{Point{Dim,T}}) where {Dim,T} =
  Ngon{length(vertices)}(vertices)

# type aliases for convenience
const Triangle   = Ngon{3}
const Quadrangle = Ngon{4}
const Pentagon   = Ngon{5}
const Hexagon    = Ngon{6}
const Heptagon   = Ngon{7}
const Octagon    = Ngon{8}
const Nonagon    = Ngon{9}
const Decagon    = Ngon{10}

issimple(::Type{<:Ngon}) = true

nvertices(::Type{<:Ngon{N}}) where {N} = N
nvertices(ngon::Ngon) = nvertices(typeof(ngon))

# measure of N-gon embedded in 2D
function signarea(ngon::Ngon{N,2}) where {N}
  v = ngon.vertices
  sum(i -> signarea(v[1], v[i], v[i+1]), 2:N-1)
end
measure(ngon::Ngon{N,2}) where {N} = abs(signarea(ngon))

# measure of N-gon embedded in higher dimension
function measure(ngon::Ngon{N}) where {N}
  areaₜ(A, B, C) = norm((B - A) × (C - A)) / 2
  v = ngon.vertices
  sum(i -> areaₜ(v[1], v[i], v[i+1]), 2:N-1)
end

hasholes(::Ngon) = false

chains(ngon::Ngon{N}) where {N} = [Chain(ngon.vertices[[1:N; 1]])]

Base.unique!(ngon::Ngon) = ngon

function Base.in(p::Point{Dim,T}, ngon::Ngon{N,Dim,T}) where {N,Dim,T}
  # decompose n-gons into triangles by
  # fan triangulation (assumes convexity)
  # https://en.wikipedia.org/wiki/Fan_triangulation
  v = ngon.vertices
  Δ(i) = Triangle(view(v, [1,i,i+1]))
  any(i -> p ∈ Δ(i), 2:N-1)
end

# ----------
# TRIANGLES
# ----------

# triangles are special
issimplex(::Type{<:Triangle}) = true
isconvex(::Type{<:Triangle}) = true

function Base.in(p::Point{2,T}, t::Triangle{2,T}) where {T}
  # given coordinates
  a, b, c = t.vertices
  x₁, y₁ = coordinates(a)
  x₂, y₂ = coordinates(b)
  x₃, y₃ = coordinates(c)
  x , y  = coordinates(p)

  # barycentric coordinates
  λ₁ = ((y₂ - y₃)*(x  - x₃) + (x₃ - x₂)*(y  - y₃)) /
       ((y₂ - y₃)*(x₁ - x₃) + (x₃ - x₂)*(y₁ - y₃))
  λ₂ = ((y₃ - y₁)*(x  - x₃) + (x₁ - x₃)*(y  - y₃)) /
       ((y₂ - y₃)*(x₁ - x₃) + (x₃ - x₂)*(y₁ - y₃))
  λ₃ = one(T) - λ₁ - λ₂

  # barycentric check
  zero(T) ≤ λ₁ ≤ one(T) &&
  zero(T) ≤ λ₂ ≤ one(T) &&
  zero(T) ≤ λ₃ ≤ one(T)
end

function Base.in(p::Point{3,T}, t::Triangle{3,T}) where {T}
  # given coordinates
  a, b, c = t.vertices

  # evaluate vectors defining geometry
  v₁ = b - a
  v₂ = c - a
  v₃ = p - a

  # calculate required dot products
  d₁₁ = v₁ ⋅ v₁
  d₁₂ = v₁ ⋅ v₂
  d₂₂ = v₂ ⋅ v₂
  d₃₁ = v₃ ⋅ v₁
  d₃₂ = v₃ ⋅ v₂

  # calculate reused denominator
  d = d₁₁ * d₂₂ - d₁₂ * d₁₂

  # barycentric coordinates
  λ₂ = (d₂₂ * d₃₁ - d₁₂ * d₃₂) / d
  λ₃ = (d₁₁ * d₃₂ - d₁₂ * d₃₁) / d

  # barycentric check
  (λ₂ ≥ zero(T)) && (λ₃ ≥ zero(T)) && ((λ₂ + λ₃) ≤ one(T))
end

"""
    normal(triangle)

Determine the normalised normal of `triangle` in three
dimensions
"""
function normal(t::Triangle{3})
  a, b, c = t.vertices
  n = (b - a) × (c - a)
  n / norm(n)
end

# ------------
# QUADRANGLES
# ------------

# Coons patch https://en.wikipedia.org/wiki/Coons_patch
function (q::Quadrangle)(u, v)
  c₀₀, c₀₁, c₁₁, c₁₀ = coordinates.(q.vertices)
  Point(c₀₀*(1-u)*(1-v) + c₀₁*u*(1-v) + c₁₀*(1-u)*v + c₁₁*u*v)
end


text:
@testset "STRtree" begin
    # front page examples
    p1 = readgeom("POLYGON((0 0,1 0,1 1,0 0))")
    p2 = readgeom("POLYGON((0 0,1 0,1 1,0 1,0 0))")
    p3 = readgeom("POLYGON((2 0,3 0,3 1,2 1,2 0))")

    # query polygon
    p4 = readgeom("POLYGON((0.5 0.5, 1.5 0.5, 1.5 1.5, 0.5 1.5, 0.5 0.5))")

    tree = STRtree((p1, p2, p3))
    res = query(tree, p4)
    @test res isa Vector{Polygon}
    @test res == [p1, p2]

    # heterogeneous tree with linestring and polygon
    ls = readgeom("LINESTRING(0 0,1 0,1 1,0 0)")
    tree = STRtree((p1, p2, p3, ls))
    res = query(tree, p4)
    @test res isa Vector{GeoInterface.AbstractGeometry}
    @test res == [p1, p2, ls]

    # Let's accidentally "drop" p2. Since it is in the tree it won't be finalized.
    p2 = nothing
    p5 = readgeom("POLYGON((0 0,1 0,1 1,0 1,0 0))") # equal to p2
    res = query(tree, p4)
    @test equals(res[1], p1)
    @test equals(res[2], p5)
end


text:
using Documenter
using SkyCoords

DocMeta.setdocmeta!(SkyCoords, :DocTestSetup, :(using SkyCoords); recursive = true)

makedocs(
    sitename = "SkyCoords",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages = ["Home" => "index.md", "API/Reference" => "api.md"],
    modules = [SkyCoords],
    strict = true,
)

deploydocs(repo = "github.com/JuliaAstro/SkyCoords.jl.git")


text:
using Protoscon
using Test

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


text:
# work around for julia issue #28948, Gurobi.jl issue #152
if VERSION ≤ v"0.7-"
    sparse_transpose(A) = sparse(transpose(A))
else
    sparse_transpose(A) = SparseMatrixCSC(transpose(A))
end

## Add Linear constraints

# add single constraint

function add_constr!(model::Model, inds::IVec, coeffs::FVec, rel::Cchar, rhs::Float64)
    length(inds) == length(coeffs) || error("Inconsistent argument dimensions.")
    ret = @grb_ccall(addconstr, Cint, (
        Ptr{Cvoid},    # model
        Cint,         # numnz
        Ptr{Cint},    # cind
        Ptr{Float64}, # cvals
        Cchar,        # sense
        Float64,      # rhs
        Ptr{UInt8}    # name
        ),
        model, length(inds), inds .- Cint(1), coeffs, rel, rhs, C_NULL)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function add_constr!(model::Model, inds::Vector, coeffs::Vector, rel::GChars, rhs::Real)
    add_constr!(model, ivec(inds), fvec(coeffs), cchar(rel), Float64(rhs))
end

function add_constr!(model::Model, coeffs::Vector, rel::GChars, rhs::Real)
    inds = findall(x->x!=0, coeffs)
    vals = coeffs[inds]
    add_constr!(model, inds, vals, rel, rhs)
end

# add multiple constraints

function add_constrs!(model::Model, cbegins::IVec, inds::IVec, coeffs::FVec, rel::CVec, rhs::FVec)
    m = length(cbegins)
    nnz = length(inds)
    (m == length(rel) == length(rhs) && nnz == length(coeffs)) || error("Inconsistent argument dimensions.")

    if m > 0
        ret = @grb_ccall(addconstrs, Cint, (
            Ptr{Cvoid},    # model
            Cint,         # num constraints
            Cint,         # num non-zeros
            Ptr{Cint},    # cbeg
            Ptr{Cint},    # cind
            Ptr{Float64}, # cval
            Ptr{Cchar},   # sense
            Ptr{Float64}, # rhs
            Ptr{UInt8}    # names
            ),
            model, m, nnz, cbegins .- Cint(1), inds .- Cint(1), coeffs,
            rel, rhs, C_NULL)

        if ret != 0
            throw(GurobiError(model.env, ret))
        end
    end
    nothing
end

function add_constrs!(model::Model, cbeg::Vector, inds::Vector, coeffs::Vector, rel::GCharOrVec, rhs::Vector)
    add_constrs!(model, ivec(cbeg), ivec(inds), fvec(coeffs), cvecx(rel, length(cbeg)), fvec(rhs))
end

function add_constrs_t!(model::Model, At::SparseMatrixCSC{Float64}, rel::GCharOrVec, b::Vector)
    n, m = size(At)
    (m == length(b) && n == num_vars(model)) || error("Incompatible argument dimensions.")
    add_constrs!(model, At.colptr[1:At.n], At.rowval, At.nzval, rel, b)
end

function add_constrs_t!(model::Model, At::Matrix{Float64}, rel::GCharOrVec, b::Vector)
    n, m = size(At)
    (m == length(b) && n == num_vars(model)) || error("Incompatible argument dimensions.")
    add_constrs_t!(model, sparse(At), rel, b)
end

function add_constrs!(model::Model, A::CoeffMat, rel::GCharOrVec, b::Vector{Float64})
    m, n = size(A)
    (m == length(b) && n == num_vars(model)) || error("Incompatible argument dimensions.")
    add_constrs_t!(model, sparse_transpose(A), rel, b)
end

# add single range constraint

function add_rangeconstr!(model::Model, inds::IVec, coeffs::FVec, lb::Float64, ub::Float64)
    ret = @grb_ccall(addrangeconstr, Cint, (
        Ptr{Cvoid},    # model
        Cint,         # numnz
        Ptr{Cint},    # cind
        Ptr{Float64}, # cvals
        Float64,      # lower
        Float64,      # upper
        Ptr{UInt8}    # name
        ),
        model, length(inds), inds .- Cint(1), coeffs, lb, ub, C_NULL)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function add_rangeconstr!(model::Model, inds::Vector, coeffs::Vector, lb::Real, ub::Real)
    add_rangeconstr!(model, ivec(inds), fvec(coeffs), Float64(lb), Float64(ub))
end

function add_rangeconstr!(model::Model, coeffs::Vector, lb::Real, ub::Real)
    inds = find(coeffs)
    vals = coeffs[inds]
    add_rangeconstr!(model, inds, vals, lb, ub)
end

# add multiple range constraints

function add_rangeconstrs!(model::Model, cbegins::IVec, inds::IVec, coeffs::FVec, lb::FVec, ub::FVec)
    m = length(cbegins)
    nnz = length(inds)
    (m == length(lb) == length(ub) && nnz == length(coeffs)) || error("Incompatible argument dimensions.")

    if m > 0
        ret = @grb_ccall(addrangeconstrs, Cint, (
            Ptr{Cvoid},    # model
            Cint,         # num constraints
            Cint,         # num non-zeros
            Ptr{Cint},    # cbeg
            Ptr{Cint},    # cind
            Ptr{Float64}, # cval
            Ptr{Float64}, # lower
            Ptr{Float64}, # upper
            Ptr{UInt8}    # names
            ),
            model, m, nnz, cbegins .- Cint(1), inds .- Cint(1), coeffs, lb, ub, C_NULL)

        if ret != 0
            throw(GurobiError(model.env, ret))
        end
    end
    nothing
end

function add_rangeconstrs!(model::Model, cbeg::Vector, inds::Vector, coeffs::Vector, lb::Vector, ub::Vector)
    add_rangeconstrs!(model, ivec(cbeg), ivec(inds), fvec(coeffs), fvec(lb), fvec(ub))
end

function add_rangeconstrs_t!(model::Model, At::SparseMatrixCSC{Float64}, lb::Vector, ub::Vector)
    add_rangeconstrs!(model, At.colptr[1:At.n], At.rowval, At.nzval, lb, ub)
end

function add_rangeconstrs_t!(model::Model, At::Matrix{Float64}, lb::Vector, ub::Vector)
    add_rangeconstrs_t!(model, sparse(At), lb, ub)
end

function add_rangeconstrs!(model::Model, A::CoeffMat, lb::Vector, ub::Vector)
    m, n = size(A)
    (m == length(lb) == length(ub) && n == num_vars(model)) || error("Incompatible argument dimensions.")

    add_rangeconstrs_t!(model, sparse_transpose(A), lb, ub)
end

function get_constrmatrix(model::Model)
    get_constrs(model::Model, 1, num_constrs(model))
end
function get_constrs(model::Model, start::Integer, len::Integer)

    m = num_constrs(model)
    @assert start >= 1
    @assert len >= 0
    @assert start + len <= m + 1
    n = num_vars(model)
    numnzP = Ref{Cint}()
    cbeg = Array{Cint}(undef, len+1)

    ret = @grb_ccall(getconstrs, Cint, (
        Ptr{Cvoid},
        Ptr{Cint},
        Ptr{Cint},
        Ptr{Cint},
        Ptr{Cdouble},
        Cint, # start
        Cint #len
        ),
        model, numnzP, cbeg, C_NULL, C_NULL, Cint(start-1), Cint(len))

    nnz = numnzP[]
    cind = Array{Cint}(undef, nnz)
    cval = Array{Cdouble}(undef, nnz)
    ret = @grb_ccall(getconstrs, Cint, (
                     Ptr{Cvoid},
                     Ptr{Cint},
                     Ptr{Cint},
                     Ptr{Cint},
                     Ptr{Cdouble},
                     Cint, # start
                     Cint #len
                     ),
                     model, numnzP, cbeg, cind, cval, Cint(start-1), Cint(len))
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    cbeg[end] = nnz
    I = Array{Int64}(undef, nnz)
    J = Array{Int64}(undef, nnz)
    V = Array{Float64}(undef, nnz)
    for i in 1:length(cbeg)-1
        for j in (cbeg[i]+1):cbeg[i+1]
            I[j] = i
            J[j] = cind[j]+1
            V[j] = cval[j]
        end
    end
    return sparse(I, J, V, m, n)
end

const GRB_SOS_TYPE1 = convert(Cint, 1)
const GRB_SOS_TYPE2 = convert(Cint, 2)

function add_sos!(model::Model, sostype::Symbol, idx::Vector{Int}, weight::Vector{Cdouble})
    ((nelem = length(idx)) == length(weight)) || error("Index and weight vectors of unequal length")
    (sostype == :SOS1) ? (typ = GRB_SOS_TYPE1) : ( (sostype == :SOS2) ? (typ = GRB_SOS_TYPE2) : error("Invalid SOS constraint type") )
    ret = @grb_ccall(addsos, Cint, (
                     Ptr{Cvoid},
                     Cint,
                     Cint,
                     Ptr{Cint},
                     Ptr{Cint},
                     Ptr{Cint},
                     Ptr{Cdouble}
                     ),
                     model, convert(Cint, 1), convert(Cint, nelem), Cint[typ], Cint[0], convert(Vector{Cint}, idx.-1), weight)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
end
function del_sos!(model::Model, idx::Vector{Cint})
    numdel = length(idx)
    ret = @grb_ccall(delsos, Cint, (
                     Ptr{Cvoid},
                     Cint,
                     Ptr{Cint}),
                     model, convert(Cint,numdel), ivec(idx.-1) )
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
end

get_sos_matrix(model::Model) = get_sos(model::Model, 1, num_sos(model))
function get_sos(model::Model, start::Integer, len::Integer)
    m = num_sos(model)
    sostype = Array{Cint}(undef, m)
    @assert m > 0
    @assert start <= m
    @assert len <= m
    n = num_vars(model)
    numnzP = Ref{Cint}()
    cbeg = Array{Cint}(undef, len+1)

    ret = @grb_ccall(getsos, Cint, (
        Ptr{Cvoid},
        Ptr{Cint},
        Ptr{Cint},
        Ptr{Cint},
        Ptr{Cint},
        Ptr{Cdouble},
        Cint, # start
        Cint #len
        ),
        model, numnzP, sostype, cbeg, C_NULL, C_NULL, Cint(start-1), Cint(len))

    nnz = numnzP[]

    cind = Array{Cint}(undef, nnz)
    cval = Array{Cdouble}(undef, nnz)
    ret = @grb_ccall(getsos, Cint, (
                     Ptr{Cvoid},
                     Ptr{Cint},
                     Ptr{Cint},
                     Ptr{Cint},
                     Ptr{Cint},
                     Ptr{Cdouble},
                     Cint, # start
                     Cint #len
                     ),
                     model, numnzP, sostype, cbeg, cind, cval, Cint(start-1), Cint(len))
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    cbeg[end] = nnz
    I = Array{Int64}(undef, nnz)
    J = Array{Int64}(undef, nnz)
    V = Array{Float64}(undef, nnz)
    for i in 1:length(cbeg)-1
        for j in (cbeg[i]+1):cbeg[i+1]
            I[j] = i
            J[j] = cind[j]+1
            V[j] = cval[j]
        end
    end
    return sparse(I, J, V, m, n), sostype#map(x-> x==Cint(1) ? :SOS1 : :SOS2 ,sostype)
end

del_constrs!(model::Model, idx::T) where {T<:Real} = del_constrs!(model, Cint[idx])
del_constrs!(model::Model, idx::Vector{T}) where {T<:Real} = del_constrs!(model, convert(Vector{Cint},idx))
function del_constrs!(model::Model, idx::Vector{Cint})
    numdel = length(idx)
    ret = @grb_ccall(delconstrs, Cint, (
                     Ptr{Cvoid},
                     Cint,
                     Ptr{Cint}),
                     model, convert(Cint,numdel), idx .- Cint(1))
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
end

function chg_coeffs!(model::Model, cidx::Real, vidx::Real, val::Real)
    ret = @grb_ccall(chgcoeffs, Cint, (
                     Ptr{Cvoid},
                     Cint,
                     Ref{Cint},
                     Ref{Cint},
                     Ref{Float64}),
                     model, 1, cidx - 1, vidx - 1, val)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
end

chg_coeffs!(model::Model, cidx::AbstractVector{<:Real}, vidx::AbstractVector{<:Real}, val::AbstractVector{<:Real}) =
    chg_coeffs!(model, ivec(cidx), ivec(vidx), fvec(val))

function chg_coeffs!(model::Model, cidx::IVec, vidx::IVec, val::FVec)
    (length(cidx) == length(vidx) == length(val)) || error("Inconsistent argument dimensions.")

    numchgs = length(cidx)
    ret = @grb_ccall(chgcoeffs, Cint, (
                     Ptr{Cvoid},
                     Cint,
                     Ptr{Cint},
                     Ptr{Cint},
                     Ptr{Float64}),
                     model, convert(Cint,numchgs), cidx .- Cint(1), vidx .- Cint(1), val)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
end

function getcoeff!(val::FVec, model::Model, cidx::Integer, vidx::Integer)
    Base.depwarn("getcoeff!(val, model, cidx, vidx) is deprecated. Instead you can retrieve a coefficient without allocating a vector by doing `val = getcoeff(model, cidx, vidx)`", :grb_getcoeff)
    @assert length(val) == 1
    ret = @grb_ccall(getcoeff, Cint, (
        Ptr{Cvoid},
        Cint,
        Cint,
        Ptr{Float64}),
        model, cidx-Cint(1), vidx-Cint(1), val)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
end

function getcoeff(model::Model, cidx::Integer, vidx::Integer)
    out = Ref{Float64}()
    ret = @grb_ccall(getcoeff, Cint, (
        Ptr{Cvoid},
        Cint,
        Cint,
        Ref{Float64}),
        model, cidx - 1, vidx - 1, out)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    out[]
end


text:
module FluxBooks

# Write your package code here.

end


text:
function create_data_for_test()
    data = zeros(Float32,2,1000)
    data[:,1:250] .= [-1,-1]
    data[:,251:500] .= [-1,1]
    data[:,501:750] .= [1,-1]
    data[:,751:1000] .= [1,1]
    return data
end

# @testset "Testing Module (Determinstic)" begin
#     data = create_data_for_test()
#     labels,clusters,weights = fit(data,100.0, iters = 200, seed = 12345, burnout = 5)
#     @test all(data[:,1:250] .== data[:,1])
#     @test all(data[:,251:500] .== data[:,251])
#     @test all(data[:,501:750] .== data[:,501])
#     @test all(data[:,751:1000] .== data[:,751])
#     @test length(clusters) == 4
#     println(weights)
#     @test all(weights .>= 0.15)
#     labels_histogram = get_labels_histogram(labels)
#     for (k,v) in labels_histogram
#         @test v == 250
#     end
# end

# @testset "Testing Module (Random mess)" begin
#     @everywhere Random.seed!(12345)
#     x,labels,clusters = generate_gaussian_data(10^6,3,10,100.0)

#     hyper_params = DPMMSubClustersStreaming.niw_hyperparams(Float32(1.0),
#                zeros(Float32,3),
#                Float32(5),
#                Matrix{Float32}(I, 3, 3)*1)

#     dp = dp_parallel(x,hyper_params,Float32(1000000000000000000000.0), 100,1,nothing,true,false,15,labels)
#     println("Full: ",mutualinfo(Int.(labels),dp[1].group.labels,normed=true))
#     @test length(dp[1].group.local_clusters) > 1
# end

# @testset "Streaming Data NIW" begin
#     @everywhere Random.seed!(12345)
#     x,labels,clusters = generate_gaussian_data(10^6,3,10,100.0)
#     x1 = x[:,1:5:end]
#     x2 = x[:,2:5:end]
#     x3 = x[:,3:5:end]
#     x4 = x[:,4:5:end]
#     x5 = x[:,5:5:end]
#     labels1 = labels[1:5:end]
#     labels2 = labels[2:5:end]
#     labels3 = labels[3:5:end]
#     labels4 = labels[4:5:end]
#     labels5 = labels[5:5:end]

#     hyper_params = DPMMSubClustersStreaming.niw_hyperparams(Float32(1.0),
#                zeros(Float32,3),
#                Float32(5),
#                Matrix{Float32}(I, 3, 3)*1)
#     dp = dp_parallel_streaming(x1,hyper_params,Float32(1000000000000000000000.0), 20,1,nothing,true,false,15,labels1)
#     labels = dp.group.labels
#     println("Part 1: ",mutualinfo(Int.(labels1),labels,normed=true))
#     run_model_streaming(dp,20,2,x2)
#     labels = dp.group.labels
#     println("Part 2: ",mutualinfo(Int.(labels2),labels,normed=true))
#     run_model_streaming(dp,20,4,x3)
#     labels = dp.group.labels
#     println("Part 3: ",mutualinfo(Int.(labels3),labels,normed=true))
#     run_model_streaming(dp,20,6,x4)
#     labels = dp.group.labels
#     println("Part 4: ",mutualinfo(Int.(labels4),labels,normed=true))
#     run_model_streaming(dp,20,8,x5)
#     labels = dp.group.labels
#     println("Part 5: ",mutualinfo(Int.(labels5),labels,normed=true))
#     @test length(dp.group.local_clusters) > 1
# end

# @testset "Streaming Data NIW 2" begin
#     @everywhere Random.seed!(12345)
#     x,labels,clusters = generate_gaussian_data(10^6,3,10,100.0)
#     parts = 10000
#     xs = [x[:,i:parts:end] for i=1:parts]
#     labelss = [labels[i:parts:end] for i=1:parts]
#     hyper_params = DPMMSubClustersStreaming.niw_hyperparams(Float32(1.0),
#                zeros(Float32,3),
#                Float32(5),
#                Matrix{Float32}(I, 3, 3)*1)
#     dp = dp_parallel_streaming(xs[1],hyper_params,Float32(1000000000000000000000.0), 1,1,nothing,true,false,15,labelss[1],0.0001)
#     labels = get_labels(dp)
#     avg_nmi = mutualinfo(Int.(labelss[1]),labels,normed=true)
#     for i=2:parts
#         run_model_streaming(dp,1,i*0.5,xs[i])
#         labels = get_labels(dp)
#         avg_nmi += mutualinfo(Int.(labelss[i]),labels,normed=true)
#     end
#     println("NMI: ",avg_nmi/parts)
#     @test length(dp.group.local_clusters) > 1
# end

@testset "Streaming Data Multinomial 2" begin
    @everywhere Random.seed!(12345)
    x,labels,clusters = generate_mnmm_data(10^4,100,20,50)
    parts = 10
    xs = [x[:,i:parts:end] for i=1:parts]
    labelss = [labels[i:parts:end] for i=1:parts]
    hyper_params = DPMMSubClustersStreaming.multinomial_hyper(ones(Float32,100))
    dp = dp_parallel_streaming(xs[1],hyper_params,Float32(10000.0), 20,1,nothing,true,false,15,labelss[1],0.0001)
    labels = get_labels(dp)
    avg_nmi = mutualinfo(Int.(labelss[1]),labels,normed=true)
    for i=2:parts
        run_model_streaming(dp,20,i*2,xs[i])
        labels = get_labels(dp)
        labels = predict(dp,xs[i])
        avg_nmi += mutualinfo(Int.(labelss[i]),labels,normed=true)
    end
    println("NMI: ",avg_nmi/parts)
    @test length(dp.group.local_clusters) > 1
end

# @testset "Streaming Data NIW 3" begin
#         @everywhere Random.seed!(12345)
#         x,labels,clusters = generate_gaussian_data(10^6,2,20,50.0)
#         labels = Int32.(labels);
#         hyper_params = DPMMSubClustersStreaming.niw_hyperparams(Float32(1.0),
#            zeros(Float32,2),
#            Float32(2),
#            Matrix{Float32}(I, 2, 2)*1)
#         parts = 1000
#         xs = [x[:,i:parts:end] for i=1:parts]
#         labelss = [labels[i:parts:end] for i=1:parts]
#         # First batch, initializing the model
#         dp = dp_parallel_streaming(xs[1],hyper_params,Float32(100.0), 30,1,nothing,true,false,10,labels,0.0001);

#         all_preds = []
#         preds = get_labels(dp)
#         preds2 = predict(dp,xs[1])
#         println(preds == preds2)
#         println(preds[1:100])
#         println(preds2[1:100])
#         push!(all_preds,preds)
#         avg_nmi = mutualinfo(Int.(labelss[1]),preds,normed=true)
#         for i=2:parts
#             run_model_streaming(dp,1,i,xs[i])
#             preds = get_labels(dp)
#             push!(all_preds,preds)
#             avg_nmi += mutualinfo(Int.(labelss[i]),preds,normed=true)
#         end
#         all_labels = reduce(hcat,labelss)
#         all_preds = reduce(hcat,all_preds)
#         println("Avg NMI: ",avg_nmi/parts)
#         println("Full NMI: ",mutualinfo(Int.(all_labels),all_preds,normed=true))
#     end

# @testset "Multinomial Module And save load" begin
#     @everywhere Random.seed!(12345)
#     x,labels,clusters = generate_mnmm_data(10^3,100,20,50)
#     @test size(x,1) == 100
#     @test size(x,2) == 10^3
#     npzwrite("save_load_test/mnm_data.npy",x')

#     dp = dp_parallel("save_load_test/multinomial_params.jl")
#     @test length(dp[1].group.local_clusters) > 1
#     dp = run_model_from_checkpoint("save_load_test/checkpoint_20.jld2")
#     @test length(dp[1].group.local_clusters) > 1
# end


text:
"""
    splitobs(n::Int, [at = 0.7]) -> Tuple

Pre-compute the indices for two disjoint subsets and return them
as a tuple of two ranges. The first range will span the first
`at` fraction of possible indices, while the second range will
cover the rest. These indices are applicable to any data
container of size `n`.

```julia
julia> splitobs(100, at = 0.7)
# (1:70,71:100)
```
"""
splitobs(n::Int; at = 0.7) = splitobs(n, at)

# partition into 2 sets
function splitobs(n::Int, at::AbstractFloat)
    0 < at < 1 || throw(ArgumentError("the parameter \"at\" must be in interval (0, 1)"))
    n1 = clamp(round(Int, at*n), 1, n)
    (1:n1, n1+1:n)
end

# has to be outside the generated function
_ispos(x) = x > 0
# partition into length(at)+1 sets
# we use @generated because we compute "N+1"
@generated function splitobs(n::Int, at::NTuple{N,AbstractFloat}) where N
    quote
        (all(map(_ispos, at)) && sum(at) < 1) || throw(ArgumentError("all elements in \"at\" must be positive and their sum must be smaller than 1"))
        nleft = n
        lst = UnitRange{Int}[]
        for (i, sz) in enumerate(at)
            ni = clamp(round(Int, sz*n), 0, nleft)
            push!(lst, n-nleft+1:n-nleft+ni)
            nleft -= ni
        end
        push!(lst, n-nleft+1:n)
        $(Expr(:tuple, (:(lst[$i]) for i in 1:N+1)...))
    end
end

"""
    splitobs(data, [at = 0.7], [obsdim]) -> Tuple

Split the `data` into multiple subsets proportional to the
value(s) of `at`.

Note that this function will perform the splits statically and
thus not perform any randomization. The function creates a
`NTuple` of data subsets in which the first N-1 elements/subsets
contain the fraction of observations of `data` that is specified
by `at`.

For example, if `at` is a `Float64` then the return-value will be
a tuple with two elements (i.e. subsets), in which the first
element contains the fracion of observations specified by `at`
and the second element contains the rest. In the following code
the first subset `train` will contain the first 70% of the
observations and the second subset `test` the rest.

```julia
train, test = splitobs(X, at = 0.7)
```

If `at` is a tuple of `Float64` then additional subsets will be
created. In this example `train` will have the first 50% of the
observations, `val` will have next 30%, and `test` the last 20%

```julia
train, val, test = splitobs(X, at = (0.5, 0.3))
```

It is also possible to call `splitobs` with multiple data
arguments as tuple, which all must have the same number of total
observations. This is useful for labeled data.

```julia
train, test = splitobs((X, y), at = 0.7)
(x_train,y_train), (x_test,y_test) = splitobs((X, y), at = 0.7)
```

If the observations should be randomly assigned to a subset,
then you can combine the function with `shuffleobs`

```julia
# This time observations are randomly assigned.
train, test = splitobs(shuffleobs((X,y)), at = 0.7)
```

When working with arrays one may want to choose which dimension
represents the observations. By default the last dimension is
assumed, but this can be overwritten.

```julia
# Here we say each row represents an observation
train, test = splitobs(X, obsdim = 1)
```

The functions also provide a type-stable API

```julia
# By avoiding keyword arguments, the compiler can infer the return type
train, test = splitobs((X,y), 0.7)
train, test = splitobs((X,y), 0.7, ObsDim.First())
```

see [`DataSubset`](@ref) for more information on data subsets.

see [`stratifiedobs`](@ref) for a related function that preserves
the target distribution.
"""
splitobs(data; at = 0.7, obsdim = default_obsdim(data)) =
    splitobs(data, at, convert(LearnBase.ObsDimension,obsdim))

# partition into 2 sets
function splitobs(data, at::AbstractFloat, obsdim=default_obsdim(data))
    allowcontainer(splitobs, data) || throw(MethodError(splitobs, (data,at,obsdim)))
    n = nobs(data, obsdim)
    idx1, idx2 = splitobs(n, at)
    datasubset(data, idx1, obsdim), datasubset(data, idx2, obsdim)
end

# partition into length(at)+1 sets
function splitobs(data, at::NTuple{N,AbstractFloat}, obsdim=default_obsdim(data)) where N
    allowcontainer(splitobs, data) || throw(MethodError(splitobs, (data,at,obsdim)))
    n = nobs(data, obsdim)
    map(idx->datasubset(data, idx, obsdim), splitobs(n, at))
end


text:
export TD3Policy, TD3Critic

struct TD3Critic
    critic_1::Flux.Chain
    critic_2::Flux.Chain
end
Flux.@functor TD3Critic
(c::TD3Critic)(s, a) = (inp = vcat(s, a); (c.critic_1(inp), c.critic_2(inp)))

mutable struct TD3Policy{
    BA<:NeuralNetworkApproximator,
    BC<:NeuralNetworkApproximator,
    TA<:NeuralNetworkApproximator,
    TC<:NeuralNetworkApproximator,
    P,
    R<:AbstractRNG,
} <: AbstractPolicy

    behavior_actor::BA
    behavior_critic::BC
    target_actor::TA
    target_critic::TC
    γ::Float32
    ρ::Float32
    batch_size::Int
    start_steps::Int
    start_policy::P
    update_after::Int
    update_freq::Int
    policy_freq::Int
    target_act_limit::Float64
    target_act_noise::Float64
    act_limit::Float64
    act_noise::Float64
    update_step::Int
    rng::R
    replay_counter::Int
    # for logging
    actor_loss::Float32
    critic_loss::Float32
end

"""
    TD3Policy(;kwargs...)

# Keyword arguments

- `behavior_actor`,
- `behavior_critic`,
- `target_actor`,
- `target_critic`,
- `start_policy`,
- `γ = 0.99f0`,
- `ρ = 0.995f0`,
- `batch_size = 32`,
- `start_steps = 10000`,
- `update_after = 1000`,
- `update_freq = 50`,
- `policy_freq = 2` # frequency in which the actor performs a gradient update_step and critic target is updated
- `target_act_limit = 1.0`, # noise added to actor target
- `target_act_noise = 0.1`, # noise added to actor target
- `act_limit = 1.0`, # noise added when outputing action
- `act_noise = 0.1`, # noise added when outputing action
- `update_step = 0`,
- `rng = Random.GLOBAL_RNG`,
"""
function TD3Policy(;
    behavior_actor,
    behavior_critic,
    target_actor,
    target_critic,
    start_policy,
    γ = 0.99f0,
    ρ = 0.995f0,
    batch_size = 64,
    start_steps = 10000,
    update_after = 1000,
    update_freq = 50,
    policy_freq = 2,
    target_act_limit = 1.0,
    target_act_noise = 0.1,
    act_limit = 1.0,
    act_noise = 0.1,
    update_step = 0,
    rng = Random.GLOBAL_RNG,
)
    copyto!(behavior_actor, target_actor)  # force sync
    copyto!(behavior_critic, target_critic)  # force sync
    TD3Policy(
        behavior_actor,
        behavior_critic,
        target_actor,
        target_critic,
        γ,
        ρ,
        batch_size,
        start_steps,
        start_policy,
        update_after,
        update_freq,
        policy_freq,
        target_act_limit,
        target_act_noise,
        act_limit,
        act_noise,
        update_step,
        rng,
        1, # keep track of numbers of replay
        0.0f0,
        0.0f0,
    )
end

# TODO: handle Training/Testing mode
function (p::TD3Policy)(env)
    p.update_step += 1

    if p.update_step <= p.start_steps
        p.start_policy(env)
    else
        D = device(p.behavior_actor)
        s = state(env)
        s = Flux.unsqueeze(s, ndims(s) + 1)
        action = p.behavior_actor(send_to_device(D, s)) |> vec |> send_to_host
        clamp(action[] + randn(p.rng) * p.act_noise, -p.act_limit, p.act_limit)
    end
end

function RLBase.update!(
    p::TD3Policy,
    traj::CircularArraySARTTrajectory,
    ::AbstractEnv,
    ::PreActStage,
)
    length(traj) > p.update_after || return
    p.update_step % p.update_freq == 0 || return
    inds, batch = sample(p.rng, traj, BatchSampler{SARTS}(p.batch_size))
    update!(p, batch)
end

function RLBase.update!(p::TD3Policy, batch::NamedTuple{SARTS})
    s, a, r, t, s′ = send_to_device(device(p.behavior_actor), batch)

    actor = p.behavior_actor
    critic = p.behavior_critic

    # !!! we have several assumptions here, need revisit when we have more complex environments
    # state is vector
    # action is scalar
    target_noise =
        clamp.(
            randn(p.rng, Float32, 1, p.batch_size) .* p.target_act_noise,
            -p.target_act_limit,
            p.target_act_limit,
        )
    # add noise and clip to tanh bounds
    a′ = clamp.(p.target_actor(s′) + target_noise, -1.0f0, 1.0f0)

    q_1′, q_2′ = p.target_critic(s′, a′)
    y = r .+ p.γ .* (1 .- t) .* (min.(q_1′, q_2′) |> vec)
    a = Flux.unsqueeze(a, 1)

    gs1 = gradient(Flux.params(critic)) do
        q1, q2 = critic(s, a)
        loss = mse(q1 |> vec, y) + mse(q2 |> vec, y)
        ignore() do
            p.critic_loss = loss
        end
        loss
    end
    update!(critic, gs1)

    if p.replay_counter % p.policy_freq == 0
        gs2 = gradient(Flux.params(actor)) do
            actions = actor(s)
            loss = -mean(critic.model.critic_1(vcat(s, actions)))
            ignore() do
                p.actor_loss = loss
            end
            loss
        end
        update!(actor, gs2)
        # polyak averaging
        for (dest, src) in zip(
            Flux.params([p.target_actor, p.target_critic]),
            Flux.params([actor, critic]),
        )
            dest .= p.ρ .* dest .+ (1 - p.ρ) .* src
        end
        p.replay_counter = 1
    end
    p.replay_counter += 1
end