text:
using GridInterpolations
using Test
using LinearAlgebra
using Random
using DelimitedFiles

# use import otherwise Interpolations.interpolate conflicts with GridInterpolations.interpolate
import Interpolations: BSpline, Linear, OnGrid, Flat
import Interpolations

function compareToGrid(testType::Symbol=:random, numDims::Int=3, pointsPerDim::Int=5, testPointsDim::Int=5; numRandomTests::Int=1000, eps::Float64=1e-10)
    # Compare interpolation results from GridInterpolations with the Grid Module:
    # create a multi-dimensional test Grid with data, cut points, etc.
    # (dataM, gridM, gridI) = compareToGrid(2, 3, 5);
    cutsDim = pointsPerDim*ones(Int,numDims)
    gridData = rand(cutsDim...)

    # Set up the data structures:
    gridI = Interpolations.interpolate(gridData, BSpline(Linear()))      # This is for the Interpolations package
    gridI = Interpolations.extrapolate(gridI, Flat())
    cutPointMat = 0:(1/(pointsPerDim-1)):1
    cutPoints = [cutPointMat for i=1:numDims]
    gridM = RectangleGrid(tuple(cutPoints...)...)         # Mykel's interpolation package

    # Package the data in a way that can be read by the interpolation package:
    dataM = Array{Float64}(undef, length(gridM))
    for i=1:length(gridM)
        # x = ind2x(gridM, i)
        x = ind2x(gridM, i)
        dataM[i] = gridData[getFractionalIndexes(gridM,x)...]
    end

    # Check a finer grid of coordinates, for comprehensiveness?

    if testType==:random
        # Interpolate in GridInterpolations with the same points, report any discrepancy
        # Select a number of random points and compare the interpolation results.
        testPoint = zeros(numDims)
        for i=1:numRandomTests
            rand!(testPoint)
            testI = gridI(((pointsPerDim.-1)*testPoint.+1)...)
            testM = GridInterpolations.interpolate(gridM,dataM,testPoint)

            if (abs(testI-testM)>eps)
                display("Failed Random Point Interpolation Test")
                display("Non-matching interpolation points:")
                display("Test point = $testPoint")
                display("Grid Package interpolation result = $testI")
                display("Interpolation Module result = $testM")
                return false
            end

        end
        return true
    end

    if testType==:extrapNeg
        # Test behavior outside the grid, i.e. extrapolation:
        testPoint = zeros(numDims)
        for i=1:numRandomTests
            testPoint = -rand!(testPoint)
            testI = gridI(((pointsPerDim-1)*testPoint+1)...)
            testM = interpolate(gridM,dataM,testPoint)

            if (abs(testI-testM)>eps)
                display("Failed Random Negative Point Extrapolation Test")
                display("Non-matching interpolation points:")
                display("Test point = $testPoint")
                display("Grid Package interpolation result = $testI")
                display("Interpolation Module result = $testM")
                return false
            end

        end
        return true
    end

    if testType==:extrapPos
        for i=1:numRandomTests
            testPoint = rand(numDims).+1
            testI = gridI(((pointsPerDim-1)*testPoint.+1)...)
            testM = interpolate(gridM,dataM,testPoint)

            if (abs(testI-testM)>eps)
                display("Failed Random Positive Point Extrapolation Test")
                display("Non-matching interpolation points:")
                display("Test point = $testPoint")
                display("Grid Package interpolation result = $testI")
                display("Interpolation Module result = $testM")
                return false
            end

        end
        return true
    end

    return true

end

function getFractionalIndexes(g::AbstractGrid, s::Array)
    # Returns the fractional index of sprime within the grid-defined discretization.

    fracInd = Array{Int64}(undef, length(g.cutPoints))

    for i=1:length(g.cutPoints)
        gridDisc = g.cutPoints[i]
        fracInd[i] = getFracIndex(gridDisc, s[i])
    end

    return fracInd

end

function getFracIndex(vararray::Array, value::Float64)
    # Searches through the array vararray and returns the fractional index position
    # that value lies in.  It CAPassumes vararray increases monotonically, and returns
    # the first or last element if the value is outside the bounds of the array.

    if value <= vararray[1]
        return 1
    end

    if value >= vararray[end]
        return argmax(vararray)
    end

    i=1
    while value > vararray[i+1]
        i+=1
    end

    return i+(value-vararray[i])÷(vararray[i+1] - vararray[i])

end

function simplexMagic(NDISC::Int=20, NPOINTS::Int=3, checkFileName::AbstractString="", eps::Float64=1e-10)
    if isempty(checkFileName)
        checkFileName = joinpath(@__DIR__, "simplexMagicTest20.txt")
    end

    val = transpose([ [8.,1,6] [3,5,7] [4,9,2] ]) # transposed magic(3) from matlab

    sGrid = SimplexGrid(1.:NPOINTS, 1.:NPOINTS)
    sInterpVal = zeros(NDISC, NDISC)

    for i = 1:NDISC, j = 1:NDISC
        valX = mapPt(i, NPOINTS, NDISC)
        valY = mapPt(j, NPOINTS, NDISC)
        sInterpVal[i, j] = interpolate(sGrid, vcat(val...), [valX, valY])
    end

    sInterpValTest = readdlm(checkFileName)

    testErr = sum(abs, sInterpVal-sInterpValTest)

    if (testErr > eps)
        display("Failed Simplex Comparison Test")
        display("Test error = $testErr")
    end

    return testErr<eps

end

# constructs a new grid from repr output and tests to see if it's the same
function reprConstruct()
    g1 = SimplexGrid([1.0, 2.0, 3.0], [1.0, 2.0, 3.0])
    g2 = eval(Meta.parse(repr(g1)))
    return g1.cutPoints == g2.cutPoints
end

function mapPt(pt::Int, NPOINTS::Int, NDISC::Int)
    return pt * (NPOINTS-1) / NDISC + 1
end

function checkCounters()
    dims = 2
    nPts = 100
    grid = RectangleGrid(1.:10, 1.:10)
    if length(grid) != nPts || dimensions(grid) != dims
        return false
    end
    data = rand(10,10)
    temp = interpolate(grid, data, [-1.5,1.5])
    temp = interpolate(grid, data, [1.,1.])
    grid = SimplexGrid(1.:10, 1.:10)
    if length(grid) != nPts || dimensions(grid) != dims
        return false
    end
    temp = interpolate(grid, data, [-1.5,1.5])
    temp = interpolate(grid, data, [1.,1.])
    temp = interpolate(grid, data, [1.5,1.5])
    #    showcompact(grid)
    #    show(grid)
    return true
end

function testMask()
    x = [1.5, 1.5]
    vals = collect(1.:9)
    mask = falses(9)
    grid = RectangleGrid(1.:3, 1.:3)
    interped = interpolate(grid, vals, x)
    masked   = maskedInterpolate(grid, vals, x, mask)
    @test interped ≈ masked atol = 1e-3
    mask[1] = true # mask the first entry
    truth  = 3.6666
    masked = maskedInterpolate(grid, vals, x, mask)
    @test truth ≈ masked atol = 1e-3
    return true
end

@testset "Comparisons" begin
    @test compareToGrid(:random) == true
    @test compareToGrid(:extrapPos) == true
    #@test compareToGrid(:extrapNeg)==true
end

@testset "Magic" begin
    @test simplexMagic() == true
end

# TODO
@testset "Repr" begin
    @test reprConstruct() == true
end

@testset "Counter" begin
    @test checkCounters() == true
end

@testset "Mask" begin
    @test testMask() == true
end

@testset "Duplicates" begin
    # Check that duplicates are disallowed
    @test_throws ErrorException RectangleGrid([1,1])
    @test_throws ErrorException SimplexGrid([1,1])
end

# Sanity check overall implementations
function test_rect_implemented()
    rgrid = RectangleGrid([2,5],[2,5])

    @test length(rgrid) == 4
    @test GridInterpolations.dimensions(rgrid) == 2

    @test ind2x(rgrid,1) == [2,2]
    @test ind2x(rgrid,2) == [5,2]
    @test ind2x(rgrid,3) == [2,5]
    @test ind2x(rgrid,4) == [5,5]
    x = [0,0]
    ind2x!(rgrid, 4, x)
    @test x == [5,5]

    @test interpolants(rgrid, [1,1]) == ([1],[1.0])
    @test interpolants(rgrid, [2,5]) == ([3],[1.0])

    indices, weights = interpolants(rgrid, [1.5,3])
    @test indices == [1,3]
    @test isapprox(weights, [0.666667, 0.333333], rtol=1e-5)

    @test interpolate(rgrid, [1,2,3,4], [1,1]) == 1.0
    @test maskedInterpolate(rgrid, [1,2,3,4], [1,1], BitArray([false, false, false, false])) == 1.0
    @test isnan(maskedInterpolate(rgrid, [1,2,3,4], [1,1], BitArray([true, false, false, false])))

    @test isapprox(interpolate(rgrid, [1,2,3,4], [1.5,3]), 1.66666666, rtol=1e-5)
    @test maskedInterpolate(rgrid, [1,2,3,4], [1.5,3], BitArray([true, false, false, false])) == 3
    @test maskedInterpolate(rgrid, [1,2,3,4], [1.5,3], BitArray([false, false, true, false])) == 1
end

@testset "Rect Impl" begin
    test_rect_implemented()
end

# Sanity check overall implementations
function test_simplex_implemented()
    grid = SimplexGrid([2,5],[2,5])

    @test length(grid) == 4
    @test GridInterpolations.dimensions(grid) == 2

    @test ind2x(grid,1) == [2,2]
    @test ind2x(grid,2) == [5,2]
    @test ind2x(grid,3) == [2,5]
    @test ind2x(grid,4) == [5,5]
    x = [0,0]
    ind2x!(grid, 4, x)
    @test x == [5,5]

    indices, weights = interpolants(grid, [1,1])
    full_weight_indices = findall(x -> x == 1, weights)
    @test  length(full_weight_indices) == 1
    @test  weights[full_weight_indices[1]] == 1.0

    indices, weights = interpolants(grid, [1.5,3])
    full_weight_indices = findall(x -> isapprox(x, 0.666667, rtol=1e-5), weights)
    @test  length(full_weight_indices) == 1
    @test  isapprox(weights[full_weight_indices[1]], 0.666667, rtol=1e-5) == true
    full_weight_indices = findall(x -> isapprox(x, 0.333333, rtol=1e-5), weights)
    @test  length(full_weight_indices) == 1
    @test  isapprox(weights[full_weight_indices[1]], 0.333333, rtol=1e-5)  == true

    @test interpolate(grid, [1,2,3,4], [1,1]) == 1.0
    @test maskedInterpolate(grid, [1,2,3,4], [1,1], BitArray([false, false, false, false])) == 1.0
    @test isnan(maskedInterpolate(grid, [1,2,3,4], [1,1], BitArray([true, false, false, false])))

    @test isapprox(interpolate(grid, [1,2,3,4], [1.5,3]), 1.66666666, rtol=1e-5)
    @test maskedInterpolate(grid, [1,2,3,4], [1.5,3], BitArray([true, false, false, false])) == 3
    @test maskedInterpolate(grid, [1,2,3,4], [1.5,3], BitArray([false, false, true, false])) == 1
end
@testset "Simplex Impl" begin
    test_simplex_implemented()
end

# check whether rectangle & simplex expect sorted order of cutpoints on dimensions
function test_ordering(grid)
    @test ind2x(grid,1) == [2,18] # 1
    @test ind2x(grid,2) == [5,18] # 2
    @test ind2x(grid,3) == [2,15] # 3
    @test ind2x(grid,4) == [5,15] # 4
    @test ind2x(grid,5) == [2,12] # 5
    @test ind2x(grid,6) == [5,12] # 6

    return interpolate(grid, [1,2,3,4,5,6], [6, 20]) == 2.0
end

# tests the order of vertices returned by vertices()
# by comparing against ind2x for each unrolled index
function test_vertices_ordering(grid)

    grid_verts = @inferred vertices(grid)

    @test length(grid_verts) == length(grid)

    for i = 1 : length(grid)
        @test grid_verts[i] == ind2x(grid,i)
    end

    return true
end
@testset "Ordering" begin
    @test_throws ErrorException test_ordering( RectangleGrid([2,5], [18,15,12]) )
    @test_throws ErrorException test_ordering( SimplexGrid([2,5], [18,15,12]) )

    @test test_vertices_ordering( RectangleGrid([2,5], [12,15,18]) ) == true
    @test test_vertices_ordering( SimplexGrid([1,4,6,10], [8,12,17]) ) == true
end

@testset "NaN" begin
    r = GridInterpolations.RectangleGrid(1:4, 1:4)
    @test_throws DomainError interpolants(r, [2.5, NaN])
end

