text:
module Hyperparameters

using FilePathsBase: AbstractPath, /
using JSON
using Memento

export HYPERPARAMETERS, hyperparam, hyperparams, report_hyperparameters, save_hyperparam

const MODULE = @__MODULE__
const LOGGER = getlogger(MODULE)

"""
    HYPERPARAMETERS::Dict{Symbol, Any}

Collection of all the hyperparameters, and their values accessed during this run.
"""
const HYPERPARAMETERS = Dict{Symbol, Any}()
const SAGEMAKER_PREFIX = "SM_HP_"

_name_to_envvar(prefix, name) = uppercase(string(prefix, name))
function _envvar_to_name(prefix, name)
    startswith(name, prefix) || throw(DomainError(name, "Should have prefix; $prefix"))
    return Symbol(lowercase(SubString(name, length(prefix) + 1)))
end

# Grabs hyperparameter from environment variables
_get_hyperparam(name::Symbol, prefix::AbstractString)= ENV[_name_to_envvar(prefix, name)]

# Records a hyperparameter into global HYPERPARAMETERS dict for later reporting
# Logs if a hyperparameter is set twice to two difference values
function _set_hyperparam(name::Symbol, value)
    if haskey(HYPERPARAMETERS, name)
        orig = HYPERPARAMETERS[name]
        if orig != value
            notice(LOGGER, "Overwriting HYPERPARAMETERS[:$name] from $(repr(orig)) to $(repr(value)).")
        elseif typeof(orig) != typeof(value)
            notice(
                LOGGER,
                "Changing type of HYPERPARAMETERS[:$name] from $(typeof(orig)) to $(typeof(value))."
            )
        end
    end
    HYPERPARAMETERS[name] = value
end

"""
    hyperparam([T::Type,] name; prefix="$SAGEMAKER_PREFIX"))

Load the hyperparameter with the given `name` from the environment variable
named with the name in uppercase, and prefixed with `prefix`.
If the type `T` is passed then it will be parsed as that type,
otherwise it will take a guess at the type chosing the first that passes successfully
from `Bool`, `Int`, then `Float` then `String`

Also stores the hyperparameter and its value in the global `HYPERPARAMETERS` dictionary.
This function is generally expected to be used with SageMaker, and supplies the default prefix for it.
```jldoctest
using Hyperparameters
ENV["HP_POWER_LEVEL"] = "9001"
hyperparam(:power_level; prefix="HP_")

# output
9001.0
```
"""
function hyperparam(name::Symbol; prefix::AbstractString=SAGEMAKER_PREFIX)
    return hyperparam(nothing, name; prefix=prefix)
end

function hyperparam(T::Type, name::Symbol; prefix::AbstractString=SAGEMAKER_PREFIX)
    value = parse(T, _get_hyperparam(name, prefix))
    _set_hyperparam(name, value)
    return value
end

function hyperparam(::Type{String}, name::Symbol; prefix::AbstractString=SAGEMAKER_PREFIX)
    value = _get_hyperparam(name, prefix)
    _set_hyperparam(name, value)
    return value
end

function hyperparam(::Nothing, name::Symbol; prefix::AbstractString=SAGEMAKER_PREFIX)
    # this one is adaptive to try and work out the type from the string.
    value = _parse_hyper(_get_hyperparam(name, prefix))
    _set_hyperparam(name, value)
    return value
end

"""
    _parse_hyper(value_string)

Guesses the type of the `value_string` and parses as that.
"""
function _parse_hyper(value_string)
    return something(
        _strict_tryparse_bool(value_string),
        tryparse(Int, value_string),
        tryparse(Float64, value_string),
        value_string,
    )
end

# tryparse(Bool, "1") == true so can't use that. This is stricter
_strict_tryparse_bool(x) = x == "true" ? true : x == "false" ? false : nothing

"""
    hyperparams(names...; prefix="$SAGEMAKER_PREFIX")

As per [`hyperparam`](@ref), but taking multiple names and returning a `NamedTuple`.
```jldoctest
using Hyperparameters
ENV["$(SAGEMAKER_PREFIX)A"] = "5"
ENV["$(SAGEMAKER_PREFIX)B"] = "1.22"
hyperparams(:a, :b, types=[Int, Float64])

# output
(a = 5, b = 1.22)
```
Also stores the hyperparameters and their values in the global `HYPERPARAMETERS`
dictionary.
"""
function hyperparams(
    names::Vararg{Symbol};
    prefix::AbstractString=SAGEMAKER_PREFIX,
    types::Union{Nothing, AbstractVector{<:Type}}=nothing
)
    parameters = if types === nothing
        (name => hyperparam(name; prefix=prefix) for name in names)
    else
        length(names) == length(types) || throw(ArgumentError(
            "Number of `names` and `types` must be equal $(length(names)) != $(length(types))"
        ))
        (name => hyperparam(type, name; prefix=prefix) for (name, type) in zip(names, types))
    end

    return (; parameters...)
end

"""
    save_hyperparam(name::Symbol, value, prefix::AbstractString="")

Save value to the enviroment variables and the global `HYPERPARAMETERS` dictionary.
"""
function save_hyperparam(name::Symbol, value; prefix::AbstractString="")
    ENV[_name_to_envvar(prefix, name)] = string(value)
    _set_hyperparam(name, value)
end

"""
    report_hyperparameters(save_dir::AbstractPath; prefix="$SAGEMAKER_PREFIX")

Saves all hyperparameters to a JSON file named "hyperparameters.json" in the `save_dir`
and prints each key-value pair to the logger.

The hyperparameters are taken from the cached `HYPERPARAMETERS` dictionary of all that were
used, combined with any enviroment variables matching the prefix.
Where things occur in both, the cached dictionary takes precedence.
Hyperparameters read from enviroment variables are all recorded as strings.
(You can overwrite this by using them via `hyperparam(type, name)` before the report)

The regex to extract the components is: `hyperparameters: (?<key>)=(?<value>)`.
"""
function report_hyperparameters(save_dir::AbstractPath; prefix=SAGEMAKER_PREFIX)
    env_hypers = Dict(
        _envvar_to_name(prefix, k) => _parse_hyper(v)
        for (k, v) in ENV if startswith(k, prefix)
    )
    all_hypers = merge(env_hypers, HYPERPARAMETERS)
    for (key, value) in all_hypers
        info(LOGGER, "hyperparameters: $key=$value")
    end

    file = save_dir / "hyperparameters.json"
    info(LOGGER, "Report: saving at $file")

    open(file, "w") do fh
        JSON.print(fh, all_hypers, 4)  # 4 space indenting
    end

    return file
end

end # module


text:
"""
`module JuLIP.Solve`

Contains a few geometry optimisation routines, for now see
the help for these:

* `minimise!`
"""
module Solve

import Optim
import LineSearches

using Optim: OnceDifferentiable, optimize, ConjugateGradient, LBFGS
using LineSearches: BackTracking

using JuLIP: AbstractAtoms, Preconditioner, update!, Identity,
            dofs, energy, gradient, set_dofs!, set_constraint!, site_energies,
            Dofs, calculator, constraint, AbstractCalculator

using JuLIP.Potentials: SitePotential
using JuLIP.Preconditioners: Exp
using JuLIP.Constraints: FixedCell

# TODO: move to Optim 0.11, and exploit the new alphaguess functionality of
#       Linesearches. This gives awful error messages.

export minimise!

Ediff(V::AbstractCalculator, at::AbstractAtoms, Es0::Vector{Float64}) =
   sum_kbn(site_energies(V, at) - Es0)

Ediff(at::AbstractAtoms, Es0::Vector{Float64}, x::Dofs) =
   Ediff(calculator(at), set_dofs!(at, x), Es0)

"""
`minimise!(at::AbstractAtoms)`: geometry optimisation

`at` must have a calculator and a constraint attached.

## Keyword arguments:
* `precond = :auto` : preconditioner; more below
* `gtol = 1e-6` : gradient tolerance (max-norm)
* `ftol = 1e-32` : objective tolerance
* `Optimiser = :auto`, `:auto` should always work, at least on the master
   branch of `Optim`; `:lbfgs` needs the `extraplbfgs2` branch, which is not
   yet merged. Other options might be introduced in the future.
* `verbose = 1`: 0 : no output, 1 : final, 2 : iteration and final
* `robust_energy_difference = false` : if true use Kahan summation of site energies
* `store_trace = false` : store history of energy and norm of forces
* `extended_trace = false`: also store full history of postions and forces
* `maxstep = Inf`: maximum step size, useful if initial gradient is very large
* `callback`: callback function to pass to `optimize()`, e.g. to use alternate convergence criteria

## Preconditioner

`precond` may be a valid preconditioner, e.g.,
`Identity()` or `Exp(at)`, or one of the following symbols

* `:auto` : the code will make the best choice it can with the avilable
   information
* `:exp` : will use `Exp(at)`
* `:id` : will use `Identity()`
"""
function minimise!(at::AbstractAtoms;
                  precond = :auto,
                  method = :auto,
                  gtol=1e-5, ftol=1e-32,
                  verbose = 1,
                  robust_energy_difference = false,
                  store_trace = false,
		  extended_trace = false,
                  maxstep = Inf,
                  callback = nothing)

   # create an objective function
   if robust_energy_difference
      Es0 = site_energies(at)
      obj_f = x -> Ediff(at, Es0, x)
   else
      obj_f = x->energy(at, x)
   end
   obj_g! = (g, x) -> copy!(g, gradient(at, x))

   # create a preconditioner
   if isa(precond, Symbol)
      if precond == :auto
         if isa(constraint(at), FixedCell)
            precond = :exp
         else
            precond = :id
         end
      end
      if precond == :exp
         if method == :lbfgs
            precond = Exp(at, energyscale = :auto)
         else
            precond = Exp(at)
         end
      elseif precond == :id
         precond = Identity()
      else
         error("unknown symbol for precond")
      end
   end

   # choose the optimisation method Optim.jl
   if method == :auto || method == :cg
      if isa(precond, Identity)
         optimiser = ConjugateGradient(linesearch = BackTracking(order=2, maxstep=maxstep))
      else
         optimiser = ConjugateGradient( P = precond,
                           precondprep = (P, x) -> update!(P, at, x),
                           linesearch = BackTracking(order=2, maxstep=maxstep) )
      end
   elseif method == :lbfgs
      optimiser = LBFGS( P = precond,
                        precondprep = (P, x) -> update!(P, at, x),
                        alphaguess = LineSearches.InitialHagerZhang(),
                        linesearch = BackTracking(order=2, maxstep=maxstep) )
   else
      error("JuLIP.Solve.minimise!: unknown `method` option")
   end

   results = optimize( obj_f, obj_g!, dofs(at), optimiser,
                        Optim.Options( f_tol = ftol, g_tol = gtol,
                                       store_trace = store_trace,
                                       extended_trace = extended_trace,
                                       callback = callback,
                                       show_trace = (verbose > 1)) )
   set_dofs!(at, Optim.minimizer(results))
   # analyse the results
   if verbose > 0
      println(results)
   end
   return results
end

end

ext: jl


text:
using DataFrames
using NamedArrays
using SchumakerSpline

bounds = (0.0, 1.0)
x_name = :alpha
evaluations = DataFrame()

function func_inc1(x)
    feval = NamedArray([x, -sqrt(x)])
    setnames!(feval, ["utility", "x_E"], 1)
    return feval
end
const_Array = NamedArray([0.1, sqrt(0.2)])
setnames!(const_Array, ["utility", "x_E"], 1)

# Testing the construction of a dataframe of functions
dd_function = make_function_dataframe([:func_inc1], Array{Any,1}([func_inc1]), [:Increasing])
dd_function2 = make_function_dataframe([:const_Array], Array{Any,1}([const_Array]), [:Constant])

#
grid                = reverse(0.0001:0.01:0.9999)
output_for_envelope = :utility
x_name              = :alpha
env                 = Envelope(dd_function, grid, output_for_envelope,x_name)
isa(env, Envelope)
env2                = Envelope(dd_function2, grid, output_for_envelope,x_name)
isa(env2, Envelope)


text:
#=
P10 (*) Run-length encoding of a list.
Use the result of problem P09 to implement the so-called run-length encoding data compression method. Consecutive duplicates of elements are encoded as lists (N E) where N is the number of duplicates of the element E.

