text:
using Distributed
@everywhere using JuMP, Gurobi
@everywhere function kmer_coverage_ilp(locus, kmer_class, allele_ids, coverages)
  # m = Model(solver=GurobiSolver(Presolve=0))
  # one minute at most per locus, gap=0.1, no need to have optimal per se;
  # One thread per solve, since every ILP is already being run in parallel.
  m = Model(solver=GurobiSolver(OutputFlag=0, MIPGap=0.1, TimeLimit=1200, Threads=1))
  n_kmer = length(kmer_class)
  n_alleles = length(allele_ids)
  # kmer decision variable:
  @variable(m, x[1:n_kmer], Bin)
  @variable(m, c[1:n_alleles], Int)
  # coverage constraints:
  kmer_list = collect(keys(kmer_class))
  kmer_cover = DefaultDict{Int,Vector{Int}}(() -> Int[]) # allele -> list of kmers that cover it
  for (i, kmer) in enumerate(kmer_list)
    for al in kmer_class[kmer]
      push!(kmer_cover[al], i)
    end
  end
  # define coverage:
  for j in 1:n_alleles
    @constraint(m, c[j] == sum(x[i] for i in kmer_cover[j]))
  end
  # min coverage:
  for j in 1:n_alleles
    @constraint(m, c[j] >= coverages[j])
    # @constraint(m, c[j] == coverages[j])
  end
  # Forcing same coverage: (takes much longer, get many more kmers; better to store coverage)
  # Update: consider the cardinality of the kmers, to minimize it also:
  # @variable(m, max_card, Int) # max cardinality of the selected kmers
  # for (i, kmer) in enumerate(kmer_list)
  #   @constraint(m, max_card >= x[i] * length(kmer_class[kmer]))
  # end
  # minimize # of selected kmers:
  @objective(m, Min, sum(x[i] for i in 1:n_kmer))
  # @objective(m, Min, sum(x[i] for i in 1:n_kmer) + 100000*max_card)
  # TODO: e parameters works well, reduces variance but increases kmer number as expected.
  status = solve(m)
  if status == :Infeasible
    return nothing
  end
  # println("Objective value: ", getobjectivevalue(m))
  # for (i,kmer) in enumerate(kmer_list)
  #   if getvalue(x[i]) == 1
  #     println("$kmer : $(kmer_class[kmer])")
  #   end
  # end
  # println(join([Int(round(getvalue(c[j]),0)) for j in 1:n_alleles],", "))
  try
    kmers = [kmer_list[i] for i in 1:n_kmer if getvalue(x[i]) == 1]
    coverages = [Int16(round(getvalue(c[j]))) for j in 1:n_alleles]
    return kmers, coverages
  catch
    println("What happened??? -> LP status: $status LOCUS:$locus.")
    exit()
  end
  # return kmers, coverages, Int16(round(getvalue(max_card)))
end

max_stars_repo_name: liberuser/MentaLiST

max_issues_repo_name: liberuser/MentaLiST

max_forks_repo_name: liberuser/MentaLiST


text:
# Poly type manipulations

module Polynomials
#todo: sparse polynomials?

using Compat

export Poly, poly
export degree, coeffs, variable
export polyval, polyint, polyder, roots, polyfit
export Pade, padeval
export truncate!

import Base: length, endof, getindex, setindex!, copy, zero, one, convert, norm, gcd
import Base: show, print, *, /, //, -, +, ==, divrem, div, rem, eltype, .*, .-, .+
import Base: promote_rule, truncate, chop
if VERSION >= v"0.4"
    import Base.call
end

eps{T}(::Type{T}) = zero(T)
eps{F<:AbstractFloat}(x::Type{F}) = Base.eps(F)
eps{T}(x::Type{Complex{T}}) = eps(T)

"""

* `Poly{T<:Number}(a::Vector)`: Construct a polynomial from its coefficients, lowest order first. That is if `p=a_n x^n + ... + a_2 x^2 + a_1 x^1 + a_0`, we construct this through `Poly([a_0, a_1, ..., a_n])`.

Example:
```
Poly([1,0,3,4])    # Poly(1 + 3x^2 + 4x^3)
```

An optional variable parameter can be added:

```
Poly([1,2,3], :s)       # Poly(1 + 2s + 3s^2)
```

The usual arithmetic operators are overloaded to work on polynomials, and combinations of polynomials and scalars.

```
p = Poly([1,2])        # Poly(1 + 2x)
q = Poly([1, 0, -1])   # Poly(1 - x^2)
2p                     # Poly(2 + 4x)
2+p                    # Poly(3 + 2x)
p - q                  # Poly(2x + x^2)
p*q                    # Poly(1 + 2x - x^2 - 2x^3)
q/2                    # Poly(0.5 - 0.5x^2)
```

Note that operations involving polynomials with different variables will error:

```j
p = Poly([1, 2, 3], :x)
q = Poly([1, 2, 3], :s)
p + q                  # ERROR: Polynomials must have same variable.
```

"""
immutable Poly{T<:Number}
    a::Vector{T}
    var::Symbol
    function Poly(a::Vector{T}, var)
        # if a == [] we replace it with a = [0]
        if length(a) == 0
            return new(zeros(T,1), @compat Symbol(var))
        else
            # determine the last nonzero element and truncate a accordingly
            a_last = max(1,findlast(a))
            new(a[1:a_last], @compat Symbol(var))
        end
    end
end

@compat Poly{T<:Number}(a::Vector{T}, var::Union{AbstractString,Symbol,Char}=:x) = Poly{T}(a, var)

# create a Poly object from its roots
"""

* `poly(r::AbstractVector)`: Construct a polynomial from its
  roots. This is in contrast to the `Poly` constructor, which
  constructs a polynomial from its coefficients.

Example:
```
## Represents (x-1)*(x-2)*(x-3)
poly([1,2,3])     # Poly(-6 + 11x - 6x^2 + x^3)
```
"""
function poly{T}(r::AbstractVector{T}, var=:x)
    n = length(r)
    c = zeros(T, n+1)
    c[1] = 1
    for j = 1:n
        for i = j:-1:1
            c[i+1] = c[i+1]-r[j]*c[i]
        end
    end
    return Poly(reverse(c), var)
end
poly(A::Matrix, var=:x) = poly(eigvals(A), var)
poly(A::Matrix, var::AbstractString) = poly(eigvals(A), @compat Symbol(var))
poly(A::Matrix, var::Char) = poly(eig(A)[1], @compat Symbol(var))

include("show.jl") # display polynomials.

convert{T}(::Type{Poly{T}}, p::Poly) = Poly(convert(Vector{T}, p.a), p.var)
convert{T, S<:Number}(::Type{Poly{T}}, x::S) = Poly(promote_type(T, S)[x])
convert{T, S<:Number,n}(::Type{Poly{T}}, x::Array{S,n}) = map(el->convert(Poly{promote_type(T,S)},el),x)
promote_rule{T, S}(::Type{Poly{T}}, ::Type{Poly{S}}) = Poly{promote_type(T, S)}
promote_rule{T, S<:Number}(::Type{Poly{T}}, ::Type{S}) = Poly{promote_type(T, S)}
eltype{T}(::Poly{T}) = T

"""

`legnth(p::Poly)`: return length of coefficient vector

"""
length(p::Poly) = length(p.a)
endof(p::Poly) = length(p) - 1

"""

`degree(p::Poly)`: return degree of polynomial `p`

"""
degree(p::Poly) = length(p) - 1

"""

`coeffs(p::Poly)`: return coefficient vector [a_0, a_1, ..., a_n]

"""
coeffs(p::Poly) = p.a

"""

Return the indeterminate of a polynomial, `x`.

* `variable(p::Poly)`: return variable of `p` as a `Poly` object.
* `variable(T<:Number, [:x])`: return poly one(T)*x
* `variable([var::Symbol])`: return polynomial 1x over `Float64`.

"""
variable{T}(p::Poly{T}) = poly(zeros(T,1), p.var)
variable{T<:Number}(::Type{T}, var=:x) = poly(zeros(T,1), var)
variable(var::Symbol=:x) = poly([0.0], var)

"""

`truncate{T}(p::Poly{T}; reltol = eps(T), abstol = eps(T))`: returns a polynomial with coefficients a_i truncated to zero if |a_i| <= reltol*maxabs(a)+abstol

"""
function truncate{T}(p::Poly{T}; reltol = eps(T), abstol = eps(T))
    a = coeffs(p)
    amax = maxabs(a)
    anew = map(ai -> abs(ai) <= amax*reltol+abstol ? zero(T) : ai, a)
    return Poly(anew, p.var)
end

"""

`chop(p::Poly{T}; kwargs...)` chop off leading values which are
approximately zero. The tolerances are passed to `isapprox`.

"""
function chop{T}(p::Poly{T}; reltol=zero(T), abstol=2 * eps(T))
    c = copy(p.a)
    for k=length(c):-1:1
        if !isapprox(c[k], zero(T); rtol=reltol, atol=abstol)
            resize!(c, k)
            return Poly(c, p.var)
        end
    end

    resize!(c,0)
    Poly(c, p.var)
end

"""

* `norm(q::Poly, [p])`: return `p` norm of polynomial `q`

"""
norm(q::Poly, args...) = norm(coeffs(q), args...)

"""

* `getindex(p::Poly, i)`: If `p=a_n x^n + a_{n-1}x^{n-1} + ... + a_1 x^1 + a_0`, then `p[i]` returns `a_i`.

"""
getindex{T}(p::Poly{T}, i) = (i+1 > length(p.a) ? zero(T) : p.a[i+1])
getindex{T}(p::Poly{T}, idx::AbstractArray) = map(i->p[i], idx)
function setindex!(p::Poly, v, i)
    n = length(p.a)
    if n < i+1
        resize!(p.a,i+1)
        p.a[n+1:i] = 0
    end
    p.a[i+1] = v
    v
end
function setindex!(p::Poly, vs, idx::AbstractArray)
    [setindex!(p, v, i) for (i,v) in zip(idx, vs)]
    p
end

copy(p::Poly) = Poly(copy(p.a), p.var)

zero{T}(p::Poly{T}) = Poly([zero(T)], p.var)
zero{T}(::Type{Poly{T}}) = Poly(T[])
one{T}(p::Poly{T}) = Poly([one(T)], p.var)
one{T}(::Type{Poly{T}}) = Poly([one(T)])

## Overload arithmetic operators for polynomial operations between polynomials and scalars
*{T<:Number,S}(c::T, p::Poly{S}) = Poly(c * p.a, p.var)
*{T<:Number,S}(p::Poly{S}, c::T) = Poly(p.a * c, p.var)
.*{T<:Number,S}(c::T, p::Poly{S}) = Poly(c * p.a, p.var)
.*{T<:Number,S}(p::Poly{S}, c::T) = Poly(p.a * c, p.var)
/(p::Poly, c::Number) = Poly(p.a / c, p.var)
-(p::Poly) = Poly(-p.a, p.var)
-{T<:Number}(p::Poly, c::T) = +(p, -c)
.-{T<:Number}(p::Poly, c::T) = +(p, -c)
.-{T<:Number}(c::T, p::Poly) = +(p, -c)
+{T<:Number}(c::T, p::Poly) = +(p, c)
.+{T<:Number}(c::T, p::Poly) = +(p, c)
.+{T<:Number}(p::Poly, c::T) = +(p, c)
function +{T<:Number}(p::Poly, c::T)
    if length(p) < 1
        return Poly([c,], p.var)
    else
        p2 = copy(p)
        p2.a[1] += c;
        return p2;
    end
end
function -{T<:Number}(c::T, p::Poly)
    if length(p) < 1
        return Poly(T[c,], p.var)
    else
        p2 = -p;
        p2.a[1] += c;
        return p2;
    end
end

function +{T,S}(p1::Poly{T}, p2::Poly{S})
    if p1.var != p2.var
        error("Polynomials must have same variable")
    end
    Poly([p1[i] + p2[i] for i = 0:max(length(p1),length(p2))], p1.var)
end
function -{T,S}(p1::Poly{T}, p2::Poly{S})
    if p1.var != p2.var
        error("Polynomials must have same variable")
    end
    Poly([p1[i] - p2[i] for i = 0:max(length(p1),length(p2))], p1.var)
end

function *{T,S}(p1::Poly{T}, p2::Poly{S})
    if p1.var != p2.var
        error("Polynomials must have same variable")
    end
    R = promote_type(T,S)
    n = length(p1)-1
    m = length(p2)-1
    a = zeros(R,m+n+1)

    for i = 0:n
        for j = 0:m
            a[i+j+1] += p1[i] * p2[j]
        end
    end
    Poly(a,p1.var)
end

function divrem{T, S}(num::Poly{T}, den::Poly{S})
    if num.var != den.var
        error("Polynomials must have same variable")
    end
    m = length(den)-1
    if m == 0 && den[0] == 0
        throw(DivideError())
    end
    R = typeof(one(T)/one(S))
    n = length(num)-1
    deg = n-m+1
    if deg <= 0
        return convert(Poly{R}, zero(num)), convert(Poly{R}, num)
    end

    aQ = zeros(R, deg)
    # aR = deepcopy(num.a)
    # @show num.a
    aR = R[ num.a[i] for i = 1:n+1 ]
    for i = n:-1:m
        quot = aR[i+1] / den[m]
        aQ[i-m+1] = quot
        for j = 0:m
            elem = den[j]*quot
            aR[i-(m-j)+1] -= elem
        end
    end
    pQ = Poly(aQ, num.var)
    pR = Poly(aR, num.var)

    return pQ, pR
end

@deprecate /(num::Poly, den::Poly)                div(num::Poly, den::Poly)
div(num::Poly, den::Poly) = divrem(num, den)[1]
rem(num::Poly, den::Poly) = divrem(num, den)[2]

function ==(p1::Poly, p2::Poly)
    if p1.var != p2.var
        return false
    else
        return p1.a == p2.a
    end
end

"""
* `polyval(p::Poly, x::Number)`: Evaluate the polynomial `p` at `x` using Horner's method.

Example:
```
polyval(Poly([1, 0, -1]), 0.1)  # 0.99
```

For `julia` version `0.4` or greater, the `call` method can be used:

```
p = Poly([1,2,3])
p(4)   # 57 = 1 + 2*4 + 3*4^2
```

"""
function polyval{T,S}(p::Poly{T}, x::S)
    R = promote_type(T,S)
    lenp = length(p)
    if lenp == 0
        return zero(R) * x
    else
        y = convert(R, p[end]) + 0*x
        for i = (endof(p)-1):-1:0
            y = p[i] + x*y
        end
        return y
    end
end

polyval(p::Poly, v::AbstractVector) = map(x->polyval(p, x), v)

if VERSION >= v"0.4"
    @compat (p::Poly)(x) = polyval(p, x)
