text:
# This package contains implementations of common shapes, which can be used to
# construct `Mesh` types from those shapes.

"""
Utility function for constructing an icosphere triangle mesh.

Adapted from https://observablehq.com/@mourner/fast-icosphere-mesh
"""
function make_icosphere_mesh(order::Int64)
    # set up a 20-triangle icosahedron
    f = (1 + 5^0.5) / 2
    T = 4^order

    vertices = zeros((10 * T + 2) * 3)
    vertices[1:36] = vec(
        [
            -1 f 0
            1 f 0
            -1 -f 0
            1 -f 0
            0 -1 f
            0 1 f
            0 -1 -f
            0 1 -f
            f 0 -1
            f 0 1
            -f 0 -1
            -f 0 1
        ]',
    )
    triangles = vec(
        [
            1 12 6
            1 6 2
            1 2 8
            1 8 11
            1 11 12
            12 11 3
            6 12 5
            2 6 10
            8 2 9
            11 8 7
            4 10 5
            4 5 3
            4 3 7
            4 7 9
            4 9 10
            10 9 2
            5 10 6
            3 5 12
            7 3 11
            9 7 8
        ]',
    )

    v = 13
    midCache = (order == 0) ? nothing : Dict() # midpoint vertices cache to avoid duplicating shared vertices

    function addMidPoint(a::Int64, b::Int64)
        key = floor((a + b - 2) * (a + b - 1) / 2) + min(a, b) # Cantor's pairing function
        if midCache !== nothing && haskey(midCache, key)
            i = midCache[key]
            delete!(midCache, key)
            return i
        end
        midCache[key] = v
        for k in [1, 2, 3]
            vertices[3*(v-1)+k] = (vertices[3*(a-1)+k] + vertices[3*(b-1)+k]) / 2
        end
        i = v
        v += 1
        return i
    end

    trianglesPrev = triangles
    for i = 1:order
        # subdivide each triangle into 4 triangles
        triangles = zeros(Int64, size(trianglesPrev)[1] * 4)
        for k = 0:3:(size(trianglesPrev)[1]-3)
            v1 = trianglesPrev[k+1]
            v2 = trianglesPrev[k+2]
            v3 = trianglesPrev[k+3]
            a = addMidPoint(v1, v2)
            b = addMidPoint(v2, v3)
            c = addMidPoint(v3, v1)
            t = k * 4
            triangles[t+1] = v1
            triangles[t+2] = a
            triangles[t+3] = c
            triangles[t+4] = v2
            triangles[t+5] = b
            triangles[t+6] = a
            triangles[t+7] = v3
            triangles[t+8] = c
            triangles[t+9] = b
            triangles[t+10] = a
            triangles[t+11] = b
            triangles[t+12] = c
            t += 12
        end
        trianglesPrev = triangles
    end
    # normalize vertices
    for i = 0:3:(size(vertices)[1]-3)
        m = 1 / (vertices[i+1]^2 + vertices[i+2]^2 + vertices[i+3]^2)^0.5
        vertices[i+1] *= m
        vertices[i+2] *= m
        vertices[i+3] *= m
    end
    reshape(vertices, 3, :), reshape(triangles, 3, :)
end

"""
Create an icosphere as a tetrahedral mesh, with potential 0 at the boundary.
"""
function make_icosphere(order::Int64, radius::Float64 = 1.0)
    verts, tris = make_icosphere_mesh(order)
    points::Matrix{Float64} = radius .* [verts [0; 0; 0]] # add the origin
    num_points = size(points, 2)
    tets::Matrix{Int64} = [tris; repeat([num_points], 1, size(tris, 2))]
    @assert size(tets, 1) == 4 "sanity check tets are in the right format"
    potentials = zeros(Float64, num_points)
    potentials[end] = radius
    Object(Mesh(points, tets, potentials))
end

"""
Creates a cuve with center at the position com and side length two.
"""
function make_cube(size::Float64 = 1.0)
    """
    returns a 2x2x2 cube with center at the origin
    """
    cube_verts =
        size .* [
            -0.5 -0.5 -0.5 -0.5 0.5 0.5 0.5 0.5 0.0
            -0.5 -0.5 0.5 0.5 -0.5 -0.5 0.5 0.5 0.0
            -0.5 0.5 -0.5 0.5 -0.5 0.5 -0.5 0.5 0.0
        ]
    cube_tets = [
        1 4 1 1 1 1 8 4 8 7 4 8
        2 2 2 5 3 5 4 3 7 6 8 6
        3 3 6 6 7 7 7 7 6 5 2 2
        9 9 9 9 9 9 9 9 9 9 9 9
    ]
    cube_pots = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, size / 2.0]
    Object(Mesh(cube_verts, cube_tets, cube_pots))
end

export make_icosphere, make_cube

max_stars_repo_name: ekzhang/hydroelastics

max_issues_repo_name: ekzhang/hydroelastics

max_forks_repo_name: ekzhang/hydroelastics


text:
using RecyclePermanent
using RNGPool
using JLD2

include("../models/gkModel.jl")

setRNGs(12345)

θ = GKθ(3.0, 1.0, 2.0, 0.5)
model = makeGKModel(θ)
ys = simulateModel(model, 1000)

@save "data/gkData.jld" ys θ

ext: jl


text:
mutable struct RandSource{F}
    f::F
    n::Int64
    RandSource(f::F, n::Int64) where {F <: Function} = new{F}(f, n)
    RandSource(; f::F = randn, n::Int64 = 1) where {F <: Function} = new{F}(f, n)
end

MultiRandSource = Tuple{Tuple, Bool}

mutable struct RandSpec{R}
    rng::R #AbstractRNG
    seed::Int64 # RNG seed (when -1, a seed will be chosen at random using the global RNG)
    init::MultiRandSource # Number of draws necessary for init function
    effects::MultiRandSource # Number of draws necessary for effects function
    derivatives::MultiRandSource # Number of draws necessary for derivatives function
    update::MultiRandSource # Number of draws necessary for update function
    RandSpec(; rng::R = MersenneTwister(0), seed=-1, init=0, effects=0, derivatives=0, update=0) where {R <: AbstractRNG} =
        new{R}(rng, seed, init, effects, derivatives, update)
end

# convert(::Type{RandSource}, ::Nothing) = RandSource(n=0)
convert(::Type{RandSource}, n::Int64) = RandSource(n=n)

# convert(::Type{MultiRandSource}, ::Nothing) = (Vector{RandSource}(), true)
convert(::Type{MultiRandSource}, r::RandSource) = ((r,), true)
convert(::Type{MultiRandSource}, n::Int64) = (n > 0 ? (RandSource(n=n),) : (), true)
convert(::Type{MultiRandSource}, t::Tuple{Function, Int64}) = ((RandSource(t...),), true)
convert(::Type{MultiRandSource}, vt::Vector{T}) where {T<:Tuple{Function, Int64}} = (collect(RandSource(t...) for t in vt), false)

convert(::Type{RandSpec}, ::Nothing) = RandSpec()

function seed(r::RandSpec, force_seed::Bool = false)
    if force_seed || r.seed < 0
        Random.seed!(r.rng, rand(1:1000000))
    else
        Random.seed!(r.rng, r.seed)
    end
end

function draw(rand_spec::RandSpec, field::Symbol)
    mrs      = getfield(rand_spec, field)
    minimize = mrs[2]
    draws    = [r.f(rand_spec.rng, r.n) for r in mrs[1]]
    return minimize ? vcat(draws...) : draws
end


text:
@testset "rationalquad" begin
    rng = MersenneTwister(123456)
    x = rand(rng) * 2
    v1 = rand(rng, 3)
    v2 = rand(rng, 3)

    @testset "RationalQuadraticKernel" begin
        α = 2.0
        k = RationalQuadraticKernel(; α=α)

        @testset "RQ ≈ EQ for large α" begin
            @test isapprox(
                RationalQuadraticKernel(; α=1e9)(v1, v2),
                SqExponentialKernel()(v1, v2);
                atol=1e-6,
                rtol=1e-6,
            )
        end

        @test metric(RationalQuadraticKernel()) == SqEuclidean()
        @test metric(RationalQuadraticKernel(; α=2.0)) == SqEuclidean()
        @test repr(k) == "Rational Quadratic Kernel (α = $(α))"

        # Standardised tests.
        TestUtils.test_interface(k, Float64)
        test_ADs(x -> RationalQuadraticKernel(; alpha=x[1]), [α])
        test_params(k, ([α],))
    end

    @testset "GammaRationalQuadraticKernel" begin
        k = GammaRationalQuadraticKernel()

        @test repr(k) == "Gamma Rational Quadratic Kernel (α = 2.0, γ = 2.0)"

        @testset "Default GammaRQ ≈ RQ with rescaled inputs" begin
            @test isapprox(
                GammaRationalQuadraticKernel()(v1 ./ sqrt(2), v2 ./ sqrt(2)),
                RationalQuadraticKernel()(v1, v2),
            )
            a = 1.0 + rand()
            @test isapprox(
                GammaRationalQuadraticKernel(; α=a)(v1 ./ sqrt(2), v2 ./ sqrt(2)),
                RationalQuadraticKernel(; α=a)(v1, v2),
            )
        end

        @testset "GammaRQ ≈ EQ for large α with rescaled inputs" begin
            v1 = randn(2)
            v2 = randn(2)
            @test isapprox(
                GammaRationalQuadraticKernel(; α=1e9)(v1 ./ sqrt(2), v2 ./ sqrt(2)),
                SqExponentialKernel()(v1, v2);
                atol=1e-6,
                rtol=1e-6,
            )
        end

        @testset "GammaRQ(γ=1) ≈ Exponential with rescaled inputs for large α" begin
            v1 = randn(4)
            v2 = randn(4)
            @test isapprox(
                GammaRationalQuadraticKernel(; α=1e9, γ=1.0)(v1, v2),
                ExponentialKernel()(v1, v2);
                atol=1e-6,
                rtol=1e-6,
            )
        end

        @testset "GammaRQ ≈ GammaExponential for same γ and large α" begin
            v1 = randn(3)
            v2 = randn(3)
            γ = rand() + 0.5
            @test isapprox(
                GammaRationalQuadraticKernel(; α=1e9, γ=γ)(v1, v2),
                GammaExponentialKernel(; γ=γ)(v1, v2);
                atol=1e-6,
                rtol=1e-6,
            )
        end

        @test metric(GammaRationalQuadraticKernel()) == Euclidean()
        @test metric(GammaRationalQuadraticKernel(; γ=2.0)) == Euclidean()
        @test metric(GammaRationalQuadraticKernel(; γ=2.0, α=3.0)) == Euclidean()

        # Standardised tests.
        TestUtils.test_interface(k, Float64)
        a = 1.0 + rand()
        test_ADs(x -> GammaRationalQuadraticKernel(; α=x[1], γ=x[2]), [a, 1 + 0.5 * rand()])
        test_params(GammaRationalQuadraticKernel(; α=a, γ=x), ([a], [x]))
    end
end


text:
import Nemo.lift

function fmpz_poly_read!(a::fmpz_poly, b::String)
  f = ccall((:fopen, :libc), Ptr{Nothing}, (Cstring, Cstring), b, "r")
  ccall((:fmpz_poly_fread, :libflint), Nothing, (Ptr{Nothing}, Ref{fmpz}), f, a)
  ccall((:fclose), Nothing, (Ptr{Nothing}, ), f)
  return a
end

function mahler_measure_bound(f::fmpz_poly)
  return root(sum([coeff(f, i)^2 for i=0:degree(f)])-1, 2)+1
end

function _val(m::fmpz, p::fmpz)  ## not used here, but fast -> Bernstein
  n, r = divrem(m, p)
  if !iszero(r)
    return 0, m
  end
  m = n
  v = 1
  p2 = p^2
  v2, m = _val(m, p2)
  v += 2*v2
  n, r = divrem(m, p)
  if iszero(r)
    v += 1
    m = n
  end
  return v, m
end

function _log(p::fmpz, m::fmpz)
  return Int(ceil(log(m)/log(p)))
end
mutable struct fmpz_poly_factor_t
  c::Int64 # actually an fmpz, the content
  p::Int64 # actually a fmpz_poly_struct *
  exp:: Int64 # slong *
  num:: Int64
  alloc::Int64
  function fmpz_poly_factor_t()
    n = new()
    ccall((:fmpz_poly_factor_init, :libflint), Nothing, (Ref{fmpz_poly_factor_t}, ), n)
    finalizer(n, fmpz_poly_factor_clear)
    return n
  end
  function fmpz_poly_factor_clear(n::fmpz_poly_factor_t)
    ccall((:fmpz_poly_factor_clear, :libflint), Nothing, (Ref{fmpz_poly_factor_t}, ), n)
  end
end

function flint_factor(f::fmpz_poly)
  n = fmpz_poly_factor_t()
  ccall((:fmpz_poly_factor, :libflint), Nothing, (Ref{fmpz_poly_factor_t}, Ref{fmpz_poly}), n, f)
  return n
end