Example:
* (encode '(a a a a b c c a a d e e e e))
((4 A) (1 B) (2 C) (2 A) (1 D)(4 E))
=#

# From P09
function _pack(coll, current_group, finished_groups)
    if isempty(coll)
        vcat(finished_groups, Vector{}[current_group])
    else
        if coll[1] == current_group[end]
            _pack(coll[2:end], vcat(current_group, coll[1]), finished_groups)
        else
            _pack(coll[2:end], [coll[1]], vcat(finished_groups, Vector{}[current_group]))
        end
    end
end

# From P09
function pack(coll)
    first, rest = coll[1], coll[2:end]
    _pack(rest, [first], [])
end

function encode(coll)
    packed = pack(coll)
    map(group -> [length(group), group[1]], packed)
end

encoded = encode([:a, :a, :a, :a, :b, :c, :c, :a, :a, :d, :e, :e, :e, :e])
@assert encoded[1] == [4, :a]
@assert encoded[2] == [1, :b]
@assert encoded[3] == [2, :c]
@assert encoded[4] == [2, :a]
@assert encoded[5] == [1, :d]
@assert encoded[6] == [4, :e]

println("Tests passed: JL-10.jl")


text:
using Astro, Test

print("37.a (1) heliocentric coordinates of Pluto")
L, B, R = heliocentric_pluto(2448908.5)
@test isapprox(round(B, digits=5), 14.58782, atol=0.1)
@test isapprox(round(L, digits=5), 232.74071, atol=0.1)
@test isapprox(round(R, digits=6), 29.71111, atol=0.1)
println(" passed")


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/BoundingSphere.jl/blob/master/LICENSE

using StaticArrays

@testset "Generic type support" begin
    inputs = []
    pts_template = [[-1,0],[1,0],[0,-1],[0,1]]
    for F in subtypes(AbstractFloat)
        F == Float16 && continue  # Float16 fires assertions
        for V in [Vector{F}, SVector{2,F}]
            pts::Vector{V} = map(V, pts_template)
            push!(inputs, pts)
        end
    end
    for A in subtypes(MB.BoundingSphereAlg)
        for pts in inputs
            alg = A()
            c, r = boundingsphere(pts, alg)
            P = eltype(pts)
            F = eltype(P)
            @test typeof(c) == P
            @test typeof(r) == F
            @inferred boundingsphere(pts, alg)
            @test norm(c) < sqrt(eps(F))
            @test r ≈ 1
        end
    end
end

@testset "support_count" begin
    alg = WelzlMTF()

    a = @SVector [0.]
    b = @SVector [1.]
    c = @SVector [2.]

    pts = [a]
    bdry = MB.create_boundary_device(pts, alg)
    ball, support_count = MB.welzl!(pts, bdry, alg)
    @test isempty(bdry)
    @test support_count == 1

    pts = [a,b]
    bdry = MB.create_boundary_device(pts, alg)
    ball, support_count = MB.welzl!(pts, bdry, alg)
    @test isempty(bdry)
    @test support_count == 2

    pts = [a,b,c]
    bdry = MB.create_boundary_device(pts, alg)
    ball, support_count = MB.welzl!(pts, bdry, alg)
    @test isempty(bdry)
    @test support_count == 2
    @test Set(pts[1:support_count]) == Set([a,c])

    vals = randn(10)
    V = SVector{1, Float64}
    pts = map(V, vals)
    v1, v2 = map(V, extrema(vals))
    bdry = MB.create_boundary_device(pts, alg)
    ball, support_count = MB.welzl!(pts, bdry, alg)
    @test isempty(bdry)
    @test support_count == 2
    @test Set(pts[1:support_count]) == Set([v1,v2])

end

@testset "random support_count" begin
    seed!(42)
    alg = WelzlMTF()
    for dim in 1:10, npoints in 1:20
        ball, pts = create_ball_points(dim, npoints, p_rep=0, p_boundary=0)
        bdry = MB.create_boundary_device(pts, alg)
        ball, support_count = MB.welzl!(pts, bdry, alg)
        @test isempty(bdry)

        # support set suffices
        pts2 = pts[1:support_count]
        c, r = boundingsphere(pts2, alg)
        @test MB.center(ball) ≈ c
        @test MB.radius(ball) ≈ r

        # support set is not too large
        pts3 = copy(pts2)
        if length(pts3) > 1
            index = rand(1:length(pts3))
            deleteat!(pts3, index)
            c, r = boundingsphere(pts3, alg)
            @test r < MB.radius(ball)
        end
    end
end

@testset "random Ritter" begin
    seed!(42)
    for dim in 1:5
        for npoints in 1:100
            random_test(Ritter(), npoints, dim,
                        rtol_radius=0.3,
                        allow_broken=false
                       )
        end
    end
end

@testset "random WelzlPivot" begin

    @testset "random WelzlPivot small" begin
        seed!(42)
        alg = WelzlPivot()
        for dim in 1:3
            for npoints in 1:10
                    random_test(alg, npoints, dim,
                                allow_broken=false,
                                codim=0)
            end
        end
    end

    @testset "non degenerate" begin
        seed!(42)
        alg = WelzlPivot()
        for dim in 1:10
            for npoints in 1:100
                    random_test(alg, npoints, dim,
                                p_boundary=0,
                                p_rep=0,
                                codim=0,
                                allow_broken=false)
            end
        end
    end

    @testset "nasty" begin
        seed!(42)
        alg = WelzlPivot()
        for _ in 1:1
            for dim in 1:10
                for npoints in 1:100
                    random_test(alg, npoints, dim,
                                p_rep=1/sqrt(npoints),
                                p_boundary=0.5,
                                allow_broken=true)

                    random_test(alg, npoints, dim,
                                p_rep=1/sqrt(npoints),
                                p_boundary=1,
                                allow_broken=true)

                    random_test(alg, npoints, dim,
                                p_rep=0,
                                p_boundary=0,
                                allow_broken=true)

                    random_test(alg, npoints, dim,
                                allow_broken=true)

                    random_test(alg, npoints, dim,
                                codim=poisson(3),
                                allow_broken=true)
                end
            end
        end
    end
end

@testset "random WelzlMTF" begin
    seed!(42)
    for dim in 1:10, npoints in 1:20
        random_test(WelzlMTF(), npoints, dim,
                    codim=0, p_rep=0, p_boundary=0,
                    allow_broken=false)
    end
end


text:
"""
The data recorded at the receivers during a seismic experiment
is a convolution of the Green's function,
the source signature and the receiver instrument
response.
The source filter `ssf` takes into account
the effective wavelet present in the data.
Similarly, the receiver filter `rf` takes the
receiver coupling into account.
The Inversion module can
can estimate both `ssf` and `rf` from the seismic data in order to achieve better data fitting.
A base source wavelet `Src` is always used to model the data
before applying the time-domain filter data type  described in this module:
* `TD` : time domain filters for both source and receivers
"""
module Coupling

import GeoPhyInv: AGeom

"""
Time-domain source and receiver filters.

# Fields

* `ssf::Array{Array{Float64,1},2}` : source filters for each supersource and recorded component
* `rf::Array{Array{Float64,2},2}` : receiver filters for each receiver, supersource and recorded component
* `fields::Vector{Symbol}` :  number of recorded components at receivers
* `tgridssf` : a  time grid for source filters with both positive and negative lags
* `tgridrf` : a  time grid for receiver filters with both positive and negative lags
* `ageom::AGeom` : acquisition ageometry
"""
mutable struct TD
	ssf::Array{Array{Float64,1},2}
	rf::Array{Array{Float64,2},2}
	fields::Vector{Symbol}
	tgridssf::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
	tgridrf::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
	ssflags::Vector{Int}
	rflags::Vector{Int}
	ageom::AGeom
	TD(ssf, rf, fields, tgridssf, tgridrf, ssflags, rflags, ageom) =
		any([
		  any([fields[ifield] ∉ [:P, :Vx, :Vz] for ifield in 1:length(fields)]),
		  length(fields) == 0,
		  broadcast(size,ssf) != [(length(tgridssf),) for iss=1:ageom.nss, ifield=1:length(fields)]
		  ]) ?
		error("error in TD construction") : new(ssf, rf, fields, tgridssf, tgridrf, ssflags, rflags, ageom)

end

"""
Initialize coupling filters `TD` with  delta functions.

# Arguments

* `tgriddata` : time grid of the data
* `tlagssf_fracs::Vector{Real}` : +ve and -ve fractions of source filter
* `tlagrf_fracs::Vector{Real}` : +ve and -ve fractions for receiver filter
* `δt:Float64` : sampling interval in time
* `fields::Vector{Symbol}` : number of components
* `ageom::AGeom` : acquisition ageometry

# Return

* time-domain coupling filters as `TD`
"""
function TD_delta(
			   tgriddata,
			   tlagssf_fracs,
			   tlagrf_fracs,
		  fields::Vector{Symbol},
		  ageom::AGeom
		 )
	δt=tgriddata.δx
	tot_t=abs(tgriddata[end]-tgriddata[1])

	error("fix these")
	#tgridssf, ssflags = lag(tot_t.*tlagssf_fracs, δt)
	#tgridrf, rflags =lag(tot_t.*tlagrf_fracs, δt)

	spiss = zeros(length(tgridssf)); spir = zeros(length(tgridrf));
	# check where to put spikes
	nspikessf=(findall(tgridssf.==0.0) == []) ? error("no t=0") : findall(tgridssf.==0.0)
	nspikerf=(findall(tgridrf.==0.0) == []) ? error("no t=0") : findall(tgridrf.==0.0)
	# put spikes
	spiss[nspikessf] = 1.0; spir[nspikerf] = 1.0;
	# number of unique receivers (implement in future)
	# one receiver filter for each unique receiver

	ssf=[spiss for iss=1:ageom.nss, ifield=1:length(fields)]

	return TD(ssf,
	   [repeat(spir, outer=[1,ageom.nr[iss]]) for iss=1:ageom.nss, ifield=1:length(fields)],
	   fields,tgridssf,tgridrf,ssflags,rflags,deepcopy(ageom))

end

end # module


text:
import Base.Iterators

abstract type Seq{T} end

struct EmptySeq{T} <: Seq{T} end

struct ConsSeq{T} <: Seq{T}
    genfun::Function
end

function Seq(::Type{T}, state, next::Function) where {T}
    ConsSeq{T}(() -> (state, next))
end

function Seq(::Type{T}, iter) where {T}
    ConsSeq{T}( () -> (
                       begin
                           nothing, state ->
                           begin
                               res = (state === nothing) ? iterate(iter) : iterate(iter, state)
                               if res === nothing
                                   (nothing::Union{T, Nothing}), state
                               else
                                   v, newstate = res
                                   (v::Union{T, Nothing}), newstate
                               end
                           end
                       end
              ))
end

function Base.eltype(x::Seq{T}) where {T}
    T
end

function Seq(range::AbstractRange{T}) where {T}
    Seq(T, range)
end

tryread(s::EmptySeq{T}) where {T} = nothing, EmptySeq{T}()

function tryread(xs::ConsSeq{T}) where {T}
    state, next = xs.genfun()
    v, newstate = next(state)
    if v === nothing
        nothing, EmptySeq{T}()
    else
        v, ConsSeq{T}(() -> (newstate, next))
    end
end

Base.map(f::Function, xs::EmptySeq{T}, ::Type{S}) where {T, S} = EmptySeq{S}()

Base.map(f::Function, xs::EmptySeq{T}) where {T} = EmptySeq{T}()

function Base.map(f::Function, xs::ConsSeq{T}, ::Type{S}) where {T, S}
    ConsSeq{S}(() ->
                begin
                    state, next = xs.genfun()
                    state, s -> begin
                                    v, newstate = next(s)
                                    if v === nothing
                                        nothing, s
                                    else
                                        f(v), newstate
                                    end
                                end
                end
              )
end

function Base.map(f::Function, xs::ConsSeq{T}) where {T}
    map(f, xs, T)
end

fold(f::Function, init, xs::EmptySeq{T}) where {T} = init

function fold(f::Function, init, xs::ConsSeq{T}) where {T}
    eof = false
    acc = init
    tail = xs
    while !eof
        v, tail = tryread(tail)
        if v === nothing
            eof = true
        else
            acc = f(acc, v)
        end
    end
    acc
end

function iter(f::Function, xs::ConsSeq{T}) where {T}
    eof = false
    tail = xs
    while !eof
        v, tail = tryread(tail)
        if v === nothing
            eof = true
        else
            f(v)
        end
    end
end

Iterators.take(xs::EmptySeq{T}, n::Integer) where {T} = EmptySeq{T}()

function Iterators.take(xs::ConsSeq{T}, n::Integer) where {T}
    ConsSeq{T}(() ->
                    begin
                        state, next = xs.genfun()
                        (state, n), s -> begin
                                             state = s[1]
                                             n = s[2]
                                             if n >= 1
                                                v, newstate = next(state)
                                                v, (newstate, n - 1)
                                             else
                                                 nothing, state
                                             end
                                    end
                    end
             )
end

function chunkbysize(xs::EmptySeq{T}, n::Integer) where {T}
    EmptySeq{Vector{T}}()
end

function chunkbysize(xs::ConsSeq{T}, n::Integer) where {T}
    return ConsSeq{Vector{T}}(() ->
        begin
            state, next = xs.genfun()
            state, s -> begin
                            chunk = Vector{T}(undef, n)
                            newstate = s
                            for k in 1:n
                                v, newstate = next(newstate)
                                if v !== nothing
                                    chunk[k] = v
                                else
                                    resize!(chunk, k - 1)
                                    break
                                end
                            end
                            if length(chunk) == 0
                                nothing, newstate
                            else
                                (chunk::Union{Vector{T}, Nothing}), newstate
                            end
                        end
        end)
end

function Base.zip(xs::NTuple{N, Seq}) where {N}
    eltypes = map(eltype, xs)
    EmptySeq{Tuple{eltypes...}}()
end

function Base.zip(xs::NTuple{N, ConsSeq}) where {N}
    eltypes = map(eltype, xs)
    return ConsSeq{Tuple{eltypes...}}(() ->
        begin
            y = map((x -> x.genfun()), xs)
            state = map((x -> x[1]), y)
            next = map((x -> x[2]), y)
            state, sarg ->
                begin
                    res = map((x -> x[2](x[1])), zip(sarg, next))
                    newstate = map((x -> x[2]), res)
                    if all([a !== nothing for (a, b) in res])
                        v = map((x -> x[1]), res)
                        (Tuple(v)::Union{Tuple{eltypes...}, Nothing}), newstate
                    else
                        nothing, newstate
                    end
                end
        end
    )
end

function Base.zip(xs::Seq...)
    Base.zip(xs)
end

function Base.zip(xs::ConsSeq...)
    Base.zip(xs)
end

function Base.foreach(f::Function, xs::Seq{T}) where {T}
    eof = false
    tail = xs
    while !eof
        v, tail = tryread(tail)
        if v === nothing
            eof = true
        else
            f(v)
        end
    end
end

function Base.collect(xs::Seq{T}) where {T}
    res = Vector{T}()
    eof = false
    tail = xs
    while !eof
        v, tail = tryread(tail)
        if v === nothing
            eof = true
        else
            push!(res, v)
        end
    end
    res
end

function concat(xs::EmptySeq{T}, ys::Seq{T}) where {T}
    ys
end

function concat(xs::Seq{T}, ys::EmptySeq{T}) where {T}
    xs
end

function concat(xs::ConsSeq{T}, ys::ConsSeq{T}) where {T}
    return ConsSeq{T}(() ->
            begin
                statex, nextx = xs.genfun()
                statey, nexty = ys.genfun()
                (statex, statey), s ->
                    begin
                        _statex, _statey = s
                        x, newstatex = nextx(_statex)
                        if x === nothing
                            y, newstatey = nexty(_statey)
                            if y === nothing
                                nothing, (newstatex, newstatey)
                            else
                                (y::Union{T, Nothing}), (newstatex, newstatey)
                            end
                        else
                            (x::Union{T, Nothing}), (newstatex, _statey)
                        end
                    end
            end
            )
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
using Test
using jInvVis
using jInv.Mesh

# tests for regular mesh
domain = [0 1.1 0 1.0 0  1.1]
n      = [8 5 3]
Mr    = getRegularMesh(domain,n)

xc    = getCellCenteredGrid(Mr)
println("=== test viewSlice2D  ===")
figure(1); clf()
subplot(1,3,1)
viewSlice2D(xc[:,1],Mr,Int(round(n[3]/2)))
xlabel("x, intensity increases --> ")
ylabel("intensity constant")
title("f(x,y,z) = x")

subplot(1,3,2)
viewSlice2D(xc[:,1],Mr,Int(round(n[2]/2)),view=:xz,addLabel=true)
xlabel("x, intensity increases --> ")
title("f(x,y,z) = x")

subplot(1,3,3)
viewSlice2D(xc[:,2],Mr,Int(round(n[1]/2)),view=:yz,addLabel=true)
xlabel("y, intensity increases --> ")
title("f(x,y,z) = y")


text:
using Test
using Comonicon

# pkg activate <string> --shared=false

module PkgCmd

using Test
using Comonicon

"""
activate the environment at `path`.

# Arguments

- `path`: the path of the environment

# Flags

- `--shared`: whether activate the shared environment
"""
@cast function activate(path; shared::Bool = false)
    @test path == "path"
    @test shared == true
end

@main

end

PkgCmd.command_main()


text:
const constant_funcs = []

const monadic_linear = [deg2rad, +, rad2deg, transpose, -, conj]

const monadic_nonlinear = [asind, log1p, acsch, acos, asec, acosh, acsc, cscd, log, tand, log10, csch, asinh, abs2, cosh, sin, cos, atan, cospi, cbrt, acosd, acoth, inv, acotd, asecd, exp, acot, sqrt, sind, sinpi, asech, log2, tan, exp10, sech, coth, asin, cotd, cosd, sinh, abs, csc, tanh, secd, atand, sec, acscd, cot, exp2, expm1, atanh]

const diadic_linear_true_true_true = [+, rem2pi, -, >, isless, <, isequal, max, min, convert]
const diadic_linear_true_true_false = [*]
const diadic_linear_true_false_false = [ / ]
const diadic_linear_false_true_false = [ \ ]
const diadic_linear_false_false_false = [hypot, atan, mod, rem, ^]

haslinearity(f, nargs) = false

# some functions strip the linearity metadata

for f in constant_funcs
    @eval begin
        haslinearity(::typeof($f), ::Val) = true
        linearity(::typeof($f), ::Val) = nothing
    end
end

# linearity of a single input function is either
# Val(true) or Val(false}
#
for f in monadic_linear
    @eval begin
        haslinearity(::typeof($f), ::Val{1}) = true
        linearity(::typeof($f), ::Val{1}) = Val(true)
    end
end

for f in monadic_nonlinear
    @eval begin
        haslinearity(::typeof($f), ::Val{1}) = true
        linearity(::typeof($f), ::Val{1}) = Val(false)
    end
end

# linearity of a 2-arg function is:
# Val((linear11, linear22, linear12))
#
# linearIJ refers to the zeroness of d^2/dxIxJ
for f in diadic_linear_true_true_true
    @eval begin
        haslinearity(::typeof($f), ::Val{2}) = true
        linearity(::typeof($f), ::Val{2}) = Val((true, true, true))
    end
end

for f in diadic_linear_true_true_false
    @eval begin
        haslinearity(::typeof($f), ::Val{2}) = true
        linearity(::typeof($f), ::Val{2}) = Val((true, true, false))
    end
end

for f in diadic_linear_true_false_false
    @eval begin
        haslinearity(::typeof($f), ::Val{2}) = true
        linearity(::typeof($f), ::Val{2}) = Val((true, false, false))
    end
end

for f in diadic_linear_false_true_false
    @eval begin
        haslinearity(::typeof($f), ::Val{2}) = true
        linearity(::typeof($f), ::Val{2}) = Val((false, true, false))
    end
end

for f in diadic_linear_false_false_false
    @eval begin
        haslinearity(::typeof($f), ::Val{2}) = true
        linearity(::typeof($f), ::Val{2}) = Val((false, false, false))
    end
end


text:
# Background:
#
# We build a novel fe space form a given fe space plus a set of linear constraints.
# We accept any single field fe space as input.  In particular, the given fe space can also be defied
# via constraints.
#
# Assumptions:
#
#  - A constrained dof depends only on master dofs, i.e., only one level of constraints.
#    This can be easily relaxed in the future.
#
#  - We can add constraints to Dirichlet dofs but their masters have to be Dirichlet as well.
#    However, free dofs can be constrained both by free and Dirichlet.
#
# Notation in this file:
#
#  - dof: a dof (either free or Dirichlet) in the original space.
#   The sign of the value tells if it is one or the other.
#
#  - DOF: a dof (either free or Dirichlet) in the original space.
#    In this case, Dirichlet dofs are also represented with a positive id "pas the end" of free dof ids.
#    Useful for indexing arrays.
#    It will refer to a free dof if 0 < DOF <= n_fdofs, and Dirichlet if n_fdofs < DOF <= n_dofs
#    with n_ddofs=n_dofs-n_fdofs
#
#  - fdof: a free dof in the original space. We have: fdof = dof = DOF
#
#  - ddof: a Dirichlet dof in the original space. We have: ddof = -dof = DOF - n_fdofs
#
#  - sDOF a slave dof (either free or Dirichlet). It is possible to tell if it is free or Dirichlet
#    thanks to the vector sDOF_to_dof.
#
#  - mdof: a master dof (i.e. a dof in the constrained space).
#    It can be free or Dirichlet. The sign of the value tells if it is one or the other.
#
#  - mDOF: a master dof (i.e. a dof in the constrained space).
#    In this case, Dirichlet dofs are also represented with a positive id "pas the end" of free dof ids.
#    It will refer to a free dof if 0 < mDOF <= n_fmdofs, and Dirichlet if n_fmdofs < mDOF <= n_mdofs
#    with n_dmdofs=n_mdofs-n_fmdofs
#
#  - fmdof: a free master dof. We have fmdof = mdof = mDOF
#
#  - dmdof: a Dirichlet master dof. We have dmdof = -mdof = mDOF - n_fmdofs
#
#  - ldof: a local dof in a cell for the original space. Always positive.
#
#  - lmdof: local master dof in a cell. That is, a local dofs in a cell for the constrained space.
#    Always positive.
#
#  - ludof: A local dof in the unconstrained interpolation space of the underlying cell. Always positive.
#
#  - coeff: A coefficient in a linear constraint

# TODO we can optimize memory by only storing info about slave DOFs
# The fields of this struct are private
struct FESpaceWithLinearConstraints{S<:SingleFieldFESpace} <: SingleFieldFESpace
  space::S
  n_fdofs::Int
  n_fmdofs::Int
  mDOF_to_DOF::Vector
  DOF_to_mDOFs::Table
  DOF_to_coeffs::Table
  cell_to_lmdof_to_mdof::Table
  cell_to_ldof_to_dof::Table
end

# Constructors

# This is the public constructor
function FESpaceWithLinearConstraints(
  sDOF_to_dof::AbstractVector{<:Integer},
  sDOF_to_dofs::Table,
  sDOF_to_coeffs::Table,
  space::SingleFieldFESpace)

  n_fdofs = num_free_dofs(space)
  n_ddofs = num_dirichlet_dofs(space)
  n_DOFs = n_fdofs + n_ddofs

  DOF_to_DOFs, DOF_to_coeffs = _prepare_DOF_to_DOFs(
    sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, n_fdofs, n_DOFs)

  FESpaceWithLinearConstraints!(DOF_to_DOFs,DOF_to_coeffs,space)

end

function _prepare_DOF_to_DOFs(
  sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, n_fdofs, n_DOFs)

  Tp = eltype(sDOF_to_dofs.ptrs)
  Td = eltype(sDOF_to_dofs.data)
  Tc = eltype(sDOF_to_coeffs.data)

  DOF_to_DOFs_ptrs = ones(Tp,n_DOFs+1)

  n_sDOFs = length(sDOF_to_dof)

  for sDOF in 1:n_sDOFs
    a = sDOF_to_dofs.ptrs[sDOF]
    b = sDOF_to_dofs.ptrs[sDOF+1]
    dof = sDOF_to_dof[sDOF]
    DOF = _dof_to_DOF(dof,n_fdofs)
    DOF_to_DOFs_ptrs[DOF+1] = b-a
  end

  length_to_ptrs!(DOF_to_DOFs_ptrs)
  ndata = DOF_to_DOFs_ptrs[end]-1
  DOF_to_DOFs_data = zeros(Td,ndata)
  DOF_to_coeffs_data = ones(Tc,ndata)

  for DOF in 1:n_DOFs
    q = DOF_to_DOFs_ptrs[DOF]
    DOF_to_DOFs_data[q] = DOF
  end

  for sDOF in 1:n_sDOFs
    dof = sDOF_to_dof[sDOF]
    DOF = _dof_to_DOF(dof,n_fdofs)
    q = DOF_to_DOFs_ptrs[DOF]-1
    pini = sDOF_to_dofs.ptrs[sDOF]
    pend = sDOF_to_dofs.ptrs[sDOF+1]-1
    for (i,p) in enumerate(pini:pend)
      _dof = sDOF_to_dofs.data[p]
      _DOF = _dof_to_DOF(_dof,n_fdofs)
      coeff = sDOF_to_coeffs.data[p]
      DOF_to_DOFs_data[q+i] = _DOF
      DOF_to_coeffs_data[q+i] = coeff
    end
  end

  DOF_to_DOFs = Table(DOF_to_DOFs_data,DOF_to_DOFs_ptrs)
  DOF_to_coeffs = Table(DOF_to_coeffs_data,DOF_to_DOFs_ptrs)

  DOF_to_DOFs, DOF_to_coeffs
end

# Private constructors

function FESpaceWithLinearConstraints(
  fdof_to_dofs::Table,
  fdof_to_coeffs::Table,
  ddof_to_dofs::Table,
  ddof_to_coeffs::Table,
  space::SingleFieldFESpace)

  DOF_to_DOFs, DOF_to_coeffs = _merge_free_and_diri_constraints(
    fdof_to_dofs, fdof_to_coeffs, ddof_to_dofs, ddof_to_coeffs)
  FESpaceWithLinearConstraints!(DOF_to_DOFs,DOF_to_coeffs,space)
end

function _merge_free_and_diri_constraints(fdof_to_dofs, fdof_to_coeffs, ddof_to_dofs, ddof_to_coeffs)
  n_fdofs = length(fdof_to_dofs)
  DOF_to_DOFs = append_tables_globally(fdof_to_dofs,ddof_to_dofs)
  for i in 1:length(DOF_to_DOFs.data)
    dof = DOF_to_DOFs.data[i]
    DOF = _dof_to_DOF(dof,n_fdofs)
    DOF_to_DOFs.data[i] = DOF
  end
  DOF_to_coeffs = Table(vcat(fdof_to_coeffs.data,ddof_to_coeffs.data),DOF_to_DOFs.ptrs)
  DOF_to_DOFs, DOF_to_coeffs
end

function FESpaceWithLinearConstraints(DOF_to_DOFs::Table,DOF_to_coeffs::Table,space::SingleFieldFESpace)
  FESpaceWithLinearConstraints!(copy(DOF_to_DOFs),copy(DOF_to_coeffs),space::SingleFieldFESpace)
end

function FESpaceWithLinearConstraints!(DOF_to_DOFs::Table,DOF_to_coeffs::Table,space::SingleFieldFESpace)

  n_fdofs = num_free_dofs(space)
  mDOF_to_DOF, n_fmdofs = _find_master_dofs(DOF_to_DOFs,n_fdofs)
  DOF_to_mDOFs = _renumber_constraints!(DOF_to_DOFs,mDOF_to_DOF)
  cell_to_ldof_to_dof = Table(get_cell_dof_ids(space))
  cell_to_lmdof_to_mdof = _setup_cell_to_lmdof_to_mdof(cell_to_ldof_to_dof,DOF_to_mDOFs,n_fdofs,n_fmdofs)

  FESpaceWithLinearConstraints(
    space,
    n_fdofs,
    n_fmdofs,
    mDOF_to_DOF,
    DOF_to_mDOFs,
    DOF_to_coeffs,
    cell_to_lmdof_to_mdof,
    cell_to_ldof_to_dof)

end

function _find_master_dofs(DOF_to_DOFs,n_fdofs)
  n_DOFs = length(DOF_to_DOFs)
  DOF_to_ismaster = fill(false,n_DOFs)
  for DOF in 1:n_DOFs
    pini = DOF_to_DOFs.ptrs[DOF]
    pend = DOF_to_DOFs.ptrs[DOF+1]-1
    for p in pini:pend
      _DOF = DOF_to_DOFs.data[p]
      @assert (DOF_to_DOFs.ptrs[_DOF+1]-DOF_to_DOFs.ptrs[_DOF]) == 1 "Rcursive constraints not allowed now"
      @assert DOF_to_DOFs.data[DOF_to_DOFs.ptrs[_DOF]] == _DOF "Rcursive constraints not allowed now"
      DOF_to_ismaster[_DOF] = true
    end
  end
  n_fmdofs = 0
  for DOF in 1:n_fdofs
    if DOF_to_ismaster[DOF]
      n_fmdofs += 1
    end
  end
  mDOF_to_DOF = findall(DOF_to_ismaster)
  mDOF_to_DOF, n_fmdofs
end

function _renumber_constraints!(DOF_to_DOFs,mDOF_to_DOF)
  DOF_to_mDOF = zeros(eltype(DOF_to_DOFs.data),length(DOF_to_DOFs))
  DOF_to_mDOF[mDOF_to_DOF] .= 1:length(mDOF_to_DOF)
  for i in 1:length(DOF_to_DOFs.data)
    DOF = DOF_to_DOFs.data[i]
    mDOF = DOF_to_mDOF[DOF]
    DOF_to_DOFs.data[i] = mDOF
  end
  DOF_to_DOFs
end

function _setup_cell_to_lmdof_to_mdof(cell_to_ldof_to_dof,DOF_to_mDOFs,n_fdofs,n_fmdofs)

  n_cells = length(cell_to_ldof_to_dof)
  cell_to_lmdof_to_mdof_ptrs = zeros(eltype(cell_to_ldof_to_dof.ptrs),n_cells+1)

  for cell in 1:n_cells
    mdofs = Set{Int}()
    pini = cell_to_ldof_to_dof.ptrs[cell]
    pend = cell_to_ldof_to_dof.ptrs[cell+1]-1
    for p in pini:pend
      dof = cell_to_ldof_to_dof.data[p]
      DOF = _dof_to_DOF(dof,n_fdofs)
      qini = DOF_to_mDOFs.ptrs[DOF]
      qend = DOF_to_mDOFs.ptrs[DOF+1]-1
      for q in qini:qend
        mDOF = DOF_to_mDOFs.data[q]
        mdof = _DOF_to_dof(mDOF,n_fmdofs)
        push!(mdofs,mdof)
      end
    end
    cell_to_lmdof_to_mdof_ptrs[cell+1] = length(mdofs)
  end

  length_to_ptrs!(cell_to_lmdof_to_mdof_ptrs)
  ndata = cell_to_lmdof_to_mdof_ptrs[end]-1
  cell_to_lmdof_to_mdof_data = zeros(eltype(cell_to_ldof_to_dof.data),ndata)

  for cell in 1:n_cells
    mdofs = Set{Int}()
    pini = cell_to_ldof_to_dof.ptrs[cell]
    pend = cell_to_ldof_to_dof.ptrs[cell+1]-1
    for p in pini:pend
      dof = cell_to_ldof_to_dof.data[p]
      DOF = _dof_to_DOF(dof,n_fdofs)
      qini = DOF_to_mDOFs.ptrs[DOF]
      qend = DOF_to_mDOFs.ptrs[DOF+1]-1
      for q in qini:qend
        mDOF = DOF_to_mDOFs.data[q]
        mdof = _DOF_to_dof(mDOF,n_fmdofs)
        push!(mdofs,mdof)
      end
    end
    o = cell_to_lmdof_to_mdof_ptrs[cell]-1
    for (lmdof, mdof) in enumerate(mdofs)
      cell_to_lmdof_to_mdof_data[o+lmdof] = mdof
    end
  end

  Table(cell_to_lmdof_to_mdof_data,cell_to_lmdof_to_mdof_ptrs)
end

@inline function _dof_to_DOF(dof,n_fdofs)
  if dof > 0
    DOF = dof
  else
    DOF = n_fdofs - dof
  end
end

@inline function _DOF_to_dof(DOF,n_fdofs)
  if DOF > n_fdofs
    dof = -(DOF-n_fdofs)
  else
    dof = DOF
  end
end

# Implementation of the SingleFieldFESpace interface

function get_cell_dof_ids(f::FESpaceWithLinearConstraints)
  f.cell_to_lmdof_to_mdof
end

function get_cell_dof_basis(f::FESpaceWithLinearConstraints)
  get_cell_dof_basis(f.space)
end

get_dirichlet_dof_ids(f::FESpaceWithLinearConstraints) = Base.OneTo(length(f.mDOF_to_DOF) - f.n_fmdofs)

num_dirichlet_tags(f::FESpaceWithLinearConstraints) = num_dirichlet_tags(f.space)

function get_dirichlet_dof_tag(f::FESpaceWithLinearConstraints)
  ddof_to_tag = get_dirichlet_dof_tag(f.space)
  dmdof_to_tag = zeros(eltype(ddof_to_tag),num_dirichlet_dofs(f))
  _setup_ddof_to_tag!(
    dmdof_to_tag,
    ddof_to_tag,
    f.mDOF_to_DOF,
    f.n_fdofs,
    f.n_fmdofs)
  dmdof_to_tag
end

function _setup_ddof_to_tag!(
  dmdof_to_tag,
  ddof_to_tag,
  mDOF_to_DOF,
  n_fdofs,
  n_fmdofs)

  for mDOF in (n_fmdofs+1):length(mDOF_to_DOF)
    mdof = _DOF_to_dof(mDOF,n_fmdofs)
    @assert mdof < 0 "Dirichlet dofs can only depend on Dirichlet dofs"
    dmdof = -mdof
    DOF = mDOF_to_DOF[mDOF]
    dof = _DOF_to_dof(DOF,n_fdofs)
    @assert dof < 0 "Dirichlet dofs can only depend on Dirichlet dofs"
    ddof = -dof
    dmdof_to_tag[dmdof] = ddof_to_tag[ddof]
  end
end

function get_dirichlet_values(f::FESpaceWithLinearConstraints)
  ddof_to_tag = get_dirichlet_values(f.space)
  dmdof_to_tag = zeros(eltype(ddof_to_tag),num_dirichlet_dofs(f))
  _setup_ddof_to_tag!(
    dmdof_to_tag,
    ddof_to_tag,
    f.mDOF_to_DOF,
    f.n_fdofs,
    f.n_fmdofs)
  dmdof_to_tag
end

function scatter_free_and_dirichlet_values(f::FESpaceWithLinearConstraints,fmdof_to_val,dmdof_to_val)
  fdof_to_val = zero_free_values(f.space)
  ddof_to_val = zero_dirichlet_values(f.space)
  _setup_dof_to_val!(
    fdof_to_val,
    ddof_to_val,
    fmdof_to_val,
    dmdof_to_val,
    f.DOF_to_mDOFs,
    f.DOF_to_coeffs,
    f.n_fdofs,
    f.n_fmdofs)
  scatter_free_and_dirichlet_values(f.space,fdof_to_val,ddof_to_val)
end

function _setup_dof_to_val!(
  fdof_to_val,
  ddof_to_val,
  fmdof_to_val,
  dmdof_to_val,
  DOF_to_mDOFs,
  DOF_to_coeffs,
  n_fdofs,
  n_fmdofs)

  T = eltype(fdof_to_val)

  for DOF in 1:length(DOF_to_mDOFs)
    pini = DOF_to_mDOFs.ptrs[DOF]
    pend = DOF_to_mDOFs.ptrs[DOF+1]-1
    val = zero(T)
    for p in pini:pend
      mDOF = DOF_to_mDOFs.data[p]
      coeff = DOF_to_coeffs.data[p]
      mdof = _DOF_to_dof(mDOF,n_fmdofs)
      if mdof > 0
        fmdof = mdof
        val += fmdof_to_val[fmdof]*coeff
      else
        dmdof = -mdof
        val += dmdof_to_val[dmdof]*coeff
      end
    end
    dof = _DOF_to_dof(DOF,n_fdofs)
    if dof > 0
      fdof = dof
      fdof_to_val[fdof] = val
    else
      ddof = -dof
      ddof_to_val[ddof] = val
    end
  end

end

function gather_free_and_dirichlet_values(f::FESpaceWithLinearConstraints,cell_to_ludof_to_val)
  fdof_to_val, ddof_to_val = gather_free_and_dirichlet_values(f.space,cell_to_ludof_to_val)
  fmdof_to_val = zero_free_values(f)
  dmdof_to_val = zero_dirichlet_values(f)
  _setup_mdof_to_val!(
    fmdof_to_val,
    dmdof_to_val,
    fdof_to_val,
    ddof_to_val,
    f.mDOF_to_DOF,
    f.n_fdofs,
    f.n_fmdofs)
  fmdof_to_val, dmdof_to_val
end

function gather_free_and_dirichlet_values!(fmdof_to_val,dmdof_to_val,f::FESpaceWithLinearConstraints,cell_to_ludof_to_val)
  fdof_to_val, ddof_to_val = gather_free_and_dirichlet_values(f.space,cell_to_ludof_to_val)
  _setup_mdof_to_val!(
    fmdof_to_val,
    dmdof_to_val,
    fdof_to_val,
    ddof_to_val,
    f.mDOF_to_DOF,
    f.n_fdofs,
    f.n_fmdofs)
  fmdof_to_val, dmdof_to_val
end

function _setup_mdof_to_val!(
  fmdof_to_val,
  dmdof_to_val,
  fdof_to_val,
  ddof_to_val,
  mDOF_to_DOF,
  n_fdofs,
  n_fmdofs)

  for mDOF in 1:length(mDOF_to_DOF)
    DOF = mDOF_to_DOF[mDOF]
    dof = _DOF_to_dof(DOF,n_fdofs)
    if dof > 0
      fdof = dof
      val = fdof_to_val[fdof]
    else
      ddof = -dof
      val = ddof_to_val[ddof]
    end
    mdof = _DOF_to_dof(mDOF,n_fmdofs)
    if mdof > 0
      fmdof = mdof
      fmdof_to_val[fmdof] = val
    else
      dmdof = -mdof
      dmdof_to_val[dmdof] = val
    end
  end

end

# Implementation of FESpace interface

function get_triangulation(f::FESpaceWithLinearConstraints)
  get_triangulation(f.space)
end

get_free_dof_ids(f::FESpaceWithLinearConstraints) = Base.OneTo(f.n_fmdofs)

function get_vector_type(f::FESpaceWithLinearConstraints)
  get_vector_type(f.space)
end

function get_cell_shapefuns(f::FESpaceWithLinearConstraints)
  get_cell_shapefuns(f.space)
end

function get_cell_shapefuns_trial(f::FESpaceWithLinearConstraints)
  get_cell_shapefuns_trial(f.space)
end

function CellField(f::FESpaceWithLinearConstraints,cellvals)
  CellField(f.space,cellvals)
end

ConstraintStyle(::Type{<:FESpaceWithLinearConstraints}) = Constrained()

function get_cell_isconstrained(f::FESpaceWithLinearConstraints)
  #TODO this can be heavily optimized
  n = length(get_cell_dof_ids(f))
  Fill(true,n)
end

function get_cell_constraints(f::FESpaceWithLinearConstraints)

  k = LinearConstraintsMap(
    f.DOF_to_mDOFs,
    f.DOF_to_coeffs,
    length(f.mDOF_to_DOF),
    f.n_fmdofs,
    f.n_fdofs)

  cell_to_mat = get_cell_constraints(f.space)
  lazy_map(k,f.cell_to_lmdof_to_mdof,f.cell_to_ldof_to_dof,cell_to_mat)

end

struct LinearConstraintsMap{A,B} <: Map
  DOF_to_mDOFs::A
  DOF_to_coeffs::B
  n_mDOFs::Int
  n_fmdofs::Int
  n_fdofs::Int
end

function return_cache(k::LinearConstraintsMap,lmdof_to_mdof,ldof_to_dof,mat)
  n_lmdofs = length(lmdof_to_mdof)
  n_ldofs = length(ldof_to_dof)
  n_ludofs = size(mat,2)
  @assert n_ldofs == size(mat,1)
  m1 = CachedArray(zeros(n_lmdofs,n_ldofs))
  m2 = CachedArray(zeros(n_lmdofs,n_ludofs))
  mDOF_to_lmdof = zeros(Int8,k.n_mDOFs)
  m1, m2, mDOF_to_lmdof
end

function evaluate!(cache,k::LinearConstraintsMap,lmdof_to_mdof,ldof_to_dof,mat)
  m1, m2, mDOF_to_lmdof = cache
  n_lmdofs = length(lmdof_to_mdof)
  n_ldofs = length(ldof_to_dof)
  n_ludofs = size(mat,2)

  setsize!(m1,(n_lmdofs,n_ldofs))
  setsize!(m2,(n_lmdofs,n_ludofs))
  a1 = m1.array
  a2 = m2.array
  fill!(a1,zero(eltype(a1)))

  for (lmdof,mdof) in enumerate(lmdof_to_mdof)
    mDOF = _dof_to_DOF(mdof,k.n_fmdofs)
    mDOF_to_lmdof[mDOF] = lmdof
  end

  for (ldof,dof) in enumerate(ldof_to_dof)
    DOF = _dof_to_DOF(dof,k.n_fdofs)
    qini = k.DOF_to_mDOFs.ptrs[DOF]
    qend = k.DOF_to_mDOFs.ptrs[DOF+1]-1
    for q in qini:qend
      mDOF = k.DOF_to_mDOFs.data[q]
      coeff = k.DOF_to_coeffs.data[q]
      lmdof = mDOF_to_lmdof[mDOF]
      a1[lmdof,ldof] = coeff
    end
  end

  #TODO this is not always needed
  mul!(a2,a1,mat)
  a2
end


text:
using Documenter, AWSTools

makedocs(;
    modules=[AWSTools],
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", nothing) == "true", assets=["assets/invenia.css"]
    ),
    pages=["Home" => "index.md"],
    repo="https://github.com/JuliaCloud/AWSTools.jl/blob/{commit}{path}#L{line}",
    sitename="AWSTools.jl",
    authors="Invenia Technical Computing",
    checkdocs=:exports,
    strict=true,
)


