text:
__precompile__()

module Dierckx

# This 'using' is required to suppress a warning about Dierckx
# not having Libdl in its dependencies (Libdl is used by BinaryProvider),
# e.g.: bicycle1885/CodecZlib.jl#26.
using Libdl

export Spline1D,
       Spline2D,
       ParametricSpline,
       evaluate,
       derivative,
       integrate,
       roots,
       evalgrid,
       get_knots,
       get_coeffs,
       get_residual

import Base: show

const depsfile = joinpath(dirname(dirname(@__FILE__)), "deps", "deps.jl")
if isfile(depsfile)
    include(depsfile)
else
    error("Dierckx is not properly installed. Please run Pkg.build(\"Dierckx\")")
end

function __init__()
    check_deps()
end

# ----------------------------------------------------------------------------
# 1-d splines

const _fit1d_messages = Dict(
2=>
"""A theoretically impossible result was found during the iteration
process for finding a smoothing spline with fp = s: s too small.
There is an approximation returned but the corresponding weighted sum
of squared residuals does not satisfy the condition abs(fp-s)/s <
tol.""",
3=>
"""The maximal number of iterations maxit (set to 20 by the program)
allowed for finding a smoothing spline with fp=s has been reached: s
too small. There is an approximation returned but the corresponding
weighted sum of squared residuals does not satisfy the condition
abs(fp-s)/s < tol.""",
10=>
"""Error on entry, no approximation returned. The following conditions
must hold:
1<=k<=5
x[1] < x[2] < ... < x[end]
w[i] > 0.0 for all i

Additionally, if spline knots are given:
length(xknots) <= length(x) + k + 1
x[1] < xknots[1] < xknots[k+2] < ... < xknots[end] < x[end]
The schoenberg-whitney conditions: there must be a subset of data points
xx[j] such that t[j] < xx[j] < t[j+k+1] for j=1,2,...,n-k-1""")

const _eval1d_messages = Dict(
1=>
"""Input point out of range""",
10=>
"""Invalid input data. The following conditions must hold:
length(x) != 0 and xb <= x[1] <= x[2] <= ... x[end] <= xe""")

_translate_bc(bc::AbstractString) = (bc == "extrapolate" ? 0 :
                             bc == "zero" ? 1 :
                             bc == "error" ? 2 :
                             bc == "nearest" ? 3 :
                             error("unknown boundary condition: \"$bc\""))
_translate_bc(bc::Int) = (bc == 0 ? "extrapolate" :
                          bc == 1 ? "zero" :
                          bc == 2 ? "error" :
                          bc == 3 ? "nearest" : "")

mutable struct Spline1D
    t::Vector{Float64}
    c::Vector{Float64}
    k::Int
    bc::Int
    fp::Float64
    wrk::Vector{Float64}
end

# add a constructor that automatically creates the `work` array
Spline1D(t, c, k, bc, fp) = Spline1D(t, c, k, bc, fp, Vector{Float64}(undef, length(t)))

get_knots(spl::Spline1D) = spl.t[spl.k+1:end-spl.k]
get_coeffs(spl::Spline1D) = spl.c[1:end-spl.k+1]
get_residual(spl::Spline1D) = spl.fp

function reallycompact(a::Vector)
    io = IOBuffer()
    io_compact = IOContext(io, :compact => true)
    if length(a) <= 5
        show(io, a)
    else
        write(io, "[")
        show(io_compact, a[1])
        write(io, ",")
        show(io_compact, a[2])
        write(io, " \u2026 ")
        show(io_compact, a[end-1])
        write(io, ",")
        show(io_compact, a[end])
        write(io, "]")
        write(io, " ($(length(a)) elements)")
    end
    seekstart(io)
    return read(io, String)
end

function show(io::IO, spl::Spline1D)
    print(io, """Spline1D(knots=$(reallycompact(get_knots(spl))), k=$(spl.k), extrapolation=\"$(_translate_bc(spl.bc))\", residual=$(spl.fp))""")
end

function Spline1D(x::AbstractVector, y::AbstractVector;
                  w::AbstractVector=ones(length(x)),
                  k::Int=3, s::Real=0.0, bc::AbstractString="nearest",
                  periodic::Bool=false)
    m = length(x)
    length(y) == m || error("length of x and y must match")
    length(w) == m || error("length of x and w must match")
    m > k || error("k must be less than length(x)")
    1 <= k <= 5 || error("1 <= k = $k <= 5 must hold")

    # ensure inputs are of correct type
    xin = convert(Vector{Float64}, x)
    yin = convert(Vector{Float64}, y)
    win = convert(Vector{Float64}, w)

    nest = 0
    if periodic
        nest = max(m + 2k, 2k + 3)
    else
        nest = max(m + k + 1, 2k + 3)
    end

    # outputs
    n = Ref{Int32}(0)
    t = Vector{Float64}(undef, nest)
    c = Vector{Float64}(undef, nest)
    fp = Ref{Float64}(0)
    ier = Ref{Int32}(0)

    # workspace
    lwrk = 0
    if periodic
        lwrk = m * (k + 1) + nest*(8 + 5k)
    else
        lwrk = m * (k + 1) + nest*(7 + 3k)
    end
    wrk = Vector{Float64}(undef, lwrk)
    iwrk = Vector{Int32}(undef, nest)

    if !periodic
        ccall((:curfit_, libddierckx), Nothing,
            (Ref{Int32}, Ref{Int32},  # iopt, m
             Ref{Float64}, Ref{Float64}, Ref{Float64},  # x, y, w
             Ref{Float64}, Ref{Float64},  # xb, xe
             Ref{Int32}, Ref{Float64},  # k, s
             Ref{Int32}, Ref{Int32}, # nest, n
             Ref{Float64}, Ref{Float64}, Ref{Float64},  # t, c, fp
             Ref{Float64}, Ref{Int32}, Ref{Int32},  # wrk, lwrk, iwrk
             Ref{Int32}),  # ier
            0, m, xin, yin, win, xin[1], xin[end], k, Float64(s),
            nest, n, t, c, fp, wrk, lwrk, iwrk, ier)
    else
        ccall((:percur_, libddierckx), Nothing,
            (Ref{Int32}, Ref{Int32}, # iopt, m
             Ref{Float64}, Ref{Float64}, Ref{Float64}, # x, y, w
             Ref{Int32}, Ref{Float64}, # k, s
             Ref{Int32}, Ref{Int32}, # nest, n
             Ref{Float64}, Ref{Float64}, Ref{Float64}, # t, c, fp
             Ref{Float64}, Ref{Int32}, Ref{Int32}, # wrk, lwrk, iwrk
             Ref{Int32}), # ier
            0, m, xin, yin, win, k, Float64(s), nest, n, t, c,
            fp, wrk, lwrk, iwrk, ier)
    end

    ier[] <= 0 || error(_fit1d_messages[ier[]])

    # resize output arrays
    resize!(t, n[])
    resize!(c, n[] - k - 1)

    return Spline1D(t, c, k, _translate_bc(bc), fp[])
end

# version with user-supplied knots
function Spline1D(x::AbstractVector, y::AbstractVector,
                  knots::AbstractVector;
                  w::AbstractVector=ones(length(x)),
                  k::Int=3, bc::AbstractString="nearest",
                  periodic::Bool=false)
    m = length(x)
    length(y) == m || error("length of x and y must match")
    length(w) == m || error("length of x and w must match")
    m > k || error("k must be less than length(x)")
    length(knots) <= m + k + 1 || error("length(knots) <= length(x) + k + 1 must hold")
    first(x) < first(knots) || error("first(x) < first(knots) must hold")
    last(x) > last(knots) || error("last(x) > last(knots) must hold")

    # ensure inputs are of correct type
    xin = convert(Vector{Float64}, x)
    yin = convert(Vector{Float64}, y)
    win = convert(Vector{Float64}, w)

    # x knots
    # (k+1) knots will be added on either end of interior knots.
    n = length(knots) + 2(k + 1)
    t = Vector{Float64}(undef, n)  # All knots
    t[k+2:end-k-1] = knots

    # outputs
    c = Vector{Float64}(undef, n)
    fp = Ref{Float64}(0)
    ier = Ref{Int32}(0)

    # workspace
    lwrk = 0
    if periodic
        lwrk = m * (k + 1) + n*(8 + 5k)
    else
        lwrk = m * (k + 1) + n*(7 + 3k)
    end
    wrk = Vector{Float64}(undef, lwrk)
    iwrk = Vector{Int32}(undef, n)

    if !periodic
        ccall((:curfit_, libddierckx), Nothing,
            (Ref{Int32}, Ref{Int32},  # iopt, m
             Ref{Float64}, Ref{Float64}, Ref{Float64},  # x, y, w
             Ref{Float64}, Ref{Float64},  # xb, xe
             Ref{Int32}, Ref{Float64},  # k, s
             Ref{Int32}, Ref{Int32}, # nest, n
             Ref{Float64}, Ref{Float64}, Ref{Float64},  # t, c, fp
             Ref{Float64}, Ref{Int32}, Ref{Int32},  # wrk, lwrk, iwrk
             Ref{Int32}),  # ier
            -1, m, xin, yin, win, xin[1], xin[end], k, -1.0,
            n, n, t, c, fp, wrk, lwrk, iwrk, ier)
    else
        ccall((:percur_, libddierckx), Nothing,
            (Ref{Int32}, Ref{Int32}, # iopt, m
             Ref{Float64}, Ref{Float64}, Ref{Float64}, # x, y, w
             Ref{Int32}, Ref{Float64}, # k, s
             Ref{Int32}, Ref{Int32}, # nest, n
             Ref{Float64}, Ref{Float64}, Ref{Float64}, # t, c, fp
             Ref{Float64}, Ref{Int32}, Ref{Int32}, # wrk, lwrk, iwrk
             Ref{Int32}), # ier
            -1, m, xin, yin, win, k, -1.0, n, n, t, c,
            fp, wrk, lwrk, iwrk, ier)
    end

    ier[] <= 0 || error(_fit1d_messages[ier[]])
    resize!(c, n - k - 1)

    return Spline1D(t, c, k, _translate_bc(bc), fp[])
end

function _evaluate(t::Vector{Float64}, c::Vector{Float64}, k::Int,
                   x::Vector{Float64}, bc::Int)
    bc in (0, 1, 2, 3) || error("bc = $bc not in (0, 1, 2, 3)")
    m = length(x)
    xin = convert(Vector{Float64}, x)
    y = Vector{Float64}(undef, m)
    ier = Ref{Int32}(0)
    ccall((:splev_, libddierckx), Nothing,
         (Ref{Float64}, Ref{Int32},
          Ref{Float64}, Ref{Int32},
          Ref{Float64}, Ref{Float64}, Ref{Int32},
          Ref{Int32}, Ref{Int32}),
         t, length(t), c, k, xin, y, m, bc, ier)

    ier[] == 0 || error(_eval1d_messages[ier[]])
    return y
end

function _evaluate(t::Vector{Float64}, c::Vector{Float64}, k::Int,
                   x::Real, bc::Int)
    bc in (0, 1, 2, 3) || error("bc = $bc not in (0, 1, 2, 3)")
    y = Ref{Float64}(0)
    ier = Ref{Int32}(0)
    ccall((:splev_, libddierckx), Nothing,
         (Ref{Float64}, Ref{Int32},
          Ref{Float64}, Ref{Int32},
          Ref{Float64}, Ref{Float64}, Ref{Int32},
          Ref{Int32}, Ref{Int32}),
         t, length(t), c, k, Float64(x), y, 1, bc, ier)

    ier[] == 0 || error(_eval1d_messages[ier[]])
    return y[]
end

evaluate(spline::Spline1D, x::AbstractVector) =
    _evaluate(spline.t, spline.c, spline.k,
              convert(Vector{Float64}, x), spline.bc)

evaluate(spline::Spline1D, x::Real) =
    _evaluate(spline.t, spline.c, spline.k, x, spline.bc)

function _derivative(t::Vector{Float64}, c::Vector{Float64}, k::Int,
                     x::Vector{Float64}, nu::Int, bc::Int, wrk::Vector{Float64})
    (1 <= nu <= k) || error("order of derivative must be positive and less than or equal to spline order")
    m = length(x)
    n = length(t)
    y = Vector{Float64}(undef, m)
    ier = Ref{Int32}(0)
    ccall((:splder_, libddierckx), Nothing,
         (Ref{Float64}, Ref{Int32}, # t, n
          Ref{Float64}, Ref{Int32}, # c, k
          Ref{Int32}, # nu
          Ref{Float64}, Ref{Float64}, Ref{Int32}, # x, y, m
          Ref{Int32}, Ref{Float64}, Ref{Int32}), # e, wrk, ier
         t, n, c, k, nu, x, y, m, bc, wrk, ier)
    ier[] == 0 || error(_eval1d_messages[ier[]])
    return y
end

function _derivative(t::Vector{Float64}, c::Vector{Float64}, k::Int,
                     x::Real, nu::Int, bc::Int, wrk::Vector{Float64})
    (1 <= nu <= k) || error("order of derivative must be positive and less than or equal to spline order")
    n = length(t)
    y = Ref{Float64}(0)
    ier = Ref{Int32}(0)
    ccall((:splder_, libddierckx), Nothing,
         (Ref{Float64}, Ref{Int32}, # t, n
          Ref{Float64}, Ref{Int32}, # c, k
          Ref{Int32}, # nu
          Ref{Float64}, Ref{Float64}, Ref{Int32}, # x, y, m
          Ref{Int32}, Ref{Float64}, Ref{Int32}), # e, wrk, ier
         t, n, c, k, nu, Float64(x), y, 1, bc, wrk, ier)
    ier[] == 0 || error(_eval1d_messages[ier[]])
    return y[]
end

# TODO: should the function name be evaluate, derivative, or grad?
#       or should it be integrated with evaluate, above?
#       (problem with that: derivative doesn't accept bc="nearest")
# TODO: should `nu` be `d`?
derivative(spline::Spline1D, x::AbstractVector, nu::Int=1) =
    _derivative(spline.t, spline.c, spline.k,
                convert(Vector{Float64}, x), nu, spline.bc, spline.wrk)

derivative(spline::Spline1D, x::Real, nu::Int=1) =
    _derivative(spline.t, spline.c, spline.k, x, nu, spline.bc, spline.wrk)

# TODO: deprecate this?
derivative(spline::Spline1D, x; nu::Int=1) = derivative(spline, x, nu)

function _integrate(t::Vector{Float64}, c::Vector{Float64}, k::Int,
                    a::Real, b::Real, wrk::Vector{Float64})
    n = length(t)
    ccall((:splint_, libddierckx), Float64,
        (Ref{Float64}, Ref{Int32},
         Ref{Float64}, Ref{Int32},
         Ref{Float64}, Ref{Float64},
         Ref{Float64}),
         t, n, c, k,
         Float64(a), Float64(b), wrk)
end

integrate(spline::Spline1D, a::Real, b::Real) =
    _integrate(spline.t, spline.c, spline.k, a, b, spline.wrk)