_last = 1
function vanHoeji(f_orig::fmpz_poly, trunk::Bool = true)

  f = f_orig
  g = gcd(f, derivative(f))
  while !isone(g)
    f = divexact(f, g)
    g = gcd(f, derivative(f))
  end
  # now f is squarefree...

  #Bounds:
  #according to Klueners et all, the coefficients of the log lattive of the factors
  # (f/f_i f_i') are bounded by
  #            choose(n-1, i) n M(f)
  # in the worst case, one has to lift up to this bound - which is maximal for
  # i=n/2 roughly.

  p = next_prime(degree(f))
  while discriminant(PolynomialRing(ResidueRing(FlintZZ, p))[1](f)) == 0
    p = next_prime(p)
  end

  n = degree(f)
  mmb = mahler_measure_bound(f)*n
  cld_bound = function(i)
    a = mmb*binomial(n-1, i)
    b = fmpz()
    ccall((:fmpz_poly_CLD_bound, :libflint), Nothing, (Ref{fmpz}, Ref{fmpz_poly}, Int64), b, f, i)
    println("Compare: $(nbits(a)) vs $(nbits(b)) as a bound for $i")
    return min(a,b)
  end

  # following Bill's paper
  # we build the matrix gradually...

  k = min(_log(fmpz(p), cld_bound(0)), _log(fmpz(p), cld_bound(n-1)))+1
  println("starting with k=k_max=$k")
  d = Int(ceil(2*nbits(n)/nbits(p))) * 1  + Int(ceil(0.2 * k))
  k_max = k

  H = factor_mod_pk_init(f, p)

  while true
    k = k_max
    println("now at precision $k_max+$d")
    pk = fmpz(p)^(d+k_max)
    R = ResidueRing(FlintZZ, pk)
    Rx = PolynomialRing(R)[1]

    println("lifting")
    @time ff = factor_mod_pk(H, d+k_max)
    ky = collect(keys(ff))
    ky = [Rx(x) for x= ky]
    Rf = Rx(f)
    r = length(ky)
    m = MatrixSpace(FlintZZ, r, r)(1)
    # gradual feeding...
    # build the complete p-adic matrix first, then worry about the version
    # for LLL:
    p_adic = MatrixSpace(R, r, n)()
    j = 1
    @time for fi in ky
      g = div(Rf, fi)
#      h = prod(ky[find(x -> x != fi, ky)])
#      @assert h == g
      @assert degree(g) == n - degree(fi)
      g *= derivative(fi)  # g = f/f[i] * f[i]'
      @assert degree(g) == n-1
      for i=0:degree(g)
        p_adic[j, i+1] = coeff(g, i)
      end
      j += 1
    end

    curr_col = -1
    bd = -1
    while true
      # need to add column(s) to m - otherwise nothing will work

      c = curr_col + 1
      if c >= n
        break
      end
      curr_col += 1
      println("adding col, now at $curr_col, $(nrows(m))")
      bd = _log(fmpz(p), cld_bound(n-1-c))+1 # bound for THIS coeff only
      if bd > k_max
        print_with_color(:red, "prec change for next iteration was $k_max is $bd for $c\n")
      end
      k_max = max(k_max, bd)
      rm = nrows(m)
      _n = MatrixSpace(R, min(r, rm), r)(sub(m, 1:min(r, rm), 1:r)) * sub(p_adic, 1:r, n-1-c+1:n-1-c+1)
      _n = lift(MatrixSpace(FlintZZ, min(r, rm), 1), _n)
      # now we need to down-scale this...
      for i=1:min(r, rm)
        if trunk && bd <= k
          _n[i,1] = div(_n[i,1] % fmpz(p)^(bd+d), fmpz(p)^bd)
        else
          _n[i,1] = div(_n[i,1], fmpz(p)^bd)
        end
      end
      m = hcat(m, vcat(_n, MatrixSpace(FlintZZ, nrows(m)-nrows(_n), 1)()))
      if k+d > bd
        m = vcat(m, MatrixSpace(FlintZZ, 1, ncols(m))())
        if trunk && bd <= k
          m[nrows(m), ncols(m)] = fmpz(p)^(d)
        else
          m[nrows(m), ncols(m)] = fmpz(p)^(k+d-bd)
        end
      else
        print_with_color(:yellow, "\n\n zero diag added \n\n")
        d *= 2
        break;
      end

      if curr_col <= 1 # *div(nbits(r), 1)
        println("not enough cols")
        continue
      end

      println("scale")
      for i = 1:r
        for j=1:nrows(m)
          m[j,i] *= r
        end
      end
      println("lll")
      @time rk, l = lll_with_removal(m, r*fmpz(r)^2 + fmpz(r)^2*curr_col)
      println("un scale")
      for i = 1:r
        for j=1:nrows(l)
          l[j,i] = div(l[j,i], r)
        end
      end

      if rk <1
        print_with_color(:red, "\n\nshit - loosing everything\n\n")
        return m
        break
      end
      @assert rk > 0
      if rk < nrows(m)
        println("loosing rows! now at $rk from $(nrows(m))")
      end
      if rk==1
        println("Poly is irreducible!")
        fact = Dict{typeof(f), Int}()
        Rx = PolynomialRing(ResidueRing(FlintZZ, p))[1]
        fact[f] = valuation(Rx(f_orig), Rx(f))
        return fact
      end
      m = sub(l, 1:rk, 1:ncols(m))

      _cmp = function(i,j)
        local k = 1
        while k< rk && l[k, i] == l[k, j]
          k += 1
        end
        return cmp(l[k, i], l[k, j])
      end
      _lt(i,j) = _cmp(i,j)<0

      s = sort(1:r, lt = _lt)

      i = 1
      de = 0
      bad = 0
      fact = Dict{typeof(f), Int}()
      cls = Set{Array{Int, 1}}()
      grp = Array{Int, 1}()
      while i <= r
        j = i+1
        push!(grp, s[i])
        while j<=r && _cmp(s[i],s[j])==0
          push!(grp, s[j])
          j += 1
        end
        push!(cls, grp)
        grp = Array{Int, 1}()
        i = j
      end

      println("found $cls of length $(length(cls)) for $(nrows(m))")
      bad = 0
      if length(cls) == nrows(m)
        println("testing")
        for grp in cls
          fa = lift(parent(f_orig), prod(ky[grp]))

          local k
          for k=0:degree(fa)
            x = coeff(fa, k)
            if x > div(pk,2)
              x -= pk
            end
            setcoeff!(fa, k, x)
          end
          if iszero(pseudorem(f, fa))
            fact[fa] = 1
            de += degree(fa)
          else
            println("crap $fa")
            bad += 1
            break
          end
        end
        if de >= n
          Rx = PolynomialRing(ResidueRing(FlintZZ, p))[1]
          Rf = Rx(f_orig)
          for k in keys(fact)
            fact[k] = valuation(Rf, Rx(k))
          end
          return fact
        end
        if bad > 1
          break
        end
      end
      if bad > 1
        break
      end
    end
    d *= 2
  end
end


text:
"""
    Variable_DR!(X::Array{Interval{Float64}},mult_lo::Vector{Float64},
                 mult_hi::Vector{Float64},LBD::Float64,UBD::Float64)

Uses duality multiplier from box constraints along with global upper bound to
tighten the size of the node, `Xt`, in place. Inputs:
* `X::Array{Interval{Float64}}`: Node over which to perform domain reduction
* `mult_lo::Vector{Float64}`: Multipliers for variable (low bound)
* `mult_hi::Vector{Float64}`: Multipliers for variable (high bound)
* `LBD::Float64` - Lower bound of problem
* `UBD::Float64` - Global upper bound for B&B algorithm
"""
function Variable_DR!(X::Array{Interval{Float64}},
                      mult_lo::Vector{Float64},
                      mult_hi::Vector{Float64},
                      LBD::Float64,
                      UBD::Float64)

  lower_cut::Float64 = 0.0
  upper_cut::Float64 = 0.0
  if (LBD<=UBD)
    for i = 1:length(X)
      if (mult_lo[i] > 0.0)
        lower_cut = X[i].hi - (UBD-LBD)/mult_lo[i]
        if (lower_cut > X[i].lo)
          X[i] = Interval(lower_cut,X[i].hi)
        end
      elseif (mult_hi[i] > 0.0)
        upper_cut = X[i].lo + (UBD-LBD)/mult_lo[i]
        if (upper_cut < X[i].hi)
          X[i] = Interval(X[i].lo,upper_cut)
        end
      end
    end
  end
end

"""
    Variable_DR_Imp!(X::Array{Interval{Float64}},mult_lo::Vector{Float64},
                     mult_hi::Vector{Float64},LBD::Float64,UBD::Float64,nx::Int64)

Uses duality multiplier from box constraints along with global upper bound to
tighten the size of the node, `Xt`, in place. Only tightens the X[(nx+1):(nx+np)]
intervals. Inputs:
* `X::Array{Interval{Float64}}`: Node over which to perform domain reduction
* `mult_lo::Vector{Float64}`: Multipliers for variable (low bound)
* `mult_hi::Vector{Float64}`: Multipliers for variable (high bound)
* `LBD::Float64` - Lower bound of problem
* `UBD::Float64` - Global upper bound for B&B algorithm
* `nx::Int64` - The number of state variables.
"""
function Variable_DR_Imp!(X::Array{Interval{Float64}},
                      mult_lo::Vector{Float64},
                      mult_hi::Vector{Float64},
                      LBD::Float64,
                      UBD::Float64,
                      nx::Int64)

  lower_cut::Float64 = 0.0
  upper_cut::Float64 = 0.0
  if (LBD<=UBD)
    for i = (nx+1):length(X)
      if (mult_lo[i-nx] > 0.0)
        lower_cut = X[i].hi - (UBD-LBD)/mult_lo[i-nx]
        if (lower_cut > X[i].lo)
          X[i] = Interval(lower_cut,X[i].hi)
        end
      elseif (mult_hi[i-nx] > 0.0)
        upper_cut = X[i].lo + (UBD-LBD)/mult_lo[i-nx]
        if (upper_cut < X[i].hi)
          X[i] = Interval(X[i].lo,upper_cut)
        end
      end
    end
  end
end


text:
import Yao: gradient
export num_gradient, diff_blocks, loss, circuit, QCOptProblem, QCOptGo!

"""
    QCOptProblem

General quantum circuit optimization problem interface
"""
abstract type QCOptProblem end

"""
    circuit(qop::QCOptProblem) -> AbstractBlock

circuit to optimize
"""
function circuit end

"""
    loss(qop::QCOptProblem) -> Real

Return the loss.
"""
function loss end

#####################################################

"""
    gradient(qop::QCOptProblem) -> Vector

the gradients with respect to `diff_blocks`.
"""
function gradient end

"""
    diff_blocks(qop::QCOptProblem) -> iterable

collection of all differentiable units.
"""
diff_blocks(qop::QCOptProblem) = collect(qop |> circuit, AbstractDiff)

"""
    num_gradient(qop::QCOptProblem) -> Vector

obtain the gradient numerically
"""
num_gradient(qop::QCOptProblem) = numdiff.(()->loss(qop), qop |> diff_blocks)

################# Optimization ###################
"""
    QCOptGo!{QT}

quantum circuit optimization problem optimizer.
"""
abstract type QCOptGo!{QT} end

include("QuGAN.jl")
include("QCBM.jl")


text:
struct SimpleTsit5 <: DiffEqBase.AbstractODEAlgorithm end
export SimpleTsit5

mutable struct SimpleTsit5Integrator{IIP, S, T, P, F} <: DiffEqBase.AbstractODEIntegrator{SimpleTsit5, IIP, S, T}
    f::F                  # eom
    uprev::S              # previous state
    u::S                  # current state
    tmp::S                # dummy, same as state
    tprev::T              # previous time
    t::T                  # current time
    t0::T                 # initial time, only for reinit
    dt::T                 # step size
    tdir::T
    p::P                  # parameter container
    u_modified::Bool
    ks::Vector{S}         # interpolants of the algorithm
    cs::SVector{6, T}     # ci factors cache: time coefficients
    as::SVector{21, T}    # aij factors cache: solution coefficients
    rs::SVector{22, T}    # rij factors cache: interpolation coefficients
end
const ST5I = SimpleTsit5Integrator

DiffEqBase.isinplace(::ST5I{IIP}) where {IIP} = IIP

#######################################################################################
# Initialization
#######################################################################################
function DiffEqBase.__init(prob::ODEProblem,alg::SimpleTsit5;
                         dt = error("dt is required for this algorithm"))
  simpletsit5_init(prob.f,DiffEqBase.isinplace(prob),prob.u0,
                   prob.tspan[1], dt, prob.p)
end

function DiffEqBase.__solve(prob::ODEProblem,alg::SimpleTsit5;
                          dt = error("dt is required for this algorithm"))
  u0 = prob.u0
  tspan = prob.tspan
  ts = Array(tspan[1]:dt:tspan[2])
  n = length(ts)
  us = Vector{typeof(u0)}(undef,n)
  @inbounds us[1] = _copy(u0)
  integ = simpletsit5_init(prob.f,DiffEqBase.isinplace(prob),prob.u0,
                           prob.tspan[1], dt, prob.p)
  # FSAL
  for i in 1:n-1
    step!(integ)
    us[i+1] = _copy(integ.u)
  end
  sol = DiffEqBase.build_solution(prob,alg,ts,us,
                                  calculate_error = false)
  DiffEqBase.has_analytic(prob.f) && DiffEqBase.calculate_solution_errors!(sol;timeseries_errors=true,dense_errors=false)
  sol
end

@inline function simpletsit5_init(f::F,
                         IIP::Bool, u0::S, t0::T, dt::T, p::P
                         ) where {F, P, T, S<:AbstractArray{T}}

    cs, as, rs = _build_tsit5_caches(T)
    ks = [zero(u0) for i in 1:7]

    !IIP && @assert S <: SArray

    integ = ST5I{IIP, S, T, P, F}(
        f, _copy(u0), _copy(u0), _copy(u0), t0, t0, t0, dt, sign(dt), p, true, ks, cs, as, rs
    )
end

