text:
module LighthouseFlux

using Zygote: Zygote
using Flux: Flux
using Lighthouse: Lighthouse, classes, log_resource_info!, log_value!

export FluxClassifier

#####
##### `FluxClassifier`
#####

struct FluxClassifier{M,O,C,P,OH,OC} <: Lighthouse.AbstractClassifier
    model::M
    optimiser::O
    classes::C
    params::P
    onehot::OH
    onecold::OC
end

"""
    FluxClassifier(model, optimiser, classes; params=Flux.params(model),
                   onehot=(label -> Flux.onehot(label, 1:length(classes))),
                   onecold=(label -> Flux.onecold(label, 1:length(classes))))

Return a `FluxClassifier <: Lighthouse.AbstractClassifier` with the given arguments:

- `model`: a Flux model. The model must additionally support LighthouseFlux's [`loss`](@ref)
  and [`loss_and_prediction`](@ref) functions.

- `optimiser`: a [Flux optimiser](https://fluxml.ai/Flux.jl/stable/training/optimisers/)

- `classes`: a `Vector` or `Tuple` of possible class values; this is the return
value of `Lighthouse.classes(::FluxClassifier)`.

- `params`: The parameters to optimise during training; generally, a `Zygote.Params`
value or a value that can be passed to `Zygote.Params`.

- `onehot`: the function used to convert hard labels to soft labels when
`Lighthouse.onehot` is called with this classifier.

- `onecold`: the function used to convert soft labels to hard labels when
`Lighthouse.onecold` is called with this classifier.
"""
function FluxClassifier(model, optimiser, classes; params=Flux.params(model),
                        onehot=(label -> Flux.onehot(label, 1:length(classes))),
                        onecold=(label -> Flux.onecold(label, 1:length(classes))))
    return FluxClassifier(Flux.testmode!(model), optimiser, classes, params, onehot,
                          onecold)
end

"""
    loss(model, batch_arguments...)

Return the scalar loss of `model` given `batch_arguments`.

This method must be implemented for all `model`s passed to [`FluxClassifier`](@ref).
"""
function loss end

"""
    loss_and_prediction(model, input_batch, other_batch_arguments...)

Return `(model_loss, model_prediction)` where:

- `model_loss` is equivalent to (and defaults to) `loss(model, input_batch, other_batch_arguments...)`.

- `model_prediction` is a matrix where the `i`th column is the soft label prediction for the `i`th
sample in `input_batch`. Thus, the numnber of columns should be `size(input_batch)[end]`, while the
number of rows is equal to the number of possible classes predicted by model. `model_prediction`
defaults to `model(input_batch)`.

This method must be implemented for all `model`s passed to [`FluxClassifier`](@ref), but
has the default return values described above, so it only needs to be overloaded if the
default definitions do not yield the expected values for a given `model` type. It
additionally may be overloaded to avoid redundant computation if `model`'s loss
function computes soft labels as an intermediate result.
"""
function loss_and_prediction(model, input_batch, other_batch_arguments...)
    return (loss(model, input_batch, other_batch_arguments...), model(input_batch))
end

#####
##### Lighthouse `AbstractClassifier` Interface
#####

Lighthouse.classes(classifier::FluxClassifier) = classifier.classes

function Lighthouse.is_early_stopping_exception(::FluxClassifier, exception)
    return exception isa Flux.Optimise.StopException
end

Lighthouse.onehot(classifier::FluxClassifier, label) = classifier.onehot(label)

Lighthouse.onecold(classifier::FluxClassifier, label) = classifier.onecold(label)

function Lighthouse.train!(classifier::FluxClassifier, batches, logger)
    Flux.trainmode!(classifier.model)
    weights = Zygote.Params(classifier.params)
    for batch in batches
        train_loss, back = log_resource_info!(logger, "train/forward_pass";
                                              suffix="_per_batch") do
            f = () -> loss(classifier.model, batch...)
            return Zygote.pullback(f, weights)
        end
        log_value!(logger, "train/loss_per_batch", train_loss)
        gradients = log_resource_info!(logger, "train/reverse_pass"; suffix="_per_batch") do
            return back(Zygote.sensitivity(train_loss))
        end
        log_resource_info!(logger, "train/update"; suffix="_per_batch") do
            Flux.Optimise.update!(classifier.optimiser, weights, gradients)
            return nothing
        end
    end
    Flux.testmode!(classifier.model)
    return nothing
end

function Lighthouse.loss_and_prediction(classifier::FluxClassifier, batch...)
    return Flux.cpu(loss_and_prediction(classifier.model, batch...))
end

#####
##### `CuIterator`
#####

Base.@deprecate_moved CuIterator CuArrays false

#####
##### miscellaneous utilities
#####

"""
    evaluate_chain_in_debug_mode(chain::Flux.Chain, input)

Evaluate `chain(input)`, printing additional debug information at each layer.
"""
function evaluate_chain_in_debug_mode(chain::Flux.Chain, input)
    for (i, layer) in enumerate(chain)
        @info "Executing layer $i / $(length(chain))..." layer size(input)
        input = layer(input)
        @info output_size = size(input)
    end
    return input
end

end # module


text:
# a collection of helper functions

dropcols(df::AbstractDataFrame, cols) = df[setdiff(names(df), cols)]

function tomissing!(d::AbstractDataFrame,x::Union{Integer,Symbol})
	allowmissing!(d,x)
	d[findall(isnan,d[x]),x] = missing
end

function alltomissing!(d::AbstractDataFrame)
	for nm in names(d)
		if eltype(d[nm]) <: Float64
			if any(isnan,d[nm])
				tomissing!(d,nm)
			end
		end
	end
end

function disallowmissing!(d::AbstractDataFrame,x::Union{Integer,Symbol})
	d[x] = convert(Vector{eltype(d[x])},d[x])
	return d
end

function disallowmissing!(d,x...)
	for y in x
		disallowmissing!(d,y)
	end
	return d
end

ext: jl


text:
using DataArrays
using Base: @get!, promote_eltype
using Base.Broadcast: bitcache_chunks, bitcache_size, dumpbitcache
using Compat: promote_eltype_op

if isdefined(Base, :OneTo)
    if VERSION < v"0.6.0-dev.1121"
        _broadcast_shape(x...) = Base.to_shape(Base.Broadcast.broadcast_shape(x...))
    else
        _broadcast_shape(x...) = Base.to_shape(Base.Broadcast.broadcast_indices(x...))
    end
else
    const _broadcast_shape = Base.Broadcast.broadcast_shape
end

# Check that all arguments are broadcast compatible with shape
# Differs from Base in that we check for exact matches
function check_broadcast_shape(shape::Dims, As::(@compat Union{AbstractArray,Number})...)
    samesize = true
    for A in As
        if ndims(A) > length(shape)
            throw(DimensionMismatch("cannot broadcast array to have fewer dimensions"))
        end
        for k in 1:length(shape)
            n, nA = shape[k], size(A, k)
            samesize &= (n == nA)
            if n != nA != 1
                throw(DimensionMismatch("array could not be broadcast to match destination"))
            end
        end
    end
    samesize
end

# Get ref for value for a PooledDataArray, adding to the pool if
# necessary
_unsafe_pdaref!(Bpool, Brefdict::Dict, val::NAtype) = 0
function _unsafe_pdaref!(Bpool, Brefdict::Dict, val)
    @get! Brefdict val begin
        push!(Bpool, val)
        length(Bpool)
    end
end

# Generate a branch for each possible combination of NA/not NA. This
# gives good performance at the cost of 2^narrays branches.
function gen_na_conds(f, nd, arrtype, outtype, daidx=find([arrtype...] .!= AbstractArray), pos=1, isna=())
    if pos > length(daidx)
        args = Any[Symbol("v_$(k)") for k = 1:length(arrtype)]
        for i = 1:length(daidx)
            if isna[i]
                args[daidx[i]] = NA
            end
        end

        # Needs to be gensymmed so that the compiler won't box it
        val = gensym("val")
        quote
            $val = $(Expr(:call, f, args...))
            $(if outtype == DataArray
                :(@inbounds unsafe_dasetindex!(Bdata, Bc, $val, ind))
            elseif outtype == PooledDataArray
                :(@inbounds (@nref $nd Brefs i) = _unsafe_pdaref!(Bpool, Brefdict, $val))
            end)
        end
    else
        k = daidx[pos]
        quote
            if $(Symbol("isna_$(k)"))
                $(gen_na_conds(f, nd, arrtype, outtype, daidx, pos+1, tuple(isna..., true)))
            else
                $(if arrtype[k] == DataArray
                    :(@inbounds $(Symbol("v_$(k)")) = $(Symbol("data_$(k)"))[$(Symbol("state_$(k)_0"))])
                else
                    :(@inbounds $(Symbol("v_$(k)")) = $(Symbol("pool_$(k)"))[$(Symbol("r_$(k)"))])
                end)
                $(gen_na_conds(f, nd, arrtype, outtype, daidx, pos+1, tuple(isna..., false)))
            end
        end
    end
end

# Broadcast implementation for DataArrays
#
# TODO: Fall back on faster implementation for same-sized inputs when
# it is safe to do so.
function gen_broadcast_dataarray(nd::Int, arrtype::@compat(Tuple{Vararg{DataType}}), outtype, f::Function)
    F = Expr(:quote, f)
    narrays = length(arrtype)
    As = [Symbol("A_$(i)") for i = 1:narrays]
    dataarrays = find([arrtype...] .== DataArray)
    abstractdataarrays = find([arrtype...] .!= AbstractArray)
    have_fastpath = outtype == DataArray && all(x->!(x <: PooledDataArray), arrtype)

    @eval let
        function _F_(B::$(outtype), $(As...))
            @assert ndims(B) == $nd

            # Set up input DataArray/PooledDataArrays
            $(Expr(:block, [
                arrtype[k] == DataArray ? quote
                    $(Symbol("na_$(k)")) = $(Symbol("A_$(k)")).na.chunks
                    $(Symbol("data_$(k)")) = $(Symbol("A_$(k)")).data
                    $(Symbol("state_$(k)_0")) = $(Symbol("state_$(k)_$(nd)")) = 1
                    @nexprs $nd d->($(Symbol("skip_$(k)_d")) = size($(Symbol("data_$(k)")), d) == 1)
                end : arrtype[k] == PooledDataArray ? quote
                    $(Symbol("refs_$(k)")) = $(Symbol("A_$(k)")).refs
                    $(Symbol("pool_$(k)")) = $(Symbol("A_$(k)")).pool
                end : nothing
            for k = 1:narrays]...))

            # Set up output DataArray/PooledDataArray
            $(if outtype == DataArray
                quote
                    Bdata = B.data
                    # Copy in case aliased
                    # TODO: check for aliasing?
                    Bna = falses(size(Bdata))
                    Bc = Bna.chunks
                    ind = 1
                end
            elseif outtype == PooledDataArray
                quote
                    Bpool = B.pool = similar(B.pool, 0)
                    Brefs = B.refs
                    Brefdict = Dict{eltype(Bpool),eltype(Brefs)}()
                end
            end)

            @nloops($nd, i, $(outtype == DataArray ? (:Bdata) : (:Brefs)),
                # pre
                d->($(Expr(:block, [
                    arrtype[k] == DataArray ? quote
                        $(Symbol("state_$(k)_")){d-1} = $(Symbol("state_$(k)_d"));
                        $(Symbol("j_$(k)_d")) = $(Symbol("skip_$(k)_d")) ? 1 : i_d
                    end : quote
                        $(Symbol("j_$(k)_d")) = size($(Symbol("A_$(k)")), d) == 1 ? 1 : i_d
                    end
                for k = 1:narrays]...))),

                # post
                d->($(Expr(:block, [quote
                    $(Symbol("skip_$(k)_d")) || ($(Symbol("state_$(k)_d")) = $(Symbol("state_$(k)_0")))
                end for k in dataarrays]...))),

                # body
                begin
                    # Advance iterators for DataArray and determine NA status
                    $(Expr(:block, [
                        arrtype[k] == DataArray ? quote
                            @inbounds $(Symbol("isna_$(k)")) = Base.unsafe_bitgetindex($(Symbol("na_$(k)")), $(Symbol("state_$(k)_0")))
                        end : arrtype[k] == PooledDataArray ? quote
                            @inbounds $(Symbol("r_$(k)")) = @nref $nd $(Symbol("refs_$(k)")) d->$(Symbol("j_$(k)_d"))
                            $(Symbol("isna_$(k)")) = $(Symbol("r_$(k)")) == 0
                        end : nothing
                    for k = 1:narrays]...))

                    # Extract values for ordinary AbstractArrays
                    $(Expr(:block, [
                        :(@inbounds $(Symbol("v_$(k)")) = @nref $nd $(Symbol("A_$(k)")) d->$(Symbol("j_$(k)_d")))
                    for k = find([arrtype...] .== AbstractArray)]...))

                    # Compute and store return value
                    $(gen_na_conds(F, nd, arrtype, outtype))

                    # Increment state
                    $(Expr(:block, [:($(Symbol("state_$(k)_0")) += 1) for k in dataarrays]...))
                    $(if outtype == DataArray
                        :(ind += 1)
                    end)
                end)

            $(if outtype == DataArray
                :(B.na = Bna)
            end)
        end
        _F_
    end