end

"""

* `polyint(p::Poly, k::Number=0)`: Integrate the polynomial `p` term
  by term, optionally adding constant term `k`. The order of the
  resulting polynomial is one higher than the order of `p`.

Examples:
```
polyint(Poly([1, 0, -1]))     # Poly(x - 0.3333333333333333x^3)
polyint(Poly([1, 0, -1]), 2)  # Poly(2.0 + x - 0.3333333333333333x^3)
```

"""
function polyint{T}(p::Poly{T}, k::Number=0)
    n = length(p)
    R = typeof(one(T)/1)
    a2 = Array(R, n+1)
    a2[1] = k
    for i = 1:n
        a2[i+1] = p[i-1] / i
    end
    return Poly(a2, p.var)
end

"""

* `polyder(p::Poly)`: Differentiate the polynomial `p` term by
  term. The order of the resulting polynomial is one lower than the
  order of `p`.

Example:
```
polyder(Poly([1, 3, -1]))   # Poly(3 - 2x)
```
"""
function polyder{T}(p::Poly{T}, order::Int=1)
    n = length(p)
    if order < 0
        error("Order of derivative must be non-negative")
    elseif n <= order
        return Poly(zeros(T,0),p.var)
    elseif order == 0
        return p
    else
        a2 = Array(T, n-order)
        for i = order:n-1
            a2[i-order+1] = p[i] * prod((i-order+1):i)
        end
        return Poly(a2, p.var)
    end
end
Base.ctranspose{T}(p::Poly{T}) = polyder(p)

polyint{T}(a::Array{Poly{T},1}, k::Number  = 0) = [ polyint(p,k) for p in a ]
polyder{T}(a::Array{Poly{T},1}, order::Int = 1) = [ polyder(p,order) for p in a ]
polyint{n,T}(a::Array{Poly{T},n}, k::Number  = 0) = map(p->polyint(p,k),a)
polyder{n,T}(a::Array{Poly{T},n}, order::Int = 1) = map(p->polyder(p,order),a)

##################################################
##
## Some functions on polynomials...

# compute the roots of a polynomial
"""

* `roots(p::Poly)`: Return the roots (zeros) of `p`, with
  multiplicity. The number of roots returned is equal to the order of
  `p`. The returned roots may be real or complex.

Examples:
```
roots(Poly([1, 0, -1]))    # [-1.0, 1.0]
roots(Poly([1, 0, 1]))     # [0.0+1.0im, 0.0-1.0im]
roots(Poly([0, 0, 1]))     # [0.0, 0.0]
roots(poly([1,2,3,4]))     # [1.0,2.0,3.0,4.0]
```
"""
function roots{T}(p::Poly{T})
    R = promote_type(T, Float64)
    length(p) == 0 && return zeros(R, 0)

    num_leading_zeros = 0
    while abs(p[num_leading_zeros]) <= 2*eps(T)
        if num_leading_zeros == length(p)-1
            return zeros(R, 0)
        end
        num_leading_zeros += 1
    end

    num_trailing_zeros = 0
    while abs(p[end - num_trailing_zeros]) <= 2*eps(T)
        num_trailing_zeros += 1
    end

    n = endof(p)-(num_leading_zeros + num_trailing_zeros)

    n < 1 && return zeros(R, length(p) - num_trailing_zeros - 1)

    companion = zeros(R, n, n)
    an = p[end-num_trailing_zeros]
    for i = 1:n-1
        companion[i,n] = -p[num_leading_zeros + i - 1] / an
        companion[i+1,i] = 1;
    end
    companion[end,end] = -p[end-num_trailing_zeros-1] / an
    D = eigvals(companion)
    r = zeros(eltype(D),length(p)-num_trailing_zeros-1)
    r[1:n] = D
    return r
end
roots{T}(p::Poly{Rational{T}}) = roots(convert(Poly{promote_type(T, Float64)}, p))

## compute gcd of two polynomials
"""

* `gcd(a::Poly, b::Poly)`: Finds the Greatest Common Denominator of
    two polynomials recursively using [Euclid's
    algorithm](http://en.wikipedia.org/wiki/Polynomial_greatest_common_divisor#Euclid.27s_algorithm).

Example:
```
gcd(poly([1,1,2]), poly([1,2,3])) # returns (x-1)*(x-2)
```
"""
function gcd{T, S}(a::Poly{T}, b::Poly{S})
    if all(abs(b.a).<=2*eps(S))
        return a
    else
        s, r = divrem(a, b)
        return gcd(b, r)
    end
end

## Fit degree n polynomial to points
"""

`polyfit(x, y, n=length(x)-1, sym=:x )`: Fit a polynomial of degree
`n` through the points specified by `x` and `y` where `n <= length(x)
- 1` using least squares fit. When `n=length(x)-1` (the default), the
interpolating polynomial is returned. The optional fourth argument can
be used to pass the symbol for the returned polynomial.

Example:

```
xs = linspace(0, pi, 5)
ys = map(sin, xs)
polyfit(xs, ys, 2)
```

Original by [ggggggggg](https://github.com/Keno/Polynomials.jl/issues/19)
"""
function polyfit(x, y, n::Int=length(x)-1, sym::Symbol=:x)
    length(x) == length(y) || throw(DomainError)
    1 <= n <= length(x) - 1 || throw(DomainError)

    A = [ float(x[i])^p for i = 1:length(x), p = 0:n ]
    Poly(A \ y, sym)
end
polyfit(x,y,sym::Symbol) = polyfit(x,y,length(x)-1, sym)

### Pull in others
include("pade.jl")

end # module Poly

ext: jl


text:
## Time stepping methods ##

@inline euler_forward!(eq::PrognosticEquation, s, dt) =
    euler_forward!(s.inc[eq.next], s.state[eq.var], s.inc[eq.inc], dt)

@inline function euler_forward!(next::T, old::T, inc::T, dt) where {T}
    n, m = size(old)
    @inbounds for j in 1:m
        for i in 1:n
            next[i, j] = dt * inc[i, j] + old[i, j]
        end
    end
end

@inline integrate!(eq::PrognosticEquation, state, dt) = eq.integrator(eq, state, dt)
@inline function integrate!(eq, state, dt) end
@inline integrate!(model, dt) = map(eq -> integrate!(eq, model.s, dt), model.equations)


text:
module CUDA

    import Base.length, Base.size
    import Base.copy!

    export

    # errors
    CuDriverError, description,

    # base
    @cucall,

    # devices
    CuDevice, CuCapability, dev_count, name, totalmem, attribute, capability,
    list_devices,

    # context
    CuContext, create_context, destroy, push, pop,

    # module
    CuModule, CuFunction, unload,

    # stream
    CuStream, synchronize,

    # execution
    launch,

    # arrays
    CuPtr, CuArray, free, to_host

    const CUDA_LIB = @windows? "nvcuda.dll" : "libcuda"
    Libdl.dlopen(CUDA_LIB)    # loads library, throws an error if not found

    include("errors.jl")
    include("funmap.jl")

    include("base.jl")
    include("devices.jl")
    include("context.jl")
    include("module.jl")
    include("stream.jl")
    include("execution.jl")

    include("arrays.jl")
end


text:
using NNlib
import NNlib: conv!, ∇conv_filter!, ∇conv_data!,
  maxpool!, meanpool!, ∇maxpool!, ∇meanpool!,
  softmax, softmax!, ∇softmax!, logsoftmax, logsoftmax!, ∇logsoftmax
using CUDAnative

# Activation functions
@cufunc σ(x) = ifelse(x < -80, zero(x), one(x) / (one(x) + exp(-x)))

@cufunc function logσ(x)
  max_v = max(zero(x), -x)
  z = exp(-max_v) + exp(-x-max_v)
  -(max_v + log(z))
end

@cufunc elu(x, α = one(x)) =
  ifelse(x ≥ 0, x/1, α * (exp(x) - one(x)))

# TODO: make @cufunc recognise its own definitions
cufunc(::typeof(swish)) = x -> x * cufunc(σ)(x)

@cufunc function selu(x)
  λ = oftype(x/1, 1.0507009873554804934193349852946)
  α = oftype(x/1, 1.6732632423543772848170429916717)
  λ * ifelse(x > 0, x/1, α * (exp(x) - 1))
end

@cufunc softplus(x) = log1p(exp(x))

if !@isdefined CUDNN
  function conv!(y::CuArray, x::CuArray, w::CuArray; kw...)
    error("CUDNN is not installed.")
  end
  function softmax!(out::CuVecOrMat, xs::CuVecOrMat)
    error("CUDNN is not installed.")
  end
  function logsoftmax!(out::CuVecOrMat, xs::CuVecOrMat)
    error("CUDNN is not installed.")
  end
end


text:
# Functions implementing the time evolution
# We only implement (which seems to work) the velocity verlet algorithm.

#TODO: implement the Extended Lagrangian based methods and integrators

function  velocity_verlet(a::Function,   # function providing the acceleration
                                            x::Array{T,1}, # current position
                                            v::Array{T,1}, # current speed
                                            vdot::Array{T,1}, # current acceleration
                                            h::T) where T <: Number           # delta time

# modification to have only one evaluation of a per iteration
    v_half = v      + h/2 * vdot
    x_next = x      + h   * v_half
    vdot_next = a(x_next)
    v_next = v_half + h/2 * vdot_next
    return (x_next, v_next, vdot_next)
end

function time_evolution(Integrator::Function, # integrator
                        a::Function,          # acceleration
                        h::T,                 # time step
                        Nsteps::Int64,        # number of stepd for evolution
                        x0::Array{T,1},       # initial guess
                        x1::Array{T,1} ) where T <: Number      # initial+1 guess

  xArray = Array{T,2}(undef, Nsteps+1, length(x0))
  vArray = Array{T,2}(undef, Nsteps+1, length(x0))
  aArray = Array{T,2}(undef, Nsteps+1, length(x0))

  println("initial conditions")
  xArray[1,:] =  x0
  xArray[2,:] =  x1

  vArray[1,:] = (x1 - x0)/(h) # is this enough?
  aArray[1,:] = a(x0)

  println("running the loop")

  for ii = 1:Nsteps
      (x, v, vdot) = Integrator(a,xArray[ii,:], vArray[ii,:], aArray[ii,:], h)
      xArray[ii+1,:] = x
      vArray[ii+1,:] = v
      aArray[ii+1,:] = vdot
  end
  println("loop finished")

  return (xArray, vArray, aArray)
end

max_stars_repo_name: Forgotten/modelscf

max_issues_repo_name: Forgotten/modelscf

max_forks_repo_name: Forgotten/modelscf


text:
abstract type AbstractDiscretizer end

struct QuantileSlicingDiscretizer <: AbstractDiscretizer
    nbins::Int64
end

function learn_discretizations(
    disc::QuantileSlicingDiscretizer,
    Ps_train,
    Xs_train::AbstractVector{<:Real},
    Xs_test::AbstractVector{<:Real},
    α,
    multiple_testing_method,
)

    @unpack nbins = disc
    qs = [-Inf; quantile([Xs_train; Xs_test], (1:(nbins-1)) ./ nbins); +Inf]
    Xs_train_d = cut(Xs_train, qs, allowempty = true)
    Xs_test_d = cut(Xs_test, qs, allowempty = true)
    (Xs_train_d, Xs_test_d, qs)
end

struct DiscretizationWeightLearner{CW<:CategoricalWeightLearner,D<:AbstractDiscretizer} <:
       WeightLearner
    wlearner::CW
    discretizer::D
end

function learn_weights(
    dw::DiscretizationWeightLearner,
    Ps_train,
    Xs_train,
    Xs_test,
    α,
    multiple_testing_method,
)
    # TODO: check levels are the same

    Xs_train_d, Xs_test_d, discr_scheme = learn_discretizations(
        dw.discretizer,
        Ps_train,
        Xs_train,
        Xs_test,
        α,
        multiple_testing_method,
    )

    ws, ws_fit = learn_weights(
        dw.wlearner,
        Ps_train,
        Xs_train_d,
        Xs_test_d,
        α,
        multiple_testing_method,
    )
    discr_ws_fit = (discr_fit = discr_scheme, ws_fit = ws_fit)
    (ws, discr_ws_fit)
end


text:
"""
#### function ```get_line_params(xi::Float64, df::DataFrame)```

    Description of ```get_line_params```
    -----------------------------------------
Returns the parameters a,b of linear approximation y = ax +b. For the segment where the random number ξ falls into.
``
"""
function get_line_params(xi::Float64, df::DataFrame)
    row = df[(df.minE .< xi) .& (df.maxE .> xi),:]
    return (row.a[1], row.b[1])
end

"""
#### function ```get_line_params(x1::Float64, y1::Float64, x2::Float64, y2::Float64)```

    Description of ```get_line_params```
    -----------------------------------------
Returns the parameters a,b of linear approximation y = ax +b.
``
"""
function get_line_params(x1::Float64, y1::Float64, x2::Float64, y2::Float64)
    a = (y1 - y2) / (x1 - x2)
    b = y1 - a*x1
    return a,b #returns the parameters of line given by y = b + ax
end

get_line_point(x::Float64,a::Float64,b::Float64) = b + a*x

"""
#### function ```get_integral_linear(minE::Real, maxE::Real, a::Real, b::Real)```

    Description of ```get_integral_linear```
    -----------------------------------------
Return the definitive integral of ``
    \\int_{minE}^{maxE} (a_iE + b_i) dE= \\frac{a_i}{2}(maxE^2 - minE^2) + b(maxE - minE)
``
"""
function get_integral_linear(minE::Real, maxE::Real, a::Real, b::Real)
    return 0.5*a * (maxE^2 - minE^2) + b * ( maxE - minE )
end

function solvequadratic(a, b, c)
    d  = sqrt(b^2 - 4*a*c)
    return (-b - d) / (2*a), (-b + d) / (2*a)
end

function plot_lines(step::Real, df::DataFrame, c::ColorPalette)
    p::Plots.Plot = plot()
    cp::Float64 = 1.0
    for e1 in 1:step:length(unique(df.E1))
        for row in eachrow(df[df.E1 .== unique(df.E1)[e1],:])
            xs = row.minE : 1e-3 : row.maxE
            line2(x) = get_line_point(x, row.a, row.b)
            plot!(xs, line2.(xs), lw = 2, alpha = 0.3, c = c[ceil(Int,cp)], legend = :false,
                    xlabel = "E2 [MeV]", ylabel ="dGdE", ylims = (minimum(df.minG), 1.1*maximum(df.maxG)),
                    title = "projection of the linear approximations, \n every $step")
        end
        if cp < 250.0
            cp += step/10.0
        end

    end
    return p
end