text:
# # A simple QR decomposition

# ## Functions used in this example

using NiLang, NiLang.AD, Test

# ## The QR decomposition
# Let us consider a naive implementation of QR decomposition from scratch.
# This implementation is just a proof of principle which does not consider reorthogonalization and other practical issues.

@i function qr(Q, R, A::Matrix{T}) where T
    @routine begin
        anc_norm ← zero(T)
        anc_dot ← zeros(T, size(A,2))
        ri ← zeros(T, size(A,1))
    end
    for col = 1:size(A, 1)
        ri .+= A[:,col]
        for precol = 1:col-1
            i_dot(anc_dot[precol], Q[:,precol], ri)
            R[precol,col] += anc_dot[precol]
            for row = 1:size(Q,1)
                ri[row] -=
                    anc_dot[precol] * Q[row, precol]
            end
        end
        i_norm2(anc_norm, ri)

        R[col, col] += anc_norm^0.5
        for row = 1:size(Q,1)
            Q[row,col] += ri[row] / R[col, col]
        end

        ~begin
            ri .+= A[:,col]
            for precol = 1:col-1
                i_dot(anc_dot[precol], Q[:,precol], ri)
                for row = 1:size(Q,1)
                    ri[row] -= anc_dot[precol] *
                        Q[row, precol]
                end
            end
            i_norm2(anc_norm, ri)
        end
    end
    ~@routine