@testset "Iteration and Indexing" begin
    @testset "2d" begin
        g2 = GridInterpolations.RectangleGrid(1:3,4:5)
        pts = ([1.,4.],[2.,4.],[3.,4.],[1.,5.],[2.,5.],[3.,5.])
        @test length(g2) == length(pts)
        @test all(pts .== [x for x in g2])
        @test iterate(g2,7) === nothing
        @test all(pts .== [g2[ci] for ci in reshape(CartesianIndices((3,2)), length(g2))])
        @test all(pts .== reshape([g2[i,j] for i=1:3, j=1:2], length(g2)))
        @test_throws BoundsError g2[CartesianIndex(0,0)]
        @test_throws BoundsError g2[CartesianIndex(4,2)]
        @test_throws BoundsError g2[CartesianIndex(3,3)]
        @test_throws BoundsError g2[0,0]
        @test_throws BoundsError g2[4,2]
        @test_throws BoundsError g2[3,3]
        @test_throws BoundsError g2[1,1,2]
    end
    @testset "3d" begin
        g3 = GridInterpolations.RectangleGrid(1:3,4:5,6:10)
        pts = reshape([[x,y,z] for x=1.:3., y=4.:5., z=6.:10.], length(g3))
        @test length(g3) == length(pts)
        @test all(pts .== [x for x in g3])
        @test iterate(g3,31) === nothing
        @test all(pts .== [g3[ci] for ci in reshape(CartesianIndices((3,2,5)), length(g3))])
        @test all(pts .== reshape([g3[i,j,k] for i=1:3, j=1:2, k=1:5], length(g3)))
        @test_throws BoundsError g3[CartesianIndex(0,0,0)]
        @test_throws BoundsError g3[CartesianIndex(4,2,5)]
        @test_throws BoundsError g3[CartesianIndex(3,3,5)]
        @test_throws BoundsError g3[CartesianIndex(3,2,6)]
        @test_throws BoundsError g3[0,0,0]
        @test_throws BoundsError g3[4,2,5]
        @test_throws BoundsError g3[3,3,5]
        @test_throws BoundsError g3[3,2,6]
        if VERSION.major > 0  # this test is known to fail for julia v0.7
            @test_throws BoundsError g3[1,1]
        end
        @test_throws BoundsError g3[1,1,1,2]
    end
end

include(joinpath(@__DIR__, "..", "bench", "interpBenchmarks.jl"))
compareBenchmarks(4, 10, 100, quiet=false)

println("All tests complete")


text:
const sx2 = [0 0.5; 0.5 0]
const isy2 = [0 0.5; -0.5 0]
const sz2 = [0.5 0; 0 -0.5]
const sp2 = sx2 + isy2
const sm2 = sx2 - isy2
const id2 = Float64[1 0; 0 1]

const sx3 = Float64[0 1 0; 1 0 1; 0 1 0] ./ sqrt(2)
const isy3 = Float64[0 1 0; -1 0 1; 0 -1 0] ./ sqrt(2)
const sz3 = Float64[1 0 0; 0 0 0; 0 0 -1]
const sp3 = sx3 + isy3
const sm3 = sx3 - isy3
const id3 = Float64[1 0 0; 0 1 0; 0 0 1]

const sdots = (kron(sp3, sm3) + kron(sm3, sp3)) / 2 + kron(sz3, sz3)
const blbq(θ) = reshape(cos(θ) .* sdots + sin(θ) .* sdots^2, 3, 3, 3, 3)

ext: jl


text:
using  Test
using SimplePosets, SimplePosetAlgorithms

P = BooleanLattice(5)
C = max_chain(P)
@test length(C) == 6
@test height(P) == 6

A = max_antichain(P)
@test length(A) == 10
@test width(P) == 10

R = realizer(P,5)
Q = realize_poset(R)
@test P==Q

P = StandardExample(4)
@test dimension(P) == 4


text:
"""
    LiftingLineSection <: AbstractModel

Lifting line section model with state variables ``v_x, v_y, v_z, \\omega_x,
\\omega_y, \\omega_z``, inputs ``F_x', F_y', F_z', M_x, M_y, M_z``, and zero
parameters.  Two-dimensional aerodynamic models may be extended to
three-dimensional models by coupling with this model.  Note that
this model has no rate equations of its own since its state variables are
defined as functions of the 3D structural model's state variables.
"""
struct LiftingLineSection <: AbstractModel end

"""
    LiftingLineSection()

Initialize an object of type [`LiftingLineSection`](@ref)
"""
LiftingLineSection()

# --- Traits --- #

number_of_states(::Type{<:LiftingLineSection}) = 6
number_of_inputs(::Type{<:LiftingLineSection}) = 6
number_of_parameters(::Type{<:LiftingLineSection}) = 0

inplaceness(::Type{LiftingLineSection}) = OutOfPlace()

# --- Internal Methods for Couplings --- #

function liftingline_velocities(vx, vz, ωy)
    u = vx
    v = vz
    ω = ωy
    return u, v, ω
end

function liftingline_accelerations(dvx, dvz, dωy)
    udot = dvx
    vdot = dvz
    ωdot = dωy
    return udot, vdot, ωdot
end


text:
include("toolbox.jl")
include("solver.jl")
include("augmentation.jl")

function warm_start(Xmat::Array{T,2}, Y::Array{T,1}, rho::T, n::Int64, d::Int64; n_partition=0, partition_size=0,max_steps=5, maxiter=50, violTOL=1.0e-3, innerTOL=1.0e-4,outerTOL=1.0e-5,block=0) where T<:AbstractFloat
    Is = Array{Int64,1}(undef,0)
    Js = Array{Int64,1}(undef,0)
    prtt_len::Int64 = 0;
    obj::Float64 = 0.;
    lamb_vec = Array{eltype(Y),1}(undef, 0);
    if partition_size != 0 && n_partition != 0
        @assert partition_size*n_partition == n
        m = partition_size
    elseif partition_size !=0
        m = partition_size
    elseif n_partition != 0
        m = cld(n, n_partition);
    else
        m = cld(n, 5);
    end
    partition = partition_n_by_m(n,m);
    for prtt in partition
        prtt_len = length(prtt);
        _,_,lamb_sub, _, I_sub, J_sub, _ ,obj_sub = active_set(Xmat[prtt,:], Y[prtt], rho; max_steps=max_steps,maxiter=maxiter,violTOL=violTOL,
            innerTOL=innerTOL, outerTOL=outerTOL, verbose=0,LBFGS=true, WS=false, block=block, aug = ActiveSetAugmentation(2,prtt_len,block,violTOL))
        append!(Is, prtt[I_sub[lamb_sub.nzval .!=0]]);
        append!(Js, prtt[J_sub[lamb_sub.nzval .!=0]]);
        append!(lamb_vec, lamb_sub[lamb_sub.!=0]);
        obj += obj_sub;
    end

    Wlen = length(lamb_vec);
    if block == 0
        W = sparse(Js, Is, fill(true, Wlen), n,n);
        lamb = sparse(Js, Is, lamb_vec, n,n);
    else
        W = sparse(Is, Js, fill(true, Wlen), n,n);
        lamb = sparse(Is, Js, lamb_vec, n,n);
    end
    return W, lamb, Wlen, obj, Is, Js;
end

function active_set(Xmat::Array{T,2}, Y::Array{T,1}, rho::T;
        max_steps= 5,maxiter=50,violTOL =1.0e-3, innerTOL=1.0e-4, outerTOL=1.0e-5,
        verbose = 1,random_state = 42,LBFGS=false, LS =true, WS=false, inexact = true,
        block = 0, aug = ActiveSetAugmentation(2,n,block,violTOL)) where {T<:AbstractFloat}
    start = time()
    if random_state != -1
        Random.seed!(random_state);
    end
    n,d = size(Xmat);
    Xflat = mat_to_flat(Xmat,n,d);
    last_obj = 0.;
    flag::Int = 0;
    if ~WS
        W = spzeros(Bool, n,n);
        Wlen::Int64 = 0;
        Wlen_prev::Int64 = 0;
        lamb = spzeros(eltype(Y),n,n);
        phi = copy(Y);
        xi = zeros(eltype(Y),n*d);
        obj = 0.;
        if block ==0
            J,I,_ = findnz(W)
        else
            I,J,_ = findnz(W)
        end
    else
        W, lamb, Wlen, obj, I, J = warm_start(Xmat, Y, rho, n, d; n_partition=5);
        Wlen_prev = Wlen;
        phi = zeros(eltype(Y), n);
        xi = zeros(eltype(Y),n*d);
        get_primal_solution!(phi,xi,Xflat,Y,rho,I,J, Wlen, lamb.nzval, n, d)
    end
    if (~LBFGS)&&LS
        L::Float64 = 0.;
    end
    if verbose >= 1
        println("Optimization Started\n")
    end
    for k in 1:maxiter
        if verbose >= 1
            println("\nIteration:  ",k)
            println("Wlen_prev:  ",Wlen_prev)
        end
        Wlen = aug(phi,xi,Xflat,Y,W,Wlen,lamb,n,d)
        if Wlen > Wlen_prev
            if block ==0
                J,I,_ = findnz(W)
            else
                I,J,_ = findnz(W)
            end
        end
        if Wlen - Wlen_prev <= 0.005*n
            flag += 1
        else
            flag = 0
        end
        Wlen_prev = Wlen;
        if verbose >= 1
            println("     Wlen:  ",Wlen)
            println("     flag:  ",flag)
        end
        if LBFGS == false
            if LS == true
                L,obj = cvxreg_qp_ineq_APG_ls(Xflat, Y, rho, I, J, Wlen, lamb, L, n, d, phi, xi; TOL=innerTOL, max_steps = max_steps)
                println("        L:  ",L)
            else
                obj = cvxreg_qp_ineq_APG(Xflat, Y, rho, I, J, Wlen, lamb, n, d, phi, xi; TOL=innerTOL, max_steps = max_steps)
            end
        else
            obj = cvxreg_qp_ineq_LBFGS(Xflat, Y, rho, I, J, Wlen, lamb, n, d, phi, xi; TOL=innerTOL, max_steps = max_steps)
        end
        if verbose >= 1
            println("Objective:  ",obj)
        end
        if k > 1&& (obj - last_obj > -outerTOL || obj > last_obj * (1+outerTOL)) && flag >= 5
            break
        else
            last_obj = obj;
        end
    end
    runtime = time() - start
    if verbose >= 1
        println("\nConvex Regression finishes in ",runtime," s.")
    end
    return phi,xi,lamb, W, I,J, Wlen,obj
end

function active_set_profiling(Xmat::Array{T,2}, Y::Array{T,1}, rho::T;
        max_steps= 5,maxiter=50,maxtime=10800,violTOL =1.0e-3, innerTOL=1.0e-4, outerTOL=1.0e-5,
        verbose = 1,random_state = 42,LBFGS=false, LS =true, WS=false, inexact = true,
        block = 0, aug = ActiveSetAugmentation(2,n,block,violTOL)) where {T<:AbstractFloat}
    start = time()
    cur_time = @elapsed begin
        if random_state != -1
            Random.seed!(random_state);
        end
        n,d = size(Xmat);
        Xflat =  mat_to_flat(Xmat,n,d);
        last_obj = 0.;
        flag::Int = 0;
        if ~WS
            W = spzeros(Bool, n,n);
            Wlen::Int64 = 0;
            Wlen_prev::Int64 = 0;
            lamb = spzeros(eltype(Y),n,n);
            phi = copy(Y);
            xi = zeros(eltype(Y),n*d);
            obj = 0.;
            if block ==0
                J,I,_ = findnz(W)
            else
                I,J,_ = findnz(W)
            end
        else
            W, lamb, Wlen, obj, I, J = warm_start(Xmat, Y, rho, n, d; n_partition=5);
            Wlen_prev = Wlen;
            phi = zeros(eltype(Y), n);
            xi = zeros(eltype(Y),n*d);
            get_primal_solution!(phi,xi,Xflat,Y,rho,I,J, Wlen, lamb.nzval, n, d)
        end
        if (~LBFGS)&&LS
            L::Float64 = 0.;
        end
    end
    searching_times = Array{Float64}([]) # len = # outer
    all_objs = Array{Float64}([]) # len = # inner + outer
    all_obj_times = Array{Float64}([]) # len = # inner + outer
    Wlen_list = Array{Int64}([]) # len = # inner + 1
    if LBFGS
        pow_times = nothing
        L_list = nothing
        f_evals_list = Array{Int64}([]) # len = # outer
    elseif ~LS
        pow_times = Array{Float64}([]) # len = # outer
        f_evals_list = nothing
        L_list = Array{Float64}([]) # len = # outer
    elseif LS
        pow_times = nothing
        L_list = Array{Float64}([]) # len = # outer
        f_evals_list = Array{Int64}([]) # len = # outer
    end

    optimization_starting_steps = Array{Int64}([]) # len = # outer
    searching_time ::Float64 = 0;
    append!(all_objs, last_obj)
    append!(all_obj_times,cur_time)
    append!(Wlen_list, Wlen)
    solvetime::Float64 = 0;
    if verbose >= 1
        println("Optimization Started\n")
    end
    for k in 1:maxiter
        if k >= 2 && searching_times[end]+cur_time >= maxtime
            println("early stopping in searching due to time limit")
            break
        end
        if verbose >= 1
            println("\nIteration:  ", k)
            println("Wlen_prev:  ",Wlen_prev)
        end
        searching_time = @elapsed Wlen = aug(phi,xi,Xflat,Y,W,Wlen,lamb,n,d)
        cur_time += searching_time
        append!(searching_times,searching_time)
        append!(Wlen_list, Wlen)
        if Wlen > Wlen_prev
            if block ==0
                J,I,_ = findnz(W)
            else
                I,J,_ = findnz(W)
            end
        end
        if Wlen - Wlen_prev <= 0.005*n
            flag += 1
        else
            flag = 0
        end
        Wlen_prev = Wlen;
        if verbose >= 1
            println("     Wlen:  ",Wlen)
            println("     flag:  ",flag)
        end
        append!(optimization_starting_steps, length(all_objs)+1)
        if LBFGS == false
            if LS == true
                solvetime = @elapsed status, L, obj, times, objs, f_evals, Ls, _ = cvxreg_qp_ineq_APG_ls_profiling(Xflat, Y, rho, I, J, Wlen, lamb, L, n, d, phi, xi; TOL=innerTOL, max_steps = max_steps,
                    start_time = cur_time, end_time = maxtime)
                append!(f_evals_list, sum(f_evals));
                append!(L_list, L);
                println("        L:  ",L)
            else
                solvetime = @elapsed status, obj, times, objs, pow_time, L, _ = cvxreg_qp_ineq_APG_profiling(Xflat, Y, rho, I, J, Wlen, lamb, n, d, phi, xi; TOL=innerTOL, max_steps = max_steps,
                    start_time = cur_time, end_time = maxtime)
                append!(pow_times, pow_time)
                append!(L_list, L)
                println("        L:  ",L)
            end
        else
            solvetime = @elapsed obj, times, objs, f_evals, _ = cvxreg_qp_ineq_LBFGS_profiling(Xflat, Y, rho, I, J, Wlen, lamb, n, d, phi, xi; TOL=innerTOL, max_steps = max_steps)
            append!(f_evals_list, sum(f_evals))
            status= 0
        end
        append!(all_objs,objs)
        append!(all_obj_times,cur_time.+times)
        cur_time += solvetime;
        if verbose >= 1
            println("Objective:  ",obj)
            println("     Time:  ",all_obj_times[end])
        end
        if status == -1
            break
        end
        if cur_time >= maxtime
            println("early stopping after optimization due to time limit")
            break
        end
        if k > 1&& (obj - last_obj > -outerTOL || obj > last_obj * (1+outerTOL)) && flag >= 5
            break
        else
            last_obj = obj;
        end
    end
    runtime = time() - start
    if verbose >= 1
        println("\nConvex Regression finishes in ",runtime," s.")
    end
    return phi,xi,lamb, W, Wlen, all_objs,all_obj_times,Wlen_list,searching_times,pow_times,f_evals_list, L_list, optimization_starting_steps