# TODO roots for parametric splines
# note: default maxn in scipy.interpolate is 3 * (length(spline.t) - 7)
function roots(spline::Spline1D; maxn::Integer=8)
    if spline.k != 3
        error("root finding only supported for cubic splines (k=3)")
    end
    n = length(spline.t)
    zeros = Vector{Float64}(undef, maxn)
    m = Vector{Int32}(undef, 1)
    ier = Vector{Int32}(undef, 1)
    ccall((:sproot_, libddierckx), Nothing,
          (Ref{Float64}, Ref{Int32},  # t, n
           Ref{Float64}, Ref{Float64},  # c, zeros
           Ref{Int32}, Ref{Int32},  # mest, m
           Ref{Int32}),  # ier
          spline.t, n, spline.c, zeros,
          maxn, m, ier)

    if ier[1] == 0
        return zeros[1:m[1]]
    elseif ier[1] == 1
        warn("number of zeros exceeded maxn; only first maxn zeros returned")
        return zeros
    elseif ier[1] == 10
        error("Invalid input data.")
    else
        error("unknown error code in sproot: $(ier[1])")
    end
end

# ----------------------------------------------------------------------------
# parametric splines

mutable struct ParametricSpline
    t::Vector{Float64}
    c::Matrix{Float64}
    k::Int
    bc::Int
    fp::Float64
    wrk::Vector{Float64}
end

ParametricSpline(t, c, k, bc, fp) =
    ParametricSpline(t, c, k, bc, fp, Vector{Float64}(undef, length(t)))

get_knots(spl::ParametricSpline) = spl.t[spl.k+1:end-spl.k]
get_coeffs(spl::ParametricSpline) = spl.c[:, 1:end-spl.k+1]
get_residual(spl::ParametricSpline) = spl.fp

function show(io::IO, spl::ParametricSpline)
    print(io, """ParametricSpline(knots=$(reallycompact(get_knots(spl))), k=$(spl.k), extrapolation=\"$(_translate_bc(spl.bc))\", residual=$(spl.fp))""")
end

function ParametricSpline(x::AbstractMatrix;
                          w::AbstractVector=ones(size(x, 2)),
                          k::Int=3, s::Real=0.0, bc::AbstractString="nearest",
                          periodic::Bool=false)
    return _ParametricSpline(nothing, x, nothing, w, k, s, bc, periodic)
end

# version with user-supplied u
function ParametricSpline(u::AbstractVector, x::AbstractMatrix;
                          ub::Real=u[1], ue::Real=u[end],
                          w::AbstractVector=ones(size(x, 2)),
                          k::Int=3, s::Real=0.0, bc::AbstractString="nearest",
                          periodic::Bool=false)
    return _ParametricSpline(u, x, nothing, w, k, s, bc, periodic)
end

# version with user-supplied knots
function ParametricSpline(x::AbstractMatrix, knots::AbstractVector;
                          w::AbstractVector=ones(size(x, 2)),
                          k::Int=3, bc::AbstractString="nearest",
                          periodic::Bool=false)
    return _ParametricSpline(nothing, x, knots, w, k, -1.0, bc, periodic)
end

# version with user-supplied u and knots
function ParametricSpline(u::AbstractVector, x::AbstractMatrix,
                          knots::AbstractVector;
                          w::AbstractVector=ones(size(x, 2)),
                          k::Int=3, bc::AbstractString="nearest",
                          periodic::Bool=false)
    return _ParametricSpline(u, x, knots, w, k, -1.0, bc, periodic)
end

function _ParametricSpline(u::Union{AbstractVector, Nothing}, x::AbstractMatrix,
                           knots::Union{AbstractVector, Nothing},
                           w::AbstractVector, k::Int, s::Real,
                           bc::AbstractString, periodic::Bool)
    idim, m = size(x)
    if periodic
        x[:, 1] == x[:, end] || error("for periodic splines x[:,1] and x[:,end] must match")
    end

    length(w) == m || error("number of data points and length of w must match")
    0 < idim < 11 || error("number of dimension must be between 1 and 10")
    m > k || error("number of data points muct be greater than k")
    1 <= k <= 5 || error("1 <= k = $k <= 5 must hold")

    local ipar, uin, ub, ue
    if u != nothing
        all(u[1:end-1] .< u[2:end]) || error("u[i] must be strictly increasing")
        ipar = 1
        uin = convert(Vector{Float64}, u)
        ub = uin[1]
        ue = uin[end]
    else
        ipar = 0
        uin = Vector{Float64}(undef, size(x, 2))
        ub = 0.0
        ue = 1.0
    end

    local iopt, nest::Int, t
    if knots != nothing
        length(knots) <= m + k + 1 || error("length(knots) <= length(x) + k + 1 must hold")
        first(x) < first(knots) || error("first(x) < first(knots) must hold")
        last(x) > last(knots) || error("last(x) > last(knots) must hold")
        iopt = -1
        nest = length(knots) + 2(k + 1)
        t = Vector{Float64}(undef, nest)
        t[k+2:end-k-1] = knots
    else
        iopt = 0
        nest = m + 2k
        if s == 0
            nest = periodic ? m + 2k : m + k + 1
        end
        nest = max(nest, 2k + 3)
        t = Vector{Float64}(undef, nest)
    end

    xin = convert(Matrix{Float64}, x)
    win = convert(Vector{Float64}, w)
    n = Ref{Int32}(nest)
    c = Vector{Float64}(undef, idim*nest)
    fp = Ref{Float64}(0)
    ier = Ref{Int32}(0)

    local lwrk::Int
    if periodic
        lwrk = m*(k + 1) + nest*(7 + idim + 5k)
    else
        lwrk = m*(k + 1) + nest*(6 + idim + 3k)
    end
    wrk = Vector{Float64}(undef, lwrk)
    iwrk = Vector{Int32}(undef, nest)

    if !periodic
        ccall((:parcur_, libddierckx), Nothing,
              (Ref{Int32}, Ref{Int32}, # iopt, ipar
              Ref{Int32}, Ref{Int32}, # idim, m
              Ref{Float64}, Ref{Int32}, # u, mx
              Ref{Float64}, Ref{Float64}, # x, w
              Ref{Float64}, Ref{Float64}, # ub, ue
              Ref{Int32}, Ref{Float64}, # k, s
              Ref{Int32}, Ref{Int32}, # nest, n
              Ref{Float64}, Ref{Int32}, # t, nc
              Ref{Float64}, Ref{Float64}, # c, fp
              Ref{Float64}, Ref{Int32}, # wrk, lwrk
              Ref{Int32}, Ref{Int32}), #iwrk, ier
              iopt, ipar,
              idim, m,
              uin, length(x),
              xin, win,
              ub, ue,
              k, s,
              nest, n,
              t, length(c),
              c, fp,
              wrk, lwrk,
              iwrk, ier)
    else
        ccall((:clocur_, libddierckx), Nothing,
              (Ref{Int32}, Ref{Int32}, # iopt, ipar
              Ref{Int32}, Ref{Int32}, # idim, m
              Ref{Float64}, Ref{Int32}, # u, mx
              Ref{Float64}, Ref{Float64}, # x, w
              Ref{Int32}, Ref{Float64}, # k, s
              Ref{Int32}, Ref{Int32}, # nest, n
              Ref{Float64}, Ref{Int32}, # t, nc
              Ref{Float64}, Ref{Float64}, # c, fp
              Ref{Float64}, Ref{Int32}, # wrk, lwrk
              Ref{Int32}, Ref{Int32}), #iwrk, ier
              iopt, ipar,
              idim, m,
              uin, length(x),
              xin, win,
              k, s,
              nest, n,
              t, length(c),
              c, fp,
              wrk, lwrk,
              iwrk, ier)
    end

    ier[] <= 0 || error(_fit1d_messages[ier[]])

    resize!(t, n[])
    c = [c[n[]*(j-1) + i] for j=1:idim, i=1:n[]-k-1]

    return ParametricSpline(t, c, k, _translate_bc(bc), fp[])
end

_evaluate(t::Vector{Float64}, c::Matrix{Float64}, k::Int,
          x::Vector{Float64}, bc::Int) =
    mapslices(v -> _evaluate(t, v, k, x, bc), c, dims=[2])

_evaluate(t::Vector{Float64}, c::Matrix{Float64}, k::Int,
          x::Real, bc::Int) =
    vec(mapslices(v -> _evaluate(t, v, k, x, bc), c, dims=[2]))

function evaluate(spline::ParametricSpline, x::AbstractVector)
    xin = convert(Vector{Float64}, x)
    _evaluate(spline.t, spline.c, spline.k, xin, spline.bc)
end

evaluate(spline::ParametricSpline, x::Real) =
    _evaluate(spline.t, spline.c, spline.k, x, spline.bc)

_derivative(t::Vector{Float64}, c::Matrix{Float64}, k::Int,
            x::Vector{Float64}, nu::Int, bc::Int, wrk::Vector{Float64}) =
    mapslices(v -> _derivative(t, v, k, x, nu, bc, wrk), c, dims=[2])

_derivative(t::Vector{Float64}, c::Matrix{Float64}, k::Int,
            x::Real, nu::Int, bc::Int, wrk::Vector{Float64}) =
    vec(mapslices(v -> _derivative(t, v, k, x, nu, bc, wrk), c, dims=[2]))

derivative(spline::ParametricSpline, x::AbstractVector, nu::Int=1) =
    _derivative(spline.t, spline.c, spline.k,
                convert(Vector{Float64}, x), nu, spline.bc, spline.wrk)

# TODO: deprecate this
derivative(spline::ParametricSpline, x; nu::Int=1) = derivative(spline, x, nu)

derivative(spline::ParametricSpline, x::Real, nu::Int=1) =
    _derivative(spline.t, spline.c, spline.k, x, nu, spline.bc, spline.wrk)

_integrate(t::Vector{Float64}, c::Matrix{Float64}, k::Int,
           a::Real, b::Real, wrk::Vector{Float64}) =
    vec(mapslices(v -> _integrate(t, v, k, a, b, wrk), c, dims=[2]))

integrate(spline::ParametricSpline, a::Real, b::Real) =
    _integrate(spline.t, spline.c, spline.k, a, b, spline.wrk)

# ----------------------------------------------------------------------------
# 2-d splines

# NOTE REGARDING ARGUMENT ORDER: In the "grid" version of the Spline2D
# constructor and evaluators, the fortran functions expects z to have
# shape (my, mx), but we'd rather have x be the fast axis in z.  So,
# in the ccall()s in these methods, all the x and y related inputs are
# swapped with regard to what the Fortran documentation says.

const _fit2d_messages = Dict(
-3=>

"""The coefficients of the spline returned have been computed as the
minimal norm least-squares solution of a (numerically) rank deficient
system (deficiency=%i). If deficiency is large, the results may be
inaccurate. Deficiency may strongly depend on the value of eps.""",

1=>

"""The required storage space exceeds the available storage space:
nxest or nyest too small, or s too small. Try increasing s.""",
# The weighted least-squares spline corresponds to the current set of knots.

2=>

"""A theoretically impossible result was found during the iteration
process for finding a smoothing spline with fp = s: s too small or
badly chosen eps.  Weighted sum of squared residuals does not satisfy
abs(fp-s)/s < tol.""",

3=>

"""the maximal number of iterations maxit (set to 20 by the program)
allowed for finding a smoothing spline with fp=s has been reached: s
too small.  Weighted sum of squared residuals does not satisfy
abs(fp-s)/s < tol.""",

4=>

"""No more knots can be added because the number of b-spline
coefficients (nx-kx-1)*(ny-ky-1) already exceeds the number of data
points m: either s or m too small.  The weighted least-squares spline
corresponds to the current set of knots.""",

5=>

"""No more knots can be added because the additional knot would
(quasi) coincide with an old one: s too small or too large a weight to
an inaccurate data point.  The weighted least-squares spline
corresponds to the current set of knots.""",

10=>

"""Error on entry, no approximation returned. The following conditions
must hold:
xb<=x[i]<=xe, yb<=y[i]<=ye, w[i]>0, i=0..m-1
If iopt==-1, then
  xb<tx[kx+1]<tx[kx+2]<...<tx[nx-kx-2]<xe
  yb<ty[ky+1]<ty[ky+2]<...<ty[ny-ky-2]<ye""")

const _eval2d_message = (
"""Invalid input data. Restrictions:
length(x) != 0, length(y) != 0
x[i-1] <= x[i] for i=2,...,length(x)
y[j-1] <= y[j] for j=2,...,length(y)
""")

mutable struct Spline2D
    tx::Vector{Float64}
    ty::Vector{Float64}
    c::Vector{Float64}
    kx::Int
    ky::Int
    fp::Float64
end

get_knots(spl::Spline2D) = (spl.tx[spl.kx+1:end-spl.kx],
                            spl.ty[spl.ky+1:end-spl.ky])
get_residual(spl::Spline2D) = spl.fp

# Helper functions for calculating required size of work arrays in surfit.
# Note that x and y here are as in the Fortran documentation.
# These are translated from scipy/interpolate/src/fitpack.pyf
function calc_surfit_lwrk1(m, kx, ky, nxest, nyest)
    u = nxest - kx - 1
    v = nyest - ky - 1
    km = max(kx, ky) + 1
    ne = max(nxest, nyest)
    bx = kx*v + ky + 1
    by = ky*u + kx + 1
    b1 = b2 = 0
    if (bx<=by)
        b1 = bx
        b2 = bx + v - ky
    else
        b1 = by
        b2 = by + u - kx
    end
    return u*v*(2 + b1 + b2) + 2*(u+v+km*(m+ne)+ne-kx-ky) + b2 + 1
end

function calc_surfit_lwrk2(m, kx, ky, nxest, nyest)
    u = nxest - kx - 1
    v = nyest - ky - 1
    bx = kx * v + ky + 1
    by = ky * u + kx + 1
    b2 = (bx <= by ? bx + v - ky : by + u - kx)
    return u * v * (b2 + 1) + b2
end