end

# Here, in order to avoid frequent uncomputing, we allocate ancillas `ri` and `anc_dot` as vectors.
# The expression in `~` is used to uncompute `ri`, `anc_dot` and `anc_norm`.
# `i_dot` and `i_norm2` are reversible functions to compute dot product and vector norm.
# One can quickly check the correctness of the gradient function

A  = randn(4,4)
q, r = zero(A), zero(A)
@i function test1(out, q, r, A)
    qr(q, r, A)
    i_sum(out, q)
end

@test check_grad(test1, (0.0, q, r, A); iloss=1)

# Here, the loss function `test1` is defined as the sum of the output unitary matrix `q`.
# The `check_grad` function is a gradient checker function defined in module `NiLang.AD`.


text:
using FakePackage, ReTest

retest(FakePackage, recursive=false)
@test FakePackage.RUN == [1, 2]

FakePackage.runtests(r"begin-end")
@test FakePackage.RUN == [1, 2, 1]
FakePackage.runtests(r" for")
@test FakePackage.RUN == [1, 2, 1, 2]

retest(dry=true) # only to update InlineTest.TESTED_MODULES
@test FakePackage.Sub1 in ReTest.TESTED_MODULES
@test FakePackage.Sub2.SubSub in ReTest.TESTED_MODULES