text:
@testset "Arithmetic" begin
    Decimals.setprecision!(precision=4)
    @test d"1" / 3 == 0.3333

    Decimals.setprecision!(precision=5)
    @test d"1" / 3 == 0.33333 == Decimals.One / Decimal(3)

    Decimals.setprecision!(precision=3, rounding=RoundUp)
    @test one(Decimal) * d"4.2" == 4.2 == d"42E-1"
    @test 4.2 * d"4.2" == 4.2 * 4.2 == d"42E-1" * d"4.2"
    @test d"1" / 3 == 0.334 == Decimals.One / Decimal(3)
    @test d"-1" / 3 == -0.334 == -Decimals.One / Decimal(3)
    @test d"1" + d"2" == d"3" == d"40E-1" - d"100E-2"
    @test Decimals.setexponent(d"3.3333", -3, RoundDown) ==  d"3.333"
    @test Decimals.setexponent(d"-3.3333", -3, RoundDown) ==  d"-3.333"
    @test Decimals.setexponent(d"3.3333", -3, RoundFloor) ==  d"3.333"
    @test Decimals.setexponent(d"-3.3333", -3, RoundFloor) ==  d"-3.334"
    @test Decimals.setexponent(d"3.3333", -3, RoundCeiling) ==  d"3.334"
    @test Decimals.setexponent(d"-3.3333", -3, RoundCeiling) ==  d"-3.333"
end


text:
using FlatBuffers
using Test

#@testset "PR 234 Arrow.jl: bugfix parsing primitive arrays" begin
#	buf = [
#	    0x14,0x00,0x00,0x00,0x00,0x00,0x0e,0x00,0x14,0x00,0x00,0x00,0x10,0x00,0x0c,0x00,0x08,
#	    0x00,0x04,0x00,0x0e,0x00,0x00,0x00,0x2c,0x00,0x00,0x00,0x38,0x00,0x00,0x00,0x38,0x00,
#	    0x00,0x00,0x38,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
#	    0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
#	    0x03,0x00,0x00,0x00,0x01,0x00,0x00,0x00,0x02,0x00,0x00,0x00,0x03,0x00,0x00,0x00,0x00,
#	    0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
#	]
#
#	struct TestData <: FlatBuffers.Table
#	    bytes::Vector{UInt8}
#	    pos::Base.Int
#	end
#
#	function Base.getproperty(x::TestData, field::Symbol)
#	    if field === :DataInt32
#	        o = FlatBuffers.offset(x, 12)
#	        o != 0 && return FlatBuffers.Array{Int32}(x, o)
#	    else
#	        @warn "field $field not supported"
#	    end
#	end
#
#	d = FlatBuffers.getrootas(TestData, buf, 0);
#	@test d.DataInt32 == UInt32[1,2,3]
#end

@testset "monsters example" begin
	codegendir = "monster/__codegen/"
	if "FLATC" ∈ keys(ENV)
		mkpath(codegendir)
		flatc = ENV["FLATC"] # path to flatc executable based on https://github.com/jonalm/flatbuffers
		run(`$flatc -o $codegendir --julia monster/monster.fbs`)
		include(joinpath(codegendir, "monster_generated.jl"))

		mana = Int16(99)
    	hp = Int16(123)
    	x, y, z = 1f0, 2f0, 3f0
    	color = Foo.Color.Red
    	w1_name = "Sword"
    	w1_damage = Int16(3)
    	w2_name = "Axe"
    	w2_damage = Int16(321)
    	b = FlatBuffers.Builder()
    	weapon_one_name = FlatBuffers.createstring!(b, w1_name)
    	weapon_two_name = FlatBuffers.createstring!(b, w2_name)
    	pos = Foo.createVec3(b, x,y,z)
    	Foo.WeaponStart(b)
    	Foo.WeaponAddName(b, weapon_one_name)
    	Foo.WeaponAddDamage(b, w1_damage)
    	weapon_one = Foo.WeaponEnd(b)
    	Foo.WeaponStart(b)
    	Foo.WeaponAddName(b, weapon_two_name)
    	Foo.WeaponAddDamage(b, w2_damage)
    	weapon_two = Foo.WeaponEnd(b)
    	weapon_vec = [weapon_one, weapon_two]
    	Foo.MonsterStartWeaponsVector(b, length(weapon_vec))
    	for w in Iterators.reverse(weapon_vec)
    	    FlatBuffers.prependoffset!(b, w)
    	end
    	weapons = FlatBuffers.endvector!(b, length(weapon_vec))

		doubles = [1.0, 2.0, 3.0]
		Foo.MonsterStartArrayofdoublesVector(b, length(doubles))
    	for w in Iterators.reverse(doubles)
    	    FlatBuffers.prependoffset!(b, w)
    	end
    	doubles_ = FlatBuffers.endvector!(b, length(doubles))

		shorts = Int16[1, 2, 3]
		Foo.MonsterStartArrayofshortsVector(b, length(shorts))
    	for w in Iterators.reverse(shorts)
    	    FlatBuffers.prependoffset!(b, w)
    	end
    	shorts_ = FlatBuffers.endvector!(b, length(shorts))

		ints = Int32[1, 2, 3]
		Foo.MonsterStartArrayofintsVector(b, length(ints))
    	for w in Iterators.reverse(ints)
    	    FlatBuffers.prependoffset!(b, w)
    	end
    	ints_ = FlatBuffers.endvector!(b, length(ints))

		floats = Float32[1.0, 2.0, 3.0]
		Foo.MonsterStartArrayoffloatsVector(b, length(floats))
    	for w in Iterators.reverse(floats)
    	    FlatBuffers.prependoffset!(b, w)
    	end
    	floats_ = FlatBuffers.endvector!(b, length(floats))

    	Foo.MonsterStart(b)
    	Foo.MonsterAddPos(b, pos)
    	Foo.MonsterAddMana(b, mana)
		Foo.MonsterAddArrayofshorts(b, shorts_)
		Foo.MonsterAddArrayofints(b, ints_)
		Foo.MonsterAddMydouble(b, 1.234)
		Foo.MonsterAddArrayofdoubles(b, doubles_)
		Foo.MonsterAddMyfloat(b, 2.34f0)
		Foo.MonsterAddArrayoffloats(b, floats_)
    	Foo.MonsterAddColor(b, color)
    	Foo.MonsterAddHp(b, hp)
    	Foo.MonsterAddWeapons(b,  weapons)
    	monster = Foo.MonsterEnd(b)
    	FlatBuffers.finish!(b, monster)
	    monsterbuf = FlatBuffers.finishedbytes(b)
	    buf = Vector{UInt8}(monsterbuf) # is this the best way of reading back?
	    monster_ = FlatBuffers.getrootas(Foo.Monster, buf, 0)
	    @test monster_.mana === mana
	    @test monster_.hp === hp
	    @test monster_.pos.x === x
	    @test monster_.pos.y === y
	    @test monster_.pos.z === z
	    @test monster_.color === color
	    @test monster_.weapons[1].name === w1_name
	    @test monster_.weapons[2].name === w2_name
	    @test monster_.weapons[1].damage === w1_damage
	    @test monster_.weapons[2].damage === w2_damage

		@test monster_.mydouble ===  1.234
		@test monster_.arrayofdoubles[1] ===  doubles_[1]
		@test monster_[Foo.MonsterProperties.arrayofdoubles][1] ===  doubles_[1]

		@test monster_.myfloat ===  2.34f0
		@test monster_.arrayoffloats[1] ===  floats_[1]
		@show monster_.arrayoffloats[1]
		@show monster_.arrayoffloats[2]
		@show monster_.arrayoffloats[3]

		@show monster_.arrayofints[1]
		@show monster_.arrayofints[2]
		@show monster_.arrayofints[3]

		@show monster_.arrayofshorts[1]
		@show monster_.arrayofshorts[2]
		@show monster_.arrayofshorts[3]

		#@test monster_[Foo.MonsterProperties.arrayofdoubles][1] ===  doubles_[1]

	else
		@info "Didn't find flatc executable in as `ENV[\"FLATC\"]`"
	end
end


text:
#Partition a dataset into K clusters.

#Finds clusters by repeatedly assigning each data point to the cluster with
#the nearest centroid and iterating until the assignments converge (meaning
#they don't change during an iteration) or the maximum number of iterations
#is reached.

#Parameters
#----------
#K : int
#    The number of clusters into which the dataset is partitioned.
#max_iters: int
#    The maximum iterations of assigning points to the nearest cluster.
#    Short-circuited by the assignments converging on their own.
#init: str, default 'random'
#    The name of the method used to initialize the first clustering.
#
#   'random' - Randomly select values from the dataset as the K centroids.
#    '++' - Select a random first centroid from the dataset, then select
#           K - 1 more centroids by choosing values from the dataset with a
#           probability distribution proportional to the squared distance
#           from each point's closest existing cluster. Attempts to create
#           larger distances between initial clusters to improve convergence
#           rates and avoid degenerate cases.

type Kmeans
    clusters::Dict{Integer,Matrix}
    clu_ind::Dict{Integer,Vector}
    k::Integer
    X::Matrix
    max_iter::Integer
    centroid::Matrix
    init::String
end

function Kmeans(;
                clusters::Dict{Integer,Matrix} = Dict{Integer,Matrix}(),
                k::Integer = 2,
                X::Matrix = zeros(2,2),
                max_iter::Integer = 150,
                centroid::Matrix = zeros(2,2),
                init::String = "++",
                clu_ind::Dict{Integer,Vector} = Dict{Integer,Vector}())

    return Kmeans(clusters, clu_ind, k, X, max_iter, centroid , init)
end

function train!(model::Kmeans, X::Matrix)
    model.X = X
end

function predict!(model::Kmeans)

    initialize_centroid(model)
    for i in 1:model.max_iter
        centroid_old = copy(model.centroid)
        assign_clusters!(model)
        update_centroid!(model)
        if is_converged(centroid_old, model.centroid)
            break
        end
    end
end

function update_centroid!(model)
    for i = 1:model.k
        if !haskey(model.clusters, i)
            r = 0
            for i = 1:model.k
                if haskey(model.clusters, i)
                    r = i
                end
            end
            model.clusters[i] = model.clusters[r]
        end
        model.centroid[i,:] = mean(model.clusters[i],1)
    end
end