# Construct spline from unstructured data
function Spline2D(x::AbstractVector, y::AbstractVector, z::AbstractVector;
                  w::AbstractVector=ones(length(x)), kx::Int=3, ky::Int=3,
                  s::Real=0.0)

    # array sizes
    m = length(x)
    (length(y) == length(z) == m) || error("lengths of x, y, z must match")
    (length(w) == m) || error("length of w must match other inputs")

    nxest = max(kx+1+ceil(Int,sqrt(m/2)), 2*(kx+1))
    nyest = max(ky+1+ceil(Int,sqrt(m/2)), 2*(ky+1))
    nmax = max(nxest, nyest)

    eps = 1.0e-16

    # bounds
    xb = minimum(x)
    xe = maximum(x)
    yb = minimum(y)
    ye = maximum(y)

    # ensure arrays are of correct type
    xin = convert(Vector{Float64}, x)
    yin = convert(Vector{Float64}, y)
    zin = convert(Vector{Float64}, z)
    win = convert(Vector{Float64}, w)

    # return values
    nx = Ref{Int32}()
    tx = Vector{Float64}(undef, nxest)
    ny = Ref{Int32}()
    ty = Vector{Float64}(undef, nyest)
    c = Vector{Float64}(undef, (nxest-kx-1) * (nyest-ky-1))
    fp = Ref{Float64}()
    ier = Ref{Int32}()

    # work arrays
    # Note: in lwrk1 and lwrk2, x and y are swapped on purpose.
    lwrk1 = calc_surfit_lwrk1(m, ky, kx, nyest, nxest)
    lwrk2 = calc_surfit_lwrk2(m, ky, kx, nyest, nxest)
    kwrk = m + (nxest - 2*kx - 1) * (nyest - 2*ky - 1)
    wrk1 = Vector{Float64}(undef, lwrk1)
    wrk2 = Vector{Float64}(undef, lwrk2)
    iwrk = Vector{Int32}(undef, kwrk)

    ccall((:surfit_, libddierckx), Nothing,
          (Ref{Int32}, Ref{Int32},  # iopt, m
           Ref{Float64}, Ref{Float64},  # y, x
           Ref{Float64}, Ref{Float64},  # z, w
           Ref{Float64}, Ref{Float64},  # yb, ye
           Ref{Float64}, Ref{Float64},  # xb, xe
           Ref{Int32}, Ref{Int32}, Ref{Float64},  # ky, kx, s
           Ref{Int32}, Ref{Int32}, Ref{Int32},  # nyest, nxest, nmax
           Ref{Float64},  # eps
           Ref{Int32}, Ref{Float64},  # ny, tx
           Ref{Int32}, Ref{Float64},  # nx, tx
           Ref{Float64}, Ref{Float64},  # c, fp
           Ref{Float64}, Ref{Int32},  # wrk1, lwrk1
           Ref{Float64}, Ref{Int32},  # wrk2, lwrk2
           Ref{Int32}, Ref{Int32}, Ref{Int32}),   # iwrk, kwrk, ier
          0, m, yin, xin, zin, win, yb, ye, xb, xe, ky, kx,
          Float64(s), nyest, nxest, nmax, eps, ny, ty, nx, tx,
          c, fp, wrk1, lwrk1, wrk2, lwrk2, iwrk, kwrk, ier)

    while ier[] > 10
        # lwrk2 is too small, i.e., there is not enough workspace
        # for computing the minimal least-squares solution of a rank
        # deficient system of linear equations. ier gives the
        # requested value for lwrk2. Rerun with that value in "continue"
        # mode, with iopt = 1.
        lwrk2 = ier[]
        resize!(wrk2, lwrk2)
        ccall((:surfit_, libddierckx), Nothing,
              (Ref{Int32}, Ref{Int32},  # iopt, m
               Ref{Float64}, Ref{Float64},  # y, x
               Ref{Float64}, Ref{Float64},  # z, w
               Ref{Float64}, Ref{Float64},  # yb, ye
               Ref{Float64}, Ref{Float64},  # xb, xe
               Ref{Int32}, Ref{Int32}, Ref{Float64},  # ky, kx, s
               Ref{Int32}, Ref{Int32}, Ref{Int32},  # nyest, nxest, nmax
               Ref{Float64},  # eps
               Ref{Int32}, Ref{Float64},  # ny, tx
               Ref{Int32}, Ref{Float64},  # nx, tx
               Ref{Float64}, Ref{Float64},  # c, fp
               Ref{Float64}, Ref{Int32},  # wrk1, lwrk1
               Ref{Float64}, Ref{Int32},  # wrk2, lwrk2
               Ref{Int32}, Ref{Int32}, Ref{Int32}),   # iwrk, kwrk, ier
              1, m, yin, xin, zin, win, yb, ye, xb, xe, ky, kx, s,
              nyest, nxest, nmax, eps, ny, ty, nx, tx, c, fp,
              wrk1, lwrk1, wrk2, lwrk2, iwrk, kwrk, ier)
    end

    if (ier[] == 0 || ier[] == -1 || ier[] == -2)
        # good values, pass.
    elseif ier[] < -2
        warn("""
        The coefficients of the spline returned have been
        computed as the minimal norm least-squares solution of a
        (numerically) rank deficient system. The rank is $(-ier[]).
        The rank deficiency is $((nx[]-kx-1)*(ny[]-ky-1)+ier[]).
        Especially if the rank deficiency is large the results may
        be inaccurate.""")
        # "The results could also seriously depend on the value of
        # eps" (not in message because eps is currently not an input)
    else
        error(_fit2d_messages[ier[]])
    end

    # Resize output arrays to the size actually used.
    resize!(tx, nx[])
    resize!(ty, ny[])
    resize!(c, (nx[] - kx - 1) * (ny[] - ky - 1))

    return Spline2D(tx, ty, c, kx, ky, fp[])
end

# Construct spline from data on a grid.
function Spline2D(x::AbstractVector, y::AbstractVector, z::AbstractMatrix;
                  kx::Int=3, ky::Int=3, s::Real=0.0)
    mx = length(x)
    my = length(y)
    @assert size(z, 1) == mx && size(z, 2) == my

    mx > kx || error("length(x) must be greater than kx")
    my > ky || error("length(y) must be greater than ky")

    # Bounds
    xb = x[1]
    xe = x[end]
    yb = y[1]
    ye = y[end]
    nxest = mx+kx+1
    nyest = my+ky+1

    # ensure arrays are of correct type
    xin = convert(Vector{Float64}, x)
    yin = convert(Vector{Float64}, y)
    zin = convert(Matrix{Float64}, z)

    # Return values
    nx = Ref{Int32}()
    tx = Vector{Float64}(undef, nxest)
    ny = Ref{Int32}()
    ty = Vector{Float64}(undef, nyest)
    c = Vector{Float64}(undef, (nxest-kx-1) * (nyest-ky-1))
    fp = Ref{Float64}()
    ier = Ref{Int32}()

    # Work arrays.
    # Note that in lwrk, x and y are swapped with respect to the Fortran
    # documentation. See "NOTE REGARDING ARGUMENT ORDER" above.
    lwrk = (4 + nyest * (mx+2*ky+5) + nxest * (2*kx+5) +
            my*(ky+1) + mx*(kx+1) + max(mx, nyest))
    wrk = Vector{Float64}(undef, lwrk)
    kwrk = 3 + mx + my + nxest + nyest
    iwrk = Vector{Int32}(undef, kwrk)

    ccall((:regrid_, libddierckx), Nothing,
          (Ref{Int32},  # iopt
           Ref{Int32}, Ref{Float64},  # my, y
           Ref{Int32}, Ref{Float64},  # mx, x
           Ref{Float64},  # z
           Ref{Float64}, Ref{Float64},  # yb, ye
           Ref{Float64}, Ref{Float64},  # xb, xe
           Ref{Int32}, Ref{Int32}, Ref{Float64},  # ky, kx, s
           Ref{Int32}, Ref{Int32},  # nyest, nxest
           Ref{Int32}, Ref{Float64},  # ny, ty
           Ref{Int32}, Ref{Float64},  # nx, tx
           Ref{Float64}, Ref{Float64},  # c, fp
           Ref{Float64}, Ref{Int32},  # wrk, lwrk
           Ref{Int32}, Ref{Int32},  # iwrk, lwrk
           Ref{Int32}),  # ier
          0f0, my, yin, mx, xin, zin, yb, ye, xb, xe, ky, kx,
          Float64(s), nyest, nxest, ny, ty, nx, tx, c, fp,
          wrk, lwrk, iwrk, kwrk, ier)

    if !(ier[] == 0 || ier[] == -1 || ier[] == -2)
        error(_fit2d_messages[ier[]])
    end

    # Resize output arrays to the size actually used.
    resize!(tx, nx[])
    resize!(ty, ny[])
    resize!(c, (nx[] - kx - 1) * (ny[] - ky - 1))

    return Spline2D(tx, ty, c, kx, ky, fp[])
end

# Evaluate spline at individual points
function evaluate(spline::Spline2D, x::AbstractVector, y::AbstractVector)
    m = length(x)
    @assert length(y) == m

    xin = convert(Vector{Float64}, x)
    yin = convert(Vector{Float64}, y)

    ier = Ref{Int32}()
    lwrk = spline.kx + spline.ky + 2
    wrk = Vector{Float64}(undef, lwrk)
    z = Vector{Float64}(undef, m)

    ccall((:bispeu_, libddierckx), Nothing,
          (Ref{Float64}, Ref{Int32},  # ty, ny
           Ref{Float64}, Ref{Int32},  # tx, nx
           Ref{Float64},  # c
           Ref{Int32}, Ref{Int32},  # ky, kx
           Ref{Float64}, Ref{Float64}, Ref{Float64},  # y, x, z
           Ref{Int32},  # m
           Ref{Float64}, Ref{Int32}, Ref{Int32}),  # wrk, lwrk, ier
          spline.ty, length(spline.ty),
          spline.tx, length(spline.tx),
          spline.c, spline.ky, spline.kx, yin, xin, z, m,
          wrk, lwrk, ier)

    ier[] == 0 || error(_eval2d_message)

    return z
end

function evaluate(spline::Spline2D, x::Real, y::Real)
    ier = Ref{Int32}()
    lwrk = spline.kx + spline.ky + 2
    wrk = Vector{Float64}(undef, lwrk)
    z = Ref{Float64}()
    ccall((:bispeu_, libddierckx), Nothing,
          (Ref{Float64}, Ref{Int32},  # ty, ny
           Ref{Float64}, Ref{Int32},  # tx, nx
           Ref{Float64},  # c
           Ref{Int32}, Ref{Int32},  # ky, kx
           Ref{Float64}, Ref{Float64}, Ref{Float64},  # y, x, z
           Ref{Int32},  # m
           Ref{Float64}, Ref{Int32}, Ref{Int32}),  # wrk, lwrk, ier
          spline.ty, length(spline.ty),
          spline.tx, length(spline.tx),
          spline.c, spline.ky, spline.kx, y, x, z, 1,
          wrk, lwrk, ier)

    ier[] == 0 || error(_eval2d_message)
    return z[]
end

# Evaluate spline on the grid spanned by the input arrays.
function evalgrid(spline::Spline2D, x::AbstractVector, y::AbstractVector)
    mx = length(x)
    my = length(y)

    xin = convert(Vector{Float64}, x)
    yin = convert(Vector{Float64}, y)

    lwrk = mx*(spline.kx + 1) + my*(spline.ky + 1)
    wrk = Vector{Float64}(undef, lwrk)
    kwrk = mx + my
    iwrk = Vector{Int32}(undef, kwrk)
    ier = Ref{Int32}()
    z = Matrix{Float64}(undef, mx, my)

    ccall((:bispev_, libddierckx), Nothing,
          (Ref{Float64}, Ref{Int32},  # ty, ny
           Ref{Float64}, Ref{Int32},  # tx, nx
           Ref{Float64},  # c
           Ref{Int32}, Ref{Int32},  # ky, kx
           Ref{Float64}, Ref{Int32},  # y, my
           Ref{Float64}, Ref{Int32},  # x, mx
           Ref{Float64},  # z
           Ref{Float64}, Ref{Int32},  # wrk, lwrk
           Ref{Int32}, Ref{Int32},  # iwrk, kwrk
           Ref{Int32}),  # ier
          spline.ty, length(spline.ty),
          spline.tx, length(spline.tx),
          spline.c, spline.ky, spline.kx, yin, my, xin, mx, z,
          wrk, lwrk, iwrk, kwrk, ier)

    ier[] == 0 || error(_eval2d_message)

    return z
end

# 2D integration
function integrate(spline::Spline2D, xb::Real, xe::Real, yb::Real, ye::Real)
        nx = length(spline.tx)
        ny = length(spline.ty)

        kx = spline.kx
        ky = spline.ky

        wrk = Vector{Float64}(undef, nx + ny -kx - ky -2)
        ccall((:dblint_, libddierckx), Float64,
              (Ref{Float64}, Ref{Int32},  # tx, nx
               Ref{Float64}, Ref{Int32},  # ty, ny
               Ref{Float64}, Ref{Int32},  # c, kx
               Ref{Int32}, #ky
               Ref{Float64}, Ref{Float64},  # xb, xe
               Ref{Float64}, Ref{Float64},  # yb, ye
               Ref{Float64}),  # wrk
              spline.tx, nx,
              spline.ty, ny,
              spline.c, spline.kx,
              spline.ky,
              xb, xe,
              yb, ye,
              wrk)
end

# call synonyms for evaluate():
(spl::Spline1D)(x::Real) = evaluate(spl, x)
(spl::Spline1D)(x::AbstractVector) = evaluate(spl, x)
(spl::ParametricSpline)(x::Real) = evaluate(spl, x)
(spl::ParametricSpline)(x::AbstractVector) = evaluate(spl, x)
(spl::Spline2D)(x::Real, y::Real) = evaluate(spl, x, y)
(spl::Spline2D)(x::AbstractVector, y::AbstractVector) =
    evaluate(spl, x, y)

end # module


text:
#####
##### Dummy runs are used to ensure the absence of runtime errors in the code
##### before launching a training session
#####

# Returned modified parameters where all num_games fields are set to 1.
# The number of iterations is set to 2.
function dummy_run_params(params, benchmark)
  self_play = SelfPlayParams(params.self_play, num_games=1)
  arena = ArenaParams(params.arena, num_games=1)
  learning = LearningParams(params.learning,
    max_batches_per_checkpoint=2,
    num_checkpoints=min(params.learning.num_checkpoints, 2))
  params = Params(params, num_iters=2,
    self_play=self_play, arena=arena, learning=learning)
  benchmark = [
    Benchmark.Duel(d.player, d.baseline,
      num_games=1, color_policy=d.color_policy)
    for d in benchmark ]
  return params, benchmark
end

function dummy_run(GameModule)
  Game = GameModule.Game
  Training = GameModule.Training
  Net = Training.Network
  netparams = Training.netparams
  params, benchmark = dummy_run_params(Training.params, Training.benchmark)
  session = Session(Game, Net{Game}, params, netparams,
    benchmark=benchmark, nostdout=true)
  resume!(session)
  return true
end

ext: jl


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