end

struct AlgorithmParameters
    LBFGS::Bool;
    LS::Bool;
    WS::Bool;
    inexact::Bool;
    max_steps::Int64;
    violTOL::Float64;
    innerTOL::Float64;
    outerTOL::Float64;
    function AlgorithmParameters(;LBFGS=false,LS=true,WS=false,inexact=true,max_steps=-1,violTOL=1.0e-3,innerTOL=-1.,outerTOL=1.0e-5)
        if inexact
            max_steps = (max_steps == -1) ? 5 : max_steps
            innerTOL = (innerTOL==-1.0) ? 1e-6 : innerTOL
        else
            max_steps = (max_steps == -1) ? 3000 : max_steps
            innerTOL = (innerTOL==-1.0) ? 1e-7 : innerTOL
        end
        new(LBFGS, LS, WS, inexact, max_steps, violTOL, innerTOL, outerTOL);
    end
end

function two_stage_active_set(Xmat::Array{T,2}, Y::Array{T,1}, rho::T; verbose = 1, random_state = 42, maxiter = 50, block = 0, dropzero = true,
        params_list::Array{AlgorithmParameters,1}, augs::Array{ActiveSetAugmentation,1}) where {T<:AbstractFloat}
    start = time()
    if random_state != -1
        Random.seed!(random_state);
    end
    n,d = size(Xmat);
    Xflat = mat_to_flat(Xmat,n,d);
    last_obj = 0.;
    flag::Int = 0;
    WS = params_list[1].WS
    LS = params_list[1].LS
    LBFGS = params_list[1].LBFGS
    if ~WS
        W = spzeros(Bool, n,n);
        Wlen::Int64 = 0;
        Wlen_prev::Int64 = 0;
        lamb = spzeros(eltype(Y),n,n);
        phi = copy(Y);
        xi = zeros(eltype(Y),n*d);
        obj = 0.;
        if block ==0
            J,I,_ = findnz(W)
        else
            I,J,_ = findnz(W)
        end
    else
        W, lamb, Wlen, obj, I, J = warm_start(Xmat, Y, rho, n, d; n_partition=5);
        Wlen_prev = Wlen;
        phi = zeros(eltype(Y), n);
        xi = zeros(eltype(Y),n*d);
        get_primal_solution!(phi,xi,Xflat,Y,rho,I,J, Wlen, lamb.nzval, n, d)
    end
    if (~LBFGS)&&LS
        L::Float64 = 0.;
    end
    stage::Int64 = 1;
    if verbose >= 1
        println("Optimization Started\n")
    end
    for k in 1:maxiter
        if verbose >= 1
            println("\nIteration:  ", k)
            println("    Stage:  ", stage)
            println("Wlen_prev:  ",Wlen_prev)
        end
        Wlen = augs[stage](phi,xi,Xflat,Y,W,Wlen,lamb,n,d)
        if Wlen > Wlen_prev
            if block ==0
                J,I,_ = findnz(W)
            else
                I,J,_ = findnz(W)
            end
        end
        if Wlen - Wlen_prev <= 0.005*n
            flag += 1
        elseif stage == 1
            flag = 0
        end
        Wlen_prev = Wlen;
        if verbose >= 1
            println("     Wlen:  ",Wlen)
            println("     flag:  ",flag)
        end
        if LBFGS == false
            if LS == true
                L,obj = cvxreg_qp_ineq_APG_ls(Xflat, Y, rho, I, J, Wlen, lamb, L, n, d, phi, xi; TOL=params_list[stage].innerTOL, max_steps = params_list[stage].max_steps)
                println("        L:  ",L)
            else
                obj = cvxreg_qp_ineq_APG(Xflat, Y, rho, I, J, Wlen, lamb, n, d, phi, xi; TOL=params_list[stage].innerTOL, max_steps = params_list[stage].max_steps)
            end
        else
            obj = cvxreg_qp_ineq_LBFGS(Xflat, Y, rho, I, J, Wlen, lamb, n, d, phi, xi; TOL=params_list[stage].innerTOL, max_steps = params_list[stage].max_steps)
        end
        if verbose >= 1
            println("Objective:  ",obj)
        end
        if stage == 1 && k > 1 && flag >= 5
            flag = 0;
            stage = 2;
            if dropzero
                lamb = dropzeros(lamb)
                if block ==0
                    J,I,_ = findnz(lamb)
                else
                    I,J,_ = findnz(lamb)
                end
                Wlen = length(lamb.nzval)
                Wlen_prev = Wlen
                W = sparse(J,I,fill(true,Wlen),n,n)
            end
        elseif stage == 2 && (k > 1&& (obj - last_obj > -params_list[stage].outerTOL || obj > last_obj * (1+params_list[stage].outerTOL)))
            if flag >= 10
                break
            else
                flag += 1
            end
        end
        last_obj = obj;
    end
    runtime = time() - start
    if verbose >= 1
        println("\nConvex Regression finishes in ",runtime," s.")
    end
    return phi,xi,lamb, W, I,J, Wlen,obj
end

function two_stage_active_set_profiling(Xmat::Array{T,2}, Y::Array{T,1}, rho::T; verbose = 1, random_state = 42, maxiter = 50, maxtime = 10800, block = 0, dropzero = true,
        params_list::Array{AlgorithmParameters,1}, augs::Array{ActiveSetAugmentation,1}) where {T<:AbstractFloat}
    start = time()
    cur_time = @elapsed begin
        if random_state != -1
            Random.seed!(random_state);
        end
        n,d = size(Xmat);
        Xflat = mat_to_flat(Xmat,n,d);
        last_obj = 0.;
        flag::Int = 0;
        WS = params_list[1].WS
        LS = params_list[1].LS
        LBFGS = params_list[1].LBFGS
        if ~WS
            W = spzeros(Bool, n,n);
            Wlen::Int64 = 0;
            Wlen_prev::Int64 = 0;
            lamb = spzeros(eltype(Y),n,n);
            phi = copy(Y);
            xi = zeros(eltype(Y),n*d);
            obj = 0.;
            if block ==0
                J,I,_ = findnz(W)
            else
                I,J,_ = findnz(W)
            end
        else
            W, lamb, Wlen, obj, I, J = warm_start(Xmat, Y, rho, n, d; n_partition=5);
            Wlen_prev = Wlen;
            phi = zeros(eltype(Y), n);
            xi = zeros(eltype(Y),n*d);
            get_primal_solution!(phi,xi,Xflat,Y,rho,I,J, Wlen, lamb.nzval, n, d)
        end
        if (~LBFGS)&&LS
            L::Float64 = 0.;
        end
        stage::Int64 = 1;
    end
    stages = Array{Int64}([]) # len = # outer
    searching_times = Array{Float64}([]) # len = # outer
    all_objs = Array{Float64}([]) # len = # inner + outer
    all_obj_times = Array{Float64}([]) # len = # inner + outer
    Wlen_list = Array{Int64}([]) # len = # inner + 1
    if LBFGS
        pow_times = nothing
        L_list = nothing
        f_evals_list = Array{Int64}([]) # len = # outer
    elseif ~LS
        pow_times = Array{Float64}([]) # len = # outer
        f_evals_list = nothing
        L_list = Array{Float64}([]) # len = # outer
    elseif LS
        pow_times = nothing
        L_list = Array{Float64}([]) # len = # outer
        f_evals_list = Array{Int64}([]) # len = # outer
    end

    optimization_starting_steps = Array{Int64}([]) # len = # outer
    searching_time ::Float64 = 0;
    append!(all_objs, last_obj)
    append!(all_obj_times,cur_time)
    append!(Wlen_list, Wlen)
    solvetime::Float64 = 0;
    if verbose >= 1
        println("Optimization Started\n")
    end
    for k in 1:maxiter
        if k >= 2 && searching_times[end]+cur_time >= maxtime
            println("early stopping in searching due to time limit")
            break
        end
        append!(stages, stage)
        if verbose >= 1
            println("\nIteration:  ", k)
            println("    Stage:  ", stage)
            println("Wlen_prev:  ",Wlen_prev)
        end
        searching_time = @elapsed Wlen = augs[stage](phi,xi,Xflat,Y,W,Wlen,lamb,n,d)
        cur_time += searching_time
        append!(searching_times,searching_time)
        append!(Wlen_list, Wlen)
        if Wlen > Wlen_prev
            if block ==0
                J,I,_ = findnz(W)
            else
                I,J,_ = findnz(W)
            end
        end
        if Wlen - Wlen_prev <= 0.005*n
            flag += 1
        elseif stage == 1
            flag = 0
        end
        Wlen_prev = Wlen;
        if verbose >= 1
            println("     Wlen:  ",Wlen)
            println("     flag:  ",flag)
        end
        append!(optimization_starting_steps, length(all_objs)+1)
        if LBFGS == false
            if LS == true
                solvetime = @elapsed status, L, obj, times, objs, f_evals, Ls, _ = cvxreg_qp_ineq_APG_ls_profiling(Xflat, Y, rho, I, J, Wlen, lamb, L, n, d, phi, xi;
                    TOL=params_list[stage].innerTOL, max_steps = params_list[stage].max_steps, start_time = cur_time, end_time = maxtime)
                append!(f_evals_list, sum(f_evals));
                append!(L_list, L);
                println("        L:  ",L)
            else
                solvetime = @elapsed status, obj, times, objs, pow_time, L, _ = cvxreg_qp_ineq_APG_profiling(Xflat, Y, rho, I, J, Wlen, lamb, n, d, phi, xi;
                    TOL=params_list[stage].innerTOL, max_steps = params_list[stage].max_steps,  start_time = cur_time, end_time = maxtime)
                append!(pow_times, pow_time)
                append!(L_list, L)
                println("        L:  ",L)
            end
        else
            solvetime = @elapsed obj, times, objs, f_evals, _ = cvxreg_qp_ineq_LBFGS_profiling(Xflat, Y, rho, I, J, Wlen, lamb, n, d, phi, xi;
                TOL=params_list[stage].innerTOL, max_steps = params_list[stage].max_steps)
            append!(f_evals_list, sum(f_evals))
            status = 0;
        end
        append!(all_objs,objs)
        append!(all_obj_times,cur_time.+times)
        cur_time += solvetime;
        if verbose >= 1
            println("Objective:  ",obj)
            println("     Time:  ",all_obj_times[end])
        end
        if status == -1
            break
        end
        if cur_time >= maxtime
            println("early stopping after optimization due to time limit")
            break
        end
        if stage == 1 && k > 1 && flag >= 5
            flag = 0;
            stage = 2;
            if dropzero
                lamb = dropzeros(lamb)
                if block ==0
                    J,I,_ = findnz(lamb)
                else
                    I,J,_ = findnz(lamb)
                end
                Wlen = length(lamb.nzval)
                Wlen_prev = Wlen
                W = sparse(J,I,fill(true,Wlen),n,n)
            end
        elseif stage == 2 &&  (k > 1&& (obj - last_obj > -params_list[stage].outerTOL || obj > last_obj * (1+params_list[stage].outerTOL)))
            if flag >= 10
                break
            else
                flag += 1
            end
        end
        last_obj = obj;
    end
    runtime = time() - start
    if verbose >= 1
        println("\nConvex Regression finishes in ",runtime," s.")
    end
    return phi,xi,lamb, W, Wlen, all_objs,all_obj_times, stages,Wlen_list,searching_times,pow_times,f_evals_list, L_list, optimization_starting_steps
end


text:
# ### m10.3t.jl

using Pkg, DrWatson

@quickactivate "StatisticalRethinkingTuring"
using Turing
using StatisticalRethinking
Turing.setprogress!(false)

delim=';'
df = CSV.read(sr_datadir("chimpanzees.csv"), DataFrame; delim);