@inline function _build_tsit5_caches(::Type{T}) where {T}

    cs = SVector{6, T}(0.161, 0.327, 0.9, 0.9800255409045097, 1.0, 1.0)

    as = SVector{21, T}(
        #=a21=# convert(T,0.161),
        #=a31=# convert(T,-0.008480655492356989),
        #=a32=# convert(T,0.335480655492357),
        #=a41=# convert(T,2.8971530571054935),
        #=a42=# convert(T,-6.359448489975075),
        #=a43=# convert(T,4.3622954328695815),
        #=a51=# convert(T,5.325864828439257),
        #=a52=# convert(T,-11.748883564062828),
        #=a53=# convert(T,7.4955393428898365),
        #=a54=# convert(T,-0.09249506636175525),
        #=a61=# convert(T,5.86145544294642),
        #=a62=# convert(T,-12.92096931784711),
        #=a63=# convert(T,8.159367898576159),
        #=a64=# convert(T,-0.071584973281401),
        #=a65=# convert(T,-0.028269050394068383),
        #=a71=# convert(T,0.09646076681806523),
        #=a72=# convert(T,0.01),
        #=a73=# convert(T,0.4798896504144996),
        #=a74=# convert(T,1.379008574103742),
        #=a75=# convert(T,-3.290069515436081),
        #=a76=# convert(T,2.324710524099774)
    )

    rs = SVector{22, T}(
        #=r11=# convert(T,1.0),
        #=r12=# convert(T,-2.763706197274826),
        #=r13=# convert(T,2.9132554618219126),
        #=r14=# convert(T,-1.0530884977290216),
        #=r22=# convert(T,0.13169999999999998),
        #=r23=# convert(T,-0.2234),
        #=r24=# convert(T,0.1017),
        #=r32=# convert(T,3.9302962368947516),
        #=r33=# convert(T,-5.941033872131505),
        #=r34=# convert(T,2.490627285651253),
        #=r42=# convert(T,-12.411077166933676),
        #=r43=# convert(T,30.33818863028232),
        #=r44=# convert(T,-16.548102889244902),
        #=r52=# convert(T,37.50931341651104),
        #=r53=# convert(T,-88.1789048947664),
        #=r54=# convert(T,47.37952196281928),
        #=r62=# convert(T,-27.896526289197286),
        #=r63=# convert(T,65.09189467479366),
        #=r64=# convert(T,-34.87065786149661),
        #=r72=# convert(T,1.5),
        #=r73=# convert(T,-4),
        #=r74=# convert(T,2.5),
    )

    return cs, as, rs
end

#######################################################################################
# Stepping
#######################################################################################
# IIP version for vectors and matrices
@inline @muladd function DiffEqBase.step!(integ::ST5I{true, S, T}) where {T, S}

    L = length(integ.u)

    c1, c2, c3, c4, c5, c6 = integ.cs;
    dt = integ.dt; t = integ.t; p = integ.p
    a21, a31, a32, a41, a42, a43, a51, a52, a53, a54,
    a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = integ.as

    k1, k2, k3, k4, k5, k6 = integ.ks; k7 = k1
    tmp = integ.tmp; f! = integ.f

    integ.uprev .= integ.u; uprev = integ.uprev

    if integ.u_modified
      f!(k1, uprev, p, t)
      integ.u_modified=false
    else
      k1 .= k7
    end

    @inbounds begin
    for i in 1:L
        tmp[i] = uprev[i]+dt*a21*k1[i]
    end
    f!(k2, tmp, p, t+c1*dt)
    for i in 1:L
        tmp[i] = uprev[i]+dt*(a31*k1[i]+a32*k2[i])
    end
    f!(k3, tmp, p, t+c2*dt)
    for i in 1:L
        tmp[i] = uprev[i]+dt*(a41*k1[i]+a42*k2[i]+a43*k3[i])
    end
    f!(k4, tmp, p, t+c3*dt)
    for i in 1:L
        tmp[i] = uprev[i]+dt*(a51*k1[i]+a52*k2[i]+a53*k3[i]+a54*k4[i])
    end
    f!(k5, tmp, p, t+c4*dt)
    for i in 1:L
        tmp[i] = uprev[i]+dt*(a61*k1[i]+a62*k2[i]+a63*k3[i]+a64*k4[i]+a65*k5[i])
    end
    f!(k6, tmp, p, t+dt)
    for i in 1:L
        integ.u[i] = uprev[i]+dt*(a71*k1[i]+a72*k2[i]+a73*k3[i]+a74*k4[i]+a75*k5[i]+a76*k6[i])
    end
    end
    f!(k7, integ.u, p, t+dt)

    integ.tprev = t
    integ.t += dt

    return  nothing
end

# OOP version for vectors and matrices
@inline @muladd function DiffEqBase.step!(integ::ST5I{false, S, T}) where {T, S}

    c1, c2, c3, c4, c5, c6 = integ.cs;
    dt = integ.dt; t = integ.t; p = integ.p
    a21, a31, a32, a41, a42, a43, a51, a52, a53, a54,
    a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = integ.as

    tmp = integ.tmp; f = integ.f
    integ.uprev = integ.u; uprev = integ.u

    if integ.u_modified
      k1 = f(uprev, p, t)
      integ.u_modified=false
    else
      @inbounds k1 = integ.ks[7];
    end

    tmp = uprev+dt*a21*k1
    k2 = f(tmp, p, t+c1*dt)
    tmp = uprev+dt*(a31*k1+a32*k2)
    k3 = f(tmp, p, t+c2*dt)
    tmp = uprev+dt*(a41*k1+a42*k2+a43*k3)
    k4 = f(tmp, p, t+c3*dt)
    tmp = uprev+dt*(a51*k1+a52*k2+a53*k3+a54*k4)
    k5 = f(tmp, p, t+c4*dt)
    tmp = uprev+dt*(a61*k1+a62*k2+a63*k3+a64*k4+a65*k5)
    k6 = f(tmp, p, t+dt)

    integ.u = uprev+dt*((a71*k1+a72*k2+a73*k3+a74*k4)+a75*k5+a76*k6)
    k7 = f(integ.u, p, t+dt)

    @inbounds begin # Necessary for interpolation
        integ.ks[1] = k7; integ.ks[2] = k2; integ.ks[3] = k3
        integ.ks[4] = k4; integ.ks[5] = k5; integ.ks[6] = k6
        integ.ks[7] = k7
    end

    integ.tprev = t
    integ.t += dt

    return  nothing
end

#######################################################################################
# Interpolation
#######################################################################################
# Interpolation function, OOP
@muladd function (integ::ST5I)(t::T) where {T}
    tnext, tprev, dt = integ.t, integ.tprev, integ.dt
    #@assert tprev ≤ t ≤ tnext
    θ = (t - tprev)/dt
    b1θ, b2θ, b3θ, b4θ, b5θ, b6θ, b7θ = bθs(integ.rs, θ)

    ks = integ.ks
    if !isinplace(integ)
        u = @inbounds integ.uprev + dt*(b1θ*ks[1] + b2θ*ks[2] + b3θ*ks[3] + b4θ*ks[4] +
                      b5θ*ks[5] + b6θ*ks[6] + b7θ*ks[7])
        return u
    else
        u = similar(integ.u)
        @inbounds for i in 1:length(u)
            u[i] = integ.uprev[i] + dt*(b1θ*ks[1][i] + b2θ*ks[2][i] + b3θ*ks[3][i] +
                   b4θ*ks[4][i] + b5θ*ks[5][i] + b6θ*ks[6][i] + b7θ*ks[7][i])
        end
        return u
    end
end
# Interpolation coefficients
@inline function bθs(rs::SVector{22, T}, θ::T) where {T}
    # θ in (0, 1) !
    r11,r12,r13,r14,r22,r23,r24,r32,r33,r34,r42,r43,r44,r52,r53,
    r54,r62,r63,r64,r72,r73,r74 = rs

    b1θ::T = @evalpoly(θ, zero(T), r11    , r12, r13, r14)
    b2θ::T = @evalpoly(θ, zero(T), zero(T), r22, r23, r24)
    b3θ::T = @evalpoly(θ, zero(T), zero(T), r32, r33, r34)
    b4θ::T = @evalpoly(θ, zero(T), zero(T), r42, r43, r44)
    b5θ::T = @evalpoly(θ, zero(T), zero(T), r52, r53, r54)
    b6θ::T = @evalpoly(θ, zero(T), zero(T), r62, r63, r64)
    b7θ::T = @evalpoly(θ, zero(T), zero(T), r72, r73, r74)

    return b1θ, b2θ, b3θ, b4θ, b5θ, b6θ, b7θ
end


text:
using FFTW, DSP, Test, DelimitedFiles

# not exported, but used for some internal testing
using DSP.Windows: makewindow
@testset "makewindow" begin
    # make sure we're evaluating the given function at the correct points for
    # all combinations of arguments
    @test makewindow(identity, 6, 0, false) ≈ [-0.5, -0.3, -0.1, 0.1, 0.3, 0.5]
    win = makewindow(identity, 6, 0, true)
    # doesn't matter whether the middle value is evaluated at +/- 0.5
    @test win ≈ [0.0, 1/6, 1/3, 1/2, -1/3, -1/6] || win ≈ [0.0, 1/6, 1/3, -1/2, -1/3, -1/6]
    # we actually only end up with one real zero here, because one of the
    # "padding" values is used to split the -1/2 from +1/2. For windows that go
    # to zero padding will actually add `padding` zeros, because evaluting at
    # +/- 0.5 will be zero.
    @test makewindow(identity, 6, 2, false) ≈ [-0.5, -0.3, -0.1, 0.1, 0.3, 0.5, 0.0, 0.0]
    @test makewindow(identity, 6, 2, true) ≈ [0.0, 1/6, 1/3, 1/2, 0.0, -1/2, -1/3, -1/6]
    @test makewindow(identity, 5, 0, false) ≈ [-0.5, -0.25, 0.0, 0.25, 0.5]
    @test makewindow(identity, 5, 0, true) ≈ [0.0, 0.2, 0.4, -0.4, -0.2]
    @test makewindow(identity, 5, 2, false) ≈ [-0.5, -0.25, 0.0, 0.25, 0.5, 0.0, 0.0]
    @test makewindow(identity, 5, 2, true) ≈ [0.0, 0.2, 0.4, 0.0, 0.0, -0.4, -0.2]

    @test makewindow(x->42.0, 1, 0, false) ≈ [42.0]
    @test makewindow(x->42.0, 1, 0, true) ≈ [42.0]
    @test makewindow(x->42.0, 1, 2, false) ≈ [42.0, 0.0, 0.0]
    @test makewindow(x->42.0, 1, 2, true) ≈ [42.0, 0.0, 0.0]
end

@testset "dspss" begin
    # Test dpss against dpss computed with MATLAB
    d1 = dpss(128, 4)
    d2 = readdlm(joinpath(dirname(@__FILE__), "data", "dpss128,4.txt"), '\t')
    @test d1 ≈ d2

    # Test dpsseig against dpss from MATLAB
    lambda = [0.9999999997159923,0.9999999731146645,0.9999988168667646,
              0.9999680890685374,0.9994167543397652,0.9925560207018469,
              0.9368556668429153]
    @test dpsseig(d1, 4) ≈ lambda
end

@testset "common windows" begin
    # Checking Hanning, Hamming, Triangular, Bartlett, Bartlett-Hann, Tukey,
    # and Blackman windows against values computed with MATLAB.
    # Lanczos and cosine are checked against values generated with DSP.jl v0.4.0
    # to test for regressions, as there's no reference MATLAB implementation
    # Gaussian is compared against DSP.jl from commit da1b195, when the
    # implementation was corrected (see GH issue #204)
    @test rect(128) == ones(128)

    hanning_jl = hanning(128)
    hanning_ml = readdlm(joinpath(dirname(@__FILE__), "data", "hanning128.txt"), '\t')
    @test hanning_jl ≈ hanning_ml

    hamming_jl = hamming(128)
    hamming_ml = readdlm(joinpath(dirname(@__FILE__), "data", "hamming128.txt"), '\t')
    @test hamming_jl ≈ hamming_ml

    triang_jl = triang(128)
    triang_ml = readdlm(joinpath(dirname(@__FILE__), "data", "triang128.txt"), '\t')
    @test triang_jl ≈ triang_ml

    # for odd `n` the `triang` window should be the middle n-2 samples of the
    # bartlett
    @test triang(5) ≈ bartlett(7)[2:6]

    bartlett_jl = bartlett(128)
    bartlett_ml = readdlm(joinpath(dirname(@__FILE__), "data", "bartlett128.txt"), '\t')
    @test bartlett_jl ≈ bartlett_ml

    barthann_jl = bartlett_hann(128)
    barthann_ml = readdlm(joinpath(dirname(@__FILE__), "data", "bartlett_hann128.txt"), '\t')
    @test bartlett_jl ≈ bartlett_ml

    blackman_jl = blackman(128)
    blackman_ml = readdlm(joinpath(dirname(@__FILE__), "data", "blackman128.txt"), '\t')
    @test blackman_jl ≈ blackman_ml

    kaiser_jl = kaiser(128, 0.4/π)
    kaiser_ml = readdlm(joinpath(dirname(@__FILE__), "data", "kaiser128,0.4.txt"), '\t')
    @test kaiser_jl ≈ kaiser_ml

    gaussian_jl = gaussian(128, 0.2)
    gaussian_ref = readdlm(joinpath(dirname(@__FILE__), "data", "gaussian128,0.2.txt"), '\t')
    @test gaussian_jl ≈ gaussian_ref

    tukey_jl = tukey(128, 0.4)
    tukey_ml = readdlm(joinpath(dirname(@__FILE__), "data", "tukey128,0.4.txt"), '\t')
    @test tukey_jl ≈ tukey_ml

    @test tukey(128, 0) == rect(128)

    lanczos_jl = lanczos(128)
    lanczos_ref = readdlm(joinpath(dirname(@__FILE__), "data", "lanczos128.txt"), '\t')
    @test lanczos_jl ≈ lanczos_ref

    cosine_jl = cosine(128)
    cosine_ref = readdlm(joinpath(dirname(@__FILE__), "data", "cosine128.txt"), '\t')
    @test cosine_jl ≈ cosine_ref
end

zeroarg_wins = [rect, hanning, hamming, cosine, lanczos,
                bartlett, bartlett_hann, blackman, triang]
