text:
# camera model

"""
Data structure for a Camera model with parameters.
Use `CameraModel(width,height,fc,cc,skew,kc)` for easy construction.
"""
struct CameraModelandParameters
    width::Int		# image width
    height::Int		# image height
    fc::Vector{Float64}	# focal length in x and y
    cc::Vector{Float64}	# camera center
    skew::Float64	    # skew value
    kc::Vector{Float64} # distortion coefficients up to fifth order
    K::Matrix{Float64} # 3x3 camera calibration matrix (Camera intrinsics)
    Ki::Matrix{Float64} # inverse of a 3x3 camera calibratio matrix
end

"""
    CameraModel(width,height,fc,cc,skew,kc)

Constructor helper for creating a camera model.
"""
function CameraModel(width,height,fc,cc,skew,kc)
    KK = [fc[1]      skew  cc[1];
             0       fc[2] cc[2];
             0		    0     1]
    # KK = [fc[1] skew*fc[1] cc[1];
    #          0       fc[2] cc[2];
    #          0		    0     1]
    Ki = inv(KK)
    CameraModelandParameters(width,height,fc,cc,skew,kc,KK,Ki)
end

# Sensor Model and Geometry
mutable struct PInt64
    i::Int64
end

"""
Data structure for holding time, gyroscope, and accelerometer data.

**Fields**
 - `utime` -- time [μs]
 - `acc`   -- accelerometer data
 - `gyro`  -- gyroscope data
"""
struct IMU_DATA
	utime::Int64
	acc::Vector{Float64}
	gyro::Vector{Float64}
	# mag::Vector{Float64}
    # IMU_DATA() = new()
end

"""
	propR!(R, gyro_sample, delta_time)

Propogate rotation matrix R with gyro data over time dt
"""
function propR!(R::Matrix{Float64}, gyro::Vector{Float64}, dt::Float64)
    # this integration can be done better
    #alpha = 0.5*(gyro + datt.dalpha)*dt
    alpha = gyro*dt
    phi = alpha
    dR = expm(skew(phi))
    R[:] = R * dR
	return nothing
end

"""
    HornAbsoluteOrientation(a::Matrix{Float64},b::Matrix{Float64})
    HornAbsoluteOrientation(a::Vector{Vector{Float64}},b::Vector{Vector{Float64}})

Compute the rotation between rows of (a and b)::Array{Float64,2}.\\
Rotate b into the frame of a\\
Returns a quaternion, aQb
"""
function HornAbsoluteOrientation(a::Matrix{Float64},b::Matrix{Float64})
  # rotate b into the frame of a

  N = zeros(4,4);

  for n = 1:size(a,1)
    li = [0;a[n,:]]
    ri = [0;b[n,:]]

    N = N + linmapR(ri)'*linmapL(li)
  end

  D,v = eig(N);
  d = diagm(D);
  # maxd = maximum(d);
  # maxind = find(maxd==d);
  maxd, maxind = findmax(D)

  q = v[:,maxind]
  if (q[1]<0)
    q = -q;
  end

  return Quaternion(q)
end

function HornAbsoluteOrientation(a::Vector{Vector{Float64}},b::Vector{Vector{Float64}})
  # rotate b into the frame of a

  N = zeros(4,4)
  valid_count = 0

  for n = 1:size(a,1)
    if isnan(a[n][1]) || isnan(b[n][1])
        continue
    end
    li = [0;a[n]]
    ri = [0;b[n]]

    N = N + linmapR(ri)'*linmapL(li)
    valid_count += 1;
  end

  if valid_count < 4
      error("Not enough vectors for a solution!")
      return Quaternion()
  end

  D,v = eig(N)
  d = diagm(D)
  # maxd = maximum(d);
  # maxind = find(maxd==d);
  maxd, maxind = findmax(D)

  q = v[:,maxind]
  if (q[1]<0)
    q = -q;
  end

  return Quaternion(q)
end

"""
    linmapL(p)
"""
function linmapL(p)
    p0 = p[1];
    px = p[2];
    py = p[3];
    pz = p[4];

    P = [p0  -px  -py -pz;
         px   p0  -pz  py;
         py   pz   p0 -px;
         pz  -py   px  p0]

    return P
end

"""
    linmapR(q)
"""
function linmapR(q)
  q0 = q[1];
  qx = q[2];
  qy = q[3];
  qz = q[4];

  Q = [q0  -qx  -qy  -qz;
       qx   q0   qz  -qy;
       qy  -qz   q0   qx;
       qz   qy  -qx   q0];

  return Q
end

"""
    integrateGyroBetweenFrames!(index, current_time, vector_data)

Estimate rotations from IMU data between time stamps.
"""
function integrateGyroBetweenFrames!(index::PInt64, ctime::Int64, imudata::Vector{IMU_DATA})

	R = eye(3)

	n = size(imudata,1)

	#end of imudata
	if (index.i > n)
 		return (false, eye(3))
	end

	#increment index to start at 2 for integration, ignore first
    (index.i < 2)? index.i+=1:nothing

	while (index.i <= n) && (imudata[index.i].utime <= ctime)

		dt = Float64(imudata[index.i].utime - imudata[index.i-1].utime)/1e6

		propR!(R, imudata[index.i].gyro, dt)

		# println("@$(index.i) : $(cR)")
		index.i += 1
	end

	return (true, R)

end

"""
	estimateRotationFromKeypoints(points_a, points_b, cameraModel, [compensate = false])

Estimate the rotation between 2 sets of Keypoints a and b using HornAbsoluteOrientation.\\
It is assumed that only possitive Keypoints are valid.\\
Set compensate to true if rotaion should be around centre off image rather than 0,0
"""
function estimateRotationFromKeypoints(points_a::Keypoints, points_b::Keypoints, cam::CameraModelandParameters; compensate = false)

	focald = cam.fc[1] #assume 1:1 aspect ratio for now TODO
	#create a Vector of 3d vectors for valid keypoints (ie. not 0) else NaN
	nps_a = map(kp -> (kp[1] > 0 < kp[2])? [kp[1], kp[2], focald] : [NaN,NaN,NaN], points_a)
	nps_b = map(kp -> (kp[1] > 0 < kp[2])? [kp[1], kp[2], focald] : [NaN,NaN,NaN], points_b)
    if compensate
        # compensate offset to centre of camera
        nps_a = map(nps -> nps - [cam.cc; 0.], nps_a)
        nps_b = map(nps -> nps - [cam.cc; 0.], nps_b)
    end
	#normalize vectors
	nps_a .= nps_a./norm.(nps_a)
	nps_b .= nps_b./norm.(nps_b)

	return HornAbsoluteOrientation(nps_a,nps_b)

end

"""
    predictHomographyIMU!(index, current_time, vector_data, CameraK, CameraK_inverse)

Estimate camera rotation from IMU data and compute predicted homography
"""
function predictHomographyIMU!(index::PInt64, ctime::Int64, imudata::Vector{IMU_DATA},K::Array{Float64,2},Ki::Array{Float64,2}; cRi::Array{Float64,2} = eye(3))

	R = eye(3)

	n = size(imudata,1)

	#end of imudata
	if (index.i >= n)
 		return (false, eye(3))
	end

	#increment index to start at 2 for integration, ignore first
    (index.i < 2)? index.i+=1:nothing

	while (index.i <= n) && (imudata[index.i].utime  <= ctime)

		dt = Float64(imudata[index.i].utime - imudata[index.i-1].utime)/1e6

		SensorFeatureTracking.propR!(R, cRi * imudata[index.i].gyro, dt)

		# println("@$(index.i) : $(cR)")
		index.i += 1
	end

	#compute homography
	# H = K*R*inv(K)
    R4 = eye(4)
    R4[1:3,1:3] = R
    K4 = eye(4)
    K4[1:2,1:2] = K[1:2,1:2]
    K4[1:2,4] = K[1:2,3]
    H = K4*R4*inv(K4)
# H./H[3,3] ??
    # H[:] ./= H[3,3]
	return (true, H);
end

function predictHomographyIMU(roll::Float64, pitch::Float64, yaw::Float64, K::Array{Float64,2}, Ki::Array{Float64,2})

    R = eye(4)
    R[1:3,1:3] = SO3(Euler(roll,pitch,yaw)).R

    K4 = eye(4)
    K4[1:2,1:2] = K[1:2,1:2]
    K4[1:2,4] = K[1:2,3]
    H = K4*R*inv(K4)
# H./H[3,3] ??
    return  K4*R*inv(K4)
end

function predictAffinity(H::Matrix{Float64})

    return AffineMap(SMatrix{2,2}(H[1:2,1:2]), SVector{2}(H[1:2,:]*[0,0,1,1]))

end

max_stars_repo_name: dehann/SensorFeatureTracking

max_issues_repo_name: dehann/SensorFeatureTracking

max_forks_repo_name: dehann/SensorFeatureTracking


text:
type ModelParams
    Zc::Float64
    R::Float64
    Rav::Float64
    C::Float64
    V0::Float64
    t1::Float64
    t2::Float64
    m1::Float64
    m2::Float64
    Emax::Float64
    Emin::Float64
    th::Vector{Float64}

    function ModelParams(old=Dict("a"=>0),restart="no")
        this = new()
        if restart == "no"
            # healthy
            # # this.Zc = 2e5; # Pa⋅s/m^3
            # # this.R = 2e7;
            # this.Zc = 1.186*2e5; # Pa⋅s/m^3
            # this.R = 0.940*3e7;
            # this.Rav = 2e5;
            # # this.Rav = 0;
            # # this.C = 1.3e-7 # m^3/Pa;
            # this.C = 0.901*1.5e-8 # m^3/Pa;
            # this.V0 = 10*cm3Tom3;
            # # this.t1 = 0.215;
            # this.t1 = 0.913*1;
            # # this.t2 = 0.362;
            # this.t2 = 1.308*0.3;
            # this.m1 = 0.821*2;
            # this.m2 = 0.389*27.4;
            # # this.m2 = 7;
            # # this.Emax = 7e7;
            # this.Emax = 0.968*1e8;
            # patient
            this.Zc = 1.007*2e5; # Pa⋅s/m^3
            this.R = 0.963*3e7;
            this.Rav = 2e5;
            this.C = 0.864*1e-8 # m^3/Pa;
            this.V0 = 10*cm3Tom3;
            this.t1 = 0.896*0.25;
            this.t2 = 0.4;
            this.m1 = 1.066*2;
            this.m2 = 1.037*20;
            this.Emax = 0.764*1.4e8;
            this.Emin = 3.77e6;
            this.th = [0.8];
        elseif restart == "yes"
            this.Zc = old["Zc"];
            this.R = old["R"];
            this.Rav = old["Rav"];
            this.C = old["C"]; # m^3/Pa
            this.V0 = old["V0"];
            this.t1 = old["t1"];
            this.t2 = old["t2"];
            this.m1 = old["m1"];
            this.m2 = old["m2"];
            this.Emax = old["Emax"];
            this.Emin = old["Emin"];
            this.th = [old["th"]];
        end
        return this
    end
end

ext: jl


text:
#
# Correctness Tests
#

using Optim
using Base.Test
using Compat

my_tests = [
    "types.jl",
    "bfgs.jl",
    "gradient_descent.jl",
    "momentum_gradient_descent.jl",
    "grid_search.jl",
    "l_bfgs.jl",
    "levenberg_marquardt.jl",
    "newton.jl",
    "newton_trust_region.jl",
    "cg.jl",
    "nelder_mead.jl",
    "optimize.jl",
    "simulated_annealing.jl",
    "particle_swarm.jl",
    "interpolating_line_search.jl",
    "api.jl",
    "golden_section.jl",
    "brent.jl",
    "type_stability.jl",
    "array.jl",
    "constrained.jl",
    "callbacks.jl",
    "precon.jl",
    "initial_convergence.jl"
]

println("Running tests:")

for my_test in my_tests
    println(" * $(my_test)")
    include(my_test)
end


text:
"""
# module FortranToJulia

Conversion of Fortran values (as strings) to equivalent Julia values.

The functions in this module are used to convert the string representation of
the values of basic Fortran data types into equivalent Julia values.

# Examples

```jldoctest
julia>
```
"""
module FortranToJulia

using Compat: isnothing

export FortranData, @f_str

struct FortranData{T <: AbstractString}
    data::T
end

macro f_str(str)
    return :(FortranData($str))
end

function Base.parse(::Type{T}, s::FortranData) where {T <: Integer}
    return parse(T, s.data)
end
function Base.parse(::Type{T}, s::FortranData) where {T <: Float32}
    return parse(T, replace(lowercase(s.data), r"(?<=[^e])(?=[+-])" => "f"))
end
function Base.parse(::Type{T}, s::FortranData) where {T <: Float64}
    return parse(T, replace(lowercase(s.data), r"d"i => "e"))
end
function Base.parse(::Type{Complex{T}}, s::FortranData) where {T <: AbstractFloat}
    str = s.data
    if first(str) == '(' && last(str) == ')' && length(split(str, ',')) == 2
        re, im = split(str[2:end - 1], ',', limit = 2)
        return Complex(parse(T, re), parse(T, im))
    else
        throw(Meta.ParseError("$str must be in complex number form (x, y)."))
    end
end
function Base.parse(::Type{Bool}, s::FortranData)
    str = lowercase(s.data)
    if str in (".true.", ".t.", "true", 't')
        return true
    elseif str in (".false.", ".f.", "false", 'f')
        return false
    else
        throw(Meta.ParseError("$str is not a valid logical constant."))
    end
end
function Base.parse(::Type{T}, s::FortranData) where {T <: AbstractString}
    str = s.data
    m = match(r"([\"'])((?:\\\1|.)*?)\1", str)
    isnothing(m) && throw(Meta.ParseError("$str is not a valid string!"))
    quotation_mark, content = m.captures
    # Replace escaped strings
    return string(replace(content, repeat(quotation_mark, 2) => quotation_mark))
end

end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

module DSP

import Base.trailingsize

export filt, filt!, deconv, conv, conv2, xcorr

_zerosi(b,a,T) = zeros(promote_type(eltype(b), eltype(a), T), max(length(a), length(b))-1)