function assign_clusters!(model::Kmeans)
    n = size(model.X,1)
    model.clusters = Dict{Integer,Matrix}()
    model.clu_ind = Dict{Integer,Vector}()
    for i = 1:n
        dist = zeros(model.k)
        for j = 1:length(dist)
            dist[j] = norm(model.centroid[j,:]-model.X[i,:])
        end
        clu = indmin(dist)
        if haskey(model.clusters,clu)
            model.clusters[clu] = vcat(model.clusters[clu], model.X[i,:]')
            model.clu_ind[clu] = vcat(model.clu_ind[clu], i)
        else
            model.clusters[clu] = model.X[i,:]'
            model.clu_ind[clu] = [i]
        end

    end
end

function initialize_centroid(model::Kmeans)
    model.centroid = zeros(model.k, size(model.X,2))
    if model.init == "random"
        model.centroid = model.X[randperm(size(model.X,1))[1:model.k],:]
    elseif model.init == "++"
        model.centroid[1,:] = model.X[rand(1:size(model.X,1)),:]
        for i = 2:model.k
            model.centroid[i,:] = find_next_centroid(model,i-1)
        end
    else
        error("you must provide a initial type")
    end

end

function find_next_centroid(model, num)
    mean_cent = vec(mean(model.centroid[1:num,:],1))
    n_sample = size(model.X,1)
    res = zeros(n_sample)
    for i = 1:n_sample
        res[i] = norm(model.X[i,:]-mean_cent)
    end
    prob = res/sum(res)
    cum = cumsum(prob, 1)
    r = rand()
    x_sel = model.X[cum .> r,:]
    x_sel = x_sel[1,:]
    return x_sel
end

function is_converged(x::Matrix,
                      y::Matrix)
    return norm(x-y) == 0
end

function plot_in_2d(model::Kmeans)
    y_ = []
    for i = 1:model.k
        push!(y_,model.clu_ind[i])
    end
    y_ = reduce(vcat,y_)
    len = size(model.X, 2)
    x_ = zeros(size(y_,1),len)
    for i = 1:size(x_,1)
        x_[i,:] = model.X[y_[i],:]
    end
    # if not 2d perform PCA

    #PCA
    if size(x_, 2) > 2
        pca_model = PCA()
        train!(pca_model, x_)
        x_ = transform(pca_model, x_)
        x_ = x_[:, 1:2]
    end

    x_sep = x_[:,1]
    y_sep = x_[:,2]
    num_ = zeros(model.k)
    for i = 1:model.k
        num_[i] = size(model.clusters[i],1)
    end
    num_ = convert(Array{Int,1},num_)
    rep = zeros(sum(num_))
    for i = 1:model.k
        for j = 1:num_[i]
            j = convert(Int,j)
            if i == 1
                rep[j] = i
            else
                kk = i-1
                rep[sum(num_[1:kk])+j] = i
            end
        end
    end
    df = DataFrame(x = x_sep,y = y_sep , cluster = rep)
    println("Computing finished")
    println("Drawing the plot.....Please Wait(Actually Gadfly is quite slow in drawing the first plot)")
    Gadfly.plot(df, x = "x", y = "y", color = "cluster",Geom.point)
end

function test_kmeans_speed()
    X, y = make_blo()
    clu = length(unique(y))
    @show clu
    model = Kmeans(k=clu)
    train!(model,X)
    predict!(model)
    plot_in_2d(model)
end

function test_kmeans_random()
    X, y = make_blo()
    clu = length(unique(y))
    @show clu
    model = Kmeans(k=clu,init="random")
    train!(model,X)
    predict!(model)
    plot_in_2d(model)
end


text:
@inline function Base.checkbounds(a::AbstractArray, i::T, j::T) where T <: NucleicAcid
  if iscertain(i) & iscertain(j)
    return true
  end
  checkbounds(a, trailing_zeros(i) + 1, trailing_zeros(j) + 1)
end

@inline function Base.checkbounds(a::AbstractArray, i::NucleicAcid)
  if iscertain(i)
      return true
  end
  checkbounds(a, trailing_zeros(i) + 1)
end

@inline function Base.getindex(a::AbstractArray, i::T, j::T) where T <: NucleicAcid
  @boundscheck checkbounds(a, i, j)
  @inbounds return a[trailing_zeros(i) + 1, trailing_zeros(j) + 1]
end

@inline function Base.getindex(a::AbstractArray, i::NucleicAcid)
  @boundscheck checkbounds(a, i)
  @inbounds return a[trailing_zeros(i) + 1]
end

@inline function Base.setindex!(a::AbstractArray, x, i::T, j::T) where T <: NucleicAcid
  @boundscheck checkbounds(a, i, j)
  @inbounds return setindex!(a, x, trailing_zeros(i) + 1, trailing_zeros(j) + 1)
end

@inline function Base.setindex!(a::AbstractArray, x, i::NucleicAcid)
  @boundscheck checkbounds(a, i)
  @inbounds return setindex!(a, x, trailing_zeros(i) + 1)
end


text:
using JSON
using Requests
import Requests: URI, get

using DandelionWebSockets
# Explicitly import the callback functions that we're going to add more methods for.
import DandelionWebSockets: on_text, on_binary,
                            state_connecting, state_open, state_closing, state_closed

type BitcoinBlockchainHandler <: WebSocketHandler
  client::WSClient
  stop_channel::Channel{Any}
  debug::Bool
  op_callbacks::Dict
  max_attempts::Int

  BitcoinBlockchainHandler() = new(WSClient(), Channel{Any}(3), false, Dict(), 6)

end

function getResult(handler, url)
  attempt = 1
  resp = get(url)

  # Retry if service temporarily not available
  while resp.status == 503 && attempt < handler.max_attempts
    sleep(attempt) #increase wait time by one second for each failed attempt
    attempt += 1
    resp = get(url)
  end
  if resp.status == 503
    error("Aborting request with url: $url after $(handler.max_attempts) attempts")
  elseif resp.status != 200
    error("$(resp.status): Error executing the request with url: $url, resp: $resp")
  else
    try
      parsedresp = Requests.json(resp)
      if "error" in keys(parsedresp)
        error("Error parsing response to request with url: $url, resp: $(parsedresp["error"])")
      end
      parsedresp
    catch e
      error("Error parsing response to request with url: $url, resp: $resp")
    end
  end
end

# These are called when you get text/binary frames, respectively.
on_text(handler::BitcoinBlockchainHandler, s::String)         = onMessage(handler, s, false)
on_binary(handler::BitcoinBlockchainHandler, data::Vector{UInt8}) = onMessage(handler, data, true)

function onMessage(handler::BitcoinBlockchainHandler, payload, isBinary)
  if handler.debug
    isBinary && println("Binary message received: $(length(payload)) bytes.")
    !isBinary && println("Text message received: $payload.")
  end

  try

    responseEventInfo = JSON.parse(payload)

    # ping response
    responseEventInfo["op"] == "pong" && return
    handle_command(handler, responseEventInfo)

  catch e
    error("Error processing payload: $payload: $e")
  end
end

function sendMessage(handler::BitcoinBlockchainHandler, payload, isBinary)
  handler.debug && println("Sending: $payload")
  isBinary && send_binary(handler.client, payload)
  !isBinary && send_text(handler.client, String(payload))
end

# These are called when the WebSocket state changes.

state_connecting(::BitcoinBlockchainHandler) = println("State: CONNECTING")

# Called when the connection is open, and ready to send/receive messages.
function state_open(handler::BitcoinBlockchainHandler)
  println("State: OPEN")
  ping(handler)
end

state_closing(handler::BitcoinBlockchainHandler) = println("State: CLOSING")
state_closed(handler::BitcoinBlockchainHandler) = println("State: CLOSED")

function handle_command(handler::BitcoinBlockchainHandler, eventInfo)
  handler.debug && println("HANDLE COMMAND")
  # Server told us to do something.
  op = eventInfo["op"]
  if haskey(handler.op_callbacks, op)
    cb = handler.op_callbacks[op]
    cb(eventInfo)
  end
end

function publish(data)
  println("Publish: $(JSON.json(data))")
end

function subscribe(handler, op, publish_key, publish_callback)
  sendMessage(handler, json(Dict("op" => op)), false)
  handler.op_callbacks[publish_key] = publish_callback
  handler.debug && println("Subscribed $op")
end

function unsubscribe(handler, op, publish_key)
  delete!(handler.op_callbacks, publish_key)
  handler.debug && println("Unsubscribed $op")
end

function ping(handler)
  @schedule begin
    while true
      sleep(30)
      sendMessage(handler, json(Dict("op" => "ping")), false)
    end
  end
end

function connect(handler::BitcoinBlockchainHandler)
  uri = URI("wss://ws.blockchain.info/inv")
  println("Connecting to $uri")
  wsconnect(handler.client, uri, handler)
end

subscribe_unconfirmed_transcations(handler, publish_callback = subscribe) = subscribe(handler, "unconfirmed_sub", "utx", publish_callback)
subscribe_new_blocks(handler, publish_callback = subscribe) = subscribe(handler, "blocks_sub", "block", publish_callback)
get_block(handler, block) = getResult(handler, "https://blockchain.info/rawblock/$block")
get_transaction(handler, tx) = getResult(handler, "https://blockchain.info/rawtx/$tx")
get_address(handler, address) = getResult(handler, "https://blockchain.info/rawaddr/$address")


text:
export # =====================PML absorbing boundaries============
       dampCoef,
       InitDampCoef,
       DampBound,
       modExpand,
       modelPadding,
       modSmooth,
       #  =====================finite difference matrix===========
       FidMtx,
       vel2lm,
       lm2vel,
       DispStable!,
       MVxxBVxx,
       MVxxBTxx,
       MVxzBVxz,
       MVxzBTxz,
       MVzxBVzx,
       MVzxBTxz,
       MVzzBVzz,
       MVzzBTzz,
       MTxxxBVx,
       MTxxxBTxxx,
       MTxxzBVz,
       MTxxzBTxxz,
       MTzzxBVx,
       MTzzxBTzzx,
       MTzzzBVz,
       MTzzzBTzzz,
       MTxzxBVz,
       MTxzxBTxzx,
       MTxzzBVx,
       MTxzzBTxzz,
       CreateFidMtx,
      #  =====================SnapShots and WFD=====================
       SnapShot,
       InitSnapShot,
       Wfd,
       InitWfd,
       CopySnapShot,
       CopySnapShot!,
       Add2SnapShots!,
       WriteSnapShot,
       ReverseOrderSnapShots,
       WriteWfd,
       extractBoundary!,
       WriteWfdBoundary,
       ReadWfdBoundary,
       ReadWfd,
       InfoWfd,
       InfoSnapShots,
       ReadSnapShot,
       ExtractCpt,
       NormOfSnapShots,
       NormOfSpt,
       ScaleWfd!,
       ReverseOrderWfd,
       IptWfd,
       partialV!,
       WritePv,
       ReadPv,
       InfoPv,
       checkPv,
       IptSpts,
       dif2spt,
      #  =====================Sources=================================
       Source,
       Ricker,
       InitSource,
       InitMultiSources,
       SrcRange,
       SourcesTimeRange,
       AddSource!,
       AddSourceBorn!,
       AddMultiSources!,
       srcs2spt,
       ConvertBornSource2Spts,
      #  =====================Records=================================
       Shot,
       InitShot,
       ShotV,
       InitShotV,
       WriteShot,
       WriteShotV,
       ReadShot,
       ReadShotV,
       InfoShotV,
       Spt2Shot!,
       Spt2ShotV!,
       AddShot2Spt!,
       AddShotV2Spt!,
       NormOfShot,
       imaging,
       SptsPv2Born,
       removeDirect,
       #  =====================Image=================================
       WriteImage,
       ReadImage,
       laplaceFilter,
       lambdamu2vpvs,
       gec,
      #  ======================SrcLoc================================
       WSrcCube,
       InitWSrcCube,
       CopyWsc,
       l112norm,
       EngDisWsc,
       Srcs2Wsc,
       AddWsc2spt!,
       spt2Wsc!

include("dampCoef.jl")
include("fidMtx.jl")
include("SnapShot.jl")
include("Source.jl")
include("Record.jl")
include("Image.jl")
include("WsrcCube.jl")

max_stars_repo_name: WneleiGao/ElasticWave


text:
#=
  Original Author: Anastasia Yendiki

  Copyright © 2022 The General Hospital Corporation (Boston, MA) "MGH"

  Terms and conditions for use, reproduction, distribution and contribution
  are found in the 'FreeSurfer Software License Agreement' contained
  in the file 'LICENSE' found in the FreeSurfer distribution, and here:

  https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense

  Reporting: freesurfer@nmr.mgh.harvard.edu
=#

using LinearAlgebra

export Tract, trk_read, trk_write

"Container for header and streamline data stored in .trk format"
mutable struct Tract
  # Header fields (.trk format version 2)
  id_string::Vector{UInt8}
  dim::Vector{Int16}
  voxel_size::Vector{Float32}
  origin::Vector{Float32}
  n_scalars::Int16
  scalar_name::Matrix{UInt8}
  n_properties::Int16
  property_name::Matrix{UInt8}
  vox_to_ras::Matrix{Float32}
  reserved::Vector{UInt8}
  voxel_order::Vector{UInt8}
  voxel_order_original::Vector{UInt8}
  image_orientation_patient::Vector{Float32}
  pad1::Vector{UInt8}
  invert_x::UInt8
  invert_y::UInt8
  invert_z::UInt8
  swap_xy::UInt8
  swap_yz::UInt8
  swap_zx::UInt8
  n_count::Int32
  version::Int32
  hdr_size::Int32

  # Streamline data fields
  npts::Vector{Int32}
  properties::Matrix{Float32}
  xyz::Vector{Any}
  scalars::Vector{Any}
end

"""
    Tract()

Return an empty `Tract` structure
"""
Tract() = Tract(
  Vector{UInt8}(undef, 0),
  Vector{Int16}(undef, 0),
  Vector{Float32}(undef, 0),
  Vector{Float32}(undef, 0),
  Int16(0),
  Matrix{UInt8}(undef, 0, 0),
  Int16(0),
  Matrix{UInt8}(undef, 0, 0),
  Matrix{Float32}(undef, 0, 0),
  Vector{UInt8}(undef, 0),
  Vector{UInt8}(undef, 0),
  Vector{UInt8}(undef, 0),
  Vector{Float32}(undef, 0),
  Vector{UInt8}(undef, 0),
  UInt8(0),
  UInt8(0),
  UInt8(0),
  UInt8(0),
  UInt8(0),
  UInt8(0),
  Int32(0),
  Int32(0),
  Int32(0),

  Vector{Int32}(undef, 0),
  Matrix{Float32}(undef, 0, 0),
  [],
  []
)

"""
    Tract(ref::MRI)

Return a `Tract` structure whose header fields are populated based on the
reference `MRI` structure `ref`
"""
function Tract(ref::MRI)

  tr = Tract()

  # In the following I must take into account that mri_read()
  # may have permuted the first 2 dimensions of volsize and volres,
  # but it has not changed the vox2ras matrix

  # Find orientation of image coordinate system
  orient = vox2ras_to_orient(ref.vox2ras)

  # Find patient-to-scanner coordinate transform:
  # Take x and y vectors from vox2RAS matrix, convert to LPS,
  # divide by voxel size
  if ref.ispermuted
    p2s = [-1 0 0; 0 -1 0; 0 0 1] * ref.vox2ras[1:3, 1:2] *
                                    Diagonal([1, 1]./ref.volres[[2,1]])
  else
    p2s = [-1 0 0; 0 -1 0; 0 0 1] * ref.vox2ras[1:3, 1:2] *
                                    Diagonal([1, 1]./ref.volres[1:2])
  end

  tr.id_string     = UInt8.(collect("TRACK\0"))

  if ref.ispermuted
    tr.dim         = Int16.(ref.volsize[[2,1,3]])
    tr.voxel_size  = Float32.(ref.volres[[2,1,3]])
  else
    tr.dim         = Int16.(ref.volsize)
    tr.voxel_size  = Float32.(ref.volres)
  end
  tr.origin        = Float32.(zeros(3))

  tr.n_scalars     = Int16(0)
  tr.scalar_name   = UInt8.(zeros(10, 20))
  tr.n_properties  = Int16(0)
  tr.property_name = UInt8.(zeros(10, 20))

  tr.vox_to_ras    = ref.vox2ras
  tr.reserved      = UInt8.(zeros(444))
  tr.voxel_order               = vcat(UInt8.(collect(orient)), UInt8(0))
  tr.voxel_order_original      = tr.voxel_order
  tr.image_orientation_patient = Float32.(p2s[:])
  tr.pad1          = UInt8.(zeros(2))

  tr.invert_x      = UInt8(0)
  tr.invert_y      = UInt8(0)
  tr.invert_z      = UInt8(0)
  tr.swap_xy       = UInt8(0)
  tr.swap_yz       = UInt8(0)
  tr.swap_zx       = UInt8(0)
  tr.n_count       = Int32(0)
  tr.version       = Int32(2)
  tr.hdr_size      = Int32(1000)

  return tr
end

"""
    trk_read(infile::String)

Read tractography streamlines from `infile` and return a `Tract` structure

Input file must be in .trk format, see:
http://trackvis.org/docs/?subsect=fileformat
"""
function trk_read(infile::String)

  tr = Tract()

  io = open(infile, "r")

  # Read .trk header
  tr.id_string     = read!(io, Vector{UInt8}(undef, 6))

  tr.dim           = read!(io, Vector{Int16}(undef, 3))
  tr.voxel_size    = read!(io, Vector{Float32}(undef, 3))
  tr.origin        = read!(io, Vector{Float32}(undef, 3))

  tr.n_scalars     = read(io, Int16)
  tr.scalar_name   = read!(io, Matrix{UInt8}(undef, 20, 10))
  tr.scalar_name   = permutedims(tr.scalar_name, [2,1])

  tr.n_properties  = read(io, Int16)
  tr.property_name = read!(io, Matrix{UInt8}(undef, 20, 10))
  tr.property_name = permutedims(tr.property_name, [2,1])

  tr.vox_to_ras    = read!(io, Matrix{Float32}(undef, 4, 4))
  tr.vox_to_ras    = permutedims(tr.vox_to_ras, [2,1])
  tr.reserved      = read!(io, Vector{UInt8}(undef, 444))
  tr.voxel_order               = read!(io, Vector{UInt8}(undef, 4))
  tr.voxel_order_original      = read!(io, Vector{UInt8}(undef, 4))
  tr.image_orientation_patient = read!(io, Vector{Float32}(undef, 6))
  tr.pad1          = read!(io, Vector{UInt8}(undef, 2))

  tr.invert_x      = read(io, UInt8)
  tr.invert_y      = read(io, UInt8)
  tr.invert_z      = read(io, UInt8)
  tr.swap_xy       = read(io, UInt8)
  tr.swap_yz       = read(io, UInt8)
  tr.swap_zx       = read(io, UInt8)
  tr.n_count       = read(io, Int32)
  tr.version       = read(io, Int32)
  tr.hdr_size      = read(io, Int32)

  # Read streamline data
  tr.npts       = Vector{Int32}(undef, tr.n_count)
  tr.properties = Matrix{Float32}(undef, tr.n_properties, tr.n_count)

  for istr in 1:tr.n_count		# Loop over streamlines
    tr.npts[istr] = read(io, Int32)

    push!(tr.xyz, Matrix{Float32}(undef, 3, tr.npts[istr]))
    push!(tr.scalars, Matrix{Float32}(undef, tr.n_scalars, tr.npts[istr]))

    for ipt in 1:tr.npts[istr]		# Loop over points on a streamline
      # Divide by voxel size and make voxel coordinates 0-based
      tr.xyz[istr][:, ipt] = read!(io, Vector{Float32}(undef, 3)) ./
                             tr.voxel_size .- .5

      tr.scalars[istr][:, ipt] = read!(io, Vector{Float32}(undef, tr.n_scalars))
    end

    tr.properties[:, istr] = read!(io, Vector{Float32}(undef, tr.n_properties))
  end

  close(io)

  return tr
end

"""
     trk_write(tr::Tract, outfile::String)

Write a `Tract` structure to a file in the .trk format

Return true if an error occurred (i.e., the number of bytes written was not the
expected based on the size of the `Tract` structure)
"""
function trk_write(tr::Tract, outfile::String)

  io = open(outfile, "w")

  nb = 0

  # Write .trk header
  nb += write(io, UInt8.(tr.id_string))

  nb += write(io, Int16.(tr.dim))
  nb += write(io, Float32.(tr.voxel_size))
  nb += write(io, Float32.(tr.origin))

  nb += write(io, Int16(tr.n_scalars))
  nb += write(io, UInt8.(permutedims(tr.scalar_name, [2,1])))
  nb += write(io, Int16(tr.n_properties))
  nb += write(io, UInt8.(permutedims(tr.property_name, [2,1])))

  nb += write(io, Float32.(permutedims(tr.vox_to_ras, [2,1])))
  nb += write(io, UInt8.(tr.reserved))
  nb += write(io, UInt8.(tr.voxel_order))
  nb += write(io, UInt8.(tr.voxel_order_original))
  nb += write(io, Float32.(tr.image_orientation_patient))
  nb += write(io, UInt8.(tr.pad1))

  nb += write(io, UInt8(tr.invert_x))
  nb += write(io, UInt8(tr.invert_y))
  nb += write(io, UInt8(tr.invert_z))
  nb += write(io, UInt8(tr.swap_xy))
  nb += write(io, UInt8(tr.swap_yz))
  nb += write(io, UInt8(tr.swap_zx))
  nb += write(io, Int32(tr.n_count))
  nb += write(io, Int32(tr.version))
  nb += write(io, Int32(tr.hdr_size))

  # Write streamline data
  for istr in 1:tr.n_count		# Loop over streamlines
    nb += write(io, Int32(tr.npts[istr]))

    for ipt in 1:tr.npts[istr]		# Loop over points on a streamline
      # Make voxel coordinates .5-based and multiply by voxel size
      nb += write(io, Float32.((tr.xyz[istr][:, ipt] .+ .5) .* tr.voxel_size))

      nb += write(io, Float32.(tr.scalars[istr][:, ipt]))
    end

    nb += write(io, Float32.(tr.properties[:, istr]))
  end

  close(io)

  err = (nb != sizeof(UInt8) * 866 +
               sizeof(Int32) * (3 + length(tr.npts)) +
               sizeof(Int16) * 5 +
               sizeof(Float32) * (28 + sum(length.(tr.xyz)) +
                                       sum(length.(tr.scalars)) +
                                       length(tr.properties)))

  return err
end

max_stars_repo_name: freesurfer/jupysurf

max_issues_repo_name: freesurfer/jupysurf

max_forks_repo_name: freesurfer/jupysurf


text:
# 3.4 plotting
tmpexpr = :(
   xlabel("Year"); grid(true);
   xlim([ idx_year2plot[1] - 1, idx_year2plot[end] + 1 ]);
)
figure(figsize = (13,8))
   subplot(2,2,1)  # gap / exp
      for tmpreform in 1:length(res_IndiInDesc["reformNames"])
         plot( Dt[:Year][idx_plot], res_IndiInDesc["gap/exp"][tmpreform], res_IndiInDesc["plotstyles"][tmpreform] )
      end
      eval(tmpexpr); ylabel("Gap / Pool Expenditure (%)");
      legend( res_IndiInDesc["plotlegends"] , loc = "lower right")
   subplot(2,2,2)  # gap / in
      for tmpreform in 1:length(res_IndiInDesc["reformNames"])
         plot( Dt[:Year][idx_plot], res_IndiInDesc["gap/in"][tmpreform], res_IndiInDesc["plotstyles"][tmpreform] )
      end
      eval(tmpexpr); ylabel("Gap / Pool Income (%)");
      # legend( string.("ζ = ",reform_zeta_levs) , loc = "upper left")
   subplot(2,2,3)  # gap / gdp
      for tmpreform in 1:length(res_IndiInDesc["reformNames"])
         plot( Dt[:Year][idx_plot], res_IndiInDesc["gap/gdp"][tmpreform], res_IndiInDesc["plotstyles"][tmpreform] )
      end
      eval(tmpexpr); ylabel("Gap / GDP (%)");
      # legend( string.("ζ = ",reform_zeta_levs) , loc = "upper left")
   subplot(2,2,4)  # gap / taxrev
      for tmpreform in 1:length(res_IndiInDesc["reformNames"])
         plot( Dt[:Year][idx_plot], res_IndiInDesc["gap/taxrev"][tmpreform], res_IndiInDesc["plotstyles"][tmpreform] )
      end
      eval(tmpexpr); ylabel("Gap / Tax Revenues (%)");
      # legend( string.("ζ = ",reform_zeta_levs) , loc = "upper left")
   tight_layout()
# 3.5 save figure
savefig( "./output/IndiAcctReform.pdf", format = "pdf" )

#


text:
include("Circuits.jl")
using Catlab.WiringDiagrams
using Catlab.CategoricalAlgebra
using Catlab.Theories
using Catlab.Syntax
using Catlab.GAT
using Catlab.Present
using Catlab.Programs
using Catlab.Graphics
using Catlab.Graphics.Graphviz
using Catlab.Programs
using Catlab.WiringDiagrams
using MLStyle
#using ..Compologicircuits: CircuitDom, Circuit, iNOT, iAND, iOR, Impl, Circuits, tfba_expr, show_diagram
#using PicoSAT
import Catlab.Theories: id, compose, otimes, braid, munit, mcopy, delete, pair, proj1, proj2, ⊗, ⋅, σ

B = Ob(FreeCartesianCategory, :B)
dup = mcopy(B)
I = id(B)
and = Hom(:AND, B⊗B, B)
or = Hom(:OR, B⊗B, B)
xor = Hom(:XOR, B⊗B,B)
not = Hom(:NOT, B,B)

to_diagram(x) = add_junctions!(to_wiring_diagram(x))

normal_form(x) = to_hom_expr(FreeCartesianCategory, to_diagram(x))

#returns set of all identical pairs of n length bool vectors
function all_match(n::Int)
    output = []
    for i in 0:(2^n)-1
        vect = .!Vector{Bool}(iszero.(digits(i, base=2, pad=n)))
        push!(output, [vect, vect])
    end
    return output
end

#binary vectors to int
function bivec_int(w)::Int
    output = 0
    v = Vector{Int}(w)

    for i in 0:length(v)-1
        output += (v[i+1])*2^i
    end

    return output
end

#brute force SAT algorithm, loops over all the possible inputs
#return true and first argument that gives true if sat., else returns false and the last argument
function brute_SAT(f::Circuit)
    sat = false
    i = 0
    dom = f.dom.n
    arr = Bool[]

    if !(f.codom == CircuitDom(1))
        error("circuit codomain needs to be 1!")
    else
        while (!sat)&&(i < 2^dom)
                arr = iszero.(digits(i, base=2, pad=dom))
                sat = f.impl(arr)[1]
                i += 1
        end
        return (sat, arr)
    end
end

brute_SAT(expr) = brute_SAT(Impl(expr))

#if dup is never used then the circtuit has all possible outputs, hence satisfiable
#returns true if contains dup, otherwise return false
function dup_check(expr)::Bool
    argss = args(expr)
    output = false
    headd = head(expr)

    if headd == :mcopy
        output = true
    elseif (headd == :generator)||(headd == :id)||(headd == :braid)||(headd == :pair)||
           (headd == :delete)||(headd == :proj1)||(headd == :proj2)
        output = false
    else
        for i in argss
            output = output||dup_check(i)
        end
    end

    return output
end

vibe_check(expr::FreeCartesianCategory.Hom)::Bool =
    @match head(expr) begin
        :mcopy => true
        :compose || :otimes => any(map(vibe_check, args(expr)))
        _ => false
    end

#replacing all instances of i and j in v with r, ignoring but keeping their signs
function copy_replace(v::Vector, i::Int, j::Int, r::Int)
    count = 1
    for a in v
        if (abs(a) == abs(i)) || (abs(a) == abs(j))
            v[count] = Int((abs(a)/a)*abs(r))
        end
        count += 1
    end
    return v
end

#replace all instances of n with m, and m with n, at the same time, ignoring but keeping signs
function braid_replace(v::Vector, n::Int, m::Int)
    count = 1
    for i in v
        if abs(i) == abs(n)
            v[count] = Int((abs(i)/i)*abs(m))
        elseif abs(i) == abs(m)
            v[count] = Int((abs(i)/i)*abs(n))
        end
        count += 1
    end
    return v
end

#replace all instances of n or -n with its negation
function not_replace(v::Vector, n::Int)
    count = 1
    for i in v
        if abs(i) == abs(n)
            v[count] = -Int((abs(i)/i)*abs(n))
        end
        count += 1
    end
    return v
end

function populate(v::Vector, anc::Int, offs::Int)
    output = []
    for i in v
        if i != anc
            push!(output, i)
        else
            push!(push!(output, i), offs)
        end
    end
    return output
end

#get rid of all instances of n in v, except the first
function cleanse(v::Vector, n::Int)
    output = []
    found = false
    for i in v
        if i != n
            push!(output, i)
        elseif !found
            found = true
            push!(output, i)
        end
    end
    return output
end

function splate(v::Vector, old::Int, new::Int)
    temp =[]
    for i in v
        if i != old
            push!(temp, i)
        else
            push!(temp, new)
        end
    end
    if v == temp
        return [v]
    else
        return [v, temp]
    end
end

#apply the expression expr to the current cnf
#only works for expressions with no nested otimes, every exression is expressible without nested otimes
function app_expr(expr, vars::Vector, cnf::Vector, begn::Int, ed::Int)
    #assuming that if head(expr) == :otimes then no argument of expr has :otimes head
    #vars = the boolean variables on each wire in order, must all be positive
    #begn = which wire to begin applying expr in vars
    #ed = the number of input wires of the previous application, might not always be relevant

    out_vars = vars #the vars on the wires after applying expr
    out_cnf = cnf #cnf after applying expr
    out_begn = begn #which wire to start the next application, should always be 1
    out_ed = ed #the input wires of expr, might not always be relevant
    new_var = maximum(vars)+1 #only needed to introduce a new variable for and, or gates
    heads = head(expr)
    argss = args(expr)

    if heads == :compose
        for i in reverse(argss)
            (out_vars, out_cnf, out_begn, out_ed) = app_expr(i, out_vars, out_cnf, out_begn, out_ed)
        end
    elseif heads == :otimes
        #@assert begn == 1
        (out_vars, out_cnf, out_begn, out_ed) = app_expr(argss[1], out_vars, out_cnf, out_begn, out_ed)
        for i in argss[2:end]
            v = app_expr(i, out_vars, out_cnf, out_ed+1, out_ed)
            out_ed += v[4]
            out_vars = v[1]
            out_cnf = v[2]
        end
    #assume to be Δ(B)
    elseif heads == :mcopy
        n = max(vars[begn], vars[begn+1])
        out_cnf = map(x -> copy_replace(x, vars[begn], vars[begn+1], n), out_cnf)
        out_vars[begn] = n
        out_vars = out_vars[1:end .!= begn+1]
        out_ed = 1
    #assumed to be σ(B,B)
    elseif heads == :braid
        n = vars[begn]
        m = vars[begn+1]
        out_vars[begn] = m
        out_vars[begn+1] = n
        out_cnf = map(x -> braid_replace(x, n, m), out_cnf)
        out_ed = 2
    elseif heads == :id
        out_ed = Impl(expr).dom.n
    elseif expr == not
        var = out_vars[begn]
        out_ed = 1
        out_cnf = map(x -> not_replace(x, var), out_cnf)
    elseif expr == and
        var = vars[begn]
        out_vars = vcat(vcat(vars[1:begn], [new_var]), vars[begn+1:end])
        temp = Vector{Int}[]
        out_cnf = map(x -> cleanse(x, var), out_cnf)
        out_cnf = map(x -> cleanse(x, -var), out_cnf)
        for v in out_cnf
            if (!(var in v) && !(-var in v))
                push!(temp, v)
            else
                temp_v = populate(v, -var, -new_var)
                temp_v = splate(temp_v, var, new_var)
                temp = vcat(temp, temp_v)
            end
        end
        out_cnf = temp
        out_ed = 2
    elseif expr == or
        var = vars[begn]
        out_vars = vcat(vcat(vars[1:begn], [new_var]), vars[begn+1:end])
        temp = Vector{Int}[]
        out_cnf = map(x -> cleanse(x, var), out_cnf)
        out_cnf = map(x -> cleanse(x, -var), out_cnf)
        for v in out_cnf
            if (!(var in v) && !(-var in v))
                push!(temp, v)
            else
                temp_v = populate(v, var, new_var)
                temp_v = splate(temp_v, -var, -new_var)
                temp = vcat(temp, temp_v)
            end
        end
        out_cnf = temp
        out_ed = 2
    end

    return (out_vars, out_cnf, out_begn, out_ed)
end

function expr_to_cnf(expr)
    @assert Impl(expr).codom.n == 1
    return app_expr(normal_form(expr), [1], [[1]], 1, 0)
end

demo1 = (dup⊗dup)⋅(not⊗σ(B,B)⊗not)⋅(and⊗and)⋅and
#=demo_cnf = expr_to_cnf(demo1)[2]

PicoSAT.solve(demo_cnf)
time brute_SAT(demo1)=#


text:
include("gumbel_softmax.jl")

using ReverseDiff

and(x, y) = Int.((x + y) .== 2)
or(x, y) = Int.((x + y) .>= 1)
xor(x, y) = Int.((x .* y) .== 0)

const OP = [and, or, xor]

X = hcat([[x1,x2] for x1 in (0,1), x2 in (0,1)]...)'

function f(ps, x1, x2)
    T = 0.3
    logits = log.(ps)
    w = gumbel_softmax(logits, T, true)
    return w' * [op(x1,x2) for op in OP]
end

function objective(ps, x1, x2, y_target)
    abs2.(y_target - f(ps, x1, x2))
end

function learn(gt::Function=xor, α=0.1, N=50)
    ps = ones(length(OP)) ./ length(OP)
    Y = [gt(X[i,1], X[i,2]) for i=1:size(X,1)]
    for i = 1:N
        @show ps
        for j = 1:4
            g, _, _, _ = ReverseDiff.gradient(objective, (ps, [X[j,1]], [X[j,2]], [Y[j]]))
            #@show g
            if !any(isnan.(g))
                ps -= α * g
            end
            ps = max.(min.(ps, 1.0), 1e-12)
            normalize!(ps,1)
        end
    end
    return indmax(ps)
end

function test(n_seeds=50)
    result = zeros(length(OP), length(OP)) #truth x result
    for (i, gt) in enumerate(OP)
        for s=1:n_seeds
            srand(s)
            imax = learn(gt)
            result[i, imax] += 1
        end
    end
    result
end

#Issues:
#1. y_hard: tell ReverseDiff to hard forward, soft back
#2. variable temperature: tell ReverseDiff that T is a constant
#3. ReverseDiff throws and error when batch_size > 1


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# ╔═╡ 63ba08cc-59a8-11eb-0a0f-27efac60d779
using Pkg, DrWatson

# ╔═╡ 6db218c6-59a8-11eb-2a8b-7107354cf590
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample, StanOptimize
	using StatisticalRethinking
end

# ╔═╡ 51fc19b8-59a8-11eb-2214-15aca59b807b
md" ## Clip-08-13-17s.jl"

# ╔═╡ 8aaa4bcc-59a8-11eb-2003-f1213b116565
begin
	df = CSV.read(sr_datadir("rugged.csv"), DataFrame)
	dropmissing!(df, :rgdppc_2000)
	dropmissing!(df, :rugged)
	df.log_gdp = log.(df[:, :rgdppc_2000])
	df.log_gdp_s = df.log_gdp / mean(df.log_gdp)
	df.rugged_s = df.rugged / maximum(df.rugged)
	df.cid = [df.cont_africa[i] == 1 ? 1 : 2 for i in 1:size(df, 1)]
	PRECIS(df[:, [:rgdppc_2000, :log_gdp, :log_gdp_s, :rugged, :rugged_s, :cid]])
end

# ╔═╡ e890c59c-6a68-11eb-0d33-4d52a21d0ccf
stan8_1 = "
data {
	int N;
	vector[N] G;
	vector[N] R;
}

parameters {
	real a;
	real b;
	real<lower=0> sigma;
}

transformed parameters {
	vector[N] mu;
	mu = a + b * (R - 0.125);
}

model {
	a ~ normal(1, 0.1);
	b ~ normal(0, 0.3);
	sigma ~ exponential(1);
	G ~ normal(mu, sigma);
}
generated quantities {
	vector[N] log_lik;
	for (i in 1:N)
		log_lik[i] = normal_lpdf(G[i] | mu[i], sigma);
}
";

# ╔═╡ 8e815d5c-6a66-11eb-21dc-99734c31f4e1
begin
	data = (N = size(df, 1), K = length(unique(df.cid)),
		G = df.log_gdp_s, R = df.rugged_s, cid=df.cid)
	m8_1s = SampleModel("m8.1s", stan8_1)
	rc8_1s = stan_sample(m8_1s; data)
	if success(rc8_1s)
		post8_1s_df = read_samples(m8_1s, :dataframe)
		PRECIS(post8_1s_df[:, [:a, :b, :sigma]])
	end
end

# ╔═╡ 7e9e3674-6aea-11eb-3206-c1f775fe2d78
stan8_2 = "
data {
	int N;
	int K;
	vector[N] G;
	vector[N] R;
	int cid[N];
}

parameters {
	vector[K] a;
	real b;
	real<lower=0> sigma;
}

transformed parameters {
	vector[N] mu;
	for (i in 1:N)
		mu[i] = a[cid[i]] + b * (R[i] - 0.215);
}

model {
	a ~ normal(1, 0.1);
	b ~ normal(0, 0.3);
	sigma ~ exponential(1);
	G ~ normal(mu, sigma);
}
generated quantities {
	vector[N] log_lik;
	for (i in 1:N)
		log_lik[i] = normal_lpdf(G[i] | mu[i], sigma);
}
";

# ╔═╡ 12db71d0-6aeb-11eb-341c-e3243ba74862
begin
	m8_2s = SampleModel("m8.2s", stan8_2)
	rc8_2s = stan_sample(m8_2s; data)
	if success(rc8_2s)
		post8_2s_df = read_samples(m8_2s, :dataframe)
		PRECIS(post8_2s_df[:, [Symbol("a.1"), Symbol("a.2"), :b, :sigma]])
	end
end

# ╔═╡ 81c10f50-6af6-11eb-291b-a9dabe8a297c
stan8_3 = "
data {
	int N;
	int K;
	vector[N] G;
	vector[N] R;
	int cid[N];
}

parameters {
	vector[K] a;
	vector[K] b;
	real<lower=0> sigma;
}

transformed parameters {
	vector[N] mu;
	for (i in 1:N)
		mu[i] = a[cid[i]] + b[cid[i]] * (R[i] - 0.215);
}

model {
	a ~ normal(1, 0.1);
	b ~ normal(0, 0.3);
	sigma ~ exponential(1);
	G ~ normal(mu, sigma);
}
generated quantities {
	vector[N] log_lik;
	for (i in 1:N)
		log_lik[i] = normal_lpdf(G[i] | mu[i], sigma);
}
";

# ╔═╡ a83bde26-6af6-11eb-29ab-83fc8ed39292
begin
	m8_3s = SampleModel("m8.3s", stan8_3)
	rc8_3s = stan_sample(m8_3s; data)
	if success(rc8_3s)
		post8_3s_df = read_samples(m8_3s, :dataframe)
		PRECIS(post8_3s_df[:, [Symbol("a.1"), Symbol("a.2"),
					Symbol("b.1"), Symbol("b.2"), :sigma]])
	end
end

# ╔═╡ ba035d2a-6b2b-11eb-382d-313cc52aeb1a
compare([m8_1s, m8_2s, m8_3s], :waic)

# ╔═╡ a27a15de-7069-11eb-277c-2f2d7b4dfa8b
plot_models([m8_1s, m8_2s, m8_3s], :waic)

# ╔═╡ c1a6492a-6b2b-11eb-0f91-796cd055a303
compare([m8_1s, m8_2s, m8_3s], :psis)

# ╔═╡ 0721effc-6af7-11eb-1149-5f8050095c92
loo8_3s, loos8_3s, pk8_3s = psisloo(m8_3s)

# ╔═╡ 246d50ae-6af7-11eb-188e-5dd9a1382b32
pk_plot(pk8_3s)

# ╔═╡ 9331b922-6b14-11eb-0fa2-b7ac480b546c
begin
	nt8_3s = read_samples(m8_3s)
	xs = -0.1:0.01:1.1
	df_africa = df[df.cid .== 1, [:country, :log_gdp_s, :rugged_s, :cid]]
	mu_africa = nt8_3s.a[1,:]' .+ nt8_3s.b[1, :]' .* xs
	scatter(xs, mu_africa[:, 1:50], leg=false, markersize=1,
		color=:blue, alpha=0.1)
end

# ╔═╡ d60aaa22-6b2c-11eb-1e68-f1609c237a37
begin
	mlu_a = meanlowerupper(mu_africa)
	fig1 = plot(;ylims=(0.7, 1.4), xlab="ruggedness (standardized)",
		ylab="log GDP (standardaized)", title="African nations")
	scatter!(df_africa.rugged_s, df_africa.log_gdp_s, leg=false)
	plot!(xs, mlu_a.mean,
		ribbon=(mlu_a.mean .- mlu_a.lower, mlu_a.upper - mlu_a.mean))
	df_a = df_africa[df_africa.rugged_s .> 0.5, :]
	for r in eachrow(df_a)
		annotate!([(r.rugged_s+0.04, r.log_gdp_s+0.015,
			Plots.text(r.country, 6, :red, :right))])
	end
	plot!()
end

# ╔═╡ a2bcceb8-6b2f-11eb-03cc-b54f771a9174
begin
	df_other = df[df.cid .== 2, [:country, :log_gdp_s, :rugged_s, :cid]]
	mu_other = nt8_3s.a[2, :]' .+ nt8_3s.b[2, :]' .* xs
	scatter(xs, mu_other[:, 1:50], leg=false, markersize=1,
		color=:blue, alpha=0.1)
end

# ╔═╡ dea1e592-6b33-11eb-2497-df10350753a0
begin
	mlu_o = meanlowerupper(mu_other, (0.055, 0.945))
	fig2 = plot(;ylims=(0.7, 1.4), xlab="ruggedness (standardized)",
		ylab="log GDP (standardaized)", title="Non-African nations")

	scatter!(df_other.rugged_s, df_other.log_gdp_s, leg=false)
	plot!(xs, mlu_o.mean,
		ribbon=(mlu_o.mean .- mlu_o.lower, mlu_o.upper - mlu_o.mean))
	df_o = df_other[df_other.rugged_s .> 0.6, :]
	for r in eachrow(df_o)
		annotate!([(r.rugged_s+0.04, r.log_gdp_s+0.015,
			Plots.text(r.country, 6, :red, :right))])
	end
	plot!()
end

# ╔═╡ 9cd48a0e-6b34-11eb-05a5-939aed1ad460
plot(fig1, fig2, layout=(1, 2))

# ╔═╡ 45767e2e-6a63-11eb-3e12-354a7e32a374
md" ## End of clip-08-13-17s.jl"

# ╔═╡ Cell order:
# ╟─51fc19b8-59a8-11eb-2214-15aca59b807b
# ╠═63ba08cc-59a8-11eb-0a0f-27efac60d779
# ╠═6db218c6-59a8-11eb-2a8b-7107354cf590
# ╠═8aaa4bcc-59a8-11eb-2003-f1213b116565
# ╠═e890c59c-6a68-11eb-0d33-4d52a21d0ccf
# ╠═8e815d5c-6a66-11eb-21dc-99734c31f4e1
# ╠═7e9e3674-6aea-11eb-3206-c1f775fe2d78
# ╠═12db71d0-6aeb-11eb-341c-e3243ba74862
# ╠═81c10f50-6af6-11eb-291b-a9dabe8a297c
# ╠═a83bde26-6af6-11eb-29ab-83fc8ed39292
# ╠═ba035d2a-6b2b-11eb-382d-313cc52aeb1a
# ╠═a27a15de-7069-11eb-277c-2f2d7b4dfa8b
# ╠═c1a6492a-6b2b-11eb-0f91-796cd055a303
# ╠═0721effc-6af7-11eb-1149-5f8050095c92
# ╠═246d50ae-6af7-11eb-188e-5dd9a1382b32
# ╠═9331b922-6b14-11eb-0fa2-b7ac480b546c
# ╠═d60aaa22-6b2c-11eb-1e68-f1609c237a37
# ╠═a2bcceb8-6b2f-11eb-03cc-b54f771a9174
# ╠═dea1e592-6b33-11eb-2497-df10350753a0
# ╠═9cd48a0e-6b34-11eb-05a5-939aed1ad460
# ╟─45767e2e-6a63-11eb-3e12-354a7e32a374


text:
using BCTRNN
using DiffEqSensitivity
using OrdinaryDiffEq
import DiffEqFlux: FastChain, FastDense
import Flux: ClipValue, ADAM

# Not in Project.toml
using Plots
gr()

include("sine_wave_dataloader.jl")

function train_sine_fc(epochs, solver=Tsit5();
  sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)),
  T=Float32, in_size=10, model_size=5,
  tgrad=false, gen_jac=false, lr=0.02, kwargs...)

  train_dl = generate_2d_data(T)

  f_in = in_size
  f_out = 1

  # input_matrices = Dict(
  #     :w_sens => ones(T, f_in, model_size),
  #     :E_sens => ones(T, f_in, model_size)
  # )
  # syns_matrices = Dict(
  #     :w_syns => ones(T, model_size, model_size),
  #     :E_syns => ones(T, model_size, model_size)
  # )
  # im = BCTRNN.InputMatrices(input_matrices)
  # sm = BCTRNN.SynsMatrices(syns_matrices)

  im = BCTRNN.InputAllToAll()
  sm = BCTRNN.SynsFullyConnected()
  om = BCTRNN.OutputAll()
  wiring = BCTRNN.WiringConfig(f_in, model_size, im,sm,om)
  cell,defs,sys = BCTRNN.LTCOrigMTK(wiring, solver, sensealg; T, tgrad, gen_jac, kwargs...)
  model = FastChain(#(x,p)->ones(eltype(x), f_in, size(x)[2]),
                    FastDense(2, in_size),
                    cell,
                    FastDense(wiring.n_out, f_out))

  hs = []
  for (k,v) in wiring.matrices
    push!(hs, heatmap(v, title=k))
  end
  display(plot(hs..., layout=length(hs)))

  cb = BCTRNN.MyCallback(T; ecb=mycb, nepochs=epochs, nsamples=length(train_dl))
  #opt = GalacticOptim.Flux.Optimiser(ClipValue(0.5), ADAM(0.02))
  opt = BCTRNN.ClampBoundOptim(BCTRNN.get_bounds(model,T)..., ClipValue(T(1.0)), ADAM(T(lr)))
  BCTRNN.optimize(model, BCTRNN.loss_seq, cb, opt, train_dl, epochs, T), model, defs, sys
