text:
struct Dist
  item_weights::Vector{Matrix{Float}}
  aux_item_weights1::Vector{Matrix{Float}}
  aux_item_weights2::Vector{Matrix{Float}}

  comp_weights::Vector{Float}
  aux_comp_weights1::Vector{Float}
  aux_comp_weights2::Vector{Float}

  aux_weight::Vector{Float}

  sorted_targets::Vector{Float}
  targets::Vector{Float}
end

const ZERO = Float(0)
const ONE = Float(1)
const TWO = Float(2)

function Dist()
  item_weights = Vector{Matrix{Float}}()
  aux_item_weights1 = Vector{Matrix{Float}}()
  aux_item_weights2 = Vector{Matrix{Float}}()

  comp_weights = Vector{Float}()
  aux_comp_weights1 = Vector{Float}()
  aux_comp_weights2 = Vector{Float}()

  sorted_targets = Vector{Float}()
  targets = Vector{Float}()

  Dist(item_weights, aux_item_weights1, aux_item_weights2, comp_weights, aux_comp_weights1, aux_comp_weights2, Float[0], sorted_targets, targets)
end

function init!(dist::Dist, M::Int, N::Int, compnum::Int; max_init::Float = (compnum == 1) ? ZERO : Float(3))

  if compnum < 1
    error("compnum must be a postive integer")
  end

  resize!(dist.comp_weights, compnum)
  fill!(resize!(dist.comp_weights, compnum), Float(0))
  fill!(resize!(dist.aux_comp_weights1, compnum), Float(0))
  fill!(resize!(dist.aux_comp_weights2, compnum), Float(0))

  resize!(dist.item_weights, compnum)
  resize!(dist.aux_item_weights1, compnum)
  resize!(dist.aux_item_weights2, compnum)

  for i = 1:compnum
    dist.item_weights[i] = (rand(Float, M, N) - Float(.5))*max_init
    dist.aux_item_weights1[i] = fill(Float(0), M, N)
    dist.aux_item_weights2[i] = fill(Float(0), M, N)
  end

  dist.aux_weight[1] = Float(0)

  resize!(dist.sorted_targets, 0)
  resize!(dist.targets, 0)

  nothing
end

compnum(dist) = length(dist.comp_weights)

function randcomp(dist::Dist)
  if length(dist.comp_weights) == 0
    error("dist isn't initialized")
  end

  return wsample(1:length(dist.comp_weights), exp.(dist.comp_weights))
end

function randseq(dist::Dist, comp::Int = randcomp(dist))
  if length(dist.item_weights) == 0
    error("dist isn't initialized")
  end

  item_weights = exp.(dist.item_weights[comp])

  M, N = size(item_weights)
  seq = zeros(Int, N)
  for i = 1:N
    seq[i] = wsample(1:M, item_weights[:, i])
  end

  return seq, comp
end

function mapseq(dist::Dist, comp::Int)
  if length(dist.item_weights) == 0
    error("dist isn't initialized")
  end

  item_weights = dist.item_weights[comp]

  M, N = size(item_weights)
  seq = zeros(Int, N)
  for i = 1:N
    seq[i] = indmax(item_weights[:, i])
  end

  return seq
end

function weights(dist::Dist)
  vcat(dist.comp_weights, map(vec, dist.item_weights)...)
end

compprobs(dist::Dist) = probmat(dist.comp_weights)

function itemprobs(dist::Dist)
  return cat(3, map(probmat, dist.item_weights)...)
end

probs(dist::Dist) = vcat(compprobs(dist), vec(itemprobs(dist)))

@enum UpdateType sgd=1 momentum=2 adagrad=3 rmsprop=4 adam=5

function update!(dist::Dist, seq::Vector{Int}, comp::Int, y::Float;
  targets_num::Int = 100,
  le::Bool = false,
  incd_last::Bool = false,
  early_updates::Bool = false,
  update::UpdateType = adagrad,
  eta::Float = Float(1e-2),
  eta_items::Float = eta,
  eta_comps::Float = eta,
  lambda::Float = Float(0),
  lambda_items::Float = lambda,
  lambda_comps::Float = lambda,
  tau1::Float = Float(.9),
  tau2::Float = Float(1e-3),
  delta::Float = Float(1e-6),
  max_val::Float = Float(88),
  min_val::Float = log(eps(Float)))

  # sizes
  M, N = size(dist.item_weights[1])
  C = length(dist.comp_weights)
  n = length(dist.sorted_targets)+1

  if targets_num < 10
    error("targets_num must be higher or equal to 10")
  end

  push!(dist.targets, y)
  update_params = (n > targets_num) || early_updates
  if n > targets_num
    del_y = splice!(dist.targets, 1)
    del_sorted_idx = searchsortedfirst(dist.sorted_targets, del_y)
    splice!(dist.sorted_targets, del_sorted_idx )
    n -= 1
  end

  if le
    idx = searchsortedlast(dist.sorted_targets, y)+1
  else
    idx = searchsortedfirst(dist.sorted_targets, y)
  end
  insert!(dist.sorted_targets, idx, y)
  if n > targets_num
    splice!(dist.sorted_targets, 1)
    idx -= 1
    n -= 1
  end
  # empirical CDF
  py = (n > 1) ? Float((idx - incd_last) / (n - incd_last)) : Float(.5)
  # example weight, distributed as uniform discrete between -1 and 1
  wy = TWO*py - ONE

  # likelihood P(X;w)
  x = falses(M, N)
  for i = 1:N
    x[seq[i], i] = true
  end
  pc = probmat(dist.comp_weights)
  pitems = map(probmat, dist.item_weights)

  # probabilites, log-probalities, and objective
  log_pc = (dist.comp_weights .- logsumexp(dist.comp_weights))
  log_pxgc = map(w -> sum(w[x] - vec(logsumexp(w))), dist.item_weights)
  log_qc = log_pc + log_pxgc
  log_px = logsumexp(log_qc)[1]
  # posterior of C given X
  qc = exp.(log_qc - log_px)
  # correct numerical inaccuracies
  qc ./= sum(qc)

  # for learning the distribution we minimize the following loss
  loss = -wy*log_px

  if update_params
    # update dist weight
    if (update == adagrad) || (update == adam)
      dist.aux_weight[1] += ONE
    elseif (update ==rmsprop)
      dist.aux_weight[1] = (1 - tau2)*dist.aux_weight[1] + tau2
    end

    # update w[c]
    comp_grad = -wy.*(qc - pc)
    update_field!(dist, comp_grad, :comp_weights, :aux_comp_weights1, :aux_comp_weights2,
                  Val{update}(), eta_comps, lambda_comps, tau1, tau2, delta, max_val, min_val)

    # update w[i,j,c]
    item_grad = [-wy.*qc[i].*(x - pitems[i]) for i in 1:C]
    update_field!(dist, item_grad, :item_weights, :aux_item_weights1, :aux_item_weights2,
                  Val{update}(), eta_items, lambda_items, tau1, tau2, delta, max_val, min_val)
  end

  return loss
end

function update_field!(dist::Dist,
  grad::Array,
  field_name::Symbol,
  aux_field_name1::Symbol,
  aux_field_name2::Symbol,
  ::Val{sgd},
  eta::Float,
  lambda::Float,
  tau1::Float,
  tau2::Float,
  delta::Float,
  max_val::Float,
  min_val::Float)

  for i in 1:length(getfield(dist, field_name))
    getfield(dist, field_name)[i] = min.(max.(
                          (ONE - eta*lambda).*getfield(dist, field_name)[i] - eta.*grad[i]
                          , min_val), max_val)
  end

  return nothing
end

function update_field!(dist::Dist,
  grad::Array,
  field_name::Symbol,
  aux_field_name1::Symbol,
  aux_field_name2::Symbol,
  ::Val{momentum},
  eta::Float,
  lambda::Float,
  tau1::Float,
  tau2::Float,
  delta::Float,
  max_val::Float,
  min_val::Float)

  for i in 1:length(getfield(dist, field_name))
    getfield(dist, aux_field_name1)[i] = tau1.*getfield(dist, aux_field_name1)[i] + eta.*grad[i]
    getfield(dist, field_name)[i] = min.(max.(
                          (ONE - eta*lambda).*getfield(dist, field_name)[i] - getfield(dist, aux_field_name1)[i]
                          , min_val), max_val)
  end

  return nothing
end

function update_field!(dist::Dist,
  grad::Array,
  field_name::Symbol,
  aux_field_name1::Symbol,
  aux_field_name2::Symbol,
  ::Val{adagrad},
  eta::Float,
  lambda::Float,
  tau1::Float,
  tau2::Float,
  delta::Float,
  max_val::Float,
  min_val::Float)

  for i in 1:length(getfield(dist, field_name))
    getfield(dist, aux_field_name1)[i] .+= grad[i].^2
    aux = sqrt.(getfield(dist, aux_field_name1)[i]./dist.aux_weight[1]) .+ delta

    getfield(dist, field_name)[i] = min.(max.(
                          (aux.*getfield(dist, field_name)[i] - eta.*grad[i])./(aux .+ eta*lambda)
                          , min_val), max_val)
  end

  return nothing
end

function update_field!(dist::Dist,
  grad::Array,
  field_name::Symbol,
  aux_field_name1::Symbol,
  aux_field_name2::Symbol,
  ::Val{rmsprop},
  eta::Float,
  lambda::Float,
  tau1::Float,
  tau2::Float,
  delta::Float,
  max_val::Float,
  min_val::Float)

  for i in 1:length(getfield(dist, field_name))
    getfield(dist, aux_field_name1)[i] = (1 - tau2).*getfield(dist, aux_field_name1)[i] + tau2.*grad[i].^2
    aux = sqrt.(getfield(dist, aux_field_name1)[i]./dist.aux_weight[1]) .+ delta

    getfield(dist, field_name)[i] = min.(max.(
                          (aux.*getfield(dist, field_name)[i] - eta.*grad[i])./(aux .+ eta*lambda)
                          , min_val), max_val)
  end

  return nothing
end

function update_field!(dist::Dist,
  grad::Array,
  field_name::Symbol,
  aux_field_name1::Symbol,
  aux_field_name2::Symbol,
  ::Val{adam},
  eta::Float,
  lambda::Float,
  tau1::Float,
  tau2::Float,
  delta::Float,
  max_val::Float,
  min_val::Float)

  for i in 1:length(getfield(dist, field_name))
    est_fix1 = ONE/(ONE - (ONE - tau1)^dist.aux_weight[1])
    est_fix2 = ONE/(ONE - (ONE - tau2)^dist.aux_weight[1])

    getfield(dist, aux_field_name1)[i] = (ONE - tau1).*getfield(dist, aux_field_name1)[i] + tau1.*grad[i]
    getfield(dist, aux_field_name2)[i] = (ONE - tau2).*getfield(dist, aux_field_name2)[i] + tau2.*grad[i].^2

    aux = sqrt.(max.(est_fix2*getfield(dist, aux_field_name2)[i],0)) .+ delta
    getfield(dist, field_name)[i] = min.(max.(
                                        (aux.*getfield(dist, field_name)[i] - (eta*est_fix1).*getfield(dist, aux_field_name1)[i])./(aux .+ eta*lambda)
                                        , min_val), max_val)
  end

  return nothing
end


text:
# TODO integration underway with Manifolds.jl, see RoME #244, also see IIF #467 regarding consolidation effort.

import DistributedFactorGraphs: getVariableType

export projectCartesian

"""
$(TYPEDEF)

XY Euclidean manifold variable node softtype.
"""
@defVariable Point2 TranslationGroup(2) [0.0;0.0]

"""
$(TYPEDEF)

XYZ Euclidean manifold variable node softtype.

Example
-------
```julia
p3 = Point3()
```
"""
@defVariable Point3 TranslationGroup(3) [0;0;0.0]

"""
$(TYPEDEF)

Pose2 is a SE(2) mechanization of two Euclidean translations and one Circular rotation, used for general 2D SLAM.
"""
@defVariable Pose2 SpecialEuclidean(2) ProductRepr([0;0.0],[1 0; 0 1.0])