# pulled_left, condition, prosoc_left
@model function ppl10_3(y, x₁, x₂)
    α ~ Normal(0, 10)
    βp ~ Normal(0, 10)
    βpC ~ Normal(0, 10)

    logits = α .+ (βp .+ βpC * x₁) .* x₂
    y .~ BinomialLogit.(1, logits)
end;

m10_3t = ppl10_3(df.pulled_left, df.condition, df.prosoc_left)
nchains = 4; sampler = NUTS(0.65); nsamples=2000
chns10_3t = mapreduce(c -> sample(m10_3t, sampler, nsamples), chainscat, 1:nchains)

# Rethinking result

m_10_03t_result = "
      Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
 a    0.05   0.13      -0.15       0.25  3284    1
 bp   0.62   0.22       0.28       0.98  3032    1
 bpC -0.11   0.26      -0.53       0.29  3184    1
";

# End of m10.03t.jl


text:
module MLJMultivariateStatsInterface

# ===================================================================
# IMPORTS
import MLJModelInterface
import MLJModelInterface: Table, Continuous, Finite, @mlj_model, metadata_pkg,
    metadata_model
import MultivariateStats
import MultivariateStats: SimpleCovariance
import StatsBase: CovarianceEstimator

using Distances
using LinearAlgebra

# ===================================================================
## EXPORTS
# Models are exported automatically by `@mlj_model` macro

# ===================================================================
## Re-EXPORTS
export SimpleCovariance, CovarianceEstimator, SqEuclidean, CosineDist

# ===================================================================
## CONSTANTS
# Define constants for easy referencing of packages
const MMI = MLJModelInterface
const MS = MultivariateStats
const PCAFitResultType = MS.PCA
const KernelPCAFitResultType = MS.KernelPCA
const ICAFitResultType = MS.ICA
const PPCAFitResultType = MS.PPCA
const FactorAnalysisResultType = MS.FactorAnalysis
const default_kernel = (x, y) -> x'y #default kernel used in KernelPCA

# Definitions of model descriptions for use in model doc-strings.
const PCA_DESCR = """
      Principal component analysis. Learns a linear transformation to
    project the data  on a lower dimensional space while preserving most of the initial
    variance.
    """
const KPCA_DESCR = "Kernel principal component analysis."
const ICA_DESCR = "Independent component analysis."
const PPCA_DESCR = "Probabilistic principal component analysis"
const FactorAnalysis_DESCR = "Factor Analysis"
const LDA_DESCR = """
      Multiclass linear discriminant analysis. The algorithm learns a
    projection matrix `P` that projects a feature matrix `Xtrain` onto a lower dimensional
    space of dimension `out_dim` such that the trace of the transformed between-class
    scatter matrix(`Pᵀ*Sb*P`) is maximized relative to the trace of the transformed
    within-class scatter matrix (`Pᵀ*Sw*P`).The projection matrix is scaled such that
    `Pᵀ*Sw*P=I` or `Pᵀ*Σw*P=I`(where `Σw` is the within-class covariance matrix) .
    Predicted class posterior probability for feature matrix `Xtest` are derived by
    applying a softmax transformationto a matrix `Pr`, such that  rowᵢ of `Pr` contains
    computed distances(based on a distance metric) in the transformed space of rowᵢ in
    `Xtest` to the centroid of each class.
    """
const BayesianLDA_DESCR = """
      Bayesian Multiclass linear discriminant analysis. The algorithm
    learns a projection matrix `P` that projects a feature matrix `Xtrain` onto a lower
    dimensional space of dimension `out_dim` such that the trace of the transformed
    between-class scatter matrix(`Pᵀ*Sb*P`) is maximized relative to the trace of the
    transformed within-class scatter matrix (`Pᵀ*Sw*P`). The projection matrix is scaled
    such that `Pᵀ*Sw*P = n` or `Pᵀ*Σw*P=I` (Where `n` is the number of training samples
    and `Σw` is the within-class covariance matrix).
    Predicted class posterior probability distibution are derived by applying Bayes rule
    with a multivariate Gaussian class-conditional distribution.
    """
const SubspaceLDA_DESCR = """
    Multiclass linear discriminant analysis. Suitable for high
    dimensional data (Avoids computing scatter matrices `Sw` ,`Sb`). The algorithm learns a
    projection matrix `P = W*L` that projects a feature matrix `Xtrain` onto a lower
    dimensional space of dimension `nc - 1` such that the trace of the transformed
    between-class scatter matrix(`Pᵀ*Sb*P`) is maximized relative to the trace of the
    transformed within-class scatter matrix (`Pᵀ*Sw*P`). The projection matrix is scaled
    such that `Pᵀ*Sw*P = mult*I` or `Pᵀ*Σw*P=mult/(n-nc)*I` (where `n` is the number of
    training samples, mult` is  one of `n` or `1` depending on whether `Sb` is normalized,
    `Σw` is the within-class covariance matrix, and `nc` is the number of unique classes
    in `y`) and also obeys `Wᵀ*Sb*p = λ*Wᵀ*Sw*p`, for every column `p` in `P`.
    Predicted class posterior probability for feature matrix `Xtest` are derived by
    applying a softmax transformation to a matrix `Pr`, such that  rowᵢ of `Pr` contains
    computed distances(based on a distance metric) in the transformed space of rowᵢ in
    `Xtest` to the centroid of each class.
    """
const BayesianSubspaceLDA_DESCR = """
       Bayesian Multiclass linear discriminant analysis. Suitable for high dimensional data
    (Avoids computing scatter matrices `Sw` ,`Sb`). The algorithm learns a projection
    matrix `P = W*L` (`Sw`), that projects a feature matrix `Xtrain` onto a lower
    dimensional space of dimension `nc-1` such that the trace of the transformed
    between-class scatter matrix(`Pᵀ*Sb*P`) is maximized relative to the trace of the
    transformed within-class scatter matrix (`Pᵀ*Sw*P`). The projection matrix is scaled
    such that `Pᵀ*Sw*P = mult*I` or `Pᵀ*Σw*P=mult/(n-nc)*I` (where `n` is the number of
    training samples, `mult` is  one of `n` or `1` depending on whether `Sb` is normalized,
    `Σw` is the within-class covariance matrix, and `nc` is the number of unique classes in
    `y`) and also obeys `Wᵀ*Sb*p = λ*Wᵀ*Sw*p`, for every column `p` in `P`.
    Posterior class probability distibution are derived by applying Bayes rule with a
    multivariate Gaussian class-conditional distribution
    """
const LinearRegressor_DESCR = """
    Linear Regression. Learns a linear combination of given
    variables to fit the response by minimizing the squared error between.
    """
const MultitargetLinearRegressor_DESCR = """
    Multitarget Linear Regression. Learns linear combinations of given
    variables to fit the responses by minimizing the squared error between.
    """
const RidgeRegressor_DESCR = """
    Ridge regressor with regularization parameter lambda. Learns a
    linear regression with a penalty on the l2 norm of the coefficients.
    """
const MultitargetRidgeRegressor_DESCR = """
    Multitarget Ridge regressor with regularization parameter lambda. Learns a
    Multitarget linear regression with a penalty on the l2 norm of the coefficients.
    """
const PKG = "MLJMultivariateStatsInterface"

# ===================================================================
# Includes
include("models/decomposition_models.jl")
include("models/discriminant_analysis.jl")
include("models/linear_models.jl")
include("utils.jl")

# ===================================================================
# List of all models interfaced
const MODELS = (
    LinearRegressor,
    MultitargetLinearRegressor,
    RidgeRegressor,
    MultitargetRidgeRegressor,
    PCA,
    KernelPCA,
    ICA,
    LDA,
    BayesianLDA,
    SubspaceLDA,
    BayesianSubspaceLDA,
    FactorAnalysis,
    PPCA
)

# ====================================================================
# PKG_METADATA
metadata_pkg.(
    MODELS,
    name = "MultivariateStats",
    uuid = "6f286f6a-111f-5878-ab1e-185364afe411",
    url = "https://github.com/JuliaStats/MultivariateStats.jl",
    license = "MIT",
    julia = true,
    is_wrapper = false
)

end


text:
"""
  `code_similarity(onestr, twostr[, code=fuzzy_soundex])`

  Produce a similarity estimate based on string comparison of the result of
  coding the two input strings.

  `code` can be any of: `soundex`, `phonex`, `phonix`, `[fuzzy_soundex]`, `nysiis`,
  `metaphone` or `match_rating_encode`. Results from using `double_metaphone` may be
  misleading due to the double encoding of that system.

  This comparison function was originally used by the authors of Fuzzy Soundex. A similar
  but non-normalised comparison scheme exists for Soundex, see `editex`, and a specific
  comparison algorithm exists for `match_rating_encode`, see `match_rating`.
"""
function code_similarity(onestr::String, twostr::String, code=fuzzy_soundex)
  #encode strings
  ionestr = code(onestr)
  itwostr = code(twostr)

  return code_similarity_internal(ionestr, itwostr, code)
end

" Internal code similarity function (not for raw strings). "
function code_similarity_internal(ionestr::String , itwostr::String, code=fuzzy_soundex)

  l1 = length(ionestr)
  l2 = length(itwostr)

  common = 0
  for i in 1:min(l1,l2)
    if ionestr[i] == itwostr[i]
      common += 1
    end
  end

  ð = (2*common)/(l1+l2)
  return ð
end

"""
  `code_match(str, array[, code=fuzzy_soundex, permissive=0.0])`

  Returns those items in `array` which are phonetic matches to `str` according
  to the algorithm `code`.

  Matches are normally on exact codes. The `permissive` option allows you to set
  a lower threshold for comparison for most `code` options, using `code_similarity()`.
  The value should be in the range 0:1.  Permissive coding is not possible for
  `double_metaphone`.

  `match_rating_encode` will be treated like other codes rather than how it is used
  in the overall approach. Those wishing to match items according to that system can
  call

      filter(item -> meets_match_rating(str, item), array)

  instead.
"""
function code_match(str::String, array::Array{T,1}, code=fuzzy_soundex::Function, permissive=0.0::Float64) where {T<:String}

  #catch weird input
  if permissive > 1 || permissive < 0
    error("Permissive value must be in range 0:1")
  end

  #If partial matching
  if permissive > 0 && permissive < 1

    #catch attempt at double_metaphone (makes code messy)
    if code == double_metaphone
      error("Cannot currently code match with double_metaphone")
    end

    #returns any better matches than the set level
    return filter(array) do item
      sim_val = code_similarity(item, str, code)
      return sim_val > permissive
    end
  end

  #encode key
  encoded_key = code(str)

  #filter straight matches
  return filter(array) do item

    encoded_val = code(item)

    #handle double_metaphone
    if code == double_metaphone && isa(encoded_val, Array{String,1})
      return (encoded_key in encoded_val)
    else
      return encoded_key == encoded_val
    end
  end

end

"""
  `code_cluster(array[, code=phonix, lower_threshold=0.7, higher_threshold=0.9, stochastic=true])`

  Uses a phonetic similarity measure to cluster an array of strings.

  This function uses a canopy clustering method to group input strings based
  on their phonetic similarity, as judged by `code_similarity`. The phonetic
  coding method to be used can be specified, as can the upper and lower thresholds
  for similarity and whether or not to use randomly sampled centroids for the clusters.

  - Setting `higher_threshold=1` and `lower_threshold=1` will result in only exact
    code matches, which may be particularly sensible for coding schemes with variable
    length codes.

  - Setting `lower_threshold` to the same value as `higher_threshold` will result
    in cleanly separated clusters -- that is, each element of the original array
    will be assigned to one and only one cluster.

  - Setting `stochastic=false` will cause centroid to be selected from left to right
    rather than randomly, meaning that the output clusters will be deterministic. Behaviour
    without this setting is to select a centroid randomly, which can result in some clusters
    appearing or disappearing between runs on the same data where either threshold < 1.
"""
function code_cluster(array::Array{T, 1}, code=phonix::Function, lower_threshold=0.7, higher_threshold=0.9, stochastic=true) where {T<:String}

  #Sanity check
  if lower_threshold > higher_threshold
    error("Lower threshold cannot be greater than higher threshold.")
  end

  clusters = Vector[]

  #local copy of input array to shrink
  coderef = copy(array)

  #encoded version of input array
  inarray = map(code, array)

  while length(inarray) > 0

    if stochastic
      centroid = rand(inarray)
    else
      centroid = inarray[1]
    end

    #similarity of all items to selected centroid
    similarities = map(x -> code_similarity_internal(x, centroid), inarray)

    #select words which are greater than the lower threshold to add to this cluster
    incluster = filter(x -> similarities[findall(y -> y == x, coderef)[1]] >= lower_threshold, coderef)
    push!(clusters, incluster)

    #remove codes and words which are higher than the higher threshold from consideration
    inarray = filter(x -> similarities[findall(y -> y == x, inarray)[1]] < higher_threshold, inarray)
    coderef = filter(x -> similarities[findall(y -> y == x, coderef)[1]] < higher_threshold, coderef)
  end
  return clusters
end

export code_similarity, code_match, code_cluster


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

using JuliaFEM, Test

X = Dict(
    1 => [2.0, 3.0, 4.0],
    2 => [6.0, 3.0, 2.0],
    3 => [2.0, 5.0, 1.0],
    4 => [4.0, 3.0, 6.0])
X[5] = 1/2*(X[1] + X[2])
X[6] = 1/2*(X[2] + X[3])
X[7] = 1/2*(X[3] + X[1])
X[8] = 1/2*(X[1] + X[4])
X[9] = 1/2*(X[2] + X[4])
X[10] = 1/2*(X[3] + X[4])