retest(FakePackage, recursive=true)
@test FakePackage.RUN == [1, 2, 1, 2, 1, 2, 3, 4]
retest(FakePackage) # equivalent to recursive=true
@test FakePackage.RUN == [1, 2, 1, 2, 1, 2, 3, 4, 1, 2, 3, 4]

module Included
using ReTest
RUN = []

@testset "include parent" begin
    include("included.jl")
    include(joinpath(@__DIR__, "included.jl"))
    if VERSION >= v"1.5"
        include(identity, "included.jl")
        include(identity, joinpath(@__DIR__, "included.jl"))
    end
end
end

retest(Included)
@test Included.RUN == (VERSION >= v"1.5" ? [0, 0, 0, 0] : [0, 0])

include("../../setup.jl")

using ReTest: process_args, and

ReTest.Test.@testset "retest: load" begin
    @test process_args(()).modules == [
        FakePackage => and(),
        FakePackage.Sub1 => and(),
        FakePackage.Sub2.SubSub => and(),
        Main.Included => and(),
    ]
    check([], dry=true, verbose=0, output = """
FakePackage
FakePackage.Sub1
FakePackage.Sub2.SubSub
Main.Included
""")

    check([], dry=true, verbose=0, load=true, output = """
FakePackage
FakePackage.Sub1
FakePackage.Sub2.SubSub
Main.Included
Main.FakePackageTests
Main.FakePackageTests.Sub
Main.FakePackageTests.Sub.SubSub
""")

    check(FakePackage, [], dry=true, verbose=0, load=true, recursive=false, output = """
FakePackage
Main.FakePackageTests
""")

    @test process_args((FakePackage, 2)).modules ==
        [
            FakePackage => and(2),
            FakePackage.Sub1 => and(2),
            FakePackage.Sub2.SubSub => and(2),
        ]
    @test process_args((FakePackage, 2), load=true, recursive=false).modules ==
        [
            FakePackage => and(2),
            FakePackageTests => and(2),
        ]
    # cumulative patterns
    @test process_args((FakePackage, 2, FakePackageTests => 3),
                       load=true, recursive=false).modules ==
        [
            FakePackage => and(2),
            FakePackageTests => and(3, 2),
        ]
    # specified pattern for test module
    @test process_args((FakePackage => 2, FakePackageTests => 3),
                       load=true, recursive=false).modules ==
        [
            FakePackage => and(2),
            FakePackageTests => and(3),
        ]
    @test process_args((FakePackage => 2, FakePackageTests => 3),
                       load=true, recursive=true).modules ==
        [
            FakePackage => and(2),
            FakePackageTests => and(3),
            FakePackage.Sub1 => and(2),
            FakePackage.Sub2.SubSub => and(2),
            FakePackageTests.Sub => and(3),
            FakePackageTests.Sub.SubSub => and(3),
        ]
end


text:
using Test
using StaticArrays
using ViewerGL
GL = ViewerGL

@testset "Frustum.jl" begin

   # projectPoint(map::FrustumMap,p3::Point3d)
   mesh = GL.GLCuboid(GL.Box3d(GL.Point3d(0,0,0),GL.Point3d(1,1,1)))
   GL.VIEW([mesh],false);
   viewport=[0,0,GL.viewer.W,GL.viewer.H];
   projection =GL.getProjection(GL.viewer);
   modelview=GL.getModelview(GL.viewer);
   themap=GL.FrustumMap(viewport,projection,modelview)
#=   @testset "FrustumMap" begin

      @test map=GL.FrustumMap(viewport,projection,modelview) ==
      ViewerGL.FrustumMap(
      [1024.0 0.0 0.0 1024.0; 0.0 768.0 0.0 768.0; 0.0 0.0 0.5 0.5; 0.0 0.0 0.0 1.0],
      [1.29904 0.0 0.0 0.0; 0.0 1.73205 0.0 0.0; 0.0 0.0 -1.00401 -0.0801603; 0.0 0.0 -1.0 0.0],
      [-0.707107 0.707107 0.0 0.0; -0.408248 -0.408248 0.816497 0.0; 0.57735 0.57735 0.57735 -5.19615; 0.0 0.0 0.0 1.0],
      [0.000976563 0.0 0.0 -1.0; 0.0 0.00130208 0.0 -1.0; 0.0 0.0 2.0 -1.0; 0.0 0.0 0.0 1.0],
      [0.7698 -0.0 -0.0 -0.0; -0.0 0.57735 0.0 -0.0; -0.0 -0.0 -0.0 -1.0; -0.0 -0.0 -12.475 12.525],
      [-0.707107 -0.408248 0.57735 3.0; 0.707107 -0.408248 0.57735 3.0; 0.0 0.816497 0.57735 3.0; 0.0 0.0 -0.0 1.0])

      @test map.viewport==GL.Matrix4([1024.0 0.0 0.0 1024.0; 0.0 768.0 0.0 768.0; 0.0 0.0 0.5 0.5; 0.0 0.0 0.0 1.0])
      @test map.projection==GL.Matrix4([1.29904 0.0 0.0 0.0; 0.0 1.73205 0.0 0.0; 0.0 0.0 -1.00401 -0.0801603; 0.0 0.0 -1.0 0.0])
      @test map.modelview==GL.Matrix4([-0.707107 0.707107 0.0 0.0; -0.408248 -0.408248 0.816497 0.0; 0.57735 0.57735 0.57735 -5.19615; 0.0 0.0 0.0 1.0])
      @test map.inv_viewport==GL.Matrix4([0.000976563 0.0 0.0 -1.0; 0.0 0.00130208 0.0 -1.0; 0.0 0.0 2.0 -1.0; 0.0 0.0 0.0 1.0])
      @test map.inv_projection==GL.Matrix4([0.7698 -0.0 -0.0 -0.0; -0.0 0.57735 0.0 -0.0; -0.0 -0.0 -0.0 -1.0; -0.0 -0.0 -12.475 12.525])
      @test map.inv_modelview=GL.Matrix4([-0.707107 -0.408248 0.57735 3.0; 0.707107 -0.408248 0.57735 3.0; 0.0 0.816497 0.57735 3.0; 0.0 0.0 -0.0 1.0])
   end

   # function FrustumMap(viewport,projection::Matrix4,modelview::Matrix4)
   @testset "FrustumMap" begin
      @test
      @test
      @test
      @test
   end
=#
   # function projectPoint(map::FrustumMap,p3::Point3d)
   @testset "projectPoint" begin
      @test GL.projectPoint(themap::GL.FrustumMap,GL.Point3d(1.5,1.5,1.5))==GL.Point3d(512.0, 384.0, 0.9865771471158417)
      @test GL.projectPoint(themap::GL.FrustumMap,GL.Point3d(-1.5,-1.5,-1.5))==GL.Point3d(512.0, 384.0,0.9968617210493019)
      @test typeof(GL.projectPoint(themap::GL.FrustumMap,GL.Point3d(1.5,1.5,1.5)))==StaticArrays.MArray{Tuple{3},Float64,1,3}
      @test typeof(themap)==ViewerGL.FrustumMap
      @test typeof(GL.Point3d(1.5,1.5,1.5))==StaticArrays.MArray{Tuple{3},Float64,1,3}
    end

   # function unprojectPoint(map::FrustumMap,x::Float32,y::Float32, z::Float32)
   @testset "unprojectPoint" begin
      x=1.5; y=1.5; z=1.5;
      p4 = (themap.inv_modelview * (themap.inv_projection * (themap.inv_viewport * GL.Point4d(x,y,z, 1.0))))
      @test typeof(themap)==ViewerGL.FrustumMap
      @test GL.Point4d(x,y,z, 1.0)==[1.5,1.5,1.5,1.0]
      @test p4 == [-37.0748323866816,-38.160305054815986,-38.32191336607133,-12.425000000000002]
   end

end


text:
# Copyright 2020 Neven Sajko. All rights reserved. See LICENSE for a license.

module FitPlotDat

using DataFrames, Gadfly, Polynomials

export main

const F = BigFloat
const P = Polynomial{F}
const V = Vector{F}
const VI = Vector{Int}

const fpprec = 8192

function readTable()::Tuple{V, V, VI}
	local r::Tuple{V, V, VI} = (F[], F[], Int[])
	while !eof(stdin)
		local latest_seed::Int = parse(Int, readline(), base = 16)
		local distinct_value_count::Int = parse(Int, readline(), base = 10)
		push!(r[1], BigFloat(latest_seed, precision=fpprec))
		push!(r[2], BigFloat(distinct_value_count, precision=fpprec))
		push!(r[3], latest_seed - distinct_value_count + 1)
	end
	r
end

function main()
	setprecision(BigFloat, fpprec)

	local x::V, y::V, d::VI = readTable()

	local line::P = fit(x, y, 1)
	local poly::P = fit(x, y, 47)

	draw(SVG("data/diffs.svg", 50cm, 50cm), plot(
	   Coord.cartesian(xmin=0, ymin=0, ymax=4),
	   layer(DataFrame(latest_seed = x, diff = d), x=:latest_seed, y=:diff)))

	draw(SVG(stdout, 50cm, 50cm), plot(
	   Coord.cartesian(xmin=0, xmax=Int(3e9), ymin=0, ymax=Int(3e9), fixed=true),
	   layer(DataFrame(latest_seed = x, distinct_value_count = y), x=:latest_seed, y=:distinct_value_count),
	   layer([x::F -> line(x), x::F -> poly(x)], 1000, x[length(x)] + 1000, color=[string("linear approximation f(x) = ", line), "approx. with a polynomial of degree 47"])))

	return nothing
end

end

FitPlotDat.main()

max_stars_repo_name: nsajko/go-prng-test

max_issues_repo_name: nsajko/go-prng-test

max_forks_repo_name: nsajko/go-prng-test


text:
using Test
using LinearAlgebra
import TensorTrains:one_body_to_matrix, one_body_mpo, two_body_to_matrix, two_body_mpo
using Random

@testset "MPO constructors" begin
    L=6
    p=rand(1:L)
    q=rand(1:L)
    #one-body MPO constructors
    H = one_body_to_matrix(p,q,L)
    Hmpo = one_body_mpo(p,q,L)
    H2 = tto_to_tensor(Hmpo)
    @test isapprox(norm(H-reshape(H2,2^L,2^L)),0.0,atol=1e-10)
    #two-body MPO constructors