"""
$(TYPEDEF)

Pose3 is currently a Euler angle mechanization of three Euclidean translations and three Circular rotation.

Future:
------
- Work in progress on AMP3D for proper non-Euler angle on-manifold operations.
- TODO the AMP upgrade is aimed at resolving 3D to Quat/SE3/SP3 -- current Euler angles will be replaced
"""
@defVariable Pose3 SpecialEuclidean(3) ProductRepr([0;0;0.0],[1 0 0; 0 1 0; 0 0 1.0])

"""
$(TYPEDEF)

Dynamic point in 2D space with velocity components: `x, y, dx/dt, dy/dt`

"""
@defVariable DynPoint2 TranslationGroup(4) zeros(4)

"""
$(TYPEDEF)

Dynamic pose variable with velocity components: `x, y, theta, dx/dt, dy/dt`

Note
- The `SE2E2_Manifold` definition used currently is a hack to simplify the transition to Manifolds.jl, see #244
- Replaced `SE2E2_Manifold` hack with `ProductManifold(SpecialEuclidean(2), TranslationGroup(2))`, confirm if it is correct.
"""
@defVariable DynPose2 Manifolds.ProductGroup(ProductManifold(SpecialEuclidean(2), TranslationGroup(2))) ProductRepr(ProductRepr([0;0.0],[1 0; 0 1.0]),[0;0.0])

"""
$SIGNATURES

Function to project only XY data onto Cartesian plane for 2D plotting.
"""
projectCartesian(pose::Union{<:Point2,<:Point3,<:Pose2,<:Pose3,<:DynPoint2,<:DynPose2},
                  x::Vector{Float64}) = [x[1]; x[2]; 0]

#

# Still experimental
# export BearingRange2
@defVariable BearingRange2 BearingRange_Manifold ProductRepr(0.0,0.0)

#

ext: jl


text:
# ---
# title: 35. Search Insert Position
# id: problem35
# author: zhwang
# date: 2022-02-12
# difficulty: Easy
# categories: Array, Binary Search
# link: <https://leetcode.com/problems/search-insert-position/description/>
# hidden: true
# ---
#
# Given a sorted array of distinct integers and a target value, return the index
# if the target is found. If not, return the index where it would be if it were
# inserted in order.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,3,5,6], target = 5
#     Output: 2
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1,3,5,6], target = 2
#     Output: 1
#
#
# **Example 3:**
#
#
#
#     Input: nums = [1,3,5,6], target = 7
#     Output: 4
#
#
# **Example 4:**
#
#
#
#     Input: nums = [1,3,5,6], target = 0
#     Output: 0
#
#
# **Example 5:**
#
#
#
#     Input: nums = [1], target = 0
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 104`
#   * `-104 <= nums[i] <= 104`
#   * `nums` contains **distinct** values sorted in **ascending** order.
#   * `-104 <= target <= 104`
#
#
## @lc code=start
using LeetCode

function search_insert(nums::Vector{Int}, target::Int)::Int
    left, right = 1, length(nums)
    while left <= right
        mid = left + ((right - left) >> 1)
        if nums[mid] >= target
            right = mid - 1
        else
            left = mid + 1
        end
    end
    left - 1
end
## @lc code=end

## @lc test=start
@testset "35.search-insert-position.jl" begin
    @test search_insert([1,3,5,6], 5) == 2
    @test search_insert([1,3,5,6], 2) == 1
    @test search_insert([1,3,5,6], 7) == 4
    @test search_insert([1,3,5,6], 0) == 0
end
## @lc test=end

## @lc info=start
# link: [solution to 35](https://leetcode-cn.com/problems/search-insert-position/solution/35-sou-suo-cha-ru-wei-zhi-by-rex-rs-8egt/)
## @lc info=end


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

const _default_float_WT = Float64 # Default type for float weights
const _default_int_WT = Int # Default type for int weights
const _default_LDT = Float64 # Default type for log-density values

"""
    struct DensitySample

A weighted sample drawn according to an statistical density,
e.g. a [`BAT.AbstractDensity`](@ref).

Fields:
    * `v`: Multivariate parameter vector
    * `logd`: log of the value of the density at `v`
    * `weight`: Weight of the sample
    * `info`: Additional info on the provenance of the sample. Content depends
       on the sampling algorithm.
    * aux: Custom user-defined information attatched to the sample.

Constructors:

```julia
DensitySample(
    v::AbstractVector{<:Real},
    logd::Real,
    weight::Real,
    info::Any,
    aux::Any
)
```

Use [`DensitySampleVector`](@ref) to store vectors of multiple samples with
an efficient column-based memory layout.
"""
struct DensitySample{
    P,
    T<:Real,
    W<:Real,
    R,
    Q
}
    v::P
    logd::T
    weight::W
    info::R
    aux::Q
end

export DensitySample

# DensitySample behaves as a scalar type under broadcasting:
@inline Base.Broadcast.broadcastable(shape::DensitySample) = Ref(shape)

import Base.==
function ==(A::DensitySample, B::DensitySample)
    A.v == B.v && A.logd == B.logd &&
        A.weight == B.weight && A.info == B.info && A.aux == B.aux
end

function Base.similar(s::DensitySample{P,T,W,R,Q}) where {P<:AbstractVector{<:Real},T,W,R,Q}
    v = fill!(similar(s.v), eltype(s.v)(NaN))
    logd = convert(T, NaN)
    weight = zero(W)
    info = R()
    aux = Q()
    P_new = typeof(v)
    DensitySample{P_new,T,W,R,Q}(v, logd, weight, info, aux)
end

function _apply_shape(shape::AbstractValueShape, s::DensitySample)
    DensitySample(
        stripscalar(shape(s.v)),
        s.logd,
        s.weight,
        s.info,
        s.aux,
    )
end

(shape::AbstractValueShape)(s::DensitySample) = _apply_shape(shape, s)

"""
    DensitySampleVector

Type alias for `StructArrays.StructArray{<:DensitySample,...}`.

Constructors:

```julia
    DensitySampleVector(
        (
            v::AbstractVector{<:AbstractVector{<:Real}},
            logd::AbstractVector{<:Real},
            weight::AbstractVector{<:Real},
            info::AbstractVector{<:Any},
            aux::AbstractVector{<:Any}
        )
    )
```

```julia
    DensitySampleVector(
            v::AbstractVector,
            logval::AbstractVector{<:Real};
            weight::Union{AbstractVector{<:Real}, Symbol} = fill(1, length(eachindex(v))),
            info::AbstractVector = fill(nothing, length(eachindex(v))),
            aux::AbstractVector = fill(nothing, length(eachindex(v)))
        )