"""
    filt(b, a, x, [si])

Apply filter described by vectors `a` and `b` to vector `x`, with an optional initial filter
state vector `si` (defaults to zeros).
"""
function filt{T,S}(b::Union{AbstractVector, Number}, a::Union{AbstractVector, Number},
                   x::AbstractArray{T}, si::AbstractArray{S}=_zerosi(b,a,T))
    filt!(Array{promote_type(eltype(b), eltype(a), T, S)}(size(x)), b, a, x, si)
end

# in-place filtering: returns results in the out argument, which may shadow x
# (and does so by default)

"""
    filt!(out, b, a, x, [si])

Same as [`filt`](:func:`filt`) but writes the result into the `out` argument, which may
alias the input `x` to modify it in-place.
"""
function filt!{T,S,N}(out::AbstractArray, b::Union{AbstractVector, Number}, a::Union{AbstractVector, Number},
                      x::AbstractArray{T}, si::AbstractArray{S,N}=_zerosi(b,a,T))
    isempty(b) && throw(ArgumentError("filter vector b must be non-empty"))
    isempty(a) && throw(ArgumentError("filter vector a must be non-empty"))
    a[1] == 0  && throw(ArgumentError("filter vector a[1] must be nonzero"))
    if size(x) != size(out)
        throw(ArgumentError("output size $(size(out)) must match input size $(size(x))"))
    end

    as = length(a)
    bs = length(b)
    sz = max(as, bs)
    silen = sz - 1
    ncols = trailingsize(x,2)

    if size(si, 1) != silen
        throw(ArgumentError("initial state vector si must have max(length(a),length(b))-1 rows"))
    end
    if N > 1 && trailingsize(si,2) != ncols
        throw(ArgumentError("initial state vector si must be a vector or have the same number of columns as x"))
    end

    size(x,1) == 0 && return out
    sz == 1 && return scale!(out, x, b[1]/a[1]) # Simple scaling without memory

    # Filter coefficient normalization
    if a[1] != 1
        norml = a[1]
        a ./= norml
        b ./= norml
    end
    # Pad the coefficients with zeros if needed
    bs<sz   && (b = copy!(zeros(eltype(b), sz), b))
    1<as<sz && (a = copy!(zeros(eltype(a), sz), a))

    initial_si = si
    for col = 1:ncols
        # Reset the filter state
        si = initial_si[:, N > 1 ? col : 1]
        if as > 1
            _filt_iir!(out, b, a, x, si, col)
        else
            _filt_fir!(out, b, x, si, col)
        end
    end
    return out
end

function _filt_iir!(out, b, a, x, si, col)
    silen = length(si)
    @inbounds for i=1:size(x, 1)
        xi = x[i,col]
        val = si[1] + b[1]*xi
        for j=1:(silen-1)
            si[j] = si[j+1] + b[j+1]*xi - a[j+1]*val
        end
        si[silen] = b[silen+1]*xi - a[silen+1]*val
        out[i,col] = val
    end
end

function _filt_fir!(out, b, x, si, col)
    silen = length(si)
    @inbounds for i=1:size(x, 1)
        xi = x[i,col]
        val = si[1] + b[1]*xi
        for j=1:(silen-1)
            si[j] = si[j+1] + b[j+1]*xi
        end
        si[silen] = b[silen+1]*xi
        out[i,col] = val
    end
end

"""
    deconv(b,a) -> c

Construct vector `c` such that `b = conv(a,c) + r`.
Equivalent to polynomial division.
"""
function deconv{T}(b::StridedVector{T}, a::StridedVector{T})
    lb = size(b,1)
    la = size(a,1)
    if lb < la
        return [zero(T)]
    end
    lx = lb-la+1
    x = zeros(T, lx)
    x[1] = 1
    filt(b, a, x)
end

"""
    conv(u,v)

Convolution of two vectors. Uses FFT algorithm.
"""
function conv{T<:Base.LinAlg.BlasFloat}(u::StridedVector{T}, v::StridedVector{T})
    nu = length(u)
    nv = length(v)
    n = nu + nv - 1
    np2 = n > 1024 ? nextprod([2,3,5], n) : nextpow2(n)
    upad = [u; zeros(T, np2 - nu)]
    vpad = [v; zeros(T, np2 - nv)]
    if T <: Real
        p = plan_rfft(upad)
        y = irfft((p*upad).*(p*vpad), np2)
    else
        p = plan_fft!(upad)
        y = ifft!((p*upad).*(p*vpad))
    end
    return y[1:n]
end
conv{T<:Integer}(u::StridedVector{T}, v::StridedVector{T}) = round(Int,conv(float(u), float(v)))
conv{T<:Integer, S<:Base.LinAlg.BlasFloat}(u::StridedVector{T}, v::StridedVector{S}) = conv(float(u), v)
conv{T<:Integer, S<:Base.LinAlg.BlasFloat}(u::StridedVector{S}, v::StridedVector{T}) = conv(u, float(v))