#    k,l,m,n=1,2,3,4
#    H = two_body_to_matrix(k,l,m,n,L) #fails but probably bug in two_body_to_matrix
#    Hmpo = two_body_mpo(k,l,m,n,L)
#    H2 = tto_to_tensor(Hmpo)
#    @test isapprox(norm(H-reshape(H2,2^L,2^L)),0.0,atol=1e-10)
end

@testset "PPP_C_NH_N" begin
    N = 6
    H_tto = PPP_C_NH_N(N,order=randperm(12))
    ψ_0 = tt_up_rks(half_filling(N),32;ϵ_wn=1e-2)
    ψ_0 = orthogonalize(ψ_0)
    E_mals, ψ_tt = mals_eigsolv(H_tto,ψ_0)
    @test isapprox(E_mals[end],-0.46752535023652453,atol=1e-6) #value in Bendazzoli, G. L., & Evangelisti, S. (1991). Full-CI calculations of alternant cyclic polyenes (CH) N, N= 2, 4, 6, ƒ 18, in the PPP approximation. Chemical physics letters, 185(1-2), 125-130.
end

max_stars_repo_name: msdupuy/Tensor-Train-Julia

max_issues_repo_name: msdupuy/Tensor-Train-Julia

max_forks_repo_name: msdupuy/Tensor-Train-Julia


text:
mutable struct ScalarKalmanState{Fl<:AbstractFloat}
    v::Fl
    F::Fl
    att::Fl
    a::Fl
    Ptt::Fl
    P::Fl
    llk::Fl
    steady_state::Bool
    # Auxiliary matrices
    P_to_check_steady_state::Fl
    function ScalarKalmanState(a1::Fl, P1::Fl) where Fl
        P_to_check_steady_state = zero(Fl)
        return new{Fl}(
            zero(Fl),
            zero(Fl),
            zero(Fl),
            a1,
            zero(Fl),
            P1,
            zero(Fl),
            false,
            P_to_check_steady_state,
        )
    end
end

function save_a1_P1_in_filter_output!(
    filter_output::FilterOutput{Fl}, kalman_state::ScalarKalmanState{Fl}
) where Fl
    filter_output.a[1] = copy(fill(kalman_state.a, 1))
    filter_output.P[1] = copy(fill(kalman_state.P, 1, 1))
    return filter_output
end

function save_kalman_state_in_filter_output!(
    filter_output::FilterOutput{Fl}, kalman_state::ScalarKalmanState{Fl}, t::Int
) where Fl
    filter_output.v[t] = copy(fill(kalman_state.v, 1))
    filter_output.F[t] = copy(fill(kalman_state.F, 1, 1))
    filter_output.a[t + 1] = copy(fill(kalman_state.a, 1))
    filter_output.att[t] = copy(fill(kalman_state.att, 1))
    filter_output.P[t + 1] = copy(fill(kalman_state.P, 1, 1))
    filter_output.Ptt[t] = copy(fill(kalman_state.Ptt, 1, 1))
    # There is no diffuse part in UnivariateKalmanState
    filter_output.Pinf[t] = copy(fill(zero(Fl), size(filter_output.Ptt[t])))
    return filter_output
end

"""
    ScalarKalmanFilter{Fl <: Real} <: KalmanFilter

Similar to the univariate Kalman filter but exploits the fact that the dimension of the
state is equal to 1.
"""
mutable struct ScalarKalmanFilter{Fl<:Real} <: KalmanFilter
    steadystate_tol::Fl
    a1::Fl
    P1::Fl
    skip_llk_instants::Int
    kalman_state::ScalarKalmanState

    function ScalarKalmanFilter(
        a1::Fl, P1::Fl, skip_llk_instants::Int=1, steadystate_tol::Fl=Fl(1e-5)
    ) where Fl
        kalman_state = ScalarKalmanState(a1, P1)
        return new{Fl}(steadystate_tol, a1, P1, skip_llk_instants, kalman_state)
    end
end

function reset_filter!(kf::ScalarKalmanFilter{Fl}) where Fl
    reset_filter!(kf.kalman_state, kf.a1, kf.P1)
    return kf
end

function reset_filter!(kalman_state::ScalarKalmanState, a1::Fl, P1::Fl) where Fl
    kalman_state.a = a1
    kalman_state.P = P1
    kalman_state.llk = zero(Fl)
    kalman_state.P_to_check_steady_state = zero(Fl)
    kalman_state.steady_state = false
    return kalman_state
end

function scalar_check_steady_state(kalman_state::ScalarKalmanState{Fl}, tol::Fl) where Fl
    if abs((kalman_state.P - kalman_state.P_to_check_steady_state) / kalman_state.P) > tol
        kalman_state.P_to_check_steady_state = kalman_state.P
        return kalman_state
    end
    kalman_state.steady_state = true
    return kalman_state
end

function scalar_update_v!(
    kalman_state::ScalarKalmanState{Fl}, y::Fl, Z::Fl, d::Fl
) where Fl
    kalman_state.v = y - Z * kalman_state.a - d
    return kalman_state
end

function scalar_update_F!(kalman_state::ScalarKalmanState{Fl}, Z::Fl, H::Fl) where Fl
    kalman_state.F = Z * kalman_state.P * Z + H
    return kalman_state
end

function scalar_update_att!(kalman_state::ScalarKalmanState{Fl}, Z::Fl) where Fl
    kalman_state.att = kalman_state.a + Z * kalman_state.P * kalman_state.v / kalman_state.F
    return kalman_state
end

function scalar_update_a!(kalman_state::ScalarKalmanState{Fl}, T::Fl, c::Fl) where Fl
    kalman_state.a = T * kalman_state.att + c
    return kalman_state
end

function scalar_update_Ptt!(kalman_state::ScalarKalmanState{Fl}, Z::Fl) where Fl
    kalman_state.Ptt =
        kalman_state.P - kalman_state.P * Z * Z * kalman_state.P / kalman_state.F
    return kalman_state
end

function scalar_update_P!(kalman_state::ScalarKalmanState{Fl}, T::Fl, RQR::Fl) where Fl
    kalman_state.P = T * kalman_state.Ptt * T + RQR
    return kalman_state
end

function update_llk!(kalman_state::ScalarKalmanState{Fl}) where Fl
    kalman_state.llk -= (
        HALF_LOG_2_PI + 0.5 * (log(kalman_state.F) + kalman_state.v^2 / kalman_state.F)
    )
    return kalman_state
end

function scalar_update_kalman_state!(
    kalman_state, y, Z, T, H, RQR, d, c, skip_llk_instants, tol, t
)
    if isnan(y)
        kalman_state.v = NaN
        scalar_update_F!(kalman_state, Z, H)
        kalman_state.att = kalman_state.a
        scalar_update_a!(kalman_state, T, c)
        kalman_state.Ptt = kalman_state.P
        scalar_update_P!(kalman_state, T, RQR)
        kalman_state.steady_state = false # Not on steadystate anymore
    elseif kalman_state.steady_state
        scalar_update_v!(kalman_state, y, Z, d)
        scalar_update_att!(kalman_state, Z)
        scalar_update_a!(kalman_state, T, c)
        if t > skip_llk_instants
            update_llk!(kalman_state)
        end
    else
        scalar_update_v!(kalman_state, y, Z, d)
        scalar_update_F!(kalman_state, Z, H)
        scalar_update_att!(kalman_state, Z)
        scalar_update_a!(kalman_state, T, c)
        scalar_update_Ptt!(kalman_state, Z)
        scalar_update_P!(kalman_state, T, RQR)
        scalar_check_steady_state(kalman_state, tol)
        if t > skip_llk_instants
            update_llk!(kalman_state)
        end
    end
    return kalman_state
end

function optim_kalman_filter(
    sys::LinearUnivariateTimeInvariant{Fl}, filter::ScalarKalmanFilter{Fl}
) where Fl
    return filter_recursions!(
        filter.kalman_state, sys, filter.steadystate_tol, filter.skip_llk_instants
    )
end

function kalman_filter!(
    filter_output::FilterOutput, sys::StateSpaceSystem, filter::ScalarKalmanFilter{Fl}
) where Fl
    filter_recursions!(
        filter_output,
        filter.kalman_state,
        sys,
        filter.steadystate_tol,
        filter.skip_llk_instants,
    )
    return filter_output
end

function filter_recursions!(
    kalman_state::ScalarKalmanState{Fl},
    sys::LinearUnivariateTimeInvariant{Fl},
    steadystate_tol::Fl,
    skip_llk_instants::Int,
) where Fl
    RQR = sys.R[1] * sys.Q[1] * sys.R[1]
    @inbounds for t in eachindex(sys.y)
        scalar_update_kalman_state!(
            kalman_state,
            sys.y[t],
            sys.Z[1],
            sys.T[1],
            sys.H,
            RQR,
            sys.d,
            sys.c[1],
            skip_llk_instants,
            steadystate_tol,
            t,
        )
    end
    return kalman_state.llk
end

function filter_recursions!(
    filter_output::FilterOutput,
    kalman_state::ScalarKalmanState{Fl},
    sys::LinearUnivariateTimeInvariant,
    steadystate_tol::Fl,
    skip_llk_instants::Int,
) where Fl
    RQR = sys.R[1] * sys.Q[1] * sys.R[1]
    save_a1_P1_in_filter_output!(filter_output, kalman_state)
    @inbounds for t in eachindex(sys.y)
        scalar_update_kalman_state!(
            kalman_state,
            sys.y[t],
            sys.Z[1],
            sys.T[1],
            sys.H,
            RQR,
            sys.d,
            sys.c[1],
            skip_llk_instants,
            steadystate_tol,
            t,
        )
        save_kalman_state_in_filter_output!(filter_output, kalman_state, t)
    end
    return filter_output
end


text:
module JuMPComplementarity

using MathProgBase, AmplNLWriter
importall MathProgBase.SolverInterface

export AmplComplementaritySolver

struct AmplComplementaritySolver <: MathProgBase.AbstractMathProgSolver
    solver_command::String
    options::Vector{String}
    filename::String
end

import AmplNLWriter.NonlinearModel
import AmplNLWriter.loadproblem!

macro complementarity(m::Model, var::Expr, g::Expr)

#=
function printModel(m::)

end
=#

end # module


text:
#
# Multiplexer.
#

function IfThenElseQuery(qcond::Query, qthen::Query, qelse::Query)
    dom = ibound(domain(input(qthen)), domain(input(qelse)))
    q = RecordQuery(ItQuery(dom), LiftQuery(c -> c ? 1 : 2, qcond))
    q >>
    Query(
        SwitchSig(),
        [qthen, qelse],
        Input(domain(output(q))),
        Output(
            obound(domain(output(qthen)), domain(output(qelse))),
            obound(mode(output(qcond)), mode(output(qthen)), mode(output(qelse)))))
end

immutable SwitchSig <: AbstractSignature
end

function ev(::SwitchSig, args::Vector{Query}, ity::Input, oty::Output, iflow::InputFlow)
    ivals = values(column(values(iflow), 1))
    imux = column(values(iflow), 2)
    ocols = Vector{Column}(length(args))
    coloty = Output(domain(fields(domain(ity))[1]), OutputMode(true, isplural(imux)))
    for k = 1:length(args)
        arg = args[k]
        ocol = switch_split_impl(ivals, imux, k)
        colflow = OutputFlow(coloty, ocol)
        icolflow = distribute(narrow(iflow, input(arg)), colflow)
        ocols[k] = column(ev(arg, icolflow))
    end
    col = switch_merge_impl(datatype(domain(oty)), imux, ocols)
    return OutputFlow(oty, col)
end

function switch_split_impl(ivals::AbstractVector, mux::Column, k::Int)
    size = 0
    for m in values(mux)
        if m == k
            size += 1
        end
    end
    offs = Vector{Int}(length(mux)+1)
    offs[1] = 1
    idxs = Vector{Int}(size)
    cr = cursor(mux)
    n = 1
    while !done(mux, cr)
        next!(mux, cr)
        for m in cr
            if m == k
                idxs[n] = cr.pos
                n += 1
            end
        end
        offs[cr.pos+1] = n
    end
    vals = ivals[idxs]
    return Column{true,isplural(mux)}(offs, vals)
end

function switch_merge_impl(T::DataType, mux::Column, cols::Vector{Column})
    optional = isoptional(mux)
    plural = isplural(mux)
    size = 0
    for col in cols
        optional |= isoptional(col)
        plural |= isplural(col)
        size += length(values(col))
    end
    offs = Vector{Int}(length(mux)+1)
    offs[1] = 1
    vals = Vector{T}(size)
    n = 1
    crs = ColumnCursor[]
    for col in cols
        push!(crs, cursor(col))
    end
    muxcr = cursor(mux)
    while !done(mux, muxcr)
        next!(mux, muxcr)
        for m in muxcr
            col = cols[m]
            cr = crs[m]
            next!(col, cr)
            copy!(vals, n, cr)
            n += length(cr)
        end
        offs[muxcr.pos+1] = n
    end
    return Column{optional,plural}(offs, vals)
end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
using MRMP
import Random: seed!

solvers =
    ["MRMP.PP", "MRMP.CBS", "MRMP.RRT", "MRMP.RRT_connect", "MRMP.PRM_direct", "MRMP.SSSP"]

TIME_LIMIT = 1.0

ins_commons_params = Dict(:N => 3, :rad => 0.025, :num_obs => 3, :rad_obs => 0.05)