```
    For `weight = :RepetitionWeight`, repeated samples will be dropped and the weight of the corresponding sample will be increased.
"""
const DensitySampleVector{
    P,T<:AbstractFloat,W<:Real,R,Q,
    PV<:AbstractVector{P},TV<:AbstractVector{T},WV<:AbstractVector{W},RV<:AbstractVector{R},QV<:AbstractVector{Q}
} = StructArray{
    DensitySample{P,T,W,R,Q},
    1,
    NamedTuple{(:v, :logd, :weight, :info, :aux), Tuple{PV,TV,WV,RV,QV}}
}

export DensitySampleVector

function StructArray{DensitySample}(
    contents::Tuple{
        AbstractVector{P},
        AbstractVector{T},
        AbstractVector{W},
        AbstractVector{R},
        AbstractVector{Q},
    }
) where {P,T<:AbstractFloat,W<:Real,R,Q}
    v, logd, weight, info, aux = contents
    StructArray{DensitySample{P,T,W,R,Q}}(contents)
end

DensitySampleVector(contents::NTuple{5,Any}) = StructArray{DensitySample}(contents)

_create_undef_vector(::Type{T}, len::Integer) where T = Vector{T}(undef, len)

function DensitySampleVector{P,T,W,R,Q}(::UndefInitializer, len::Integer, npar::Integer) where {
    PT<:Real, P<:AbstractVector{PT}, T<:AbstractFloat, W<:Real, R, Q
}
    contents = (
        VectorOfSimilarVectors(ElasticArray{PT}(undef, npar, len)),
        Vector{T}(undef, len),
        Vector{W}(undef, len),
        _create_undef_vector(R, len),
        _create_undef_vector(Q, len)
    )

    DensitySampleVector(contents)
end

DensitySampleVector(::Type{S}, varlen::Integer) where {P<:AbstractVector{<:Real},T<:AbstractFloat,W<:Real,R,Q,S<:DensitySample{P,T,W,R,Q}} =
    DensitySampleVector{P,T,W,R,Q}(undef, 0, varlen)

function DensitySampleVector(
    v::AbstractVector,
    logval::AbstractVector{<:Real};
    weight::Union{AbstractVector{<:Real}, Symbol} = fill(1, length(eachindex(v))),
    info::AbstractVector = fill(nothing, length(eachindex(v))),
    aux::AbstractVector = fill(nothing, length(eachindex(v)))
)
    if weight == :RepetitionWeight
        idxs, weight = repetition_to_weights(v)
        return DensitySampleVector((ArrayOfSimilarArrays(v[idxs]), logval[idxs], weight, info[idxs], aux[idxs]))
    else
        return DensitySampleVector((ArrayOfSimilarArrays(v), logval, weight, info, aux))
    end
end

# Specialize getindex to properly support ArraysOfArrays and similar, preventing
# conversion to exact element type:

@inline Base.getindex(A::StructArray{<:DensitySample}, I::Int...) =
    DensitySample(A.v[I...], A.logd[I...], A.weight[I...], A.info[I...], A.aux[I...])

@inline Base.getindex(A::StructArray{<:DensitySample}, I::Union{Int,AbstractArray,Colon}...) =
    DensitySampleVector((A.v[I...], A.logd[I...], A.weight[I...], A.info[I...], A.aux[I...]))

# Specialize IndexStyle, current default for StructArray seems to be IndexCartesian()
Base.IndexStyle(::StructArray{<:DensitySample, 1}) = IndexLinear()

# Specialize comparison, currently StructArray seems fall back to `(==)(A::AbstractArray, B::AbstractArray)`
import Base.==
function(==)(A::DensitySampleVector, B::DensitySampleVector)
    A.v == B.v &&
    A.logd == B.logd &&
    A.weight == B.weight &&
    A.info == B.info &&
    A.aux == B.aux
end

function Base.merge!(X::DensitySampleVector, Xs::DensitySampleVector...)
    for Y in Xs
        append!(X, Y)
    end
    X
end

Base.merge(X::DensitySampleVector, Xs::DensitySampleVector...) = merge!(deepcopy(X), Xs...)

Base.@propagate_inbounds function _bcasted_apply_to_params(f, A::DensitySampleVector)
    DensitySampleVector((
        f.(A.v),
        A.logd,
        A.weight,
        A.info,
        A.aux
    ))
end

Base.copy(
    instance::Base.Broadcast.Broadcasted{
        <:Base.Broadcast.AbstractArrayStyle{1},
        <:Any,
        <:Union{AbstractValueShape,typeof(unshaped)},
        <:Tuple{DensitySampleVector}
    }
) = _bcasted_apply_to_params(instance.f, instance.args[1])

ValueShapes.varshape(A::DensitySampleVector) = elshape(A.v)

function _get_statw(f::Function, samples::DensitySampleVector, resultshape::AbstractValueShape)
    shape = varshape(samples)
    X = unshaped.(samples.v)
    w = FrequencyWeights(samples.weight)
    r_unshaped = f(X, w)
    resultshape(r_unshaped)
end

Statistics.mean(samples::DensitySampleVector) = _get_statw(mean, samples, varshape(samples))
Statistics.var(samples::DensitySampleVector) = _get_statw(var, samples, replace_const_shapes(ValueShapes.const_zero_shape, varshape(samples)))
Statistics.std(samples::DensitySampleVector) = _get_statw(std, samples, replace_const_shapes(ValueShapes.const_zero_shape, varshape(samples)))

function Statistics.median(samples::DensitySampleVector)
    shape = varshape(samples)
    flat_samples = flatview(unshaped.(samples.v))
    n_params = size(flat_samples)[1]
    median_params = Vector{Float64}()

    for param in Base.OneTo(n_params)
        median_param = quantile(flat_samples[param,:], FrequencyWeights(samples.weight), 0.5)
        push!(median_params, median_param)
    end
    shape(median_params)
end

function _get_stat(f::Function, samples::DensitySampleVector)
    shape = varshape(samples)
    X = unshaped.(samples.v)
    r_unshaped = f(X)
    shape(r_unshaped)
end

Base.minimum(samples::DensitySampleVector) = _get_stat(minimum, samples)
Base.maximum(samples::DensitySampleVector) = _get_stat(maximum, samples)

Statistics.cov(samples::DensitySampleVector{<:AbstractVector{<:Real}}) = cov(samples.v, FrequencyWeights(samples.weight))
Statistics.cor(samples::DensitySampleVector{<:AbstractVector{<:Real}}) = cor(samples.v, FrequencyWeights(samples.weight))

function _get_mode(samples::DensitySampleVector)
    shape = varshape(samples)
    i = findmax(samples.logd)[2]
    v_unshaped = unshaped.(samples.v)[i]
    v = copy(shape(v_unshaped))
    (v, i)
end

StatsBase.mode(samples::DensitySampleVector) = _get_mode(samples)[1]

"""
    drop_low_weight_samples(
        samples::DensitySampleVector,
        fraction::Real = 10^-4
    )

*BAT-internal, not part of stable public API.*

Drop `fraction` of the total probability mass from samples to filter out the
samples with the lowest weight.
"""
function drop_low_weight_samples(samples::DensitySampleVector, fraction::Real = 10^-5; threshold::Real=10^-2)
    W = float(samples.weight)
    if minimum(W) / maximum(W) > threshold
        samples
    else
        W_s = sort(W)
        Q = cumsum(W_s)
        Q ./= maximum(Q)
        @assert last(Q) ≈ 1
        ind = searchsortedlast(Q, fraction)
        if ind !== 0
            thresh = W_s[searchsortedlast(Q, fraction)]
            idxs = findall(x -> x >= thresh, samples.weight)
            samples[idxs]
        else
            samples
        end
    end
end

"""
    repetition_to_weights(v::AbstractVector)

*BAT-internal, not part of stable public API.*

Drop (subsequently) repeated samples by adding weights.
"""
function repetition_to_weights(v::AbstractVector)
    idxs = Vector{Int}()
    counts = Vector{Int}()
    push!(idxs, 1)
    push!(counts, 1)
    for i in 2:length(v)
        if v[i] == v[i-1]
            counts[end] += 1
        else
            push!(idxs, i)
            push!(counts, 1)
        end
    end
    return (idxs, counts)
end


text:
using MySQL
using BioMedQuery.DBUtils

dbname="test"

con = init_mysql_database("127.0.0.1", dbname, "root", "", true)
init_mysql_database(con, dbname, true)
init_mysql_database(con, dbname, false)

all_tables = select_all_tables(con)

pritntln(all_tables)


text:
###################################################
## network/solvers/socp.jl
## SOCP that finds new traveltimes to optimize cost function
###################################################

"""
    socpTimes :
    optimize travel times to minimize L1 error from data with given paths

    `continuityConstraint` can be
    - "none" : no constraint
    - "simple": simple ϵ-based continuity
    - "neighborhoods": continuity constraint for each neighborhood
    - "regularized": soft continuity constraint
"""
function socpTimes(s::IterativeState;
                   inequalityConstraints::Bool = true,
                   continuityConstraint::AbstractString="none",
                   continuityParam::Float64 = 0.1,
                   args...)
    g = s.data.network.graph
    paths = s.paths
    tripData = s.trips
    roads = s.data.network.roads

    #Create the model (will be changed to avoid hard-coded parameters)
    # !BarConvTol needs to be changed
    m = Model(solver = MosekSolver(; args...))

    # DECISION VARIABLES
    # Road times
    @variable(m, t[i=vertices(g), j=outneighbors(g,i)] >= s.data.minTimes[i,j])
    # Absolute difference between tripData times and computed times
    @variable(m, epsilon[d=eachindex(tripData)] >= 0)
    @variable(m, T[d=eachindex(tripData)] >= 0)

    if continuityConstraint == "neighborhoods"
        edgesAndNeighbors = flatten(
        [
            [
                (src(edge), dst(edge), src(nearEdge), dst(nearEdge))
                 for nearEdge in findNearEdgesSameType(s.data.network, edge)
            ]
         for edge in edges(g)
        ])
        @variable(m, velocity[(i,j,p,q) = edgesAndNeighbors] >= 0)
    end
    if continuityConstraint == "regularized"
        edgesAndNeighbors = flatten(
        [
            [
                (src(edge), dst(edge), src(nearEdge), dst(nearEdge))
                 for nearEdge in findNearEdgesSameType(s.data.network, edge)
            ]
         for edge in edges(g)
        ])
        @variable(m, velocityDiff[(i,j,p,q) = edgesAndNeighbors] >= 0)
    end
    # OBJECTIVE
    if continuityConstraint == "regularized" # we add a soft constraint for regularized velocities, weighted by road distances
        @objective(m, Min, sum(tripData[d].weight * epsilon[d] for d in eachindex(tripData))
        + continuityParam *sum(velocityDiff[(i,j,p,q)]*2/(roads[i,j].distance + roads[p,q].distance) for (i,j,p,q) in edgesAndNeighbors)
        )
    else
        @objective(m, Min, sum(tripData[d].weight * epsilon[d] for d in eachindex(tripData)))
    end
    # CONSTRAINTS
    # big T constraints
    @constraint(m, pathTime[d=eachindex(tripData)],
        T[d] == sum(paths[d][1][edge] * t[src(edge), dst(edge)] for edge in keys(paths[d][1])))
    # second order cone constraints (define epsilon), equal to time of first path
    @constraint(m, epsLower[d=eachindex(tripData)],
        norm([2 * sqrt(tripData[d].time), T[d] - epsilon[d]])
        <= sum(paths[d][1][edge] * t[src(edge), dst(edge)] for edge in keys(paths[d][1])) + epsilon[d]
        )
    @constraint(m, epsUpper[d=eachindex(tripData)],
        sum(paths[d][1][edge] * t[src(edge), dst(edge)] for edge in keys(paths[d][1])) <=
        epsilon[d] * tripData[d].time
        )

    # inequality constraints
    if inequalityConstraints
        @constraint(m, inequalityPath[d=eachindex(tripData), p=1:(length(paths[d])-1)],
            sum(paths[d][p+1][edge] * t[src(edge), dst(edge)] for edge in keys(paths[d][p+1])) >=
            sum(paths[d][1][edge] * paths[d][1][edge] * t[src(edge), dst(edge)] for edge in keys(paths[d][1]))
            )
    end
    if continuityConstraint == "simple"
        # continuity constraints
        for i in vertices(g), j in outneighbors(g,i)
            for edge in findNearEdgesSameType(s.data.network, Edge(i,j))
                p = src(edge)
                q = dst(edge)
                @constraint(m, t[i,j]/roads[i,j].distance - t[p,q]/roads[p,q].distance
                    <= continuityParam * (roads[i,j].distance + roads[p,q].distance) / 2)
                @constraint(m, t[i,j]/roads[i,j].distance - t[p,q]/roads[p,q].distance
                    >= - continuityParam * (roads[i,j].distance + roads[p,q].distance) / 2)
            end
        end
    elseif continuityConstraint == "regularized"
        # continuity constraints
        for (i,j,p,q) = edgesAndNeighbors
            @constraint(m, t[i,j]/roads[i,j].distance - t[p,q]/roads[p,q].distance
                <= velocityDiff[(i,j,p,q)])
            @constraint(m, t[i,j]/roads[i,j].distance - t[p,q]/roads[p,q].distance
                >= - velocityDiff[(i,j,p,q)])
        end
    elseif continuityConstraint == "neighborhoods"
        # continuity constraints
        clusters = clusterEdges(s.data.network)
        for (k, cluster) in enumerate(clusters)
            for edge in cluster
                i=src(edge); j=dst(edge);
                for nearEdge in findNearEdgesSameType(s.data.network, edge)
                    p=src(nearEdge); q=dst(nearEdge);
                    @constraint(m, t[i,j]/roads[i,j].distance - t[p,q]/roads[p,q].distance
                        <= velocity[(i,j,p,q)] * (roads[i,j].distance + roads[p,q].distance) / 2)
                    @constraint(m, t[i,j]/roads[i,j].distance - t[p,q]/roads[p,q].distance
                        >= - velocity[(i,j,p,q)] * (roads[i,j].distance + roads[p,q].distance) / 2)
                end
            end
            vList = flatten(
            [
                [
                    (src(edge), dst(edge), src(nearEdge), dst(nearEdge))
                     for nearEdge in findNearEdgesSameType(s.data.network, edge)
                ]
                 for edge in cluster
            ])
            @constraint(m, sum(velocity[(i,j,p,q)] for (i,j,p,q) in vList) <= length(vList) * continuityParam)
        end
    end
    # SOLVE SOCP
    status = solve(m)
    # if infeasible, remove possible causes of infeasibility
    if status == :Infeasible && continuityConstraint=="none"
        return socpTimes(s, inequalityConstraints = false)
    end
    times = getvalue(t)

    # Export result as sparse matrix
    result = spzeros(Float64, nv(g), nv(g))
    for i in vertices(g), j in outneighbors(g,i)
        result[i,j] = times[i,j]
    end

    return result
end


text:
using Distributions, Plots, LaTeXStrings; pyplot()

prior(lam) = pdf(TriangularDist(0, 10, 3), lam)
data = [2,1,0,0,1,0,2,2,5,2,4,0,3,2,5,0]

like(lam) = *([pdf(Poisson(lam),x) for x in data]...)
posteriorUpToK(lam) = like(lam)*prior(lam)

delta = 10^-4.
lamRange = 0:delta:10
K = sum([posteriorUpToK(lam)*delta for lam in lamRange])
posterior(lam) = posteriorUpToK(lam)/K

bayesEstimate = sum([lam*posterior(lam)*delta for lam in lamRange])
println("Bayes estimate: ",bayesEstimate)

plot(lamRange, prior.(lamRange),
	c=:blue, label="Prior distribution")
plot!(lamRange, posterior.(lamRange),
	c=:red, label="Posterior distribution",
	xlims=(0, 10), ylims=(0, 1.2),
	xlabel=L"\lambda",ylabel="Density")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using Pkg

Pkg.instantiate()

using PkgDev

try
    version_arg = ARGS[3]
    new_version = nothing

    if version_arg=="Next"
        new_version = nothing
    elseif version_arg=="Major"
        new_version = :major
    elseif version_arg=="Minor"
        new_version = :minor
    elseif version_arg=="Patch"
        new_version = :patch
    else
        new_version = Base.VersionNumber(version_arg)
    end

    PkgDev.tag(PkgDev.FilePathsBase.cwd(), new_version, credentials=ARGS[1], github_username=ARGS[2])

catch err
    Base.display_error(err, catch_backtrace())
end

println("FINISHED")
readline()

max_stars_repo_name: pgunnink/julia-vscode

max_issues_repo_name: pgunnink/julia-vscode

max_forks_repo_name: pgunnink/julia-vscode


text:
"""
    BundleMaster

Bundle method implementation of Lagrangian master method.
This simply uses BundleMethod.jl package.

# Arguments
- `constructor`: struct for bundle method
- `optimizer`: optimization solver for bundle master
- `inner`: bundle method object
"""
mutable struct BundleMaster <: AbstractLagrangeMaster
    constructor::DataType
    optimizer::Union{Nothing,DataType,MOI.OptimizerWithAttributes}
    inner::Union{Nothing,BM.AbstractMethod}
    params::BM.Parameters
    obj_limit::Float64

    function BundleMaster(constructor, optimizer, params = BM.Parameters())
        bm = new()
        bm.constructor = constructor
        bm.optimizer = optimizer
        bm.inner = nothing
        bm.params = params
        bm.obj_limit = Inf
        return bm
    end
end

function load!(method::BundleMaster, num_coupling_variables::Int, num_blocks::Int, eval_function::Function, init_sol::Vector{Float64})
    method.inner = method.constructor(num_coupling_variables, num_blocks, eval_function; init = init_sol, params = method.params)

    # Set optimizer to bundle method
    model = BM.get_jump_model(method.inner)
    JuMP.set_optimizer(model, method.optimizer)

    # This builds the bunlde model.
    BM.build_bundle_model!(method.inner)

    # set objective limit
    BM.set_obj_limit(method.inner, -method.obj_limit)
end

function add_constraints!(LD::AbstractLagrangeDual, method::BundleMaster)
    model = BM.get_jump_model(method.inner)
    λ = model[:x]
    for (id, vars) in LD.block_model.variables_by_couple
        @constraint(model, sum(λ[index_of_λ(LD, v)] for v in vars) == 0)
    end
end

run!(method::BundleMaster) = BM.run!(method.inner)
get_objective(method::BundleMaster) = -BM.get_objective_value(method.inner)
get_solution(method::BundleMaster) = copy(BM.get_solution(method.inner))
get_times(method::BundleMaster)::Vector{Float64} = method.inner.model.time
function set_obj_limit!(method::BundleMaster, val::Float64)
    method.obj_limit = val
end


text:
#Function which compute the Energy and Force of the system

function force_pair(p1,p2,r,Data)

 r6 = (r^2)^3
 r12 = (r^2)^6
 r13 = r12*r
 r7 = r6*r

 dudr1 = -12*(Data.sig12/r13)
 dudr2 = 6*(Data.sig6/r7)

 drdx = -(p1[1]-p2[1])/r
 drdy = -(p1[2]-p2[2])/r

  upair = Data.eps4*(Data.sig12/r12 - Data.sig6/r6)
  fx = (Data.eps4*(dudr1*drdx + dudr2*drdx), Data.eps4*(dudr1*drdy + dudr2*drdy))

 return upair,fx
end


text:
"""
FWSim_Julia contains the created julia code for the FW-Sim
"""
mutable struct FWSim_Julia <: FWSim
    code::String
    fluxes::Vector{String} # defines the order of fluxes which have to be passed to the function, is equal to the flux order in the supplied network
    semu_mass_ratios::Vector{ASTEMUMassRatio} # defines the order of defined semus
    target_mass_ratios::Vector{ASTEMUMassRatio} # defines the order of result semus , is equal to the supplied vector
end

""""
function ast_to_juliacode_01( data::FluxiData , ast::ASTEMUSim , target_smrs::Vector{ASTEMUMassRatio} ; fast_lu=true)

 compiles the ast to juliacode
 creates a julia_fw_sim object with:
-------------------------------------------------------------------------------
 .code has structure:
 "function( semu_mass_ratios::Vector{Float64} ,  v::Vector{Float64} )
     return md
   end"
-------------------------------------------------------------------------------
 the other fields define the order of the two supplied vectors and the result
 vector (see julia_fw_sim struct documentation for details)

 fast_lu : if true, then we use lu-decomp without pivoting. this SHOULD BE
           NO PROBLEM, IF supplied fluxes are NOT VERY CLOSE TO ZERO..
"""
function ast_to_juliacode_01( net::Network , ast::ASTEMUSim , target_smrs::Vector{ASTEMUMassRatio} ; fast_lu=true)
    flag_printDebugOutput = false

    code_header::IOBuffer = IOBuffer()

    # 1. create function "header"
    print(code_header , "function( semu_mass_ratios::Vector{Float64} ,  v::Vector{Float64} )\n\n")

    # here we add all that we find in the X matrix, in any step
    semu_mass_ratios_computed::Set{ASTEMUMassRatio} = Set{ASTEMUMassRatio}()

    # 3. step in, start with processing the sim steps..
    code_computation::IOBuffer = IOBuffer()

    req_constants::Dict{ASTNumNode,String} = Dict{ASTNumNode,String}()
    for siz in sort( collect( keys(ast.sim_steps) ) )
        ss_cnt = 0
        for ss in ast.sim_steps[siz]
            ss_cnt+=1
            print(code_computation, "#----Init A---------------------------------------------------------------------------------\n")
            print( code_computation , to_code_ast_matrix( "sA_$(siz)_$(ss_cnt)" , ss.A , req_constants ) )
            print(code_computation, "#----Init B---------------------------------------------------------------------------\n")
            print( code_computation , to_code_ast_matrix( "sB_$(siz)_$(ss_cnt)" , ss.B , req_constants ) )
            print(code_computation, "#----Init Y---------------------------------------------------------------------------\n")
            print( code_computation , to_code_ast_matrix( "sY_$(siz)_$(ss_cnt)" , ss.Y , req_constants ) )

            print(code_computation, "\n#------COMPUTE!-----------------------------------------------------------------------\n")
            # now add the computation of the matrix Y, and then read out the newly computed values
            if(!fast_lu)
                print( code_computation , "sX_$(siz)_$(ss_cnt) = sA_$(siz)_$(ss_cnt) \\ ( sB_$(siz)_$(ss_cnt) * sY_$(siz)_$(ss_cnt) )\n" )
            end
            if(fast_lu)
                print( code_computation , "( sA_L_$(siz)_$(ss_cnt) , sA_U_$(siz)_$(ss_cnt) ) = lu(sA_$(siz)_$(ss_cnt),Val{false})\n" )
                print( code_computation , "sX_pre_$(siz)_$(ss_cnt) = sA_L_$(siz)_$(ss_cnt) \\ ( sB_$(siz)_$(ss_cnt) * sY_$(siz)_$(ss_cnt) )\n" )
                print( code_computation , "sX_$(siz)_$(ss_cnt) = sA_U_$(siz)_$(ss_cnt) \\ sX_pre_$(siz)_$(ss_cnt) \n\n" )
            end

            if(flag_printDebugOutput)
                s_ma = "sA_$(siz)_$(ss_cnt)"
                s_mb = "sB_$(siz)_$(ss_cnt)"
                s_my = "sY_$(siz)_$(ss_cnt)"
                s_mx = "sX_$(siz)_$(ss_cnt)"

                println( code_computation , "print( \"STEP: A = \\n \") \n" )
                println( code_computation , "print($(s_ma))")
                println( code_computation , "print(\"\\n\")")
                println( code_computation , "print( \"STEP: B = \\n \") \n" )
                println( code_computation , "print($(s_mb))")
                println( code_computation , "print(\"\\n\")")
                println( code_computation , "print( \"STEP: Y = \\n \") \n")
                println( code_computation , "print($(s_my))")
                println( code_computation , "print(\"\\n\")")
                println( code_computation , "print( \"STEP: X = \\n \") \n")
                println( code_computation , "print($(s_mx))")
                println( code_computation , "print(\"\\n\")")
            end

            print(code_computation, "#-------------------------------------------------------------------------------------\n")
            # now read out the newly computed values from X:
            for zib = 1:size(ss.X.m,2)
                for zia = 1:size(ss.X.m,1)
                    smri = ss.X.m[zia,zib]
                    @assert typeof(smri) == ASTEMUMassRatio
                    var_name_si = create_var_name(smri)
                    req_constants[smri] = var_name_si
                    push!( semu_mass_ratios_computed , smri )
                    print(code_computation, "$(var_name_si) = sX_$(siz)_$(ss_cnt)[$(zia),$(zib)]\n")
                end
            end
            print(code_computation, "#-------------------------------------------------------------------------------------\n")
        end
    end

    # 4. now we add the required variable initializations
    # 4.1 collect the required mass distribution ratios:
    code_varinit::IOBuffer = IOBuffer()

    rc_flux = filter( x1 -> typeof(x1)==ASTFluxRate , collect(keys(req_constants)) )
    rc_emr  = filter( x1 -> typeof(x1)==ASTEMUMassRatio , collect(keys(req_constants)) )
    rc_flux = convert(Vector{ASTFluxRate},collect(rc_flux))
    rc_emr  = convert(Vector{ASTEMUMassRatio},collect(rc_emr))

    rc_emr_needed::Vector{ASTEMUMassRatio} = collect( setdiff( Set(rc_emr) , Set(semu_mass_ratios_computed) ) )
    rc_emr_needed                          = sort(rc_emr_needed)
    print("\n\nRequired EMRs: \n");
    print(code_varinit, "#-----Init EMRs-------------------------------------------------------------------------------------\n")
    emri_cnt = 0
    for emri in rc_emr_needed
        emri_cnt += 1
        println("$(emri)")
        print(code_varinit,"$(req_constants[emri]) = semu_mass_ratios[$(emri_cnt)]\n")
    end

    rc_flux = sort(rc_flux) # NOTE: this does not change anything, just makes the resulting code look nicer..
    print("\n\nRequired Fluxes: \n");
    print(code_varinit, "#-----Init Fluxes-----------------------------------------------------------------------------------\n")
    for fluxii in rc_flux
        println("$(fluxii)")
        #fluxpos = find( xfi -> xfi==fluxii.flux , data.net.fluxes )
        fluxpos = find( xfi -> xfi==fluxii.flux , net.fluxes )
        @assert length(fluxpos)==1
        print(code_varinit,"$(req_constants[fluxii]) = v[$(fluxpos[1])]\n")
    end
    print(code_varinit, "#---------------------------------------------------------------------------------------------------\n")

    # 5. collect the results
    code_results::IOBuffer = IOBuffer()
    print(code_results, "#---RESULTS-----------------------------------------------------------------------------------------\n")
    cnt_tsr = 0
    print(code_results , "target_smrs = Array{Float64,1}($(length(target_smrs)))\n")
    for t_semu in target_smrs
        cnt_tsr +=1
        print(code_results , "target_smrs[$(cnt_tsr)] = $(req_constants[t_semu])\n")
    end
    print(code_results, "#---------------------------------------------------------------------------------------------------\n")
    print(code_results,"return target_smrs\n")
    print(code_results, "#---------------------------------------------------------------------------------------------------\n")

    #str_total = string( takebuf_string(code_header) ,  takebuf_string(code_header) ,  takebuf_string(code_computation) )
    code_total = string( String(take!(code_header)) ,  String(take!(code_varinit)) ,  String(take!(code_computation) ) , String(take!(code_results) ) , "end" )

    #return FWSim_Julia(code_total,data.net.fluxes,rc_emr_needed,[])
    return FWSim_Julia(code_total,net.fluxes,rc_emr_needed,target_smrs)
end

function to_code_ast_matrix(matrix_name::String , m::ASTMatrix ,  req_constants::Dict{ASTNumNode,String} )
    sb::IOBuffer = IOBuffer()
    #print(sb,"$(matrix_name)::Array{Float64,2} = Array{Float64,2}($(size(m.m,1)),$(size(m.m,2))) \n")
    print(sb,"$(matrix_name)::Array{Float64,2} = zeros(Float64,$(size(m.m,1)),$(size(m.m,2))) \n")

    # fill the matrix..
    for zj=1:size(m.m,2)
        for zi=1:size(m.m,1)
            # evaluate cell..
            if(isassigned(m.m,zi,zj))
                str_cell = to_code_ast_num(m.m[zi,zj],req_constants)
                print(sb,"$(matrix_name)[$(zi),$(zj)] = $(str_cell) \n")
            end
        end
    end
    return String( take!( sb ) )
end

function to_code_ast_num( n::ASTSymmOp , req_constants::Dict{ASTNumNode,String})
    str = "( "
    op_str = "$(n.op)"
    for zi=1:length(n.args)
        str = string( str , to_code_ast_num(n.args[zi] , req_constants) , " " , op_str )
    end
    str = str[1:length(str)-1]
    str = string(str," ) ")
    return str
end

function to_code_ast_num( n::ASTRealValue , req_constants::Dict{ASTNumNode,String} )
    return " $(n.val) "
end

function to_code_ast_num( n::ASTFluxRate , req_constants::Dict{ASTNumNode,String} )
    if( ! haskey( req_constants , n ) )
        req_constants[n] = create_var_name(n)
    end
    return " $(req_constants[n]) "
end
function to_code_ast_num( n::ASTEMUMassRatio , req_constants::Dict{ASTNumNode,String} )
    if( ! haskey( req_constants , n ) )
        req_constants[n] = create_var_name(n)
    end
    return " $(req_constants[n]) "
end

function create_var_name( n::ASTFluxRate)
    return string("v_",n.flux)
end

function create_var_name( n::ASTEMUMassRatio)
    @assert isinsinglespecies(n.semu)
    ss::String = string("mr_",collect(n.semu)[1].species,"_")
    n_semu_pos = sort( collect( map( x1 -> x1.pos , n.semu ) ) )
    for seip in n_semu_pos
        ss = string( ss , seip , "_" )
    end
    #ss = ss[1:length(ss)-1]
    mass_i =  round(Int , n.mass )
    ss = string( ss , "M_$(mass_i)")
    return String(ss)
end

function benchmark_matrix_solving_01()

    # create lots of toy problems to solve..
    # i.e. diagonal-heavy matrix to invert

    m_sizes        = [100;20;6]
    m_density      = 0.25
    m_diagheavy_a  = 0.1
    m_diagheavy_b  = 1.5
    # A: ms[1] x ms[1] , B: ,s[1] x ms[2] , Y: s[2] x x[3]
    # diagheavy: to the diagonal of m we add
    # dh_a and dh_b times the abs value of the row

    n_tests = 1000

    test_A = Vector{Array{Float64,2}}(n_tests)
    test_B = Vector{Array{Float64,2}}(n_tests)
    test_Y = Vector{Array{Float64,2}}(n_tests)

    test_X_A = Vector{Array{Float64,2}}(n_tests)
    test_X_B = Vector{Array{Float64,2}}(n_tests)
    test_X_C = Vector{Array{Float64,2}}(n_tests)
    test_X_D = Vector{Array{Float64,2}}(n_tests)

    failed_problems = Vector{Array{Array{Float64,2}}}()

    # prepare..
    println("\nPrepare Test:\n")
    for zi=1:n_tests
        Ai = rand(m_sizes[1],m_sizes[1])
        Ai[Ai.<(1.-m_density)] = 0
        for zx=1:size(Ai,1); Ai[zx,zx] = Ai[zx,zx] + m_diagheavy_a  + m_diagheavy_b*sum(abs.(Ai[zx,:])); end
        Bi = rand(m_sizes[1],m_sizes[2])
        Bi[Bi.<(1.-m_density)] = 0
        Yi = rand(m_sizes[2],m_sizes[3])
        Yi[Yi.<(1.-m_density)] = 0

        Ai = Ai * 20 * rand()
        Bi = Bi * 20 * rand()
        Yi = Yi * 20 * rand()

        test_A[zi] = Ai
        test_B[zi] = Bi
        test_Y[zi] = Yi
        if(rem(zi,n_tests/20)==0); print('.'); end
    end

    println("\nRun Test:\n")

    if(false)
        print("Test A: ")
        tic()
        for zi=1:n_tests
            test_X_A[zi] = inv(test_A[zi]) * test_B[zi] * test_Y[zi]
            if(rem(zi,n_tests/20)==0); print('.'); end
        end
        println("")
        toc()
        println("DONE!")

        print("Test B: ")
        tic()
        for zi=1:n_tests
            test_X_A[zi] = test_A[zi] \ ( test_B[zi] * test_Y[zi] )
            if(rem(zi,n_tests/20)==0); print('.'); end
        end
        println("")
        toc()
        println("DONE!")
    end
    if(true)
        print("Test B: ")
        tic()
        for zi=1:n_tests
            test_X_B[zi] = test_A[zi] \ ( test_B[zi] * test_Y[zi] )
            if(rem(zi,n_tests/20)==0); print('.'); end
        end
        println("\n")
        toc()
        println("DONE!")
    end
    if(true)
        print("Test C: ")
        tic()
        for zi=1:n_tests
            (A_L,A_U,A_P) = lu(test_A[zi])
            x = Array{Float64,2}(1,1)

            x_pre = A_L \ ( ( test_B[zi] * test_Y[zi] )[A_P,:] )
            x     = A_U \ x_pre
            test_X_C[zi] = x
            #test_X_B[zi] = test_A[zi] \ ( test_B[zi] * test_Y[zi] )

            if(rem(zi,n_tests/20)==0); print('.'); end
        end
        println("")
        toc()
        println("DONE!")
    end
    if(true)
        print("\nTest D: ")
        tic()
        for zi=1:n_tests
            (A_L,A_U,A_P) = lu(test_A[zi],Val{false})
            x_pre = A_L \ ( test_B[zi] * test_Y[zi] )
            x     = A_U \ x_pre
            test_X_D[zi] = x
            #test_X_B[zi] = test_A[zi] \ ( test_B[zi] * test_Y[zi] )

            if(rem(zi,n_tests/20)==0); print('.'); end
        end
        println("")
        toc()
        println("DONE!")
    end

    # compute max error in D:
    f_dist = (x,y) -> sum( ((x.-y).^2)[:] )
    errors_d = Array{Float64,1}(length(test_X_B))
    for zi=1:length(test_X_B)
        errors_d[zi] = f_dist(test_X_B[zi],test_X_D[zi])
    end
    println("Max Error in D: $(maximum(errors_d))")
end

function test_lu_solving_01()
    a1 = rand(1000,1000)
    b1 = rand(1000,12)

    (a_l,a_u,a_p) = lu(a1)
    x_pre = a_l \ b1[a_p,:]; x = a_u \ x_pre;

    (a2_l,a2_u) = lu(a1,Val{false})
    x_pre_2 = a2_l \ b1; x_2 = a2_u \ x_pre_2;

    print("\n\n")
    #println(a1\b1)
    #println(x)
    #println(x_2)
    x_ref = a1\b1
    f_dist = (x,y) -> sum( ((x.-y).^2)[:] )

    println("Diff 1: $(f_dist(x_ref,x))")
    println("Diff 2: $(f_dist(x_ref,x_2))")
end


text:
module ReinLearn
    abstract type Element end
    abstract type Action end

    struct State{T<:Element}
        values::Array{T}
        State() = new(Array{T,1})
    end

    struct Environment
        state::State
        actions::Vector{Action}
    end

    struct Agent
        policy::Function
        env::Environment
    end
end


text:
module Dec2Part2
function get_lower_upper_bounds(x::String)
    bounds = split(x, '-')
    bounds = parse.(Int64, bounds)
    return bounds
end

function get_char(x::String)
    return x[1]
end

function get_counts(x::String, y::Char)
    local counter = 0
    for c in x
        if c == y
            counter += 1
        end
    end
    return counter
end

function isinbounds(n::Int64, bounds::Array{Int64})
    lower = bounds[1]
    upper = bounds[2]
    return lower ≤ n ≤ upper
end

function get_chars(x::String, bounds::Array{Int64})
    return x[bounds[1]], x[bounds[2]]
end

function only_one_matches(c::Char, lower_c::Char, upper_c::Char)
    if (lower_c == c) & (upper_c != c)
        return true
    elseif (upper_c == c) & (lower_c != c)
        return true
    else
        return false
    end

end

function valid_password(password::String, c::Char, bounds::Array{Int64})
    lower_c, upper_c = get_chars(password, bounds)
    return only_one_matches(c, lower_c, upper_c)
end

function is_valid_password(x::String)
    arr = split(x, ' ')
    arr = String.(arr)
    bounds = get_lower_upper_bounds(arr[1])
    char = get_char(arr[2])
    flag_password = valid_password(arr[3], char, bounds)
    return flag_password
end

function main(input_file::String)::Int64
    lines = readlines(input_file)
    valid_passwords = is_valid_password.(lines)
    return sum(valid_passwords)
end
end


text:
bingos = []

scores = []
deleteBingos = []

function removeValue(val)
    for bingo in bingos
        for i in eachindex(bingo)
            if bingo[i] == val
                bingo[i] = 0
            end
        end
    end
end

function checkComplete(num)
    rows = [0 0 0 0 0]
    cols = [0 0 0 0 0]
    for (i, bingo) in enumerate(bingos)
        if any(x -> x == 0, [sum(bingo, dims=1)..., sum(bingo, dims=2)...])
            s = sum(bingo) - sum((x -> x != 0).(bingo))
            push!(scores, s * (num - 1))
            push!(deleteBingos, i)
        end
    end
    return -1
end

function main()
    toInt(x) = parse(Int, x)
    lines = readlines()
    drawn = toInt.(split(lines[1], ",")) .+ 1
    y = 1
    for line in lines[2:end]
        ints = toInt.(split(line)) .+ 1
        if line == ""
            push!(bingos, zeros(Int, (5, 5)))
            y = 1
            continue
        end
        for i in eachindex(ints)
            bingos[end][y, i] = ints[i]
        end
        y += 1
    end
    for num in drawn
        removeValue(num)
        checkComplete(num)
        for i in reverse(deleteBingos)
            deleteat!(bingos, i)
        end
        empty!(deleteBingos)
    end
    println(scores[1])
    println(scores[end])
end

main()

max_stars_repo_name: LiquidFun/adventofcode

max_issues_repo_name: LiquidFun/adventofcode

max_forks_repo_name: LiquidFun/adventofcode


text:
@testset "AffineConstraint" begin
    v = [_SV(_VI(i)) for i in 1:3]
    ac = Cerberus.AffineConstraint(v[1] + 2.0 * v[2] + 3.0 * v[3], _ET(3.0))
    @test _is_equal(
        ac.f,
        _CSAF([1.0, 2.0, 3.0], [_CVI(1), _CVI(2), _CVI(3)], 0.0),
    )
    @test ac.s == _ET(3.0)
end

function _test_polyhedron(p::Cerberus.Polyhedron)
    @test Cerberus.num_constraints(p) == 2
    @test Cerberus.num_constraints(p, _LT) == length(p.lt_constrs) == 1
    @test Cerberus.num_constraints(p, _GT) == length(p.gt_constrs) == 0
    @test Cerberus.num_constraints(p, _ET) == length(p.et_constrs) == 1
    et_constr = @inferred Cerberus.get_constraint(p, _ET, 1)
    @test et_constr === p.et_constrs[1]
    _is_equal(
        et_constr.f,
        _CSAF([1.0, 2.1, 3.0], [_CVI(1), _CVI(2), _CVI(3)], 0.0),
    )
    @test et_constr.s == _ET(3.0)

    lt_constr = @inferred Cerberus.get_constraint(p, _LT, 1)
    @test lt_constr === p.lt_constrs[1]
    @test _is_equal(lt_constr.f, _CSAF([-3.5, 1.2], [_CVI(1), _CVI(2)], 0.0))
    @test lt_constr.s == _LT(4.0)

    @test p.bounds == [_IN(0.5, 1.0), _IN(-1.3, 2.3), _IN(0.0, 1.0)]

    return nothing
end

@testset "Polyhedron" begin
    p = @inferred _build_polyhedron()
    _test_polyhedron(p)

    # TODO: Test throws on malformed Polyhedron
    @test_throws AssertionError Cerberus.Polyhedron(
        [
            Cerberus.AffineConstraint(
                _CSAF([1.0, 2.0], [_CVI(1), _CVI(2)], 0.0),
                _ET(1.0),
            ),
        ],
        [_IN(0.0, 1.0)],
    )
    @testset "ambient_dim" begin
        @test Cerberus.ambient_dim(p) == 3
    end
    @testset "num_constraints" begin
        @test Cerberus.num_constraints(p) == 2
    end
    @testset "add_variable" begin
        Cerberus.add_variable(p)
        @test Cerberus.ambient_dim(p) == 4
    end
    @testset "empty constructor" begin
        p = @inferred Cerberus.Polyhedron()
        @test Cerberus.ambient_dim(p) == 0
        @test Cerberus.num_constraints(p) == 0
    end
end

@testset "DMIPFormulation" begin
    fm = @inferred _build_dmip_formulation()
    _test_polyhedron(fm._feasible_region)
    @test _is_equal(fm.obj, _CSAF([1.0, -1.0], [_CVI(1), _CVI(2)], 0.0))
    @test isempty(fm.disjunction_formulaters)
    @test [Cerberus.get_variable_kind(fm, _CVI(i)) for i in 1:3] == [_ZO(), nothing, _ZO()]

    @testset "empty constructor" begin
        empty_fm = @inferred Cerberus.DMIPFormulation()
        @test Cerberus.num_variables(empty_fm) == 0
        @test Cerberus.ambient_dim(empty_fm._feasible_region) == 0
        @test Cerberus.num_constraints(empty_fm._feasible_region) == 0
        @test _is_equal(empty_fm.obj, _CSAF())
        @test isempty(empty_fm.disjunction_formulaters)
        @test isempty(empty_fm._variable_kind)
    end

    @testset "num_variables" begin
        @test @inferred Cerberus.num_variables(fm) == 3
    end

    @testset "add_variable" begin
        @inferred Cerberus.add_variable(fm)
        @test Cerberus.num_variables(fm) == 4
    end

    @testset "all_variables" begin
        av = @inferred Cerberus.all_variables(fm)
        @test av == [_CVI(i) for i in 1:4]
        Cerberus.add_variable(fm)
        @test Cerberus.all_variables(fm) == [_CVI(i) for i in 1:5]
    end

    @testset "get_variable_kind" begin
        @test Cerberus.get_variable_kind(fm, _CVI(1)) == _ZO()
        @test Cerberus.get_variable_kind(fm, _CVI(2)) === nothing
        @test Cerberus.get_variable_kind(fm, _CVI(3)) == _ZO()
        @test Cerberus.get_variable_kind(fm, _CVI(4)) === nothing
    end

    @testset "set_variable_kind!" begin
        @inferred Cerberus.set_variable_kind!(fm, _CVI(4), _GI())
        @test Cerberus.get_variable_kind(fm, _CVI(4)) == _GI()
    end

    @testset "num_constraints" begin
        @test @inferred Cerberus.num_constraints(fm) == 2
        @test @inferred Cerberus.num_constraints(fm, _CCI{_LT}) == 1
        @test @inferred Cerberus.num_constraints(fm, _CCI{_GT}) == 0
        @test @inferred Cerberus.num_constraints(fm, _CCI{_ET}) == 1
    end

    @testset "get_constraint(s)" begin
        lt_true = @inferred Cerberus.get_constraint(fm, _CCI{_LT}(1))
        @test _is_equal(lt_true.f, _CSAF([-3.5, 1.2], [_CVI(1), _CVI(2)], 0.0))
        @test lt_true.s == _LT(4.0)
        et_true = @inferred Cerberus.get_constraint(fm, _CCI{_ET}(1))
        @test _is_equal(
            et_true.f,
            _CSAF([1.0, 2.1, 3.0], [_CVI(1), _CVI(2), _CVI(3)], 0.0),
        )
        @test et_true.s == _ET(3.0)
        @test length(Cerberus.get_constraints(fm, _CCI{_LT})) == 1
        @test Cerberus.get_constraints(fm, _CCI{_LT}) == [lt_true]
        @test length(Cerberus.get_constraints(fm, _CCI{_GT})) == 0
        @test length(Cerberus.get_constraints(fm, _CCI{_ET})) == 1
        @test Cerberus.get_constraints(fm, _CCI{_ET}) == [et_true]
    end

    @testset "add_constraint" begin
        f = _CSAF([2.1, 5.4], [_CVI(3), _CVI(4)], 1.2)
        s = _GT(7.2)
        ac = Cerberus.AffineConstraint(f, s)
        @inferred Cerberus.add_constraint(fm, ac)
        @test Cerberus.num_constraints(fm, _CCI{_GT}) == 1
        @test Cerberus.get_constraint(fm, _CCI{_GT}(1)) == ac
        @test Cerberus.get_constraints(fm, _CCI{_GT})[1] == ac
    end

    @testset "isempty" begin
        @test !isempty(fm)
        empty_form = Cerberus.DMIPFormulation()
        @test isempty(empty_form)
        Cerberus.add_constraint(
            empty_form,
            Cerberus.AffineConstraint(_CSAF(), _LT(1.0)),
        )
        @test !isempty(empty_form)
    end

    @testset "get_bounds" begin
        l, u = @inferred Cerberus.get_bounds(fm, _CVI(1))
        @test Cerberus.get_bounds(fm, _CVI(1)) == (0.5, 1.0)
        @test Cerberus.get_bounds(fm, _CVI(2)) == (-1.3, 2.3)
        @test Cerberus.get_bounds(fm, _CVI(3)) == (0.0, 1.0)
        @test Cerberus.get_bounds(fm, _CVI(4)) == (-Inf, Inf)
    end

    @testset "set_bounds!" begin
        @inferred Cerberus.set_bounds!(fm, _CVI(4), _IN(5.7, 8.4))
        @test Cerberus.get_bounds(fm, _CVI(4)) == (5.7, 8.4)
    end

    # TODO: Test throws on malformed DMIPFormulation
end

function _test_gi_polyhedron(p::Cerberus.Polyhedron)
    @test Cerberus.num_constraints(p) == 1
    @test Cerberus.num_constraints(p, _LT) == 1
    @test Cerberus.num_constraints(p, _GT) == 0
    @test Cerberus.num_constraints(p, _ET) == 0
    lt_constr = @inferred Cerberus.get_constraint(p, _LT, 1)
    @test lt_constr === p.lt_constrs[1]
    _is_equal(
        lt_constr.f,
        _CSAF([1.3, 3.7, 2.4], [_CVI(1), _CVI(2), _CVI(3)], 0.0),
    )
    @test lt_constr.f.constant == 0.0
    @test lt_constr.s == _LT(5.5)

    @test p.bounds == [_IN(0.0, 4.5), _IN(0.0, 1.0), _IN(0.0, 3.0)]
end

@testset "General integer polyhedron/formulation" begin
    p = @inferred _build_gi_polyhedron()
    _test_gi_polyhedron(p)

    fm = @inferred _build_gi_dmip_formulation()
    _test_gi_polyhedron(fm._feasible_region)
    @test isempty(fm.disjunction_formulaters)
    @test [Cerberus.get_variable_kind(fm, _CVI(i)) for i in 1:3] == [nothing, _ZO(), _GI()]
end

struct DummyMethod <: DisjunctiveConstraints.AbstractDisjunctiveFormulation end
function Cerberus.new_variables_to_attach(
    ::Cerberus.DisjunctiveFormulater{DummyMethod},
)
    return [_ZO(), _GI(), nothing, nothing, nothing, nothing]
end

@testset "isempty(::DMIPFormulation)" begin
    form = Cerberus.DMIPFormulation()
    @test isempty(form)
    disjunction = Cerberus.Disjunction(
        _CSAF[],
        DisjunctiveConstraints.DisjunctiveSet(
            Matrix{Float64}(undef, 0, 0),
            Matrix{Float64}(undef, 0, 0),
        ),
    )
    formulater = Cerberus.DisjunctiveFormulater(disjunction, DummyMethod())
    Cerberus.attach_formulater!(form, formulater)
    @test !isempty(form)
    form = _build_dmip_formulation()
    @test !isempty(form)
end

@testset "attach_formulater!" begin
    form = _build_dmip_formulation()
    @test Cerberus.num_disjunctive_constraints(form) == 0
    @assert Cerberus.num_variables(form) == 3
    disjunction = Cerberus.Disjunction(
        _CSAF[],
        DisjunctiveConstraints.DisjunctiveSet(
            Matrix{Float64}(undef, 0, 0),
            Matrix{Float64}(undef, 0, 0),
        ),
    )
    formulater = Cerberus.DisjunctiveFormulater(disjunction, DummyMethod())
    @inferred Cerberus.attach_formulater!(form, formulater)
    @test Cerberus.num_variables(form) == 3 + 6
    @test length(form.disjunction_formulaters) == 1
    @test haskey(form.disjunction_formulaters, formulater)
    @test form.disjunction_formulaters[formulater] ==
          [_CVI(i) for i in (3+1):(3+6)]
    @test Cerberus.num_disjunctive_constraints(form) == 1

    @test_throws ArgumentError Cerberus.attach_formulater!(form, formulater)
end


text:
using CUSPARSE
using CUDArt
using Base.Test

m = 25
n = 35
k = 10
blockdim = 5

#################
# test_bsrilu02 #
#################

function test_bsrilu02!(elty)
    A = rand(elty,m,m)
    A += transpose(A)
    A += m * eye(elty,m)
    d_A = CudaSparseMatrixCSR(sparse(A))
    d_A = CUSPARSE.switch2bsr(d_A, convert(Cint,5))
    d_A = CUSPARSE.ilu02!(d_A,'O')
    h_A = to_host(CUSPARSE.switch2csr(d_A))
    Alu = lufact(full(A),Val{false})
    Ac = sparse(Alu[:L]*Alu[:U])
    h_A = ctranspose(h_A) * h_A
    @test_approx_eq(h_A.rowval,Ac.rowval)
    @test reduce(&,isfinite(h_A.nzval))
    d_A = CudaSparseMatrixCSR(sparse(rand(elty,m,n)))
    d_A = CUSPARSE.switch2bsr(d_A, convert(Cint,5))
    @test_throws(DimensionMismatch,CUSPARSE.ilu02!(d_A,'O'))
end

function test_bsrilu02(elty)
    A = rand(elty,m,m)
    A += transpose(A)
    A += m * eye(elty,m)
    d_A = CudaSparseMatrixCSR(sparse(A))
    d_A = CUSPARSE.switch2bsr(d_A, convert(Cint,5))
    d_B = CUSPARSE.ilu02(d_A,'O')
    h_A = to_host(CUSPARSE.switch2csr(d_B))
    Alu = lufact(full(A),Val{false})
    Ac = sparse(Alu[:L]*Alu[:U])
    h_A = ctranspose(h_A) * h_A
    @test_approx_eq(h_A.rowval,Ac.rowval)
    @test reduce(&,isfinite(h_A.nzval))
end

types = [Float32,Float64,Complex64,Complex128]
for elty in types
    tic()
    test_bsrilu02!(elty)
    println("bsrilu02! took ", toq(), " for ", elty)
    tic()
    test_bsrilu02(elty)
    println("bsrilu02 took ", toq(), " for ", elty)
end


text:
# A tokenizer of octal, decimal, hexadecimal and floating point numbers
# =====================================================================

import Automa
import Automa.RegExp: @re_str
const re = Automa.RegExp

# Describe patterns in regular expression.
oct      = re"0o[0-7]+"
dec      = re"[-+]?[0-9]+"
hex      = re"0x[0-9A-Fa-f]+"
prefloat = re"[-+]?([0-9]+\.[0-9]*|[0-9]*\.[0-9]+)"
float    = prefloat | re.cat(prefloat | re"[-+]?[0-9]+", re"[eE][-+]?[0-9]+")
number   = oct | dec | hex | float
numbers  = re.cat(re.opt(number), re.rep(re" +" * number), re" *")

# Register action names to regular expressions.
number.actions[:enter] = [:mark]
oct.actions[:exit]     = [:oct]
dec.actions[:exit]     = [:dec]
hex.actions[:exit]     = [:hex]
float.actions[:exit]   = [:float]

# Compile a finite-state machine.
machine = Automa.compile(numbers)

# This generates a SVG file to visualize the state machine.
# write("numbers.dot", Automa.machine2dot(machine))
# run(`dot -Tpng -o numbers.png numbers.dot`)

# Bind an action code for each action name.
actions = Dict(
    :mark  => :(mark = p),
    :oct   => :(emit(:oct)),
    :dec   => :(emit(:dec)),
    :hex   => :(emit(:hex)),
    :float => :(emit(:float)),
)

# Generate a tokenizing function from the machine.
context = Automa.CodeGenContext()
@eval function tokenize(data::String)
    tokens = Tuple{Symbol,String}[]
    mark = 0
    $(Automa.generate_init_code(context, machine))
    p_end = p_eof = lastindex(data)
    emit(kind) = push!(tokens, (kind, data[mark:p-1]))
    $(Automa.generate_exec_code(context, machine, actions))
    return tokens, cs == 0 ? :ok : cs < 0 ? :error : :incomplete
end

tokens, status = tokenize("1 0x0123BEEF 0o754 3.14 -1e4 +6.022045e23")


text:
module ANTsRegistration

using Images, Glob, Random, Unitful, Suppressor

export register, motioncorr, warp, Global, SyN, MeanSquares, CC, MI, Stage

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("ANTsRegistration not installed properly, run `] build ANTsRegistration', restart Julia and try again")
end
include(depsjl_path)