end

@time res, model, defs, sys = train_sine_fc(200, Tsit5(); lr=0.02, in_size=3, model_size=5,
tgrad=true, gen_jac=false,
save_everystep=false, save_start=false, save_end=true,
abstol=1e-4, reltol=1e-3#, dtmin=1e-3, dtmax=1.0
#dt=0.2f0, adaptive=false,
)


text:
module ThreadedStockCorr

# Threaded implementation of test case from Issue #445

using Base.Threads
if VERSION >= v"0.7.0-DEV.3406"
    using Random
end
if VERSION >= v"0.7.0-DEV.3449"
    using LinearAlgebra
else
    using Base.LinAlg
    const LinearAlgebra = Base.LinAlg
end

# Run paths in parallel (has to be in its own function due to #10718)
if VERSION >= v"0.7.0-DEV.3204"
    function runpath!(n, Wiener, CorrWiener, SA, SB, T, UpperTriangle, k11, k12, k21, k22, rngs)
        @threads for i = 1:n
            randn!(rngs[threadid()], Wiener)
            LinearAlgebra.mul!(CorrWiener, Wiener, UpperTriangle)
            @simd for j = 2:T
                @inbounds SA[j, i] = SA[j-1, i] * exp(k11 + k12*CorrWiener[j-1, 1])
                @inbounds SB[j, i] = SB[j-1, i] * exp(k21 + k22*CorrWiener[j-1, 2])
            end
        end
    end