function test_all_solvers(ins_gen_fn)
    seed!(1)
    config_init, config_goal, obstacles, ins_params... = ins_gen_fn(; ins_commons_params...)
    connect = gen_connect(config_init[1], obstacles, ins_params...)
    collide = gen_collide(config_init[1], ins_params...)
    check_goal = MRMP.gen_check_goal(config_goal)

    for solver in solvers
        @info solver
        seed!(1)
        @time eval(Meta.parse(solver))(
            config_init,
            config_goal,
            connect,
            collide,
            check_goal;
            TIME_LIMIT = TIME_LIMIT,
            VERBOSE = 1,
        )
    end
end

@testset "point2d" begin
    test_all_solvers(MRMP.gen_random_instance_StatePoint2D)
end
@testset "point3d" begin
    test_all_solvers(MRMP.gen_random_instance_StatePoint3D)
end
@testset "line2d" begin
    test_all_solvers(MRMP.gen_random_instance_StateLine2D)
end
@testset "arm22" begin
    test_all_solvers(MRMP.gen_random_instance_StateArm22)
end
@testset "arm33" begin
    test_all_solvers(MRMP.gen_random_instance_StateArm33)
end
@testset "snake2d" begin
    test_all_solvers(MRMP.gen_random_instance_StateSnake2D)
end
@testset "dubins2d" begin
    test_all_solvers(MRMP.gen_random_instance_StateDubins)
end
@testset "capsule3d" begin
    test_all_solvers(MRMP.gen_random_instance_StateDubins)
end


text:
# ## Data

import CSV
import Random
import TuringModels

using DataFrames

Random.seed!(1)

data_path = joinpath(TuringModels.project_root, "data", "UCBadmit.csv")
df = CSV.read(data_path, DataFrame; delim=';')

# ## Model

using StatsFuns: logistic
using Turing

@model function m11_5(admit, applications)
  θ ~ truncated(Exponential(1), 0, Inf)
  α ~ Normal(0, 2)

  ## alpha and beta for the BetaBinomial must be provided.
  ## The two parameterizations are related by
  ## alpha = prob * theta, and beta = (1-prob) * theta.
  ## See https://github.com/rmcelreath/rethinking/blob/master/man/dbetabinom.Rd

  prob = logistic(α)
  alpha = prob * θ
  beta = (1 - prob) * θ
  admit .~ BetaBinomial.(applications, alpha, beta)
end

model = m11_5(df.admit, df.applications);

# ## Output

chns = sample(model, NUTS(), 1000)

# \defaultoutput{}

# ## Original output

"""
        mean   sd  5.5% 94.5% n_eff Rhat
theta   2.74 0.96  1.43  4.37  3583    1
a       -0.37 0.31 -0.87  0.12  3210    1
""";


text:
function hc4(X::IntervalBox{N,T}, constraints::Vector{Constraint{T}}, tol=1e-5) where{N, T}
    n = length(constraints)
    while true
        X_temp = X
        for i in 1:n
            X = invokelatest(constraints[i].C,  constraints[i].bound, X)
        end
        if isempty(X) || sum(dist.(X, X_temp)) < tol
            break
        end
    end

    #new_constraints = Constraint{T}[]

    #for i in 1:n
    #    if !(invokelatest(constraints[i].C, X) ⊆ constraints[i].bound)
    #        push!(new_constraints, constraints[i])
    #    end
    #end
    #return new_constraints, X
    return constraints, X
end

function contraction(f::Function, C, global_min::Float64, X::IntervalBox{N,T}, constraints::Vector{Constraint{T}}, tol=1e-5) where {N, T}
    lb = -Inf
    while true
        X_temp = X
        X = invokelatest(C, -Inf..global_min, X)
        lb = inf(f(X))
        constraints, X = hc4(X, constraints)

        if isempty(X) || sum(dist.(X, X_temp)) < tol
            break
        end
    end

    return lb, X, constraints
end

function generate_random_feasible_point(X::IntervalBox{N, T}, constraints::Vector{Constraint{T}}) where{N, T}
    for i in 1:30
        point = [X[j].lo + (1-rand())*(X[j].hi - X[j].lo) for j in 1:length(X)]      # discover a random point in interval box X
        for j in 1:length(constraints)
            if !(invokelatest(constraints[j].C, point) ⊆ constraints[j].bound)
                break
            end
            if j == length(constraints)
                return (true, point)    # return a feasible point
            end
        end
        if i == 30
            return (false, point)    # returns a infeasible point
        end
    end

end

function convex_hull(vec::Vector{IntervalBox{N,T}}) where{N, T}

    x_big = Interval{T}[]
    num_variables = N
    num_boxes = length(vec)

    for i in 1:num_variables
        lower_bound = Inf
        upper_bound = -Inf
        for j in 1:num_boxes
            if vec[j][i].lo < lower_bound
                lower_bound = vec[j][i].lo
            end
            if upper_bound < vec[j][i].hi
                upper_bound = vec[j][i].hi
            end
        end
        append!(x_big, Interval(lower_bound, upper_bound))
    end
    return IntervalBox(x_big...)
end

function check_feasiblity(point::SVector{N,T}, constraints::Vector{Constraint{T}}) where{N, T}
    for constraint in constraints
        if !(invokelatest(constraint.C, point) ∈ constraint.bound)
            return false
        end
    end
    return true
end

function ibc_minimise(f::Function , X::IntervalBox{N,T}, constraints::Vector{Constraint{T}}; ibc_chnl = RemoteChannel(()->Channel{Tuple{IntervalBox{N,T}, T}}(0)),
               diffevol_chnl = RemoteChannel(()->Channel{Tuple{SVector{N, T}, T, Union{Nothing, IntervalBox{N, T}}}}(0)), structure = SortedVector, debug = false, tol=1e-6, ibc_ind = true) where{N, T}

    vars = [Variable(Symbol(:(x[$i])))() for i in 1:N]
    g(x...) = f(x)
    C = BasicContractor(vars, g)

    working = structure([(X, inf(f(X)))], x->x[2]) # list of boxes with corresponding lower bound, arranged according to selected structure :
    minimizers = IntervalBox{N,T}[]
    global_min = ∞  # upper bound

    info = Information(0, 0, 0)
    num_bisections = 0

    while !isempty(working)

        info.iterations= info.iterations + 1

        if !ibc_ind
            if isready(ibc_chnl)
                from_diff = take!(ibc_chnl)     # Receiving best individual from ibc_minimise
                if debug
                    println("Box recevied from DifferentialEvolution: ", from_diff)
                end
                global_min = min(from_diff[2], global_min)
                info.de_to_ibc = info.de_to_ibc + 1
            end
        end

        (X, X_min) = popfirst!(working)

        gra = gradient(f, X.v)
        (lbb, fcb, cb) = bauman_form(X, f, gra)      # ----- second order form

        if global_min > lbb
            X_min = max(X_min, lbb)
            if fcb < global_min && check_feasiblity(cb, constraints)
                global_min = fcb
                x_best = cb
            end
        else
            continue
        end

        if debug
            println("New search-space : ", X)
        end

        X_min, X, constraints = contraction(f, C, global_min, X, constraints)

        if debug
            println("Contracted search_space: ", X)
        end

        if X_min > global_min    # X_min == inf(f(X))
            continue
        end

        status, output = generate_random_feasible_point(X, constraints)   # finding a feasible point in the interval box if there present any

        if status                           # true if we got any feasible point
            feas_point = output
            m = sup(f(Interval.(feas_point)))  # find candidate for upper bound of global minimum by evaluating f on a feasible point
            if m < global_min
                global_min = m
                x_best = SVector(mid(X))
                if !ibc_ind
                    if debug
                        println("Box send to DifferentialEvolution: ", x_best )
                    end
                    if info.iterations % 200 == 0
                       x_big = convex_hull(working.data)
                        put!(diffevol_chnl, (x_best, global_min, x_big))  # sending best individual to diffevol
                    else
                        put!(diffevol_chnl, (x_best, global_min, nothing))
                    end
                    info.ibc_to_de = info.ibc_to_de + 1
                end
            end

            filter_elements!(working , (X, global_min) )   # Remove all boxes whose lower bound is greater than the current one:

            if diam(X) < tol
                push!(minimizers, X)
            else
                X1, X2 = bisect(X)
                push!( working, (X1, inf(f(X1))) )
                push!( working, (X2, inf(f(X2))) )
                num_bisections += 1
            end
        else
            if diam(X) < tol
                push!(minimizers, X)
            else
                X1, X2 = bisect(X)
                push!( working, (X1, inf(f(X1))) )
                push!( working, (X2, inf(f(X2))) )
                num_bisections += 1
            end
        end

    end

    if debug
        println("IBC search is ended")
    end

    if !ibc_ind
        if isready(diffevol_chnl)
            take!(diffevol_chnl)
            put!(diffevol_chnl,(SVector(mid(X)), Inf, nothing) )
        else
            put!(diffevol_chnl,(SVector(mid(X)), Inf, nothing) )
        end
        if debug
            println("DifferentialEvolution is also terminated")
        end
        if isready(ibc_chnl)
            take!(ibc_chnl)
        end
    end

    lower_bound = minimum(inf.(f.(minimizers)))

    return Interval(lower_bound, global_min), minimizers, info

end

function ibc_maximise(f::Function , X::IntervalBox{N,T}, constraints::Vector{Constraint{T}}; ibc_chnl = RemoteChannel(()->Channel{Tuple{IntervalBox{N,T}, T}}(0)),
               diffevol_chnl = RemoteChannel(()->Channel{Tuple{SVector{N, T}, T, Union{Nothing, IntervalBox{N, T}}}}(0)), structure = SortedVector, debug = false, tol=1e-6, ibc_ind = true) where{N, T}
    bound, minimizer, info = ibc_minimise(x -> -f(x), X, constraints, ibc_chnl = ibc_chnl, diffevol_chnl = diffevol_chnl, debug = debug, tol = tol, ibc_ind = ibc_ind)
    return -bound, minimizer, info
end

max_issues_repo_name: yashcodes/Charibde

max_forks_repo_name: yashcodes/Charibde


text:
module DeconvOptim

export gpu_or_cpu

 # to check wether CUDA is enabled
using Requires
# for fast array regularizers
using Tullio
# optional CUDA dependency
include("requires.jl")

 # for optimization
using Optim
 #mean
using Statistics
using StatsBase
using FFTW
FFTW.set_num_threads(12)
using LineSearches

# possible up_sampling
using Interpolations

# for defining custom derivatives
using ChainRulesCore
using LinearAlgebra

using FillArrays

gpu_or_cpu(x) = Array

include("forward_models.jl")
include("lossfunctions.jl")
include("mappings.jl")
# special CUDA regularizers
include("regularizer_cuda.jl")
include("regularizer.jl")
include("utils.jl")
include("conv.jl")
include("generic_invert.jl")
include("lucy_richardson.jl")
include("deconvolution.jl")
include("analysis_tools.jl")

# refresh Zygote to load the custom rrules defined with ChainRulesCore
using Zygote: gradient

# end module
end


text:
module TestVariableMap

using Test

using VecMathOptInterface
const MOI = VecMathOptInterface

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

struct VariableDummyBridge <: MOI.Bridges.Variable.AbstractBridge
    id::Int
end

function MOI.Bridges.Variable.unbridged_map(
    ::VariableDummyBridge,
    ::MOI.VariableIndex,
)
    return nothing
end

function test_empty()
    map = MOI.Bridges.Variable.Map()
    @test isempty(map)
    @test !MOI.Bridges.Variable.has_bridges(map)
    @test isempty(keys(map))
    return
end