# Module initialization function
function __init__()
    check_deps()
end

# Per-user working path
function userpath()
    td = tempdir()
    user = ENV["USER"]
    up = joinpath(td, user, "ANTs")
    if !ispath(up)
        mkpath(up)
    end
    return up
end

abstract type AbstractTransformation end

const globaltransformations = ("Translation", "Rigid", "Similarity", "Affine")

struct Global <: AbstractTransformation
    mode::String
    gradientstep::Float64

    function Global(mode::AbstractString, gradientstep::Real)
        cmode = uppercase(mode[1])*lowercase(mode[2:end])
        cmode ∈ globaltransformations || error("mode $mode not recognized")
        return new(cmode, gradientstep)
    end
end
"""
    Global(mode, gradientstep=0.1)

Sets the mode of a global (affine) transformation. `mode` is one of the following strings:
  - "Translation"  (for image shifts only)
  - "Rigid"        (for shifts + rotations)
  - "Similarity"   (for rigid + scaling)
  - "Affine"       (for similarity + shear)
"""
Global(mode::AbstractString) = Global(mode, 0.1)

Base.show(io::IO, st::Global) = print(io, st.mode, '[', st.gradientstep, ']')

struct SyN <: AbstractTransformation
    gradientstep::Float64
    updateFV::Float64
    totalFV::Float64