"""
    conv2(u,v,A)

2-D convolution of the matrix `A` with the 2-D separable kernel generated by
the vectors `u` and `v`.
Uses 2-D FFT algorithm.
"""
function conv2{T}(u::StridedVector{T}, v::StridedVector{T}, A::StridedMatrix{T})
    m = length(u)+size(A,1)-1
    n = length(v)+size(A,2)-1
    B = zeros(T, m, n)
    B[1:size(A,1),1:size(A,2)] = A
    u = fft([u;zeros(T,m-length(u))])
    v = fft([v;zeros(T,n-length(v))])
    C = ifft(fft(B) .* (u * v.'))
    if T <: Real
        return real(C)
    end
    return C
end

"""
    conv2(B,A)

2-D convolution of the matrix `B` with the matrix `A`. Uses 2-D FFT algorithm.
"""
function conv2{T}(A::StridedMatrix{T}, B::StridedMatrix{T})
    sa, sb = size(A), size(B)
    At = zeros(T, sa[1]+sb[1]-1, sa[2]+sb[2]-1)
    Bt = zeros(T, sa[1]+sb[1]-1, sa[2]+sb[2]-1)
    At[1:sa[1], 1:sa[2]] = A
    Bt[1:sb[1], 1:sb[2]] = B
    p = plan_fft(At)
    C = ifft((p*At).*(p*Bt))
    if T <: Real
        return real(C)
    end
    return C
end
conv2{T<:Integer}(A::StridedMatrix{T}, B::StridedMatrix{T}) = round(Int,conv2(float(A), float(B)))
conv2{T<:Integer}(u::StridedVector{T}, v::StridedVector{T}, A::StridedMatrix{T}) = round(Int,conv2(float(u), float(v), float(A)))

"""
    xcorr(u,v)

Compute the cross-correlation of two vectors.
"""
function xcorr(u, v)
    su = size(u,1); sv = size(v,1)
    if su < sv
        u = [u;zeros(eltype(u),sv-su)]
    elseif sv < su
        v = [v;zeros(eltype(v),su-sv)]
    end
    flipdim(conv(flipdim(u, 1), v), 1)
end

end # module

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
"Construction of the Kravchuk functions for at most 256 points (j ≤ 128)"

using HypergeometricFunctions
using SpecialFunctions

export kp, binom, Φ, Φ2

"`kp(n, q, N)` gives the symmetric Kravchuk polynomia for n ∈ [0, N] and q ∈ [-N,N]"
kp(n, q, N) = HypergeometricFunctions._₂F₁general2(-BigFloat(n), -q, -N, 2)

A(n, j) = ((-1.0)^n) / (2.0^j)

"`binom(n, k)` gives the binomial coefficient using the Beta function definition"
binom(n, k) = (1/n) * (1/beta(k + 1, n - k + 1))

B(n, q, j) = √(binom(2 * j, n) * binom(2 * j, j + q))

"Fundamental definition of the Kravchuk function"
ϕ(n, q, j) = Float64(A(n, j) * B(n, q, j) * kp(n, j + q, 2 * j))

"`Φ(n, q, j)` gives the Kravchuk function within the ranges n ∈ [0, 2j], q ∈ [-j, j]. It is tested to give correct values for j ≤ 128."
function Φ(n, q, j)
    if 0 ≤ n ≤ j
        return ϕ(n, q, j)
    elseif j < n ≤ 2 * j
        return (-1.0)^q * ϕ(2 * j - n, q, j)
    end
end

"`Φ2(n₁, n₂, q₁, q₂, j)` evaluates the two-dimensional Kravchuk functions over the set of lineal independent eigenmodes"
Φ2(n₁, n₂, q₁, q₂, j) = Φ(n₁, q₁, j) * Φ(n₂, q₂, j)


text:
"""
ptcsolsc_test

Make sure ptcsolsc finds the stable steady state sqrt(2)/2

f(u) = u^3 - lambda u = 0 with lambda = 1/2

The answer(s!!) are u=0, -sqrt(2)/2, and sqrt(2)/2.

The initial iterate is u0=.1, so the correct answer is sqrt(2)/2.

I'm testing for correctness and a match of the iterations statistics

"""
function ptcsolsc_test()
    u0 = 0.1
    ustable = 0.5 * sqrt(2.0)
    uunstable = 0.0
    lambda = 0.5
    #
    # Convergence to the right solution
    #
    ptcdata1 = ptcsolsc(sptest, u0, sptestp; delta0 = 1.0, rtol = 1.e-12)
    ptcdata2 = ptcsolsc(
        spitchfork,
        u0;
        delta0 = 1.0,
        rtol = 1.e-12,
        pdata = lambda,
        keepsolhist = false,
    )
    ptcdatasec = secant(spitchfork, u0; rtol = 1.e-12, pdata = lambda, keepsolhist = false)
    dh = ptcdata1.history - ptcdata2.history
    ndh = norm(dh[:, 1], Inf)
    fdok = (ndh < 1.e-7)
    ptcerr = ptcdata1.solhist .- ustable
    ptcfun = ptcdata1.history
    secok = abs(ptcdatasec.solution) < 1.e-10
    solok = (abs(ptcdata1.solution - ustable) < 1.e-10)
    funok = (abs(ptcdata1.functionval) < 1.e-12)
    histok = (length(ptcfun) == 18)
    ptcdataf = ptcsolsc(sptest, u0; delta0 = 0.1, rtol = 1.e-12)
    errcode = ptcdataf.errcode
    failok = ~ptcdataf.idid && (errcode == 10)
    ptcok = fdok && solok && funok && histok && failok && secok
    if ~ptcok
        println("Failure in Scalar PTC")
        println(ptcdata1)
    end
    return ptcok
end


text:
# Run collapsed Jain-Neal sampler on a gene expression data set.
module GeneExpression

using DelimitedFiles
using Statistics
using BayesianMixtures
B = BayesianMixtures

# Load data
collection = "deSouto"
data_ID = "armstrong-2002-v1"
# collection = "Siri"
# data_ID = "Blood1"
if collection=="deSouto"
    dataset = readdlm("datasets/gene-deSouto/$(data_ID)_database.txt",'\t')
    x = convert(Array{Float64,2}, dataset[3:end,2:end])
    x = log2.(x)
    x = [x[:,i]::Array{Float64,1} for i = 1:size(x,2)]
elseif collection=="Siri"
    dataset = readdlm("datasets/gene-Siri/$(data_ID).csv",',')
    x = convert(Array{Float64,2}, dataset[2:end,2:end])
    x = [x[:,i]::Array{Float64,1} for i = 1:size(x,2)]
else
    error("Unknown collection: $collection")
end
# Normalize to zero mean, unit variance
mu = mean(x)
v = mean([xi.*xi for xi in x]) - mu.*mu  # sample variance
x = [((xi-mu)./sqrt.(v))::Array{Float64,1} for xi in x]

# Run sampler
n_total = 1000  # number of MCMC sweeps to run the algorithm
options = B.options("MVNaaC","MFM",x,n_total; n_keep=1000,t_max=20)
result = B.run_sampler(options)

# Posteriors on t and k
B.open_figure(1)
B.plot_t_posterior(result; color="b", marker="s", label="p(t|data)")
B.plot_k_posterior(result; color="g", marker="o", label="p(k|data)")
B.labels("Posteriors on t and k","t (clusters) or k (components)","")
B.PyPlot.legend(loc="upper right",numpoints=1)
B.PyPlot.xlim(0,15)

# Posterior similarity matrix (probability that i and j are in same cluster)
B.open_figure(2; figure_size=(5,4))
B.plot_similarity_matrix(result)

end


text:
include("v0.2.jl")
include("v0.4.jl")


text:
using HashFunctions
using Test

@testset "HashFunctions.jl" begin
    include("fnv.jl")
end


text:
"""
        initGigaSOM(train, xdim, ydim = xdim;  norm = :none, toroidal = false)

Initialises a SOM.

# Arguments:
- `train`: training data
- `xdim, ydim`: geometry of the SOM
           If DataFrame, the column names will be used as attribute names.
           Codebook vectors will be sampled from the training data.
- `norm`: optional normalisation; one of :`minmax, :zscore or :none`
- `toroidal`: optional flag; if true, the SOM is toroidal.
"""
function initGigaSOM( train, xdim, ydim = xdim;
             norm::Symbol = :none, toroidal = false)

    if typeof(train) == DataFrame
        colNames = [String(x) for x in names(train)]
    else
        colNames = ["x$i" for i in 1:size(train,2)]
    end

    train = convertTrainingData(train)

    numCodes = xdim * ydim

    # normalise training data:
    train, normParams = normTrainData(train, norm)

    # initialise the codes with random samples
    codes = train[rand(1:size(train,1), numCodes),:]
    grid = gridRectangular(xdim, ydim)

    normParams = convert(DataFrame, normParams)
    names!(normParams, Symbol.(colNames))

    # create X,Y-indices for neurons:
    #
    x = y = collect(1:numCodes)
    indices = DataFrame(X = x, Y = y)

    # make SOM object:
    som = Som(codes = codes, colNames = colNames,
           normParams = normParams, norm = norm,
           xdim = xdim, ydim = ydim,
           numCodes = numCodes,
           grid = grid, indices = indices,
           toroidal = toroidal,
           population = zeros(Int, numCodes))
    return som
end

"""
    trainGigaSOM(som::Som, train::DataFrame, kernelFun = gaussianKernel,
                        r = 0.0, epochs = 10)

# Arguments:
- `som`: object of type Som with an initialised som
- `train`: training data
- `kernel::function`: optional distance kernel; one of (`bubbleKernel, gaussianKernel`)
            default is `gaussianKernel`
- `rStart`: optional training radius.
       If r is not specified, it defaults to (xdim+ydim)/3
- `rFinal`: target radius at the last epoch, defaults to 0.1
- `knnTreeFun`: Constructor of the KNN-tree (e.g. from NearestNeighbors package)
- `metric`: Passed as metric argument to the KNN-tree constructor
- `radiusFun`: Function that generates radius decay, e.g. `linearRadius` or `expRadius(10.0)`
- `epochs`: number of SOM training iterations (default 10)

Training data must be convertable to Array{Float34,2} with `convert()`.
Training samples are row-wise; one sample per row. An alternative kernel function
can be provided to modify the distance-dependent training. The function must fit
to the signature fun(x, r) where x is an arbitrary distance and r is a parameter
controlling the function and the return value is between 0.0 and 1.0.
"""
function trainGigaSOM(som::Som, train::DataFrame;
                      kernelFun::Function = gaussianKernel,
                      metric = Euclidean(),
                      knnTreeFun = BruteTree,
                      rStart = 0.0, rFinal=0.1, radiusFun=linearRadius,
                      epochs = 10)

    train = convertTrainingData(train)

    # set default radius:

    if rStart == 0.0
        rStart = √(som.xdim^2 + som.ydim^2) / 2
        @info "The radius has been determined automatically."
    end

    dm = distMatrix(som.grid, som.toroidal)

    codes = som.codes

    nWorkers = nprocs()
    dTrain = distribute(train)

    for j in 1:epochs

        globalSumNumerator = zeros(Float64, size(codes))
        globalSumDenominator = zeros(Float64, size(codes)[1])

        tree = knnTreeFun(Array{Float64,2}(transpose(codes)), metric)

        if nworkers() > 1

            R = Vector{Any}(undef,nworkers())

            @sync begin
                for (idx, pid) in enumerate(workers())
                    @async begin
                        R[idx] =  fetch(@spawnat pid begin doEpoch(localpart(dTrain), codes, tree) end)
                        globalSumNumerator += R[idx][1]
                        globalSumDenominator += R[idx][2]
                    end
                end
            end
        else
            # only batch mode
            sumNumerator, sumDenominator = doEpoch(localpart(dTrain), codes, tree)
            globalSumNumerator += sumNumerator
            globalSumDenominator += sumDenominator
        end

        r = radiusFun(rStart, rFinal, j, epochs)
        println("Radius: $r")
        if r <= 0
            error("Sanity check: radius must be positive")
        end

        wEpoch = kernelFun(dm, r)
        codes = (wEpoch*globalSumNumerator) ./ (wEpoch*globalSumDenominator)
    end

    som.codes[:,:] = codes[:,:]

    return som
end

"""
    doEpoch(x::Array{Float64}, codes::Array{Float64}, tree)

vectors and the adjustment in radius after each epoch.

# Arguments:
- `x`: training Data
- `codes`: Codebook
- `tree`: knn-compatible tree built upon the codes
"""
function doEpoch(x::Array{Float64, 2}, codes::Array{Float64, 2}, tree)

     # initialise numerator and denominator with 0's
     sumNumerator = zeros(Float64, size(codes))
     sumDenominator = zeros(Float64, size(codes)[1])

     # for each sample in dataset / trainingsset
     for s in 1:size(x, 1)

         (bmuIdx, bmuDist) = knn(tree, x[s, :], 1)

         target = bmuIdx[1]

         sumNumerator[target, :] .+= x[s, :]
         sumDenominator[target] += 1
     end

     return sumNumerator, sumDenominator
end

"""
    mapToGigaSOM(som::Som, data)

Return a DataFrame with X-, Y-indices and index of winner neuron for
every row in data.

# Arguments
- `som`: a trained SOM
- `data`: Array or DataFrame with training data.
- `knnTreeFun`: Constructor of the KNN-tree (e.g. from NearestNeighbors package)
- `metric`: Passed as metric argument to the KNN-tree constructor

Data must have the same number of dimensions as the training dataset
and will be normalised with the same parameters.
"""
function mapToGigaSOM(som::Som, data::DataFrame;
                      knnTreeFun = BruteTree,
                      metric = Euclidean())

    data::Array{Float64,2} = convertTrainingData(data)
    if size(data,2) != size(som.codes,2)
        println("    data: $(size(data,2)), codes: $(size(som.codes,2))")
        error(SOM_ERRORS[:ERR_COL_NUM])
    end

    nWorkers = nprocs()
    vis = Int64[]
    tree = knnTreeFun(Array{Float64,2}(transpose(som.codes)), metric)

    if nWorkers > 1
        # distribution across workers
        dData = distribute(data)

        R = Array{Future}(undef,nWorkers, 1)
        @sync for (p, pid) in enumerate(workers())
            @async R[p] = @spawnat pid begin
                # knn() returns a tuple of 2 arrays of arrays (one with indexes
                # that we take out, the second with distances that we discard
                # here). vcat() nicely squashes the arrays-in-arrays into a
                # single vector.
                vcat(knn(tree, transpose(localpart(dData)), 1)[1]...)
            end
        end

        @sync begin
            for (p, pid) in enumerate(sort!(workers()))
                append!(vis, fetch(R[p]))
            end
        end
    else
        vis = vcat(knn(tree, transpose(data), 1)[1]...)
    end

    return DataFrame(index = vis)
end

"""
    scaleEpochTime(iteration::Int64, epochs::Int64)

Convert iteration ID and epoch number to relative time in training.
"""
function scaleEpochTime(iteration::Int64, epochs::Int64)
    # prevent division by zero on 1-epoch training
    if epochs>1
        epochs -= 1
    end

    return Float64(iteration-1) / Float64(epochs)
end

"""
    linearRadius(initRadius::Float64, iteration::Int64, decay::String, epochs::Int64)

Return a neighbourhood radius. Use as the `radiusFun` parameter for `trainGigaSOM`.

# Arguments
- `initRadius`: Initial Radius
- `finalRadius`: Final Radius
- `iteration`: Training iteration
- `epochs`: Total number of epochs
"""
function linearRadius(initRadius::Float64, finalRadius::Float64,
                      iteration::Int64, epochs::Int64)

    scaledTime = scaleEpochTime(iteration,epochs)
    return initRadius*(1-scaledTime) + finalRadius*scaledTime
end

"""
    expRadius(steepness::Float64)

Return a function to be used as a `radiusFun` of `trainGigaSOM`, which causes
exponencial decay with the selected steepness.

Use: `trainGigaSOM(..., radiusFun = expRadius(0.5))`

# Arguments
- `steepness`: Steepness of exponential descent. Good values range
  from -100.0 (almost linear) to 100.0 (really quick decay).

"""
function expRadius(steepness::Float64 = 1.0)
    return (initRadius::Float64, finalRadius::Float64,
            iteration::Int64, epochs::Int64) -> begin

        scaledTime = scaleEpochTime(iteration,epochs)

        if steepness < -100.0
            # prevent floating point underflows
            error("Sanity check: steepness too low, use linearRadius instead.")
        end

        # steepness is simulated by moving both points closer to zero
        adjust = finalRadius * (1 - 1.1^(-steepness))

        if initRadius <= 0 || (initRadius-adjust) <= 0 || finalRadius <= 0
            error("Radii must be positive. (Possible alternative cause: steepness is too high.)")
        end

        initRadius -= adjust
        finalRadius -= adjust

        return adjust + initRadius * ((finalRadius/initRadius)^scaledTime)
    end
end


text:
# ---
# title: 801. Minimum Swaps To Make Sequences Increasing
# id: problem801
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/minimum-swaps-to-make-sequences-increasing/description/>
# hidden: true
# ---
#
# We have two integer sequences `A` and `B` of the same non-zero length.
#
# We are allowed to swap elements `A[i]` and `B[i]`.  Note that both elements
# are in the same index position in their respective sequences.
#
# At the end of some number of swaps, `A` and `B` are both strictly increasing.
# (A sequence is _strictly increasing_ if and only if `A[0] < A[1] < A[2] < ...
# < A[A.length - 1]`.)
#
# Given A and B, return the minimum number of swaps to make both sequences
# strictly increasing.  It is guaranteed that the given input always makes it
# possible.
#
#
#
#     **Example:**
#     Input: A = [1,3,5,4], B = [1,2,3,7]
#     Output: 1
#     Explanation:
#     Swap A[3] and B[3].  Then the sequences are:
#     A = [1, 3, 5, 7] and B = [1, 2, 3, 4]
#     which are both strictly increasing.
#
#
# **Note:**
#
#   * `A, B` are arrays with the same length, and that length will be in the range `[1, 1000]`.
#   * `A[i], B[i]` are integer values in the range `[0, 2000]`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module CVRPSEP

include("main.jl")

export LIB_CVRPSEP, ConstraintType,
        CnstrRecord, CnstrPointer, CnstrPointerList,
        CnstrMgrRecord, CnstrMgrPointer,
        CutRecord, CutManager,
        convert_CMP, init_CMP,
        rounded_capacity_cuts!, strengthened_comb_inequalities!,
        framed_capacity_inequalities!, hypotour_inequalities!,
        multistar_inequalities!

end


text:
using Cubature
using BasicInterpolators: CubicSplineInterpolator, NoBoundaries
using Base.Threads: @threads
using Test

using ClearSky
using ClearSky: 𝓇stream, dIdω, P2ω

const R = 8.31446262

##

totalopticaldepth(σ, g, μ, Pₛ, Pₜ=0) = dτdP(σ, g, μ)*(Pₛ - Pₜ)

function analyticaloutgoing(σ, g, μ, cₚ, Pₛ, Tₛ)
    #whole atmosphere optical depth, to a pressure of 1 Pa
    τ∞ = totalopticaldepth(σ, g, μ, Pₛ)
    #adiabat exponent
    γ = R/(μ*cₚ)
    #blackbody power
    P = stefanboltzmann(Tₛ)
    #analytical OLR solution (Principles of Planetary Climate equation 4.32)
    P*(exp(-τ∞) + τ∞^(-4γ)*hquadrature(τ′->exp(-τ′)*τ′^(4γ), 0, τ∞)[1])
end

function numericaloutgoing(σ, g, μ, cₚ, Pₛ, Tₛ, Pₜ=1e-3, tol=1e-9)
    #nice big, wide wavenumber sample
    ν = [logrange(1e-6, 1e5, 10000, 4); 1e6]
    #absorber wrapper
    U = UnifiedAbsorber(GrayGas(σ, ν))
    #dry adiabat, whole atmosphere
    Γ = DryAdiabat(Tₛ, Pₛ, cₚ, μ)
    # 1/cos(θ)
    m = 1.0
    #mean molar mass
    fμ(T,P) = μ
    #transformed coords
    ω₁, ω₂ = P2ω(Pₛ), P2ω(Pₜ)
    #compute OLR
    olr = zeros(length(ν))
    @threads for i ∈ eachindex(ν)
        I₀ = planck(ν[i], Tₛ)
        olr[i] = π*stream(dIdω, I₀, ω₁, ω₂, U, i, g, m, Γ, fμ, tol)
    end
    #make an interpolator for smooth integration
    ϕ = CubicSplineInterpolator([0.0; ν], [0.0; olr], NoBoundaries())
    #integrate it
    hquadrature(ν->ϕ(ν), 0.0, maximum(ν))[1]
end

##

#define atmospheric parameters
g = 10 #gravity [m/s^2]
μ = 0.01 #molar mass [kg/molecule]
cₚ = 1e3 #heat capacity [J/kg/K]
Pₛ = 1e5 #surface pressure [Pa]
Tₛ = 300 #surface temperature [K]
σ = 10 .^ range(-29, -23, length=10) #gray gas absorption coefs [cm^2/molecule]

#test accuracy with different optical depths
abserr = similar(σ)
relerr = similar(σ)
τ = similar(σ)
for (i,σᵢ) ∈ enumerate(σ)
    τ[i] = totalopticaldepth(σᵢ, g, μ, Pₛ)
    OLRₐ = analyticaloutgoing(σᵢ, g, μ, cₚ, Pₛ, Tₛ)
    OLRₙ = numericaloutgoing(σᵢ, g, μ, cₚ, Pₛ, Tₛ, 1e-6)
    err = OLRₙ - OLRₐ
    abserr[i] = abs(err)
    relerr[i] = abserr[i] / OLRₐ
    @test relerr[i] < 0.01
end

##

"""using Plots

p = plot(τ, relerr,
    xaxis=:log,
    yaxis=:log,
    xlabel="Whole Atmosphere Optical Depth",
    ylabel="OLR [W/m^2] Relative Error",
    legend=false)"""


text:
#= Copyright (c) 2017 Dave Crist

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
=#

function uniquetol_dict(itr,tol=1e-6)
    ut_vals_dict = Dict()
    ut_idxs_dict = Dict()
    if length(itr) == 0
        return ut_vals_dict, ut_idxs_dict
    end

    _itr = copy(itr)
    while(length(_itr) > 0)
        idxs = find(abs.(itr - _itr[1]) .< tol)
        utval = mean(itr[idxs])
        ut_vals_dict[ utval ] = itr[idxs]
        ut_idxs_dict[ utval ] =     idxs
        deleteat!(_itr,find(abs.(_itr - _itr[1]) .< tol))
    end

    return ut_vals_dict, ut_idxs_dict
end

function uniquetol_list(itr::Array{T,1},tol=1e-6) where {T<:Real}
    preal = sortperm(real(itr))     # buffer for permutation idxs of reals
    dreal = diff(real(itr[preal]))  # buffer for differences of reals

    out = [ ];
    pri1=1;
    for pri2 = 1:length(preal)
        # look for tol-large gap in reals
        if pri2==length(preal) || dreal[pri2] > tol
            push!(out, Array{T,1}(itr[preal[pri1:pri2]]))
            pri1=pri2+1;
        end # real tol check
    end # pri2, tol-large gap in reals
    return out
end

function uniquetol_list(itr::Array{T,1},tol=1e-6) where {T<:Complex}
    preal = sortperm(real(itr))     # buffer for permutation idxs of reals
    dreal = diff(real(itr[preal]))  # buffer for differences of reals

    pimag = Array{Int}(size(preal))     # buffer for permutation idxs of imags
    dimag = Array{Float64}(size(dreal)) # buffer for differences of imags

    out = [ ];
    pri1=1;
    for pri2 = 1:length(preal)
        # look for tol-wide gap in reals
        if pri2==length(preal) || dreal[pri2] > tol
            pridxs = preal[pri1:pri2]
            sortperm!( pimag[1:length(pridxs)], imag(itr[pridxs]) )
            dimag[1:(length(pridxs)-1)] = diff(imag(itr[preal[pridxs]]))

            pii1=1;
            for pii2 = 1:length(pridxs)
                if pii2==length(pridxs) || dimag[pii2] > tol
                    push!(out, Array{T,1}(itr[pridxs[pii1:pii2]]))
                    pii1=pii2+1;
                end # imag tol check
            end # pii2, tol-large gap in imags

            pri1=pri2+1;
        end # real tol check
    end # pri2, tol-large gap in reals
    return out
end

# residue.jl
using Polynomials
function residue(pn::Poly, pd::Poly, root_tol=1e-6)
    # calculate degree and roots of denominator
    degn = degree(pn)
    degd = degree(pd)

    # Check if numerator degree >= denominator degree
    # If degn >= degd, divide out polynomial GCD, and continue with remainder
    if degn >= degd
        pq, pr = divrem(pn, pd)
    else
        pq=Poly(0); pr=copy(pn)
    end

    # get roots, and find roots that are numerically close
    # rootsd = roots of denominator polynomial
    rootsd = roots(pd)
    # rd_ut_idxs = groups of idxs in to rootsd,
    #     that represent numerically close roots
    rd_ut_vals, rd_ut_idxs = uniquetol_dict(rootsd,1e-6)

    # Loop over denominator roots, in groups, determined by uniquetol_map
    # uniquetol_map is used to find roots that are numerically close but uneqal
    # abc_polys are the polynomials multiplying A, B, C, etc.
    # in the form
    # (Numerator-Poly) = A*abc_poly[1] + B*abc_poly[2] + C*abc_poly[3] + ...
    abc_polys = []

    for rdidxs = values(rd_ut_idxs)
        # @show rdidxs
        # Account for all other-valued roots
        other_roots = rootsd[setdiff( 1:length(rootsd), rdidxs)]
        abc_poly = prod( [Poly([-rL,1]) for rL in other_roots] )
        # print("rD=$(rootsd[rdidxs[1]]) , (normal): $abc_poly\n")
        push!(abc_polys, abc_poly)

        # Account for repeated roots numerically close to to rD
        for repidx = 2:length(rdidxs) # start at 2 to ignore non-repeated roots
            rD = rootsd[rdidxs[repidx]]
            abc_poly *= Poly([-rD,1])
            push!(abc_polys, abc_poly)
        end
    end

    # With abc_polys, set up a nice set of linear equations for
    # coefficients of NumeratorCoef_i * x^pow = (A+B+C+...) * x^pow

    A = [ p[pidx] for pidx=0:(degree(pd)-1), p in abc_polys ]
    b = coeffs(pr)

    return A \ b, rootsd

end # residue function


text:
struct Numerics
    T::Float64
    dt::Float64
    S::Float64
    ds::Float64
end

struct Simulation
    equation::String #("GPE", "NPSE")
    potential::String #("barrier", "ellipse")
end

struct Potential
    type::String
    width::Float64
    position::Float64
    energy::Float64
    ϵ::Float64
    a::Float64
end

struct Apparatus
    m::Float64 # kg ω
    as::Float64 # m
    ω_perp::Float64 # rad/s
    N::Int64
    γ::Float64
    ω_z::Float64
end

struct InitialState
    type::String
    width::Float64
    v0::Float64
    x0::Float64
end

struct Coefficients
    α::ComplexF64
    β::Function
    γ::Function
    initial::Function
end

function potential(sim::Simulation, app::Apparatus, pot::Potential, s::Float64)
    hbar = 6.62607015e-34 / (2 * pi)
    l_perp = sqrt(hbar / (app.m * app.ω_perp))
    l_z = sqrt(hbar / (app.m * app.ω_z))
    if pot.type == "barrier"
        value = pot.energy * exp(-(s - pot.position)^2 / (2 * (pot.width))^2)
    elseif pot.type == "ellipse"
        target_error = 0.001 * pot.a
        lower_ϕ = -pi
        upper_ϕ = 3 * pi
        let error
            let midpoint
                error = target_error + 1
                while error > target_error
                    midpoint = (lower_ϕ + upper_ϕ) / 2
                    if s - pot.a * E(midpoint, pot.ϵ) > 0
                        lower_ϕ = midpoint
                    else
                        upper_ϕ = midpoint
                    end
                    error = abs(s - pot.a * E(midpoint, pot.ϵ))
                end
                value = -hbar^2 * (1 / pot.a * (sqrt(1 - pot.ϵ^2)) / (sin(midpoint)^2 + sqrt(1 - pot.ϵ^2) * cos(midpoint)^2)^(3 / 2))^2 / (8 * app.m)
            end
        end
    end
    return value
end

function β(sim::Simulation, app::Apparatus, pot::Potential, s::Float64)
    hbar = 6.62607015e-34 / (2 * pi)
    l_perp = sqrt(hbar / (app.m * app.ω_perp))
    l_z = sqrt(hbar / (app.m * app.ω_z))
    if (sim.equation == "NPSE")
        value = im * hbar / (8 * app.m)
    elseif (sim.equation == "GPE")
        value = -im * app.ω_perp
    end
    value += (-im / hbar) * potential(sim::Simulation, app::Apparatus, pot::Potential, s::Float64)
    return value
end

function γ(sim::Simulation, app::Apparatus, state::InitialState, ψ::ComplexF64)
    hbar = 6.62607015e-34 / (2 * pi)
    l_perp = sqrt(hbar / (app.m * app.ω_perp))
    l_z = sqrt(hbar / (app.m * app.ω_z))
    npse_gamma = app.γ * hbar * app.ω_perp / l_z^6 * 1e-4
    σ(ψ::ComplexF64) = l_perp * sqrt(sqrt(1 + 2 * app.as * (app.N - 1) * abs.(ψ)^2))
    if sim.equation == "NPSE"
        if 1 + 2 * app.as * (app.N - 1) * abs.(ψ)^2 >= 0
            value = -im * hbar / (2 * app.m * σ(ψ)^2) - im * app.m * app.ω_perp^2 / (2 * hbar) * σ(ψ)^2 - im * 2 * hbar * app.as * (app.N - 1) / (app.m * σ(ψ)^2) * abs.(ψ)^2 - npse_gamma * app.N^2 * abs.(ψ)^4
        else
            #print("\nCollapse detected!")
            return NaN
        end
    elseif sim.equation == "GPE"
        value = -im * 2 * hbar * app.as * (app.N - 1) / (app.m * l_perp^2) * abs.(ψ)^2
    end
    return value
end

function wave(sim::Simulation, app::Apparatus, state::InitialState, s::Float64)
    hbar = 6.62607015e-34 / (2 * pi)
    l_perp = sqrt(hbar / (app.m * app.ω_perp))
    interaction_g = abs(2*hbar^2 * app.as / app.m / l_perp^2)
    ggg = 2*hbar*app.ω_perp*abs(app.as)
    GSEnergy = -app.N^3/24*ggg^2 * app.m/hbar^2
    Energy = GSEnergy
    k0 = app.m/hbar * state.v0

    if state.type == "gaussian" # Gaussian Pulse
        return sqrt(1 / (sqrt(2 * pi) * state.width)) * exp.(-(s-state.x0) .^ 2 / (4 * state.width^2)) * exp(im * (s-state.x0) * k0)
    elseif state.type == "sech"
        return sqrt(1 / (2 * state.width)) * 2 ./ (exp.(-((s-state.x0) / (state.width))) .+ exp.((s-state.x0) / (state.width))) * exp(im * (s-state.x0) * k0)
    end
end

function get_coefficients(sim::Simulation, app::Apparatus, pot::Potential, state::InitialState)
    # implement with PhysicalConstants.CODATA2018.h
    hbar = 6.62607015e-34 / (2 * pi)
    l_perp = sqrt(hbar / (app.m * app.ω_perp))
    l_z = sqrt(hbar / (app.m * app.ω_z))
    npse_gamma = app.γ * hbar * app.ω_perp / l_z^6 * 1e-4
    @assert(sim.equation in ["NPSE", "GPE"])
    @assert(pot.type in ["barrier", "ellipse", ""])

    α = -im * hbar / (2 * app.m)
    beta(s::Float64) = β(sim::Simulation, app::Apparatus, pot::Potential, s)
    gamma(ψ::ComplexF64) = γ(sim::Simulation, app::Apparatus, state::InitialState, ψ)
    initial_state(s::Float64) = wave(sim::Simulation, app::Apparatus, state::InitialState, s)
    return Coefficients(α, beta, gamma, initial_state)
end

function run_ground_state(num::Numerics, sim::Simulation, app::Apparatus, pot::Potential, state::InitialState)

    coeffs = get_coefficients(sim, app, pot, state)

    display("Running ground-state simulation")
    results = ground_state_solve(num, coeffs)
    plot_ground_state(results...)
end

function run_dynamics(num::Numerics, sim::Simulation, app::Apparatus, pot::Potential, state::InitialState)

    coeffs = get_coefficients(sim, app, pot, state)

    display("Running Dynamic simulation")
    results = ssfm_solve(num, coeffs)
    plot_dynamics(results...)
end

max_stars_repo_name: lorenzifrancesco/soliton-BEC

max_issues_repo_name: lorenzifrancesco/soliton-BEC

max_forks_repo_name: lorenzifrancesco/soliton-BEC


text:
module GridDensities

using Distributions

export GridDensity, pdf

struct GridDensity
    lo::Vector{Float64}
    hi::Vector{Float64}
    bins::Vector{Int64}
    dist::Distributions.Categorical
    volume::Float64
    cart
    linear
end

function GridDensity(data, lo, hi, bins)
    n = prod(bins)
    @assert length(data) == n
    @assert length(lo) == length(hi) == length(bins)
    @assert all(hi .≥ lo)
    volume = prod(hi .- lo)/n # volume for each cell
    p = copy(data) ./ sum(data)
    dist = Categorical(p)
    cart = CartesianIndices(Tuple(bins))
    linear = LinearIndices(Tuple(bins))
    return GridDensity(lo, hi, bins, dist, volume, cart, linear)
end

function _bin(x, lo, hi, bins)
    if x < lo || x > hi
        return 0
    elseif x < hi
        return Int64(floor(bins*(x - lo)/(hi - lo)) + 1)
    else
        return bins
    end
end

function pdf(d::GridDensity, x)
    dims = length(d.bins)
    @assert length(x) == dims
    @assert all(d.lo .≤ x .≤ d.hi)
    inds = [_bin(x[i], d.lo[i], d.hi[i], d.bins[i]) for i = 1:dims]
    ind = d.linear[inds...]
    return Distributions.pdf(d.dist, ind) / d.volume
end

function _rand(bin::Int64, lo, hi, bins)
    return lo + (bin-1+Base.rand())*(hi-lo)/bins
end

function Base.rand(d::GridDensity)
    ind = rand(d.dist)
    dims = length(d.bins)
    coords = d.cart[ind]
    return [_rand(coords[i], d.lo[i], d.hi[i], d.bins[i]) for i = 1:dims]
end

function Base.rand(d::GridDensity, n::Int64)
    return [Base.rand(d) for i = 1:n]
end

end # module


text:
import Distributed

@inline function Distributed.manage(manager::ExistingProcessManager,
                                    id::Integer,
                                    config::Distributed.WorkerConfig,
                                    op::Symbol)
    return nothing
end


text:
@testset "Serialized models" begin
    m = test_simpleLP()

    sm = serialize_model(m, tmpfile("toy1.serialized"))
    sm2 = serialize_model(sm, tmpfile("toy2.serialized"))

    @test typeof(sm) == Serialized{CoreModel} # expected type
    @test typeof(sm2) == Serialized{CoreModel} # no multi-layer serialization

    precache!(sm)

    @test isequal(m, sm.m) # the data is kept okay
    @test sm2.m == nothing # nothing is cached here
    @test isequal(m, COBREXA.Serialization.deserialize(tmpfile("toy2.serialized"))) # it was written as-is
    @test issetequal(
        reactions(convert(StandardModel, sm)),
        reactions(convert(StandardModel, sm2)),
    )
    sm.m = nothing
    @test issetequal(
        metabolites(convert(CoreModelCoupled, sm)),
        metabolites(convert(CoreModelCoupled, sm2)),
    )
end


text:
"""
    Multivector3{T <: Real} <: AbstractMultivector{T}

An immutable 8-tuple of real numbers that represents a member
of a 3-dimensional multivector algebra.

Each `Multivector3` has the form
```math
    p+qC
```
where ``p`` and ``q`` are 2-dimensional multivectors, or
```math
    a+bA+cB+dAB+fC+gAC+hBC+j(AB)C
```
where ``a``, ``b``, ``c``, ``d``, ``f``, ``g``, ``h``, and ``j`` are real
(and of the same type), and ``A * A = 0``, ``B * B = 0``, `` C * C = 0``;
``AB = A * B = -B * A``, ``AC = A * C = -C * A``, ``BC = B * C = -C * B``;
and ``(AB)C = (A * B) * C``.
Here ``*`` is the wedge product.
"""
struct Multivector3{T <: Real} <: AbstractMultivector{T}
    l::Multivector2{T}
    r::Multivector2{T}

    Multivector3{U}(l::Multivector2{U}, r::Multivector2{U}) where {U <: Real} = new(l, r)
end

function Multivector3(x::Multivector2{T}, y::Multivector2{T}) where {T <: Real}
    Multivector3{T}(x, y)
end

function Multivector3(a::T, b::T, c::T, d::T, f::T, g::T, h::T, j::T) where {T <: Real}
    Multivector3{T}(Multivector2(a, b, c, d), Multivector2(f, g, h, j))
end

function Multivector3(a::T, b::T, c::T, d::T, f::T, g::T, h::T) where {T <: Real}
    Multivector3{T}(Multivector2(a, b, c, d), Multivector2(f, g, h))
end

function Multivector3(a::T, b::T, c::T, d::T, f::T, g::T) where {T <: Real}
    Multivector3{T}(Multivector2(a, b, c, d), Multivector2(f, g))
end

function Multivector3(a::T, b::T, c::T, d::T, f::T) where {T <: Real}
    Multivector3{T}(Multivector2(a, b, c, d), Multivector2(f))
end

function Multivector3(a::T, b::T, c::T, d::T) where {T <: Real}
    Multivector3{T}(Multivector2(a, b, c, d), zero(Multivector2{T}))
end

function Multivector3(a::T, b::T, c::T) where {T <: Real}
    Multivector3{T}(Multivector2(a, b, c), zero(Multivector2{T}))
end

function Multivector3(a::T, b::T) where {T <: Real}
    Multivector3{T}(Multivector2(a, b), zero(Multivector2{T}))
end

function Multivector3(a::T) where {T <: Real}
    Multivector3{T}(Multivector2(a), zero(Multivector2{T}))
end

function Multivector3(a::Real, b::Real, c::Real, d::Real, f::Real, g::Real, h::Real, j::Real)
    Multivector3(promote(a, b, c, d, f, g, h, j)...)
end

function Multivector3(a::Real, b::Real, c::Real, d::Real, f::Real, g::Real, h::Real)
    Multivector3(promote(a, b, c, d, f, g, h)...)
end

function Multivector3(a::Real, b::Real, c::Real, d::Real, f::Real, g::Real)
    Multivector3(promote(a, b, c, d, f, g)...)
end

function Multivector3(a::Real, b::Real, c::Real, d::Real, f::Real)
    Multivector3(promote(a, b, c, d, f)...)
end

function Multivector3(a::Real, b::Real, c::Real, d::Real)
    Multivector3(promote(a, b, c, d)...)
end

function Multivector3(a::Real, b::Real, c::Real)
    Multivector3(promote(a, b, c)...)
end

function Multivector3(a::Real, b::Real)
    Multivector3(promote(a, b)...)
end

function Multivector3(x::Multivector2, y::Multivector2)
    Multivector3(promote(asarray(x)..., asarray(y)...)...)
end

function Multivector3(x::Multivector2)
    Multivector3(x, zero(x))
end

function Multivector3(r::Multivector1, s::Multivector1, t::Multivector1, u::Multivector1)
    Multivector3(Multivector2(r, s), Multivector2(t, u))
end

function Multivector3(r::Multivector1, s::Multivector1, t::Multivector1)
    Multivector3(Multivector2(r, s), Multivector2(t))
end

function Multivector3(r::Multivector1, s::Multivector1)
    Multivector3(Multivector2(r, s), Multivector2(zero(s)))
end

function Multivector3(r::Multivector1)
    Multivector3(Multivector2(r))
end

function show(io::IO, z::Multivector3)
    print(io, "[1: ")
    print(io, z.l.l.l)
    print(io, ", A: ")
    print(io, z.l.l.r)
    print(io, ", B: ")
    print(io, z.l.r.l)
    print(io, ", AB: ")
    print(io, z.l.r.r)
    print(io, ", C: ")
    print(io, z.r.l.l)
    print(io, ", AC: ")
    print(io, z.r.l.r)
    print(io, ", BC: ")
    print(io, z.r.r.l)
    print(io, ", (AB)C: ")
    print(io, z.r.r.r)
    print(io, "]")
end

function real(z::Multivector3)
    real(z.l)
end

function unreal(z::Multivector3)
    vcat(unreal(z.l), asarray(z.r))
end

function zero(z::Multivector3{T}) where {T <: Real}
    Multivector3{T}(zero(z.l), zero(z.r))
end

function zero(::Type{Multivector3{T}}) where {T <: Real}
    Multivector3{T}(zero(T))
end

function one(z::Multivector3{T}) where {T <: Real}
    Multivector3{T}(one(z.l))
end

function one(::Type{Multivector3{T}}) where {T <: Real}
    Multivector3{T}(one(T))
end

"""
    conj{T <: Real}(z::Multivector3{T})

The `Multivector3` conjugate. If ``z=a+bA+cB+dAB+fC+gAC+hBC+j(AB)C``,
then `conj(z)` gives
```math
    a-bA-cB-dAB-fC-gAC-hBC-j(AB)C
```
This operation is an involution.
"""
function conj(z::Multivector3{T}) where {T <: Real}
    Multivector3{T}(conj(z.l), -z.r)
end

"""
    cloak{T <: Real}(z::Multivector3{T})

The cloak conjugate changes the sign of even blades. If ``z=a+bA+cB+dAB+fC+gAC+hBC+j(AB)C``,
then `cloak(z)` gives
```math
    -a+bA+cB-dAB+fC-gAC-hBC+j(AB)C
```
This operation is equivalent to `-dagger(z)` and thus is also an involution.
"""
function cloak(z::Multivector3{T}) where {T <: Real}
    Multivector3{T}(cloak(z.l), dagger(z.r))
end

"""
    dagger{T <: Real}(z::Multivector3{T})

The dagger conjugate changes the sign of odd blades. If ``z=a+bA+cB+dAB+fC+gAC+hBC+j(AB)C``, then `dagger(z)` gives
```math
    a-bA-cB+dAB-fC+gAC+hBC-j(AB)C
```
This operation is an involution.
"""
function dagger(z::Multivector3{T}) where {T <: Real}
    Multivector3{T}(dagger(z.l), cloak(z.r))
end

"""
    star{T <: Real}(z::Multivector3{T})

Returns the Hodge star conjugate. If ``z=a+bA+cB+dAB+fC+gAC+hBC+j(AB)C``, then `star(z)` gives
```math
    j+hA-gB+fAB+dC-cAC+bBC+a(AB)C
```
This operation is an involution.
"""
function star(z::Multivector3{T}) where {T <: Real}
    Multivector3{T}(dagger(star(z.r)), star(z.l))
end

"""
    selfstar(z::Multivector3)

The self-star-conjugate part.
In odd number of dimensions, this operation is idempotent.
"""
function selfstar(z::Multivector3)
    (z + star(z)) / 2
end

"""
    antiselfstar(z::Multivector3)

The anti-self-star-conjugate part.
In odd number of dimensions, this operation is idempotent.
"""
function antiselfstar(z::Multivector3)
    (z - star(z)) / 2
end

function (+)(x::Multivector3, y::Multivector3)
    Multivector3(x.l + y.l, x.r + y.r)
end

function (+)(x::Multivector3, y::Multivector2)
    Multivector3(x.l + y, x.r)
end

function (+)(x::Multivector2, y::Multivector3)
    Multivector3(x + y.l, y.r)
end

function (+)(x::Multivector3, y::Multivector1)
    Multivector3(x.l + y, x.r)
end

function (+)(x::Multivector1, y::Multivector3)
    Multivector3(x + y.l, y.r)
end

function (+)(z::Multivector3, a::Real)
    Multivector3(z.l + a, z.r)
end

function (+)(a::Real, z::Multivector3)
    Multivector3(z.l + a, z.r)
end

function (-)(z::Multivector3{T}) where {T <: Real}
    Multivector3{T}(-z.l, -z.r)
end

function (-)(x::Multivector3, y::Multivector3)
    Multivector3(x.l - y.l, x.r - y.r)
end

function (-)(x::Multivector3, y::Multivector2)
    Multivector3(x.l - y, x.r)
end

function (-)(x::Multivector2, y::Multivector3)
    Multivector3(x - y.l, -y.r)
end

function (-)(x::Multivector3, y::Multivector1)
    Multivector3(x.l - y, x.r)
end

function (-)(x::Multivector1, y::Multivector3)
    Multivector3(x - y.l, -y.r)
end

function (-)(z::Multivector3, a::Real)
    Multivector3(z.l - a, z.r)
end

function (-)(a::Real, z::Multivector3)
    Multivector3(a - z.l, -z.r)
end

"""
    (*)(x::Multivector3, y::Multivector3)

Wedge product of two 3-dimensional multivectors.
This operation is non-commutative and non-associative.
"""
function (*)(x::Multivector3, y::Multivector3)
    Multivector3(x.l * y.l,  (y.r * x.l) + (x.r * conj(y.l)))
end

function (*)(x::Multivector3, y::Multivector2)
    Multivector3(x.l * y, x.r * conj(y))
end

function (*)(x::Multivector2, y::Multivector3)
    Multivector3(x * y.l,  y.r * x)
end

function (*)(x::Multivector3, y::Multivector1)
    Multivector3(x.l * y, x.r * conj(y))
end

function (*)(x::Multivector1, y::Multivector3)
    Multivector3(x * y.l,  y.r * x)
end

"""
    (*)(z::Multivector3, a::Real)
    (*)(a::Real, z::Multivector3)

Scaling and/or reflection of a `Multivector3` by a real number.
"""
function (*)(z::Multivector3, a::Real)
    Multivector3(z.l * a, z.r * a)
end

function (*)(a::Real, z::Multivector3)
    Multivector3(a * z.l, a * z.r)
end

function (/)(x::Multivector3, y::Multivector3)
    if iszerodivisor(y)
        error(ZeroDivisorDenominator)
    end

    x * inv(y)
end

function (\)(y::Multivector3, x::Multivector3)
    if iszerodivisor(y)
        error(ZeroDivisorDenominator)
    end

    inv(y) * x
end

function (/)(z::Multivector3, a::Real)
    Multivector3(z.l / a, z.r / a)
end

function (\)(a::Real, z::Multivector3)
    Multivector3(a \ z.l, a \ z.r)
end

function random(::Type{Multivector3{T}}) where T <: Real
    Multivector3{T}(random(Multivector2{T}), random(Multivector2{T}))
end


text:
using MixedModelsSim
using Test

include("utilities.jl")
include("simdat.jl")
# there's no point in testing the big power stuff
# until all the stuff it depends on has been tested
include("power.jl")


text:
# test the entropy stable interface calculation in parallel
"""
  Test the entropy stable integrals in parallel
"""
function test_ESS_parallel()
  @testset "----- testing ESS parallel -----" begin

    if MPI.Comm_rank(MPI.COMM_WORLD) == 1
      rmfile("./entropy.dat")
    end

    fname = "input_vals_ESS_parallel.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)
    data = readdlm("entropy.dat")

    npts, ncols = size(data)

    s_start = data[1, 3]
    s_end = data[end, 3]

    # check total change in entropy function
    @test  abs( (s_start - s_end)/s_start)  < 1e-12

    # check w^T * res at each timestep
    for i=1:npts
      @test  abs(data[i, 4])  < 1e-13
    end

  end

  return nothing
end

#test_ESS_parallel()
add_func1!(EulerTests, test_ESS_parallel, [TAG_LONGTEST])


text:
import Base: +, -
using Dates
using Printf

const minutes_in_day = 24*60
const minutes_in_hour = 60

to_relative_minutes(minutes::Int64) = mod(minutes, minutes_in_day)

to_relative_minutes(hour::Int64, minutes::Int64) = to_relative_minutes(hour*60 + minutes)

struct Clock
    relative_minutes::Int64
    Clock(minutes::Int64) = new(to_relative_minutes(minutes))
    function Clock(hour::Int64, minutes::Int64)
        new(to_relative_minutes(hour, minutes))
    end
end

function (+)(clock::Clock, period::Dates.Period)
    Clock(clock.relative_minutes + convert(Dates.Minute, period).value)
end

function (-)(clock::Clock, period::Dates.Period)
    Clock(clock.relative_minutes - convert(Dates.Minute, period).value)
end

function Base.show(io::IO, clock::Clock)
    minutes = mod(clock.relative_minutes, minutes_in_hour)
    hours = div(clock.relative_minutes, minutes_in_hour)
    timepad(num) = lpad(num, 2, "0")
    print(io, "\"$(timepad(hours)):$(timepad(minutes))\"")
end

max_stars_repo_name: koladilip/julia-exercism

max_issues_repo_name: koladilip/julia-exercism

max_forks_repo_name: koladilip/julia-exercism


text:
create(c::Client, ::Type{Integration}, g::AbstractGuild; kwargs...) = create_guild_integration(c, g.id; kwargs...)

retrieve(c::Client, ::Type{Integration}, g::AbstractGuild) = get_guild_integrations(c, g.id)

update(c::Client, i::Integration, g::AbstractGuild; kwargs...) = modify_guild_integration(c, g.id, i.id; kwargs...)

delete(c::Client, i::Integration, g::AbstractGuild) = delete_guild_integration(c, g.id, i.id)


text:
module Materials

using Images
using FileIO

using ..GfxBase

export Material, Texture
export ShadingModel, Flat, Normal
export PhysicalShadingModel, Lambertian, BlinnPhong

export get_diffuse

## Shading model type hierarchy ##
# Types are as follows:
#
# ShadingModel
#   - Flat - simply the diffuse color of the object
#   - Normal - color-code a pixel according to its surface normal
#   - PhysicalShadingModel
#       - Lambertian - diffuse shading
#       - BlinnPhong - diffuse+specular shading

abstract type ShadingModel end
struct Flat<:ShadingModel end
struct Normal <: ShadingModel end

abstract type PhysicalShadingModel <: ShadingModel end

struct Lambertian <: PhysicalShadingModel end

mutable struct BlinnPhong <: PhysicalShadingModel
    specular_color::RGB{Float32} # color of the highlight
    specular_exp::Float64 # "sharpness" exponent
end

## Texture struct definition ##
mutable struct Texture
    image_data::Array{RGB{Float32},2}
    repeat::Bool
end

""" Texture constructor - loads image data from a file"""
function Texture(image_fn::String, repeat::Bool)
    image_data = load(image_fn)
    Texture(image_data, repeat)
end

## Material struct definition ##
mutable struct Material
    shading_model::ShadingModel
    mirror_coeff::Float64
    texture::Union{Texture,Nothing}
    diffuse_color::Union{RGB{Float32},Nothing}
end

""" Get the diffuse color of a material; if the material is textured,
provide the uv coordinate on the object of that material. """
function get_diffuse(material::Material, uv::Union{Vec2,Nothing})
    if isnothing(material.texture) || isnothing(uv)
        return material.diffuse_color
    else
        return get_texture_value(material.texture, uv)
    end
end

""" Look up a texture value given the uv coordinates """
function get_texture_value(texture::Texture, uv::Vec2)
    # this equation is from the book
    i = convert(Int, round(uv[1]*size(texture.image_data, 1)-.5))
    j = convert(Int, round(size(texture.image_data,2)-uv[2]*size(texture.image_data,2)-.5))
    i =  i%size(texture.image_data,1)
    j =  j%size(texture.image_data,2)

    # make sure we dont index a value that isnt valid
    if i < 1
        i += size(texture.image_data,1)
    end
    if j < 1
        j += size(texture.image_data, 2)
    end
    # get rgb value from uv coord
    rgb = texture.image_data[j, i]
    return rgb
end

end # module Materials


text:
#= TODO: remove
function findfiles(dir::AbstractPath; exts = ("md", "ipynb"), includehidden = false)
    it = filter(collect((relative(p, dir) for p in walkpath(dir)))) do p
        (extension(p) in exts) && (includehidden || !startswith(filename(p), '.'))
    end
    return it
end
=#

mutable struct Project
    sources::Dict{AbstractPath, XTree}
    outputs::Dict{AbstractPath, XTree}
    rewriters::Vector{<:Rewriter}
end

function Project(rewriters)
    sources = merge(Dict{AbstractPath, XTree}(), [createsources!(rewriter) for rewriter in rewriters]...)
    outputs = Dict{AbstractPath, XTree}()
    return Project(sources, outputs, rewriters)
end

Base.show(io::IO, project::Project) = print(io, "Project($(length(project.sources)) documents, $(length(project.rewriters)) rewriters)")

"""
    rewritesources!(project, paths) -> paths

Rewrite source documents named by `paths` as well as new source documents
recursively created by rewriters. Return a list of all rewritten paths.
"""
function rewritesources!(project, paths = Set(keys(project.sources)))
    rewritesources!(project.sources, project.outputs, project.rewriters, paths)
end

function rewritesources!(sources::Dict, outputs::Dict, rewriters::Vector{<:Rewriter}, paths)
    dirtypaths = []
    docs = filter(((k, v),) -> k in paths, sources)

    while !isempty(docs)
        merge!(outputs, rewritedocs(docs, rewriters))
        docs = createsources!(rewriters)
        paths = union(paths, keys(docs))
    end

    return paths
end

"""
    rewritedocs(sources, rewriters) -> outputs

Applies `rewriters` to a collection of `sources`.
"""
function rewritedocs(sources, rewriters)
    outputs = Dict{AbstractPath, XTree}()
    paths = collect(keys(sources))
    Threads.@threads for i in 1:length(paths)
        p = paths[i]
        xtree = sources[p]
        for rewriter in rewriters
            xtree = rewritedoc(rewriter, p, xtree)
        end
        outputs[p] = xtree
    end
    return outputs
end

"""
    createsources!(rewriters) -> sources

Creates new source documents from `rewriters`
"""
function createsources!(rewriters::Vector{<:Rewriter})
    docs = Vector{Dict{AbstractPath, XTree}}(undef, length(rewriters))
    docs = []
    Threads.@threads for i in 1:length(rewriters)
        #docs[i] = createsources!(rewriters[i])
        push!(docs, createsources!(rewriters[i]))
    end
    return merge(docs...)
end
#=

function addfiles(
        sources::Dict,
        outputs::Dict,
        rewriters,
        newsources;
        dirtypaths = Set())
    sources = copy(sources)
    outputs = copy(outputs)
    dirtypaths = addfiles!(sources, outputs, rewriters, newsources; dirtypaths = dirtypaths)
    return sources, outputs, dirtypaths
end

"""
    addfiles(sources, outputs, rewriters, newsources) -> (sources', outputs', dirtypaths)
    addfiles!(sources, outputs, rewriters, newsources) -> dirtypaths

Updates `sources` and `outputs` based on new or updated `changedsources`
using `rewriters`.
"""
function addfiles!(
        sources::Dict,
        outputs::Dict,
        rewriters,
        newsources::Dict;
        dirtypaths = Set())

    isempty(newsources) && return dirtypaths

    # Process new/changed files on document-level
    for (p, xtree) in newsources
        sources[p] = xtree
        for rewriter in rewriters
            xtree = rewritedoc(rewriter, p, xtree)
        end
        outputs[p] = xtree
        push!(dirtypaths, p)
    end

    # Apply project-level changes like updating tree elements,
    # and creating new files.
    newersources = Dict{AbstractPath, XNode}()
    for rewriter in rewriters
        new = createdocs(rewriter)
        newersources = merge(newersources, new)
    end

    # Run recursively with new documents
    return addfiles!(sources, outputs, rewriters, newersources; dirtypaths = dirtypaths)
end

function addfiles!(project::Project, newsources)
    dirtypaths = addfiles!(project.sources, project.outputs, project.rewriters, newsources)
    return dirtypaths
end

=#
function reset!(project::Project)
    foreach(k -> delete!(project.sources, k), keys(project.sources))
    foreach(k -> delete!(project.outputs, k), keys(project.outputs))
    foreach(reset!, project.rewriters)
end


text:
mutable struct VariableFidelitySurrogate{X,Y,L,U,N,F,E} <: AbstractSurrogate
    x::X
    y::Y
    lb::L
    ub::U
    num_high_fidel::N
    low_fid_surr::F
    eps_surr::E
end

function VariableFidelitySurrogate(x,y,lb,ub;
                                   num_high_fidel = Int(floor(length(x)/2)),
                                   low_fid_structure = RadialBasisStructure(radial_function = linearRadial, scale_factor=1.0, sparse = false),
                                   high_fid_structure = RadialBasisStructure(radial_function = cubicRadial ,scale_factor=1.0,sparse=false))

            x_high = x[1:num_high_fidel]
            x_low = x[num_high_fidel+1:end]
            y_high = y[1:num_high_fidel]
            y_low = y[num_high_fidel+1:end]

            #Fit low fidelity surrogate:
            if low_fid_structure[1] == "RadialBasis"
                    #fit and append to local_surr
                low_fid_surr = RadialBasis(x_low,y_low,lb,ub,rad = low_fid_structure.radial_function, scale_factor = low_fid_structure.scale_factor, sparse = low_fid_structure.sparse)

            elseif low_fid_structure[1] == "Kriging"
                low_fid_surr = Kriging(x_low, y_low,lb,ub, p = low_fid_structure.p, theta = low_fid_structure.theta)

            elseif low_fid_structure[1] == "GEK"
                low_fid_surr = GEK(x_low, y_low,lb,ub, p = low_fid_structure.p, theta = low_fid_structure.theta)

            elseif low_fid_structure == "LinearSurrogate"
                low_fid_surr = LinearSurrogate(x_low, y_low,lb,ub)

            elseif low_fid_structure[1] == "InverseDistanceSurrogate"
                low_fid_surr = InverseDistanceSurrogate(x_low, y_low,lb,ub, p = low_fid_structure.p)

            elseif low_fid_structure[1] == "LobacheskySurrogate"
                low_fid_surr = LobacheskySurrogate(x_low, y_low,lb,ub,alpha = low_fid_structure.alpha , n = low_fid_structure.n, sparse = low_fid_structure.sparse)

            elseif low_fid_structure[1] == "NeuralSurrogate"
                low_fid_surr = NeuralSurrogate(x_low, y_low,lb,ub, model = low_fid_structure.model , loss = low_fid_structure.loss ,opt = low_fid_structure.opt ,n_echos = low_fid_structure.n_echos)

            elseif low_fid_structure[1] == "RandomForestSurrogate"
                low_fid_surr = RandomForestSurrogate(x_low, y_low,lb,ub, num_round = low_fid_structure.num_round)

            elseif low_fid_structure == "SecondOrderPolynomialSurrogate"
                low_fid_surr = SecondOrderPolynomialSurrogate(x_low, y_low,lb,ub)

            elseif low_fid_structure[1] == "Wendland"
                low_fid_surr = Wendand(x_low, y_low,lb,ub, eps = low_fid_surr.eps, maxiters = low_fid_surr.maxiters, tol = low_fid_surr.tol)
            else
                throw("A surrogate with the name provided does not exist or is not currently supported with VariableFidelity")
            end

            #Fit surrogate eps on high fidelity data with objective function y_high - low_find_surr
            y_eps = zeros(eltype(y),num_high_fidel)
            @inbounds for i = 1:num_high_fidel
                y_eps[i] = y_high[i]-low_fid_surr(x_high[i])
            end

            if high_fid_structure[1] == "RadialBasis"
                    #fit and append to local_surr
                eps = RadialBasis(x_high,y_eps,lb,ub,rad = high_fid_structure.radial_function, scale_factor = high_fid_structure.scale_factor, sparse = high_fid_structure.sparse)

            elseif high_fid_structure[1] == "Kriging"
                eps = Kriging(x_high, y_eps,lb,ub, p = high_fid_structure.p, theta = high_fid_structure.theta)

            elseif high_fid_structure == "LinearSurrogate"
                eps = LinearSurrogate(x_high, y_eps,lb,ub)

            elseif high_fid_structure[1] == "InverseDistanceSurrogate"
                eps = InverseDistanceSurrogate(x_high, y_eps,lb,ub, high_fid_structure.p)

            elseif high_fid_structure[1] == "LobacheskySurrogate"
                eps = LobacheskySurrogate(x_high, y_eps,lb,ub,alpha = high_fid_structure.alpha , n = high_fid_structure.n, sparse = high_fid_structure.sparse)

            elseif high_fid_structure[1] == "NeuralSurrogate"
                eps = NeuralSurrogate(x_high, y_eps,lb,ub, model = high_fid_structure.model , loss = high_fid_structure.loss ,opt = high_fid_structure.opt ,n_echos = high_fid_structure.n_echos)

            elseif high_fid_structure[1] == "RandomForestSurrogate"
                eps = RandomForestSurrogate(x_high, y_eps,lb,ub, num_round = high_fid_structure.num_round)

            elseif high_fid_structure == "SecondOrderPolynomialSurrogate"
                eps = SecondOrderPolynomialSurrogate(x_high, y_eps,lb,ub)

            elseif high_fid_structure[1] == "Wendland"
                eps = Wendand(x_high, y_eps,lb,ub, eps = high_fid_structure.eps, maxiters = high_fid_structure.maxiters, tol = high_fid_structure.tol)
            else
                throw("A surrogate with the name provided does not exist or is not currently supported with VariableFidelity")
            end
            return VariableFidelitySurrogate(x,y,lb,ub,num_high_fidel,low_fid_surr,eps)
end

#=
function (varfid::VariableFidelitySurrogate)(val::Number)
    return varfid.eps_surr(val) + varfid.low_fid_surr(val)
end

"""
VariableFidelitySurrogate(x,y,lb,ub;
                                   num_high_fidel = Int(floor(length(x)/2))
                                   low_fid = RadialBasisStructure(radial_function = linearRadial, scale_factor=1.0, sparse = false),
                                   high_fid = RadialBasisStructure(radial_function = cubicRadial ,scale_factor=1.0,sparse=false))
First section (1:num_high_fidel) of samples are high fidelity, second section are low fidelity
"""
function VariableFidelitySurrogate(x,y,lb,ub;
                                   num_high_fidel = Int(floor(length(x)/2))
                                   low_fid = RadialBasisStructure(radial_function = linearRadial, scale_factor=1.0, sparse = false),
                                   high_fid = RadialBasisStructure(radial_function = cubicRadial ,scale_factor=1.0,sparse=false))

end
=#

function (varfid::VariableFidelitySurrogate)(val)
    return varfid.eps_surr(val) + varfid.low_fid_surr(val)
end

"""
add_point!(varfid::VariableFidelitySurrogate,x_new,y_new)

I expect to add low fidelity data to the surrogate.

"""
function add_point!(varfid::VariableFidelitySurrogate,x_new,y_new)
    if length(varfid.x[1]) == 1
        #1D
        varfid.x = vcat(varfid.x,x_new)
        varfid.y = vcat(varfid.y,y_new)

        #I added a new lowfidelity datapoint, I need to update the low_fid_surr:
        add_point!(varfid.low_fid_surr,x_new,y_new)
    else
        #ND
        varfid.x = vcat(varfid.x,x_new)
        varfid.y = vcat(varfid.y,y_new)

        #I added a new lowfidelity datapoint, I need to update the low_fid_surr:
        add_point!(varfid.low_fid_surr,x_new,y_new)
    end

end


text:
using Libdl

nautyver = "nauty27rc4"
nautydir = "$(@__DIR__)/$nautyver/"

run(`tar xfz $(@__DIR__)/$nautyver.tar.gz -C $(@__DIR__)`)
cd(()->run(`./configure --enable-tls`), nautydir)
cd(()->run(`make nauty.a CCOBJ='${CC} -fPIC -c ${CFLAGS} -o $@'`), nautydir)

nautyobj = nautydir .* ["traces.o", "gtools.o", "naututil.o", "nautinv.o", "gutil1.o", "gutil2.o", "gtnauty.o", "naugroup.o", "nautycliquer.o", "nauty.o", "nautil.o", "nausparse.o", "naugraph.o", "schreier.o", "naurng.o"]

run(`gcc -shared -fPIC -I $nautydir -o $(@__DIR__)/nauty.$dlext $(@__DIR__)/extra.c $nautyobj`)


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

using Markdown
using InteractiveUtils

# ╔═╡ 828b5fdc-2f57-11ec-2ba2-bd6c394912fa
begin
	using PlutoUI
	using LinearAlgebra
	using PrettyTables
end

# ╔═╡ 8c65af87-71e1-45ed-a635-d97a6e6b7a3e
md"""$(PlutoUI.LocalResource("./figs/Fig-ReactorFilter-F21.png"))

__Fig 1: Reactor-Filter problem schematic.__ Starting material A is converted to product B by enzyme E in a well-mixed reactor with volume V. The reactor is well insulated. The output from the reactor is fed into a filtration unit that has no moving parts and is NOT well insulated."""

# ╔═╡ b56695f3-68c1-4b76-8cc8-bb57d7e34699
md"""
Consider the reaction/separation process shown in Fig 1.
    A starting material $A$ is converted to product $B$ by enzyme $E$ in a well-mixed and well-insulated reactor.
    Mixing is achieved in the reactor using a mechanical stirring device.
    The enzyme $E$ is immobilized in the reactor (does not flow out), and is stable. Downstream of the reactor,
    a filtration device separates unreacted starting material $A$ from product $B$.
    The filtration unit has no moving parts, and is NOT well insulated.

 __Assumptions:__
* the reactor and filtration units operate at steady-state;
* neglect changes in the kinetic and potential energy in the systems and streams;
* the volume of the reactor is $V$ = 14 L;
* the density of the systems and streams is constant $\rho$ = 1.2 kg/L.

"""

# ╔═╡ 59a9a683-e2da-4690-9950-80e2bcedbb18
md"""

__Compute:__
* (a) Compute the missing values for the mol and volumetric flow rates if the specific reaction rate $\hat{r}_{1}$ = 6.0 mmol L$^{-1}$ min$^{-1}$.
* (b) Compute the missing mass flow rate values.
* (c) Compute the rate of work $\dot{W}_{sh}$ (kJ/min) required to keep the reactor well mixed.
* (d) Compute the rate of heat flow $\dot{Q}$ (kJ/min) into (or from) the filtration unit.
"""

# ╔═╡ f198f507-e5ca-4ec1-a2bd-648965f94780
md"""
##### a) Compute the missing mol and volumetric flow rates in the table.

###### Volumetric flow rates:
Let's start with the volumetric flow rates. You can compute the volumetic flow rate from the _total mass balance_ equation. At steady state, for _any process unit_ we know that:

$$\sum_{s}v_{s}\dot{m}_{s} = 0$$
However, we know that $\dot{m}_{s}$ = $\rho_{s}\dot{F}_{s}$, where $\rho_{s}$ denotes the density of stream $s$ (units: kg/L), and $\dot{F}_{s}$ denotes the volumetric flow rate for stream $s$ (units: L/min). Thus, we can re-write the total mass balance in terms of the volumetric flow rates:

$$\sum_{s}v_{s}\rho_{s}\dot{F}_{s} = 0$$
In this problem (from the assumptions) we have assumed that the density of all streams was the same (and constant) i.e., $\rho_{1}=\rho_{2}\dots\rho_{\mathcal{S}}\equiv\rho$ which means we can pull the density out of the sum, and divide by $\rho$ to give:

$$\sum_{s}v_{s}\dot{F}_{s} = 0$$

For the __reactor__ we have a single in and a single out, so $\dot{F}_{1}$ = $\dot{F}_{2}$. However, for the __filter__ we have a single in, but two exit streams, thus: $\dot{F}_{2} = \dot{F}_{3}+\dot{F}_{4}$
"""

# ╔═╡ 3e17b4fe-1025-471b-ba7e-f257d3f9ef0b
md"""
###### Mol flow rates:

To compute the missing mol flow rates, we start with the steady-state species mol balance equation (which is true for any process unit):

$$\sum_{s}v_{s}\dot{n}_{k,s}+\dot{n}_{gen,k} = 0$$

where $\dot{n}_{k,s}$ denotes the flow of mols of species $k$ in stream $s$, and
$\dot{n}_{gen,k}$ denotes the generation (reaction) terms for species $k$ that are occuring inside the process unit (system). Lastly, we known (in the absence of cells) that the generation (reaction) term(s) are given by:

$$\dot{n}_{gen,k} = \left(\sum_{r}\sigma_{kr}\hat{r}_{r}\right)V$$

where $\sigma_{kr}$ denotes the stoichiometric coefficient for species $k$ in reaction $r$, and $\hat{r}_{r}$ denotes the _specific reaction rate_ for reaction $r$ (units: mmol/min-L) and $V$ denotes volume (units: L). For this problem, we have a single reaction whose rate $\hat{r}_{1}$ which was given in the problem, the volume $V$ was given, and $\sigma_{11}$ = -1 and $\sigma_{21}$ = +1, where $A$ = 1 and $B$ = 2.

For the __reactor__, the species mol balances are given by:

$$\begin{eqnarray}
\dot{n}_{11} -\dot{n}_{12} -\hat{r}_{1}V &=& 0\\
\dot{n}_{21} -\dot{n}_{22} +\hat{r}_{1}V &=& 0\\
\end{eqnarray}$$

For the __filtration unit__, the species mol balances are given by:

$$\begin{eqnarray}
\dot{n}_{12} -\dot{n}_{13} - \dot{n}_{14} &=& 0\\
\dot{n}_{22} -\dot{n}_{23} - \dot{n}_{24} &=& 0\\
\end{eqnarray}$$
"""

# ╔═╡ 75e41abc-bf6b-42b6-a400-3e68b49e53a1
begin

	# data given in the problem setup -
	n11 = 100.0 	# units: mol/min
	n21 = 0.0 		# units: mol/min
	n24 = 0.0 		# units: mol/min
	F1_dot = 10.0 	# units: L/min
	F3_dot = 4.0 	# units: L/min

	# reaction terms -> change for γ and δ versions of prelim
	r1_hat = 4.0    # units: mmol/min-L
	V = 14.0 		# units: L
	n14 = 8.0 		# units: mol/min

	# show -
	nothing
end

# ╔═╡ 51cf0922-fa58-4335-9f34-ae033203f4df
# matrix vector around reactor and filtration -
begin

	# setup both units - unknowns: n12, n22, n13, n23, F2_dot and F4_dot -
	AT = [-1 0 0 0 0 0; 0 -1 0 0 0 0; 1 0 -1 0 0 0; 0 1 0 -1 0 0; 0 0 0 0 -1 0; 0 0 0 0 1 -1];
	bT = [-n11 + r1_hat*V ; -n21 - r1_hat*V ; n14 ; n24 ; -F1_dot ; F3_dot]
	x_dot = inv(AT)*bT
	n12 = x_dot[1]
	n22 = x_dot[2]
	n13 = x_dot[3]
	n23 = x_dot[4]
	F2_dot = x_dot[5]
	F4_dot = x_dot[6]

	# show -
	nothing
end

# ╔═╡ 99d53a10-dc5e-43ae-8095-4ef594f9c420
begin
	# put all the data into a table, so we can display the solution using the PrettyTable.j package -
	data_table = zeros(4,4)

	# row 1 -
	data_table[1,1] = 1
	data_table[1,2] = n11
	data_table[1,3] = n21
	data_table[1,4] = F1_dot

	# row 2 -
	data_table[2,1] = 2
	data_table[2,2] = n12
	data_table[2,3] = n22
	data_table[2,4] = F2_dot

	# row 3 -
	data_table[3,1] = 3
	data_table[3,2] = n13
	data_table[3,3] = n23
	data_table[3,4] = F3_dot

	# row 4 -
	data_table[4,1] = 4
	data_table[4,2] = n14
	data_table[4,3] = n24
	data_table[4,4] = F4_dot

	with_terminal() do

		# setup the table -
		header_row = (["Stream","n1,i","n2,i","Vol Flow Rate"],["","[mmol/min]","[mmol/min]","[L/min]"]);

		# write the table -
		pretty_table(data_table; header=header_row)
	end
end

# ╔═╡ ccee2491-9a2e-45d2-a607-f27b54ad66ae
md"""
###### b) Compute the missing mass flow rates in the table.

To compute the missing mass flow rates, we apply a total mass balance to each process unit. At steady-state we know that:

$$\sum_{s}v_{s}\dot{m}_{s} = 0$$

For the __reactor__ we have a single in (given) and a single out, thus: $\dot{m}_{1}=\dot{m}_{2}$. However, for the filter we have a single in $\dot{m}_{2}$ but two output streams $\dot{m}_{3}$ and $\dot{m}_{4}$, thus:

$$\dot{m}_{2} = \dot{m}_{3} + \dot{m}_{4}$$

"""

# ╔═╡ c3967197-6c0c-4f7f-ac02-20a95ba921ca
begin

	# setup problem -
	m1_dot = 12.0 	# units: kg/min
	m3_dot = 4.8 	# units: kg/min

	# reactor -
	m2_dot = m1_dot

	# filter -
	m4_dot = m2_dot - m3_dot

	# show -
	nothing
end

# ╔═╡ 230abf7b-db83-49b9-a0d0-124aa2b328e5
begin

	# initialize -
	mass_data_table = zeros(4,4)

	# row 1
	mass_data_table[1,1] = 1
	mass_data_table[1,2] = 1.2
	mass_data_table[1,3] = m1_dot
	mass_data_table[1,4] = 120.0

	# row 2
	mass_data_table[2,1] = 2
	mass_data_table[2,2] = 1.2
	mass_data_table[2,3] = m2_dot
	mass_data_table[2,4] = 200.0

	# row 3
	mass_data_table[3,1] = 3
	mass_data_table[3,2] = 1.2
	mass_data_table[3,3] = m3_dot
	mass_data_table[3,4] = 125.0

	# row 3
	mass_data_table[4,1] = 4
	mass_data_table[4,2] = 1.2
	mass_data_table[4,3] = m4_dot
	mass_data_table[4,4] = 260.0

	with_terminal() do

		# setup the table -
		mass_header_row = (["Stream","density","mdot_i","H"],["","[kg/L]","[kg/min]","[kJ/kg]"]);

		# write the table -
		pretty_table(mass_data_table; header=mass_header_row)

	end
end

# ╔═╡ 00a50f4a-838e-4f5a-b01d-0a754ca5a256
md"""
###### c) Compute the rate of work $\dot{W}_{sh}$ (kJ/min) required to keep the reactor well mixed.

To compute the rate of work, we start from the steady-state open energy balance equation:

$$\dot{Q}+\dot{W}_{sh}+\sum_{s}v_{s}\dot{m}_{s}H_{s} = 0$$

and throw out (or simplify) terms that don't apply to reactor process unit. First,
for the reactor we have a single input and a single output, thus $\dot{m}_{1}=\dot{m}_{2}$. Next, the reactor is well-insulated, thus, we know that $\dot{Q}$ = 0. Putting these ideas together gives:

$$\dot{W}_{sh} + \dot{m}_{1}\left(H_{1}-H_{2}\right) = 0$$ or (after we solve for work):

$$\dot{W}_{sh} = \dot{m}_{1}\left(H_{2}-H_{1}\right)$$.
"""

# ╔═╡ 1f41d0db-abb1-48c6-8651-655ba4f73686
begin
	# alias the H's
	H₁ = mass_data_table[1,4];
	H₂ = mass_data_table[2,4];
	H₃ = mass_data_table[3,4];
	H₄ = mass_data_table[4,4];

	# compute the ΔH -
	ΔH = (H₂ - H₁)

	# compute the work -
	W = m1_dot*ΔH
end

# ╔═╡ ddaabaaa-f1dd-4853-9412-f7972dc4dc6d
md"""
###### d) Compute the rate of heat flow $\dot{Q}$ (kJ/min) into (or from) the filtration unit.

To compute the rate of heat transfer, we start from the steady-state open energy balance equation:

$$\dot{Q}+\dot{W}_{sh}+\sum_{s}v_{s}\dot{m}_{s}H_{s} = 0$$

and throw out (or simplify) terms that don't apply to filtration process unit.
First, the filtration unit has no moving parts, thus, the rate of shaft work $\dot{W}_{sh}$. However, unlike the reactor, we have a single input but multiple outputs, thus, the energy balance for the filtration unit is given by:

$$\dot{Q} + \dot{m}_{2}H_{2} - \dot{m}_{3}H_{3} - \dot{m}_{4}H_{4} = 0$$

All the enthalpies were given in the problem, and we know the $\dot{m}_{s}$ from part b). To compute $\dot{Q}$:

$$\dot{Q} = \dot{m}_{3}H_{3} + \dot{m}_{4}H_{4} - \dot{m}_{2}H_{2}$$

"""

# ╔═╡ 8196a83d-b705-479a-b8fd-551f1ef95bd3
begin
	# compute the Qdot -
	Q_dot = m3_dot*H₃ + m4_dot*H₄ - m2_dot*H₂
end

# ╔═╡ db72171c-3964-4992-8d00-c9f8217a829d
html"""
<style>
main {
    max-width: 1200px;
    width: 95%;
    margin: auto;
    font-family: "Roboto, monospace";
}

a {
    color: blue;
    text-decoration: none;
}

.H1 {
    padding: 0px 30px;
}
</style>"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"

[compat]
PlutoUI = "~0.7.16"
PrettyTables = "~1.2.2"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[HypertextLiteral]]
git-tree-sha1 = "f6532909bf3d40b308a0f360b6a0e626c0e263a8"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.1"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "98f59ff3639b3d9485a03a72f3ab35bab9465720"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.6"

[[PlutoUI]]
deps = ["Base64", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "4c8a7d080daca18545c56f1cac28710c362478f3"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.16"

[[PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "69fd065725ee69950f3f58eceb6d144ce32d627d"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.2.2"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "fed34d0e71b91734bf0a7e10eb1bb05296ddbcd0"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.6.0"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
"""

# ╔═╡ Cell order:
# ╟─8c65af87-71e1-45ed-a635-d97a6e6b7a3e
# ╟─b56695f3-68c1-4b76-8cc8-bb57d7e34699
# ╟─59a9a683-e2da-4690-9950-80e2bcedbb18
# ╟─f198f507-e5ca-4ec1-a2bd-648965f94780
# ╟─3e17b4fe-1025-471b-ba7e-f257d3f9ef0b
# ╠═75e41abc-bf6b-42b6-a400-3e68b49e53a1
# ╠═99d53a10-dc5e-43ae-8095-4ef594f9c420
# ╟─51cf0922-fa58-4335-9f34-ae033203f4df
# ╟─ccee2491-9a2e-45d2-a607-f27b54ad66ae
# ╠═c3967197-6c0c-4f7f-ac02-20a95ba921ca
# ╠═230abf7b-db83-49b9-a0d0-124aa2b328e5
# ╟─00a50f4a-838e-4f5a-b01d-0a754ca5a256
# ╠═1f41d0db-abb1-48c6-8651-655ba4f73686
# ╟─ddaabaaa-f1dd-4853-9412-f7972dc4dc6d
# ╠═8196a83d-b705-479a-b8fd-551f1ef95bd3
# ╠═828b5fdc-2f57-11ec-2ba2-bd6c394912fa
# ╟─db72171c-3964-4992-8d00-c9f8217a829d
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
@testset "High-level functions" begin
    test_ex(VulkanGen.promote_hl(APIFunction(func_by_name(:vkEnumeratePhysicalDevices), false)), :(
        function enumerate_physical_devices(instance)::ResultTypes.Result{Vector{PhysicalDevice}, VulkanError}
            val = @propagate_errors(_enumerate_physical_devices(instance))
            val
        end
    ))

    test_ex(VulkanGen.promote_hl(APIFunction(func_by_name(:vkGetPerformanceParameterINTEL), false)), :(
        function get_performance_parameter_intel(device, parameter::PerformanceParameterTypeINTEL)::ResultTypes.Result{PerformanceValueINTEL, VulkanError}
            val = @propagate_errors _get_performance_parameter_intel(device, parameter)
            PerformanceValueINTEL(val)
        end
    ))

    test_ex(VulkanGen.promote_hl(APIFunction(func_by_name(:vkDestroyInstance), false)), :(
        function destroy_instance(instance; allocator = C_NULL)
            _destroy_instance(instance; allocator)
        end
    ))

    test_ex(VulkanGen.promote_hl(APIFunction(func_by_name(:vkGetPhysicalDeviceProperties2), false)), :(
        function get_physical_device_properties_2(physical_device, next_types::Type...)
            next_types_hl = next_types
            next_types = intermediate_type.(next_types)
            PhysicalDeviceProperties2(_get_physical_device_properties_2(physical_device, next_types...), next_types_hl...)
        end
    ))

    test_ex(VulkanGen.promote_hl(APIFunction(func_by_name(:vkGetPhysicalDeviceImageFormatProperties2), false)), :(
        function get_physical_device_image_format_properties_2(physical_device, image_format_info::PhysicalDeviceImageFormatInfo2, next_types::Type...)::ResultTypes.Result{ImageFormatProperties2, VulkanError}
            next_types_hl = next_types
            next_types = intermediate_type.(next_types)
            val = @propagate_errors(_get_physical_device_image_format_properties_2(physical_device, convert(_PhysicalDeviceImageFormatInfo2, image_format_info), next_types...))
            ImageFormatProperties2(val, next_types_hl...)
        end
    ))

    test_ex(VulkanGen.promote_hl(APIFunction(create_func(:vkCreateInstance), false)), :(
        function create_instance(enabled_layer_names::AbstractArray, enabled_extension_names::AbstractArray; allocator = C_NULL, next = C_NULL, flags = 0, application_info = C_NULL)::ResultTypes.Result{Instance, VulkanError}
            val = @propagate_errors(_create_instance(enabled_layer_names, enabled_extension_names; allocator, next, flags, application_info))
            val
        end
    ))
end

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
bitstype 8 Element

import Base: box, unbox

convert(::Type{Element}, no::Uint8) = box(Element, unbox(Uint8, no))
convert(::Type{Uint8}, elm::Element) = box(Uint8, unbox(Element, elm))

# List of elements (118 elements)
const elements = [
    # 1
    :H,  :He, :Li, :Be, :B,  :C,  :N,  :O,  :F,  :Ne,
    :Na, :Mg, :Al, :Si, :P,  :S,  :Cl, :Ar, :K,  :Ca,
    :Sc, :Ti, :V,  :Cr, :Mn, :Fe, :Co, :Ni, :Cu, :Zn,
    :Ga, :Ge, :As, :Se, :Br, :Kr, :Rb, :Sr, :Y,  :Zr,
    :Nb, :Mo, :Tc, :Ru, :Rh, :Pd, :Ag, :Cd, :In, :Sn,
    # 51
    :Sb, :Te, :I,  :Xe, :Cs, :Ba, :La, :Ce, :Pr, :Nd,
    :Pm, :Sm, :Eu, :Gd, :Tb, :Dy, :Ho, :Er, :Tm, :Yb,
    :Lu, :Hf, :Ta, :W,  :Re, :Os, :Ir, :Pt, :Au, :Hg,
    :Tl, :Pb, :Bi, :Po, :At, :Rn, :Fr, :Ra, :Ac, :Th,
    :Pa, :U,  :Np, :Pu, :Am, :Cm, :Bk, :Cf, :Es, :Fm,
    # 101
    :Md, :No, :Lr, :Rf, :Db, :Sg, :Bh, :Hs, :Mt, :Ds,
    :Rg, :Cn, :Uut,:Fl, :Uup,:Lv, :Uus,:Uuo
]

for (no, elm) in enumerate(elements)
    @eval begin
        const $(symbol(string("Elm_", elm))) = convert(Element, uint8($no))
    end
end

show(io::IO, elm::Element) = print(io, string(elements[uint8(elm)]))

function parse_element(s::String)
    # TODO: slow?
    sym = symbol(ucfirst(lowercase(strip(s))))
    for (no, elm) in enumerate(elements)
        if elm === sym
            return convert(Element, uint8(no))
        end
    end
    error("invalid element: '$s'")
end


text:
module CreateTableTasks

import SearchLight.Migrations: create_table, column, primary_key, add_index, drop_table

function up()
  create_table(:tasks) do
    [
      primary_key()
      column(:content, :string)
      column(:done, :bool)
    ]
  end

  add_index(:tasks, :content)
end

function down()
  drop_table(:tasks)
end

end

max_stars_repo_name: yshutaro/GenieTodoList

max_issues_repo_name: yshutaro/GenieTodoList

max_forks_repo_name: yshutaro/GenieTodoList


text:
@testset "quadratures" begin
    for T in (Float32, Float64, BigFloat)
        x, w = legendregauss(T, 1)
        @test iszero(x)
        @test w ≈ [2one(T)]

        x, w = legendregauss(T, 1, OneDimensionalNodes.LeftEndPoint())
        @test x ≈ [-one(T)]
        @test w ≈ [2one(T)]

        x, w = legendregauss(T, 1, OneDimensionalNodes.RightEndPoint())
        @test x ≈ [one(T)]
        @test w ≈ [2one(T)]

        x, w = legendregauss(T, 2, OneDimensionalNodes.LeftEndPoint())
        @test x ≈ [-one(T); T(1 // 3)]
        @test w ≈ [T(1 // 2); T(3 // 2)]

        x, w = legendregauss(T, 2, OneDimensionalNodes.RightEndPoint())
        @test x ≈ [T(-1 // 3); one(T)]
        @test w ≈ [T(3 // 2); T(1 // 2)]

        x, w = legendregauss(T, 2)
        @test x ≈ [-1 / sqrt(T(3)); 1 / sqrt(T(3))]
        @test w ≈ ones(T, 2)

        x, w = legendregauss(T, 3)
        @test x ≈ [-sqrt(T(3 // 5)); T(0); sqrt(T(3 // 5))]
        @test w ≈ [T(5 // 9); T(8 // 9); T(5 // 9)]

        err = ErrorException("Must have at least two points for both ends.")
        @test_throws err legendregauss(T, 1, OneDimensionalNodes.BothEndPoint())

        a, b = OneDimensionalNodes.legendrecoefficients(T, 100)
        err = ErrorException("No convergence after 1 iterations " *
                             "(try increasing maxiterations)")
        @test_throws err OneDimensionalNodes.gaussrule(-one(T), one(T), a, b, OneDimensionalNodes.BothEndPoint(), 1)
    end

    x1, w1 = legendregauss(23)
    x2, w2 = FastGaussQuadrature.gausslegendre(23)

    @test x1 ≈ x2
    @test w1 ≈ w2

    x1, w1 = legendregausslobatto(33)
    x2, w2 = FastGaussQuadrature.gausslobatto(33)

    @test x1 ≈ x2
    @test w1 ≈ w2

    x1, w1 = legendregauss(13, OneDimensionalNodes.LeftEndPoint())
    x2, w2 = FastGaussQuadrature.gaussradau(13)

    @test x1 ≈ x2
    @test w1 ≈ w2
end


text:
function EuropeanCallEFD(K, T, S, sig, r, div, N, Nj, dx)

  dt = T/N
  nu = r - div - 0.5 * sig^2
  edx = exp(dx)
  pu = 0.5 * dt * ((sig/dx)^2 + nu/dx)
  pm = 1 - dt * (sig/dx)^2 - r * dt
  pd = 0.5 * dt * ((sig/dx)^2 - nu/dx)

  ## Initialise asset prices at maturity

  St = zeros(Float64, 2N+1)

  ## creating a shifted ranged for addressing the Array in the usual i,j notation

  J = sub(1:(2N+1), (1+(N+1):((2N+1)+(N+1))))

  ## Asset price at node N,-N
  St[J[-Nj]] = S * exp(-Nj * dx)

  ## Asset prices from bottom -N, to N, at time step N
  for j = (-Nj+1):Nj
    St[J[j]] = St[J[j]-1] * edx
  end

  ## initialise Array for the tree

  C = zeros(Float64, (2N+1, N+1))

  ## Option Values at expiry, time step = N, in rows -N to N

  for j = -Nj:Nj
    C[J[j], Nj+1] = max(0, St[J[j]] - K)
  end

  ## Stepping back through the lattice

  for i in range(N-1, -1, N)
    for j = -Nj+1:Nj-1
      C[J[j], i+1] = pu * C[J[j]+1, i+2] + pm * C[J[j], i+2] + pd * C[J[j]-1, i+2]
    end
    C[J[-Nj], i+1] = C[J[-Nj]+1, i+1]
    C[J[Nj], i+1] = C[J[Nj]-1, i+1] + (St[J[Nj]]-St[J[Nj]-1])

  end

C
end


text:
# using Revise

using DistributedFactorGraphs
using IncrementalInference
using Test

@testset "deprecation of old api" begin

LinearRelative(Normal(2.0, 0.1))

end

@testset "testing compare functions for variables and factors..." begin

fg = initfg()

v1 = addVariable!(fg, :x0, ContinuousScalar)
f1 = addFactor!(fg, [:x0;], Prior(Normal()))

@test compareVariable(v1,v1)
@test compareFactor(f1,f1)

v2 = addVariable!(fg, :x1, ContinuousScalar)
f2 = addFactor!(fg, [:x0;:x1], LinearRelative(Normal(2.0, 0.1)))

fg2 = deepcopy(fg)

@test !compareVariable(v1,v2)

# not testing different factors in this way
# @test !compareFactor(f1,f2)

@test compareAllVariables(fg, fg)

@test compareAllVariables(fg, fg2)

@test compareSimilarVariables(fg, fg)
@test compareSimilarVariables(fg, fg2)

@test compareSimilarFactors(fg, fg)
@test compareSimilarFactors(fg, fg2)

@test compareFactorGraphs(fg, fg)
@test compareFactorGraphs(fg, fg2)

# easier error messages
getSolverParams(fg).multiproc = false

tree, smt, hist = solveTree!(fg)

x1a = getVariable(fg, :x0)
x1b = getVariable(fg2, :x0)

@test !compareVariable(x1a, x1b, skipsamples=false)

@test !compareSimilarVariables(fg, fg2, skipsamples=false)
@test !compareSimilarFactors(fg, fg2, skipsamples=false)

@test compareFactorGraphs(fg, fg)
@test !compareFactorGraphs(fg, fg2, skipsamples=false)

ensureAllInitialized!(fg2)

@test compareSimilarVariables(fg, fg2, skipsamples=true, skip=Symbol[:initialized;:inferdim;:ppeDict;:solvedCount])
# fg2 has been solved, so it should fail on the estimate dictionary
@test !compareSimilarVariables(fg, fg2, skipsamples=true, skip=Symbol[:initialized;:inferdim])

tree = buildTreeReset!(fg2)

@test compareSimilarFactors(fg, fg2, skipsamples=true, skipcompute=true)

@test !compareSimilarFactors(fg, fg2, skipsamples=true, skipcompute=false)

@error "Suppressing one specific factor graph compare test post DFG v0.6.0 due to unknown (likely false) compare failure"
# @test compareFactorGraphs(fg, fg2, skipsamples=true, skipcompute=true, skip=[:initialized;:inferdim;:ppeDict; :solvedCount; :fncargvID])

end

@testset "test subgraph functions..." begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal()))

addVariable!(fg, :x1, ContinuousScalar)
addFactor!(fg, [:x0;:x1], LinearRelative(Normal(2.0, 0.1)))

addVariable!(fg, :x2, ContinuousScalar)
addFactor!(fg, [:x1;:x2], LinearRelative(Normal(4.0, 0.1)))

addVariable!(fg, :l1, ContinuousScalar)
addFactor!(fg, [:x1;:l1], LinearRelative(Rayleigh()))

sfg = buildSubgraph(fg, [:x0;:x1], 1) # distance=1 to include factors

#FIXME JT - this doesn't make sense to pass, it is a subgraph so should it not rather be ⊂ [subset]?
# compareDFG(fg1, fg2, by=⊂, skip=...)
@test fg.sessionId == sfg.sessionId[1:length(fg.sessionId)]
@test compareFactorGraphs(fg, sfg, skip=[:labelDict;:addHistory;:logpath;:sessionId])

# drawGraph(sfg)

end

#