else
    function runpath!(n, Wiener, CorrWiener, SA, SB, T, UpperTriangle, k11, k12, k21, k22, rngs)
        @threads for i = 1:n
            randn!(rngs[threadid()], Wiener)
            A_mul_B!(CorrWiener, Wiener, UpperTriangle)
            @simd for j = 2:T
                @inbounds SA[j, i] = SA[j-1, i] * exp(k11 + k12*CorrWiener[j-1, 1])
                @inbounds SB[j, i] = SB[j-1, i] * exp(k21 + k22*CorrWiener[j-1, 2])
            end
        end
    end
end

function perf_pstockcorr(n)
    ## Correlated asset information
    const CurrentPrice = [78. 102.]     # Initial Prices of the two stocks
    const Corr = [1. 0.4; 0.4 1.]       # Correlation Matrix
    const T = 500                       # Number of days to simulate = 2years = 500days
    const dt = 1/250                    # Time step (1year = 250days)
    const Div = [0.01 0.01]               # Dividend
    const Vol = [0.2 0.3]                 # Volatility

    ## Market Information
    const r = 0.03                      # Risk-free rate

    ## Define storages
    SimulPriceA = zeros(T,n)            # Simulated Price of Asset A
    SimulPriceA[1,:] = CurrentPrice[1]
    SimulPriceB = zeros(T,n)            # Simulated Price of Asset B
    SimulPriceB[1,:] = CurrentPrice[2]

    ## Generating the paths of stock prices by Geometric Brownian Motion
    const UpperTriangle = full(chol(Corr))    # UpperTriangle Matrix by Cholesky decomposition

    # Optimization: pre-allocate these for performance
    # NOTE: the new GC will hopefully fix this, but currently GC time
    # kills performance if we don't do in-place computations
    Wiener = Matrix{Float64}(undef, T-1, 2)
    CorrWiener = Matrix{Float64}(undef, T-1, 2)

    # Runtime requirement: need per-thread RNG since it stores state
    rngs = [MersenneTwister(777+x) for x in 1:nthreads()]

    # Optimization: pre-computable factors
    # NOTE: this should be automatically hoisted out of the loop
    k11 = (r-Div[1]-Vol[1]^2/2)*dt
    k12 = Vol[1]*sqrt(dt)
    k21 = (r-Div[2]-Vol[2]^2/2)*dt
    k22 = Vol[2]*sqrt(dt)

    runpath!(n, Wiener, CorrWiener, SimulPriceA, SimulPriceB, T, UpperTriangle, k11, k12, k21, k22, rngs)

    return (SimulPriceA, SimulPriceB)