end
"""
    Syn(gradientstep=0.1, updateFieldVariance=3, totalFieldVariance=0)

Sets the mode to diffeomorphic (warping) registration.
"""
SyN() = SyN(0.1)
SyN(gradientstep) = SyN(gradientstep, 3, 0)

Base.show(io::IO, syn::SyN) = print(io, "SyN[", syn.gradientstep, ',', syn.updateFV, ',', syn.totalFV, ']')

abstract type AbstractMetric end

struct MeanSquares <: AbstractMetric
    metricweight::Float64
end
"""
    MeanSquares(metricweight=1.0)

Compare images using mean-squared difference of voxel values.
"""
MeanSquares() = MeanSquares(1)
metric(cmd::Cmd, fixedname, movingname, m::MeanSquares) =
    `$cmd --metric MeanSquares\[$fixedname,$movingname,$(m.metricweight)\]`

struct MI <: AbstractMetric
    metricweight::Float64
    nhist::Int
    sampling::String
    percentage::Float64
end

"""
    MI(metricweight=1.0, nhist=32, sampling="Regular", percentage=0.25)

Compare images using mutual information.
"""
MI() = MI(1, 32, "Regular", 0.25)
metric(cmd::Cmd, fixedname, movingname, m::MI) =
    `$cmd --metric MI\[$fixedname,$movingname,$(m.metricweight),$(m.nhist),$(m.sampling),$(m.percentage)\]`