end

datype(A_1::PooledDataArray, As...) = tuple(PooledDataArray, datype(As...)...)
datype(A_1::DataArray, As...) = tuple(DataArray, datype(As...)...)
datype(A_1, As...) = tuple(AbstractArray, datype(As...)...)
datype() = ()

datype_int(A_1::PooledDataArray, As...) = (@compat(UInt64(2)) | (datype_int(As...) << 2))
datype_int(A_1::DataArray, As...) = (@compat(UInt64(1)) | (datype_int(As...) << 2))
datype_int(A_1, As...) = (datype_int(As...) << 2)
datype_int() = @compat UInt64(0)

# The following four methods are to avoid ambiguity warnings on 0.4
Base.map!(f::Base.Callable, B::DataArray) =
    invoke(map!, Tuple{Base.Callable, AbstractArray}, f, B)
Base.map!(f::Base.Callable, B::PooledDataArray) =
    invoke(map!, Tuple{Base.Callable, AbstractArray}, f, B)
Base.broadcast!(f::Base.Function, B::DataArray) =
    invoke(map!, Tuple{Base.Callable, AbstractArray}, f, B)
Base.broadcast!(f::Base.Function, B::PooledDataArray) =
    invoke(map!, Tuple{Base.Callable, AbstractArray}, f, B)

for bsig in (DataArray, PooledDataArray), asig in (Union{Array,BitArray,Number},DataArray, PooledDataArray,)
    @eval let cache = Dict{Function,Dict{UInt64,Dict{Int,Function}}}()
        function Base.map!(f::Base.Callable, B::$bsig, As::$asig...)
            nd = ndims(B)
            length(As) <= 8 || throw(ArgumentError("too many arguments"))
            samesize = check_broadcast_shape(size(B), As...)
            samesize || throw(DimensionMismatch("Argument dimensions must match"))
            arrtype = datype_int(As...)

            cache_f    = @get! cache      f        Dict{UInt64,Dict{Int,Function}}()
            cache_f_na = @get! cache_f    arrtype  Dict{Int,Function}()
            func       = @get! cache_f_na nd       gen_broadcast_dataarray(nd, datype(As...), $bsig, f)

            func(B, As...)
            B
        end

        function Base.broadcast!(f::Function, B::$bsig, As::$asig...)
            nd = ndims(B)
            length(As) <= 8 || throw(ArgumentError("too many arguments"))
            samesize = check_broadcast_shape(size(B), As...)
            arrtype = datype_int(As...)

            cache_f    = @get! cache      f        Dict{UInt64,Dict{Int,Function}}()
            cache_f_na = @get! cache_f    arrtype  Dict{Int,Function}()
            func       = @get! cache_f_na nd       gen_broadcast_dataarray(nd, datype(As...), $bsig, f)

            # println(code_typed(func, typeof(tuple(B, As...))))
            func(B, As...)
            B
        end
    end
end

databroadcast(f::Function, As...) = broadcast!(f, DataArray(promote_eltype(As...), _broadcast_shape(As...)), As...)
pdabroadcast(f::Function, As...) = broadcast!(f, PooledDataArray(promote_eltype(As...), _broadcast_shape(As...)), As...)

function exreplace!(ex::Expr, search, rep)
    for i = 1:length(ex.args)
        if ex.args[i] == search
            splice!(ex.args, i, rep)
            break
        else
            exreplace!(ex.args[i], search, rep)
        end
    end
    ex
end
exreplace!(ex, search, rep) = ex

macro da_broadcast_vararg(func)
    if (func.head != :function && func.head != :(=)) ||
       func.args[1].head != :call || !isa(func.args[1].args[end], Expr) ||
       func.args[1].args[end].head != :...
        throw(ArgumentError("@da_broadcast_vararg may only be applied to vararg functions"))
    end

    va = func.args[1].args[end]
    defs = Any[]
    for n = 1:4, aa = 0:n-1
        def = deepcopy(func)
        rep = Any[Symbol("A_$(i)") for i = 1:n]
        push!(rep, va)
        exreplace!(def.args[2], va, rep)
        rep = Vector{Any}(n+1)
        for i = 1:aa
            rep[i] = Expr(:(::), Symbol("A_$i"), AbstractArray)
        end
        for i = aa+1:n
            rep[i] = Expr(:(::), Symbol("A_$i"), (@compat Union{DataArray, PooledDataArray}))
        end
        rep[end] = Expr(:..., Expr(:(::), va.args[1], AbstractArray))
        exreplace!(def.args[1], va, rep)
        push!(defs, def)
    end
    esc(Expr(:block, defs...))
end

macro da_broadcast_binary(func)
    if (func.head != :function && func.head != :(=)) ||
       func.args[1].head != :call ||
       length(func.args[1].args) != 3
        throw(ArgumentError("@da_broadcast_binary may only be applied to two-argument functions"))
    end
    (ff, A, B) = func.args[1].args
    f = esc(ff)
    body = func.args[2]
    quote
        ($f)($A::(@compat Union{DataArray, PooledDataArray}), $B::(@compat Union{DataArray, PooledDataArray})) = $(body)
        ($f)($A::(@compat Union{DataArray, PooledDataArray}), $B::AbstractArray) = $(body)
        ($f)($A::AbstractArray, $B::(@compat Union{DataArray, PooledDataArray})) = $(body)
    end
end

# Broadcasting DataArrays returns a DataArray
@da_broadcast_vararg Base.broadcast(f::Function, As...) = databroadcast(f, As...)

# Definitions for operators,
(.*)(A::BitArray, B::(@compat Union{DataArray{Bool}, PooledDataArray{Bool}})) = databroadcast(*, A, B)
(.*)(A::(@compat Union{DataArray{Bool}, PooledDataArray{Bool}}), B::BitArray) = databroadcast(*, A, B)
@da_broadcast_vararg (.*)(As...) = databroadcast(*, As...)
@da_broadcast_binary (.%)(A, B) = databroadcast(%, A, B)
@da_broadcast_vararg (.+)(As...) = broadcast!(+, DataArray(promote_eltype_op(@functorize(+), As...), _broadcast_shape(As...)), As...)
@da_broadcast_binary (.-)(A, B) =
    broadcast!(-, DataArray(promote_op(@functorize(-), eltype(A), eltype(B)),
                            _broadcast_shape(A,B)), A, B)
@da_broadcast_binary (./)(A, B) =
    broadcast!(/, DataArray(promote_op(@functorize(/), eltype(A), eltype(B)),
                            _broadcast_shape(A, B)), A, B)
@da_broadcast_binary (.\)(A, B) =
    broadcast!(\, DataArray(promote_op(@functorize(\), eltype(A), eltype(B)),
                            _broadcast_shape(A, B)), A, B)
(.^)(A::(@compat Union{DataArray{Bool}, PooledDataArray{Bool}}), B::(@compat Union{DataArray{Bool}, PooledDataArray{Bool}})) = databroadcast(>=, A, B)
(.^)(A::BitArray, B::(@compat Union{DataArray{Bool}, PooledDataArray{Bool}})) = databroadcast(>=, A, B)
(.^)(A::AbstractArray{Bool}, B::(@compat Union{DataArray{Bool}, PooledDataArray{Bool}})) = databroadcast(>=, A, B)
(.^)(A::(@compat Union{DataArray{Bool}, PooledDataArray{Bool}}), B::BitArray) = databroadcast(>=, A, B)
(.^)(A::(@compat Union{DataArray{Bool}, PooledDataArray{Bool}}), B::AbstractArray{Bool}) = databroadcast(>=, A, B)
@da_broadcast_binary (.^)(A, B) =
    broadcast!(^, DataArray(promote_op(@functorize(^), eltype(A), eltype(B)),
                            _broadcast_shape(A, B)), A, B)

# XXX is a PDA the right return type for these?
Base.broadcast(f::Function, As::PooledDataArray...) = pdabroadcast(f, As...)
(.*)(As::PooledDataArray...) = pdabroadcast(*, As...)
(.%)(A::PooledDataArray, B::PooledDataArray) = pdabroadcast(%, A, B)
(.+)(As::PooledDataArray...) =
    broadcast!(+, PooledDataArray(promote_eltype_op(@functorize(+), As...), _broadcast_shape(As...)), As...)
(.-)(A::PooledDataArray, B::PooledDataArray) =
    broadcast!(-, PooledDataArray(promote_op(@functorize(-), eltype(A), eltype(B)),
                                  _broadcast_shape(A,B)), A, B)
(./)(A::PooledDataArray, B::PooledDataArray) =
    broadcast!(/, PooledDataArray(promote_op(@functorize(/), eltype(A), eltype(B)),
                                  _broadcast_shape(A, B)), A, B)
(.\)(A::PooledDataArray, B::PooledDataArray) =
    broadcast!(\, PooledDataArray(promote_op(@functorize(\), eltype(A), eltype(B)),
                                  _broadcast_shape(A, B)), A, B)
(.^)(A::PooledDataArray{Bool}, B::PooledDataArray{Bool}) = databroadcast(>=, A, B)
(.^)(A::PooledDataArray, B::PooledDataArray) =
    broadcast!(^, PooledDataArray(promote_op(@functorize(^), eltype(A), eltype(B)),
                                  _broadcast_shape(A, B)), A, B)

for (sf, vf) in zip(scalar_comparison_operators, array_comparison_operators)
    @eval begin
        # ambiguity
        $(vf)(A::(@compat Union{PooledDataArray{Bool},DataArray{Bool}}), B::(@compat Union{PooledDataArray{Bool},DataArray{Bool}})) =
            broadcast!($sf, DataArray(Bool, _broadcast_shape(A, B)), A, B)
        $(vf)(A::(@compat Union{PooledDataArray{Bool},DataArray{Bool}}), B::AbstractArray{Bool}) =
            broadcast!($sf, DataArray(Bool, _broadcast_shape(A, B)), A, B)
        $(vf)(A::AbstractArray{Bool}, B::(@compat Union{PooledDataArray{Bool},DataArray{Bool}})) =
            broadcast!($sf, DataArray(Bool, _broadcast_shape(A, B)), A, B)

        @da_broadcast_binary $(vf)(A, B) = broadcast!($sf, DataArray(Bool, _broadcast_shape(A, B)), A, B)
    end
end


text:
"""
  is_inside_outline(fo,outline)

Checks whether filling object `fo` is inside the `outline` object. Special methods
called for circle, rectangle, and polygon outlines.
"""
function is_inside_outline(fo::FillingCircle{T},outline::OutlineCircle{T}) where T<:Real
  furthest_point = fo.radius + sqrt((fo.center.x-outline.center.x)^2 + (fo.center.y-outline.center.y)^2)
  return (furthest_point < outline.radius ? true : false)
end
function is_inside_outline(fo::FillingCircle{T},outline::OutlineRectangle{T}) where T<:Real
  if (fo.center.y + fo.radius) > outline.p2.y         # check top
      return false
  elseif (fo.center.y - fo.radius) < outline.p1.y     # check bottom
      return false
  end

  if (fo.center.x + fo.radius) > outline.p2.x         # check right
      return false
  elseif (fo.center.x - fo.radius) < outline.p1.x     # check left
      return false
  end

  return true
end
function is_inside_outline(fo::FillingCircle{T},outline::OutlinePolygon{T}) where T<:Real
  return is_inside_polygon(outline.pList, fo.center)
end

function is_inside_polygon(polygon::Vector{P}, p::P; extreme = Point(100_000.0, p.y)) where P<:Point
  n = length(polygon)
  # There must be at least 3 vertices in polygon
  if (n < 3); return false; end

  # Count intersections of the above line with sides of polygon
  count = 0
  for i=1:n
    next = mod(i,n)+1
    # Check if the line segment from 'p' to 'extreme' intersects
    # with the line segment from 'polygon[i]' to 'polygon[next]'
    if (do_intersect(polygon[i], polygon[next], p, extreme))
      # If the point 'p' is colinear with line segment 'i-next', then check if it lies on segment.
      # If it lies, return true, otherwise false
      if (get_orientation(polygon[i], p, polygon[next]) == 0)
         return on_segment(polygon[i], p, polygon[next])
      end
      count += 1
    end
  end

  return (isodd(count) ? true : false)
end

function is_intersecting_walls(outline::O,fo::F) where {O<:AbstractOutlineObject,F<:AbstractFillingObject}
  for tw=1:length(outline.wlist)
    pointOnWall = get_nearest_point(fo.center,outline.wlist[tw])
    if is_point_in_wall(pointOnWall,outline.wlist[tw])
      Δx = pointOnWall.x - fo.center.x
      Δy = pointOnWall.y - fo.center.y
      dist = sqrt(Δx^2 + Δy^2)

      if dist < fo.radius
        return true
      end
    end
  end

  return false