onearg_wins = [gaussian, kaiser, tukey]
@testset "zero-phase windows" begin
    for winf in zeroarg_wins
        if winf == triang
            # triang needs to be special-cased here because it has different
            # definitions for odd and even `n` (and also the underlying
            # continuous function changes with `n`), so it doesn't match the
            # zerophase assumption below
            @test triang(6, zerophase=true) ≈ [1.0, 0.75, 0.5, 0.25, 0.5, 0.75]
            @test triang(7, zerophase=true) ≈ [1.0, 0.75, 0.5, 0.25, 0.25, 0.5, 0.75]
        else
            @test winf(8, zerophase=true) ≈ ifftshift(winf(9)[1:8])
            # the `zerophase=false` version doesn't hit the center point, but the
            # zerophase one does, so this ends up introducing a 1/2-sample shift.
            @test winf(9, zerophase=true) ≈ ifftshift(winf(19)[2:2:end])
        end
    end

    # test the window functions that need extra args
    for winf in onearg_wins
        @test winf(8, 0.5, zerophase=true) == ifftshift(winf(9, 0.5)[1:8])
        @test winf(9, 0.5, zerophase=true) == ifftshift(winf(19, 0.5)[2:2:end])
    end

    @test dpss(8, 2, 1, zerophase=true)[:] == ifftshift(dpss(9, 2, 1)[1:8])
    # odd-length zerophase dpss windows not currently supported
    @test_throws ArgumentError dpss(9, 2, 1, zerophase=true)
end

@testset "window return types" begin
    # make sure return types are correct
    n = 12
    ft = typeof(1.0)
    for W in(:rect, :hanning, :hamming, :cosine, :lanczos, :triang, :bartlett, :bartlett_hann, :blackman)
        @eval @test Array{$ft,1} == typeof($W($n)) && length($W($n)) == $n
    end

    @test Array{ft,1} == typeof(tukey(n, 0.4)) && length(tukey(n, 0.4)) == n
    @test Array{ft,1} == typeof(gaussian(n, 0.4)) && length(gaussian(n, 0.4)) == n
    @test Array{ft,1} == typeof(kaiser(n, 0.4)) && length(kaiser(n, 0.4)) == n
    @test Array{ft,2} == typeof(dpss(n, 1.5)) && size(dpss(n, 1.5),1) == n  # size(,2) depends on the parameters
end

@testset "tensor product windows" begin
    # test all combinations of arguments. Each arg and kwarg can be not present,
    # a single value, or a 2-tuple
    for winf in [zeroarg_wins; onearg_wins]
        for arg in (nothing, 0.4, (0.4, 0.5))
            # skip invalid combinations
            winf in zeroarg_wins && arg !== nothing && continue
            winf in onearg_wins && arg === nothing && continue
            for padding in (nothing, 4, (4,5))
                for zerophase in (nothing, true, (true,false))
                    w1_expr = :($winf(15))
                    w2_expr = :($winf(20))
                    w3_expr = :($winf((15,20)))

                    if arg isa Real
                        push!(w1_expr.args, arg)
                        push!(w2_expr.args, arg)
                        push!(w3_expr.args, arg)
                    elseif arg isa Tuple
                        push!(w1_expr.args, arg[1])
                        push!(w2_expr.args, arg[2])
                        push!(w3_expr.args, arg)
                    end

                    if padding isa Integer
                        push!(w1_expr.args, Expr(:kw, :padding, padding))
                        push!(w2_expr.args, Expr(:kw, :padding, padding))
                        push!(w3_expr.args, Expr(:kw, :padding, padding))
                    elseif padding isa Tuple
                        push!(w1_expr.args, Expr(:kw, :padding, padding[1]))
                        push!(w2_expr.args, Expr(:kw, :padding, padding[2]))
                        push!(w3_expr.args, Expr(:kw, :padding, padding))
                    end

                    if zerophase isa Bool
                        push!(w1_expr.args, Expr(:kw, :zerophase, zerophase))
                        push!(w2_expr.args, Expr(:kw, :zerophase, zerophase))
                        push!(w3_expr.args, Expr(:kw, :zerophase, zerophase))
                    elseif zerophase isa Tuple
                        push!(w1_expr.args, Expr(:kw, :zerophase, zerophase[1]))
                        push!(w2_expr.args, Expr(:kw, :zerophase, zerophase[2]))
                        push!(w3_expr.args, Expr(:kw, :zerophase, zerophase))
                    end

                    w1 = eval(w1_expr)
                    w2 = eval(w2_expr)
                    w3 = eval(w3_expr)
                    @test w3 ≈ w1 * w2'
                end
            end
        end
    end
end


text:
"""
    bar(x[, y])

Compute the Bar index between `x` and `y`. If `y` is missing, compute
the Bar index between all pairs of columns of `x`.

# Examples
```julia-repl
julia> bar([1, 2], [1, 2])
1
```
"""
function bar(x, y) ...

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# MathOptInterface Variables

##################################################
## Getters
MOI.get(model::Optimizer, ::MOI.NumberOfVariables) = length(model.variable_info)

function MOI.get(model::Optimizer, ::MOI.ListOfVariableIndices)
    return [MOI.VariableIndex(i) for i in 1:length(model.variable_info)]
end

##################################################
## Bridges with KNITRO.
function MOI.add_variable(model::Optimizer)
    # If model has been optimized, KNITRO does not support adding
    # another variable.
    (model.number_solved >= 1) && throw(AddVariableError())
    push!(model.variable_info, VariableInfo())
    KN_add_var(model.inner)
    return MOI.VariableIndex(length(model.variable_info))
end
function MOI.add_variables(model::Optimizer, n::Int)
    return [MOI.add_variable(model) for i in 1:n]
end

function check_inbounds(model::Optimizer, vi::MOI.VariableIndex)
    num_variables = length(model.variable_info)
    if !(1 <= vi.value <= num_variables)
        error("Invalid variable index $vi. ($num_variables variables in the model.)")
    end
    return
end

##################################################
## Check inbounds for safety.
function check_inbounds(model::Optimizer, var::MOI.SingleVariable)
    return check_inbounds(model, var.variable)
end

function check_inbounds(model::Optimizer, aff::MOI.ScalarAffineFunction)
    for term in aff.terms
        check_inbounds(model, term.variable_index)
    end
    return
end

function check_inbounds(model::Optimizer, quad::MOI.ScalarQuadraticFunction)
    for term in quad.affine_terms
        check_inbounds(model, term.variable_index)
    end
    for term in quad.quadratic_terms
        check_inbounds(model, term.variable_index_1)
        check_inbounds(model, term.variable_index_2)
    end
    return
end

function has_upper_bound(model::Optimizer, vi::MOI.VariableIndex)
    return model.variable_info[vi.value].has_upper_bound
end

function has_lower_bound(model::Optimizer, vi::MOI.VariableIndex)
    return model.variable_info[vi.value].has_lower_bound
end

function is_fixed(model::Optimizer, vi::MOI.VariableIndex)
    return model.variable_info[vi.value].is_fixed
end

##################################################
## PrimalStart
function MOI.supports(::Optimizer, ::MOI.VariablePrimalStart,
                      ::Type{MOI.VariableIndex})
    return true
end
function MOI.set(model::Optimizer, ::MOI.VariablePrimalStart,
                 vi::MOI.VariableIndex, value::Union{Real, Nothing})
    check_inbounds(model, vi)
    if isa(value, Real)
        KN_set_var_primal_init_values(model.inner, vi.value - 1, Cdouble(value))
    else
        # By default, initial value is set to 0.
        KN_set_var_primal_init_values(model.inner, vi.value - 1, Cdouble(0.))
    end
    return
end

##################################################
## Naming
# WIP: this part aims at supporting JuMP directly, without caching optimizer.
MOI.supports(::Optimizer, ::MOI.VariableName, ::Type{MOI.VariableIndex}) = true
function MOI.set(model::Optimizer, ::MOI.VariableName, vi::MOI.VariableIndex, name::String)
    model.variable_info[vi.value].name = name
    return
end
function MOI.get(model::Optimizer, ::MOI.VariableName, vi::MOI.VariableIndex)
    return model.variable_info[vi.value].name
end


text:
# accurate accumlulation

function SumOfSquares(x::Vector{Float64})
    if length(x) >= 375299968947538
        throw(DomainError)
    end

    SumSquares(x)
end

function SumOfSquares(x::Vector{Float32})
    if length(x) >= 699047
        throw(DomainError)
    end

    SumSquares(x)
end

function SumSquares{T<:AbstractFloat}(x::Vector{T})
    ssqHiLo = zero(FF{T})
    for i in 1:length(x)
        xsqHiLo = eftMul(x[i], x[i])
        ssqHiLo = SumNonNeg(ssqHiLo, xsqHiLo)
    end
    ssqHiLo
end

function SumNonNeg{T<:AbstractFloat}(a::FF{T}, b::FF{T})
    hi,md = eftAdd(a.hi,b.hi)
    c = a.lo + b.lo
    d = lo + c
    hi, md = eftAddGTE(hi,d)
    FF(hi, md)
end


text:
using TrajectoryOptimization, RobotDynamics
using StaticArrays, LinearAlgebra
using TrajectoryOptimization: QuadraticObjective, cost_expansion!
using Altro
using Test
const TO = TrajectoryOptimization

n,m,N = 12,6,11
dt = 0.1
tf = dt*(N-1)
A = rand(n,n)
B = rand(n,m)
model = LinearModel(A,B,dt=dt)

# Objective
Q = Diagonal(@SVector rand(n))
R = Diagonal(@SVector rand(m))
Qf = Q*N
x0 = rand(n)
xf = 10*x0
obj = LQRObjective(Q,R,Qf,xf,N)
U0 = [randn(m) for k = 1:N-1]

# Linear Constraints
cons0 = ConstraintList(n,m,N)
add_constraint!(cons0, GoalConstraint(xf), N)

cons_linear = copy(cons0)
add_constraint!(cons_linear, BoundConstraint(n,m, x_min=-2, x_max=2, u_min=-rand(m), u_max=rand(m)), 1:N-1)
add_constraint!(cons_linear, LinearConstraint(n,m,rand(3,n+m), rand(3),Inequality()), 1:N-1)

cons_nonlinear = copy(cons_linear)
add_constraint!(cons_nonlinear, CircleConstraint(n, [1,1],[1,2],[0.1,0.1]), 1:N-1)

prob = Problem(model, obj, xf, tf, x0=x0, U0=U0, integration=PassThrough)
rollout!(prob)
Z = prob.Z

E = QuadraticObjective(n,m,N)
cost_expansion!(E, obj, Z)
@test TO.is_quadratic(E)

# AL Objective
alobj = Altro.ALObjective(obj, cons0, model)
cost_expansion!(E, alobj, Z)
@test TO.is_quadratic(E)

alobj = Altro.ALObjective(obj, cons_linear, model)
cost_expansion!(E, alobj, Z)
@test !TO.is_quadratic(E)

# with with nonlinear constraint, should no longer be quadratic
alobj = Altro.ALObjective(obj, cons_nonlinear, model)
cost_expansion!(E, alobj, Z)
@test !TO.is_quadratic(E)


text:
using DashBootstrapComponents, DashHtmlComponents

list_group = html_div([
    dbc_listgroup(
        [
            dbc_listgroupitem("Item 1"),
            dbc_listgroupitem("Item 2"),
            dbc_listgroupitem("Item 3"),
        ],
        horizontal = true,
        className = "mb-2",
    ),
    dbc_listgroup(
        [
            dbc_listgroupitem("Item 1"),
            dbc_listgroupitem("Item 2"),
            dbc_listgroupitem("Item 3"),
        ],
        horizontal = "lg",
    ),
]);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
# ---
# title: 961. N-Repeated Element in Size 2N Array
# id: problem961
# author: Indigo
# date: 2021-01-25
# difficulty: Easy
# categories: Hash Table
# link: <https://leetcode.com/problems/n-repeated-element-in-size-2n-array/description/>
# hidden: true
# ---
#
# In a array `A` of size `2N`, there are `N+1` unique elements, and exactly one
# of these elements is repeated `N` times.
#
# Return the element repeated `N` times.
#
#
#
# **Example 1:**
#
#
#
#     Input: [1,2,3,3]
#     Output: 3
#
#
# **Example 2:**
#
#
#
#     Input: [2,1,2,5,3,2]
#     Output: 2
#
#
# **Example 3:**
#
#
#
#     Input: [5,1,5,2,5,3,5,4]
#     Output: 5
#
#
#
#
# **Note:**
#
#   * `4 <= A.length <= 10000`
#   * `0 <= A[i] < 10000`
#   * `A.length` is even
#
#
## @lc code=start
using LeetCode

function repeated_n_times(A::Vector{Int})
    for i in 2:length(A)
        if A[i] == A[i - 1]
            return A[i]
        end
    end

    ## 3X3X OR 3XX3 OR X3X3
    return (A[1] == A[3] || A[1] == A[4]) ? A[1] : A[2]
end
## @lc code=end


text:
# ---
# title: 504. Base 7
# id: problem504
# author: zhwang
# date: 2022-03-10
# difficulty: Easy
# categories:
# link: <https://leetcode.com/problems/base-7/description/>
# hidden: true
# ---
#
# Given an integer, return its base 7 string representation.
#
# **Example 1:**
#
#
#
#     Input: 100
#     Output: "202"
#
#
# **Example 2:**
#
#
#
#     Input: -7
#     Output: "-10"
#
#
# **Note:** The input will be in range of [-1e7, 1e7].
#
#
## @lc code=start
using LeetCode

function convert_to_base7(num::Int)::String
    num == 0 && return "0"
    sign = num < 0 ? -1 : 1
    num *= sign
    res = Int[]
    while num > 0
        push!(res, num % 7)
        num ÷= 7
    end
    reverse!(res)
    return sign != -1 ? join(res) : "-" * join(res)
end

## @lc code=end


text:
using Test

include("cas_infer.jl")

ex1 = :(cos(1 + 3.0) + 4 + (4-4im))
ex2 = :("ciao" * 2)
ex3 = :("ciao" * " mondo")

@test ComplexF64 == infer(ex1)
@test_throws MethodError infer(ex2)
@test String == infer(ex3)


text:
# https://github.com/JuliaLang/julia/pull/42585#discussion_r727423060
sample_stackpointer() = ccall("llvm.frameaddress", llvmcall, Ptr{Cvoid}, (Int,), 0)


text:
using PyPlot
using Statistics

using SubmodularMaximization

num_trials = 100

num_agents = 20
num_sensors = 6
nominal_area = 2.0

sensor_radius = sqrt(nominal_area / (num_agents * pi))
station_radius = 2 * sensor_radius