struct CC <: AbstractMetric
    metricweight::Float64
    radius::Int
end
"""
    CC(metricweight=1.0, radius=4)

Compare images using neighborhood cross correlation.
"""
CC() = CC(1, 4)
metric(cmd::Cmd, fixedname, movingname, m::CC) =
    `$cmd --metric CC\[$fixedname,$movingname,$(m.metricweight),$(m.radius)\]`

struct Stage{N,T}
    transform::AbstractTransformation
    metric::AbstractMetric
    shrink::NTuple{N,Int}
    smoothing::NTuple{N,T}
    iterations::NTuple{N,Int}
    convergencethresh::Float64
    convergencewindow::Int
end

Stage(img::AbstractArray, transform::AbstractTransformation, metric::AbstractMetric,
      shrink::NTuple{N,Int}, smooth::NTuple{N,Any}, iterations::NTuple{N,Int}) where N =
          Stage(transform, metric, shrink, smooth, iterations, default_convergence(img, transform)[2:3]...)
Stage(img::AbstractArray, transform::AbstractTransformation, metric::AbstractMetric,
      shrink::NTuple{N,Int}, smooth::NTuple{N,Any}) where N =
          Stage(transform, metric, shrink, smooth, default_convergence(img, transform)...)
Stage(img::AbstractArray, transform::AbstractTransformation, metric::AbstractMetric,
      shrink::NTuple{N,Int}) where N =
          Stage(img, transform, metric, shrink, default_smooth(img, transform))