end

"""
  do_intersect(p1,q1,p2,q2)

Checks if the line segment p1 to q1 intersects line from p2 to q2. Returns
a bool.
"""
function do_intersect(p1::Point{T},q1::Point{T},p2::Point{T},q2::Point{T}) where T<:Real
    # Find the four orientations needed for general case
    o1 = get_orientation(p1, q1, p2)
    o2 = get_orientation(p1, q1, q2)
    o3 = get_orientation(p2, q2, p1)
    o4 = get_orientation(p2, q2, q1)

    # General case
    if (o1 != o2 && o3 != o4)
        return true
    else
        return false
    end
end

"""
  get_orientation(p,q,r)

Checks orientation of polygon p-q-r. Returns 0 for colinear,
1 for clockwise, and 2 for counter-clockwise.
"""
function get_orientation(p::Point{T},q::Point{T},r::Point{T}) where T<:Real
    val = (q.y-p.y)*(r.x-q.x) - (q.x-p.x)*(r.y-q.y);

    if isapprox(val, 0; atol=eps(Float64), rtol=0)
      return zero(Int)
    end

    return (val > 0) ? one(Int) : 2*one(Int) # clock- or counterclock-wise
end

function get_nearest_point!(point::Point{T},node::Point{T},wall::LineWall) where T<:Real
  px=node.x; py=node.y

  Ax=wall.nodes[1].x; Ay=wall.nodes[1].y
  Bx=wall.nodes[2].x; By=wall.nodes[2].y

  bx=px-Ax; by=py-Ay
  ax=Bx-Ax; ay=By-Ay

  ℓ2 = ax^2+ay^2

  dotprod = ax*bx + ay*by

  point.x = dotprod*ax/ℓ2 + Ax
  point.y = dotprod*ay/ℓ2 + Ay

  nothing
end
function get_nearest_point!(point::Point{T},node::Point{T},wall::CircleWall) where T<:Real
  px=node.x; py=node.y
  cx=wall.center.x; cy=wall.center.y
  r = wall.radius
  θ = atan(py-cy,px-cx)

  point.x = cx + r*cos(θ)
  point.y = cy + r*sin(θ)

  nothing
end
function get_nearest_point(node::Point{T},wall::W) where {T<:Real,W<:AbstractWall}
  point = Point(0.0,0.0)
  get_nearest_point!(point,node,wall)
  return point
end

"""
  is_point_in_wall(point,wall)

Checks whether `point` is on `wall`. Different methods for handling
circle and line walls. Returns a bool.
"""
function is_point_in_wall(point::Point{T},wall::LineWall{T}) where T<:Real
  return on_segment(wall.nodes[1],point,wall.nodes[2])
end
function is_point_in_wall(point::Point{T},wall::CircleWall{T}) where T<:Real
  cx = wall.center.x; cy=wall.center.y
  val = abs(wall.radius - sqrt((cx-point.x)^2 + (cy-point.y)^2))
  return isapprox(val, 0; atol=eps(Float64), rtol=0) ? true : false
end

"""
  on_segment(p,q,r)

Checks whether the point `q` is on the line segment from `p` to `r`.
Returns a bool.
"""
function on_segment(p::Point{T},q::Point{T},r::Point{T}) where T<:Real
  if abs(distance(p, q) + distance(q, r) - distance(p, r)) < eps()
    return true
  else
    return false
  end
end

function distance(a::Point{T},b::Point{T}) where T<:Real
  return sqrt((a.x - b.x)^2 + (a.y - b.y)^2)
end

max_stars_repo_name: pseastham/CircGeometry

max_issues_repo_name: pseastham/CircGeometry


text:
function tensordecomp(Λ::AbstractArray{T,3}, p::Int=5;
            scaleout::Bool=true) where T
    𝒯=TensorToolbox.cp_als(Λ,p)
    𝕏=TensorToolbox.full(𝒯)
    A=mean(𝕏[:,:,i] for i=1:size(𝕏,3))
    if scaleout
        A ./=maximum(abs.(A))
    end
    return A
end

max_stars_repo_name: jamesjcai/ScTenifold

max_issues_repo_name: jamesjcai/ScTenifold

max_forks_repo_name: jamesjcai/ScTenifold


text:
"""
Authors: Adrien BLASSIAU Corentin JUVIGNY
TP RODD - 2019-2020
"""

include("../helpers/generate.jl")
include("../helpers/problem.jl")

using .G_ULS
using .Generator

@doc """
Prend un fichier en entrée et exécute un algorithme avec.
""" ->
function main(input_file::String)

	inst = generate(input_file)
	sol = Solution(inst)

	cpu_time = @elapsed(others = run(inst, sol))

	post_process(cpu_time, inst, sol, others)
end

max_stars_repo_name: AdrienBlassiau/RODD

max_issues_repo_name: AdrienBlassiau/RODD

max_forks_repo_name: AdrienBlassiau/RODD


text:
using Temporal
using Base.Test
using Base.Dates

n = 100                                                 # number of observations
k = 4                                                   # number of variables
data = cumsum(randn(n,k), 1)                            # toy random data
dates = collect(today()-Week(n-1):Week(1):today())      # range of Date
times = collect(now():Hour(1):now()+Hour(n-1))          # range of DateTime
fields = ["Field $i" for i=1:k]                         # array of String field names
X = ts(data, dates)                                     # auto-generate field names
Y = ts(data, times, fields)                             # specify field names

# Basic functionality
@test isa(X, TS)
@test isa(Y, TS)
@test size(X) == size(Y) == (n, k)
@test X.values == Y.values == data

# Indexing
i = 1
r = 1:2
a = collect(r)
dsr = dates[1]:Week(1):dates[2]
tsr = times[1]:Hour(1):times[2]
# Single row
@test X[i].values == data[i,:]'
@test Y[i].values == data[i,:]'
@test X[i,i].values == [data[i,i]]
@test Y[i,i].values == [data[i,i]]
@test X[i,r].values == data[i,r]'
@test Y[i,r].values == data[i,r]'
@test X[i,a].values == data[i,r]'
@test Y[i,a].values == data[i,r]'
# Single date/time
@test X[dates[i]].values == data[i,:]'
@test Y[times[i]].values == data[i,:]'
@test X[dates[i],i].values == [data[i,i]]
@test Y[times[i],i].values == [data[i,i]]
@test X[dates[i],r].values == data[i,r]'
@test Y[times[i],r].values == data[i,r]'
@test X[dates[i],a].values == data[i,r]'
@test Y[times[i],a].values == data[i,r]'
# Range of rows
@test X[r].values == data[r,:]
@test Y[r].values == data[r,:]
@test X[r,i].values == data[r,i]
@test Y[r,i].values == data[r,i]
@test X[r,r].values == data[r,r]
@test Y[r,r].values == data[r,r]
@test X[r,a].values == data[r,r]
@test Y[r,a].values == data[r,r]
# Range of dates/times
@test X[dates[r]].values == data[r,:]
@test Y[times[r]].values == data[r,:]
@test X[dates[r],i].values == data[r,i]
@test Y[times[r],i].values == data[r,i]
@test X[dates[r],r].values == data[r,r]
@test Y[times[r],r].values == data[r,r]
@test X[dates[r],a].values == data[r,r]
@test Y[times[r],a].values == data[r,r]
# Array of rows
@test X[a].values == data[r,:]
@test Y[a].values == data[r,:]
@test X[a,i].values == data[r,i]
@test Y[a,i].values == data[r,i]
@test X[a,r].values == data[r,r]
@test Y[a,r].values == data[r,r]
@test X[a,a].values == data[r,r]
@test Y[a,a].values == data[r,r]
# Array of dates/times
@test X[collect(dates[r])].values == data[r,:]
@test Y[collect(times[r])].values == data[r,:]
@test X[collect(dates[r]),i].values == data[r,i]
@test Y[collect(times[r]),i].values == data[r,i]
@test X[collect(dates[r]),r].values == data[r,r]
@test Y[collect(times[r]),r].values == data[r,r]
@test X[collect(dates[r]),a].values == data[r,r]
@test Y[collect(times[r]),a].values == data[r,r]
# Step ranges
@test X[dsr].values == data[r,:]
@test Y[tsr].values == data[r,:]
@test X[dsr,i].values == data[r,i]
@test Y[tsr,i].values == data[r,i]
@test X[dsr,r].values == data[r,r]
@test Y[tsr,r].values == data[r,r]


text:
###############################################################################
#
# Plot stoma on axis
#
###############################################################################
"""
    plot_stoma!(
               ax::PyObject,
               xy::Tuple{Number,Number};
               width::Number = 10,
               height::Number = 10,
               stoma::Number = 0.2,
               angle::Number = 0
    )

Plot a stoma on the axis, given
- `ax` Axis to plot on
- `xy` Center of the stoma
- `width` Width of the stoma
- `height` Height of the stoma
- `stoma` Stomatal pore width ratio
- `angle` Rotation angle of the stoma
"""
function plot_stoma!(
            ax::PyObject,
            xy::Tuple{Number,Number};
            width::Number = 10,
            height::Number = 10,
            stoma::Number = 0.2,
            angle::Number = 0
)
    patch = pyimport("matplotlib.patches");
    outl  = patch.Ellipse(xy=xy,
                          width=width,
                          height=height,
                          color="green",
                          angle=angle,
                          clip_on=false);
    pore  = patch.Ellipse(xy=xy,
                          width=width*stoma,
                          height=height*0.9,
                          color="white",
                          angle=angle,
                          clip_on=false);
    ax.add_artist(outl);
    ax.add_artist(pore);

    return nothing
end


text:
using DataFrames

function sxtevt_df2vec(df::DataFrame)
    ti=df[:TIME]
    ccdframe = df[:CCDFrame]
    x= df[:X]
    y=df[:Y]
    rawx= df[:RAWX]
    rawy= df[:RAWY]
    detx= df[:DETX]
    dety= df[:DETY]
    phas= [df[:PHAS1] df[:PHAS2] df[:PHAS3] df[:PHAS4] df[:PHAS5] df[:PHAS6] df[:PHAS7] df[:PHAS8] df[:PHAS9]]
    stts= df[:STATUS]
    pha= df[:PHA]
    grade= df[:GRADE]
    pixsabove= df[:PixsAbove]
    pi_channel= df[:PI]
    ra= df[:RA]
    decl= df[:Dec]
    phaso= [df[:PHASO1] df[:PHASO2] df[:PHASO3] df[:PHASO4] df[:PHASO5] df[:PHASO6] df[:PHASO7] df[:PHASO8] df[:PHASO9]]
    evt_col_names=["TIME", "CCDFrame", "X", "Y", "RAWX", "RAWY", "DETX", "DETY", "PHAS", "STATUS","PHA", "GRADE","PixsAbove", "PI", "RA","Dec", "PHASO"]
    evt_col_vals=[ti,ccdframe, x,y,rawx,rawy,detx,dety,phas,stts,pha,grade,pixsabove, pi_channel,ra,decl,phaso]
    return evt_col_names, evt_col_vals
    end

    function sxtbadpix_df2vec(df::DataFrame)
        rawx= df[:RAWX]
        rawy= df[:RAWY]
        typ= df[:TYPE]
        yextnt= df[:YEXTENT]
        badflag= df[:BADFLAG]
        badpix_col_names=["RAWX","RAWY","TYPE","YEXTENT","BADFLAG"]
        badpix_col_vals=[rawx,rawy,typ,yextnt,badflag]
        return badpix_col_names, badpix_col_vals
    end

    function gti_df2vec(df::DataFrame)
        gti_start=df[:START]
        gti_stop=df[:STOP]
        return gti_start, gti_stop
    end

max_stars_repo_name: gulabd/SXTMerger

max_issues_repo_name: gulabd/SXTMerger


text:
convert_html(content::AbstractUIHTMLElement) = content.html

function convert_html(content::NTuple{N, E}) where {N, E <: AbstractUIHTMLElement}
    join(map(c -> convert_html(c), content), "")
end

function default_header(title)
    UIHeader("""

<v-app-bar
color="indigo darken-4"
dark>

        <v-toolbar-title>
            $title
        </v-toolbar-title>
        <v-spacer></v-spacer>
</v-app-bar>

    """)
end

function generate_template(title, ui, server_module)
    unrolled = unroll(ui())
    header = extract_uitype(UIHeader, unrolled)
    footer = extract_uitype(UIFooter, unrolled)
    content = extract_uitype(UIElement, unrolled)
    models = extract_uitype(UIModel, unrolled)
    watch = extract_uitype(UIWatch, unrolled)

    dep_tree = dependency_tree(server_module)
    rev_dep = reverse_dependency_tree(server_module)

    if length(header) == 0
        header = default_header(title)
    end

    header = convert_html(header)
    footer = convert_html(footer)
    content = convert_html(content)