agent_specification = CircleAgentSpecification(sensor_radius, station_radius,
                                         num_sensors)

f(x) = mean_area_coverage(x, floor(Integer, 20 * num_agents / 4))

solvers = Any[]

partitions = [8]
for num_partitions in partitions
  solve_n(p) = solve_n_partitions(num_partitions, p)
  push!(solvers, (solve_n, "Partition-$num_partitions"))
end

push!(solvers, (solve_sequential, "Sequential"))

results = zeros(num_trials, length(solvers))

for trial_num in 1:num_trials
  println("Trial: $trial_num")

  agents = generate_agents(agent_specification, num_agents)
  problem = ExplicitPartitionProblem(f, agents)

  for solver_num in 1:length(solvers)
    solver = solvers[solver_num][1]
    name = solvers[solver_num][2]

    solution = solver(problem)

    results[trial_num, solver_num] = solution.value
  end
end

mean_value = mean(results; dims=2)

figure()
boxplot(results, notch=false, vert=false)
yticks(1:length(solvers), map(x->x[2], solvers))

nothing

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
# distritbutions
@rule typeof(-)(:in2, Marginalisation) (m_out::UnivariateNormalDistributionsFamily, m_in1::UnivariateNormalDistributionsFamily) = begin
    mout, vout = mean_var(m_out)
    min1, vin1 = mean_var(m_in1)
    return NormalMeanVariance(- mout + min1, vout + vin1)
end

@rule typeof(-)(:in2, Marginalisation) (m_out::MultivariateNormalDistributionsFamily, m_in1::MultivariateNormalDistributionsFamily) = begin
    mout, vout = mean_cov(m_out)
    min1, vin1 = mean_cov(m_in1)
    return MvNormalMeanCovariance(- mout + min1, vout + vin1)
end

# PointMass
@rule typeof(-)(:in2, Marginalisation) (m_out::PointMass, m_in1::PointMass) = PointMass(- mean(m_out) + mean(m_in1))

@rule typeof(-)(:in2, Marginalisation) (m_out::PointMass, m_in1::NormalMeanPrecision) = begin
    min1, win1 = mean_precision(m_in1)
    return NormalMeanPrecision(- mean(m_out) + min1, win1)
end

@rule typeof(-)(:in2, Marginalisation) (m_out::NormalMeanPrecision, m_in1::PointMass) = begin
    mout, wout = mean_precision(m_out)
    return NormalMeanPrecision(- mout + mean(m_in1), wout)
end

@rule typeof(-)(:in2, Marginalisation) (m_out::PointMass, m_in1::MvNormalMeanPrecision) = begin
    min1, win1 = mean_precision(m_in1)
    return MvNormalMeanPrecision(- mean(m_out) + min1, win1)
end

@rule typeof(-)(:in2, Marginalisation) (m_out::MvNormalMeanPrecision, m_in1::PointMass) = begin
    mout, wout = mean_precision(m_out)
    return MvNormalMeanPrecision(- mout + mean(m_in1), wout)
end

@rule typeof(-)(:in2, Marginalisation) (m_out::PointMass, m_in1::UnivariateNormalDistributionsFamily) = begin
    min1, vin1 = mean_var(m_in1)
    return NormalMeanVariance(- mean(m_out) + min1, vin1)
end

@rule typeof(-)(:in2, Marginalisation) (m_out::UnivariateNormalDistributionsFamily, m_in1::PointMass) = begin
    mout, vout = mean_var(m_out)
    return NormalMeanVariance(- mout + mean(m_in1), vout)
end

@rule typeof(-)(:in2, Marginalisation) (m_out::PointMass, m_in1::MultivariateNormalDistributionsFamily) = begin
    min1, vin1 = mean_cov(m_in1)
    return MvNormalMeanCovariance(- mean(m_out) + min1, vin1)
end

@rule typeof(-)(:in2, Marginalisation) (m_out::MultivariateNormalDistributionsFamily, m_in1::PointMass) = begin
    mout, vout = mean_cov(m_out)
    return MvNormalMeanCovariance(- mout + mean(m_in1), vout)
end


text:
############################### Reserve Variables` #########################################
"""
This function checks if the variables for reserves were created
"""
function check_activeservice_variables(
    container::OptimizationContainer,
    contributing_services::Vector{T},
) where {T <: PSY.Service}
    for service in contributing_services
        get_variable(
            container,
            ActivePowerReserveVariable(),
            typeof(service),
            PSY.get_name(service),
        )
    end
    return
end

################################## Reserve Requirement Constraint ##########################
"""
This function creates the requirement constraint that will be attained by the apropriate services
"""
function add_constraints!(
    container::OptimizationContainer,
    ::Type{RequirementConstraint},
    service::SR,
    contributing_services::Vector{<:PSY.Service},
    model::ServiceModel{SR, GroupReserve},
) where {SR <: PSY.StaticReserveGroup}
    time_steps = get_time_steps(container)
    service_name = PSY.get_name(service)
    add_constraints_container!(
        container,
        RequirementConstraint(),
        SR,
        [service_name],
        time_steps;
        meta=service_name,
    )
    constraint = get_constraint(container, RequirementConstraint(), SR, service_name)
    use_slacks = get_use_slacks(model)
    reserve_variables = [
        get_variable(container, ActivePowerReserveVariable(), typeof(r), PSY.get_name(r)) for r in contributing_services
    ]

    requirement = PSY.get_requirement(service)
    for t in time_steps
        resource_expression = JuMP.GenericAffExpr{Float64, JuMP.VariableRef}()
        for reserve_variable in reserve_variables
            JuMP.add_to_expression!(resource_expression, sum(reserve_variable[:, t]))
        end
        if use_slacks
            resource_expression += slack_vars[t]
        end
        constraint[service_name, t] =
            JuMP.@constraint(container.JuMPmodel, resource_expression >= requirement)
    end

    return
end


text:
#Demo 5: Generating DataRS & DataHR Datsets
#-------------------------------------------------------------------------------

using CMDimData
using CMDimData.MDDatasets
using CMDimData.EasyPlot

#==Attributes
===============================================================================#
LBL_AXIS_TIME = "Time (s)"
LBL_AXIS_AMPLITUDE = "Amplitude"

#==Input data
===============================================================================#
tfund = 1e-9 #Fundamental
osr = 20 #(fundamental) oversampling ratio
nfund = 20 #cycles of the fundamental

#==Computations
===============================================================================#
t = DataF1(0:(tfund/osr):(nfund*tfund))

#Generate parameter sweeps:
sweeplist = PSweep[
	PSweep("period", [1,2,3]*tfund)
	PSweep("amp", [1, 1.5, 2])
]

#Generate data:
tonesHR = fill(DataHR{DataF1}, sweeplist) do T, amp
	return amp*sin(t*(2pi/T))
end

tonesRS = DataRS(tonesHR)

#==Generate plot
===============================================================================#
plot = cons(:plot, nstrips=2, legend=true,
	ystrip1 = set(axislabel=LBL_AXIS_AMPLITUDE, striplabel="DataHR"),
	ystrip2 = set(axislabel=LBL_AXIS_AMPLITUDE, striplabel="DataRS"),
	xaxis = set(label=LBL_AXIS_TIME),
)
push!(plot,
	cons(:wfrm, tonesHR, label="tones", strip=1),
	cons(:wfrm, tonesRS, label="tones", strip=2),
)

pcoll = cons(:plot_collection, title="DataHR & DataMD", ncolumns=1)
	push!(pcoll, plot)

#==Return pcoll to user (call evalfile(...))
===============================================================================#
pcoll


text:
using StructTypes
using JSON3

StructTypes.StructType(::Type{SolSummary}) = StructTypes.Struct()
StructTypes.StructType(::Type{MutableResult}) = StructTypes.Struct()
StructTypes.StructType(::Type{Result}) = StructTypes.Struct()
StructTypes.StructType(::Type{Density}) = StructTypes.Struct()
StructTypes.StructType(::Type{Volume}) = StructTypes.Struct()
StructTypes.StructType(::Type{Options}) = StructTypes.Struct()

"""

```
save(R::Result, filename::String)
```

Function to write the result data structure to a json file.

"""
function save(R::Result, filename::String)
    f = open(filename, "w")
    JSON3.write(f, R)
    close(f)
end


text:
##############################################################################
#
# REFERENCES: Using definition from Devroye, IX.7
#  This definition differs from definitions on Wikipedia and other sources,
#  where distribution is over [0, 1] rather than [-1, 1].
#
##############################################################################

immutable Arcsine <: ContinuousUnivariateDistribution
end

function cdf(d::Arcsine, x::Number)
    if x < -1.0
        return 0.0
    elseif x > 1.0
        return 1.0
    else
        return (2.0 / pi) * asin(sqrt((x + 1.0) / 2.0))
    end
end

entropy(d::Arcsine) = -log(2.0) / pi

function insupport(d::Arcsine, x::Number)
    if -1.0 <= x <= 1.0
        return true
    else
        return false
    end
end

mean(d::Arcsine) = 0.0

median(d::Arcsine) = 0.0

function pdf(d::Arcsine, x::Number)
    if insupport(d, x)
        return 1.0 / (pi * sqrt(1.0 - x^2))
    else
        return 0.0
    end
end

quantile(d::Arcsine, p::Real) = 2.0 * sin((pi / 2.0) * p)^2 - 1.0

rand(d::Arcsine) = sin(2.0 * pi * rand())

skewness(d::Arcsine) = 0.0

var(d::Arcsine) = 1.0 / 2.0


text:
########################################################################
# Client
########################################################################\

using LinearAlgebra
using SparseArrays
using Random
using Printf

# Abstract class for client
abstract type AbstractClient end

# Client upload model to server
function update! end

# # Client download model to server
# function download end

include("fedProxClient.jl")
include("ScaffoldClient.jl")
include("fedDCDClient.jl")


text:
module problem118

description = """
Pandigital prime sets
Using all of the digits 1 through 9 and concatenating them freely to form decimal integers, different sets can be formed.
Interestingly with the set {2,5,47,89,631}, all of the elements belonging to it are prime.

How many distinct sets containing each of the digits one through nine exactly once contain only prime elements?
"""

# Convert an array of digits to an integer
function digitsToInt(digits :: Array)
  n = 0
  for d in digits
    n = n*10 + d
  end
  return n
end

# calculate the number of permutations of q that are prime
# Stores results in memo, for faster lookup for the same digits in future calls
function primePermutations(q :: Vector, memo :: Dict)
  total = 0
  sort!(q)
  key = digitsToInt(q)
  if haskey(memo, key)
    return memo[key]
  end

  for p in permutations(q)
    if isprime(digitsToInt(p))
      total += 1
    end
  end
  memo[key] = total
  return total
end

# Calculate the number of sets containing permutations of digits 1-9 that are all prime
# Julia's built in partitions, permutations and isprime do all the work :)
# For a given partition of [1:9] the number of sets is the product of prime permutations of each sub-set of digits.
function run()
  total = 0
  memo = (Int=>Int)[]
  for p in partitions([1:9])
    total += prod([primePermutations(q, memo) for q in p])
  end
  return total
end

end

max_stars_repo_name: mbuhot/mbuhot-euler-solutions

max_issues_repo_name: mbuhot/mbuhot-euler-solutions

max_forks_repo_name: mbuhot/mbuhot-euler-solutions


text:
module HenonMap
export henon_map

using Parameters: @with_kw, @unpack
using ..ExampleBase: LEDemo, DiscreteExample

@with_kw struct HenonMapParam
    a::Float64 = 1.4
    b::Float64 = 0.3
end

@inline function phase_dynamics!(u_next, u, p, t)
    @unpack a, b = p
    u_next[1] = 1 + u[2] - a * u[1]^2
    u_next[2] = b * u[1]
end

@inline @views function tangent_dynamics!(u_next, u, p, t)
    @unpack a, b = p
    phase_dynamics!(u_next[:, 1], u[:, 1], p, t)
    u_next[1, 2:end] .= u[2, 2:end] .- a * 2 * u[1, 1] .* u[1, 2:end]
    u_next[2, 2:end] .= b .* u[1, 2:end]
end

"""
Return a [`LEDemo`](@ref) for the Hénon map.

* M. Hénon, Commun. Math. Phys. Phys. 50, 69-77 (1976)
* <http://sprott.physics.wisc.edu/chaos/comchaos.htm>
* <https://en.wikipedia.org/wiki/H%C3%A9non_map>
"""
function henon_map(;
        u0=[-0.3332501594195585, 0.2743725483437717],
        t_renorm=10,
        t_attr=100000,
        atol=0, rtol=1e-2,
        terminator_options = [:atol => atol, :rtol => rtol],
        kwargs...)
    LEDemo(DiscreteExample(
        "Hénon map",
        phase_dynamics!, u0, t_renorm, HenonMapParam(),
        tangent_dynamics!,
        t_attr,
        [0.41922, -1.62319],   # known_exponents
        atol, rtol,
        terminator_options,
    ); kwargs...)
end

end


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

# assert that all columns of table are Continuous
function assert_continuous(table)
  types = schema(table).scitypes
  @assert all(T <: Continuous for T in types) "columns must hold continuous variables"
end


text:
# Translating
# https://proplot.readthedocs.io/en/latest/subplots.html#Shared-and-spanning-labels
# to Julia

# Import proplot using PyCall
using PyCall
proplot = pyimport("proplot")

N = 50
M = 40
colors = proplot.Colors("grays_r", M, left=0.1, right=0.8)
datas = []
for scale in (1, 3, 7, 0.2)
    data = scale * cumsum(rand(N, M) .- 0.5, dims=1)[N÷2+1:end, :]
    push!(datas, data)
end

# Same plot with different sharing and spanning settings
for share in (0, 1, 2, 3)
    f, axs = proplot.subplots(
        ncols=4, aspect=1, axwidth=1.2,
        sharey=share, spanx=share≥2
    )
    for (ax, data) in zip(axs, datas)
        ax.plot(data, cycle=colors)
        ax.format(
            suptitle="Axis-sharing level: $share, spanning labels $(share≥2 ? :on : :off)",
            grid=false, xlabel="spanning", ylabel="shared"
        )
    end
    f.savefig("figs/Shared-and-spanning-labels-1$(Char(97+share)).svg", transparent=false)