Stage(img::AbstractArray, transform::AbstractTransformation, metric::AbstractMetric) =
    Stage(img, transform, metric, default_shrink(img, transform))

"""
    stage = Stage(fixedimg, transform, metric=MI(), shrink=(8,4,2,1), smooth=(3,2,1,0), iterations=(1000,500,250,0))

Create a registration stage. `transform` is an AbstractTransformation;
currently [`Global`](@ref) and [`Syn`](@ref) are supported. `metric`
is one of [`MI`](@ref), [`MeanSquares`](@ref), or [`CC`](@ref). The
remaining arguments are tuple arguments for multilevel optimization;
each corresponding entry indicates the shrinking factor (8 means 8×
reduced in size), smoothing size (in voxels), and number of iterations
allowed when optimizing.
"""
Stage(img::AbstractArray, transform::AbstractTransformation) =
    Stage(img, transform, MI())

function shcmd(cmd::Cmd, stage::Stage, fixedname, movingname; ismc::Bool=false)
    cmd = `$cmd --transform $(stage.transform)`
    cmd = metric(cmd, fixedname, movingname, stage.metric)
    if ismc
        cmd = `$cmd --iterations $(xstring(stage.iterations))`
        cmd = `$cmd --shrinkFactors $(xstring(stage.shrink)) --smoothingSigmas $(xstring(stage.smoothing))`
    else
        cmd = `$cmd --convergence \[$(xstring(stage.iterations)),$(stage.convergencethresh),$(stage.convergencewindow)\]`
        cmd = `$cmd --shrink-factors $(xstring(stage.shrink)) --smoothing-sigmas $(xqstring(stage.smoothing))$(ustring(stage.smoothing))`
    end
    return cmd
end

default_shrink(img::AbstractArray, transform) = default_shrink(size_spatial(img), transform)
function default_shrink(sz::Dims, transform::Global)
    islarge = any(s->s>=256, sz)
    return islarge ? (12,8,4,2) : (8,4,2,1)
end

function default_shrink(sz::Dims, transform::SyN)
    islarge = any(s->s>=256, sz)
    return islarge ? (10,6,4,2,1) : (8,4,2,1)
end

default_smooth(img::AbstractArray, transform) = default_smooth(size_spatial(img), transform)
function default_smooth(sz::Dims, transform::Global)
    islarge = any(s->s>=256, sz)
    return islarge ? (4,3,2,1) : (3,2,1,0)
end

function default_smooth(sz::Dims, transform::SyN)
    islarge = any(s->s>=256, sz)
    return islarge ? (5,3,2,1,0) : (3,2,1,0)
end

default_convergence(img::AbstractArray, transform) = default_convergence(size_spatial(img), transform)
function default_convergence(sz::Dims, transform::Global)
    islarge = any(s->s>=256, sz)
    return ((1000,500,250,0), 1e-6, 10)
end

function default_convergence(sz::Dims, transform::SyN)
    islarge = any(s->s>=256, sz)
    return islarge ? ((100,100,70,50,0), 1e-6, 10) : ((100,70,50,0), 1e-6, 10)
end

function register(output, nd::Int, fixedname::AbstractString, movingname::AbstractString, pipeline::AbstractVector{<:Stage}; histmatch::Bool=false, winsorize=nothing, verbose::Bool=false, suppressout::Bool=false)
    cmd = `$antsRegistration -d $nd`
    if verbose
        cmd = `$cmd -v`
    end
    if histmatch
        cmd = `$cmd --use-histogram-matching`
    end
    if isa(winsorize, Tuple{Real,Real})
        cmd = `$cmd --winsorize-image-intensities \[$(winsorize[1]),$(winsorize[2])\]`
    end
    for pipe in pipeline
        cmd = shcmd(cmd, pipe, fixedname, movingname)
    end
    cmd = `$cmd --output `
    cmd = shcmd(cmd, output)
    if verbose
        @show cmd
    end
    if suppressout
    @suppress_out run(cmd)
    else
        run(cmd)
    end
end

function register(output, fixed::AbstractArray, moving::AbstractArray, pipeline::AbstractVector{<:Stage}; kwargs...)
    maskfile = ""
    if any(isnan, fixed)
        # Create a mask
        error("not sure how a mask file should be implemented")
    end
    fixedname = write_nrrd(fixed)
    movingname = write_nrrd(moving)
    imgw = register(output, sdims(fixed), fixedname, movingname, pipeline; kwargs...)
    rm(movingname)
    rm(fixedname)
end

function register(fixed::AbstractArray, moving, pipeline::AbstractVector{<:Stage}; kwargs...)
    up = userpath()
    outname = randstring(10)
    tfmname = joinpath(up, outname*"_warp")
    warpedname = joinpath(up, outname*".nrrd")
    output = (tfmname, warpedname)
    register(output, fixed, moving, pipeline; kwargs...)
    imgw = load(warpedname)
    rm(warpedname)
    for tfmfile in glob(outname*"_warp"*"*.mat", up)
        rm(tfmfile)
    end
    for tfmfile in glob(outname*"_warp"*"*.nii.gz", up)
        rm(tfmfile)
    end
    return imgw
end

register(output, fixed::AbstractArray, moving, pipeline::Stage; kwargs...) =
    register(output, fixed, moving, [pipeline]; kwargs...)

"""
    imgw = register(fixed, moving, pipeline; kwargs...)

Return a version of `moving` that has been warped to match
`fixed`. `pipeline` is a single [`Stage`](@ref) or a vector of stages.

# Example

    ## Create some images
    using Images, TestImages, Rotations, CoordinateTransformations
    img = testimage("cameraman")
    tfm = Translation(125,250) ∘ LinearMap(RotMatrix(pi/50)) ∘ Translation(-125,-250)
    img_rotated = warp(img, tfm)
    fixed = img[50:300, 50:300]
    moving = img_rotated[50:300, 50:300]

    ## Perform the registration
    rigid = Stage(fixed, Global("Rigid"))
    syn = Stage(fixed, SyN())
    imgw = register(fixed, moving, [rigid,syn])
"""
register(fixed::AbstractArray, moving, pipeline::Stage; kwargs...) =
    register(fixed, moving, [pipeline]; kwargs...)

"""
    motioncorr(output, fixed, movingname, pipeline; kwargs...)

Perform motion correction in an image series. All images are
registered to `fixed`.  `movingname` is the name of the file that
contains all the individual images, and must be stored in a format
that ITK can handle (see the README). See [`register`](@ref) for more
information about the other arguments.

`output` is generally a tuple of the form `(deformationinfoprefix,
warpedimagefilename)`.

# Example

    output = (joinpath(mydir, "seriesinfo"), joinpath(mydir, "series_warped.nrrd"))
    # Set the number of iterations to be pretty modest (this example uses 2 levels)
    stage = Stage(fixed, Global("Rigid"), CC(), (2,1), (1,0), (15,3))
    motioncorr(output, fixed, movingfilename, stage; verbose=true)
"""
function motioncorr(output, fixed::AbstractArray, movingname::AbstractString, pipeline::AbstractVector{<:Stage}; kwargs...)
    maskfile = ""
    if any(isnan, fixed)
        # Create a mask
        error("not sure how a mask file should be implemented")
    end
    fixedname = write_nrrd(fixed)
    motioncorr(output, sdims(fixed), fixedname, movingname, pipeline; kwargs...)
    rm(fixedname)