"""
    test_map()

TODO(odow): this test is too intertwined. It should be refactored into something
simpler.
"""
function test_map()
    map = MOI.Bridges.Variable.Map()
    b1 = VariableDummyBridge(1)
    set1 = MOI.EqualTo(0.0)
    F1 = MOI.VariableIndex
    S1 = typeof(set1)
    v1, c1 = MOI.Bridges.Variable.add_key_for_bridge(map, () -> b1, set1)
    cannot_unbridge_err = ErrorException(
        "Cannot unbridge function because some variables are bridged by variable" *
        " bridges that do not support reverse mapping, e.g., `ZerosBridge`.",
    )
    @test v1.value == c1.value == -1
    @test MOI.Bridges.Variable.constraint(map, v1) == c1
    @test haskey(map, v1)
    @test map[v1] == b1
    @test MOI.Bridges.Variable.constrained_set(map, v1) == S1
    @test MOI.Bridges.Variable.number_with_set(map, S1) == 1
    @test MOI.Bridges.Variable.constraints_with_set(map, S1) == [c1]
    @test MOI.Bridges.Variable.function_for(map, c1) == v1
    @test_throws cannot_unbridge_err MOI.Bridges.Variable.unbridged_function(
        map,
        MOI.VariableIndex(1),
    )

    @test MOI.Bridges.Variable.number_of_variables(map) == 1
    @test MOI.Bridges.Variable.list_of_constraint_types(map) == Set([(F1, S1)])
    @test length(map) == 1
    @test !isempty(map)
    @test MOI.Bridges.Variable.has_bridges(map)
    @test collect(keys(map)) == [v1]
    @test collect(values(map)) == [b1]

    b2 = VariableDummyBridge(2)
    set2 = MOI.Zeros(4)
    F2 = MOI.VectorOfVariables
    S2 = typeof(set2)
    v2, c2 = MOI.Bridges.Variable.add_keys_for_bridge(map, () -> b2, set2)
    @test v2[1].value == c2.value == -2
    @test MOI.Bridges.Variable.has_keys(map, v2)
    @test !MOI.Bridges.Variable.has_keys(map, v2[4:-1:1])
    for i in 1:4
        @test MOI.Bridges.Variable.constraint(map, v2[i]) == c2
        @test haskey(map, v2[i])
        @test map[v2[i]] == b2
        @test MOI.Bridges.Variable.constrained_set(map, v2[i]) == S2
        @test MOI.Bridges.Variable.length_of_vector_of_variables(map, v2[i]) ==
              4
        @test MOI.Bridges.Variable.index_in_vector_of_variables(map, v2[i]) ==
              MOI.Bridges.IndexInVector(i)
    end
    @test MOI.Bridges.Variable.number_with_set(map, S2) == 1
    @test MOI.Bridges.Variable.constraints_with_set(map, S2) == [c2]
    @test MOI.Bridges.Variable.function_for(map, c2) ==
          MOI.VectorOfVariables(v2)
    @test_throws cannot_unbridge_err MOI.Bridges.Variable.unbridged_function(
        map,
        MOI.VariableIndex(1),
    )

    @test MOI.Bridges.Variable.number_of_variables(map) == 5
    @test MOI.Bridges.Variable.list_of_constraint_types(map) ==
          Set([(F1, S1), (F2, S2)])
    @test length(map) == 2
    @test sprint(MOI.Bridges.print_num_bridges, map) ==
          "\nwith 2 variable bridges"
    @test !isempty(map)
    @test MOI.Bridges.Variable.has_bridges(map)
    @test collect(keys(map)) == [v1; v2]
    @test collect(values(map)) == [b1, b2]

    b3 = VariableDummyBridge(3)
    set3 = MOI.Zeros(0)
    v3, c3 = MOI.Bridges.Variable.add_keys_for_bridge(map, () -> b3, set3)

    @test isempty(v3)
    @test c3.value == 0
    bridges = collect(values(map))
    @test sort([b.id for b in bridges]) == 1:2
    elements = sort(collect(map), by = el -> el.second.id)
    @test elements[1].first == v1
    @test elements[1].second == b1
    @test elements[2].first == v2[1]
    @test elements[2].second == b2

    delete!(map, v1)
    @test MOI.Bridges.Variable.number_of_variables(map) == 4
    @test length(map) == 1
    @test sprint(MOI.Bridges.print_num_bridges, map) ==
          "\nwith 1 variable bridge"
    @test !isempty(map)
    @test MOI.Bridges.Variable.has_bridges(map)
    elements = collect(map)
    @test elements[1].first == v2[1]
    @test elements[1].second == b2
    @test MOI.Bridges.Variable.number_with_set(map, S1) == 0
    @test MOI.Bridges.Variable.number_with_set(map, S2) == 1
    @test isempty(MOI.Bridges.Variable.constraints_with_set(map, S1))
    @test MOI.Bridges.Variable.constraints_with_set(map, S2) == [c2]
    @test MOI.Bridges.Variable.list_of_constraint_types(map) == Set([(F2, S2)])
    @test collect(keys(map)) == v2
    @test collect(values(map)) == [b2]
    @test !haskey(map, v1)
    @test MOI.Bridges.Variable.has_keys(map, v2)

    rev_v2 = reverse(v2)
    err = ArgumentError(
        "`$rev_v2` is not a valid key vector as returned by `add_keys_for_bridge`.",
    )
    @test_throws err delete!(map, rev_v2)

    delete!(map, v2[3])
    left = [1, 2, 4]
    @test MOI.Bridges.Variable.number_of_variables(map) == 3
    @test length(map) == 1
    @test !isempty(map)
    @test MOI.Bridges.Variable.has_bridges(map)
    elements = collect(map)
    @test elements[1].first == v2[1]
    @test elements[1].second == b2
    @test MOI.Bridges.Variable.number_with_set(map, S1) == 0
    @test MOI.Bridges.Variable.number_with_set(map, S2) == 1
    @test isempty(MOI.Bridges.Variable.constraints_with_set(map, S1))
    @test MOI.Bridges.Variable.constraints_with_set(map, S2) == [c2]
    @test MOI.Bridges.Variable.list_of_constraint_types(map) == Set([(F2, S2)])
    @test collect(keys(map)) == v2[left]
    @test collect(values(map)) == [b2]
    @test !haskey(map, v1)
    @test !haskey(map, v2[3])
    @test !MOI.Bridges.Variable.has_keys(map, v2)
    @test MOI.Bridges.Variable.has_keys(map, v2[left])
    for (j, i) in enumerate(left)
        @test MOI.Bridges.Variable.constraint(map, v2[i]) == c2
        @test haskey(map, v2[i])
        @test map[v2[i]] == b2
        @test MOI.Bridges.Variable.constrained_set(map, v2[i]) == S2
        @test MOI.Bridges.Variable.length_of_vector_of_variables(map, v2[i]) ==
              3
        @test MOI.Bridges.Variable.index_in_vector_of_variables(map, v2[i]) ==
              MOI.Bridges.IndexInVector(j)
    end
    @test MOI.Bridges.Variable.function_for(map, c2) ==
          MOI.VectorOfVariables(v2[left])

    delete!(map, v2[1])
    left = [2, 4]
    @test MOI.Bridges.Variable.number_of_variables(map) == 2
    @test length(map) == 1
    @test !isempty(map)
    @test MOI.Bridges.Variable.has_bridges(map)
    elements = collect(map)
    @test elements[1].first == v2[1]
    @test elements[1].second == b2
    @test MOI.Bridges.Variable.number_with_set(map, S1) == 0
    @test MOI.Bridges.Variable.number_with_set(map, S2) == 1
    @test isempty(MOI.Bridges.Variable.constraints_with_set(map, S1))
    @test MOI.Bridges.Variable.constraints_with_set(map, S2) == [c2]
    @test MOI.Bridges.Variable.list_of_constraint_types(map) == Set([(F2, S2)])
    @test collect(keys(map)) == v2[left]
    @test collect(values(map)) == [b2]
    @test !haskey(map, v1)
    @test !haskey(map, v2[1])
    @test !haskey(map, v2[3])
    @test !MOI.Bridges.Variable.has_keys(map, v2)
    @test !MOI.Bridges.Variable.has_keys(map, v2[[2, 3, 4]])
    @test !MOI.Bridges.Variable.has_keys(map, v2[[1, 2, 4]])
    @test MOI.Bridges.Variable.has_keys(map, v2[left])
    for (j, i) in enumerate(left)
        @test MOI.Bridges.Variable.constraint(map, v2[i]) == c2
        @test haskey(map, v2[i])
        @test map[v2[i]] == b2
        @test MOI.Bridges.Variable.constrained_set(map, v2[i]) == S2
        @test MOI.Bridges.Variable.length_of_vector_of_variables(map, v2[i]) ==
              2
        @test MOI.Bridges.Variable.index_in_vector_of_variables(map, v2[i]) ==
              MOI.Bridges.IndexInVector(j)
    end
    @test MOI.Bridges.Variable.function_for(map, c2) ==
          MOI.VectorOfVariables(v2[left])

    delete!(map, v2[left])
    @test MOI.Bridges.Variable.number_of_variables(map) == 0
    @test length(map) == 0
    @test isempty(map)
    @test MOI.Bridges.Variable.has_bridges(map)
    @test isempty(collect(map))
    @test MOI.Bridges.Variable.number_with_set(map, S1) == 0
    @test MOI.Bridges.Variable.number_with_set(map, S2) == 0
    @test isempty(MOI.Bridges.Variable.constraints_with_set(map, S1))
    @test isempty(MOI.Bridges.Variable.constraints_with_set(map, S2))
    @test isempty(MOI.Bridges.Variable.list_of_constraint_types(map))
    @test isempty(keys(map))
    @test isempty(values(map))
    @test !haskey(map, v1)
    @test !MOI.Bridges.Variable.has_keys(map, v2)
    @test !MOI.Bridges.Variable.has_keys(map, v2[left])
    return
end

function test_EmptyMap()
    S1 = MOI.EqualTo{Float64}
    S2 = MOI.Zeros
    map = MOI.Bridges.Variable.EmptyMap()
    empty!(map)
    @test isempty(map)
    @test length(map) == 0
    @test isempty(keys(map))
    @test isempty(values(map))
    @test !MOI.Bridges.Variable.has_bridges(map)
    @test iszero(MOI.Bridges.Variable.number_of_variables(map))
    @test iszero(MOI.Bridges.Variable.number_with_set(map, S1))
    @test iszero(MOI.Bridges.Variable.number_with_set(map, S2))
    @test isempty(MOI.Bridges.Variable.constraints_with_set(map, S1))
    @test isempty(MOI.Bridges.Variable.constraints_with_set(map, S2))
    @test sprint(MOI.Bridges.print_num_bridges, map) == ""
    return
end

end  # module

TestVariableMap.runtests()


text:
# TODO: create test functions for the methods

# """
#     Left hand side of eqn(6) in https://arxiv.org/pdf/1703.00443.pdf
# """
# function create_LHS_matrix(z, λ, Q, G, h, A=nothing)
#     if A == nothing || size(A)[1] == 0
#         return [Q                G';
#                 Diagonal(λ) * G    Diagonal(G * z - h)]
#     else
#         @assert size(A)[2] == size(G)[2]
#         p, n = size(A)
#         m    = size(G)[1]
#         return [Q                  G'                    A';
#                 Diagonal(λ) * G    Diagonal(G * z - h)   zeros(m, p);
#                 A                  zeros(p, m)           zeros(p, p)]
#     end
# end

"""
    create_LHS_matrix(z, λ, Q, G, h, A=nothing)

Inverse matrix specified on RHS of eqn(7) in https://arxiv.org/pdf/1703.00443.pdf

Helper method while calling [`backward!`](@ref)
"""
function create_LHS_matrix(z, λ, Q, G, h, A=nothing)
    if A === nothing || size(A)[1] == 0
        return [Q         G' * Diagonal(λ);
                G         Diagonal(G * z - h)]
    else
        p, n = size(A)
        m    = size(G, 1)
        if n != size(G, 2)
            throw(DimensionError("Sizes of $A and $G do not match"))
        end
        return [Q         G' * Diagonal(λ)       A';
                G         Diagonal(G * z - h)    spzeros(m, p);
                A         spzeros(p, m)          spzeros(p, p)]
    end
end

# """
#     Right hand side of eqn(6) in https://arxiv.org/pdf/1703.00443.pdf
# """
# function create_RHS_matrix(z, dQ, dq, λ, dG, dh, ν=nothing, dA=nothing, db=nothing)
#     if dA == nothing || size(dA)[1] == 0
#         return -[dQ * z + dq + dG' * λ      ;
#                  Diagonal(λ) * (dG * z - dh)]
#     else
#         return -[dQ * z + dq + dG' * λ + dA' * ν;
#                  Diagonal(λ) * (dG * z - dh)    ;
#                  dA * z - db                    ]
#     end
# end

is_equality(set::MOI.AbstractSet) = false
is_equality(set::MOI.EqualTo) = true

"""
    get_problem_data(model::MOI.AbstractOptimizer)

Return problem parameters as matrices along with other program info such as number of constraints, variables, etc
"""
function get_problem_data(model::MOI.AbstractOptimizer)
    var_list = MOI.get(model, MOI.ListOfVariableIndices())
    nz = size(var_list, 1)

    # handle inequality constraints
    ineq_con_idx = MOI.get(
                        model,
                        MOI.ListOfConstraintIndices{
                            MOI.ScalarAffineFunction{Float64},
                            MOI.LessThan{Float64}
                        }())
    nineq = length(ineq_con_idx)

    G = spzeros(nineq, nz)
    h = spzeros(nineq)

    for i in 1:nineq
        con = ineq_con_idx[i]

        func = MOI.get(model, MOI.ConstraintFunction(), con)
        set = MOI.get(model, MOI.ConstraintSet(), con)

        for (j, var_idx) in enumerate(var_list)
            for term in func.terms
                if term.variable_index == var_idx
                    G[i,j] = MOI.coefficient(term)
                end
            end
        end
        h[i] = set.upper - func.constant
    end

    # handle equality constraints
    eq_con_idx = MOI.get(
                        model,
                        MOI.ListOfConstraintIndices{
                            MOI.ScalarAffineFunction{Float64},
                            MOI.EqualTo{Float64}
                        }())
    neq = length(eq_con_idx)

    A = spzeros(neq, nz)
    b = spzeros(neq)

    for i in 1:neq
        con = eq_con_idx[i]

        func = MOI.get(model, MOI.ConstraintFunction(), con)
        set = MOI.get(model, MOI.ConstraintSet(), con)

        for x in func.terms
            A[i, x.variable_index.value] = x.coefficient
        end
        b[i] = set.value - func.constant
    end

    # handle objective
    # works both for affine and quadratic objective functions
    objective_function = MOI.get(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}())
    Q = spzeros(nz, nz)
    q = spzeros(nz)

    if objective_function isa MathOptInterface.ScalarAffineFunction{Float64}
        for x in objective_function.terms
            q[x.variable_index.value] = x.coefficient
        end
    elseif objective_function isa MathOptInterface.ScalarQuadraticFunction{Float64}

        var_to_id = Dict(var_list .=> 1:nz)

        for quad in objective_function.quadratic_terms
            i = var_to_id[quad.variable_index_1]
            j = var_to_id[quad.variable_index_2]
            Q[i,j] = quad.coefficient
            Q[j,i] = quad.coefficient
        end

        q = MOI.coefficient.(objective_function.affine_terms)
    end

    return Q, q, G, h, A, b, nz, var_list, nineq, ineq_con_idx, neq, eq_con_idx
end

# might slow down computation
# need to find a faster way
function CSRToCSC(B::MatOI.SparseMatrixCSRtoCSC{Int64})
    A = spzeros(B.m, B.n)
    last = 0
    for i in 1:B.n
        rnge = (last+1):B.colptr[i]
        A[(1 .+ B.rowval[rnge]), i] = B.nzval[rnge]
        last = B.colptr[i]
    end
    return A
end