end

max_stars_repo_name: natgeo-wong/JuliaProPlotExamples

max_issues_repo_name: natgeo-wong/JuliaProPlotExamples

max_forks_repo_name: natgeo-wong/JuliaProPlotExamples


text:
export PointRobotNDBoxes

# ---------- Point Robot (amongst N-d boxes) ----------

struct BoxBounds{N,T<:AbstractFloat}
    lo::SVector{N,T}
    hi::SVector{N,T}
end
BoxBounds(lo::AbstractVector, hi::AbstractVector) = BoxBounds(SVector(lo), SVector(hi))
BoxBounds(lohi::Matrix) = BoxBounds(lohi[:,1], lohi[:,2])
inflate{N,T}(BB::BoxBounds{N,T}, ε) = ε > 0 ? BoxBounds(BB.lo - T(ε), BB.hi + T(ε)) : BB
changeprecision{T<:AbstractFloat}(::Type{T}, BB::BoxBounds) =
    BoxBounds(changeprecision(T, BB.lo), changeprecision(T, BB.hi))

mutable struct PointRobotNDBoxes{N,T} <: SweptCollisionChecker
    boxes::Vector{BoxBounds{N,T}}
    count::Int
end
PointRobotNDBoxes{N,T}(boxes::Vector{BoxBounds{N,T}}) = PointRobotNDBoxes(boxes, 0)
PointRobotNDBoxes{T}(box_list::Vector{Matrix{T}}) = PointRobotNDBoxes(map(BoxBounds, box_list))
PointRobotNDBoxes{T}(boxes::Array{T,3}) = PointRobotNDBoxes(vec(mapslices(BoxBounds, boxes, [1,2])))
PointRobotNDBoxes(boxhcat::Matrix) = PointRobotNDBoxes(reshape(boxhcat, size(boxhcat, 1), 2, div(size(boxhcat, 2), 2)))
changeprecision{T<:AbstractFloat}(::Type{T}, CC::PointRobotNDBoxes) = PointRobotNDBoxes(changeprecision(T, CC.boxes))

is_free_state(v::AbstractVector, CC::PointRobotNDBoxes) = is_free_state(v, CC.boxes)
is_free_motion(v::AbstractVector, w::AbstractVector, CC::PointRobotNDBoxes) = (CC.count += 1;
                                                                               is_free_motion(v, w, CC.boxes))
is_free_path(P::Path, CC::PointRobotNDBoxes) = is_free_path(P, CC.boxes)

inflate{N,T}(CC::PointRobotNDBoxes{N,T}, ε) = ε > 0 ? PointRobotNDBoxes([inflate(B, T(ε)) for B in CC.boxes]) : CC
addobstacle(CC::PointRobotNDBoxes, o) = PointRobotNDBoxes(vcat(CC.boxes, BoxBounds(o)))
addblocker(CC::PointRobotNDBoxes, v::AbstractVector, r) = addobstacle(CC, BoxBounds(v-r, v+r))
closest(p::AbstractVector, CC::PointRobotNDBoxes, W::AbstractMatrix) = closest(p, CC.boxes, W)
closeR(p::AbstractVector, CC::PointRobotNDBoxes, W::AbstractMatrix, r2) = closeR(p, CC.boxes, W, r2)

plot(CC::PointRobotNDBoxes, lo = zeros(2), hi = ones(2); kwargs...) =
    map(o -> plot_rectangle(o.lo, o.hi, color = "red", edgecolor = "none",
                            xmin = lo[1], xmax = hi[1], ymin = lo[2], ymax = hi[2]; kwargs...), CC.boxes)

### Point/box obstacle checking in R^n (side note: SAT-style might be faster)

is_free_state{N}(v::AbstractVector, BB::BoxBounds{N}) = @any [!(BB.lo[i] <= v[i] <= BB.hi[i]) for i in 1:N]
is_free_state{N,T}(v::AbstractVector, BL::Vector{BoxBounds{N,T}}) = @all [is_free_state(v, BL[k]) for k in 1:length(BL)]
is_free_motion_broadphase{N}(l::AbstractVector, h::AbstractVector, BB::BoxBounds{N}) =
    @any [BB.hi[i] < l[i] || BB.lo[i] > h[i] for i in 1:N]
function is_free_motion{N}(v::AbstractVector, w::AbstractVector, BB::BoxBounds{N})
    v_to_w = w - v
    corner = blend(map(<, v, BB.lo), BB.lo, BB.hi)
    lambdas = (corner - v) ./ v_to_w
    !(@any [(@all [i == j || BB.lo[j] <= v[j] + v_to_w[j]*lambdas[i] <= BB.hi[j] for j in 1:N]) for i in 1:N])
end
function is_free_motion{N,T}(v::AbstractVector, w::AbstractVector, BL::Vector{BoxBounds{N,T}})
    bb_min = map(min, v, w)    # TODO: check back when StaticArrays has implemented min(v,w)
    bb_max = map(max, v, w)
    @all [is_free_motion_broadphase(bb_min, bb_max, BL[k]) || is_free_motion(v, w, BL[k]) for k in 1:length(BL)]
end
is_free_path{B<:BoxBounds}(P::Path, BL::Vector{B}) = @all [is_free_motion(P[i], P[i+1], BL) for i in 1:length(BL)-1]

### Closest Point

function closest(p::AbstractVector, BB::BoxBounds, W::AbstractMatrix)
    # v = Variable(length(p))
    # problem = minimize(quad_form(v-p, W), view(o,:,1) <= v, v <= view(o,:,2))
    # solve!(problem, ECOS.ECOSSolver(verbose=false))
    # return (problem.optval, vec(v.value))
    L = chol(W)
    vmin = bvls(L, L*p, BB.lo, BB.hi)
    d2min = dot(vmin - p, W*(vmin - p))
    d2min, vmin
end

function closest{N,T}(p::AbstractVector, BL::Vector{BoxBounds{N,T}}, W::AbstractMatrix)
    d2min, vmin = T(Inf), p
    for k = 1:length(BL)
        (d2, v) = closest(p, BL[k], W)
        if d2 < d2min
            d2min, vmin = d2, v
        end
    end
    d2min, vmin
end

function closeR{N,T}(p::AbstractVector, BL::Vector{BoxBounds{N,T}}, W::AbstractMatrix, r2)
    cps = [closest(p, BL[k], W) for k in 1:length(BL)]
    sort!(cps[[c[1] for c in cps] .< r2], by=first)
end


text:
function niche_network(S::Int, C::Float64)

    cond = false
    adj = fill(0, S, S)
    n = fill(0.0, S)
    r = fill(0.0, S)
    c = fill(0.0, S)

    while !cond
        n = sort(rand(S))
        β = 1 / (2 * C) - 1
        r = n .* rand(Beta(1, β), S)

        # set r[position of min[η]] = 0, so that every web has at least one basal species
        r[1] = 0.0

        for pred in 1:S

            c[pred] = rand(Uniform(r[pred]/2, n[pred]))

            for prey in 1:S
                if ((c[pred] - r[pred] / 2) <= n[prey] <= (c[pred] + r[pred] / 2))
                    adj[prey, pred] = 1
                else
                    adj[prey, pred] = 0
                end
            end
        end

        cond = is_connected(DiGraph(adj))
    end

    return (fw = adj, n = n, r = r, c = c)
end

function prob_network(S, C)

    cond = false
    adj = fill(0, S, S)

    n = fill(0.0, S)
    r = fill(0.0, S)
    c = fill(0.0, S)

    while !cond

        n = sort(rand(S))
        β = 1 / (2 * C) - 1
        r = n .* rand(Beta(1, β), S)

        # set r[position of min[η]] = 0, so that every web has at least one basal species
        r[1] = 0.0

        @inbounds for pred in 1:S

            c[pred] = rand(Uniform(r[pred] / 2, n[pred]))

            for prey in 1:S
                ϕ = exp(-((n[prey] - c[pred]) / (r[pred] / 2)) ^ 2)
                adj[prey, pred] = rand(Bernoulli(ϕ), 1)[1]
            end
        end

        cond = is_connected(DiGraph(adj))
    end

    return (fw = adj, n = n, r = r, c = c)
end

function ϕ_mat(θ, m)

    ϕ = fill(0.0, m.S, m.S)

    @inbounds for pred in 1:m.nonbasal
        for prey in 1:m.S
            ϕ[prey, m.basal + pred] = 0.99 * exp(-((θ.n[prey] - θ.c[pred] * θ.n[m.basal + pred]) / (θ.r[pred] / 2)) ^ 2)
        end
    end

    return ϕ

end

max_stars_repo_name: clint-leach/niche-model

max_issues_repo_name: clint-leach/niche-model

max_forks_repo_name: clint-leach/niche-model


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

@with_kw mutable struct MaximumIterationParameters
    maximum::Int = 1000
end

"""
    MaximumIterations

Functor object for using maximum number of iterations as termination criterion in a quasigradient algorithm. Create by supplying a [`AfterMaximumIterations`](@ref) object through `terminate` to `QuasiGradient.Optimizer` or by setting the [`Termination`](@ref) attribute.

...
# Parameters
- `maximum::Integer = 1000`: Maximum number of iterations
...
"""
struct MaximumIterations <: AbstractTerminationCriterion
    parameters::MaximumIterationParameters

    function MaximumIterations(; kw...)
        return new(MaximumIterationParameters(; kw...))
    end
end

function Progress(termination::MaximumIterations, str::AbstractString)
    return Progress(termination.parameters.maximum, 0.0, str)
end

function progress_value(::MaximumIterations, k::Integer, f::AbstractFloat, ∇f_norm::AbstractFloat)
    return k
end

function terminate(termination::MaximumIterations, k::Integer, f::Float64, x::AbstractVector, ∇f::AbstractVector)
    return k >= termination.parameters.maximum
end

# API
# ------------------------------------------------------------
"""
    AfterMaximumIterations

Factory object for [`MaximumIterations`](@ref). Pass to `terminate` in `Quasigradient.Optimizer` or set the [`Termination`](@ref) attribute. See ?MaximumIterations for parameter descriptions.

"""
struct AfterMaximumIterations <: AbstractTermination
    parameters::MaximumIterationParameters
end
AfterMaximumIterations(maximum::Integer) = AfterMaximumIterations(MaximumIterationParameters(; maximum = maximum))
AfterMaximumIterations(; kw...) = AfterMaximumIterations(MaximumIterationParameters(; kw...))

function (criteria::AfterMaximumIterations)(::Type{T}) where T <: AbstractFloat
    return MaximumIterations(; type2dict(criteria.parameters)...)
end


text:
import NCDatasets

filename = tempname()

dimlen = 10

T = Int32
data = Vector{Vector{T}}(undef,dimlen)
for i = 1:length(data)
    data[i] = T.(collect(1:i) .+ 100 * i)
end

ncdata = Vector{NCDatasets.nc_vlen_t{T}}(undef,dimlen)

for i = 1:length(data)
    ncdata[i] = NCDatasets.nc_vlen_t{T}(length(data[i]), pointer(data[i]))
end

varname = "varname"
vlentypename = "name-vlen"

ncid = NCDatasets.nc_create(filename,NCDatasets.NC_CLOBBER | NCDatasets.NC_NETCDF4)

dimid = NCDatasets.nc_def_dim(ncid, "casts", dimlen)

typeid = NCDatasets.nc_def_vlen(ncid, vlentypename, NCDatasets.ncType[T])
varid = NCDatasets.nc_def_var(ncid, varname, typeid, [dimid])

for i = 1:dimlen
    NCDatasets.nc_put_var1(ncid, varid, [i-1], data[i])
end

typeids = NCDatasets.nc_inq_typeids(ncid)
typeid = typeids[1]
name2,datum_size2,base_nc_type2 = NCDatasets.nc_inq_vlen(ncid,typeid)

@test name2 == vlentypename
@test NCDatasets.jlType[base_nc_type2] == T

NCDatasets.nc_close(ncid)

# Reopen

ncid = NCDatasets.nc_open(filename,NCDatasets.NC_NOWRITE)

varid = NCDatasets.nc_inq_varid(ncid,varname)

typeids = NCDatasets.nc_inq_typeids(ncid)
#@show typeids

xtype = NCDatasets.nc_inq_vartype(ncid,varid)

@test xtype == typeids[1]

if xtype >= NCDatasets.NC_FIRSTUSERTYPEID
    #@show xtype,NCDatasets.NC_VLEN

    typename,shape,base_nc_type,numfields,class = NCDatasets.nc_inq_user_type(ncid,xtype)

    #@show typename,shape,base_nc_type,numfields,class
    @test base_nc_type == NCDatasets.NC_INT

    T2 = NCDatasets.jlType[base_nc_type]

    @test T == T2
    if class == NCDatasets.NC_VLEN
        data2 = Vector{Vector{T}}(undef,dimlen)
        NCDatasets.nc_get_var!(ncid,varid,data2)
        @test data == data2

        for i = 1:dimlen
            tmp2 = NCDatasets.nc_get_var1(Vector{T},ncid,varid,[i-1])
            @test data[i] == tmp2
        end
    end
end

NCDatasets.nc_close(ncid)


text:
using Pkg; Pkg.instantiate()

using GeoStats
using DensityRatioEstimation
using LossFunctions
using ProgressMeter
using DataFrames
using DataDeps
using MLJ, CSV
using Statistics
using Random

# reproducible results
Random.seed!(2020)

# estimators of generalization error
error_cv( m, p, k, ℒ) = error(PointwiseLearn(m), p, CrossValidation(k, loss=ℒ))
error_bcv(m, p, r, ℒ) = error(PointwiseLearn(m), p, BlockCrossValidation(r, loss=ℒ))
error_drv(m, p, k, ℒ) = error(PointwiseLearn(m), p, DensityRatioValidation(k, loss=ℒ, estimator=LSIF(σ=2.0,b=10)))