"""
<!DOCTYPE html>

  <link href="https://cdn.jsdelivr.net/npm/@mdi/font@4.x/css/materialdesignicons.min.css" rel="stylesheet">
  <meta name="viewport" content="width=device-width, initial-scale=1, maximum-scale=1, user-scalable=no, minimal-ui">
  <link rel="icon" href="favicon.ico" type="image/x-icon" />
  <title>$title</title>

<v-app>
$header
<v-main>
$content
</v-main>
$footer
    <v-snackbar
        v-model="error_snackbar"
        color = "error"
        :timeout="10000">
            {{ matte_error_msg }}
            <v-btn
            text
            @click="error_snackbar = false">
                Close
            </v-btn>
    </v-snackbar>
    <v-overlay :value="matte_notconnected_overlay">

        <v-progress-circular
          indeterminate
          color="primary"
        ></v-progress-circular>

    </v-overlay>
</v-app>

"""
end


text:
using Powersense, GLPK
using JuMP

"""
The attributes to the Powersense and the external solver can be added using the following method with JuMP
"""
#optimizer_with_attributes( Powersense.Optimizer, "external_optimizer" => optimizer_with_attributes((GLPK.Optimizer,"msg_lev"=>GLPK.MSG_OFF)))

model = Model(Powersense.Optimizer);
set_optimizer_attribute(model, "external_optimizer", GLPK.Optimizer)

@variable(model, X);
@variable(model, Y);
@objective(model, Min, X^2 + X);
@NLconstraint(model, X^2 - X == 2);
@NLconstraint(model, X*Y == 1);
@NLconstraint(model, X*Y >= 0);
@constraint(model, X >= -2);

println("________________________________________");
print(model);
println("________________________________________");
JuMP.optimize!(model);

xsol = JuMP.value.(X)
ysol = JuMP.value.(Y)
status = termination_status(model)

println("Xsol = ", xsol);
println("Ysol = ", ysol);

println("Status: ", status);


text:
module Day13

using ..ReTest, ..DocStringExtensions

import ..solve, ..PuzzleInfo, ..TEST_DATA_DIRECTORY

abstract type Day13PuzzleInfo <: PuzzleInfo end

"""
Day 13 part 1 puzzle information.

$(FIELDS)
"""
struct Part1PuzzleInfo <: Day13PuzzleInfo
    """Path to input data file."""
    datapath::String
end

"""
Day 13 part 2 puzzle information.

$(FIELDS)
"""
struct Part2PuzzleInfo <: Day13PuzzleInfo
    """Path to input data file."""
    datapath::String
end

"""
    $(FUNCTIONNAME)(input)

Read fold instructions from input stream `input`.
"""
function readfoldinstructions(input::IO)
    dotindices = [
        CartesianIndex(map(s -> parse(Int, s), split(pairstring, ','))...)
        for pairstring in split(readuntil(input, "\n\n"), '\n')
    ]
    function parsefoldinstruction(foldinstruction)
        instruction, coordinate = split(foldinstruction, '=')
        (instruction[end], parse(Int, coordinate))
    end
    folds = map(parsefoldinstruction, readlines(input))
    dotindices, folds
end

"""
    $(FUNCTIONNAME)(dotindices)

Get bit matrix with true values corresponding to indexes in `dotindices` minus an index
offset `offset`.
"""
function getdotsandoffset(dotindices::Vector{CartesianIndex{2}})
    min_x = minimum(p -> p[1], dotindices)
    max_x = maximum(p -> p[1], dotindices)
    min_y = minimum(p -> p[2], dotindices)
    max_y = maximum(p -> p[2], dotindices)
    dots = falses(
        max_x - min_x + (isodd(max_x - min_x) ? 2 : 1),
        max_y - min_y + (isodd(max_y - min_y) ? 2 : 1)
    )
    offset = CartesianIndex(min_x - 1, min_y - 1)
    for index in dotindices
        dots[index - offset] = true
    end
    dots, offset
end

"""
    $(FUNCTIONNAME)(dots)

Convert bit matrix `dots` to a corresponding Unicode 'block' string.
"""
function dotstostring(dots::BitMatrix)
    join(
        (join(dots[i, j] ? '█' : ' ' for i in 1:size(dots, 1)) for j in 1:size(dots, 2)),
        '\n'
    )
end

function solve(info::Day13PuzzleInfo)
    dotindices, folds = open(readfoldinstructions, info.datapath)
    dots, offset = getdotsandoffset(dotindices)
    for (axis, coordinate) in folds[1:(isa(info, Part1PuzzleInfo) ? 1 : end)]
        if axis == 'x'
            leftdots = dots[1:coordinate-offset[1]-1, :]
            rightdots = dots[end:-1:coordinate-offset[1]+1, :]
            dots = leftdots .| rightdots
        else
            topdots = dots[:, 1:coordinate-offset[2]-1]
            bottomdots = dots[:, end:-1:coordinate-offset[2]+1]
            dots = topdots .| bottomdots
        end
    end
    if isa(info, Part1PuzzleInfo)
        count(dots)
    else
        print(dotstostring(dots))
        dots
    end
end