element = Element(Tet10, (1, 2, 3, 4, 5, 6, 7, 8, 9, 10))
update!(element, "youngs modulus", 480.0)
update!(element, "poissons ratio", 1/3)
update!(element, "geometry", X)
update!(element, "density", 105.0)

problem = Problem(Heat, "tet10", 1)
add_element!(problem, element)
time = 0.0
assemble_mass_matrix!(problem, time)

M_expected = [
    6   1   1   1  -4  -6  -4  -4  -6  -6
    1   6   1   1  -4  -4  -6  -6  -4  -6
    1   1   6   1  -6  -4  -4  -6  -6  -4
    1   1   1   6  -6  -6  -6  -4  -4  -4
   -4  -4  -6  -6  32  16  16  16  16   8
   -6  -4  -4  -6  16  32  16   8  16  16
   -4  -6  -4  -6  16  16  32  16   8  16
   -4  -6  -6  -4  16   8  16  32  16  16
   -6  -4  -6  -4  16  16   8  16  32  16
   -6  -6  -4  -4   8  16  16  16  16  32]

@test isapprox(problem.assembly.M, M_expected)


text:
# ---
# title: 3. Longest Substring Without Repeating Characters
# id: problem3
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Hash Table, Two Pointers, String, Sliding Window
# link: <https://leetcode.com/problems/longest-substring-without-repeating-characters/description/>
# hidden: true
# ---
#
# Given a string `s`, find the length of the **longest substring** without
# repeating characters.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "abcabcbb"
#     Output: 3
#     Explanation: The answer is "abc", with the length of 3.
#
#
# **Example 2:**
#
#
#
#     Input: s = "bbbbb"
#     Output: 1
#     Explanation: The answer is "b", with the length of 1.
#
#
# **Example 3:**
#
#
#
#     Input: s = "pwwkew"
#     Output: 3
#     Explanation: The answer is "wke", with the length of 3.
#     Notice that the answer must be a substring, "pwke" is a subsequence and not a substring.
#
#
# **Example 4:**
#
#
#
#     Input: s = ""
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `0 <= s.length <= 5 * 104`
#   * `s` consists of English letters, digits, symbols and spaces.
#
#
## @lc code=start
using LeetCode

function length_of_longest_substring(s::String)::Int
    seen = Set{Char}()
    res = 0

    for x in s
        if x in seen
            empty!(seen)
        end
        push!(seen, x)
        res = max(res, length(seen))
    end
    return res
end
## @lc code=end


text:
using DataDeps

register(DataDep(
    "IMDB movie reviews dataset",
    """
    author    : Maas, Andrew L.  and  Daly, Raymond E.  and  Pham, Peter T.  and  Huang, Dan  and  Ng, Andrew Y.  and  Potts, Christopher
    title     : Learning Word Vectors for Sentiment Analysis
    booktitle : Proceedings of the 49th Annual Meeting of the Association for Computational Linguistics: Human Language Technologies
    month     : June
    year      : 2011
    address   : Portland, Oregon, USA
    publisher : Association for Computational Linguistics

    This is a dataset for binary sentiment classification.:
    25,000 highly polar reviews for training, 25,000 reviews for testing and additional unlabelled examples.

    Website:https://ai.stanford.edu/~amaas/data/sentiment/

    When using this dataset, please cite this paper
    http://www.aclweb.org/anthology/P11-1015
    """,
    "http://ai.stanford.edu/~amaas/data/sentiment/aclImdb_v1.tar.gz",
    "c40f74a18d3b61f90feba1e17730e0d38e8b97c05fde7008942e91923d1658fe",
    post_fetch_method = unpack
))


text:
using Distributed
@everywhere begin
    using Pkg
    # Pkg.activate("C:/users/boltg/.julia/dev/")
    Pkg.activate("/home/boltg/.julia/dev/")
end
using Distributed, JLD2, StatsBase, Distances
@everywhere using InteractionNetworkModels, Distributions

# Script inputs
if (length(ARGS)==0)
    n_samples=10
    n_reps = nworkers()
    n_train = 20
elseif (length(ARGS)==1)
    n_samples = parse(Int, ARGS[1])
    n_reps = nworkers()
    n_train = 20
elseif (length(ARGS)==2)
    n_train = parse(Int, ARGS[1])
    n_samples = parse(Int, ARGS[2])
    n_reps = nworkers()
else
    n_train = parse(Int, ARGS[1])
    n_samples = parse(Int, ARGS[2])
    n_reps = parse(Int, ARGS[3])
end

println("Making model...")

E = [[1,2,1,2],
    [1,2,1],
    [3,4,3],
    [3,4],
    [1,2],
    [1,2,1],
    [1,2,3],
    [4,5],
    [7,8]
]
γ = 2.6
V = 1:10
d = MatchingDist(FastLCS(101))
K_inner, K_outer = (DimensionRange(2,100), DimensionRange(1,25))

model = SIM(
    E, γ,
    d,
    V,
    K_inner, K_outer
)

mcmc_sampler = SimMcmcInsertDelete(
    ν_ed=6, ν_td=1, β=0.7,
    len_dist=TrGeometric(0.8, 1, model.K_inner.u),
    burn_in=5000, lag=50
)

println("Sampling data...")
@time mcmc_out = mcmc_sampler(
    model,
    desired_samples=50,
    lag=500,
    burn_in=10000
)

println("Constructing posterior...")
data = mcmc_out.sample
E_prior = SIM(E, 0.1, model.dist, model.V, model.K_inner, model.K_outer)
γ_prior = Uniform(0.5,7.0)
posterior = SimPosterior(data, E_prior, γ_prior)

# Construct posterior sampler
posterior_sampler = SimIexInsertDelete(
    mcmc_sampler,
    ν_ed=1, ν_td=1, len_dist=TrGeometric(0.8, K_inner.l, K_inner.u),
    ε=0.1,
    K=200,
    α=0.0,
    desired_samples=n_samples, burn_in=0, lag=1
)

println("Setting-up mapper to workers....")

@everywhere function f(
    posterior::SimPosterior,
    mcmc_sampler::SimPosteriorSampler,
    S_init::InteractionSequence,
    γ_init::Float64
    )

    # Small run for precompile
    @time mcmc_sampler(
        posterior,
        S_init=S_init,
        γ_init=γ_init,
        desired_samples=5,
        burn_in=0, lag=1,
        loading_bar=false,
        aux_init_at_prev=true
    )
    out = @timed mcmc_sampler(
        posterior,
        S_init=S_init,
        γ_init=γ_init,
        loading_bar=true,
        aux_init_at_prev=true
    )
    return (chain=out[1], time=out[2])
end

Ê = sample_frechet_mean(posterior.data, posterior.dist)[1]
pars = [(
    post=posterior,
    post_sampler=posterior_sampler,
    S_init=Ê,
    γ_init=γ
    ) for i in 1:n_reps
]

println("Sampling....")
out = pmap(x->f(x.post, x.post_sampler, x.S_init, x.γ_init), pars)

chains = [x.chain for x in out]
times = [x.time/60/60 for x in out]
# cd("Z:/simulations/")
cd("/luna/simulations")
save(
    "iex_sim_$(n_train)_match_dist_multi.jld2",
    "chains", chains,
    "times", times
)


text:
include("../src/tDATAsets.jl")

using Plots
default(size=(400, 400))

InftySymbol = rand(Infty(), 100)
x,y = eachrow(InftySymbol)
scatter(x,y, size = (400, 200))

S1 = rand(Sphere(), 100)
x,y = eachrow(S1)
scatter(x,y)

Plots.plotly()
default(msw = 0, ma = 0.5, ms = 1, mc = :black)

S2 = rand(Sphere(3), 1000)
x,y,z = eachrow(S2)
scatter(x,y,z)

SR = rand(SwissRoll(), 1000)
x,y,z = eachrow(SR)
scatter(x,y,z)

T = rand(Torus(), 1000, noise = 0.0)
x,y,z = eachrow(T)
scatter(x,y,z, lims = (-6,6), legend = :none)


text:
import Colors

export show_plan

### DAG-based graphing

global _part_labels = Dict()

function write_node(ctx, io, t::Chunk, c)
    _part_labels[t]="part_$c"
    c+1
end

function node_id(t::Thunk)
    t.id
end

function node_id(t)
    dec(hash(t))
end

function node_id(t::Chunk)
    _part_labels[t]
end

function write_dag(io, t::Thunk)
    !istask(t) && return
    deps = dependents(t)
    c=1
    for k in keys(deps)
        c = write_node(nothing, io, k, c)
    end
    for (k, v) in deps
        for dep in v
            if isa(k, Union{Chunk, Thunk})
                println(io, "$(node_id(k)) -> $(node_id(dep))")
            end
        end
    end
end

### Timespan-based graphing

pretty_time(ts::Timespan) = pretty_time(ts.finish-ts.start)
function pretty_time(t)
    r(t) = round(t; digits=3)
    if t > 1000^3
        "$(r(t/(1000^3))) s"
    elseif t > 1000^2
        "$(r(t/(1000^2))) ms"
    elseif t > 1000
        "$(r(t/1000)) us"
    else
        "$(r(t)) ns"
    end
end
function pretty_size(sz)
    if sz > 1024^4
        "$(sz/(1024^4)) TB (terabytes)"
    elseif sz > 1024^3
        "$(sz/(1024^3)) GB (gigabytes)"
    elseif sz > 1024^2
        "$(sz/(1024^2)) MB (megabytes)"
    elseif sz > 1024
        "$(sz/1024) KB (kilobytes)"
    else
        "$sz B (bytes)"
    end
end

node_label(x) = repr(x)
node_label(x::T) where T<:AbstractArray =
    "$T\nShape: $(size(x))\nSize: $(pretty_size(sizeof(x)))"
node_label(x::Chunk) = "Chunk on $(x.processor)"

node_proc(x) = nothing
node_proc(x::Chunk) = x.processor

_proc_color(ctx, proc::Processor) = get!(ctx.proc_to_color, proc) do
    _color = ctx.proc_colors[ctx.proc_color_idx[]]
    ctx.proc_color_idx[] = clamp(ctx.proc_color_idx[]+1, 0, 128)
    "#$(Colors.hex(_color))"
end
_proc_color(ctx, id::Int) = _proc_color(ctx, ctx.id_to_proc[id])
_proc_color(ctx, ::Nothing) = "black"
_proc_shape(ctx, proc::Processor) = get!(ctx.proc_to_shape, typeof(proc)) do
    _shape = ctx.proc_shapes[ctx.proc_shape_idx[]]
    ctx.proc_shape_idx[] = clamp(ctx.proc_shape_idx[]+1, 0, length(ctx.proc_shapes))
    _shape
end
_proc_shape(ctx, ::Nothing) = "ellipse"

function write_node(ctx, io, t::Thunk, c)
    f = isa(t.f, Function) ? "$(t.f)" : "fn"
    println(io, "n_$(t.id) [label=\"$f - $(t.id)\"];")
    c
end

dec(x) = Base.dec(x, 0, false)
function write_node(ctx, io, t, c, id=dec(hash(t)))
    l = replace(node_label(t), "\""=>"")
    proc = node_proc(t)
    color = _proc_color(ctx, proc)
    shape = _proc_shape(ctx, proc)
    println(io, "n_$id [label=\"$l\",color=\"$color\",shape=\"$shape\",penwidth=5];")
    c
end

function write_node(ctx, io, ts::Timespan, c)
    f, proc = ts.timeline
    f = isa(f, Function) ? "$f" : "fn"
    t_comp = pretty_time(ts)
    color = _proc_color(ctx, proc)
    shape = _proc_shape(ctx, proc)
    # TODO: t_log = log(ts.finish - ts.start) / 5
    ctx.id_to_proc[ts.id] = proc
    println(io, "n_$(ts.id) [label=\"$f\n$t_comp\",color=\"$color\",shape=\"$shape\",penwidth=5];")
    # TODO: "\n Thunk $(ts.id)\nResult Type: $res_type\nResult Size: $sz_comp\",
    c
end

function write_edge(ctx, io, ts_move::Timespan, logs, inputname=nothing, inputarg=nothing)
    f, id = ts_move.timeline
    t_move = pretty_time(ts_move)
    if id > 0
        print(io, "n_$id -> n_$(ts_move.id[1]) [label=\"Move: $t_move")
        color_src = _proc_color(ctx, id)
    else
        @assert inputname !== nothing
        @assert inputarg !== nothing
        print(io, "n_$inputname -> n_$(ts_move.id[1]) [label=\"Move: $t_move")
        proc = node_proc(inputarg)
        color_src = _proc_color(ctx, proc)
    end
    color_dst = _proc_color(ctx, ts_move.id[1])
    # TODO: log_t = log(ts_move.finish-ts_move.start) / 5
    println(io, "\",color=\"$color_src;0.5:$color_dst\",penwidth=2];")
end

write_edge(ctx, io, from::String, to::String) = println(io, "n_$from -> n_$to;")

getargs!(d, t) = nothing
function getargs!(d, t::Thunk)
    d[t.id] = [filter(x->!istask(x[2]), collect(enumerate(t.inputs)))...,]
    foreach(i->getargs!(d, i), t.inputs)