end
motioncorr(output, fixed::AbstractArray, movingname::AbstractString, stage::Stage; kwargs...) =
    motioncorr(output, fixed, movingname, [stage]; kwargs...)

function motioncorr(output, nd::Int, fixedname::AbstractString, movingname::AbstractString, pipeline::AbstractVector{<:Stage}; verbose::Bool=false)
    cmd = `$antsMotionCorr -u 1 -e 1 -d $nd`
    if verbose
        cmd = `$cmd -v`
    end
    for pipe in pipeline
        cmd = shcmd(cmd, pipe, fixedname, movingname; ismc=true)
    end
    cmd = `$cmd --output `
    cmd = shcmd(cmd, output)
    if verbose
        @show cmd
    end
    run(cmd)
end

shcmd(cmd::Cmd, str::AbstractString) = `$cmd $str`
function shcmd(cmd::Cmd, strs::Tuple{Vararg{S}}) where S <: AbstractString
    io = IOBuffer()
    print(io, strs[1])
    for i = 2:length(strs)
        print(io, ',', strs[i])
    end
    cmdstr = String(take!(io))
    return `$cmd \[$cmdstr\]`
end

function write_nrrd(img::AbstractArray)
    fn = joinpath(userpath(), randstring(10)*".nrrd")
    save(fn, img)
    return fn
end

function xstring(io::IO, t::Tuple)
    print(io, t[1])
    for i = 2:length(t)
        print(io, 'x', t[i])
    end
end

function xstring(t::Tuple)
    io = IOBuffer()
    xstring(io, t)
    return String(take!(io))
end

xqstring(t) = xstring(to_mm_or_vox.(t))

to_mm_or_vox(x::Number) = x
# to_mm_or_vox(x::Quantity) = float(ustrip(uconvert(u"mm", x)))
to_mm_or_vox(x::Quantity) = float(ustrip(x))  # it seems that ITK treats any physical units as "mm"

ustring(t) = ustring(t[1])
ustring(x::Number) = "vox"
ustring(x::Quantity) = "mm"   # seems to require mm

end # module


text:
x = range(0, stop=π, length=100)
y1 = 0
y2 = @. cos(5x)*exp(-x)
plot(x, y2, lw=0.05)
fill_between!(x, y1, y2, fill="aliceblue")
hline(y1)
xlim(0, π)


text:
function getLinesFromInput(input)
    lines = zeros(Int32, 0)
    splitInput = split.(input, " ")
    for i in splitInput
        (x1, y1) = parse.(Int32, split(i[1],","))
        (x2, y2) = parse.(Int32, split(i[3],","))
        append!(lines, [x1, y1, x2, y2,])
    end
    return transpose(reshape(lines, 4, :))
end

function processInput(filename)
    input = readlines(filename)
    return getLinesFromInput(input)
end

function getSequence(n1, n2)
    if n1 < n2
        return range(n1,n2)
    else
        return range(start=n1, step=-1, stop=n2)
    end
end

function getPointsInLine(line)
    points = zeros(typeof(l[1]),0)
    if line[1] == line[3]
        for i in getSequence(line[2],line[4])
            append!(points, line[1])
            append!(points, i)
        end
    elseif line[2] == line[4]
        for i in getSequence(line[1],line[3])
            append!(points, i)
            append!(points, line[2])
        end
    else # only for part 2, delete the else case for part 1
        xs = getSequence(line[1],line[3])
        ys = getSequence(line[2],line[4])
        for i in 1:length(xs)
            append!(points, xs[i])
            append!(points, ys[i])
        end
     end
     return transpose(reshape(points, 2, :))
end

function fillGrid(lines)
    size = convert(Int64,maximum(lines)) + 1
    grid = reshape(zeros(typeof(lines[1]), size^2), size, :)
    for i in 1:length(lines[:,1])
        points = getPointsInLine(lines[i,:])
        for point in eachrow(points)
            grid[point[1]+1,point[2]+1] += 1
        end
    end
    return grid
end

function findOverlaps(filename)
    lines = processInput(filename)
    return count(x->x>1,fillGrid(lines))
end


text:
# push!(LOAD_PATH,"../src/")

using MIRTio
using Documenter

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

makedocs(;
    modules = [MIRTio],
    authors = "Jeff Fessler <fessler@umich.edu> and contributors",
    repo = "https://github.com/JeffFessler/MIRTio.jl/blob/{commit}{path}#{line}",
    sitename = "MIRTio.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
#       canonical = "https://JeffFessler.github.io/MIRTio.jl/stable",
#       assets = String[],
    ),
    pages = [
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo = "github.com/JeffFessler/MIRTio.jl.git",
    devbranch = "main",
#   devurl = "dev",
#   versions = ["stable" => "v^", "dev" => "dev"],
#   push_preview = true, # only for debug
)


text:
module Merly
import HttpServer.mimetypes
import Base.|

using HttpServer,
      HttpCommon,
      JSON,
      XMLDict

export app, @page, @route, GET,POST,PUT,DELETE,Get,Post,Put,Delete
|(x::ASCIIString, y::ASCIIString)="$x|$y"

q=Dict()

POST="POST"
PUT="PUT"
DELETE="DELETE"
GET="GET"
b=[]
root=pwd()
nfmessage=""
#metodof=1
exten="\"\""::AbstractString

type Q
  query
  params
  body
end

type Pag
  method::Regex
  route::ASCIIString
  code::Function
end

type Fram
  notfound::Function
  start::Function
  use::Function
end

function pag(url::ASCIIString,code::Function)
    Pag(Regex("GET"),url,code)
end

function pag(met::ASCIIString,url::ASCIIString,cod::Function)
    Pag(Regex(met),url,cod)
end

macro page(exp1,exp2)
  quote
    push!(b,pag($exp1,(q,r)->$exp2))
    #nothing
  end
end

macro route(exp1,exp2,exp3)
  quote
    push!(b,pag($exp1,$exp2,(q,r)->$exp3 ))
    #nothing
  end
end

function Get(URL::ASCIIString, fun::Function)
  push!(b,pag("GET",URL,fun))
end

function Post(URL::ASCIIString, fun::Function)
  push!(b,pag("POST",URL,fun))
end

#function Post(URL::ASCIIString, fun::Function)
#  push!(b,pag("POST",URL,fun))
#end

function Put(URL::ASCIIString, fun::Function)
  push!(b,pag("PUT",URL,fun))
end

function Delete(URL::ASCIIString, fun::Function)
  push!(b,pag("DELETE",URL,fun))
end

function _url(ruta::ASCIIString, resource::UTF8String)
  global q
  q.params=Dict()
  q.query=Dict()

  resource = split(resource,"?")
  try
    if length(resource[2])>=1
      q.query=parsequerystring(resource[2])
    end
  end

  resource = split(resource[1],"/")
  resource =resource[2:end]

  ruta = split(ruta,"/")
  ruta =ruta[2:end]

  lruta=length(ruta)
  lresource=length(resource)

  try
    if ruta[end]==""
      lruta=lruta-1
    end

    if resource[end]==""
      lresource=lresource-1
    end
  end

  s=true
  if(lruta==lresource)
    for i=1:lruta
      if length(ruta[i])>=1
        if !(ruta[i][1]==':')
          s=s && (ismatch(Regex(ruta[i]),resource[i]))
        else
          r=ruta[i][2:end]
          q.params[r]=resource[i]
        end
      end
    end
    return s,q
  end
  return false,q
end

function _body(data::Array{UInt8,1},ty::ASCIIString)
  global q
  bo="{}"
  ld=length(data)
  if(ld>=1)
    bo=""
    for i=1:ld
      bo*="$(Char(data[i]))"
    end
    q.body=bo
  end

  if ismatch(Regex("application/json"),ty)
    q.body= JSON.parse(bo)
    return q
  end
  if ismatch(Regex("application/xml"),ty)
    q.body= xml_dict(bo)
    return q
  end

  return q
end

function files(arch::Array{Any,1})
  global root
  for i=1:length(arch)
    roop=replace(replace(arch[i],root,""),"\\","/")
    extencion=split(roop,".")[end]
    if !ismatch(r"(/\.)",roop)
      @page roop begin
      try
        r.headers["Content-Type"]=mimetypes[extencion]
      end
      File(roop[2:end], r)
      end
    end
  end
end

function WebServer(rootte::ASCIIString)
  cd(rootte)
  ls= readdir()
  arrfile=[]
  arrdir=[]
  for i=1:length(ls)
    if isfile(ls[i])
        if (ismatch(Regex("((.)*\\.(?!($exten)))"),ls[i])) && !ismatch(r"^(\.)",ls[i])
          push!(arrfile,normpath(rootte,ls[i]))
        end
    end
    if isdir(ls[i])
      push!(arrdir,normpath(rootte,ls[i]))
    end
  end
  files(arrfile)
  for i=1:length(arrdir)
    WebServer(arrdir[i])
  end
end

function NotFound(res)
  global nfmessage
  res.data = nfmessage
  res.status = 404
end

function File(file::ASCIIString, res::HttpCommon.Response)
  global root
  path = normpath(root, file)
  if isfile(path)
    res.data = readall(path)
  else
    NotFound(res)
  end
end

function process(element::Merly.Pag,q,res::HttpCommon.Response,req::HttpCommon.Request)

  if !(ismatch(Regex("GET"),req.method))
    #println("interpetando los Bytes de req.data como caracteres: ")
    global cors
    try
    body= _body(req.data,req.headers["Accept"])
    catch
    body= _body(req.data,"")
    end
  end
  #h["Content-Type"]="text/html"
  if cors
    res.headers["Access-Control-Allow-Origin"]="*"
    res.headers["Access-Control-Allow-Methods"]="POST,GET,OPTIONS"
  end
  res.status = 200
  #----------aqui escribe el programador-----------
  respond = element.code(q,res)
  #----------------------------------------------
  sal=[]
  try
    sal = matchall(Regex("{{([A-Z]|[a-z]|[0-9])*}}"),respond)
    d=length(sal)
    if d>0
      for i=1:d
        respond= replace(respond,Regex(sal[i]),q.params["$(sal[i][3:end-2])"])
      end
    end
  end
  if length(respond)>0
    try
    res.data= respond
    end
  end
  #res.headers=h
  return res
end

function handler(b::Array{Any,1},req::HttpCommon.Request,res::HttpCommon.Response)
  tam= length(b)
  if tam>0
    for s=1:tam
      pm=ismatch(b[s].method,req.method)
      pasa,q=_url(b[s].route,req.resource)
      if pm && pasa
        resp=process(b[s],q,res,req)
        return resp
      end
    end
  end
  NotFound(res)
  return res
end

function app(r=pwd()::AbstractString,load=""::AbstractString)
global root
global exten
global q
global cors
cors=false::Bool
root=r::AbstractString
q=Q("","","")
if root[end]=='/'
  root=root[1:end-1]
end

if OS_NAME==:Windows
  if root[end]=='\\'
    root=root[1:end-1]
  end
end
if length(load)>0
  if load=="*"
    WebServer(root)
  else
    exten=load::AbstractString
    WebServer(root)
  end
end
cd(root)

  function use(y::AbstractString)
    if(y=="CORS")
      cors=true
    else
      cors=false
    end
  end

  function notfound(x::AbstractString)
    global nfmessage
    nfmessage=x

    try
      path = normpath(root, x)
      if isfile(path)
          nfmessage = readall(path)
      end
    end

    nothing
  end

  function start(host="localhost"::AbstractString,port=8000::Integer)
    http = HttpHandler((req, res)-> handler(b,req,res))
    http.events["error"]  = (client, error) -> println(error)
    http.events["listen"] = (port)          -> println("Listening on $port...")
    server = Server(http)
    if host=="localhost"
      host="127.0.0.1"
    end
    try
      IPv4(host)
      #@async run(server, host=IPv4(host), port=port)
      run(server, host=IPv4(host), port=port)
    catch
      "only IPv4 addresses"
    end
  end
  return Fram(notfound,start,use)
end
end # module