@testset begin
    @test solve(Part1PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day13.txt"))) == 17
end

end # module


text:
mutable struct ScaffoldServer{T1<:Int64, T2<:Float64, T3<:SparseMatrixCSC{Float64, Int64}, T4<:Vector{Int64}, T5<:Vector{ScaffoldClient}, T6<:Matrix{Float64}} <:AbstractServer
    Xtest::T3                        # testing data
    Ytest::T4                        # testing label
    num_classes::T1                  # number of classes
    num_clients::T1                  # number of clients
    participation_rate::T2           # participation rate
    clients::T5                      # set of clients
    W::T6                           # global models
    C::T6                           # Control variate
    τ::T1                           # number of particiating clients
    lr::T2                          # global learning rate
    selectedIndices::T4             # selected clients
    function ScaffoldServer(Xtest::SparseMatrixCSC{Float64, Int64}, Ytest::Vector{Int64}, config::Dict{String, Real})
        num_classes = config["num_classes"]
        num_clients = config["num_clients"]
        participation_rate = config["participation_rate"]
        lr = config["learning_rate"]
        τ = floor(Int64, num_clients * participation_rate)
        _, d = size(Xtest)
        clients = Vector{ScaffoldClient}(undef, num_clients)
        W = zeros(Float64, d, num_classes)
        C = zeros(Float64, d, num_classes)
        new{Int64, Float64, SparseMatrixCSC{Float64, Int64}, Vector{Int64}, Vector{ScaffoldClient}, Matrix{Float64}}(Xtest, Ytest, num_classes, num_clients, participation_rate, clients, W, C, τ, lr)
    end
end

# Select clients
function select!(
    server::ScaffoldServer
    )
    server.selectedIndices = randperm(server.num_clients)[1:server.τ]
    return nothing
end

# Send global model
function sendModel!(
    server::ScaffoldServer
)
    # Only send model to selected clients
    for idx in server.selectedIndices
        server.clients[idx].W .= copy(server.W)
    end
    return nothing
end

# Send global model to all clients
function sendModelToAllClients!(
    server::ScaffoldServer
)
    # Only send model to selected clients
    for client in server.clients
        client.W .= copy(server.W)
    end
    return nothing
end

# Aggregation
function aggregate!(
    server::ScaffoldServer
)
    # Take the average of the selected clients' model
    d, K = size( server.W )
    ΔW = zeros(Float64, d, K)
    ΔC = zeros(Float64, d, K)
    for i = 1:server.τ
        idx = server.selectedIndices[i]
        ΔW += server.clients[idx].ΔW
        ΔC += server.clients[idx].ΔC
    end
    ΔW ./= server.τ
    ΔC ./= server.num_clients
    # Update
    server.W .+= server.lr*ΔW
    server.C .+= server.lr*ΔC
    return nothing
end

# Calculate objective value
function getObjValue(
    server::ScaffoldServer
)
    objValue = 0.0
    totalNumData = 0
    for i = 1:server.num_clients
        n = size(server.clients[i].Xtrain, 1)
        totalNumData += n
        objValue += n * getObjValue(server.clients[i])
    end
    return objValue / totalNumData
end

########################################################################################################
mutable struct ScaffoldServerNN{T1<:Int64, T2<:Float64, T3<:SparseMatrixCSC{Float64, Int64}, T4<:Zygote.Params, T5<:Vector{ScaffoldClientNN}, T6<:Flux.Chain, T7<:Vector{Int64}, T8<:Flux.OneHotArray} <:AbstractServer
    Xtest::T3                        # testing data
    Ytest::T8                        # testing label
    num_classes::T1                  # number of classes
    num_clients::T1                  # number of clients
    participation_rate::T2           # participation rate
    clients::T5                      # set of clients
    W::T6                           # global models
    C::T4                           # Control variate
    τ::T1                           # number of particiating clients
    lr::T2                          # global learning rate
    selectedIndices::T7             # selected clients
    function ScaffoldServerNN(Xtest::SparseMatrixCSC{Float64, Int64}, Ytest::Vector{Int64}, dim::Int64, config::Dict{String, Real})
        num_classes = config["num_classes"]
        Ytest = Flux.onehotbatch(Ytest, 1:num_classes)
        num_clients = config["num_clients"]
        participation_rate = config["participation_rate"]
        lr = config["learning_rate"]
        τ = floor(Int64, num_clients * participation_rate)
        # Set up clients and initialize models and control variates.
        clients = Vector{ScaffoldClient}(undef, num_clients)
        W = Chain( Dense(780, dim, relu, bias=false), Dense(dim, 10, bias=false), NNlib.softmax)
        C = deepcopy( params(W) )
        for j = 1:length(C)
            fill!( C[j], 0.0 )
        end
        new{Int64, Float64, SparseMatrixCSC{Float64, Int64}, Zygote.Params, Vector{ScaffoldClientNN}, Flux.Chain, Vector{Int64}, Flux.OneHotArray}(Xtest, Ytest, num_classes, num_clients, participation_rate, clients, W, C, τ, lr)
    end
end

# Select clients
function select!(
    server::ScaffoldServerNN
    )
    server.selectedIndices = randperm(server.num_clients)[1:server.τ]
    return nothing
end

# Send global model
function sendModel!(
    server::ScaffoldServerNN
)
    # Only send model to selected clients
    for i in server.selectedIndices
        c = server.clients[i]
        for j = 1:length(params(server.W))
            params(c.W)[j] .= params(server.W)[j]
            c.globalC[j] .= server.C[j]
        end
    end
    return nothing
end

# Send global model to all clients
function sendModelToAllClients!(
    server::ScaffoldServerNN
)
    for c in server.clients
        for j = 1:length(params(server.W))
            params(c.W)[j] .= params(server.W)[j]
            c.globalC[j] .= server.C[j]
        end
    end
    return nothing
end

# Aggregation
function aggregate!(
    server::ScaffoldServerNN
)
    # Take the average of the selected clients' model
    l = length(params(server.W))
    for i = 1:server.τ
        idx = server.selectedIndices[i]
        c = server.clients[idx]
        for j = 1:l
            params(server.W)[j] .+= (server.lr/server.τ) * c.ΔW[j]
            server.C[j] .+= (server.lr/server.num_clients) * c.ΔC[j]
        end
    end
    return nothing
end

# Calculate objective value
function getObjValue(
    server::ScaffoldServerNN
)
    objValue = 0.0
    totalNumData = 0
    for i = 1:server.num_clients
        c = server.clients[i]
        n = size(c.Xtrain, 1)
        totalNumData += n
        objValue += n * obj( c.XtrainT, c.Ytrain, server.W, c.λ)
    end
    return objValue / totalNumData
end


text:
############################################################
# joNumber methods #########################################
############################################################

## constructors
"""
joNumber outer constructor

    joNumber(num)

Create joNumber with types matching given number

"""
joNumber(num::NT) where {NT<:Number} = joNumber{NT,NT}(num,num)
"""
joNumber outer constructor

    joNumber(num,A::joAbstractLinearOperator{DDT,RDT})

Create joNumber with types matching the given operator.

"""
joNumber(num::NT,A::joAbstractLinearOperator{DDT,RDT}) where {NT<:Number,DDT,RDT} =
    joNumber(jo_convert(DDT,num),jo_convert(RDT,num)) # where {DDT,RDT}

## other functions
-(n::joNumber{DDT,RDT}) where {DDT,RDT} = joNumber(-n.ddt,-n.rdt) # where {DDT,RDT}
inv(n::joNumber{DDT,RDT}) where {DDT,RDT} = joNumber(inv(n.ddt),inv(n.rdt)) # where {DDT,RDT}


text:
import Pkg; Pkg.activate(joinpath(@__DIR__, "../../../"))

using Turing
using Distributions
using CytofDensityEstimation
using StatsPlots
import LinearAlgebra
const CDE = CytofDensityEstimation
include("scenarios.jl")

@model M1(yC_finite, yT_finite, K) = begin
  etaC ~ Dirichlet(K, 1/K)
  etaT ~ Dirichlet(K, 1/K)

  mu ~ filldist(Normal(0, 3), K)
  sigma ~ filldist(LogNormal(0, 1), K)
  nu ~ filldist(LogNormal(0, 1), K)
  phi ~ filldist(Normal(0, 6), K)

  yC_finite .~ MixtureModel(CDE.Util.SkewT.(mu, sigma, nu, phi), etaC)
  yT_finite .~ MixtureModel(CDE.Util.SkewT.(mu, sigma, nu, phi), etaT)
end
numparams(K, beta) = beta > 0 ? 6K : 5K
m1params = [:etaC, :etaT, :mu, :sigma, :nu, :phi]
get_last_state(sym, chain) = group(chain[end], sym).value.data
get_last_state(chain) = [get_last_state(sym, chain) for sym in m1params]
get_last_state(chain1)
getparam(sym, chain) = group(chain, sym).value.data[:, :, 1]

simdata = scenarios(1)
yC_finite = filter(isfinite, simdata[:yC])
yT_finite = filter(isfinite, simdata[:yT])
K = 6
m1 = M1(yC_finite, yT_finite, K)

function getinit(; s=[.1, .01, .01], nmcmc=fill(100, 3))
  init = nothing
  chain = nothing
  for i in 1:length(s)
    mh = MH(LinearAlgebra.I(numparams(K, 1)) * s[i])
    if i == 1
      chain = sample(m1, mh, nmcmc[i], save_state=true)
    else
      chain = sample(m1, mh, nmcmc[i], save_state=true, init_theta=init)
    end
    init = get_last_state(chain)
  end
  index_lp = findall(names(chain) .== :lp)[1]

  mat = hcat([getparam(sym, chain) for sym in m1params]...)
  return init, cov(unique(mat, dims=1))
end

init, _ = getinit(s=[.1, .01, .001], nmcmc=[100, 1000, 1000])
mh = MH(LinearAlgebra.I(numparams(K, 1)) * 0.001)
chain1 = sample(m1, mh, 3000, init_theta=init);
plot(vec(get(chain1, :lp)[1].data), label=nothing)

plot(group(chain1, :mu).value.data[:, :, 1], label=nothing)
plot(group(chain1, :etaC).value.data[:, :, 1], label=nothing)
plot(group(chain1, :etaT).value.data[:, :, 1], label=nothing)

param = Dict((sym, getparam(sym, chain1)) for sym in m1params)

function postpred(param; ygrid)
  B = size(param[:mu], 1)
  comps(b) = CDE.Util.SkewT.(param[:mu][b, :], param[:sigma][b, :],
                             param[:nu][b, :], param[:phi][b, :])
  etaC(b) = param[:etaC][b, :]
  etaT(b) = param[:etaT][b, :]

  dC = hcat([pdf.(MixtureModel(comps(b), etaC(b)), ygrid) for b in 1:B]...)
  dT = hcat([pdf.(MixtureModel(comps(b), etaT(b)), ygrid) for b in 1:B]...)

  return (dC=dC, dT=dT)
end

ygrid = range(-6, 6, length=100)
pp = postpred(param, ygrid=ygrid)

pplower = vec(CDE.MCMC.quantiles(pp[:dC], .025, dims=2))
ppupper = vec(CDE.MCMC.quantiles(pp[:dC], .975, dims=2))
plot(ygrid, pplower, fillrange=ppupper, alpha=.3, color=:blue, label="C", legend=:topleft)
plot!(ygrid, pdf.(simdata[:mmC], ygrid), ls=:dot, color=:black, label=nothing)

pplower = vec(CDE.MCMC.quantiles(pp[:dT], .025, dims=2))
ppupper = vec(CDE.MCMC.quantiles(pp[:dT], .975, dims=2))
plot!(ygrid, pplower, fillrange=ppupper, alpha=.3, color=:red, label="T")
plot!(ygrid, pdf.(simdata[:mmT], ygrid), ls=:dot, color=:black, label=nothing)

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
# This is working but should live in a @recipe instead.
# I am waiting for a more stable version of the recipes in Makie, or at least one
# that doesn't need the full Makie thing as a dependency.
# Needs to do: `using GLMakie, Color, ColorSchemes`
function plot_mtg_Makie(mtg; color_var = :YY)

    df_coordinates = mtg_coordinates_df(mtg)
    colouring_var = df_coordinates[:,color_var]
    color = colouring_var ./ maximum(colouring_var)

    fig, ax, p = scatter(df_coordinates.XX, df_coordinates.YY, df_coordinates.ZZ, color = color, colormap = :viridis)
    # ?Note: could use meshscatter! instead here

    for i in 2:size(df_coordinates, 1)
        lines!(
            [df_coordinates.XX_from[i], df_coordinates.XX[i]],
            [df_coordinates.YY_from[i], df_coordinates.YY[i]],
            [df_coordinates.ZZ_from[i], df_coordinates.ZZ[i]],
            color = get(ColorSchemes.viridis, [color[i - 1],color[i]])
        )
    end
    # DataInspector(fig)
    hidedecorations!(ax);
    hidespines!(ax)
    # ax.aspect = DataAspect()
    # p.nlabels_offset[] = Point2f(0.2, - 0.2)
    # autolimits!(ax)
    fig
end

plot_mtg_Makie(mtg)


text:
#
# Replace default random number generator by other that aleviates the
# mutiple-threading problems
#

import Random
using StableRNGs
using Future:randjump

function init_random(options)
  if options.seed > 0
    if options.StableRNG == true
      RNG = StableRNGs.StableRNG(options.seed)
    else
      RNG = [randjump(Random.MersenneTwister(options.seed),big(10)^20)]
      foreach(_ -> push!(RNG, randjump(last(RNG),big(10)^20)), 2:Threads.nthreads())
    end
  else
    if options.StableRNG == true
      seed = abs(rand(Int))
      RNG = StableRNGs.StableRNG(seed)
    else
      RNG = [randjump(Random.MersenneTwister(),big(10)^20)]
      foreach(_ -> push!(RNG, randjump(last(RNG),big(10)^20)), 2:Threads.nthreads())
    end
  end
  return RNG
end

random(RNG::StableRNGs.LehmerRNG) = rand(RNG)
random(RNG::StableRNGs.LehmerRNG, arg) = rand(RNG,arg)

random(RNG::Array{T}) where T = rand(RNG[Threads.threadid()])
random(RNG::Array{T}, arg) where T = rand(RNG[Threads.threadid()],arg)


text:
export getRotate3D,rotateAndMove3D,rotateAndMove3DTranspose,smooth3D,rotateSensitivity3D

function loc2cs3D(loc1::Union{Int64,Array{Int64}},loc2::Union{Int64,Array{Int64}},loc3::Union{Int64,Array{Int64}},n::Array{Int64,1})
@inbounds cs = loc1 .+ (loc2.-1)*n[1] .+ (loc3.-1)*n[1]*n[2];
return cs;
end

function getRotate3D(theta_rad::Float64,phi_rad::Float64)
## theta rotates on the x-y plain (or x1-x2)
## phi rotates on the x-z plain (or x1-x3).
ct = cos(theta_rad); st = sin(theta_rad);
cp = cos(phi_rad); sp = sin(phi_rad);
Rxy = [ct -st 0.0; st ct 0.0; 0.0 0.0 1.0];
Rxz = [cp 0.0 sp; 0.0 1.0 0.0; -sp 0.0 cp];
R   = Rxz*Rxy;
return R;
end

function rotateAndMove3D(u::Array{Float64,3},theta_phi_rad::Vector{Float64},b::Vector{Float64},doTranspose=false,v::Array = copy(u),XT::Array = zeros(0),XTT = zeros(0))
###
# b here is in pixels - not in Mesh units (no mesh here...)
###
n = collect(size(u));
v[:] .= 0.0;
R = getRotate3D(theta_phi_rad[1],theta_phi_rad[2]);
invR = inv(R);
midphalf = n./2.0 .+ 0.5;

midphalf_plus_b = midphalf .+ b;

maxU = maximum(u);
thresh = 1e-4*maxU;

xt = zeros(3);
x = zeros(3);

if length(XT)==0
	XT = zeros(3,prod(n));
	XTT = zeros(3,prod(n));
	# println("Building XT");
else
	XT[:].=0.0;
	XTT[:].=0.0;
end
for i3 = 1:size(u,3)
	for i2 = 1:size(u,2)
		for i1 = 1:size(u,1)
			@inbounds xt[1] = i1; xt[2] = i2; xt[3] = i3;
			# x = xt - midphalf_plus_b;
			BLAS.axpy!(-1.0,midphalf_plus_b,xt); # xt = xt - midphalf_plus_b;
			ii = loc2cs3D(i1,i2,i3,n);
			@inbounds XT[:,ii] = xt;
		end
	end
end
# gemm!(tA, tB, alpha, A, B, beta, C)
XTT[1,:] .= midphalf[1];
XTT[2,:] .= midphalf[2];
XTT[3,:] .= midphalf[3];
BLAS.gemm!('N','N',1.0,invR,XT,1.0,XTT);
# I think that here we assume that the domain is square to simplify.
II = findall(vec(sum((XTT .> 0.5) .& (XTT .< n[1]),dims=1).>=3));

LOC = round.(Int64,XTT[:,II]);
LOC = loc2cs3D(LOC[1,:],LOC[2,:],LOC[3,:],n);

if !doTranspose
	v[II] = u[LOC];
else
	# The following two options are actually different in terms of repeatitions:
	# v[LOC] .+= u[II];
	for ii=1:length(LOC)
		v[LOC[ii]] += u[II[ii]];
	end
end

# the code below is an altenative to the code above.
# iii=1;
# for i3 = 1:size(u,3)
	# for i2 = 1:size(u,2)
		# for i1 = 1:size(u,1)
			# @inbounds xt[1] = i1; xt[2] = i2; xt[3] = i3;
			# x = xt - midphalf_plus_b;
			# BLAS.axpy!(-1.0,midphalf_plus_b,xt); # xt = xt - midphalf_plus_b;
			# BLAS.gemv!('N',1.0,invR,xt,0.0,x); # x = 0.0*x + invR*xt;
			# BLAS.axpy!(1.0,midphalf,x); # x = xt + midphalf;
			# @inbounds x[:] = XTT[:,iii];
			# iii+=1;

			# @inbounds if x[1] < 0.5 || (x[1] > n[1]) || (x[2] < 0.5) || (x[2] > n[2])|| x[3] < 1 || x[3] > n[3]
				# continue;
			# end
			# ii = round.(Int64,x);
			# @inbounds v[i1,i2,i3] = u[ii[1],ii[2],ii[3]];
		# end
	# end
# end
return v,XT,XTT;
end

# function rotateAndMove3DTranspose(v::Array,theta_phi_rad::Array{Float64},b::Array{Float64})
# ###
# # b here is in pixels - not in Mesh units (no mesh here...)
# ###
# n = collect(size(v));
# u = zeros(eltype(v),size(v));
# R = getRotate3D(theta_phi_rad[1],theta_phi_rad[2]);
# invR = inv(R);
# mid = n/2.0;
# maxU = maximum(v);
# thresh = 1e-4*maxU;
# for i3 = 1:size(v,3)
	# for i2 = 1:size(v,2)
		# for i1 = 1:size(v,1)
			# @inbounds vi = v[i1,i2,i3];
			# if abs(vi) > thresh
				# x = [i1;i2;i3] .- (mid .+ 0.5 ) .- b; ### here we actually do the inverse operation in a forward way
				# x = invR*x;
				# x .+= mid;
				# @inbounds if x[1] <=  0.5 || (x[1] >= n[1]-0.5) || (x[2] <= 0.5) || (x[2] >= n[2]-0.5)|| x[3] <= 0.5 || x[3] >= n[3]-0.5
					# warn("Rotation3D: rotation went out of bounds. Please increase zero padding.");
					# continue;
				# end
				# x .+= 0.5;
				# xloc = convert(Array{Int64},round.(x));
				# @inbounds u[xloc[1],xloc[2],xloc[3]] += vi;
			# end
		# end
	# end
# end
# return u;
# end

# function smooth3D(v::Array,numFilters = 1)
# N = collect(-1:1);
# E3 = ones(3,3,3);
# E2 = [1.0 ; 2.0 ; 1.0]*[1.0 ; 2.0 ; 1.0]';
# E3[:,:,1] = E2;
# E3[:,:,2] = 2*E2;
# E3[:,:,3] = E2;
# E = E ./ sum(E[:]);

# for k = 1:2
	# for i3 = 2:size(u,3)-1
		# for i2 = 2:size(u,2)-1
			# for i1 = 2:size(u,1)-1
				# v[i1, i2, i3] = sum(E.*v[i1 + N, i2 + N, i3 + N]);
			# end
		# end
	# end
# end
# return v;
# end


text:
# SPDX-License-Identifier: X11
# 2020-11-21 (Timed)
# Quadruple (400pt)

function main()
  N,K = parse.(Int, (split ∘ readline)())
  K = abs(K)

  sum = 0
  for i ∈ K + 2:2N
    sum += pairs(N, i) * pairs(N, i - K)
  end

  println(sum)
end

pairs(n::Integer, k::Integer) = min(k - 1, 2n + 1 - k)

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
module Polynomials

export Polynomial, deg, X, Y

abstract Poly

immutable Polynomial{T <: Union(Number, Poly)} <: Poly
  coeffs :: Vector{T}
end

Base.zero{T}(::Polynomial{T}) = zero(Polynomial{T})
Base.zero{T}(::Type{Polynomial{T}}) = Polynomial{T}([])

Base.one{T}(::Polynomial{T}) = one(Polynomial{T})
Base.one{T}(::Type{Polynomial{T}}) = Polynomial{T}([one(T)])

Base.copy{T}(p::Polynomial{T}) = Polynomial(copy(p.coeffs))

Base.convert{T}(::Type{Polynomial}, x::T) = convert(Polynomial{T}, x)
Base.convert{T}(::Type{Polynomial{T}}, x) = Polynomial{T}([convert(T,x)])
Base.convert{T}(::Type{Polynomial{T}}, p::Polynomial) = Polynomial{T}([convert(T,a) for a in p.coeffs])

Base.promote_rule{T,R}(::Type{Polynomial{T}}, ::Type{Polynomial{R}}) = Polynomial{promote_type(T,R)}

Base.promote_rule{T,R}(::Type{T}, ::Type{Polynomial{R}}) = Polynomial{promote_type(T,R)}

function +{T}(p::Polynomial{T}, q::Polynomial{T})
  n = min(length(p.coeffs), length(q.coeffs))
  Polynomial{T}([p.coeffs[1:n] .+ q.coeffs[1:n], p.coeffs[n+1:end], q.coeffs[n+1:end]])
end

function *{T}(p::Polynomial{T}, q::Polynomial{T})
  if isempty(p.coeffs)
    p
  else
    Polynomial(T[p.coeffs[1] * x for x in q.coeffs]) + Polynomial([zero(T), (Polynomial{T}(p.coeffs[2:end]) * q).coeffs])
  end
end

-(p::Polynomial) = Polynomial(map(-, p.coeffs))
-(p::Polynomial, q::Polynomial) = +(p,-q)

X = Polynomial([0,1])
Y = Polynomial([X])

+(x::Union(Number,Poly), y::Union(Number,Poly)) = +(promote(x,y)...)
-(x::Union(Number,Poly), y::Union(Number,Poly)) = -(promote(x,y)...)
*(x::Union(Number,Poly), y::Union(Number,Poly)) = *(promote(x,y)...)

iszero(n::Number) = n == zero(n)

iszero(p::Polynomial) = deg(p) == -1

==(p::Polynomial, q::Polynomial) = iszero(p-q)

==(x::Union(Number,Poly), y::Union(Number,Poly)) = ==(promote(x,y)...)

function deg(p::Polynomial)
  d = length(p.coeffs)
  while d > 0 && iszero(p.coeffs[d])
    d -= 1
  end
  d - 1
end

euclid(a::Integer, b::Integer) = divrem(a,b)
euclid(a::Real, b::Real) = (a/b,0)
euclid(a::Complex, b::Complex) = (a/b,0)

function euclid{T}(p::Polynomial{T}, q::Polynomial{T})
  if iszero(q)
    error("Divide by zero")
  elseif deg(q) > deg(p)
    (zero(p), p)
  else
    (d,r) = euclid(p.coeffs[1+deg(p)], q.coeffs[1+deg(q)])
    (d,r) = map(x -> Polynomial([x]), (d,r))
    d = d * X^(deg(p)-deg(q))
    r = r * X^deg(p)
    (dd,rr) = euclid(p-(d*q+r), q)
    (d+dd,r+rr)
  end
end

Base.divrem(x::Union(Number,Poly), y::Union(Number,Poly)) = euclid(promote(x,y)...)
Base.divrem(x::Number, y::Number) = (div(x,y), rem(x,y))

Base.div(p::Polynomial, q::Polynomial) = divrem(p,q)[1]
Base.rem(p::Polynomial, q::Polynomial) = divrem(p,q)[2]

end


text:
#############################################################################
# card.jl
# Handles functions of the form f(A) = g(length(A)), where g is a concave function
#############################################################################

export card, compose, CardBasedAtom
export in
export sign, monotonicity, modularity, evaluate

mutable struct CardBasedAtom{T} <: SubmodFunc
  head::Symbol
  id_hash::UInt64
  children::Tuple{AbstractExpr, SetVariable}
  size::Tuple{Int, Int}
  setvariables::Array{SetVariable}
  func::T
end

function CardBasedAtom(F::Function, S::CombiSet)
  z = Variable(1)
  children = (F(z), S)                           # create an expression
  if F(0) != 0
    error("A function has to return 0 at 0 in order to derive a cardinality-based function.")
  end
  setvariables = get_sv(S)
  return(CardBasedAtom(:card, hash(children), children, (1, 1), setvariables, F))
end

function CardBasedAtom(w::Array{Float64}, S::CombiSet)
  children = (Constant(w), S)                           # create an expression
  setvariables = get_sv(S)
  if length(setvariables) != 1
    error("The expression should have exactly one set variable.")
  elseif size(w) != (setvariables[1].cardinality, )
    error("The sizes of the set variable and the array mismatch.")
  end
  w₁ = copy(w)
  for i = 2:length(w)
    w₁[i] = w[i] - w[i-1]
  end
  return(CardBasedAtom(:card, hash(children), children, (1, 1), setvariables, w₁))
end

card(f::Function, S::CombiSet) = CardBasedAtom(f, S)

card(w::Array{Float64}, S::CombiSet) = CardBasedAtom(w, S)

card(w::Array{Int}, S::CombiSet) = CardBasedAtom(Float64.(w), S)

card(S::CombiSet) = card(x -> x, S)

function sign(F::CardBasedAtom)
  return sign(F.children[1])
end

function sign(F::CardBasedAtom{Array{Float64, 1}})
  if all(x -> x>=0, F.func)
    return Positive()
  elseif all(x -> x<=0, F.func)
    return Negative()
  else
    return NoSign()
  end
end

function monotonicity(F::CardBasedAtom)
  return (monotonicity(F.children[1]), )
end

function monotonicity(F::CardBasedAtom{Array{Float64, 1}})
  if all(x -> x>=0, F.func)
    return Nondecreasing()
  elseif all(x -> x<=0, F.func)
    return Nonincreasing()
  elseif all(x -> x==0, F.func)
    return ConstMonotonicity()
  else
    return NoMonotonicity()
  end
end

function modularity(F::CardBasedAtom)
  if vexity(F.children[1]) == ConcaveVexity()
    return SubModularity()
  elseif vexity(F.children[1]) == ConvexVexity()
    return SuperModularity()
  elseif vexity(F.children[1]) == AffineVexity()
    return Modularity()
  elseif vexity(F.children[1]) == ConstVexity()
    return ConstModularity()
  else
    return NotDetermined()
  end
end

function modularity(F::CardBasedAtom{Array{Float64, 1}})
  w₁ = sort(F.func, rev=true)
  if w₁ == F.func
    return SubModularity()
  else
    w₂ = sort(F.func)
    if w₂ == F.func
      return SuperModularity()
    elseif all(x -> x==0, F.func)
      return ConstModularity()
    else
      return NotDetermined()
    end
  end
end

function evaluate(F::CardBasedAtom)
  var = F.setvariables[1]
  var1 = get_v(F.children[1])
  var1[1].value = length(Set(get_elements(var)))
  evaluate(F.children[1])
end

function evaluate(F::CardBasedAtom{Array{Float64, 1}})
  return sum(F.func[1:length(get_elements(F.children[2]))])
end

function in(w::AbstractArray, p::SubmodPoly{CardBasedAtom})
  n = length(w)
  @assert length(p.V) == n
  checker = true
  S = sort(p.V)
  ordering = sortperm(w, rev = true)
  h = 0
  ind = []
  y = zeros(n + 1)
  for i = 1:length(w)
    y[i + 1] = evaluate(p.F, S[push!(ind, ordering[i])])[1]
    h += - y[i] - w[i] + y[i + 1]
    if h < 0
      checker = false
      break
    end
  end
  return checker
end

function compose(f::Function, g::CardBasedAtom)
  return card(f, g.setvariables[1])
end


text:
using Distributed

# launch worker processes
num_cores = parse(Int, ENV["SLURM_CPUS_PER_TASK"])
addprocs(num_cores)

println("Number of cores: ", nprocs())
println("Number of workers: ", nworkers())

# each worker gets its id, process id and hostname
for i in workers()
    id, pid, host = fetch(@spawnat i (myid(), getpid(), gethostname()))
    println(id, " " , pid, " ", host)
end

# remove the workers
for i in workers()
    rmprocs(i)
end


text:
@doc raw"""
    SymmetricPositiveDefinite{N} <: AbstractDecoratorManifold{𝔽}

The manifold of symmetric positive definite matrices, i.e.

````math
\mathcal P(n) =
\bigl\{
p ∈ ℝ^{n × n}\ \big|\ a^\mathrm{T}pa > 0 \text{ for all } a ∈ ℝ^{n}\backslash\{0\}
\bigr\}
````

The tangent space at ``T_p\mathcal P(n)`` reads

```math
    T_p\mathcal P(n) =
    \bigl\{
        X \in \mathbb R^{n×n} \big|\ X=X^\mathrm{T}
    \bigr\},
```
i.e. the set of symmetric matrices,

# Constructor

    SymmetricPositiveDefinite(n)

generates the manifold $\mathcal P(n) \subset ℝ^{n × n}$
"""
struct SymmetricPositiveDefinite{N} <: AbstractDecoratorManifold{ℝ} end

SymmetricPositiveDefinite(n::Int) = SymmetricPositiveDefinite{n}()

function active_traits(f, ::SymmetricPositiveDefinite, args...)
    return merge_traits(IsEmbeddedManifold(), IsDefaultMetric(LinearAffineMetric()))
end

@doc raw"""
    check_point(M::SymmetricPositiveDefinite, p; kwargs...)

checks, whether `p` is a valid point on the [`SymmetricPositiveDefinite`](@ref) `M`, i.e. is a matrix
of size `(N,N)`, symmetric and positive definite.
The tolerance for the second to last test can be set using the `kwargs...`.
"""
function check_point(M::SymmetricPositiveDefinite{N}, p; kwargs...) where {N}
    if !isapprox(norm(p - transpose(p)), 0.0; kwargs...)
        return DomainError(
            norm(p - transpose(p)),
            "The point $(p) does not lie on $(M) since its not a symmetric matrix:",
        )
    end
    if !all(eigvals(p) .> 0)
        return DomainError(
            eigvals(p),
            "The point $p does not lie on $(M) since its not a positive definite matrix.",
        )
    end
    return nothing
end

"""
    check_vector(M::SymmetricPositiveDefinite, p, X; kwargs... )

Check whether `X` is a tangent vector to `p` on the [`SymmetricPositiveDefinite`](@ref) `M`,
i.e. atfer [`check_point`](@ref)`(M,p)`, `X` has to be of same dimension as `p`
and a symmetric matrix, i.e. this stores tangent vetors as elements of the corresponding
Lie group.
The tolerance for the last test can be set using the `kwargs...`.
"""
function check_vector(M::SymmetricPositiveDefinite{N}, p, X; kwargs...) where {N}
    if !isapprox(norm(X - transpose(X)), 0.0; kwargs...)
        return DomainError(
            X,
            "The vector $(X) is not a tangent to a point on $(M) (represented as an element of the Lie algebra) since its not symmetric.",
        )
    end
    return nothing
end

function get_embedding(M::SymmetricPositiveDefinite)
    return Euclidean(representation_size(M)...; field=ℝ)
end

embed(::SymmetricPositiveDefinite, p) = p
embed(::SymmetricPositiveDefinite, p, X) = X

@doc raw"""
    injectivity_radius(M::SymmetricPositiveDefinite[, p])
    injectivity_radius(M::MetricManifold{SymmetricPositiveDefinite,LinearAffineMetric}[, p])
    injectivity_radius(M::MetricManifold{SymmetricPositiveDefinite,LogCholeskyMetric}[, p])

Return the injectivity radius of the [`SymmetricPositiveDefinite`](@ref).
Since `M` is a Hadamard manifold with respect to the [`LinearAffineMetric`](@ref) and the
[`LogCholeskyMetric`](@ref), the injectivity radius is globally $∞$.
"""
injectivity_radius(::SymmetricPositiveDefinite) = Inf
injectivity_radius(::SymmetricPositiveDefinite, p) = Inf
injectivity_radius(::SymmetricPositiveDefinite, ::AbstractRetractionMethod) = Inf
injectivity_radius(::SymmetricPositiveDefinite, p, ::AbstractRetractionMethod) = Inf

@doc raw"""
    manifold_dimension(M::SymmetricPositiveDefinite)

returns the dimension of
[`SymmetricPositiveDefinite`](@ref) `M`$=\mathcal P(n), n ∈ ℕ$, i.e.
````math
\dim \mathcal P(n) = \frac{n(n+1)}{2}.
````
"""
@generated function manifold_dimension(::SymmetricPositiveDefinite{N}) where {N}
    return div(N * (N + 1), 2)
end

"""
    mean(
        M::SymmetricPositiveDefinite,
        x::AbstractVector,
        [w::AbstractWeights,]
        method = GeodesicInterpolation();
        kwargs...,
    )

Compute the Riemannian [`mean`](@ref mean(M::AbstractManifold, args...)) of `x` using
[`GeodesicInterpolation`](@ref).
"""
mean(::SymmetricPositiveDefinite, ::Any)

function default_estimation_method(::SymmetricPositiveDefinite, ::typeof(mean))
    return GeodesicInterpolation()
end

@doc raw"""
    project(M::SymmetricPositiveDefinite, p, X)

project a matrix from the embedding onto the tangent space $T_p\mathcal P(n)$ of the
[`SymmetricPositiveDefinite`](@ref) matrices, i.e. the set of symmetric matrices.
"""
project(::SymmetricPositiveDefinite, p, X)

project!(::SymmetricPositiveDefinite, Y, p, X) = (Y .= Symmetric((X + X') / 2))

@doc raw"""
    rand(M::SymmetricPositiveDefinite; σ::Real=1)

Generate a random symmetric positive definite matrix on the
`SymmetricPositiveDefinite` manifold `M`.
"""
rand(M::SymmetricPositiveDefinite; σ::Real=1)

function Random.rand!(
    M::SymmetricPositiveDefinite{N},
    pX;
    vector_at=nothing,
    σ::Real=one(eltype(pX)) / (vector_at === nothing ? 1 : norm(vector_at)),
    tangent_distr=:Gaussian,
) where {N}
    if vector_at === nothing
        D = Diagonal(1 .+ rand(N)) # random diagonal matrix
        s = qr(σ * randn(N, N)) # random q
        pX .= Symmetric(s.Q * D * transpose(s.Q))
    elseif tangent_distr === :Gaussian
        # generate ONB in TxM
        I = one(vector_at)
        B = get_basis(M, vector_at, DiagonalizingOrthonormalBasis(I))
        Ξ = get_vectors(M, vector_at, B)
        Ξx =
            vector_transport_to.(
                Ref(M),
                Ref(I),
                Ξ,
                Ref(vector_at),
                Ref(ParallelTransport()),
            )
        pX .= sum(σ * randn(length(Ξx)) .* Ξx)
    elseif tangent_distr === :Rician
        C = cholesky(Hermitian(vector_at))
        R = C.L + sqrt(σ) * triu(randn(size(vector_at, 1), size(vector_at, 2)), 0)
        pX .= R * R'
    end
    return pX
end
function Random.rand!(
    rng::AbstractRNG,
    M::SymmetricPositiveDefinite{N},
    pX;
    vector_at=nothing,
    σ::Real=one(eltype(pX)) / (vector_at === nothing ? 1 : norm(vector_at)),
    tangent_distr=:Gaussian,
) where {N}
    if vector_at === nothing
        D = Diagonal(1 .+ rand(rng, N)) # random diagonal matrix
        s = qr(σ * randn(rng, N, N)) # random q
        pX .= Symmetric(s.Q * D * transpose(s.Q))
    elseif tangent_distr === :Gaussian
        # generate ONB in TxM
        I = one(vector_at)
        B = get_basis(M, vector_at, DiagonalizingOrthonormalBasis(I))
        Ξ = get_vectors(M, vector_at, B)
        Ξx =
            vector_transport_to.(
                Ref(M),
                Ref(I),
                Ξ,
                Ref(vector_at),
                Ref(ParallelTransport()),
            )
        pX .= sum(σ * randn(rng, length(Ξx)) .* Ξx)
    elseif tangent_distr === :Rician
        C = cholesky(Hermitian(vector_at))
        R = C.L + sqrt(σ) * triu(randn(rng, size(vector_at, 1), size(vector_at, 2)), 0)
        pX .= R * R'
    end
    return pX
end

@doc raw"""
    representation_size(M::SymmetricPositiveDefinite)

Return the size of an array representing an element on the
[`SymmetricPositiveDefinite`](@ref) manifold `M`, i.e. $n × n$, the size of such a
symmetric positive definite matrix on $\mathcal M = \mathcal P(n)$.
"""
@generated representation_size(::SymmetricPositiveDefinite{N}) where {N} = (N, N)

function Base.show(io::IO, ::SymmetricPositiveDefinite{N}) where {N}
    return print(io, "SymmetricPositiveDefinite($(N))")
end

@doc raw"""
    zero_vector(M::SymmetricPositiveDefinite,x)

returns the zero tangent vector in the tangent space of the symmetric positive
definite matrix `x` on the [`SymmetricPositiveDefinite`](@ref) manifold `M`.
"""
zero_vector(::SymmetricPositiveDefinite, ::Any)

zero_vector!(::SymmetricPositiveDefinite{N}, v, ::Any) where {N} = fill!(v, 0)


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

struct VirtualNetworksApi <: SwaggerApi
    client::Swagger.Client
end

"""

Creates or updates a virtual network in the specified resource group.
Param: resourceGroupName::String (required)
Param: virtualNetworkName::String (required)
Param: parameters::VirtualNetwork (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: VirtualNetwork
"""
function _swaggerinternal_virtualNetworksCreateOrUpdate(_api::VirtualNetworksApi, resourceGroupName::String, virtualNetworkName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", VirtualNetwork, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/virtualNetworks/{virtualNetworkName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkName", virtualNetworkName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function virtualNetworksCreateOrUpdate(_api::VirtualNetworksApi, resourceGroupName::String, virtualNetworkName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworksCreateOrUpdate(_api, resourceGroupName, virtualNetworkName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function virtualNetworksCreateOrUpdate(_api::VirtualNetworksApi, response_stream::Channel, resourceGroupName::String, virtualNetworkName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworksCreateOrUpdate(_api, resourceGroupName, virtualNetworkName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Deletes the specified virtual network.
Param: resourceGroupName::String (required)
Param: virtualNetworkName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function _swaggerinternal_virtualNetworksDelete(_api::VirtualNetworksApi, resourceGroupName::String, virtualNetworkName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/virtualNetworks/{virtualNetworkName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkName", virtualNetworkName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function virtualNetworksDelete(_api::VirtualNetworksApi, resourceGroupName::String, virtualNetworkName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworksDelete(_api, resourceGroupName, virtualNetworkName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function virtualNetworksDelete(_api::VirtualNetworksApi, response_stream::Channel, resourceGroupName::String, virtualNetworkName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworksDelete(_api, resourceGroupName, virtualNetworkName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets the specified virtual network by resource group.
Param: resourceGroupName::String (required)
Param: virtualNetworkName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: expand::String
Return: VirtualNetwork
"""
function _swaggerinternal_virtualNetworksGet(_api::VirtualNetworksApi, resourceGroupName::String, virtualNetworkName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", VirtualNetwork, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/virtualNetworks/{virtualNetworkName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkName", virtualNetworkName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_param(_ctx.query, "\$expand", expand)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function virtualNetworksGet(_api::VirtualNetworksApi, resourceGroupName::String, virtualNetworkName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworksGet(_api, resourceGroupName, virtualNetworkName, api_version, subscriptionId; expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function virtualNetworksGet(_api::VirtualNetworksApi, response_stream::Channel, resourceGroupName::String, virtualNetworkName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworksGet(_api, resourceGroupName, virtualNetworkName, api_version, subscriptionId; expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets all virtual networks in a resource group.
Param: resourceGroupName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: VirtualNetworkListResult
"""
function _swaggerinternal_virtualNetworksList(_api::VirtualNetworksApi, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", VirtualNetworkListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/virtualNetworks", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function virtualNetworksList(_api::VirtualNetworksApi, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworksList(_api, resourceGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function virtualNetworksList(_api::VirtualNetworksApi, response_stream::Channel, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworksList(_api, resourceGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets all virtual networks in a subscription.
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: VirtualNetworkListResult
"""
function _swaggerinternal_virtualNetworksListAll(_api::VirtualNetworksApi, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", VirtualNetworkListResult, "/subscriptions/{subscriptionId}/providers/Microsoft.Network/virtualNetworks", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function virtualNetworksListAll(_api::VirtualNetworksApi, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworksListAll(_api, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function virtualNetworksListAll(_api::VirtualNetworksApi, response_stream::Channel, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworksListAll(_api, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Updates a virtual network tags.
Param: resourceGroupName::String (required)
Param: virtualNetworkName::String (required)
Param: parameters::TagsObject (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: VirtualNetwork
"""
function _swaggerinternal_virtualNetworksUpdateTags(_api::VirtualNetworksApi, resourceGroupName::String, virtualNetworkName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PATCH", VirtualNetwork, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/virtualNetworks/{virtualNetworkName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkName", virtualNetworkName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function virtualNetworksUpdateTags(_api::VirtualNetworksApi, resourceGroupName::String, virtualNetworkName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworksUpdateTags(_api, resourceGroupName, virtualNetworkName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function virtualNetworksUpdateTags(_api::VirtualNetworksApi, response_stream::Channel, resourceGroupName::String, virtualNetworkName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworksUpdateTags(_api, resourceGroupName, virtualNetworkName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export virtualNetworksCreateOrUpdate, virtualNetworksDelete, virtualNetworksGet, virtualNetworksList, virtualNetworksListAll, virtualNetworksUpdateTags


text:
"""
    readop4(filename::String)

Searc file and find matrix names (no matname) in Nastran binary OP4 file.
"""
function readop4(filename::String)
    matrices = readop4(filename, MAT_NAMES)

    return collect(keys(matrices))
end

"""
    readop4(filename::String, matname::String)

Read specified matrix `matname` from Nastran binary OP4 file into sparse matrices.
"""
function readop4(filename::String, matname::String)
    file = FortranFiles.FortranFile(filename)

    matrices = Dict{String, Union{SparseMatrixCSC, Nothing}}()

    while !eof(file.io)
        (ncol, nrow, mform, mtype, mname) = FortranFiles.read(file, Int32, Int32, Int32, Int32, FortranFiles.FString{8})
        mname = FortranFiles.trimstring(mname)

        colptr = zeros(Int32, ncol + 1)#Vector{Int32}(ncol + 1)
        rowind = Vector{Int32}(undef, 0)
        data = Vector{Float64}(undef, 0)

        # pre allocate and then resize!, maybe not worth the effort
        sizehint!(rowind, 2ncol)
        sizehint!(data, 2ncol)

        icol_count::Int32 = 1
        if mname == matname || matname == MAT_ALL
            icol_count, matrices[mname] = readop4!(file, data, icol_count, colptr, rowind, ncol, nrow)
        else
            # not sough matrix, keep looking
            spool!(file, ncol)

            if matname == MAT_NAMES
                matrices[mname] = nothing
            end
        end
    end

    if isempty(matrices)
        println("No matrices found, try searching the file for available matrices.")
    end

    return matrices
end

"""
    readop4!(
        file::FortranFiles.FortranFile,
        data::Vector,
        icol_count::Int32,
        colptr::Vector,
        rowind::Vector,
        ncol::Int32,
        nrow::Int32,
    )

Read Nastran binary OP4 file and populate sparse matrix.
"""
function readop4!(
    file::FortranFiles.FortranFile,
    data::Vector,
    icol_count::Int32,
    colptr::Vector,
    rowind::Vector,
    ncol::Int32,
    nrow::Int32,
)
    while true
        @fread file icol::Int32 irow::Int32 nword::Int32 iline::Array{Int32}(undef, nword)

        p::Int32 = 1
        m::Int32 = 1

        if icol > ncol
            # read one extra line after end to terminate matrix
            colptr[ncol + 1] = length(rowind) + 1
            fillcolptr!(colptr, icol, icol_count)
            return icol_count, SparseMatrixCSC(nrow, ncol, convert(Vector{Int}, colptr), convert(Vector{Int}, rowind), data)
        end

        if ncol < COL #|| nrow > 0 # small matrix, comment for now version diff.
            p = readsmallop4!(data, iline, p, nword, colptr, icol, rowind)
        else # big matrix
            p, m = readbigop4!(data, iline, p, nword, colptr, icol, rowind)
        end

        # correct col_ptr
        fillcolptr!(colptr, icol, icol_count)
        icol_count = icol + 1
    end
end

"""
    spool(file::FortranFiles.FortranFile, ncol::Int32)

If nothing found, keep looking for data.
"""
function spool!(file::FortranFiles.FortranFile, ncol::Int32)
    icol2::Int32 = 1
    while icol2 < ncol + 1
        # spool forward
        icol2 = FortranFiles.read(file, Int32)
    end

    return nothing
end

"""
    readsmallop4!(data, iline, p, nword, colptr, icol, rowind)

Read small OP4 matrix.
"""
function readsmallop4!(data, iline, p, nword, colptr, icol, rowind)
    while p <= nword
        if p == 1
            colptr[icol] = length(rowind) + 1
        end
        slength = getlength(iline[p])
        wsize = Int32(slength/2)
        irow = getirow(iline[p])
        #println("p: ", p, ", irow: ", irow, ", slength: ", slength, ", wsize: ", wsize, ", nword: ", nword)

        append!(rowind, irow:(irow + wsize - 1))
        append!(data, reinterpret(Float64, iline[(p + 1):(p + slength)]))

        p = p + slength + 1
    end

    return p
end

"""
    readbigop4!(data, iline, p, nword, colptr, icol, rowind)

Read big OP4 matrix.
"""
function readbigop4!(data, iline, p, nword, colptr, icol, rowind)
    line = reinterpret(Float64, iline)
    while p <= nword
        if p == 1
            colptr[icol] = length(rowind) + 1
        end
        slength = iline[p] - 1
        irow = iline[p+1]
        wsize = Int32(slength/2)
        #println("p: ", p, ", m: " , m, ", irow: ", irow, ", slength: ", slength, ", wsize: ", wsize)

        append!(rowind, irow:(irow + wsize - 1))
        append!(data, line[(m + 1):(m + wsize)])

        p = p + slength + 2
        m = m + wsize + 1
    end

    return p, m
end

"""
    getlength(head::Int32)

Get record length, see Nastran documentation.
"""
getlength(head::Int32) = trunc(Int32, head*iCOL) - 1

"""
    getirow(head::Int32)

Get row index, see Nastran documentation.
"""
getirow(head::Int32) = head - COL*(getlength(head) + 1)

"""
    fillcolptr!(colptr, icol, icol_count)

Fill column pointer, see Nastran documentation.
"""
function fillcolptr!(colptr, icol, icol_count)
    while icol_count < icol
        colptr[icol_count] = colptr[icol]
        icol_count = icol_count + 1
    end

    return nothing
end

max_stars_repo_name: mgcth/SDNastran

max_issues_repo_name: mgcth/SDNastran

max_forks_repo_name: mgcth/SDNastran


text:
using Cerberus, Gurobi, GZip, Justitia, MathOptInterface
const MOI = MathOptInterface

using Cerberus, Gurobi, Justitia

const INSTANCES = Dict(
    name => Justitia.MIPLIBInstance(name) for name in [
        "30n20b8",
        "50v-10",
        "academictimetablessmall",
        "air05",
        "app1-1",
        "app1-2",
        "assign1-5-8",
        "atlanta-ip",
        "b1c1s1",
        "beasleyC3",
        "binkar10_1",
        "blp-ar98",
        "blp-ic98",
        "bnatt400",
        "bppc4-08",
        "brazil3",
        "buildingenergy",
    ]
)

const TIME_LIMIT_SEC = 60.0

const APPROACHES = Dict(
    "cerberus-no-incrementalism" =>
        Justitia.MOIBasedApproach{Cerberus.Optimizer}(
            () -> begin
                model = Cerberus.Optimizer(
                    Cerberus.AlgorithmConfig(
                        Cerberus.AlgorithmConfig(
                            warm_start_strategy = Cerberus.NO_WARM_STARTS,
                            model_reuse_strategy = Cerberus.NO_MODEL_REUSE,
                        ),
                    ),
                )
                MOI.set(model, MOI.TimeLimitSec(), TIME_LIMIT_SEC)
                return model
            end,
        ),
    "cerberus-only-warm-start" =>
        Justitia.MOIBasedApproach{Cerberus.Optimizer}(
            () -> begin
                model = Cerberus.Optimizer(
                    Cerberus.AlgorithmConfig(
                        warm_start_strategy = Cerberus.WHENEVER_POSSIBLE,
                        model_reuse_strategy = Cerberus.NO_MODEL_REUSE,
                    ),
                )
                MOI.set(model, MOI.TimeLimitSec(), TIME_LIMIT_SEC)
                return model
            end,
        ),
    "cerberus-full-incrementalism" =>
        Justitia.MOIBasedApproach{Cerberus.Optimizer}(
            () -> begin
                model = Cerberus.Optimizer(
                    Cerberus.AlgorithmConfig(
                        Cerberus.AlgorithmConfig(
                            warm_start_strategy = Cerberus.WHENEVER_POSSIBLE,
                            model_reuse_strategy = Cerberus.USE_SINGLE_MODEL,
                        ),
                    ),
                )
                MOI.set(model, MOI.TimeLimitSec(), TIME_LIMIT_SEC)
                return model
            end,
        ),
    "gurobi" => Justitia.MOIBasedApproach{Gurobi.Optimizer}(
        () -> begin
            model = Gurobi.Optimizer()
            MOI.set(model, MOI.TimeLimitSec(), TIME_LIMIT_SEC)
            return model
        end,
    ),
)

result_table = Justitia.CSVRecord("results.csv", Justitia.MILPResult)
Justitia.run_experiments!(
    result_table,
    INSTANCES,
    APPROACHES,
    Justitia.MILPResult,
)


text:
# Plotting functions
using Random

function interval_plot(results, name)
  #comp = results.model.component_dict[name]
  ### TODO fix this
  idx = 0
  for (i, p) in enumerate(results.model.params)
    if p.name == name
      idx = i
    end
  end
  comp = results.model.params[idx]
  # Check if these are already computed, if not compute
  if length(comp.likelihood_x) == 0
    compute_profiled_uncertainties!(results, 0.68)
  end
  x, y = comp.likelihood_x, comp.likelihood_y
  low, high = comp.low, comp.high
  mid = comp.fit
  ns, ps = abs(mid - low), abs(mid - high)
  labeluse = "$name: " *
    " \$\\mu = $(floor(mid*100)/100)" *
    "^{+$(floor(ps*100)/100)}" *
    "_{-$(ns)}\$"
  plot(x, y, label=labeluse)
  vlines(low, 0, y[x .>= low][1], color="red")
  vlines(high, 0, y[x .>= high][1], color="red")
  ylim(0, 1.2)
  xlim(minimum(x), maximum(x))
  legend(loc=2, frameon=false)
end

function _mean_comp(x, y)
  sum(x.*y)
end

function _variance(x, y, μ)
  sum( (x .- μ).^2 .*y )
end

function _covariance(x, y, μx, μy, z)
  N = sum(z)
  k = transpose(x .- μx).*(y .- μy).*z
  sum( k )/N
end

function correlation_plots(results; kwargs...)
  kwargs = Dict(kwargs)
  steps = get(kwargs, :steps, 200)
  figwidth = get(kwargs, :figwidth, 10)
  ## Everything stored in results
  M = length(results.model.params)
  fig, ax = subplots(nrows=M, ncols=M, figsize=(M*figwidth, M*figwidth))

  ## We will be lazy for now and robust later. Lazy way, assume we know the range
  # What does a slice in space look like?
  # We can optimize again, with a constraint in place (bounds)
  levels = exp.(-[Inf,11.83,6.18,2.3,0].^0.5./2)

  for (i, v) in enumerate(results.model.params)
    for (j, w) in enumerate(results.model.params)
      if j < M
        ax[j, i].set_xticklabels([])
      end
      if i > 1
        ax[j, i].set_yticklabels([])
      end
      if i > j
        ax[j, i].axis("off")
        continue
      end
      # profile this parameter
      nlow = copy(results.lower_bounds)
      nhigh = copy(results.upper_bounds)
      p0 = copy(results.min_parameters)
      nll = []
      if i != j
        A = range(minimum(v.likelihood_x), length=steps, stop=maximum(v.likelihood_x)) |> collect
        B = range(minimum(w.likelihood_x), length=steps+5, stop=maximum(w.likelihood_x)) |> collect
        function getnll(a, b)
          nlow[i], nlow[j] = a, b
          nhigh[i], nhigh[j] = a, b
          results.opt.lower_bounds = nlow
          results.opt.upper_bounds = nhigh
          #p0 = copy(results.min_parameters) << This is slow
          p0[i], p0[j] = a, b
          minf, minx, ret = optimize!(results.opt, p0)
          if minf == 0
            minf = 10000
          end
          exp(-minf)
        end
        Z = [getnll(a, b) for b in B, a in A]
        Z = Z / maximum(Z)

        ## Given Z, A, B
        yb = sum(Z, dims=2)
        ya = transpose( sum(Z, dims=1) )
        ya = ya / sum(ya)
        yb = yb / sum(yb)

        mean_a = _mean_comp(A, ya)
        mean_b = _mean_comp(B, yb)
        var_a = sqrt( _variance(A, ya, mean_a) )
        var_b = sqrt( _variance(B, yb, mean_b) )
        cov = _covariance(A, B, mean_a, mean_b, Z)
        cor = cov / var_a / var_b
        ax[i, j].text(0,0.5, @sprintf("Cor: %0.3f", cor))
        ax[j,i].contourf(A, B, Z, levels=levels)
        continue
      end
      ax[i, i].plot( v.likelihood_x, v.likelihood_y )
      ax[i, i].set_title(v.name)
    ax[i,i].set_ylim(0, 1.2)
    end
  end
end

function correlation_plots2(results; kwargs...)
  kwargs = Dict(kwargs)
  steps = get(kwargs, :steps, 200)
  figwidth = get(kwargs, :figwidth, 10)
  ## Everything stored in results
  M = length(results.model.params)
  fig, ax = subplots(nrows=M, ncols=M, figsize=(M*figwidth, M*figwidth))

  levels = exp.(-[Inf,11.83,6.18,2.3,0].^0.5./2)

  for (i, v) in enumerate(results.model.params)
    for (j, w) in enumerate(results.model.params)
      if j < M
        ax[j, i].set_xticklabels([])
      end
      if i > 1
        ax[j, i].set_yticklabels([])
      end
      if i > j
        ax[j, i].axis("off")
        continue
      end
      # profile this parameter
      nlow = copy(results.lower_bounds)
      nhigh = copy(results.upper_bounds)
      p0 = copy(results.min_parameters)
      nll = []
      if i != j
        ## Here we need to go quadrant by quadrant
        xl = range( v.fit, length=steps, stop=minimum(v.likelihood_x) ) |> collect
        xr = range( v.fit, length=steps, stop=maximum(v.likelihood_x) ) |> collect
        yu = range( w.fit, length=steps, stop=maximum(w.likelihood_x) ) |> collect
        yd = range( w.fit, length=steps, stop=minimum(w.likelihood_x) ) |> collect

        function getz(xset, yset)
          # We always start at the best fit, so lets reset p0 to best fit
          coords = []
          p0 = copy(results.min_parameters)
          minf0 = copy(results.min_objective)
          pHold = copy(p0)
          for y in yset
            rowz = Array{Float64}(undef, 0)
            p0 = copy(pHold)
            reset = true
            for x in xset
              nlow[i], nlow[j] = x, y
              nhigh[i], nhigh[j] = x, y
              results.opt.lower_bounds = nlow
              results.opt.upper_bounds = nhigh
              p0[i], p0[j] = x, y
              minf, minx, ret = optimize!(results.opt, p0)
              if minf == 0
                minf = 10000
              end
              push!(rowz, exp(-(minf - minf0)) )
              if reset
                reset = false
                pHold = copy(p0)
              end
            end
            push!(coords, rowz)
          end
          return coords
        end
        z1 = getz(xr, yu)
        z2 = getz(xl, yu)
        z3 = getz(xl, yd)
        z4 = getz(xr, yd)
        # Reorder
        z2 = [reverse(a) for a in z2]
        z3 = reverse([reverse(a) for a in z3])
        z4 = reverse(z4)
        zleft = vcat(z3, z2)
        zright = vcat(z4, z1)
        X = vcat(reverse(xl), xr)
        Y = vcat(reverse(yd), yu)
        Z = Array([vcat(i, j) for (i,j) in zip(zleft, zright)])
        Z = hcat(Z...)'
        Z = Z / maximum(Z)

        ## Given Z, A, B
        yb = sum(Z, dims=2)
        ya = transpose( sum(Z, dims=1) )
        ya = ya / sum(ya)
        yb = yb / sum(yb)

        mean_a = _mean_comp(X, ya)
        mean_b = _mean_comp(Y, yb)
        var_a = sqrt( _variance(X, ya, mean_a) )
        var_b = sqrt( _variance(Y, yb, mean_b) )
        cov = _covariance(X, Y, mean_a, mean_b, Z)
        cor = cov / var_a / var_b
        ax[i, j].text(0,0.5, @sprintf("Cor: %0.3f", cor))
        ax[j,i].contourf(X, Y, Z, levels=levels)
        continue
      end
      ax[i, i].plot( v.likelihood_x, v.likelihood_y )
      ax[i, i].set_title(v.name)
    ax[i,i].set_ylim(0, 1.2)
    end
  end
end