end
function write_dag(io, logs::Vector, t=nothing)
    ctx = (proc_to_color = Dict{Processor,String}(),
           proc_colors = Colors.distinguishable_colors(128),
           proc_color_idx = Ref{Int}(1),
           proc_to_shape = Dict{Type,String}(),
           proc_shapes = ("ellipse","box","triangle"),
           proc_shape_idx = Ref{Int}(1),
           id_to_proc = Dict{Int,Processor}())
    argmap = Dict{Int,Vector}()
    getargs!(argmap, t)
    c = 1
    # Compute nodes
    for ts in filter(x->x.category==:compute, logs)
        c = write_node(ctx, io, ts, c)
    end
    # Argument nodes
    argnodemap = Dict{Int,Vector{String}}()
    argids = IdDict{Any,String}()
    for id in keys(argmap)
        nodes = String[]
        arg_c = 1
        for (argidx,arg) in argmap[id]
            name = "arg_$(argidx)_to_$(id)"
            if !isimmutable(arg)
                if arg in keys(argids)
                    name = argids[arg]
                else
                    argids[arg] = name
                    c = write_node(ctx, io, arg, c, name)
                end
                push!(nodes, name)
            else
                c = write_node(ctx, io, arg, c, name)
                push!(nodes, name)
            end
            # Arg-to-compute edges
            for ts in filter(x->x.category==:move &&
                                x.id[1]==id &&
                                x.timeline[2]==-argidx, logs)
                write_edge(ctx, io, ts, logs, name, arg)
            end
            arg_c += 1
        end
        argnodemap[id] = nodes
    end
    # Move edges
    for ts in filter(x->x.category==:move && x.timeline[2]>0, logs)
        write_edge(ctx, io, ts, logs)
    end
    #= FIXME: Legend (currently it's laid out horizontally)
    println(io, """
    subgraph  {
        graph[style=dotted,newrank=true,rankdir=TB];
        edge [style=invis];
        Legend [shape=box];""")
    cur = "Legend"
    for id in keys(ctx.id_to_proc)
        color = _proc_color(ctx, id)
        shape = _proc_shape(ctx, ctx.id_to_proc[id])
        name = "proc_$id"
        println(io, "$name [color=\"$color\",shape=\"$shape\"];")
        println(io, "$cur -> $name;")
        cur = name
    end
    println(io, "}")
    =#
end

function show_plan(io::IO, t)
    println(io, """strict digraph {
    graph [layout=dot,rankdir=LR];""")
    write_dag(io, t)
    println(io, "}")
end
function show_plan(io::IO, logs::Vector{Timespan}, t::Thunk)
    println(io, """strict digraph {
    graph [layout=dot,rankdir=LR];""")
    write_dag(io, logs, t)
    println(io, "}")
end

function show_plan(t::Union{Thunk,Vector{Timespan}})
    io = IOBuffer()
    show_plan(io, t)
    return String(take!(io))
end
function show_plan(logs::Vector{Timespan}, t::Thunk)
    io = IOBuffer()
    show_plan(io, logs, t)
    return String(take!(io))
end

function show_plan(c)
    ctx = Context(global_context())
    t = thunkize(ctx, stage(ctx, c))
    show_plan(t)
end

function show_plan(t::Tuple)
    show_plan(TupleCompute(t))
end

# function printing

argname(x::Symbol) = x

function argname(x)
    @assert x.head == :(::)
    x.args[1]
end

function show_statement(io, x)
    if x.head == :return
        x = x.args[1]
    end
    print(io, x)
end

function fnbody(io, x)
    body = x.args[3]
    statements = filter(x -> !isa(x, LineNumberNode), body.args)
    for s in statements
        show_statement(io, s)
    end
end

function show_ast(io, f)
    ast = Base.uncompressed_ast(f.code)
    args = map(x -> string(argname(x)), ast.args[1])
    write(io, "(", join(args, ','), ") -> ")
    fnbody(io, ast)
end

function showfn(io, f::Function)
    if isa(f, Function)
        show(io, f)
    else
        show_ast(io, f)
    end
end
showfn(io, f) = show(io, f)


text:
val = 42

using Jive

@onlyonce begin
    sleep(0.1)
    val = 0
end


text:
import ..TestSuite

@testset "Precomputed coefficients" begin
    TestSuite.runtests(LegendreUnitCoeff{Float64}(TestSuite.LMAX))
    TestSuite.runtests(LegendreUnitCoeff{Float32}(TestSuite.LMAX))
    TestSuite.runtests(LegendreUnitCoeff{BigFloat}(TestSuite.LMAX))

    @testset "Coefficient table conversion" begin
        dtab = LegendreSphereCoeff{Float64}(10)
        ftab = convert(LegendreSphereCoeff{Float32}, dtab)

        @test ftab isa LegendreSphereCoeff{Float32}
        @test ftab.α == Float32.(dtab.α)
        @test_throws MethodError convert(LegendreUnitCoeff{Float64}, dtab)
    end

    @testset "Setting mmax" begin
        LMAX = 10
        MMAX = 2
        impltab = LegendreUnitCoeff{Float64}(LMAX)
        fulltab = LegendreUnitCoeff{Float64}(LMAX, LMAX)
        mmaxtab = LegendreUnitCoeff{Float64}(LMAX, MMAX)

        # Equality of implicit and explicit [maximal] mmax
        @test impltab.α == fulltab.α
        @test impltab.α == fulltab.α
        @test impltab.β == fulltab.β
        @test impltab.μ == fulltab.μ
        @test impltab.ν == fulltab.ν

        # Equality of coefficients up to mmax for mmax limited
        @test mmaxtab.α == @view fulltab.α[:, 1:(MMAX+1)]
        @test mmaxtab.β == @view fulltab.β[:, 1:(MMAX+1)]
        @test mmaxtab.μ == @view fulltab.μ[1:(MMAX+1)]
        @test mmaxtab.ν == @view fulltab.ν[1:(MMAX+1)]
    end
end


text:
function prime(n)
c=0
if n == 1
    println("1 is Neither Prime nor Composite")
else
    for i=2:sqrt(n)
        if n%i == 0
            c=c+1
        end
    end
    if c==0
        println(n, " is a Prime")
    else
        println(n," is Not a Prime")
    end
end
end

i=735997
(prime(i))

max_stars_repo_name: landrumb/julia

max_issues_repo_name: landrumb/julia

max_forks_repo_name: landrumb/julia


text:
using Base.Sys
using DFControl.Utils: searchdir
cd(assetfile("Fe/"))

big_files = searchdir(".", "bz2")
if islinux()
    for f in big_files
        run(`bunzip2 $f`)
    end
end

job = Job(".")
abgrid = DFW.AbInitioKGrid(job)

big_files = map(x -> splitext(x)[1], big_files)
if islinux()
    for f in big_files
        run(`bzip2 $f`)
    end
end


text:
abstract GraphLoader

function loadmissingdata(loader::GraphLoader, items::Set{Container})
	# given a partially loaded set of containers (edges or vertices), loading missing data
	throw(NotImplementedException()) # must be implemented for specific graph loader
end

function loadmissingdata(loader::GraphLoader, results::Set{CompositeQueryResultItem})
	# given a partially loaded set of result items, loading missing data
	throw(NotImplementedException()) # must be implemented for specific graph loader
end

function loadvertices(loader::GraphLoader)
	throw(NotImplementedException()) # must be implemented for specific graph loader
end

function loadvertices(loader::GraphLoader, where::Function)
	throw(NotImplementedException()) # must be implemented for specific graph loader
end

function loadvertices(loader::GraphLoader, whereExp::Expr)
	throw(NotImplementedException()) # must be implemented for specific graph loader
end

function loadedges(loader::GraphLoader)
	throw(NotImplementedException()) # must be implemented for specific graph loader
end

function loadedges(loader::GraphLoader, where::Function)
	throw(NotImplementedException()) # must be implemented for specific graph loader
end

function loadedges(loader::GraphLoader, whereExp::Expr)
	throw(NotImplementedException()) # must be implemented for specific graph loader
end

function vertices(loader::GraphLoader)
	# Operation that selects the vertices of a graph.  A CompositeQuery which begins with this operation is constructed.

	cq = vertices(loader.graph)
	cq.loader = loader

	return cq
end

function vertices(loader::GraphLoader, where::Function)
	# Operation that selects the vertices of a graph that match a supplied where condition.  A CompositeQuery which begins with this operation is constructed.
	cq = vertices(loader.graph, where)
	cq.loader = loader

	return cq
end

function vertices(loader::GraphLoader, whereExp::Expr)
	# Operation that selects the vertices of a graph that match a supplied where condition.  A CompositeQuery which begins with this operation is constructed.

	cq = vertices(loader.graph, whereExp)
	cq.loader = loader

	return cq
end

function edges(loader::GraphLoader)
	# Operation that selects the edges of a graph.  A CompositeQuery which begins with this operation is constructed.

	cq = edges(loader.graph)
	cq.loader = loader

	return cq
end

function edges(loader::GraphLoader, where::Function)
	# Operation that selects the edges of a graph that match a specified where condition.  A CompositeQuery which begins with this operation is constructed.

	cq = edges(loader.graph, where)
	cq.loader = loader

	return cq
end