# actual error (known labels)
function error_empirical(m, p, ℒ)
  t  = task(p)
  Ωs = sourcedata(p)
  Ωt = targetdata(p)

  # learn task on source domain and perform
  # it on both source and target domains
  lm = learn(t, Ωs, m)
  ŷs = perform(t, Ωs, lm)
  ŷt = perform(t, Ωt, lm)

  # error on source
  ϵs = map(outputvars(t)) do var
    var => LossFunctions.value(ℒ[var], Ωs[var], ŷs[var], AggMode.Mean())
  end

  # error on target
  ϵt = map(outputvars(t)) do var
    var => LossFunctions.value(ℒ[var], Ωt[var], ŷt[var], AggMode.Mean())
  end

  Dict(ϵs), Dict(ϵt)
end

function experiment(m, p, r, k, ℒ)
  # try different error estimates
  cv  = error_cv( m, p, k, ℒ)
  bcv = error_bcv(m, p, r, ℒ)
  drv = error_drv(m, p, k, ℒ)

  # actual error (unhide labels)
  ϵs, ϵt = error_empirical(m, p, ℒ)

  # model name without suffix
  model = replace(info(m).name, r"(.*)(Regressor|Classifier)" => s"\g<1>")

  map(outputvars(task(p))) do var
    (MODEL=model, SOURCE=ϵs[var], TARGET=ϵt[var],
     CV=cv[var], BCV=bcv[var], DRV=drv[var])
  end
end

# -------------
# MAIN SCRIPT
# -------------

# download dataset if needed
register(DataDep("NewZealand",
         "Taranaki Basin Curated Well Logs",
         "https://zenodo.org/record/3832955/files/taranaki-basin-curated-well-logs.tar.gz",
         "608f7aad5a4e9fded6441fd44f242382544d3f61790446175f5ede83f15f4d11",
         post_fetch_method=DataDeps.unpack))

# name of the CSV file in the dataset
csv = joinpath(datadep"NewZealand","taranaki-basin-curated-well-logs","logs.csv")

# logs used in the experiment
logs = [:GR,:SP,:DENS,:NEUT,:DTC]

# read/clean raw data
df = CSV.read(csv)
df = df[:,[logs...,:X,:Y,:Z,:FORMATION,:ONSHORE]]
dropmissing!(df)
categorical!(df, :FORMATION)
categorical!(df, :ONSHORE)
for log in logs
  x = df[!,log]
  μ = mean(x)
  σ = std(x, mean=μ)
  df[!,log] .= (x .- μ) ./ σ
end

# define spatial data
wells = GeoDataFrame(df, [:X,:Y,:Z])

# select the two most frequent formations
formations = groupby(wells, :FORMATION)
frequency = sortperm(npoints.(formations), rev=true)
𝒞 = DataCollection(formations[frequency[1:2]])

# eliminate duplicate coordinates
Ω = uniquecoords(𝒞)

# split onshore (True) vs. offshore (False)
onoff = groupby(Ω, :ONSHORE)
order = sortperm(onoff[:values], rev=true)
Ωs, Ωt = onoff[order]

# we are left with two formations onshore and offshore
# make sure that these two are balanced for classification
fs, ft = Ωs[:FORMATION], Ωt[:FORMATION]

# formation counts
ms = count(isequal("Urenui"), fs)
ns = count(isequal("Manganui"), fs)
mt = count(isequal("Urenui"), ft)
nt = count(isequal("Manganui"), ft)

# formation proportions
ps = ms / (ms + ns)
pt = mt / (mt + nt)

# weighted sampling
ws = [f == "Urenui" ? 0.5/ps : 0.5/(1-ps) for f in fs]
wt = [f == "Urenui" ? 0.5/pt : 0.5/(1-pt) for f in ft]
Ωs = sample(Ωs, 300000, ws, replace=false)
Ωt = sample(Ωt,  50000, wt, replace=false)

# drop levels to avoid known downstream issues in MLJ
fs, ft = Ωs[:FORMATION], Ωt[:FORMATION]
levels!(fs, ["Urenui","Manganui"])
levels!(ft, ["Urenui","Manganui"])
𝒫s = georef(OrderedDict(:FORMATION => fs), domain(Ωs))
𝒫t = georef(OrderedDict(:FORMATION => ft), domain(Ωt))
Ωs = join(view(Ωs, logs), 𝒫s)
Ωt = join(view(Ωt, logs), 𝒫t)

# additional configuration without shift
𝒞 = DataCollection(Ωs, Ωt)
fraction = npoints(Ωs) / (npoints(Ωs) + npoints(Ωt))
Γs, Γt = split(𝒞, fraction)

# -------------------
# PROBLEM DEFINITION
# -------------------
# predict formation from well logs
t = ClassificationTask(logs, :FORMATION)

# onshore -> offshore problem
p = LearningProblem(Ωs, Ωt, t)

# problem without shift
q = LearningProblem(Γs, Γt, t)

# -----------
# EXPERIMENT
# -----------
@load DummyClassifier pkg="ScikitLearn"
@load RidgeClassifier pkg="ScikitLearn"
@load LogisticClassifier pkg="ScikitLearn"
@load KNeighborsClassifier pkg="ScikitLearn"
@load GaussianNBClassifier pkg="ScikitLearn"
@load BayesianLDA pkg="ScikitLearn"
@load PerceptronClassifier pkg="ScikitLearn"
@load DecisionTreeClassifier pkg="DecisionTree"

# list of models
mrange = [RidgeClassifier(), LogisticClassifier(), KNeighborsClassifier(),
          GaussianNBClassifier(), BayesianLDA(), PerceptronClassifier(),
          DecisionTreeClassifier(), DummyClassifier()]

# block sides and number of folds for error estimators
r = (10000., 10000., 500.)
k = length(GeoStats.partition(Ωs, BlockPartitioner(r)))

# misclassification loss
ℒ = Dict(:FORMATION => MisclassLoss())

# experiment iterator and progress
iterator  = Iterators.product(mrange)
pprogress = Progress(length(iterator), "ONSHORE → OFFSHORE ")
qprogress = Progress(length(iterator), "NO COVARIATE SHIFT ")

# return missing in case of failure
skip = e -> (println("Skipped: $e"); missing)

# perform experiments
presults = progress_pmap(iterator, progress=pprogress, on_error=skip) do (m,)
  experiment(m, p, r, k, ℒ)
end
qresults = progress_pmap(iterator, progress=qprogress, on_error=skip) do (m,)
  experiment(m, q, r, k, ℒ)
end

# merge all results into a single table
pres = Iterators.flatten(skipmissing(presults))
qres = Iterators.flatten(skipmissing(qresults))
pres = [(SHIFT="YES", r...) for r in pres]
qres = [(SHIFT="NO",  r...) for r in qres]
ares = DataFrame(vcat(pres, qres))

# save all results to disk
fname = joinpath(@__DIR__,"results","newzealand.csv")
CSV.write(fname, ares)

max_issues_repo_name: IBM/geostats-gen-error

max_forks_repo_name: IBM/geostats-gen-error


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

using Test, Tensors, LinearAlgebra

# Kronecker delta
@test delta(1, 1) == 1
@test delta(1, 2) == 0
@test_throws MethodError delta(1.0, 2.0)  # indices must be integers
@test_throws MethodError delta(1, BigInt(2))  # both must have the same type
@test delta(1, 2) isa Int  # the output type matches the input
@test delta(BigInt(1), BigInt(2)) isa BigInt