end

end # module


text:
"""
    toeplitz{T}(c::AbstractArray{T},r::AbstractArray{T})

Returns a Toeplitz matrix where `c` is the first column and `r` is the first row.
"""
function toeplitz(c::AbstractArray{T},r::AbstractArray{T}) where T
    nc = length(c)
    nr = length(r)
    A = zeros(T, nc, nr)
    A[:,1] = c
    A[1,:] = r
    @views for i in 2:nr
        A[2:end,i] = A[1:end-1,i-1]
    end
    A
end

"""
    y,A = getARregressor(y::AbstractVector,na::Integer)

Returns a shortened output signal `y` and a regressor matrix `A` such that the least-squares
AR model estimate of order `na` is `y\\A`
"""
function getARregressor(y::AbstractVector,na)
    A = toeplitz(y[na+1:end],y[na+1:-1:1])
    y = copy(A[:,1])
    A = A[:,2:end]
    return y,A
end

"""
    getARXregressor(y::AbstractVector,u::AbstractVecOrMat, na, nb)

Returns a shortened output signal `y` and a regressor matrix `A` such that the least-squares
ARX model estimate of order `na,nb` is `y\\A`

Return a regressor matrix used to fit an ARX model on, e.g., the form
`A(z)y = B(z)f(u)`
with output `y` and input `u` where the order of autoregression is `na` and
the order of input moving average is `nb`

# Example
Here we test the model with the Function `f(u) = √(|u|)`
```julia
A     = [1,2*0.7*1,1] # A(z) coeffs
B     = [10,5] # B(z) coeffs
u     = randn(100) # Simulate 100 time steps with Gaussian input
y     = filt(B,A,sqrt.(abs.(u)))
yr,A  = getARXregressor(y,u,3,2) # We assume that we know the system order 3,2
bfe   = MultiUniformRBFE(A,[2,2,4,4,4], normalize=true)
bfa   = BasisFunctionApproximation(yr,A,bfe, 1e-3)
e_bfe = √(mean((yr - bfa(A)).^2)) # (0.005174261451622258)
plot([yr bfa(A)], lab=["Signal" "Prediction"])
```
See README (`?BasisFunctionExpansions`) for more details
"""
function getARXregressor(y::AbstractVector,u::AbstractVecOrMat, na, nb)
    @assert(length(nb) == size(u,2))
    m    = max(na+1,maximum(nb))
    n    = length(y) - m+1
    offs = m-na-1
    A    = toeplitz(y[offs+na+1:n+na+offs],y[offs+na+1:-1:1])
    y    = copy(A[:,1])
    A    = A[:,2:end]
    for i = 1:length(nb)
        offs = m-nb[i]
        A = [A toeplitz(u[nb[i]+offs:n+nb[i]+offs-1,i],u[nb[i]+offs:-1:1+offs,i])]
    end
    return y,A
end

"size(y) = (T-1, n)"
function matricesn(x,u)
    A = [x[1:end-1,:] u[1:end-1,:]]
    y = x[2:end,:]
    y,A
end

"Convenience tyoe for estimation of LPV state-space models"
struct LPVSS
    bfe::BasisFunctionExpansion
    params
    cov
    σ
end

Base.show(io::IO, model::LPVSS) = print(io,"LPVSS model with $(typeof(model.bfe)) and $(length(model.params)*length(model.params[1])) parameters")

"""
    LPVSS(x, u, nc; normalize=true, λ = 1e-3)

Linear Parameter-Varying State-space model. Estimate a state-space model with
varying coefficient matrices `x(t+1) = A(v)x(t) + B(v)u(t)`. Internally a `MultiRBFE` spanning
the space of `X × U` is used. `x` and `u` should have time in first dimension. Centers are
found automatically using k-means, see `MultiRBFE`.

# Examples
```jldoctest
using Plots, BasisFunctionExpansions
x,xm,u,n,m = BasisFunctionExpansions.testdata(1000)
nc         = 10 # Number of centers
model      = LPVSS(x, u, nc; normalize=true, λ = 1e-3) # Estimate a model
xh         = model(x,u) # Form prediction

eRMS       = √(mean((xh[1:end-1,:]-x[2:end,:]).^2))

plot(xh[1:end-1,:], lab="Prediction", c=:red, layout=2)
plot!(x[2:end,:], lab="True", c=:blue); gui()
eRMS <= 0.37

# output

true
```
"""
function LPVSS(x, u, nc; normalize=true, λ = 1e-3)
    y,A         = matricesn(x,u)
    bfe         = MultiRBFE(A, nc; normalize=normalize) # A is sched/regressor matrix
    params, cov, σ = fit_ss(x,u,A,bfe,λ)
    LPVSS(bfe, params, cov, σ)
end

"""
    LPVSS(x, u, v, nc; normalize=true, λ = 1e-3)

Linear Parameter-Varying State-space model. Estimate a state-space model with
varying coefficient matrices `x(t+1) = A(v)x(t) + B(v)u(t)`. Internally a `MultiRBFE` or
`UniformRBFE` spanning the space of `v` is used, depending on the dimensionality of
`v`. `x`, `u` and `v` should have time in first dimension. Centers are found automatically
using k-means, see `MultiRBFE`.

# Examples
```jldoctest
using Plots, BasisFunctionExpansions
T          = 1000
x,xm,u,n,m = BasisFunctionExpansions.testdata(T)
nc         = 4
v          = 1:T
model      = LPVSS(x, u, v, nc; normalize=true, λ = 1e-3)
xh         = model(x,u,v)

eRMS       = √(mean((xh[1:end-1,:]-x[2:end,:]).^2))

plot(xh[1:end-1,:], lab="Prediction", c=:red, layout=(2,1))
plot!(x[2:end,:], lab="True", c=:blue); gui()
eRMS <= 0.26

# output

true
```
"""
function LPVSS(x, u, v::AbstractVecOrMat, nc; normalize=true, λ = 1e-3)
    if isa(v,AbstractMatrix)
        bfe  = MultiRBFE(v, nc; normalize=normalize)
    else
        bfe  = UniformRBFE(v, nc; normalize=normalize)
    end
    params, cov, σ = fit_ss(x,u,v,bfe,λ)
    LPVSS(bfe, params, cov, σ)
end

function mega_regressor(bfe,v,A)
    nc = length(bfe.μ) ÷ (ndims(v) > 1 ? size(v,2) : 1)
    ϕ = bfe(v)
    if isa(ϕ, Vector)
        ϕ = ϕ'
    end
    ϕ = repeat(ϕ,1,size(A,2)) # Extend activations from nc to nc×(n+m)
    ϕ = ϕ.* repeat(A,1,nc)  # Extend regressor from n+m to nc×(n+m)
end

shorten_v(v::AbstractVector) = v[1:end-1]
shorten_v(v::AbstractMatrix) = v[1:end-1,:]

function fit_ss(x,u,v,bfe,λ)
    y,A  = matricesn(x,u) # A is sched matrix [x[1:end-1,:] u[1:end-1,:]]
    if size(v,1) > size(A,1)
        v = shorten_v(v)
    end
    ϕ = mega_regressor(bfe,v,A)
    p = size(ϕ,2)
    B = qr(λ == 0 ? ϕ : [ϕ; λ*I])
    params = mapslices(y,dims=1) do y
        if λ == 0
            x = B\y
        else
            x = B\[y;zeros(p)]
        end
    end # We now have n vectors of nc(n+m) parameters = nc(n+m)×n
    σ = [std(y[:,i] - ϕ*params[:,i]) for i =1:size(params,2)]
    ATA = Cholesky(B.R,:U,0) # TODO: Isn't R already a cholesky factor of R'R??
    cov = λ == 0 ? inv(ATA) : ATA\(ϕ'ϕ)/Matrix(ATA) # / not defined for factorizations https://github.com/JuliaLang/julia/issues/12436

    # icov = cholfact(Hermitian(q))
    # icov = cholfact(Symmetric(q[:R]'q[:R]))
    return params, cov, σ
end

"""
    predict(model::LPVSS, x::AbstractMatrix, u, v=[x u])

If no `v` provided, return a prediction of the output `x'` given the state `x` and input `u`
This function is called when a `model::LPVSS` object is called like `model(x,u)`

Provided `v`, return a prediction of the output `x'` given the state `x`, input `u` and
scheduling parameter `v`
"""
function predict(model::LPVSS, x::AbstractMatrix, u, v=[x u])
    A = [x u]
    ϕ = mega_regressor(model.bfe,v,A)
    y = similar(x)
    for i = 1:size(y,2)
        y[:,i] = ϕ*model.params[:,i]
    end
    y
end

function predict(model::LPVSS, x::AbstractVector, u, v)
    A = [x' u']
    ϕ = mega_regressor(model.bfe,v,A)
    y = map(1:length(x)) do i
        vecdot(ϕ,model.params[:,i])
    end
end

(model::LPVSS)(x,u)                 = predict(model, x, u)
(model::LPVSS)(x,u,v)               = predict(model, x, u, v)
(model::LPVSS)(x::AbstractMatrix,u) = predict(model, x, u, [x u])
(model::LPVSS)(x::AbstractVector,u) = predict(model, x, u, [x' u'])

function AB_covariance(model::LPVSS)
    ATA = model.cov
    σ²  = model.σ^2
    # ATA has size (n+m)×(n+m)
    # Full covariance matrix should have size n(n+m)×n(n+m), hence ATA should be composed in
    # a block form and multiplied with the corresponding σ²

end

"""
    output_variance(model::LPVSS, x::AbstractVector, u::AbstractVector, v=[x u])

Return a vector of prediction variances. Note, no covariance between dimensions in output is
provided
"""
function output_variance(model::LPVSS, x::AbstractVector, u::AbstractVector, v=[x u])
    A = [x' u']
    ϕ = mega_regressor(model.bfe,v,A)[:]
    mean_var = [σ^2*ϕ'model.cov*ϕ for σ in model.σ]
    pred_var = [model.σ[i] + sqrt(mean_var[i]) for i = 1:length(model.σ)]
end


text:
module GaussianTest

using Test
using ForneyLab
using LinearAlgebra: Diagonal

@testset "sample" begin
    # Univariate
    @test isa(sample(ProbabilityDistribution(Univariate, GaussianMeanVariance, m=1.0, v=2.0)), Float64)

    # Multivariate
    @test isa(sample(ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[1.2, 2.7], v=[2.0 -0.5; -0.5 1.5])), Vector{Float64})
    @test isa(sample(ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[1.2, 2.7], v=Diagonal([2.0, 1.5]))), Vector{Float64})
end

@testset "prod!" begin
    # Univariate
    @test ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=0.0, w=1.0) * ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=0.0, w=1.0) == ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=0.0, w=2.0)
    @test ProbabilityDistribution(Univariate, GaussianMeanVariance, m=0.0, v=1.0) * ProbabilityDistribution(Univariate, GaussianMeanVariance, m=0.0, v=1.0) == ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=0.0, w=2.0)
    @test ProbabilityDistribution(Univariate, GaussianMeanVariance, m=0.0, v=1.0) * ProbabilityDistribution(Univariate, PointMass, m=1.0) == ProbabilityDistribution(Univariate, PointMass, m=1.0)
    @test ProbabilityDistribution(Univariate, PointMass, m=1.0) * ProbabilityDistribution(Univariate, GaussianMeanVariance, m=0.0, v=1.0) == ProbabilityDistribution(Univariate, PointMass, m=1.0)

    # Multivariate
    @test ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=zeros(2), w=diageye(2)) * ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=zeros(2), w=diageye(2)) == ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=zeros(2), w=2.0*diageye(2))
    @test ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=zeros(2), v=diageye(2)) * ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=zeros(2), v=diageye(2)) == ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=zeros(2), w=2.0*diageye(2))
    @test ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=zeros(2), v=diageye(2)) * ProbabilityDistribution(Multivariate, PointMass, m=ones(2)) == ProbabilityDistribution(Multivariate, PointMass, m=ones(2))
    @test ProbabilityDistribution(Multivariate, PointMass, m=ones(2)) * ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=zeros(2), v=diageye(2)) == ProbabilityDistribution(Multivariate, PointMass, m=ones(2))
end

end #module


text:
# using Glob
import FileIO
using JLD2

using DataFrame