function edges(loader::GraphLoader, whereExpr::Function)
	# Operation that selects the edges of a graph that match a specified where condition.  A CompositeQuery which begins with this operation is constructed.

	cq = edges(loader.graph, whereExpr)
	cq.loader = loader

	return cq
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct Method <: ProtoType
    name::AbstractString
    request_type_url::AbstractString
    request_streaming::Bool
    response_type_url::AbstractString
    response_streaming::Bool
    options::Base.Vector{Option}
    syntax::Int32
    Method(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Method

mutable struct Mixin <: ProtoType
    name::AbstractString
    root::AbstractString
    Mixin(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Mixin

mutable struct Api <: ProtoType
    name::AbstractString
    methods::Base.Vector{Method}
    options::Base.Vector{Option}
    version::AbstractString
    source_context::SourceContext
    mixins::Base.Vector{Mixin}
    syntax::Int32
    Api(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Api

export Api, Method, Mixin


text:
"""
    Buffer

Struct used for reading out vertices from Action struct.
Contains symmetry related flags, (asymptotic) kernel specification and interpolation parameters.
* `exchange_flag :: Bool`  : flag for site [(i0, j) <-> (j, i0)] and spin [(μ, ν) <-> (ν, μ)] exchange
* `map_flag      :: Bool`  : flag for channel mapping s <-> u
* `sgn_μν        :: Bool`  : sign function for combined spin indices
* `sgn_μ         :: Bool`  : sign function for first spin index
* `sgn_ν         :: Bool`  : sign function for second spin index
* `kernel        :: Int64` : specification of asymptotic kernel to be interpolated
* `p1            :: Param` : interpolation parameters for bosonic frequency argument
* `p2            :: Param` : interpolation parameters for first fermionic frequency argument
* `p3            :: Param` : interpolation parameters for second fermionic frequency argument
"""
struct Buffer
    exchange_flag :: Bool
    map_flag      :: Bool
    sgn_μν        :: Bool
    sgn_μ         :: Bool
    sgn_ν         :: Bool
    kernel        :: Int64
    p1            :: Param
    p2            :: Param
    p3            :: Param
end

# function to invert flag
function set_flag(
    flag :: Bool
    )    :: Bool

    if flag
        flag = false
    else
        flag = true
    end

    return flag
end

# generate buffer dummy
function get_buffer_empty() :: Buffer

    b = Buffer(false, false, false, false, false, 0, get_param_empty(), get_param_empty(), get_param_empty())

    return b
end

# generate generic access buffer for Action struct given flags
function get_buffer(
    w             :: Float64,
    v             :: Float64,
    vp            :: Float64,
    Ω             :: Vector{Float64},
    ν             :: Vector{Float64},
    exchange_flag :: Bool,
    map_flag      :: Bool,
    sgn_μν        :: Bool,
    sgn_μ         :: Bool,
    sgn_ν         :: Bool
    )             :: Buffer

    if Ω[end] < abs(w)
        return get_buffer_empty()
    else
        if ν[end] < abs(v)
            if ν[end] < abs(vp)
                # interpolation for q1
                return Buffer(exchange_flag, map_flag, sgn_μν, sgn_μ, sgn_ν, 1, get_param(w, Ω), get_param_empty(), get_param_empty())
            else
                # interpolation for q2_2
                return Buffer(exchange_flag, map_flag, sgn_μν, sgn_μ, sgn_ν, 3, get_param(w, Ω), get_param_empty(), get_param(vp, ν))
            end
        else
            if ν[end] < abs(vp)
                # interpolation for q2_1
                return Buffer(exchange_flag, map_flag, sgn_μν, sgn_μ, sgn_ν, 2, get_param(w, Ω), get_param(v, ν), get_param_empty())
            else
                # interpolation for q3
                return Buffer(exchange_flag, map_flag, sgn_μν, sgn_μ, sgn_ν, 4, get_param(w, Ω), get_param(v, ν), get_param(vp, ν))
            end
        end
    end
end

# generate access buffer for s channel of Action struct
function get_buffer_s(
    w  :: Float64,
    v  :: Float64,
    vp :: Float64,
    m  :: Mesh
    )  :: Buffer

    # init flags
    exchange_flag = false
    map_flag      = false
    sgn_μν        = false
    sgn_μ         = false
    sgn_ν         = false

    # do -w -> w
    if w < 0.0
        w             *= -1.0
        exchange_flag  = set_flag(exchange_flag)
    end

    # do -v -> v
    if v < 0.0
        v             *= -1.0
        exchange_flag  = set_flag(exchange_flag)
        map_flag       = set_flag(map_flag)
        sgn_μ          = set_flag(sgn_μ)
    end

    # do -vp -> vp
    if vp < 0.0
        vp       *= -1.0
        map_flag  = set_flag(map_flag)
        sgn_ν     = set_flag(sgn_ν)
    end

    return get_buffer(w, v, vp, m.Ωs, m.νs, exchange_flag, map_flag, sgn_μν, sgn_μ, sgn_ν)
end

# generate access buffer for t channel of Action struct
function get_buffer_t(
    w  :: Float64,
    v  :: Float64,
    vp :: Float64,
    m  :: Mesh
    )  :: Buffer

    # init flags
    exchange_flag = false
    map_flag      = false
    sgn_μν        = false
    sgn_μ         = false
    sgn_ν         = false

    # do -w -> w
    if w < 0.0
        w      *= -1.0
        sgn_μν  = set_flag(sgn_μν)
    end

    # do -v -> v
    if v < 0.0
        v     *= -1.0
        sgn_μ  = set_flag(sgn_μ)
    end

    # do -vp -> vp
    if vp < 0.0
        vp    *= -1.0
        sgn_ν  = set_flag(sgn_ν)
    end

    return get_buffer(w, v, vp, m.Ωt, m.νt, exchange_flag, map_flag, sgn_μν, sgn_μ, sgn_ν)
end

# generate access buffer for u channel of Action struct
function get_buffer_u(
    w  :: Float64,
    v  :: Float64,
    vp :: Float64,
    m  :: Mesh
    )  :: Buffer

    # init flags
    exchange_flag = false
    map_flag      = false
    sgn_μν        = false
    sgn_μ         = false
    sgn_ν         = false

    # do -w -> w
    if w < 0.0
        w             *= -1.0
        exchange_flag  = set_flag(exchange_flag)
        sgn_μν         = set_flag(sgn_μν)
    end

    # do -v -> v
    if v < 0.0
        v             *= -1.0
        exchange_flag  = set_flag(exchange_flag)
        map_flag       = set_flag(map_flag)
        sgn_ν          = set_flag(sgn_ν)
    end

    # do -vp -> vp
    if vp < 0.0
        vp       *= -1.0
        map_flag  = set_flag(map_flag)
        sgn_ν     = set_flag(sgn_ν)
    end

    return get_buffer(w, v, vp, m.Ωs, m.νs, exchange_flag, map_flag, sgn_μν, sgn_μ, sgn_ν)
end


text:
function LinearAlgebra.mul!(
    C::Matrix{T},
    blkA::BlockedSparse{T},
    adjB::Adjoint{T,<:BlockedSparse{T}},
    α::Number,
    β::Number,
) where {T}
    A = blkA.cscmat
    B = adjB.parent.cscmat
    B.m == size(C, 2) && A.m == size(C, 1) && A.n == B.n || throw(DimensionMismatch(""))
    anz = nonzeros(A)
    arv = rowvals(A)
    bnz = nonzeros(B)
    brv = rowvals(B)
    isone(β) || rmul!(C, β)
    @inbounds for j = 1:A.n
        for ib in nzrange(B, j)
            αbnz = α * bnz[ib]
            jj = brv[ib]
            for ia in nzrange(A, j)
                C[arv[ia], jj] += anz[ia] * αbnz
            end
        end
    end
    C
end

LinearAlgebra.mul!(
    C::StridedVecOrMat{T},
    A::StridedVecOrMat{T},
    adjB::Adjoint{T,<:BlockedSparse{T}},
    α::Number,
    β::Number,
) where {T} = mul!(C, A, adjB.parent.cscmat', α, β)

LinearAlgebra.mul!(
    C::StridedVector{T},
    adjA::Adjoint{T,<:BlockedSparse{T}},
    B::StridedVector{T},
    α::Number,
    β::Number,
) where {T} = mul!(C, adjA.parent.cscmat', B, α, β)

function LinearAlgebra.ldiv!(
    adjA::Adjoint{T,<:LowerTriangular{T,UniformBlockDiagonal{T}}},
    B::StridedVector{T},
) where {T}
    A = adjA.parent
    length(B) == size(A, 2) || throw(DimensionMismatch(""))
    Adat = A.data.data
    m, n, k = size(Adat)
    bb = reshape(B, (n, k))
    for j in axes(Adat, 3)
        ldiv!(adjoint(LowerTriangular(view(Adat, :, :, j))), view(bb, :, j))
    end
    B
end

function LinearAlgebra.rdiv!(
    A::Matrix{T},
    adjB::Adjoint{T,<:LowerTriangular{T,UniformBlockDiagonal{T}}},
) where {T}
    m, n = size(A)
    Bd = adjB.parent.data
    Bdd = Bd.data
    r, s, blk = size(Bdd)
    n == size(Bd, 1) && r == s || throw(DimensionMismatch())
    for b = axes(Bd.data, 3)
        coloffset = (b - 1) * s
        rdiv!(view(A, :, coloffset+1:coloffset+s), adjoint(LowerTriangular(view(Bdd, :, :, b))))
    end
    A
end

function LinearAlgebra.rdiv!(
    A::BlockedSparse{T,S,P},
    B::Adjoint{T,<:LowerTriangular{T,UniformBlockDiagonal{T}}},
) where {T,S,P}
    Bpd = B.parent.data
    Bdat = Bpd.data
    j, k, l = size(Bdat)
    cbpt = A.colblkptr
    nzv = A.cscmat.nzval
    P == j == k && length(cbpt) == (l + 1) || throw(DimensionMismatch(""))
    for j in axes(Bdat, 3)
        rdiv!(
            reshape(view(nzv,cbpt[j]:(cbpt[j+1]-1)),:,P),
            adjoint(LowerTriangular(view(Bdat,:,:,j)))
            )
    end
    A
end


text:
module CobylaTests

using Printf
using OptimPack.Powell

function runtests(;revcom::Bool = false, scale::Real = 1.0)
    # Beware that order of operations may affect the result (whithin
    # rounding errors).  I have tried to keep the same ordering as F2C
    # which takes care of that, in particular when converting expressions
    # involving powers.
    prt(s) = println("\n       "*s)
    for nprob in 1:10
        if nprob == 1
            # Minimization of a simple quadratic function of two variables.
            prt("Output from test problem 1 (Simple quadratic)")
            n = 2
            m = 0
            xopt = Array{Cdouble}(undef, n)
            xopt[1] = -1.0
            xopt[2] = 0.0
            ftest = (x::DenseVector{Cdouble}) -> begin
                r1 = x[1] + 1.0
                r2 = x[2]
                fc = 10.0*(r1*r1) + (r2*r2)
                return fc
            end
        elseif nprob == 2
            # Easy two dimensional minimization in unit circle.
            prt("Output from test problem 2 (2D unit circle calculation)")
            n = 2
            m = 1
            xopt = Array{Cdouble}(undef, n)
            xopt[1] = sqrt(0.5)
            xopt[2] = -xopt[1]
            ftest = (x::DenseVector{Cdouble}, con::DenseVector{Cdouble}) -> begin
                fc = x[1]*x[2]
                con[1] = 1.0 - x[1]*x[1] - x[2]*x[2]
                return fc
            end
        elseif nprob == 3
            # Easy three dimensional minimization in ellipsoid.
            prt("Output from test problem 3 (3D ellipsoid calculation)")
            n = 3
            m = 1
            xopt = Array{Cdouble}(undef, n)
            xopt[1] = 1.0/sqrt(3.0)
            xopt[2] = 1.0/sqrt(6.0)
            xopt[3] = -0.33333333333333331
            ftest = (x::DenseVector{Cdouble}, con::DenseVector{Cdouble}) -> begin
                fc = x[1]*x[2]*x[3]
                con[1] = 1.0 - (x[1]*x[1]) - 2.0*(x[2]*x[2]) - 3.0*(x[3]*x[3])
                return fc
            end
        elseif nprob == 4
            # Weak version of Rosenbrock's problem.
            prt("Output from test problem 4 (Weak Rosenbrock)")
            n = 2
            m = 0
            xopt = Array{Cdouble}(undef, n)
            xopt[1] = -1.0
            xopt[2] = 1.0
            ftest = (x::DenseVector{Cdouble}) -> begin
                r2 = x[1]
                r1 = r2*r2 - x[2]
                r3 = x[1] + 1.0
                fc = r1*r1 + r3*r3
                return fc
            end
        elseif nprob == 5
            # Intermediate version of Rosenbrock's problem.
            prt("Output from test problem 5 (Intermediate Rosenbrock)")
            n = 2
            m = 0
            xopt = Array{Cdouble}(undef, n)
            xopt[1] = -1.0
            xopt[2] = 1.0
            ftest = (x::DenseVector{Cdouble}) -> begin
                r2 = x[1]
                r1 = r2*r2 - x[2]
                r3 = x[1] + 1.0
                fc = r1*r1*10.0 + r3*r3
                return fc
            end
        elseif nprob == 6
            # This problem is taken from Fletcher's book Practical Methods
            # of Optimization and has the equation number (9.1.15).
            prt("Output from test problem 6 (Equation (9.1.15) in Fletcher)")
            n = 2
            m = 2
            xopt = Array{Cdouble}(undef, n)
            xopt[1] = sqrt(0.5)
            xopt[2] = xopt[1]
            ftest = (x::DenseVector{Cdouble}, con::DenseVector{Cdouble}) -> begin
                fc = -x[1] - x[2]
                r1 = x[1]
                con[1] = x[2] - r1*r1
                r1 = x[1]
                r2 = x[2]
                con[2] = 1.0 - r1*r1 - r2*r2
                return fc
            end
        elseif nprob == 7
            # This problem is taken from Fletcher's book Practical Methods
            # of Optimization and has the equation number (14.4.2).
            prt("Output from test problem 7 (Equation (14.4.2) in Fletcher)")
            n = 3
            m = 3
            xopt = Array{Cdouble}(undef, n)
            xopt[1] = 0.0
            xopt[2] = -3.0
            xopt[3] = -3.0
            ftest = (x::DenseVector{Cdouble}, con::DenseVector{Cdouble}) -> begin
                fc = x[3]
                con[1] = x[1]*5.0 - x[2] + x[3]
                r1 = x[1]
                r2 = x[2]
                con[2] = x[3] - r1*r1 - r2*r2 - x[2]*4.0
                con[3] = x[3] - x[1]*5.0 - x[2]
                return fc
            end
        elseif nprob == 8
            # This problem is taken from page 66 of Hock and Schittkowski's
            # book Test Examples for Nonlinear Programming Codes. It is
            # their test problem Number 43, and has the name Rosen-Suzuki.
            prt("Output from test problem 8 (Rosen-Suzuki)")
            n = 4
            m = 3
            xopt = Array{Cdouble}(undef, n)
            xopt[1] = 0.0
            xopt[2] = 1.0
            xopt[3] = 2.0
            xopt[4] = -1.0
            ftest = (x::DenseVector{Cdouble}, con::DenseVector{Cdouble}) -> begin
                r1 = x[1]
                r2 = x[2]
                r3 = x[3]
                r4 = x[4]
                fc = (r1*r1 + r2*r2 + r3*r3*2.0 + r4*r4 - x[1]*5.0
                      - x[2]*5.0 - x[3]*21.0 + x[4]*7.0)
                r1 = x[1]
                r2 = x[2]
                r3 = x[3]
                r4 = x[4]
                con[1] = (8.0 - r1*r1 - r2*r2 - r3*r3 - r4*r4 - x[1]
                          + x[2] - x[3] + x[4])
                r1 = x[1]
                r2 = x[2]
                r3 = x[3]
                r4 = x[4]
                con[2] = (10.0 - r1*r1 - r2*r2*2.0 - r3*r3 - r4*r4*2.0
                          + x[1] + x[4])
                r1 = x[1]
                r2 = x[2]
                r3 = x[3]
                con[3] = (5.0 - r1*r1*2.0 - r2*r2 - r3*r3 - x[1]*2.0
                          + x[2] + x[4])
                return fc
            end
        elseif nprob == 9
            # This problem is taken from page 111 of Hock and
            # Schittkowski's book Test Examples for Nonlinear Programming
            # Codes. It is their test problem Number 100.
            prt("Output from test problem 9 (Hock and Schittkowski 100)")
            n = 7
            m = 4
            xopt = Array{Cdouble}(undef, n)
            xopt[1] =  2.330499
            xopt[2] =  1.951372
            xopt[3] = -0.4775414
            xopt[4] =  4.365726
            xopt[5] = -0.624487
            xopt[6] =  1.038131
            xopt[7] =  1.594227
            ftest = (x::DenseVector{Cdouble}, con::DenseVector{Cdouble}) -> begin
                r1 = x[1] - 10.0
                r2 = x[2] - 12.0
                r3 = x[3]
                r3 *= r3
                r4 = x[4] - 11.0
                r5 = x[5]
                r5 *= r5
                r6 = x[6]
                r7 = x[7]
                r7 *= r7
                fc = (r1*r1 + r2*r2*5.0 + r3*r3 + r4*r4*3.0
                      + r5*(r5*r5)*10.0 + r6*r6*7.0 + r7*r7
                      - x[6]*4.0*x[7] - x[6]*10.0 - x[7]*8.0)
                r1 = x[1]
                r2 = x[2]
                r2 *= r2
                r3 = x[4]
                con[1] = (127.0 - r1*r1*2.0 - r2*r2*3.0 - x[3]
                          - r3*r3*4.0 - x[5]*5.0)
                r1 = x[3]
                con[2] = (282.0 - x[1]*7.0 - x[2]*3.0 - r1*r1*10.0
                          - x[4] + x[5])
                r1 = x[2]
                r2 = x[6]
                con[3] = (196.0 - x[1]*23.0 - r1*r1 - r2*r2*6.0
                          + x[7]*8.0)
                r1 = x[1]
                r2 = x[2]
                r3 = x[3]
                con[4] = (r1*r1*-4.0 - r2*r2 + x[1]*3.0*x[2]
                          - r3*r3*2.0 - x[6]*5.0 + x[7]*11.0)
                return fc
            end
        elseif nprob == 10
            # This problem is taken from page 415 of Luenberger's book
            # Applied Nonlinear Programming. It is to maximize the area of
            # a hexagon of unit diameter.
            prt("Output from test problem 10 (Hexagon area)")
            n = 9
            m = 14
            xopt = fill!(Array{Cdouble}(undef, n), 0.0)
            ftest = (x::DenseVector{Cdouble}, con::DenseVector{Cdouble}) -> begin
                fc = -0.5*(x[1]*x[4] - x[2]*x[3] + x[3]*x[9] - x[5]*x[9]
                           + x[5]*x[8] - x[6]*x[7])
                r1 = x[3]
                r2 = x[4]
                con[1] = 1.0 - r1*r1 - r2*r2
                r1 = x[9]
                con[2] = 1.0 - r1*r1
                r1 = x[5]
                r2 = x[6]
                con[3] = 1.0 - r1*r1 - r2*r2
                r1 = x[1]
                r2 = x[2] - x[9]
                con[4] = 1.0 - r1*r1 - r2*r2
                r1 = x[1] - x[5]
                r2 = x[2] - x[6]
                con[5] = 1.0 - r1*r1 - r2*r2
                r1 = x[1] - x[7]
                r2 = x[2] - x[8]
                con[6] = 1.0 - r1*r1 - r2*r2
                r1 = x[3] - x[5]
                r2 = x[4] - x[6]
                con[7] = 1.0 - r1*r1 - r2*r2
                r1 = x[3] - x[7]
                r2 = x[4] - x[8]
                con[8] = 1.0 - r1*r1 - r2*r2
                r1 = x[7]
                r2 = x[8] - x[9]
                con[9] = 1.0 - r1*r1 - r2*r2
                con[10] = x[1]*x[4] - x[2]*x[3]
                con[11] = x[3]*x[9]
                con[12] = -x[5]*x[9]
                con[13] = x[5]*x[8] - x[6]*x[7]
                con[14] = x[9]
                return fc
            end
        else
            error("bad problem number ($nprob)")
        end

        x = Array{Cdouble}(undef, n)
        for icase in 1:2
            fill!(x, 1.0)
            rhobeg = 0.5
            rhoend = (icase == 2 ? 1e-4 : 0.001)
            if revcom
                # Test the reverse communication variant.
                c = Array{Cdouble}(undef, max(m, 0))
                ctx = Cobyla.create(n, m, rhobeg, rhoend;
                                    verbose = 1, maxeval = 2000)
                status = getstatus(ctx)
                while status == Cobyla.ITERATE
                    if m > 0
                        # Some constraints.
                        fx = ftest(x, c)
                        status = iterate(ctx, fx, x, c)
                    else
                        # No constraints.
                        fx = ftest(x)
                        status = iterate(ctx, fx, x)
                    end
                end
                if status != Cobyla.SUCCESS
                    println("Something wrong occured in COBYLA: ",
                            getreason(status))
                end
            elseif scale == 1
                cobyla!(ftest, x, m, rhobeg, rhoend;
                        verbose = 1, maxeval = 2000)
            else
                Cobyla.minimize!(ftest, x, m, rhobeg/scale, rhoend/scale;
                                 scale = fill!(Array{Cdouble}(undef, n), scale),
                                 verbose = 1, maxeval = 2000)
            end
            if nprob == 10
                tempa = x[1] + x[3] + x[5] + x[7]
                tempb = x[2] + x[4] + x[6] + x[8]
                tempc = 0.5/sqrt(tempa*tempa + tempb*tempb)
                tempd = tempc*sqrt(3.0)
                xopt[1] = tempd*tempa + tempc*tempb
                xopt[2] = tempd*tempb - tempc*tempa
                xopt[3] = tempd*tempa - tempc*tempb
                xopt[4] = tempd*tempb + tempc*tempa
                for i in 1:4
                    xopt[i + 4] = xopt[i]
                end
            end
            temp = 0.0
            for i in 1:n
                r1 = x[i] - xopt[i]
                temp += r1*r1
            end
            @printf("\n     Least squares error in variables =%16.6E\n", sqrt(temp))
        end
        @printf("  ------------------------------------------------------------------\n")
    end
end

end # module


text:
using FstFiles
using FileIO
using DataFrames
using IterableTables
using Base.Test

@testset "FstFiles" begin

df = DataFrame(a=[1,2,3], b=[1.,2.,3.], c=["A", "B", "C"])

filename = joinpath(@__DIR__, "test.fst")

save(filename, df)

try

    fstfile = load(filename)

    df2 = DataFrame(fstfile)

    @test df == df2

finally
    rm(filename)
end

end


text:
import Base: isempty, push!, length

struct Empty <: Exception end
struct Full <: Exception end

struct Queue
    maxsize::Integer
    uselifo::Bool
    """
    mutex must be held whenever the queue is mutating.  All methods
    that acquire mutex must release it before returning.  mutex
    is shared between the two conditions, so acquiring and
    releasing the conditions also acquires and releases mutex.
    """
    mutex::ReentrantLock
    """
    Notify not_empty whenever an item is added to the queue; a
    thread waiting to get is notified then.
    """
    not_empty::Threads.Condition
    """
    Notify not_full whenever an item is removed from the queue;
    a thread waiting to put is notified then
    """
    not_full::Threads.Condition
    queue::Array
end

"""
   function Queue(maxsize::Integer=0, use_lifo::Bool=false)

Initialize a queue object with a given maximum size.

If `maxsize` is <= 0, the queue size is infinite.

If `use_lifo` is True, this Queue acts like a Stack (LIFO).
"""
function Queue(maxsize::Integer=0, uselifo::Bool=false)
    mutex = ReentrantLock()
    return Queue(maxsize, uselifo, mutex, Threads.Condition(mutex), Threads.Condition(mutex), [])
end

"""
    function length(q::Queue) :: Integer

Return the approximate size of the queue (not reliable!).
"""
function length(q::Queue) :: Integer
    lock(q.mutex)
    n = _length(q)
    unlock(q.mutex)
    return n
end

"""
    function isempty(q::Queue) :: Bool

Return True if the queue is empty, False otherwise (not reliable!).
"""
function isempty(q::Queue) :: Bool
    lock(q.mutex)
    n = _isempty(q)
    unlock(q.mutex)
    return n
end

"""
    function full(q::Queue) :: Bool

Return True if the queue is full, False otherwise (not reliable!).
"""

function full(q::Queue) :: Bool
    lock(q.mutex)
    n = _full(q)
    unlock(q.mutex)
    return n
end

"""

    function push!(q::Queue, item, block::Bool=true, timeout::Union{Integer, Nothing} = nothing)

Put an item into the queue.

        If optional args `block` is True and `timeout` is None (the
        functionault), block if necessary until a free slot is
        available. If `timeout` is a positive number, it blocks at
        most `timeout` seconds and raises the ``Full`` exception if no
        free slot was available within that time.  Otherwise (`block`
        is false), put an item on the queue if a free slot is
        immediately available, else raise the ``Full`` exception
        (`timeout` is ignored in that case).
        """

function push!(q::Queue, item, block::Bool=true, timeout::Union{Integer, Nothing} = nothing)
    lock(q.not_full)
    try
        if !block
            if _isfull(q)
                throw(Full())
            end
        elseif timeout == nothing
            while _isfull(q)
                wait(q.not_full)
            end
        else
            if timeout < 0
                throw(ArgumentError("'timeout' must be a positive number"))
              end
            endtime = _time() + timeout
            while _isfull(q)
                remaining = endtime - _time()
                if remaining <= 0.0
                    throw(Full())
                end
                wait(q.not_full)
                #self.not_full.wait(remaining)
            end
        end
        _push!(q, item)
        notify(q.not_empty)
     finally
        unlock(q.not_full)
     end
end

"""

Put an item into the queue without blocking.

Only enqueue the item if a free slot is immediately available.
Otherwise raise the ``Full`` exception.
        """
push!_nowait(q::Queue, item) = push!(q, item, false)
        """Remove and return an item from the queue.

        If optional args `block` is True and `timeout` is None (the
        functionault), block if necessary until an item is available. If
        `timeout` is a positive number, it blocks at most `timeout`
        seconds and raises the ``Empty`` exception if no item was
        available within that time.  Otherwise (`block` is false),
        return an item if one is immediately available, else raise the
        ``Empty`` exception (`timeout` is ignored in that case).
        """
    function get(q::Queue, block::Bool=true, timeout::Union{Integer, Nothing}=nothing)
        lock(q.not_empty)
        try
            if !block
                if _isempty(q)
                    raise(Empty)
                end
            elseif timeout === None
                while _isempty(q)
                    wait(q.not_empty)
                end
            else
                if timeout < 0
                    raise(ValueError("'timeout' must be a positive number"))
                end
                endtime = _time() + timeout
                while _isempty(q)
                    remaining = endtime - _time()
                    if remaining <= 0.0
                        raise(Empty)
                    end
                    wait(q.not_empty)
                    #self.not_empty.wait(remaining)
                end
            end
            item = _get(q)
            notify(q.not_full)
            return item
        finally
            unlock(q.not_empty)
        end
    end

      """Remove and return an item from the queue without blocking.

        Only get an item if one is immediately available. Otherwise
        raise the ``Empty`` exception.
        """

get_nowait(q::Queue) = get(q)

#Override these methods to implement other queue organizations
    # (e.g. stack or priority queue).
    # These will only be called with appropriate locks held

    # Initialize the queue representation

function _qsize(q::Queue) :: Integer
    return length(q.queue)
end

function _isempty(q::Queue) :: Bool
    return isempty(q.queue)
end

function _isfull(q::Queue) :: Bool
    return q.maxsize > 0 && length(q.queue) == q.maxsize
end

function _push!(q::Queue, item)
    push!(q.queue, item)
end

function _get(q::Queue)
    if q.use_lifo
        # LIFO
        return pop!(q.queue)
    else
        # FIFO
        return popfirst!(q.queue)
    end
end

max_stars_repo_name: lucianolorenti/StructDatabaseMapping

max_issues_repo_name: lucianolorenti/StructDatabaseMapping

max_forks_repo_name: lucianolorenti/StructDatabaseMapping


text:
###############################################################################
#
# Volume scattering functions adatpted from Volscatt version 2 by W. Verhoef
#
###############################################################################
"""
    volscatt!(cache::Array{FT,1},
              sza::FT,
              vza::FT,
              raa::FT,
              ttl::FT
    ) where {FT<:AbstractFloat}

Calculate interception parameters (`chi_s` and `chi_s`) and leaf reflectance
    multiplier (`frho`) and transmittance multiplier (`ftau`), given
- `cache` Array cache for results
- `sza` Solar zenith angle
- `vza` Viewing zenith angle
- `raa` Relative azimuth angle
- `ttl` Leaf inclination angle
"""
function volscatt!(
            cache::Array{FT,1},
            sza::FT,
            vza::FT,
            raa::FT,
            ttl::FT
) where {FT<:AbstractFloat}
    psi_rad = deg2rad(raa);
    cos_raa = cosd(raa);
    cos_ttl = cosd(ttl);
    sin_ttl = sind(ttl);
    cos_sza = cosd(sza);
    sin_sza = sind(sza);
    cos_vza = cosd(vza);
    sin_vza = sind(vza);
    Cs      = cos_ttl*cos_sza;
    Ss      = sin_ttl*sin_sza;
    Co      = cos_ttl*cos_vza;
    So      = sin_ttl*sin_vza;

    cosbts = FT(1.0);
    cosbto = FT(1.0);
    if (abs(Ss)>1e-6)
        cosbts = -Cs / Ss;
    end
    if (abs(So)>1e-6)
        cosbto = -Co / So;
    end

    if (abs(cosbts)<1)
        bts = acos(cosbts);
        ds  = Ss;
    else
        bts = FT(pi);
        ds  = Cs;
    end

    if abs(cosbto)<1
        bto = acos(cosbto);
        doo = So;
    elseif vza<90
        bto = FT(pi);
        doo = Co;
    else
        bto = 0;
        doo = -Co;
    end

    chi_s = 2 / FT(pi) * ( (bts - FT(pi)/2) * Cs + sin(bts) * Ss );
    chi_o = 2 / FT(pi) * ( (bto - FT(pi)/2) * Co + sin(bto) * So );

    # Computation of auxiliary azimut angles bt1, bt2, bt3 used
    # for the computation of the bidirectional scattering coefficient w

    btran1 = abs(bts - bto);
    btran2 = 2FT(pi) - bts - bto;
    #btran_=pi-abs(bts+bto-FT(pi))

    if psi_rad <= btran1
        bt1 = psi_rad;
        bt2 = btran1;
        bt3 = btran2;
    else
        bt1 = btran1;
        if psi_rad <= btran2
            bt2 = psi_rad;
            bt3 = btran2;
        else
            bt2 = btran2;
            bt3 = psi_rad;
        end
    end

    t1 = 2Cs * Co + Ss * So * cos_raa;
    t2 = 0;
    if bt2 > 0
        t2 = sin(bt2) * ( 2ds * doo + Ss * So * cos(bt1) * cos(bt3) );
    end

    denom = 2 * FT(pi)^2;
    frho  = ((pi-bt2) * t1 + t2) / denom;
    ftau  = (-bt2     * t1 + t2) / denom;

    # fill the values to cache
    cache[1] = chi_s;
    cache[2] = abs(chi_o);
    cache[3] = frho;
    cache[4] = ftau;

    return nothing
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land