@doc raw"""DeploymentStrategy describes how to replace existing pods with new ones.

    IoK8sApiAppsV1beta2DeploymentStrategy(;
        rollingUpdate=nothing,
        type=nothing,
    )

    - rollingUpdate::IoK8sApiAppsV1beta2RollingUpdateDeployment : Rolling update config params. Present only if DeploymentStrategyType = RollingUpdate.
    - type::String : Type of deployment. Can be \"Recreate\" or \"RollingUpdate\". Default is RollingUpdate.
"""
mutable struct IoK8sApiAppsV1beta2DeploymentStrategy <: SwaggerModel
    rollingUpdate::Any # spec type: Union{ Nothing, IoK8sApiAppsV1beta2RollingUpdateDeployment } # spec name: rollingUpdate
    type::Any # spec type: Union{ Nothing, String } # spec name: type

    function IoK8sApiAppsV1beta2DeploymentStrategy(;rollingUpdate=nothing, type=nothing)
        o = new()
        validate_property(IoK8sApiAppsV1beta2DeploymentStrategy, Symbol("rollingUpdate"), rollingUpdate)
        setfield!(o, Symbol("rollingUpdate"), rollingUpdate)
        validate_property(IoK8sApiAppsV1beta2DeploymentStrategy, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        o
    end
end # type IoK8sApiAppsV1beta2DeploymentStrategy

const _property_map_IoK8sApiAppsV1beta2DeploymentStrategy = Dict{Symbol,Symbol}(Symbol("rollingUpdate")=>Symbol("rollingUpdate"), Symbol("type")=>Symbol("type"))
const _property_types_IoK8sApiAppsV1beta2DeploymentStrategy = Dict{Symbol,String}(Symbol("rollingUpdate")=>"IoK8sApiAppsV1beta2RollingUpdateDeployment", Symbol("type")=>"String")
Base.propertynames(::Type{ IoK8sApiAppsV1beta2DeploymentStrategy }) = collect(keys(_property_map_IoK8sApiAppsV1beta2DeploymentStrategy))
Swagger.property_type(::Type{ IoK8sApiAppsV1beta2DeploymentStrategy }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAppsV1beta2DeploymentStrategy[name]))}
Swagger.field_name(::Type{ IoK8sApiAppsV1beta2DeploymentStrategy }, property_name::Symbol) =  _property_map_IoK8sApiAppsV1beta2DeploymentStrategy[property_name]

function check_required(o::IoK8sApiAppsV1beta2DeploymentStrategy)
    true
end

function validate_property(::Type{ IoK8sApiAppsV1beta2DeploymentStrategy }, name::Symbol, val)
end


text:
using DataDeps