function ic_to_df(ic::Reproduce.ItemCollection; filter=["_GIT_INFO", "_SAVE", "save_dir"])
    dlist = [merge(item.parsed_args, Dict("folder_str"=>item.folder_str)) for item in ic]
    doa = Dict()
    for k in keys(dlist[1])
	if filter isa Nothing || !(k in filter)
	    doa[k] = [datum[k] for datum in dlist]
	end
    end
    DataFrame(doa)
end

function files_to_ic_to_df(files::Vector{String}; kwargs...)
    ic_arr = [ItemCollection(at(file)) for file in files]
    ic_to_df(ItemCollection(vcat(getfield.(ic_arr, :items)...)); kwargs...)
end

function load_settings_file(settings_file)
    dict = FileIO.load(settings_file)
    if "settings_dict" ∈ keys(dict)
        dict = dict["settings_dict"]
    end
    parsed_args = dict["parsed_args"]
    used_keys = dict["used_keys"]
    (parsed_args=parsed_args[used_keys], folder_str=basename(dirname(settings_file)))
end

function settings_dataframe(folders::Vector{<:AbstractString}; kwargs...)
    dfs = [build_dataframe(folder; kwargs...) for folder in folders]
    for i in 2:length(dfs)
	append!(dfs[1], dfs[i])
    end
    dfs[1]
end

function settings_dataframe(folder::AbstractString;
			 filter=nothing,
			 force=false,
			 settings_file="settings.jld2")

    dir = splitpath(folder)[end] == "data" ? folder : joinpath(folder, "data")
    dir_list = readdir(dir)

    cache_loc = joinpath(dir, "data_frame.jld2")
    id = hash(string(dir_list))
    if isfile(cache_loc)
        data = FileIO.load(cache_loc)
        if id == data["id"] && !force
            return data["data"]
        end
    end

    # df = DataFrame()
    doa = Dict()
    item = load_settings_file(joinpath(dir, dir_list[1], settings_file))
    dlist = merge(item.parsed_args, Dict("folder_str"=>item.folder_str))

    for (k,v) in dlist
        if isnothing(filter) || k ∉ filter
            doa[k] = [v]
        end
    end

    for p in dir_list
	if basename(p) ∈ ["item_col.jld2", "data_frame.jld2"]
	    continue
	end
        item = load_settings_file(joinpath(dir, p, settings_file))
        dlist = merge(item.parsed_args, Dict("folder_str"=>item.folder_str))
        for k in keys(doa)
           push!(doa[k], dlist[k])
        end
    end

    df = DataFrame(doa)
    FileIO.save(cache_loc, "data", df, "id", id)

    return df
end

export ItemCollection, search, details

"""
    Item
An Item in the experiment. Contains the parsed arguments.
"""
struct Item
    folder_str::String
    parsed_args::Dict
    hash_keys
end

function Item(settings_file::AbstractString)
    dict = FileIO.load(settings_file)
    if "settings_dict" ∈ keys(dict)
        dict = dict["settings_dict"]
    end
    parsed_args = dict["parsed_args"]
    used_keys = dict["used_keys"]
    return Item(dirname(settings_file), parsed_args, used_keys)
end

"""
    ItemCollection
A collection of items. Mostly helpful, but not really used yet.
"""
struct ItemCollection
    items::Array{Item,1}
    dir_hash::UInt64
end

ItemCollection(items::Array{Item, 1}) = ItemCollection(items, 0x0)

@deprecate ItemCollection(::AbstractString) settings_dataframe(::AbstractString)

function ItemCollection(dir::AbstractString; settings_file="settings.jld2", data_folder="data")

    dir = splitpath(dir)[end] == data_folder ? dir : joinpath(dir, data_folder)
    dir_list = readdir(dir)

    d = joinpath(dir, "item_col.jld2")
    id = hash(string(dir_list))
    if isfile(d)
        ic = FileIO.load(d)["ic"]
        if ic.dir_hash == id
            return ic
        end
    end

    items = Array{Item,1}()
    for p in dir_list
        if isfile(joinpath(dir, p, settings_file))
            push!(items, Item(joinpath(dir, p, settings_file)))
        end
    end

    ic = ItemCollection(items, id)
    FileIO.save(d, "ic", ic)

    return ic
end

Base.show(io::IO, ic::ItemCollection) =
    print(io, "ItemCollection(Size: ", length(ic), ", Dir Hash: ", ic.dir_hash, ")")

# Iterator
Base.eltype(::Type{ItemCollection}) = Item
Base.length(ic::ItemCollection) = length(ic.items)
Base.getindex(ic::ItemCollection, idx) = ic.items[idx]
Base.firstindex(ic::ItemCollection, idx) = firstindex(ic.items)
Base.lastindex(ic::ItemCollection, idx) = lastindex(ic.items)

Base.iterate(ic::ItemCollection, state=1) = state > length(ic) ? nothing : (ic[state], state + 1)

"""
    search

Search for specific entries, or a number of entries.
"""
function search(itemCollection::ItemCollection, search_dict)

    dict_keys = keys(search_dict)
    found_items = Array{Item, 1}()
    for (item_idx, item) in enumerate(itemCollection.items)
        if search_dict == filter(k->((k[1] in dict_keys)), item.parsed_args)
            push!(found_items, item)
        end
    end
    return ItemCollection(found_items)

end

function search(f::Function, ic::ItemCollection)
    found_items = Vector{Reproduce.Item}()
    for (item_idx, item) in enumerate(ic.items)
        if f(item)
            push!(found_items, item)
        end
    end
    return ItemCollection(found_items)
end

"""
    details

get details of the pointed directory
"""
function details(itemCollection::ItemCollection)
    # itemCollection = ItemCollection(dir; settings_file=settings_file)

    for item in itemCollection.items
        args = ["$(k):$(item.parsed_args[k])" for k in item.hash_keys]
        println(item.parsed_args["_SAVE"], " ", join(args, " "))
    end
end

details(dir::AbstractString; settings_file="settings.jld2") =
    details(ItemCollection(dir; settings_file=settings_file))

import Base.-
"""
    -(l::Dict{K,T}, r::Dict{K,T}) where {K, T}

Get the difference between two dictionaries. Helper function for diff
"""
function -(l::Dict{K, V}, r::Dict{K, V}) where {K, V}
    ret = Dict{K, V}()
    for key in keys(l)
        if l[key] != r[key]
            ret[key] = (l[key], r[key])
        end
    end
    return ret
end

-(l::Item, r::Item) = l.parsed_args - r.parsed_args

"""
    diff

get difference of the list of items.
"""
function Base.diff(items::Array{Item, 1};
                   exclude_keys::Union{Array{String,1}, Array{Symbol,1}} = Array{String, 1}(),
                   exclude_parse_values::Bool=true)

    kt = keytype(items[1].parsed_args)
    exclude_keys = kt.(exclude_keys)
    if exclude_parse_values == true
        exclude_keys = [exclude_keys; kt.([Reproduce.HASH_KEY, Reproduce.SAVE_NAME_KEY, Reproduce.GIT_INFO_KEY])]
    end
    diff_parsed = Dict{kt, Array}()
    for item in items
        tmp_dict = items[1] - item
        for key in filter((k)->k ∉ exclude_keys, keys(tmp_dict))
            if key ∉ keys(diff_parsed)
		diff_parsed[key] = Array{Any, 1}()
            end

            if tmp_dict[key][1] ∉ diff_parsed[key]
                push!(diff_parsed[key], tmp_dict[key][1])
				new_type = typeof(tmp_dict[key])
            end
            if tmp_dict[key][2] ∉ diff_parsed[key]
                push!(diff_parsed[key], tmp_dict[key][2])
            end
        end
    end
    for key in keys(diff_parsed)
	diff_parsed[key] = collect(promote(diff_parsed[key]...))
        sort!(diff_parsed[key])
    end
    return diff_parsed
end

Base.diff(itemCollection::ItemCollection; kwargs...) =
    Base.diff(itemCollection.items; kwargs...)


text:
using Test
import MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities

include("../dummy.jl")

remove_moi(x::String) = replace(x, "MathOptInterface." => "")
function compare_without_moi(x::String, y::String)
    @test remove_moi(x) == remove_moi(y)
end

@testset "IndexMap" begin
    map = MOIU.IndexMap()
    @test length(map) == 0
    compare_without_moi(sprint(show, map), "Utilities.IndexMap()")
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    cx = MOI.ConstraintIndex{MOI.SingleVariable, MOI.Integer}(1)
    cy = MOI.ConstraintIndex{MOI.SingleVariable, MOI.Integer}(2)
    map = MOIU.IndexMap(Dict(x => y), Dict(cx => cy))
    @test length(map) == 2
    # `x=>y` in Julia <= 1.1 and `x => y` in Julia >= 1.2
    if VERSION < v"1.2"
        x_y = string(x) *  "=>" * string(y)
    else
        x_y = string(x => y)
    end
    compare_without_moi(sprint(show, map), "Utilities.IndexMap($x_y,Pair{ConstraintIndex,ConstraintIndex}($cx, $cy))")
end

@testset "AUTOMATIC" begin
    src = DummyModel()
    dest = DummyModel()
    @test_throws ErrorException MOIU.automatic_copy_to(dest, src)
    try
        @test_throws ErrorException MOIU.automatic_copy_to(dest, src)
    catch err
        @test sprint(showerror, err) == "Model DummyModel does not" *
        " support copy with names."
    end
end
@testset "Default" begin
    @test !MOIU.supports_default_copy_to(DummyModel(), false)
    @test !MOIU.supports_default_copy_to(DummyModel(), true)
    model = MOIU.Model{Float64}()
    MOIT.failcopytestc(model)
    MOIT.failcopytestia(model)
    MOIT.failcopytestva(model)
    MOIT.failcopytestca(model)
    MOIT.copytest(model, MOIU.Model{Float64}())
end
@testset "Allocate-Load" begin
    @test !MOIU.supports_allocate_load(DummyModel(), false)
    @test !MOIU.supports_allocate_load(DummyModel(), true)
    mock = MOIU.MockOptimizer(MOIU.Model{Float64}(), needs_allocate_load=true)
    MOIT.failcopytestc(mock)
    MOIT.failcopytestia(mock)
    MOIT.failcopytestva(mock)
    MOIT.failcopytestca(mock)
    MOIT.copytest(mock, MOIU.Model{Float64}())
end

struct DummyEvaluator <: MOI.AbstractNLPEvaluator end

@testset "Create variables in same ordering when NLPBlock is used (#849)" begin
    model = MOIU.UniversalFallback(MOIU.Model{Float64}())
    a = MOI.add_variable(model)
    b, c = MOI.add_variables(model, 2)
    x, cx = MOI.add_constrained_variable(model, MOI.GreaterThan(0.0))
    y, cy = MOI.add_constrained_variables(model, MOI.Nonnegatives(1))
    nlp_data = MOI.NLPBlockData(
        [MOI.NLPBoundsPair(0.0, 1.0) for i in 1:5],
        DummyEvaluator(), false)
    MOI.set(model, MOI.NLPBlock(), nlp_data)
    copy = MOIU.UniversalFallback(MOIU.Model{Float64}())
    index_map = MOIU.default_copy_to(copy, model, true)
    for vi in [a, b, c, x, y[1]]
        @test index_map[vi] == vi
    end
    for ci in [cx, cy]
        @test index_map[ci] == ci
    end
end

# Model that distinguish variables created constrained
struct ConstrainedVariablesModel <: MOI.ModelLike
    allocate_load::Bool
    added_constrained::Vector{Bool}
end
function ConstrainedVariablesModel(allocate_load::Bool)
    return ConstrainedVariablesModel(allocate_load, Bool[])
end
function MOI.empty!(model::ConstrainedVariablesModel)
    empty!(model.added_constrained)
end

function MOIU.supports_default_copy_to(model::ConstrainedVariablesModel, ::Bool)
    return !model.allocate_load
end
function MOIU.supports_allocate_load(model::ConstrainedVariablesModel, ::Bool)
    return model.allocate_load
end
function MOI.copy_to(dest::ConstrainedVariablesModel, src::MOI.ModelLike; kws...)
    MOIU.automatic_copy_to(dest, src; kws...)
end

function MOI.add_variables(model::ConstrainedVariablesModel, n)
    m = length(model.added_constrained)
    for i in 1:n
        push!(model.added_constrained, false)
    end
    return MOI.VariableIndex.(m .+ (1:n))
end
function MOIU.allocate_variables(model::ConstrainedVariablesModel, n)
    return MOI.add_variables(model, n)
end
function MOIU.load_variables(model::ConstrainedVariablesModel, n)
end
function MOI.add_constrained_variables(model::ConstrainedVariablesModel,
                                       set::MOI.AbstractVectorSet)
    m = length(model.added_constrained)
    for i in 1:MOI.dimension(set)
        push!(model.added_constrained, true)
    end
    ci = MOI.ConstraintIndex{MOI.VectorOfVariables, typeof(set)}(m + 1)
    return MOI.VariableIndex.(m .+ (1:MOI.dimension(set))), ci
end
function MOIU.allocate_constrained_variables(model::ConstrainedVariablesModel, set::MOI.AbstractVectorSet)
    return MOI.add_constrained_variables(model, set)
end
function MOIU.load_constrained_variables(model::ConstrainedVariablesModel, vis::Vector{MOI.VariableIndex}, ci::MOI.ConstraintIndex{MOI.VectorOfVariables}, set::MOI.AbstractVectorSet)
end
function MOI.add_constraint(model::ConstrainedVariablesModel,
                            func::MOI.VectorOfVariables,
                            set::MOI.AbstractVectorSet)
    return MOI.ConstraintIndex{typeof(func), typeof(set)}(func.variables[1].value)
end
function MOIU.allocate_constraint(model::ConstrainedVariablesModel, func::MOI.VectorOfVariables, set::MOI.AbstractVectorSet)
    return MOI.add_constraint(model, func, set)
end
function MOIU.load_constraint(model::ConstrainedVariablesModel, ::MOI.ConstraintIndex{MOI.VectorOfVariables}, ::MOI.VectorOfVariables, ::MOI.AbstractVectorSet)
end

@testset "Duplicates in VectorOfVariables: $allocate_load" for allocate_load in [false, true]
    src = MOIU.Model{Int}()
    x = MOI.add_variables(src, 3)
    cx = MOI.add_constraint(src, [x[1], x[3], x[1], x[2]], MOI.Nonnegatives(4))
    y, cy = MOI.add_constrained_variables(src, MOI.Nonpositives(3))
    dest = ConstrainedVariablesModel(allocate_load)
    idxmap = MOI.copy_to(dest, src)
    for vi in x
        @test !dest.added_constrained[idxmap[vi].value]
    end
    for vi in y
        @test dest.added_constrained[idxmap[vi].value]
    end
end