# Various tensors
let Z3 = zeros(3, 3),
    O3 = ones(3, 3)
    @test isapprox(tovoigt(II()), I)
    @test isapprox(tovoigt(IT()), [I   Z3;
                                   Z3  Z3])

    @test isapprox(tovoigt(IS()), [I   Z3;
                                   Z3  1//2*I])
    @test isapprox(tovoigt(IA()), [Z3  Z3;
                                   Z3  1//2*I])

    @test isapprox(tovoigt(IV()), [1//3*O3  Z3;
                                   Z3       Z3])
    @test isapprox(tovoigt(ID()), [(I - 1//3*O3)  Z3;
                                   Z3             1//2*I])

    @test let lambda = 10.0,
              mu = 1.0
        isapprox(tovoigt(isotropic_elasticity_tensor(lambda, mu)), [(lambda*O3 + 2*mu*I)  Z3;
                                                                    Z3                    mu*I])
    end
end

# Lamé parameters for isotropic solids
@test all(isapprox(result, expected)
        for (result, expected) in zip(lame(1e11, 0.3), (5.769230769230769e10, 3.846153846153846e10)))
@test all(isapprox(result, expected)
        for (result, expected) in zip(delame(lame(1e11, 0.3)...), (1e11, 0.3)))

# Mutating function to non-mutating function conversion
let  # introduce a local scope so the name `f!` is only defined locally for this test.
    function f!(out, x)
        out[:] = [sin(elt) for elt in x]
        return nothing
    end

    let
        out = [0.0]
        @test all([f!(out, [pi/4]) == nothing,
                 isapprox(out, [1/sqrt(2)])])
    end

    let
        out = [0.0]
        f = debang(f!)
        @test f isa Function
        @test all([isapprox(f([pi/4]), [1/sqrt(2)]),
                 out == [0.0]])
    end
end

# Newton root finder
let g(x) = [(1 - x[1]^2) + x[2]],
    x0 = [0.8, 0.2]
    @test !isapprox(g(x0), [0.0], atol=1e-15)  # find_root should have to actually do something
    @test isapprox(g(find_root(g, x0)), [0.0], atol=1e-15)
end


text:
using Pages
using Test

methods = [GET,HEAD,POST,PUT,DELETE,CONNECT,OPTIONS,TRACE,PATCH]
for m in methods
    Endpoint("/hello",m) do request::HTTP.Request
        string(m)
    end
    e = Pages.endpoints["/hello"]
    @test HTTP.handle(e.handlers[Pages.method(m)],HTTP.Request()) === string(m)
end


text:
##############################################################################
#
# VortexFilament.jl
#
# Part of CVortex.jl
# Representation of a vortex filament.
#
# Copyright 2019-2020 HJA Bird
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
# IN THE SOFTWARE.
##############################################################################

"""
CVortex internal representation of a straight singular vortex filament

The filament starts at coord1, ends at coord2 and has vorticity per unit
length of vorticity_density.

You only need to use this if you plan on calling CVortex's underlying
library directly.
"""
struct VortexFilament
    coord1 :: Vec3f
    coord2 :: Vec3f
    vorticity_density :: Float32
end

function VortexFilament(
    coord1::Vector{<:Real}, coord2::Vector{<:Real}, vort::Real)
    return VortexFilament(Vec3f(coord1), Vec3f(coord2), Float32(vort))
end

"""
    filament_induced_velocity(
        filament_start_coord :: Vector{<:Real},
        filament_end_coord :: Vector{<:Real},
        filament_strength :: Real,
        measurement_point :: Vector{<:Real})

    filament_induced_velocity(
        filament_start_coords :: Matrix{<:Real},
        filament_end_coords :: Matrix{<:Real},
        filament_strengths :: Vector{<:Real},
        measurement_point :: Vector{<:Real})

    filament_induced_velocity(
        filament_start_coords :: Matrix{<:Real},
        filament_end_coords :: Matrix{<:Real},
        filament_strengths :: Vector{<:Real},
        measurement_points :: Matrix{<:Real})

Compute the velocity induced in the flow field by vortex filaments. The third
multiple-multiple method may be GPU accelerated.

# Arguments
- `filament_start_coord` : The start coordinates of vortex filaments
- `filament_end_coord` : The end coordinates of vortex filaments
- `filament_strengths` : The vorticity per unit length of vortex filaments.
- `measurement_points` : The points where induced velocity is measured.

Vector arguments are expected to have length 3. Filament matrix arguments are
expected to have size N by 3. Measurement matrix arguments are expected to
have size M by 3. Returns an M by 3 matrix representing velocities.
"""
function filament_induced_velocity(
    filament_start_coord :: Vector{<:Real},
    filament_end_coord :: Vector{<:Real},
    filament_strength :: Real,
    measurement_point :: Vector{<:Real})

    check_filament_definition(
        filament_start_coord, filament_end_coord, filament_strength)
    convertable_to_Vec3f_vect(measurement_point, "measurement_point")

    inducing_filament = VortexFilament(filament_start_coord, filament_end_coord,
        filament_strength)
    mes_pnt = Vec3f(measurement_point)

    ret = Vec3f(0., 0., 0.)
    #=
    CVTX_EXPORT bsv_V3f cvtx_F3D_S2S_vel(
        const cvtx_F3D *self,
        const bsv_V3f mes_point);
    =#
    ret = ccall(
            ("cvtx_F3D_S2S_vel", libcvortex),
            Vec3f,
            (Ref{VortexFilament}, Vec3f),
            inducing_filament, mes_pnt)
    return Vector{Float32}(ret)
end

function filament_induced_velocity(
    filament_start_coords :: Matrix{<:Real},
    filament_end_coords :: Matrix{<:Real},
    filament_strengths :: Vector{<:Real},
    measurement_point :: Vector{<:Real})

    check_filament_definition(
        filament_start_coords, filament_end_coords, filament_strengths)
    convertable_to_Vec3f_vect(measurement_point, "measurement_point")

    ni = size(filament_end_coords)[1]
    inducing_filaments = map(
        i->VortexFilament(filament_start_coords[i,:], filament_end_coords[i, :],
        filament_strengths[i]),
        1:ni)
    mes_pnt = Vec3f(measurement_point)

    pargarr = Vector{Ptr{VortexFilament}}(undef, ni)
    for i = 1 : ni
        pargarr[i] = Base.pointer(inducing_filaments, i)
    end
    ret = Vec3f(0., 0., 0.)
	GC.@preserve pargarr inducing_filaments begin
		#=
		CVTX_EXPORT bsv_V3f cvtx_F3D_M2S_vel(
			const cvtx_F3D **array_start,
			const int num_filaments,
		=#
		ret = ccall(
				("cvtx_F3D_M2S_vel", libcvortex),
				Vec3f,
				(Ref{Ptr{VortexFilament}}, Cint, Vec3f),
				pargarr, ni, mes_pnt
				)
	end
    return Vector{Float32}(ret)
end

function filament_induced_velocity(
    filament_start_coords :: Matrix{<:Real},
    filament_end_coords :: Matrix{<:Real},
    filament_strengths :: Vector{<:Real},
    measurement_points :: Matrix{<:Real})

    check_filament_definition(
        filament_start_coords, filament_end_coords, filament_strengths)
    convertable_to_Vec3f_vect(measurement_points, "measurement_points")

    ni = size(filament_end_coords)[1]
    np = size(measurement_points)[1]
    inducing_filaments = map(
        i->VortexFilament(filament_start_coords[i,:], filament_end_coords[i, :],
            filament_strengths[i]),
        1:ni)
    mes_pnts = map(i->Vec3f(measurement_points[i, :]), 1:np)

    pargarr = Vector{Ptr{VortexFilament}}(undef, ni)
    for i = 1 : length(pargarr)
        pargarr[i] = Base.pointer(inducing_filaments, i)
    end
    ret = Vector{Vec3f}(undef, np)
	GC.@preserve pargarr inducing_filaments begin
		#=
		CVTX_EXPORT void cvtx_F3D_M2M_vel(
			const cvtx_F3D **array_start,
			const int num_filaments,
			const bsv_V3f *mes_start,
			const int num_mes,
			bsv_V3f *result_array);
		=#
		ccall(
			("cvtx_F3D_M2M_vel", libcvortex),
			Cvoid,
			(Ptr{Ptr{VortexFilament}}, Cint, Ptr{Vec3f},
				Cint, Ref{Vec3f}),
			pargarr, ni, mes_pnts, np, ret)
	end
    return Matrix{Float32}(ret)
end

"""
    filament_induced_dvort(
        filament_start_coord :: Vector{<:Real},
        filament_end_coord :: Vector{<:Real},
        filament_strength :: Real,
        induced_particle_position :: Vector{<:Real},
        induced_particle_vorticity :: Vector{<:Real})

    filament_induced_dvort(
        filament_start_coords :: Matrix{<:Real},
        filament_end_coords :: Matrix{<:Real},
        filament_strengths :: Vector{<:Real},
        induced_particle_position :: Vector{<:Real},
        induced_particle_vorticity :: Vector{<:Real})

    filament_induced_dvort(
        filament_start_coords :: Matrix{<:Real},
        filament_end_coords :: Matrix{<:Real},
        filament_strengths :: Vector{<:Real},
        induced_particle_position :: Matrix{<:Real},
        induced_particle_vorticity :: Matrix{<:Real})

Compute the rate of change of vorticity of vortex particles induced by
vortex filaments. Multiple-multiple method may be GPU accelerated. Modelled
as singular vortex particles and filaments.

# Arguments
- `filament_start_coord` : The start coordinates of vortex filaments
- `filament_end_coord` : The end coordinates of vortex filaments
- `filament_strengths` : The vorticity per unit length of vortex filaments.
- `induced_particle_position` : The positions of vortex particles.
- `induced_particle_vorticity` : The vorticities of vortex particles

Vector arguments are expected to have length 3. Filament matrix arguments are
expected to have size N by 3. Vortex particles matrix arguments are expected to
have size M by 3. Returns an M by 3 matrix representing particle vorticity
derivatives.
"""
function filament_induced_dvort(
    filament_start_coord :: Vector{<:Real},
    filament_end_coord :: Vector{<:Real},
    filament_strength :: Real,
    induced_particle_position :: Vector{<:Real},
    induced_particle_vorticity :: Vector{<:Real})

    check_filament_definition(
        filament_start_coord, filament_end_coord, filament_strength)
    check_particle_definition_3D(induced_particle_position,
        induced_particle_vorticity)

    inducing_filament = VortexFilament(filament_start_coord, filament_end_coord,
        filament_strength)
    induced_particle = VortexParticle3D(
        induced_particle_position, induced_particle_vorticity, 0.0)

    ret = Vec3f(0., 0., 0.)
    #=
    CVTX_EXPORT bsv_V3f cvtx_F3D_S2S_dvort(
        const cvtx_F3D *self,
        const cvtx_P3D *induced_particle);
    =#
    ret = ccall(
            ("cvtx_F3D_S2S_dvort", libcvortex),
            Vec3f,
            (Ref{VortexFilament}, Ref{VortexParticle3D}),
            inducing_filament, induced_particle
            )
    return Vector{Float32}(ret)
end

function filament_induced_dvort(
    filament_start_coords :: Matrix{<:Real},
    filament_end_coords :: Matrix{<:Real},
    filament_strengths :: Vector{<:Real},
    induced_particle_position :: Vector{<:Real},
    induced_particle_vorticity :: Vector{<:Real})

    check_filament_definition(
        filament_start_coords, filament_end_coords, filament_strengths)
    check_particle_definition_3D(induced_particle_position,
        induced_particle_vorticity)

    ni = size(filament_end_coords)[1]
    inducing_filaments = map(
        i->VortexFilament(filament_start_coords[i,:], filament_end_coords[i, :],
        filament_strengths[i]),
        1:ni)
    induced_particle = VortexParticle3D(
        induced_particle_position, induced_particle_vorticity, 0.0)

    pargarr = Vector{Ptr{VortexFilament}}(undef, length(inducing_filaments))
    for i = 1 : length(pargarr)
        pargarr[i] = Base.pointer(inducing_filaments, i)
    end
	GC.@preserve pargarr inducing_filaments begin
		#=
		CVTX_EXPORT bsv_V3f cvtx_F3D_M2S_dvort(
			const cvtx_F3D **array_start,
			const int num_filaments,
			const cvtx_P3D *induced_particle);
		=#
		ret = ccall(
				("cvtx_F3D_M2S_dvort", libcvortex),
				Vec3f,
				(Ref{Ptr{VortexFilament}}, Cint, Ref{VortexParticle3D}),
				pargarr, length(inducing_filaments), induced_particle
				)
	end
    return Vector{Float32}(ret)
end

function filament_induced_dvort(
    filament_start_coords :: Matrix{<:Real},
    filament_end_coords :: Matrix{<:Real},
    filament_strengths :: Vector{<:Real},
    induced_particle_position :: Matrix{<:Real},
    induced_particle_vorticity :: Matrix{<:Real})

    check_filament_definition(
        filament_start_coords, filament_end_coords, filament_strengths)
    check_particle_definition_3D(induced_particle_position,
        induced_particle_vorticity)

    ni = size(filament_end_coords)[1]
    np = size(induced_particle_position)[1]
    inducing_filaments = map(
        i->VortexFilament(filament_start_coords[i,:], filament_end_coords[i, :],
        filament_strengths[i]),
        1:ni)
    induced_particles = map(
        i->VortexParticle3D(
            induced_particle_position[i, :],
            induced_particle_vorticity[i, :], 0.0),
        1:np)

    pargarr = Vector{Ptr{VortexFilament}}(undef, ni)
    for i = 1 : length(pargarr)
        pargarr[i] = Base.pointer(inducing_filaments, i)
    end
    indarg = Vector{Ptr{VortexParticle3D}}(undef, np)
    for i = 1 : length(indarg)
        indarg[i] = Base.pointer(induced_particles, i)
    end
    ret = Vector{Vec3f}(undef, length(induced_particles))
	GC.@preserve pargarr inducing_filaments indarg induced_particles begin
		#=
		CVTX_EXPORT void cvtx_F3D_M2M_dvort(
			const cvtx_F3D **array_start,
			const int num_filaments,
			const cvtx_P3D **induced_start,
			const int num_induced,
			bsv_V3f *result_array);
		=#
		ccall(
			("cvtx_F3D_M2M_dvort", libcvortex),
			Cvoid,
			(Ptr{Ptr{VortexFilament}}, Cint, Ptr{Ptr{VortexParticle3D}}, Cint,
				Ptr{Vec3f}),
			pargarr, ni, indarg, np, ret)
	end
    return Matrix{Float32}(ret)
end

"""
    induced_velocity_influence_matrix(
        filament_start_coords :: Matrix{<:Real},
        filament_end_coords :: Matrix{<:Real},
        measurement_points :: Matrix{<:Real},
        measurement_directions :: Matrix{<:Real})

The influence of vortex filaments on normal velocities at points in the
domain.

# Arguments
- `filament_start_coord` : The start coordinates of vortex filaments. Matrix
of size N by 3.
- `filament_end_coord` : The end coordinates of vortex filaments. Matrix of
size N by 3.
- `measurement_points` : The positions at which induced velocity is evaluated.
Matrix of size M by 3
- `measurement_directions` : The directions for which velocity compents are
computed. Matrix of size M by 3

Returns a Matrix{Float32} with size (M, N).
"""
function filament_induced_velocity_influence_matrix(
    filament_start_coords :: Matrix{<:Real},
    filament_end_coords :: Matrix{<:Real},
    measurement_points :: Matrix{<:Real},
    measurement_directions :: Matrix{<:Real})

    check_filament_definition(
        filament_start_coords, filament_end_coords)
    @assert(size(measurement_points)[2]==3, "The size of the measurement "*
        "point vector must be N by 3. Actual size is ",
        size(measurement_points), ".")
    @assert(size(measurement_directions)[2]==3, "The size of the measurement "*
        "direction vector must be N by 3. Actual size is ",
        size(measurement_directions), ".")
    @assert(size(measurement_directions)==size(measurement_points),
        "The number of points defined by the measurement points and "*
        "measurement direction vectors should match, but instead measurement"*
        " directions define ", size(measurement_directions)[1], " directions "*
        "and measurement points define ", size(measurement_points), "points.")

    ni = size(filament_end_coords)[1]
    np = size(measurement_directions)[1]
    inducing_filaments = map(
        i->VortexFilament(filament_start_coords[i,:], filament_end_coords[i, :],
        1.0),
        1:ni)

    mes_pnts = map(i->Vec3f(measurement_points[i, :]), 1:np)
    mes_dirs = map(i->Vec3f(measurement_directions[i, :]), 1:np)

    pargarr = Vector{Ptr{VortexFilament}}(undef, length(inducing_filaments))
    for i = 1 : length(pargarr)
        pargarr[i] = Base.pointer(inducing_filaments, i)
    end
    # Julia is column major, C is row major.
    ret = Matrix{Float32}(undef, length(inducing_filaments),
        length(mes_pnts))
	GC.@preserve pargarr inducing_filaments begin
		#=
		CVTX_EXPORT void cvtx_F3D_inf_mtrx(
			const cvtx_F3D **array_start,
			const int num_filaments,
			const bsv_V3f *mes_start,
			const bsv_V3f *dir_start,
			const int num_mes,
			float *result_matrix);
		=#
		ccall(
			("cvtx_F3D_inf_mtrx", libcvortex),
			Cvoid,
			(Ptr{Ptr{VortexFilament}}, Cint, Ptr{Vec3f}, Ptr{Vec3f},
				Cint, Ptr{Float32}),
			pargarr, length(inducing_filaments), mes_pnts, mes_dirs, np, ret)
	end
    return transpose(ret) # Fix row major -> column major
end


text:
struct NaturalExponentialFamily{P,X}
    logh :: Function # :: X -> Real
    t    :: Function # :: X -> StaticVector{N, Real}
    a    :: Function # :: P -> Real
end

struct NatExpFamDist{P,X}
    fam :: NaturalExponentialFamily{P,X}
    # a :: Real # may add this back in later
    η :: P
end

export partype
partype(::Type{NaturalExponentialFamily{P,X}}) where {P,X} = P
partype(::Type{NatExpFamDist{P,X}}) where {P,X} = P

export eltype
eltype(::Type{NaturalExponentialFamily{P,X}}) where {P,X} = X
eltype(::Type{NatExpFamDist{P,X}}) where {P,X} = X

export logpdf
function Distributions.logpdf(d::NatExpFamDist{P,X}, x::X)::Real where {P,X}
    d.fam.logh(x) + d.η ⋅ d.fam.t(x) - d.fam.a(d.η)
end

function Distributions.logpdf(d::NatExpFamDist{P,X}, x)::Real where {P,X}
    d.fam.logh(x) + d.η ⋅ d.fam.t(x) - d.fam.a(d.η)
end

function (fam::NaturalExponentialFamily{P,X})(η::P) where {P,X}
    NatExpFamDist{P,X}(fam, η)
end

export iid
function iid(d::NatExpFamDist{P,X} ) where {P,X}
    logh(xs :: AbstractArray{<:X}) = sum(d.fam.logh.(xs))
    t(xs :: AbstractArray{<:X}) = reduce((a,b) -> a .+ b,d.fam.t.(xs))
    a(η::P) = d.fam.a(η)

    fam = NaturalExponentialFamily{P,AbstractArray{X}}(logh, t, a)
    NatExpFamDist(fam,d.η)
end


text:
using Test
@testset "IBox commands test" begin
data = getData(localIP,8888,"WidthQRS",0,0,10,10)
att = helloIBox.getAttr(localIP,8888,"WidthQRS")
@test att["dstype"]=="interval"
@test att["offsetdata"]=="QPoint"
data = getData(localIP,8888,"QPoint",0,0,1,1)
@test data == [101]

data = getData(localIP,8888,"WidthQRS",0,0,1,30)
@test data[1] == 0x0017

tag = getTag(localIP,8888,"Ecg")
@test tag=="L"

dataCount = helloIBox.getCountDataInInterval(localIP,8888,"QPoint",0,0,0,1000)
@test dataCount==6

amp = getData(localIP,8888,"AmpECG",0,0,1,1)
@test amp[1]==1063
Fr = getData(localIP,8888,"Freq",0,1,1)
@test Fr[1]==257
end

#
# fields = "QPoint,WidthQRS"
# allFields = split(fields,",")
# alldataType = [] #собираем типы данных, чтобы работать с ними
# for f in allFields
#     push!(alldataType, helloIBox.getDataType(localIP,8888,string(f))) #узнаем тип данных
# end