register(DataDep(
  "statlog-project-german-credit",
  "http://archive.ics.uci.edu/ml/datasets/Statlog+Project",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/statlog/german/german.data",
  "b21f3d81db8071257d5ff1deaeba1fd4303b62712e6fcc9715c7a86202cb5871",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=21,
        feature_indices=1:20,
        categoric_indices=[1; 3:4; 6:7; 9:10; 12; 14:15; 17; 19:20],
        delim=' ',
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
function create_bibfile(savelocation)
  bibfile = joinpath(savelocation, "internalbibtex.bib")
  if isfile(bibfile)
    return bibfile
  end
  open(bibfile, "w")
  return bibfile
end

function load_database(savelocation)
  maindb = sqlite.DB(joinpath(savelocation, "notes2graphdb.sqlite"))
  return maindb
end

function initialize_database(savelocation)
  maindb = sqlite.DB(joinpath(savelocation, "notes2graphdb.sqlite"))
  # t1 is a table that keeps a list of all nodes, i.e. concepts
  stmt = "CREATE TABLE IF NOT EXISTS t1 (nodeid INTEGER NOT NULL UNIQUE, nodename TEXT NOT NULL UNIQUE COLLATE NOCASE, PRIMARY KEY (nodeid, nodename))"
  sqlite.query(maindb, stmt)

  # t2 is a table that keeps a list of descriptions and references for each node or relation
  stmt = "CREATE TABLE IF NOT EXISTS t2 (descrid INTEGER NOT NULL UNIQUE PRIMARY KEY, descr TEXT NOT NULL UNIQUE COLLATE NOCASE, ref TEXT)"
  sqlite.query(maindb, stmt)

  # t3 is a table that keeps a list of associations between node ids and description ids
  stmt = "CREATE TABLE IF NOT EXISTS t3 (nodeid INTEGER NOT NULL, descrid INTEGER NOT NULL, PRIMARY KEY (nodeid, descrid))"
  sqlite.query(maindb, stmt)

  # t4 is a table that keeps a list of associations between related node ids
  stmt = "CREATE TABLE IF NOT EXISTS t4 (nodeid INTEGER NOT NULL, relatedid INTEGER NOT NULL, PRIMARY KEY (nodeid, relatedid))"
  sqlite.query(maindb, stmt)

  # t5 is a table that keeps a list of hierarchies between node ids.
  stmt = "CREATE TABLE IF NOT EXISTS t5 (nodeid INTEGER NOT NULL, upid INTEGER NOT NULL, PRIMARY KEY (nodeid, upid))"
  sqlite.query(maindb, stmt)

  # t6 is a table that keeps a list of all the derived words for each stem in the table 1 as `nodename`
  stmt = "CREATE TABLE IF NOT EXISTS t6 (nodeid INTEGER NOT NULL, derivative TEXT NOT NULL COLLATE NOCASE, PRIMARY KEY (nodeid, derivative))"
  sqlite.query(maindb, stmt)

  return maindb
end

function add_t1!(maindb::SQLite.DB, nodename)
  nodeid = table_length(maindb, "t1") + 1
  nodename = replace(nodename, "'" => "''")
  try
    stmt = "INSERT INTO t1 (nodeid, nodename) VALUES ($nodeid, '$nodename')"
    sqlite.query(maindb, stmt)
  catch
    stmt = "SELECT nodeid FROM t1 WHERE nodename = '$nodename' "
    nodeid = DataFrame(sqlite.query(maindb, stmt))[1,1]
    # nodeid -= 1
  end
  return nodeid
end

function add_t2!(maindb::SQLite.DB, descr, ref)
  descr = replace(descr, "'" => "''")
  stmt = "INSERT OR IGNORE INTO t2 (descr, ref) VALUES ('$descr', '$ref')"
  sqlite.query(maindb, stmt)
  descrid = DataFrame(sqlite.query(maindb, "SELECT MAX(descrid) FROM t2"))[1,1]
  return descrid
end

function add_t3!(maindb::SQLite.DB, nodeid::Integer, descrid::Integer)
  stmt = "INSERT OR IGNORE INTO t3 (nodeid, descrid) VALUES ($nodeid, $descrid)"
  sqlite.query(maindb, stmt)
end

function add_t4!(maindb::SQLite.DB, nodeid::Integer, relatedid::Integer)
  stmt = "INSERT OR IGNORE INTO t4 (nodeid, relatedid) VALUES ($nodeid, $relatedid)"
  sqlite.query(maindb, stmt)
end

function add_t5!(maindb::SQLite.DB, nodeid::Integer, upid::Integer)
  stmt = "INSERT OR IGNORE INTO t5 (nodeid, upid) VALUES ($nodeid, $upid)"
  sqlite.query(maindb, stmt)
end

function add_t6!(maindb::SQLite.DB, nodeid::Integer, derivative)
  derivative = replace(derivative, "'" => "''")
  stmt = "INSERT OR IGNORE INTO t6 (nodeid, derivative) VALUES ($nodeid, '$derivative')"
  sqlite.query(maindb, stmt)
end

"populate all tables"
function add_proposition!(maindb::SQLite.DB, proposition, nodes, stems, refs)
  nnodes = length(nodes)
  all_nodeids = Array{Integer}(undef, nnodes)
  for nn in 1:nnodes
    nodeid = add_t1!(maindb, stems[nn])
    all_nodeids[nn] = nodeid
    if stems[nn] != nodes[nn]
      add_t6!(maindb, nodeid, nodes[nn])
    end
    descrid = add_t2!(maindb, proposition, refs)
    add_t3!(maindb, nodeid, descrid)
  end
  if length(all_nodeids) > 1
    combs = combinations(all_nodeids, 2)
    for comb in combs
      add_t4!(maindb, comb[1], comb[2])
    end
  end
end

# TODO update t5

function table_length(maindb::SQLite.DB, tablename)
  # nrows = sqlite.execute!(maindb, "SELECT Count(1) FROM $tablename")
  nrows = DataFrame(sqlite.query(maindb, "SELECT Count(1) FROM $tablename"))[1,1]
  return nrows
end

## Functions to retrieve information about a concept##

function derivatives(maindb, nodeid)
  stmt = "SELECT derivative FROM t6 WHERE nodeid = $nodeid"
  result = DataFrame(sqlite.query(maindb, stmt))
  if size(result)[1] == 0
    return ""
  else
    ddd = result[1]
    return join(ddd, ", ")
  end
end

"Returns node ID of a query"
function find_nodeid(maindb, query)
  stmt = "SELECT nodeid FROM t1 WHERE nodename = '$query'"
  result = DataFrame(sqlite.query(maindb, stmt))
  if size(result)[1] == 0 # if not found search t6
    stmt = "SELECT nodeid FROM t6 WHERE derivative = '$query'"
    result = DataFrame(sqlite.query(maindb, stmt))
    if size(result)[1] == 0
      @info "Query is not in the database"
      return
    end
  end
  nodeid = result[1,1]
  return nodeid
end

"Returns IDs if related nodes to a given ID"
function related_concepts(maindb, nodeid)
  stmt = "SELECT relatedid FROM t4 WHERE nodeid = $nodeid"
  result = DataFrame(sqlite.query(maindb, stmt))
  if size(result)[1] == 0
    return ""
  else
    nmatches = length(result[1])
    nodenames = Array{String}(undef, nmatches)
    for m in 1:nmatches
      idd = result[1][m]
      stmt = "SELECT nodename FROM t1 WHERE nodeid = $idd"
      match = DataFrame(sqlite.query(maindb, stmt))[1,1]
      nodenames[m] = match
    end
    return join(nodenames, ", ")
  end
end

function find_descrid(maindb, nodeid)
  stmt = "SELECT descrid FROM t3 WHERE nodeid = $nodeid"
  result = DataFrame(sqlite.query(maindb, stmt))[1]
  return result
end

function descriptions(maindb, nodeid)
  descrids = find_descrid(maindb, nodeid)
  ndescrs = length(descrids)
  descrs = Array{String}(undef, ndescrs)
  refs = Array{String}(undef, ndescrs)
  for ind in 1:ndescrs
    stmt = "SELECT descr, ref FROM t2 WHERE descrid = $(descrids[ind])"
    result = DataFrame(sqlite.query(maindb, stmt))
    descrs[ind] = result[1,1]
    refs[ind] = result[1,2]
  end
  return descrs, refs
end

"Returns descrids of descriptions of relations between node1id and node2id"
function relationid(maindb, node1id::Integer, node2id::Integer)
  stmt = "SELECT descrid FROM t3 WHERE nodeid = $node1id INTERSECT SELECT descrid FROM t3 WHERE nodeid = $node2id"
  results = DataFrame(sqlite.query(maindb, stmt))[1]
  return results
end

function relation_descr(maindb, node1::AbstractString, node2::AbstractString)
  node1id = find_nodeid(maindb, node1)
  node2id = find_nodeid(maindb, node2)
  relationids = relationid(maindb, node1id, node2id)
  ndescr = length(relationids)
  descrs = Array{Tuple}(undef, ndescr)
  for nn in 1:ndescr
    descr = descriptions(maindb, relationids[nn])
    descrs[nn] =  descr
  end
  return descrs
end


text:
include("../../../../examples/monte_carlo_pi.jl")
using Test
@test abs(πₐₚₚᵣₒₓ - π) / π < 1e-2


text:
using OrdinaryDiffEq, DiffEqCallbacks, Test

function some_dynamics(u,p,t)
    du = zeros(length(u))
    for i in 1:length(u)
        for j in 1:length(u)
            du[i] += p[i,j]*u[j]*u[i]
        end
    end
    return du
end

p = rand(4,4)
startp = copy(p)

u0 = zeros(length(p[1,:])) .+ 0.1
tspan = (0.,1.)
prob = ODEProblem(some_dynamics,u0,tspan,p)
cb = PresetTimeCallback(0.5,integrator -> integrator.p .= rand(4,4))
sol = solve(prob,Tsit5(),callback = cb)
@test 0.5 ∈ sol.t
@test p != startp

p = rand(4,4)
startp = copy(p)

prob = ODEProblem(some_dynamics,u0,tspan,p)
cb = PresetTimeCallback([0.3,0.6],integrator -> integrator.p .= rand(4,4))
sol = solve(prob,Tsit5(),callback = cb)
@test 0.3 ∈ sol.t
@test 0.6 ∈ sol.t
@test p != startp


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

using Markdown
using InteractiveUtils

# ╔═╡ 6fd616f6-1ea1-11eb-3814-8bfb4a096c49
using Pkg, DrWatson

# ╔═╡ 6ffe2628-1ea1-11eb-24ea-57f985146a72
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Distributions, Random, StatsBase, StatsPlots
	Random.seed!(1)
end;

# ╔═╡ ed174bc4-1ea0-11eb-1e2f-a32874cec549
md"## Listing 3.11"

# ╔═╡ 70179b62-1ea1-11eb-13ce-7fb4d54bad4b
dist1 = TriangularDist(0,10,5)

# ╔═╡ 75c334f0-2134-11eb-132d-dfcb59e14474
dist2 = DiscreteUniform(1,5)

# ╔═╡ 75c37638-2134-11eb-150e-a5788c7f72a7
theorMean1, theorMean2 = mean(dist1), mean(dist2)

# ╔═╡ 75c3fa78-2134-11eb-0259-699e1c89a023
begin
	N = 10^6
	data1 = rand(dist1,N)
	data2 = rand(dist2,N)
	estMean1, estMean2 = mean(data1), mean(data2)
end

# ╔═╡ de66cfee-2134-11eb-2c7a-add4bee67341
density(data1, lab="Triangukar density")

# ╔═╡ a315bbee-2134-11eb-0249-17f3c6af3543
density!(data2, lab="Discrete Uniform density")

# ╔═╡ 75cedfbc-2134-11eb-3b3c-5bdc69e4735a
Text("Symmetric Triangular Distiribution on [0,10] has mean $theorMean1
	(estimated: $estMean1)")

# ╔═╡ 75cf627a-2134-11eb-0dc9-e5e3d869f4a8
Text("Discrete Uniform Distiribution on {1,2,3,4,5} has mean $theorMean2
	(estimated: $estMean2)")

# ╔═╡ 70182cd8-1ea1-11eb-094d-8d0d49cf15f3
md"## End of listing 3.11"

# ╔═╡ Cell order:
# ╟─ed174bc4-1ea0-11eb-1e2f-a32874cec549
# ╠═6fd616f6-1ea1-11eb-3814-8bfb4a096c49
# ╠═6ffe2628-1ea1-11eb-24ea-57f985146a72
# ╠═70179b62-1ea1-11eb-13ce-7fb4d54bad4b
# ╠═75c334f0-2134-11eb-132d-dfcb59e14474
# ╠═75c37638-2134-11eb-150e-a5788c7f72a7
# ╠═75c3fa78-2134-11eb-0259-699e1c89a023
# ╠═de66cfee-2134-11eb-2c7a-add4bee67341
# ╠═a315bbee-2134-11eb-0249-17f3c6af3543
# ╠═75cedfbc-2134-11eb-3b3c-5bdc69e4735a
# ╠═75cf627a-2134-11eb-0dc9-e5e3d869f4a8
# ╟─70182cd8-1ea1-11eb-094d-8d0d49cf15f3


text:
# Generate gsw_check_data.jl from the NetCDF file gsw_data_v3_0.nc
# downloaded from https://raw.githubusercontent.com/TEOS-10/GSW-C/master/gsw_data_v3_0.nc
# if it is not aleady present

# The code is adapted from https://github.com/TEOS-10/GSW-C/blob/master/make_check_data.py

using NCDatasets

datafile = joinpath(dirname(@__FILE__),"gsw_check_data.jl")

fname = joinpath(dirname(@__FILE__),"gsw_data_v3_0.nc")
if !isfile(fname)
    fname = download("https://raw.githubusercontent.com/TEOS-10/GSW-C/master/gsw_data_v3_0.nc")
end

ds = Dataset(fname)

work_vars = [
    ["ct", "CT_chck_cast"],
    ["rt", "Rt_chck_cast"],
    ["sa", "SA_chck_cast"],
    ["sk", "SK_chck_cast"],
    ["sp", "SP_chck_cast"],
    ["t", "t_chck_cast"],
    ["p", "p_chck_cast"],
    ["delta_p", "delta_p_chck_cast"],
    ["p_shallow", "p_chck_cast_shallow"],
    ["p_deep", "p_chck_cast_deep"],
    ["lat_cast", "lat_chck_cast"],
    ["long_cast", "long_chck_cast"],
    ["ct_arctic", "CT_Arctic"],
    ["sa_arctic", "SA_Arctic"],
    ["t_arctic", "t_Arctic"],
    ["p_arctic", "p_Arctic"],
    ["sa_seaice", "SA_seaice"],
    ["t_seaice", "t_seaice"],
    ["w_seaice", "w_seaice"],
    ["t_ice", "t_ice"],
    ["w_ice", "w_ice"],
    ["sa_bulk", "SA_bulk"],
    ["h_pot_bulk", "h_pot_bulk"],
    ["h_bulk", "h_bulk"],
    ["pref", "pr"]]

vars = [
    ["C_from_SP", ""],
    ["SP_from_C", ""],
    ["SP_from_SK", ""],
    ["SA_from_SP", ""],
    ["Sstar_from_SP", ""],
    ["CT_from_t", ""],
    ["deltaSA_from_SP", ""],
    ["SR_from_SP", ""],
    ["SP_from_SR", ""],
    ["SP_from_SA", ""],
    ["Sstar_from_SA", ""],
    ["SA_from_Sstar", ""],
    ["SP_from_Sstar", ""],
    ["pt_from_CT", ""],
    ["t_from_CT", ""],
    ["CT_from_pt", ""],
    ["pt0_from_t", ""],
    ["pt_from_t", ""],
    ["z_from_p", ""],
    ["entropy_from_pt", ""],
    ["pt_from_entropy", ""],
    ["CT_from_entropy", ""],
    ["entropy_from_CT", ""],
    ["entropy_from_t", ""],
    ["adiabatic_lapse_rate_from_CT", ""],
    ["specvol", ""],
    ["alpha", ""],
    ["beta", ""],
    ["alpha_on_beta", ""],
    ["v_vab", ""],
    ["alpha_vab", ""],
    ["beta_vab", ""],
    ["v_SA", ""],
    ["v_CT", ""],
    ["v_P", ""],
    ["v_SA_SA", ""],
    ["v_SA_CT", ""],
    ["v_CT_CT", ""],
    ["v_SA_P", ""],
    ["v_CT_P", ""],
    ["v_SA_wrt_h", ""],
    ["v_h", ""],
    ["v_SA_SA_wrt_h", ""],
    ["v_SA_h", ""],
    ["v_h_h", ""],
    ["specvol_anom_standard", ""],
    ["rho", ""],
    ["rho_rab", ""],
    ["alpha_rab", ""],
    ["beta_rab", ""],
    ["rho_SA", ""],
    ["rho_CT", ""],
    ["rho_P", ""],
    ["rho_SA_SA", ""],
    ["rho_SA_CT", ""],
    ["rho_CT_CT", ""],
    ["rho_SA_P", ""],
    ["rho_CT_P", ""],
    ["rho_SA_wrt_h", ""],
    ["rho_h", ""],
    ["rho_SA_SA_wrt_h", ""],
    ["rho_SA_h", ""],
    ["rho_h_h", ""],
    ["sigma0", ""],
    ["sigma1", ""],
    ["sigma2", ""],
    ["sigma3", ""],
    ["sigma4", ""],
    ["sound_speed", ""],
    ["kappa", ""],
    ["cabbeling", ""],
    ["thermobaric", ""],
    ["SA_from_rho", ""],
    ["CT_from_rho", ""],
    ["CT_maxdensity", ""],
    ["internal_energy", ""],
    ["enthalpy", ""],
    ["enthalpy_diff", ""],
    ["CT_from_enthalpy", ""],
    ["dynamic_enthalpy", ""],
    ["h_SA", ""],
    ["h_CT", ""],
    ["h_SA_SA", ""],
    ["h_SA_CT", ""],
    ["h_CT_CT", ""],
    ["CT_SA", ""],
    ["CT_pt", ""],
    ["CT_SA_SA", ""],
    ["CT_SA_pt", ""],
    ["CT_pt_pt", ""],
    ["eta_SA", ""],
    ["eta_CT", ""],
    ["eta_SA_SA", ""],
    ["eta_SA_CT", ""],
    ["eta_CT_CT", ""],
    ["pt_SA", ""],
    ["pt_CT", ""],
    ["pt_SA_SA", ""],
    ["pt_SA_CT", ""],
    ["pt_CT_CT", ""],
    ["CT_freezing", ""],
    ["CT_freezing_poly", ""],
    ["t_freezing", ""],
    ["t_freezing_poly", ""],
    ["pot_enthalpy_ice_freezing", ""],
    ["pot_enthalpy_ice_freezing_poly", ""],
    ["SA_freezing_from_CT", ""],
    ["SA_freezing_from_CT_poly", ""],
    ["SA_freezing_from_t", ""],
    ["SA_freezing_from_t_poly", ""],
    ["CTfreezing_SA", ""],
    ["CTfreezing_P", ""],
    ["CTfreezing_SA_poly", ""],
    ["CTfreezing_P_poly", ""],
    ["tfreezing_SA", ""],
    ["tfreezing_P", ""],
    ["tfreezing_SA_poly", ""],
    ["tfreezing_P_poly", ""],
    ["pot_enthalpy_ice_freezing_SA", ""],
    ["pot_enthalpy_ice_freezing_P", ""],
    ["pot_enthalpy_ice_freezing_SA_poly", ""],
    ["pot_enthalpy_ice_freezing_P_poly", ""],
    ["latentheat_melting", ""],
    ["latentheat_evap_CT", ""],
    ["latentheat_evap_t", ""],
    ["grav", ""],
    ["enthalpy_CT_exact", ""],
    ["h_SA_CT_exact", ""],
    ["h_CT_CT_exact", ""],
    ["h_SA_SA_CT_exact", ""],
    ["h_SA_CT_CT_exact", ""],
    ["h_CT_CT_CT_exact", ""],
    ["rho_t_exact", ""],
    ["pot_rho_t_exact", ""],
    ["alpha_wrt_t_exact", ""],
    ["beta_const_t_exact", ""],
    ["specvol_t_exact", ""],
    ["sound_speed_t_exact", ""],
    ["kappa_t_exact", ""],
    ["enthalpy_t_exact", ""],
    ["CT_SA_wrt_t", ""],
    ["CT_T_wrt_t", ""],
    ["CT_P_wrt_t", ""],
    ["chem_potential_water_t_exact", ""],
    ["t_deriv_chem_potential_water_t_exact", ""],
    ["dilution_coefficient_t_exact", ""],
    ["deltaSA_atlas", ""],
    ["Fdelta", ""],
    ["n2", ""],
    ["p_mid_n2", ""],
    ["Tu", ""],
    ["Rsubrho", ""],
    ["p_mid_TuRsr", ""],
    ["IPVfN2", ""],
    ["p_mid_IPVfN2", ""],
    ["geo_strf_dyn_height", ""],
    ["geo_strf_dyn_height_pc", ""],
    ["geo_strf_dyn_height_pc_p_mid", ""],
    ["rho_ice", ""],
    ["alpha_wrt_t_ice", ""],
    ["specvol_ice", ""],
    ["pressure_coefficient_ice", ""],
    ["sound_speed_ice", ""],
    ["kappa_ice", ""],
    ["kappa_const_t_ice", ""],
    ["internal_energy_ice", ""],
    ["enthalpy_ice", ""],
    ["entropy_ice", ""],
    ["cp_ice", ""],
    ["chem_potential_water_ice", ""],
    ["Helmholtz_energy_ice", ""],
    ["adiabatic_lapse_rate_ice", ""],
    ["pt0_from_t_ice", ""],
    ["pt_from_t_ice", ""],
    ["t_from_pt0_ice", ""],
    ["pot_enthalpy_from_pt_ice", ""],
    ["pt_from_pot_enthalpy_ice", ""],
    ["pot_enthalpy_from_pt_ice_poly", ""],
    ["pt_from_pot_enthalpy_ice_poly", ""],
    ["pressure_freezing_CT", ""],
    ["melting_ice_SA_CT_ratio", ""],
    ["melting_ice_SA_CT_ratio_poly", ""],
    ["melting_ice_equilibrium_SA_CT_ratio", ""],
    ["melting_ice_equilibrium_SA_CT_ratio_poly", ""],
    ["melting_ice_into_seawater_SA_final", ""],
    ["melting_ice_into_seawater_CT_final", ""],
    ["ice_fraction_to_freeze_seawater_SA_freeze", ""],
    ["ice_fraction_to_freeze_seawater_CT_freeze", ""],
    ["ice_fraction_to_freeze_seawater_w_Ih", ""],
    ["dSA_dCT_frazil", ""],
    ["dSA_dP_frazil", ""],
    ["dCT_dP_frazil", ""],
    ["dSA_dCT_frazil_poly", ""],
    ["dSA_dP_frazil_poly", ""],
    ["dCT_dP_frazil_poly", ""],
    ["frazil_properties_potential_SA_final", ""],
    ["frazil_properties_potential_CT_final", ""],
    ["frazil_properties_potential_w_Ih_final", ""],
    ["frazil_properties_potential_poly_SA_final", ""],
    ["frazil_properties_potential_poly_CT_final", ""],
    ["frazil_properties_potential_poly_w_Ih_final", ""],
    ["frazil_properties_SA_final", ""],
    ["frazil_properties_CT_final", ""],
    ["frazil_properties_w_Ih_final", ""],
    ["melting_seaice_SA_CT_ratio", ""],
    ["melting_seaice_SA_CT_ratio_poly", ""],
    ["melting_seaice_equilibrium_SA_CT_ratio", ""],
    ["melting_seaice_equilibrium_SA_CT_ratio_poly", ""],
    ["melting_seaice_into_seawater_SA_final", ""],
    ["melting_seaice_into_seawater_CT_final", ""],
    ["seaice_fraction_to_freeze_seawater_SA_freeze", ""],
    ["seaice_fraction_to_freeze_seawater_CT_freeze", ""],
    ["seaice_fraction_to_freeze_seawater_w_Ih", ""]
]

open(datafile,"w") do f
    println(f,"# This file is automatically generated by $(basename(@__FILE__))\n\n")

    for (codename,varname) in vcat(work_vars,vars)
        if varname == ""
            varname = codename
        end

        codename = lowercase(codename)
        data = ds[varname][:];

        print(f,"# $(varname)\n\n")
        print(f,"$(codename) = ")
        show(f,data)
        println(f,"\n")

        ca = get(ds[varname].attrib,"computation_accuracy",nothing)
        if ca != nothing
            print(f,"$(codename)_ca = $(ca)")
            println(f)
        end

        println(f,"\n\n")
    end
end

close(ds)


text:
export agg_loglikelihood

function agg_loglikelihood(icnf::AbstractICNF{T}, xs::AbstractMatrix{T}; agg::Function=mean)::T where {T <: AbstractFloat}
    logp̂x = inference(icnf, TestMode(), xs)
    agg(logp̂x)
end

function agg_loglikelihood(icnf::AbstractCondICNF{T}, xs::AbstractMatrix{T}, ys::AbstractMatrix{T}; agg::Function=mean)::T where {T <: AbstractFloat}
    logp̂x = inference(icnf, TestMode(), xs, ys)
    agg(logp̂x)
end


text:
#HW2 Problem 3
#Minhao Yan my497, Wentong Chen wc422

#Firm's Profit Function
function profit(η,c1,c2,q)
    profit1=sum(q)^(-1/η)*q[1]-c1*q[1] # Firm 1's profit
    profit2=sum(q)^(-1/η)*q[2]-c2*q[2] # Firm 2's profit
    profit=[profit1; profit2]
    return profit
end

#Function that solves q_prime
function q_prime(η,c1,c2,q)
    q_prime11=(-1/η)*sum(q)^(-1/η-1)-(1/η)*sum(q)^(-1/η-1)+(1/η^2+1/η)*q[1]*sum(q)^(-1/η-2)-c1
    q_prime12=(-1/η)*sum(q)^(-1/η-1)+(1/η^2+1/η)*sum(q)^(-1/η-2)*q[1]
    q_prime21=(-1/η)*sum(q)^(-1/η-1)+(1/η^2+1/η)*sum(q)^(-1/η-2)*q[2]
    q_prime22=(-1/η)*sum(q)^(-1/η-1)-(1/η)*sum(q)^(-1/η-1)+(1/η^2+1/η)*q[2]*sum(q)^(-1/η-2)-c2
    q_prime_m=[q_prime11 q_prime12; q_prime21 q_prime22]
    return q_prime_m
end

#Newton's method
function cournot_newton(η,c1,c2,initial_guess,tolerance)
    diff=[Inf; Inf]
    q_old=initial_guess
    q=[1e10;1e10]

    while abs.(diff)>tolerance
        q = q_old - transpose(transpose(profit(η,c1,c2,q_old))*inv(q_prime(η,c1,c2,q_old)))
        diff= q- q_old
        q_old= q
    end
 #  Print out final output
 println("The solution is $q.")
end

#   Parameters and Inputs
η=1.6
c1=0.15
c2=0.2
initial_guess=[20;20]
tolerance=[1e-5;1e-5]

#   Optimization
q_opt=cournot_newton(η,c1,c2,initial_guess,tolerance)


text:
module DEACfitness
using Statistics
using LinearAlgebra
const functionNames = ["isfonly",
                       "isfMn1",
                       "isfM0",
                       "isfMn1M0",
                       "isfMn1M0M1",
                       "isfMn1M0M1M3",
                       "Mn1M0M1",
                       "Mn1M0M1M3"]
const fitMethodNames = ["prodFit",
                        "sumFit",
                        "maxFit"]
function prodabsFit(fT::Float64)
    1.0 + fT;
end

function prodabsFit(fT::Float64,moment::Float64)
    1.0 + abs(moment - fT);
end

function absFit(fT::Float64,moment::Float64)
    abs(moment - fT);
end

function prodFit(f1::Float64,f2::Float64,f3::Float64,f4::Float64,f5::Float64,momentn1::Float64,moment0::Float64,moment1::Float64,moment3::Float64)
    prodabsFit(f1) * prodabsFit(f2,momentn1) * prodabsFit(f3,moment0) * prodabsFit(f4,moment1) * prodabsFit(f5,moment3) - 1.0;
end

function sumFit(f1::Float64,f2::Float64,f3::Float64,f4::Float64,f5::Float64,momentn1::Float64,moment0::Float64,moment1::Float64,moment3::Float64)
    f1 + absFit(f2,momentn1) + absFit(f3,moment0) + absFit(f4,moment1) + absFit(f5,moment3);
end

function sumFitDebug(f1::Float64,f2::Float64,f3::Float64,f4::Float64,f5::Float64,momentn1::Float64,moment0::Float64,moment1::Float64,moment3::Float64)
    t1 = f1;
    t2 = absFit(f2,momentn1);
    t3 = absFit(f3,moment0);
    t4 = absFit(f4,moment1);
    t5 = absFit(f5,moment3);
    println("t1: $(t1)");
    println("t2: $(t2)");
    println("t3: $(t3)");
    println("t4: $(t4)");
    println("t5: $(t5)");
    t1 + t2 + t3 + t4 + t5;
end

function maxFit(f1::Float64,f2::Float64,f3::Float64,f4::Float64,f5::Float64,momentn1::Float64,moment0::Float64,moment1::Float64,moment3::Float64)
    max(absFit(f1),absFit(f2),absFit(f3),absFit(f4),absFit(f5));
end

function isfonly(fmfunc::Function,fitnessP::Array{Float64,1},P::Array{Float64,2},isf::Array{Float64,1},momentn1::Float64,moment0::Float64,moment1::Float64,moment3::Float64,F_term_P1::Array{Float64,2},F_term_P2::Array{Float64,2},F_term_avg::Array{Float64,1},mn1_term_P1::Array{Float64,1},mn1_term_P2::Array{Float64,1},m0_term_P1::Array{Float64,1},m0_term_P2::Array{Float64,1},m1_term_P1::Array{Float64,1},m1_term_P2::Array{Float64,1},m3_term_P1::Array{Float64,1},m3_term_P2::Array{Float64,1},isfmT_P1::Array{Float64,2},mn1T_P1::Array{Float64,1},m0T_P1::Array{Float64,1},m1T_P1::Array{Float64,1},m3T_P1::Array{Float64,1},isfmT_P2::Array{Float64,2},mn1T_P2::Array{Float64,1},m0T_P2::Array{Float64,1},m1T_P2::Array{Float64,1},m3T_P2::Array{Float64,1})
    fitInd = trues(size(P,2));
    isfonly(fmfunc,fitnessP,fitInd,P,isf,momentn1,moment0,moment1,moment3,F_term_P1,F_term_P2,F_term_avg,mn1_term_P1,mn1_term_P2,m0_term_P1,m0_term_P2,m1_term_P1,m1_term_P2,m3_term_P1,m3_term_P2,isfmT_P1,mn1T_P1,m0T_P1,m1T_P1,m3T_P1,isfmT_P2,mn1T_P2,m0T_P2,m1T_P2,m3T_P2)
end

function isfonly(fmfunc::Function,fitnessP::Array{Float64,1},fitInd::BitArray{1},P::Array{Float64,2},isf::Array{Float64,1},momentn1::Float64,moment0::Float64,moment1::Float64,moment3::Float64,F_term_P1::Array{Float64,2},F_term_P2::Array{Float64,2},F_term_avg::Array{Float64,1},mn1_term_P1::Array{Float64,1},mn1_term_P2::Array{Float64,1},m0_term_P1::Array{Float64,1},m0_term_P2::Array{Float64,1},m1_term_P1::Array{Float64,1},m1_term_P2::Array{Float64,1},m3_term_P1::Array{Float64,1},m3_term_P2::Array{Float64,1},isfmT_P1::Array{Float64,2},mn1T_P1::Array{Float64,1},m0T_P1::Array{Float64,1},m1T_P1::Array{Float64,1},m3T_P1::Array{Float64,1},isfmT_P2::Array{Float64,2},mn1T_P2::Array{Float64,1},m0T_P2::Array{Float64,1},m1T_P2::Array{Float64,1},m3T_P2::Array{Float64,1})
    mul!(F_term_P1,isfmT_P1,P[1:end-1,:]);
    mul!(F_term_P2,isfmT_P2,P[2:end,:]);
    mean!(F_term_avg,abs.(isf .- (F_term_P1 .+ F_term_P2))');
    f1 = 0.0;
    f2 = momentn1;
    f3 = moment0;
    f4 = moment1;
    f5 = moment3;
    @inbounds for j = 1:size(fitInd,1)
        if fitInd[j]
            f1 = F_term_avg[j];
            fitnessP[j] = fmfunc(f1,f2,f3,f4,f5,momentn1,moment0,moment1,moment3); #testnewfit4
        end
    end
    nothing
end

function isfMn1(fmfunc::Function,fitnessP::Array{Float64,1},P::Array{Float64,2},isf::Array{Float64,1},momentn1::Float64,moment0::Float64,moment1::Float64,moment3::Float64,F_term_P1::Array{Float64,2},F_term_P2::Array{Float64,2},F_term_avg::Array{Float64,1},mn1_term_P1::Array{Float64,1},mn1_term_P2::Array{Float64,1},m0_term_P1::Array{Float64,1},m0_term_P2::Array{Float64,1},m1_term_P1::Array{Float64,1},m1_term_P2::Array{Float64,1},m3_term_P1::Array{Float64,1},m3_term_P2::Array{Float64,1},isfmT_P1::Array{Float64,2},mn1T_P1::Array{Float64,1},m0T_P1::Array{Float64,1},m1T_P1::Array{Float64,1},m3T_P1::Array{Float64,1},isfmT_P2::Array{Float64,2},mn1T_P2::Array{Float64,1},m0T_P2::Array{Float64,1},m1T_P2::Array{Float64,1},m3T_P2::Array{Float64,1})
    fitInd = trues(size(P,2));
    isfMn1(fmfunc,fitnessP,fitInd,P,isf,momentn1,moment0,moment1,moment3,F_term_P1,F_term_P2,F_term_avg,mn1_term_P1,mn1_term_P2,m0_term_P1,m0_term_P2,m1_term_P1,m1_term_P2,m3_term_P1,m3_term_P2,isfmT_P1,mn1T_P1,m0T_P1,m1T_P1,m3T_P1,isfmT_P2,mn1T_P2,m0T_P2,m1T_P2,m3T_P2)
end

function isfMn1(fmfunc::Function,fitnessP::Array{Float64,1},fitInd::BitArray{1},P::Array{Float64,2},isf::Array{Float64,1},momentn1::Float64,moment0::Float64,moment1::Float64,moment3::Float64,F_term_P1::Array{Float64,2},F_term_P2::Array{Float64,2},F_term_avg::Array{Float64,1},mn1_term_P1::Array{Float64,1},mn1_term_P2::Array{Float64,1},m0_term_P1::Array{Float64,1},m0_term_P2::Array{Float64,1},m1_term_P1::Array{Float64,1},m1_term_P2::Array{Float64,1},m3_term_P1::Array{Float64,1},m3_term_P2::Array{Float64,1},isfmT_P1::Array{Float64,2},mn1T_P1::Array{Float64,1},m0T_P1::Array{Float64,1},m1T_P1::Array{Float64,1},m3T_P1::Array{Float64,1},isfmT_P2::Array{Float64,2},mn1T_P2::Array{Float64,1},m0T_P2::Array{Float64,1},m1T_P2::Array{Float64,1},m3T_P2::Array{Float64,1})
    mul!(F_term_P1,isfmT_P1,P[1:end-1,:]);
    mul!(F_term_P2,isfmT_P2,P[2:end,:]);
    mean!(F_term_avg,abs.(isf .- (F_term_P1 .+ F_term_P2))');
    mul!(mn1_term_P1,P[1:end-1,:]',mn1T_P1);
    mul!(mn1_term_P2,P[2:end,:]',mn1T_P2);
    f1 = 0.0;
    f2 = 0.0;
    f3 = moment0;
    f4 = moment1;
    f5 = moment3;
    @inbounds for j = 1:size(fitInd,1)
        if fitInd[j]
            f1 = F_term_avg[j];
            f2 = mn1_term_P1[j]+mn1_term_P2[j];
            fitnessP[j] = fmfunc(f1,f2,f3,f4,f5,momentn1,moment0,moment1,moment3); #testnewfit4
        end
    end
    nothing
end

function isfM0(fmfunc::Function,fitnessP::Array{Float64,1},P::Array{Float64,2},isf::Array{Float64,1},momentn1::Float64,moment0::Float64,moment1::Float64,moment3::Float64,F_term_P1::Array{Float64,2},F_term_P2::Array{Float64,2},F_term_avg::Array{Float64,1},mn1_term_P1::Array{Float64,1},mn1_term_P2::Array{Float64,1},m0_term_P1::Array{Float64,1},m0_term_P2::Array{Float64,1},m1_term_P1::Array{Float64,1},m1_term_P2::Array{Float64,1},m3_term_P1::Array{Float64,1},m3_term_P2::Array{Float64,1},isfmT_P1::Array{Float64,2},mn1T_P1::Array{Float64,1},m0T_P1::Array{Float64,1},m1T_P1::Array{Float64,1},m3T_P1::Array{Float64,1},isfmT_P2::Array{Float64,2},mn1T_P2::Array{Float64,1},m0T_P2::Array{Float64,1},m1T_P2::Array{Float64,1},m3T_P2::Array{Float64,1})
    fitInd = trues(size(P,2));
    isfM0(fmfunc,fitnessP,fitInd,P,isf,momentn1,moment0,moment1,moment3,F_term_P1,F_term_P2,F_term_avg,mn1_term_P1,mn1_term_P2,m0_term_P1,m0_term_P2,m1_term_P1,m1_term_P2,m3_term_P1,m3_term_P2,isfmT_P1,mn1T_P1,m0T_P1,m1T_P1,m3T_P1,isfmT_P2,mn1T_P2,m0T_P2,m1T_P2,m3T_P2)
end

function isfM0(fmfunc::Function,fitnessP::Array{Float64,1},fitInd::BitArray{1},P::Array{Float64,2},isf::Array{Float64,1},momentn1::Float64,moment0::Float64,moment1::Float64,moment3::Float64,F_term_P1::Array{Float64,2},F_term_P2::Array{Float64,2},F_term_avg::Array{Float64,1},mn1_term_P1::Array{Float64,1},mn1_term_P2::Array{Float64,1},m0_term_P1::Array{Float64,1},m0_term_P2::Array{Float64,1},m1_term_P1::Array{Float64,1},m1_term_P2::Array{Float64,1},m3_term_P1::Array{Float64,1},m3_term_P2::Array{Float64,1},isfmT_P1::Array{Float64,2},mn1T_P1::Array{Float64,1},m0T_P1::Array{Float64,1},m1T_P1::Array{Float64,1},m3T_P1::Array{Float64,1},isfmT_P2::Array{Float64,2},mn1T_P2::Array{Float64,1},m0T_P2::Array{Float64,1},m1T_P2::Array{Float64,1},m3T_P2::Array{Float64,1})
    mul!(F_term_P1,isfmT_P1,P[1:end-1,:]);
    mul!(F_term_P2,isfmT_P2,P[2:end,:]);
    mean!(F_term_avg,abs.(isf .- (F_term_P1 .+ F_term_P2))');
    mul!(m0_term_P1,P[1:end-1,:]',m0T_P1);
    mul!(m0_term_P2,P[2:end,:]',m0T_P2);
    f1 = 0.0;
    f2 = momentn1;
    f3 = 0.0;
    f4 = moment1;
    f5 = moment3;
    @inbounds for j = 1:size(fitInd,1)
        if fitInd[j]
            f1 = F_term_avg[j];
            f3 = m0_term_P1[j]+m0_term_P2[j];
            fitnessP[j] = fmfunc(f1,f2,f3,f4,f5,momentn1,moment0,moment1,moment3); #testnewfit4
        end
    end
    nothing
end

function isfMn1M0(fmfunc::Function,fitnessP::Array{Float64,1},P::Array{Float64,2},isf::Array{Float64,1},momentn1::Float64,moment0::Float64,moment1::Float64,moment3::Float64,F_term_P1::Array{Float64,2},F_term_P2::Array{Float64,2},F_term_avg::Array{Float64,1},mn1_term_P1::Array{Float64,1},mn1_term_P2::Array{Float64,1},m0_term_P1::Array{Float64,1},m0_term_P2::Array{Float64,1},m1_term_P1::Array{Float64,1},m1_term_P2::Array{Float64,1},m3_term_P1::Array{Float64,1},m3_term_P2::Array{Float64,1},isfmT_P1::Array{Float64,2},mn1T_P1::Array{Float64,1},m0T_P1::Array{Float64,1},m1T_P1::Array{Float64,1},m3T_P1::Array{Float64,1},isfmT_P2::Array{Float64,2},mn1T_P2::Array{Float64,1},m0T_P2::Array{Float64,1},m1T_P2::Array{Float64,1},m3T_P2::Array{Float64,1})
    fitInd = trues(size(P,2));
    isfMn1M0(fmfunc,fitnessP,fitInd,P,isf,momentn1,moment0,moment1,moment3,F_term_P1,F_term_P2,F_term_avg,mn1_term_P1,mn1_term_P2,m0_term_P1,m0_term_P2,m1_term_P1,m1_term_P2,m3_term_P1,m3_term_P2,isfmT_P1,mn1T_P1,m0T_P1,m1T_P1,m3T_P1,isfmT_P2,mn1T_P2,m0T_P2,m1T_P2,m3T_P2)
end

function isfMn1M0(fmfunc::Function,fitnessP::Array{Float64,1},fitInd::BitArray{1},P::Array{Float64,2},isf::Array{Float64,1},momentn1::Float64,moment0::Float64,moment1::Float64,moment3::Float64,F_term_P1::Array{Float64,2},F_term_P2::Array{Float64,2},F_term_avg::Array{Float64,1},mn1_term_P1::Array{Float64,1},mn1_term_P2::Array{Float64,1},m0_term_P1::Array{Float64,1},m0_term_P2::Array{Float64,1},m1_term_P1::Array{Float64,1},m1_term_P2::Array{Float64,1},m3_term_P1::Array{Float64,1},m3_term_P2::Array{Float64,1},isfmT_P1::Array{Float64,2},mn1T_P1::Array{Float64,1},m0T_P1::Array{Float64,1},m1T_P1::Array{Float64,1},m3T_P1::Array{Float64,1},isfmT_P2::Array{Float64,2},mn1T_P2::Array{Float64,1},m0T_P2::Array{Float64,1},m1T_P2::Array{Float64,1},m3T_P2::Array{Float64,1})
    mul!(F_term_P1,isfmT_P1,P[1:end-1,:]);
    mul!(F_term_P2,isfmT_P2,P[2:end,:]);
    mean!(F_term_avg,abs.(isf .- (F_term_P1 .+ F_term_P2))');
    mul!(mn1_term_P1,P[1:end-1,:]',mn1T_P1);
    mul!(mn1_term_P2,P[2:end,:]',mn1T_P2);
    mul!(m0_term_P1,P[1:end-1,:]',m0T_P1);
    mul!(m0_term_P2,P[2:end,:]',m0T_P2);
    f1 = 0.0;
    f2 = 0.0;
    f3 = 0.0;
    f4 = moment1;
    f5 = moment3;
    @inbounds for j = 1:size(fitInd,1)
        if fitInd[j]
            f1 = F_term_avg[j];
            f2 = mn1_term_P1[j]+mn1_term_P2[j];
            f3 = m0_term_P1[j]+m0_term_P2[j];
            fitnessP[j] = fmfunc(f1,f2,f3,f4,f5,momentn1,moment0,moment1,moment3); #testnewfit4
        end
    end
    nothing
end

function isfMn1M0M1(fmfunc::Function,fitnessP::Array{Float64,1},P::Array{Float64,2},isf::Array{Float64,1},momentn1::Float64,moment0::Float64,moment1::Float64,moment3::Float64,F_term_P1::Array{Float64,2},F_term_P2::Array{Float64,2},F_term_avg::Array{Float64,1},mn1_term_P1::Array{Float64,1},mn1_term_P2::Array{Float64,1},m0_term_P1::Array{Float64,1},m0_term_P2::Array{Float64,1},m1_term_P1::Array{Float64,1},m1_term_P2::Array{Float64,1},m3_term_P1::Array{Float64,1},m3_term_P2::Array{Float64,1},isfmT_P1::Array{Float64,2},mn1T_P1::Array{Float64,1},m0T_P1::Array{Float64,1},m1T_P1::Array{Float64,1},m3T_P1::Array{Float64,1},isfmT_P2::Array{Float64,2},mn1T_P2::Array{Float64,1},m0T_P2::Array{Float64,1},m1T_P2::Array{Float64,1},m3T_P2::Array{Float64,1})
    fitInd = trues(size(P,2));
    isfMn1M0M1(fmfunc,fitnessP,fitInd,P,isf,momentn1,moment0,moment1,moment3,F_term_P1,F_term_P2,F_term_avg,mn1_term_P1,mn1_term_P2,m0_term_P1,m0_term_P2,m1_term_P1,m1_term_P2,m3_term_P1,m3_term_P2,isfmT_P1,mn1T_P1,m0T_P1,m1T_P1,m3T_P1,isfmT_P2,mn1T_P2,m0T_P2,m1T_P2,m3T_P2)
end

function isfMn1M0M1(fmfunc::Function,fitnessP::Array{Float64,1},fitInd::BitArray{1},P::Array{Float64,2},isf::Array{Float64,1},momentn1::Float64,moment0::Float64,moment1::Float64,moment3::Float64,F_term_P1::Array{Float64,2},F_term_P2::Array{Float64,2},F_term_avg::Array{Float64,1},mn1_term_P1::Array{Float64,1},mn1_term_P2::Array{Float64,1},m0_term_P1::Array{Float64,1},m0_term_P2::Array{Float64,1},m1_term_P1::Array{Float64,1},m1_term_P2::Array{Float64,1},m3_term_P1::Array{Float64,1},m3_term_P2::Array{Float64,1},isfmT_P1::Array{Float64,2},mn1T_P1::Array{Float64,1},m0T_P1::Array{Float64,1},m1T_P1::Array{Float64,1},m3T_P1::Array{Float64,1},isfmT_P2::Array{Float64,2},mn1T_P2::Array{Float64,1},m0T_P2::Array{Float64,1},m1T_P2::Array{Float64,1},m3T_P2::Array{Float64,1})
    mul!(F_term_P1,isfmT_P1,P[1:end-1,:]);
    mul!(F_term_P2,isfmT_P2,P[2:end,:]);
    mean!(F_term_avg,abs.(isf .- (F_term_P1 .+ F_term_P2))');
    mul!(mn1_term_P1,P[1:end-1,:]',mn1T_P1);
    mul!(mn1_term_P2,P[2:end,:]',mn1T_P2);
    mul!(m0_term_P1,P[1:end-1,:]',m0T_P1);
    mul!(m0_term_P2,P[2:end,:]',m0T_P2);
    mul!(m1_term_P1,P[1:end-1,:]',m1T_P1);
    mul!(m1_term_P2,P[2:end,:]',m1T_P2);
    f1 = 0.0;
    f2 = 0.0;
    f3 = 0.0;
    f4 = 0.0;
    f5 = moment3;
    @inbounds for j = 1:size(fitInd,1)
        if fitInd[j]
            f1 = F_term_avg[j];
            f2 = mn1_term_P1[j]+mn1_term_P2[j];
            f3 = m0_term_P1[j]+m0_term_P2[j];
            f4 = m1_term_P1[j]+m1_term_P2[j];
            fitnessP[j] = fmfunc(f1,f2,f3,f4,f5,momentn1,moment0,moment1,moment3); #testnewfit4
        end
    end
    nothing
end

function isfMn1M0M1M3(fmfunc::Function,fitnessP::Array{Float64,1},P::Array{Float64,2},isf::Array{Float64,1},momentn1::Float64,moment0::Float64,moment1::Float64,moment3::Float64,F_term_P1::Array{Float64,2},F_term_P2::Array{Float64,2},F_term_avg::Array{Float64,1},mn1_term_P1::Array{Float64,1},mn1_term_P2::Array{Float64,1},m0_term_P1::Array{Float64,1},m0_term_P2::Array{Float64,1},m1_term_P1::Array{Float64,1},m1_term_P2::Array{Float64,1},m3_term_P1::Array{Float64,1},m3_term_P2::Array{Float64,1},isfmT_P1::Array{Float64,2},mn1T_P1::Array{Float64,1},m0T_P1::Array{Float64,1},m1T_P1::Array{Float64,1},m3T_P1::Array{Float64,1},isfmT_P2::Array{Float64,2},mn1T_P2::Array{Float64,1},m0T_P2::Array{Float64,1},m1T_P2::Array{Float64,1},m3T_P2::Array{Float64,1})
    fitInd = trues(size(P,2));
    isfMn1M0M1M3(fmfunc,fitnessP,fitInd,P,isf,momentn1,moment0,moment1,moment3,F_term_P1,F_term_P2,F_term_avg,mn1_term_P1,mn1_term_P2,m0_term_P1,m0_term_P2,m1_term_P1,m1_term_P2,m3_term_P1,m3_term_P2,isfmT_P1,mn1T_P1,m0T_P1,m1T_P1,m3T_P1,isfmT_P2,mn1T_P2,m0T_P2,m1T_P2,m3T_P2)
end

function isfMn1M0M1M3(fmfunc::Function,fitnessP::Array{Float64,1},fitInd::BitArray{1},P::Array{Float64,2},isf::Array{Float64,1},momentn1::Float64,moment0::Float64,moment1::Float64,moment3::Float64,F_term_P1::Array{Float64,2},F_term_P2::Array{Float64,2},F_term_avg::Array{Float64,1},mn1_term_P1::Array{Float64,1},mn1_term_P2::Array{Float64,1},m0_term_P1::Array{Float64,1},m0_term_P2::Array{Float64,1},m1_term_P1::Array{Float64,1},m1_term_P2::Array{Float64,1},m3_term_P1::Array{Float64,1},m3_term_P2::Array{Float64,1},isfmT_P1::Array{Float64,2},mn1T_P1::Array{Float64,1},m0T_P1::Array{Float64,1},m1T_P1::Array{Float64,1},m3T_P1::Array{Float64,1},isfmT_P2::Array{Float64,2},mn1T_P2::Array{Float64,1},m0T_P2::Array{Float64,1},m1T_P2::Array{Float64,1},m3T_P2::Array{Float64,1})
    mul!(F_term_P1,isfmT_P1,P[1:end-1,:]);
    mul!(F_term_P2,isfmT_P2,P[2:end,:]);
    mean!(F_term_avg,abs.(isf .- (F_term_P1 .+ F_term_P2))');
    mul!(mn1_term_P1,P[1:end-1,:]',mn1T_P1);
    mul!(mn1_term_P2,P[2:end,:]',mn1T_P2);
    mul!(m0_term_P1,P[1:end-1,:]',m0T_P1);
    mul!(m0_term_P2,P[2:end,:]',m0T_P2);
    mul!(m1_term_P1,P[1:end-1,:]',m1T_P1);
    mul!(m1_term_P2,P[2:end,:]',m1T_P2);
    mul!(m3_term_P1,P[1:end-1,:]',m3T_P1);
    mul!(m3_term_P2,P[2:end,:]',m3T_P2);
    f1 = 0.0;
    f2 = 0.0;
    f3 = 0.0;
    f4 = 0.0;
    f5 = 0.0;
    @inbounds for j = 1:size(fitInd,1)
        if fitInd[j]
            f1 = F_term_avg[j];
            f2 = mn1_term_P1[j]+mn1_term_P2[j];
            f3 = m0_term_P1[j]+m0_term_P2[j];
            f4 = m1_term_P1[j]+m1_term_P2[j];
            f5 = m3_term_P1[j]+m3_term_P2[j];
            fitnessP[j] = fmfunc(f1,f2,f3,f4,f5,momentn1,moment0,moment1,moment3); #testnewfit4
        end
    end
    nothing
end

function Mn1M0M1(fmfunc::Function,fitnessP::Array{Float64,1},P::Array{Float64,2},isf::Array{Float64,1},momentn1::Float64,moment0::Float64,moment1::Float64,moment3::Float64,F_term_P1::Array{Float64,2},F_term_P2::Array{Float64,2},F_term_avg::Array{Float64,1},mn1_term_P1::Array{Float64,1},mn1_term_P2::Array{Float64,1},m0_term_P1::Array{Float64,1},m0_term_P2::Array{Float64,1},m1_term_P1::Array{Float64,1},m1_term_P2::Array{Float64,1},m3_term_P1::Array{Float64,1},m3_term_P2::Array{Float64,1},isfmT_P1::Array{Float64,2},mn1T_P1::Array{Float64,1},m0T_P1::Array{Float64,1},m1T_P1::Array{Float64,1},m3T_P1::Array{Float64,1},isfmT_P2::Array{Float64,2},mn1T_P2::Array{Float64,1},m0T_P2::Array{Float64,1},m1T_P2::Array{Float64,1},m3T_P2::Array{Float64,1})
    fitInd = trues(size(P,2));
    Mn1M0M1(fmfunc,fitnessP,fitInd,P,isf,momentn1,moment0,moment1,moment3,F_term_P1,F_term_P2,F_term_avg,mn1_term_P1,mn1_term_P2,m0_term_P1,m0_term_P2,m1_term_P1,m1_term_P2,m3_term_P1,m3_term_P2,isfmT_P1,mn1T_P1,m0T_P1,m1T_P1,m3T_P1,isfmT_P2,mn1T_P2,m0T_P2,m1T_P2,m3T_P2)
end

function Mn1M0M1(fmfunc::Function,fitnessP::Array{Float64,1},fitInd::BitArray{1},P::Array{Float64,2},isf::Array{Float64,1},momentn1::Float64,moment0::Float64,moment1::Float64,moment3::Float64,F_term_P1::Array{Float64,2},F_term_P2::Array{Float64,2},F_term_avg::Array{Float64,1},mn1_term_P1::Array{Float64,1},mn1_term_P2::Array{Float64,1},m0_term_P1::Array{Float64,1},m0_term_P2::Array{Float64,1},m1_term_P1::Array{Float64,1},m1_term_P2::Array{Float64,1},m3_term_P1::Array{Float64,1},m3_term_P2::Array{Float64,1},isfmT_P1::Array{Float64,2},mn1T_P1::Array{Float64,1},m0T_P1::Array{Float64,1},m1T_P1::Array{Float64,1},m3T_P1::Array{Float64,1},isfmT_P2::Array{Float64,2},mn1T_P2::Array{Float64,1},m0T_P2::Array{Float64,1},m1T_P2::Array{Float64,1},m3T_P2::Array{Float64,1})
    mul!(mn1_term_P1,P[1:end-1,:]',mn1T_P1);
    mul!(mn1_term_P2,P[2:end,:]',mn1T_P2);
    mul!(m0_term_P1,P[1:end-1,:]',m0T_P1);
    mul!(m0_term_P2,P[2:end,:]',m0T_P2);
    mul!(m1_term_P1,P[1:end-1,:]',m1T_P1);
    mul!(m1_term_P2,P[2:end,:]',m1T_P2);
    f1 = 0.0;
    f2 = 0.0;
    f3 = 0.0;
    f4 = 0.0;
    f5 = moment3;
    @inbounds for j = 1:size(fitInd,1)
        if fitInd[j]
            f2 = mn1_term_P1[j]+mn1_term_P2[j];
            f3 = m0_term_P1[j]+m0_term_P2[j];
            f4 = m1_term_P1[j]+m1_term_P2[j];
            fitnessP[j] = fmfunc(f1,f2,f3,f4,f5,momentn1,moment0,moment1,moment3); #testnewfit4
        end
    end
end

function Mn1M0M1M3(fmfunc::Function,fitnessP::Array{Float64,1},P::Array{Float64,2},isf::Array{Float64,1},momentn1::Float64,moment0::Float64,moment1::Float64,moment3::Float64,F_term_P1::Array{Float64,2},F_term_P2::Array{Float64,2},F_term_avg::Array{Float64,1},mn1_term_P1::Array{Float64,1},mn1_term_P2::Array{Float64,1},m0_term_P1::Array{Float64,1},m0_term_P2::Array{Float64,1},m1_term_P1::Array{Float64,1},m1_term_P2::Array{Float64,1},m3_term_P1::Array{Float64,1},m3_term_P2::Array{Float64,1},isfmT_P1::Array{Float64,2},mn1T_P1::Array{Float64,1},m0T_P1::Array{Float64,1},m1T_P1::Array{Float64,1},m3T_P1::Array{Float64,1},isfmT_P2::Array{Float64,2},mn1T_P2::Array{Float64,1},m0T_P2::Array{Float64,1},m1T_P2::Array{Float64,1},m3T_P2::Array{Float64,1})
    fitInd = trues(size(P,2));
    Mn1M0M1M3(fmfunc,fitnessP,fitInd,P,isf,momentn1,moment0,moment1,moment3,F_term_P1,F_term_P2,F_term_avg,mn1_term_P1,mn1_term_P2,m0_term_P1,m0_term_P2,m1_term_P1,m1_term_P2,m3_term_P1,m3_term_P2,isfmT_P1,mn1T_P1,m0T_P1,m1T_P1,m3T_P1,isfmT_P2,mn1T_P2,m0T_P2,m1T_P2,m3T_P2)
end

function Mn1M0M1M3(fmfunc::Function,fitnessP::Array{Float64,1},fitInd::BitArray{1},P::Array{Float64,2},isf::Array{Float64,1},momentn1::Float64,moment0::Float64,moment1::Float64,moment3::Float64,F_term_P1::Array{Float64,2},F_term_P2::Array{Float64,2},F_term_avg::Array{Float64,1},mn1_term_P1::Array{Float64,1},mn1_term_P2::Array{Float64,1},m0_term_P1::Array{Float64,1},m0_term_P2::Array{Float64,1},m1_term_P1::Array{Float64,1},m1_term_P2::Array{Float64,1},m3_term_P1::Array{Float64,1},m3_term_P2::Array{Float64,1},isfmT_P1::Array{Float64,2},mn1T_P1::Array{Float64,1},m0T_P1::Array{Float64,1},m1T_P1::Array{Float64,1},m3T_P1::Array{Float64,1},isfmT_P2::Array{Float64,2},mn1T_P2::Array{Float64,1},m0T_P2::Array{Float64,1},m1T_P2::Array{Float64,1},m3T_P2::Array{Float64,1})
    mul!(mn1_term_P1,P[1:end-1,:]',mn1T_P1);
    mul!(mn1_term_P2,P[2:end,:]',mn1T_P2);
    mul!(m0_term_P1,P[1:end-1,:]',m0T_P1);
    mul!(m0_term_P2,P[2:end,:]',m0T_P2);
    mul!(m1_term_P1,P[1:end-1,:]',m1T_P1);
    mul!(m1_term_P2,P[2:end,:]',m1T_P2);
    mul!(m3_term_P1,P[1:end-1,:]',m3T_P1);
    mul!(m3_term_P2,P[2:end,:]',m3T_P2);
    f1 = 0.0;
    f2 = 0.0;
    f3 = 0.0;
    f4 = 0.0;
    f5 = 0.0;
    @inbounds for j = 1:size(fitInd,1)
        if fitInd[j]
            f2 = mn1_term_P1[j]+mn1_term_P2[j];
            f3 = m0_term_P1[j]+m0_term_P2[j];
            f4 = m1_term_P1[j]+m1_term_P2[j];
            f5 = m3_term_P1[j]+m3_term_P2[j];
            fitnessP[j] = fmfunc(f1,f2,f3,f4,f5,momentn1,moment0,moment1,moment3); #testnewfit4
        end
    end
    nothing
end
end


text:
using Tools
using Documenter

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

makedocs(;
    modules=[Tools],
    authors="KronosTheLate",
    repo="https://github.com/KronosTheLate/Tools.jl/blob/{commit}{path}#{line}",
    sitename="Tools.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://KronosTheLate.github.io/Tools.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/KronosTheLate/Tools.jl",
    devbranch="master",
)


text:
mutable struct FileCachedCovariate{T<:AbstractFloat, S<:AbstractCovariate{T}} <: AbstractCovariate{T}
    cachepath::String
    basecovariate::S
    lockobj::Threads.TatasLock

    FileCachedCovariate{T, S}(cachepath, basecovariate, lockobj) where {S<:AbstractCovariate{T}} where {T<:AbstractFloat} = begin
        res = new(cachepath, basecovariate, lockobj)
        finalizer(res, x -> begin
            try
                if isfile(x.cachepath)
                    rm(x.cachepath)
                end
            finally
            end
            end)
        res
    end
end

Base.length(var::FileCachedCovariate) = length(var.basecovariate)

getname(var::FileCachedCovariate) = getname(var.basecovariate)

function FileCachedCovariate(covariate::S) where {S<:AbstractCovariate{T}} where {T<:AbstractFloat}
    FileCachedCovariate{T, S}("", covariate, Threads.TatasLock())
end

function filecache(covariate::S) where {S<:AbstractCovariate{T}} where {T<:AbstractFloat}
    FileCachedCovariate{T, S}("", covariate, Threads.TatasLock())
end

function slice(covariate::FileCachedCovariate{T, S}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T<:AbstractFloat} where {S<:Union{FileCachedCovariate, CachedCovariate, Covariate, FileCovariate}}
    basecovariate = covariate.basecovariate
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    slice(basecovariate, fromobs, toobs, slicelength)
end

function slice(covariate::FileCachedCovariate{T, S}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {S<:AbstractCovariate{T}} where {T<:AbstractFloat}
    basecovariate = covariate.basecovariate
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    lockobj = covariate.lockobj
    lock(lockobj)
    try
        if !isfile(covariate.cachepath)
            v = convert(Vector{T}, basecovariate)
            tmpjuml = joinpath(tempdir(), JUMLDIR)
            mkpath(tmpjuml)
            outdatpath = joinpath(tmpjuml, "$(randstring(10)).dat")
            open(outdatpath, "w") do f
                write(f, v)
            end
            covariate.cachepath = outdatpath
        end
    finally
        unlock(lockobj)
    end
    sizeT = sizeof(T)
    iostream = open(covariate.cachepath, "r")
    seek(iostream, sizeT * (fromobs - 1))
    buffer = Vector{T}(undef, slicelength)
    map(Seq(Vector{T}, (buffer, iostream, fromobs, toobs, slicelength), nextdatachunk), SubArray{T,1,Array{T,1},Tuple{UnitRange{Int64}},true}) do slice
        view(slice, 1:length(slice))
    end
end

function Base.map(covariate::FileCachedCovariate, dataframe::AbstractDataFrame)
    FileCachedCovariate(map(covariate.basecovariate, dataframe))
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
# -*- mode: julia; -*-
module PulseAudio
using Reexport
using FixedPointNumbers
# reexport so users automatically get SampledSignals API
@reexport using SampledSignals

export PulseAudioSink

# typedef enum pa_sample_format
const PA_SAMPLE_U8        =  0 # Unsigned 8 Bit PCM
const PA_SAMPLE_ALAW      =  1 # 8 Bit a-Law
const PA_SAMPLE_ULAW      =  2 # 8 Bit mu-Law
const PA_SAMPLE_S16LE     =  3 # Signed 16 Bit PCM, little endian (PC)
const PA_SAMPLE_S16BE     =  4 # Signed 16 Bit PCM, big endian
const PA_SAMPLE_FLOAT32LE =  5 # 32 Bit IEEE floating point, little endian (PC), range -1.0 to 1.0
const PA_SAMPLE_FLOAT32BE =  6 # 32 Bit IEEE floating point, big endian, range -1.0 to 1.0
const PA_SAMPLE_S32LE     =  7 # Signed 32 Bit PCM, little endian (PC)
const PA_SAMPLE_S32BE     =  8 # Signed 32 Bit PCM, big endian
const PA_SAMPLE_S24LE     =  9 # Signed 24 Bit PCM packed, little endian (PC). \since 0.9.15
const PA_SAMPLE_S24BE     = 10 # Signed 24 Bit PCM packed, big endian. \since 0.9.15
const PA_SAMPLE_S24_32LE  = 11 # Signed 24 Bit PCM in LSB of 32 Bit words, little endian (PC). \since 0.9.15
const PA_SAMPLE_S24_32BE  = 12 # Signed 24 Bit PCM in LSB of 32 Bit words, big endian. \since 0.9.15

# map julia types onto PulseAudio stream types
# yep, we're just assuming little-endianness
specformat(::Type{Float32}) = PA_SAMPLE_FLOAT32LE
specformat(::Type{UInt8}) = PA_SAMPLE_U8
specformat(::Type{Int16}) = PA_SAMPLE_S16LE
specformat(::Type{Fixed{Int16, 15}}) = PA_SAMPLE_S16LE
specformat(::Type{Int32}) = PA_SAMPLE_S32LE
specformat(::Type{Fixed{Int32, 31}}) = PA_SAMPLE_S32LE
specformat(T) = error("Element type $T not supported")

immutable pa_sample_spec
    format::Int32
    rate::UInt32
    channels::UInt8
end

immutable pa_channel_map
    channels::UInt8

    # map data (max 32 channels)
    map0::Cint
    map1::Cint
    map2::Cint
    map3::Cint
    map4::Cint
    map5::Cint
    map6::Cint
    map7::Cint
    map8::Cint
    map9::Cint
    map10::Cint
    map11::Cint
    map12::Cint
    map13::Cint
    map14::Cint
    map15::Cint
    map16::Cint
    map17::Cint
    map18::Cint
    map19::Cint
    map20::Cint
    map21::Cint
    map22::Cint
    map23::Cint
    map24::Cint
    map25::Cint
    map26::Cint
    map27::Cint
    map28::Cint
    map29::Cint
    map30::Cint
    map31::Cint

    pa_channel_map() = new(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
end

immutable pa_buffer_attr
    maxlength::UInt32
    tlength::UInt32
    prebuf::UInt32
    minreq::UInt32
    fragsize::UInt32
end

const pa_simple = Ptr{Void}
const LibPulseSimple = "libpulse-simple"
const PA_STREAM_PLAYBACK = 1
const PA_CHANNEL_MAP_AIFF = 0
const PA_CHANNEL_MAP_DEFAULT = PA_CHANNEL_MAP_AIFF

mutable struct PulseAudioSink{T} <: SampleSink
    pulsesink::pa_simple
    spec::pa_sample_spec

    function PulseAudioSink{T}(ch, fs, name, description) where T
        eltype = specformat(T)

        spec = pa_sample_spec(eltype, fs, ch)
        s = ccall((:pa_simple_new, LibPulseSimple),
                  pa_simple,
                  (Cstring, Cstring, Cint, Cstring, Cstring,
                   Ptr{pa_sample_spec}, Ptr{pa_channel_map}, Ptr{pa_buffer_attr},
                   Ptr{Cint}),
                  C_NULL, # Use the default server
                  name,  # Application name
                  PA_STREAM_PLAYBACK,
                  C_NULL, # Use the default device
                  description, # description of stream
                  Ref(spec),
                  C_NULL, # Use default channel map
                  C_NULL, # Use default buffering attributes
                  C_NULL) # Ignore error code
        if s == C_NULL
            error("pa_simple_new failed")
        end

        instance = new(s, spec)
        finalizer(instance, close)

        instance
    end
end

"""
    PulseAudioSink(; channels=2, samplerate=48000, eltype=Float32,
                   name="Julia", description="PulseAudioSink")

Create an audio sink that you can write to. The sink is samplerate-aware and
interoperates with the SampledSignals ecosystem, so channel mapping and
samplerate conversion are handled automatically.

Example:

```
sink = PulseAudioSink()
write(sink, sin.(2pi*220*linspace(0,1,48000))*0.2)
```
"""
function PulseAudioSink(; channels=2, samplerate=48000, eltype=Float32,
                        name="Julia", description="PulseAudioSink")
    PulseAudioSink{eltype}(channels, samplerate, name, description)
end

SampledSignals.nchannels(sink::PulseAudioSink) = Int(sink.spec.channels)
SampledSignals.samplerate(sink::PulseAudioSink) = Float64(sink.spec.rate)
Base.eltype(sink::PulseAudioSink{T}) where T = T

function Base.close(sink::PulseAudioSink)
    if sink.pulsesink != C_NULL
        ccall((:pa_simple_free, LibPulseSimple), Void, (pa_simple,), sink.pulsesink)
        sink.pulsesink = C_NULL
    end

    nothing
end

function SampledSignals.unsafe_write(sink::PulseAudioSink, buf::Array,
                                     frameoffset, framecount)
    # pulseaudio wants interleaved data, so transpose
    data = buf[(1:framecount) + frameoffset, :]'
    write_ret = ccall((:pa_simple_write, LibPulseSimple),
                      Cint,
                      (pa_simple, Ptr{Void}, Csize_t, Ptr{Cint}),
                      sink.pulsesink, data, sizeof(data), C_NULL)
    if write_ret != 0
        error("pa_simple_write failed with $write_ret")
    end

    return framecount
end

end # module


text:
using Blink

using MeshCat
using CoordinateTransformations
using Rotations
using GeometryBasics: HyperRectangle, Vec, Point, Mesh, Point3f0
using Colors: RGBA, RGB

if (false)
    vis = MeshCat.Visualizer()
    win = Blink.Window()
    open(vis, win)
    Blink.AtomShell.opentools(win)
end

# box = HyperRectangle(Vec(0., 0, 0), Vec(1., 1, 1))
# setobject!(vis, box)
# settransform!(vis, Translation(0., 1, 0))

# green_box_vis = setobject!(vis["group1"]["greenbox"], box, green_material)
# settransform!(green_box_vis, Translation(0, 0, 1))
# group1 = vis["group1"]
# settransform!(group1, Translation(0, 0, -1))

delete!(vis)

v = vis[:lines]
settransform!(v, Translation(-1, -1, 0))
points = [
    Point3f0(1, 1, 1),
    Point3f0(2, 2, 2),
    Point3f0(1, 1, 3),
    Point3f0(2, 2, 3),
]
colors = [
    RGBA(1, 0, 0, 1),
    RGBA(1, 0, 0, 1),
    RGBA(0, 1, 0, 1),
    RGBA(0, 1, 0, 1),
]
mat = LineBasicMaterial(
    color=RGBA(1.0, 1.0, 0.0, 1.0),
    linewidth = 11,
    # vertexColors = 1,
)
setobject!(v[:line_segments], LineSegments(points, colors))

# v = vis[:lines]
# settransform!(v, Translation(-1, -1, 0))
# θ = range(0, stop=2π, length=10)
# mat = LineBasicMaterial(
#     color=RGBA(1.0, 1.0, 0.0, 1.0),
#     linewidth = 11,
#     # vertexColors = 1,
# )
# setobject!(v[:line_segments], LineSegments(Point.(0.5 .* sin.(θ), 0, 0.5 .* cos.(θ)), mat))

# box = HyperRectangle(Vec(0., 0, 0), Vec(1., 1, 1))
# setobject!(vis["box"], box)
# settransform!(vis, Translation(0., 2, 0))

# vbox = vis["box"]

# va1 = vbox["arrow1"]
# arrow1 = ArrowVisualizer(va1)
# setobject!(arrow1)
# settransform!(arrow1, Point(5.0, 5.0, 5.0), Vec(-4.0, -4.0, -4.0), shaft_radius=0.1)

# va2 = vbox["arrow2"]
# arrow2 = ArrowVisualizer(va2)
# setobject!(arrow2)
# settransform!(arrow2, Point(-1.0, 2.0, 2.0), Vec(1.0, -1.0, -1.0))

# va3 = vbox["arrow3"]
# arrow3 = ArrowVisualizer(va3)
# setobject!(arrow3)
# settransform!(arrow3, Point(-1.0, -1.0, 2.0), Vec(1.0, 1.0, -1.0))

# va4 = vbox["arrow4"]
# arrow4 = ArrowVisualizer(va4)
# setobject!(arrow4)
# settransform!(arrow4, Point(2.0, -1.0, 2.0), Vec(-1.0, 1.0, -1.0))