text:
"""
    function csi()
---
## INPUT
- `K` : kernel matrix n × n
- `Y` : target vector n × d
- `m` : maximal rank
- `κ` : trade-off between approximation of `K` and prediction of `Y` (suggested: `.99`)
- `centering` : `1` if centering, `0` otherwise (suggested: 1)
- `δ` : number of columns of cholesky performed in advance (suggested: `40`)
- `tol` : minimum gain at iteration (suggested: `1e-4`)

## OUTPUT
- `G` : Cholesky decomposition -> K(P,P) is approximated by G*G'
- `P` : permutation matrix
- `Q`, `R` : QR decomposition of `G` (or `center(G)` if centering)
- `error1` : `tr(K-G*G')/tr(K)` at each step of the decomposition
- `error2` : `||Y-Q*Q'*Y||_F^2 / ||Y||_F^2` at each step of the decomposition
- `predicted_gain` : predicted gain before adding each column
- `true_gain` : actual gain after adding each column
"""
function csi(K, Y, m; centering = true, κ = 0.99, δ = 40, tol = 1e-4)
    n = size(K, 1)
    d = size(Y, 2)
    size(Y, 1) == n && @error "The targets have the wrong number of data points"
    tdm = min(m + δ, n)         # try dimension
    G = zeros(n, tdm)           # Cholesky factor
    diagK = diag(K)
    P = 1:n
    Q = zeros(n, tdm)           # Q of QR factors
    R = zeros(tdm, tdm)         # R of QR factors
    traceK = tr(K)
    λ = (1 - κ) / traceK
    centering && (Y .-= mean(Y, dims=1))
    sumY² = sum(Y.^2)
    μ = κ / sumY²
    # error1(1) = traceK  # to be confirmed if a vector
    # error2(1) = sumY2;
    k = 0               # current index of the Cholesky decomposition
    kadv = 0            # current index of the look ahead steps
    Dadv = diagK
    D = diagK
    # approximation cost cached quantities
    A1 = zeros(n, 1)
    A2 = zeros(n, 1)
    A3 = zeros(n, 1)
    GTG = zeros(m)
    QTY = zeros(m, d)
    QTYYTQ = zeros(m, m)
    # makes sure that delta is smaller than n
    δ = min(δ, n);

    for i=1:delta # first performs δ steps of Cholesky and QR decomposition
        kadv = kadv + 1
    end
        #=
        # select best index
        diagmax = Dadv(kadv)
        jast = 1
        for j in 1 : n - kadv + 1
            if Dadv(j + kadv - 1) > diagmax / .99
                diagmax = Dadv(j + kadv - 1)
                jast = j
            end
        end
        if diagmax<1e-12
            # all pivots are too close to zero, stops
            # this can only happen if the matrix has rank less than delta
            kadv = kadv - 1;
            break;
        else
            jast = jast + kadv - 1;

            # permute indices
            P( [kadv jast] )        = P( [jast kadv] );
            Dadv( [kadv jast] )     = Dadv( [jast kadv] );
            D( [kadv jast] )        = D( [jast kadv] );
            A1( [kadv jast] )        = A1( [jast kadv] );
            G([kadv jast],1:kadv-1) = G([ jast kadv],1:kadv-1);
            Q([kadv jast],1:kadv-1) = Q([ jast kadv],1:kadv-1);

            # compute new Cholesky column
            G(kadv,kadv)=Dadv(kadv);
            G(kadv,kadv)=sqrt(G(kadv,kadv));
            newKcol = K(P(kadv+1:n),P(kadv));
            G(kadv+1:n,kadv)=1/G(kadv,kadv)*( newKcol - G(kadv+1:n,1:kadv-1)*(G(kadv,1:kadv-1))');

            # update diagonal
            Dadv(kadv+1:n) =  Dadv(kadv+1:n) - G(kadv+1:n,kadv).^2;
            Dadv(kadv) = 0;

            # performs QR
            if centering
                Gcol = G(:,kadv) - 1/n * repmat( sum(G(:,kadv) ,1),n,1 );
            else
                Gcol = G(:,kadv);
            end
            R(1:kadv-1,kadv) = Q(:,1:kadv-1)' * Gcol;
            Q(:,kadv) = Gcol - Q(:,1:kadv-1) * R(1:kadv-1,kadv);
            R(kadv,kadv) = norm(Q(:,kadv));
            Q(:,kadv) = Q(:,kadv) / R(kadv,kadv);

            # update cached quantities
            if centering
                GTG(1:kadv,kadv) = G(:,1:kadv)' * G(:,kadv);
            else
                GTG(1:kadv,kadv) = R(1:kadv,1:kadv)' * R(1:kadv,kadv);
            end
            GTG(kadv,1:kadv) = GTG(1:kadv,kadv)';
            QTY(kadv,:) = Q(:,kadv)' * Y(P,:);
            QTYYTQ(kadv,1:kadv) = QTY(kadv,:) * QTY(1:kadv,:)';
            QTYYTQ(1:kadv,kadv) = QTYYTQ(kadv,1:kadv)';

            # update costs
            A1(kadv:n) = A1(kadv:n) + GTG(kadv,kadv) * ( G(kadv:n,kadv).^2 );
            A1(kadv:n) = A1(kadv:n) + 2 * G(kadv:n,kadv) .* ( G(kadv:n,1:kadv-1) * GTG(1:kadv-1,kadv) );

        end
    end

% compute remaining costs for all indices
A2 = sum( ( G(:,1:kadv) * ( R(1:kadv,1:kadv)' * QTY(1:kadv,:) ) ).^2 , 2);
A3 = sum( ( G(:,1:kadv) * R(1:kadv,1:kadv)' ).^2 , 2);

% start main loop
while k<m,

    k = k +1;

    % compute the gains in approximation for all remaining indices
    dJK = zeros(n-k+1,1);
    for i=1:n-k+1
        kast = k+i-1;
        if D(kast)<1e-12,
            % this column is already generated by already
            % selected columns -> cannot be selected
            dJK(i)=-1e100;
        else
            dJK(i) = A1(kast);
            if kast > kadv,
                % add eta
                dJK(i) = dJK(i) + D(kast)^2 - (D(kast) - Dadv(kast))^2;
            end
            dJK(i) = dJK(i) / D(kast);
        end
    end

    dJY = zeros(n-k+1,1);
    if kadv>k
        for i=1:n-k+1
            kast = k+i-1;
            if A3(kast) < 1e-12
                dJY(i) = 0;
            else
                dJY(i) = A2(kast) / A3(kast);
            end
        end
    end

    % select the best column
    dJ = lambda * dJK + mu * dJY;
    diagmax = -1;
    jast = 0;
    for j=1:n-k+1
        if D(j+k-1)>1e-12
            if dJ(j) > diagmax / 0.9
                diagmax = dJ(j);
                jast = j;
            end
        end
    end
    if jast==0,
        % no more good indices, exit
        k=k-1;
        break;
    end
    jast = jast + k - 1;
    predicted_gain(k) = diagmax;

    % performs one cholesky + QR step:
    % if new pivot not already selected, use pivot
    % otherwise, select new look ahead index that maximize Dadv

    if jast > kadv,
        newpivot = jast;
        jast = kadv + 1;
    else
        a = 1e-12;
        b = 0;
        for j=1:n-kadv
            if Dadv(j+kadv)>a/.99
                a=Dadv(j+kadv);
                b = j+kadv;
            end
        end
        if b==0, newpivot = 0;
        else newpivot = b;
        end
    end

    if newpivot > 0
        % performs steps

        kadv = kadv + 1;

        % permute
        P( [kadv newpivot] ) = P( [newpivot kadv] );
        Dadv( [kadv newpivot] ) = Dadv( [newpivot kadv] );
        D( [kadv newpivot] ) = D( [newpivot kadv] );
        A1( [kadv newpivot] ) = A1( [newpivot kadv] );
        A2( [kadv newpivot] ) = A2( [newpivot kadv] );
        A3( [kadv newpivot] ) = A3( [newpivot kadv] );
        G([kadv newpivot],1:kadv-1)=G([ newpivot kadv],1:kadv-1);
        Q([kadv newpivot],1:kadv-1)=Q([ newpivot kadv],1:kadv-1);

        % compute new Cholesky column
        G(kadv,kadv)=Dadv(kadv);
        G(kadv,kadv)=sqrt(G(kadv,kadv));
        newKcol = K(P(kadv+1:n),P(kadv));
        G(kadv+1:n,kadv)=1/G(kadv,kadv)*( newKcol - G(kadv+1:n,1:kadv-1)*(G(kadv,1:kadv-1))');

        % update diagonal
        Dadv(kadv+1:n) =  Dadv(kadv+1:n) - G(kadv+1:n,kadv).^2;
        Dadv(kadv) = 0;

        % performs QR
        if centering
            Gcol = G(:,kadv) - 1/n * repmat( sum(G(:,kadv) ,1),n,1 );
        else
            Gcol = G(:,kadv);
        end
        R(1:kadv-1,kadv) = Q(:,1:kadv-1)' * Gcol;
        Q(:,kadv) = Gcol - Q(:,1:kadv-1) * R(1:kadv-1,kadv);
        R(kadv,kadv) = norm(Q(:,kadv));
        Q(:,kadv) = Q(:,kadv) / R(kadv,kadv);

        % update the cached quantities
        if centering
            GTG(k:kadv,kadv) = G(:,k:kadv)' * G(:,kadv);
        else
            GTG(k:kadv,kadv) = R(1:kadv,k:kadv)' * R(1:kadv,kadv);
        end
        GTG(kadv,k:kadv) = GTG(k:kadv,kadv)';
        QTY(kadv,:) = Q(:,kadv)' * Y(P,:);
        QTYYTQ(kadv,k:kadv) = QTY(kadv,:) * QTY(k:kadv,:)';
        QTYYTQ(k:kadv,kadv) = QTYYTQ(kadv,k:kadv)';

        % update costs
        A1(kadv:n) = A1(kadv:n) + GTG(kadv,kadv) * ( G(kadv:n,kadv).^2 );
        A1(kadv:n) = A1(kadv:n) + 2 * G(kadv:n,kadv) .* ( G(kadv:n,k:kadv-1) * GTG(k:kadv-1,kadv) );
        A3(kadv:n) = A3(kadv:n) + G(kadv:n,kadv).^2 * sum( R(k:kadv,kadv).^2, 1 );
        temp = R(k:kadv,kadv)' * R(k:kadv,k:kadv-1);
        A3(kadv:n) = A3(kadv:n) + 2 *  G(kadv:n,kadv) .* ( G(kadv:n,k:kadv-1) * temp' );
        temp = R(k:kadv,kadv)' *  QTYYTQ(k:kadv,k:kadv);
        temp1 = temp * R(k:kadv,kadv) ;
        A2(kadv:n) = A2(kadv:n) + G(kadv:n,kadv).^2 * temp1;
        temp2 = temp * R(k:kadv,k:kadv-1);
        A2(kadv:n) = A2(kadv:n) + 2 * G(kadv:n,kadv) .* ( G(kadv:n,k:kadv-1) * temp2' );

    end

    % permute pivots in the Cholesky and QR decomposition between p,q
    p = k;
    q = jast;
    if p<q

        % store some quantities
        Gbef = G(:,p:q);
        Gbeftotal = G(:,k:kadv);
        GTGbef = GTG(p:q,p:q);
        QTYYTQbef = QTYYTQ(p:q,k:kadv);
        Rbef = R(p:q,p:q);
        Rbeftotal = R(k:kadv,k:kadv);
        tempG = eye(q-p+1,q-p+1);
        tempQ = eye(q-p+1,q-p+1);

        for s=q-1:-1:p
            % permute indices
            P( [s s+1] ) = P( [s+1 s] );
            Dadv( [s s+1] ) = Dadv( [s+1 s] );
            D( [s s+1] ) = D( [s+1 s] );
            A1( [s s+1] ) = A1( [s+1 s] );
            A2( [s s+1] ) = A2( [s+1 s] );
            A3( [s s+1] ) = A3( [s+1 s] );
            G([s s+1],1:kadv)=G([s+1 s], 1:kadv);
            Gbef([s s+1],:)=Gbef([s+1 s], :);
            Gbeftotal([s s+1],:)=Gbeftotal([s+1 s], :);
            Q([s s+1],1:kadv)=Q([s+1 s] ,1:kadv);

            % update decompositions
            [Q1,R1] = qr2( G(s:s+1,s:s+1)');
            G(:,s:s+1) = G(:,s:s+1) * Q1;
            G(s,s+1)=0;
            R(1:kadv,s:s+1) = R(1:kadv,s:s+1) * Q1;
            [Q2,R2] = qr2( R(s:s+1,s:s+1) );
            R(s:s+1,1:kadv) = Q2' * R(s:s+1,1:kadv);
            Q(:,s:s+1) = Q(:,s:s+1) * Q2;
            R(s+1,s)=0;

            % update relevant quantities
            nonchanged = [k:s-1 s+2:kadv ];
            GTG(nonchanged,s:s+1) = GTG(nonchanged,s:s+1) * Q1;
            GTG(s:s+1,nonchanged) = GTG(nonchanged,s:s+1)';
            GTG(s:s+1,s:s+1) = Q1' * GTG(s:s+1,s:s+1) * Q1;
            QTY(s:s+1,:) = Q2' * QTY(s:s+1,:);
            QTYYTQ(nonchanged,s:s+1) = QTYYTQ(nonchanged,s:s+1) * Q2;
            QTYYTQ(s:s+1,nonchanged) = QTYYTQ(nonchanged,s:s+1)';
            QTYYTQ(s:s+1,s:s+1) = Q2' * QTYYTQ(s:s+1,s:s+1) * Q2;
            tempG(:,s-p+1:s-p+2) = tempG(:,s-p+1:s-p+2) * Q1;
            tempQ(:,s-p+1:s-p+2) = tempQ(:,s-p+1:s-p+2) * Q2;
        end

        % update costs
        tempG = tempG(:,1);
        tempGG = GTGbef * tempG ;
        A1(k:n) = A1(k:n) - 2 * G(k:n,k) .* ( Gbef(k:n,:) * tempGG );                % between p and q -> different
        A1(k:n) = A1(k:n) - 2 * G(k:n,k) .* ( G(k:n,k:p-1) * GTG(k:p-1,k) );         % below p
        A1(k:n) = A1(k:n) - 2 * G(k:n,k) .* ( G(k:n,q+1:kadv) * GTG(q+1:kadv,k) );   % above q
        tempQ = tempQ(:,1);
        temp = G(k:n,q+1:kadv) * R(k,q+1:kadv)';
        temp = temp + G(k:n,k:p-1) * R(k,k:p-1)';
        temp = temp + Gbef(k:n,:) * Rbef' * tempQ;
        A3(k:n) = A3(k:n) - temp.^2 ;
        A2(k:n) = A2(k:n) + temp.^2 * QTYYTQ(k,k);
        temp2 = ( tempQ' * QTYYTQbef ) * Rbeftotal;
        A2(k:n) = A2(k:n) - 2 * temp .* ( Gbeftotal(k:n,:) * temp2' );

    else

                % update costs
        A1(k:n) = A1(k:n) - 2 * G(k:n,k) .* ( G(k:n,k:kadv) * GTG(k:kadv,k) );
        A3(k:n) = A3(k:n) - ( G(k:n,k:kadv) * R(k,k:kadv)' ).^2;
        temp = G(k:n,k:kadv) * R(k,k:kadv)';
        A2(k:n) = A2(k:n) + ( temp ).^2 * QTYYTQ(k,k);
        temp2 = QTYYTQ(k,k:kadv) * R(k:kadv,k:kadv);
        A2(k:n) = A2(k:n) - 2 * temp .* ( G(k:n,k:kadv) * temp2' );
    end

    % update diagonal and other quantities (A1,B1)
    D(k+1:n) =  D(k+1:n) - G(k+1:n,k).^2;
    D(k) = 0 ;
    A1(k:n) = A1(k:n) + GTG(k,k) * ( G(k:n,k).^2 );

    % compute errors and true gains
    temp2 =  Q(:,k)' * Y(P,:);
    temp2 = sum(temp2.^2);
    temp1 = sum( G(:,k).^2 );
    true_gain(k) = temp1 * lambda + temp2 * mu;
    error1(k+1) = error1(k) - temp1;
    error2(k+1) = error2(k) - temp2;
    if true_gain(k) < tol, break; end
end

% reduce dimensions of decomposition
G=G(:,1:k);
Q=Q(:,1:k);
R=R(1:k,1:k);

% compute and normalize errors
error = lambda * error1 + mu * error2;
error1 = error1 / traceK;
error2 = error2 / sumY2;
    =#
end

"""
    function csi_gaussian(x, α, Y, m; centering=1, κ=.99, δ=40, tol=1e-4)
---
## INPUT
- `x` : input data `nvar` × `n`
- `α` : kernel parameter of the RBF kernel
- `Y` : target vector `n` × `d`
- `m` : maximal rank
- `κ` : trade-off between approximation of K and prediction of Y (suggested: .99)
- `centering` : 1 if centering, 0 otherwise (suggested: 1)
- `δ` : number of columns of cholesky performed in advance (suggested: 40)
- `tol` : minimum gain at iteration (suggested: 1e-4)
## OUTPUT
- `G` : Cholesky decomposition → `K(P,P)` is approximated by `G × G'`
- `P` : permutation matrix
- `Q`, `R` : QR decomposition of `G` (or center(`G`) if centering)
- `error1` : `tr(K - G × G') / tr(K)` at each step of the decomposition
- `error2` : `||Y - Q × Q' × Y||_F^2 / || Y || _F^2` at each step of the decomposition
- `predicted_gain` : predicted gain before adding each column
- `true_gain` : actual gain after adding each column

"""
function csi_gaussian(x, α, Y, m;
                      centering=1, κ=.99, δ=40, tol=1e-4)
    K = exp(-α * sqdist(x, x))
    G, P, Q, R, error1, error2, error, predicted_gain, true_gain =
        csi(K, centering, κ, δ, tol)
    # return G, P, Q, R, error1, error2, error, predicted_gain, true_gain
end

function center_m(M)
    n = size(M, 1)
    for i in 1:size(M, 2)
        M[:, i] .-= mean(M[:, i])
    end
    M
end

"""
    function sqdist(a, b)
---
Computes squared Euclidean distance matrix,
a rectangular matrix of pairwise distances
between points in `a` (given in columns) and points in `b`.

- NB: very fast implementation taken from Roland Bunschoten
- XJ: my tests showed that `aa` and `bb` need to be centered.
"""
function sqdist(a, b)
    aa = sum(a .* a, dims = 1)
    bb = sum(b .* b, dims = 1)
    ab = a'b
    d  = abs.(repeat(aa', outer = (1, size(bb, 2))) +
              repeat(bb, outer = (size(aa, 2), 1)) -
              2ab)
end

"""
    function demo_csi()
---
This is a Julia version of [`csi`](https://www.di.ens.fr/~fbach/csi/csi-1.0.zip).
Here, I try to reproduce the method they described for [**Predictive low-rank decomposition for kernel methods**](https://www.di.ens.fr/~fbach/ICML_2005_CSI_3.pdf).
"""
function demo_csi()
    n, m = 500, 100
    n *= 2                      # why not just n = 1000?
    @debug "Dimensions" n m

    tm = [1 -.5; -.5 1]         # tiny matrix
    x1 = hcat(2tm * randn(2, n ÷ 4) + repeat([-3; 5], outer = (1, n ÷ 4)),
              tm * randn(2, n ÷ 4) + repeat([3; -.5], outer = (1, n ÷ 4)))
    x2 = hcat(tm * randn(2, n ÷ 4) + repeat([2; 3], outer = (1, n ÷ 4)),
              tm * randn(2, n ÷ 4) + repeat([6; 0], outer = (1, n ÷ 4)))
    x = [x1 x2]
    y = [ones(1, n ÷ 2) zeros(1, n ÷ 2)
         zeros(1, n ÷ 2) ones(1, n ÷ 2)]

    rp = randperm(n)

    xtest = x[:, rp[n÷2:end]]
    ytest = y[:, rp[n÷2:end]]
    x = x[:, rp[1:n÷2]]
    y = y[:, rp[1:n÷2]]

    ix = y[1, :] .== 1
    scatter(x[1, ix], x[2, ix])
    ix = y[2, :] .== 1
    scatter!(x[1, ix], x[2, ix])
    savefig("t.pdf")

    α = .5
    κ = 1e-4
    @debug "par" κ
end


text:
println(raw"                ______ _   _______ _   __")
println(raw"                | ___ \ | / /_   _| | / /   Version  2.2")
println(raw"                | |_/ / |/ /  | | | |/ /")
println(raw"                |    /|    \  | | |    \   David K. Zhang")
println(raw"                | |\ \| |\  \ | | | |\  \     (c) 2019")
println(raw"                |_| \_\_| \_/ |_| |_| \_/")
println()
println("RKTK is free software distributed under the terms of the MIT license.")
println()
flush(stdout)

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

using Base.Threads: @threads, nthreads, threadid
using Printf: @sprintf
using Statistics: mean, std
using UUIDs: UUID, uuid4

using DZLinearAlgebra
using DZOptimization
using MultiFloats
using RungeKuttaToolKit

push!(LOAD_PATH, @__DIR__)
using DZMisc

set_zero_subnormals(true)
use_standard_multifloat_arithmetic()

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

function precision_type(prec::Int)::Type
    if     prec <= 32;  Float32
    elseif prec <= 64;  Float64
    elseif prec <= 128; Float64x2
    elseif prec <= 192; Float64x3
    elseif prec <= 256; Float64x4
    elseif prec <= 320; Float64x5
    elseif prec <= 384; Float64x6
    elseif prec <= 448; Float64x7
    elseif prec <= 512; Float64x8
    else
        setprecision(prec)
        BigFloat
    end
end

approx_precision(::Type{Float32  }) = 32
approx_precision(::Type{Float64  }) = 64
approx_precision(::Type{Float64x2}) = 128
approx_precision(::Type{Float64x3}) = 192
approx_precision(::Type{Float64x4}) = 256
approx_precision(::Type{Float64x5}) = 320
approx_precision(::Type{Float64x6}) = 384
approx_precision(::Type{Float64x7}) = 448
approx_precision(::Type{Float64x8}) = 512
approx_precision(::Type{BigFloat }) = precision(BigFloat)

function Base.show(io::IO, ::Type{MultiFloat{Float64,N}}) where {N}
    write(io, "Float64x")
    show(io, N)
end

function Base.show(io::IO, ::Type{BigFloat})
    write(io, "BigFloat(")
    show(io, precision(BigFloat))
    write(io, ')')
end

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

RKOCBFGSOptimizer{T} = BFGSOptimizer{
    RKOCExplicitBackpropObjectiveFunctor{T},
    RKOCExplicitBackpropGradientFunctor{T}, T}

struct RKTKID
    order::Int
    num_stages::Int
    uuid::UUID
end

function Base.show(io::IO, id::RKTKID)
    write(io, "RKTK-")
    write(io, lpad(id.order, 2, '0'))
    write(io, lpad(id.num_stages, 2, '0'))
    write(io, '-')
    write(io, uppercase(string(id.uuid)))
end

const RKTKID_REGEX = Regex(
    "RKTK-([0-9]{2})([0-9]{2})-([0-9A-Fa-f]{8}-[0-9A-Fa-f]{4}-" *
    "[0-9A-Fa-f]{4}-[0-9A-Fa-f]{4}-[0-9A-Fa-f]{12})")
const RKTK_FILENAME_REGEX = Regex(
    "^[0-9]{4}-[0-9]{4}-[0-9]{4}-RKTK-([0-9]{2})([0-9]{2})-([0-9A-Fa-f]{8}-" *
    "[0-9A-Fa-f]{4}-[0-9A-Fa-f]{4}-[0-9A-Fa-f]{4}-[0-9A-Fa-f]{12})\\.txt\$")

function find_rktkid(str::String)::Union{RKTKID,Nothing}
    m = match(RKTKID_REGEX, str)
    if m != nothing
        RKTKID(parse(Int, m[1]), parse(Int, m[2]), UUID(m[3]))
    else
        nothing
    end
end

function find_filename_by_id(dir::String, id::RKTKID)::Union{String,Nothing}
    result = nothing
    for filename in readdir(dir)
        m = match(RKTKID_REGEX, filename)
        if ((m != nothing) && (parse(Int, m[1]) == id.order) &&
                (parse(Int, m[2]) == id.num_stages) && (UUID(m[3]) == id.uuid))
            if result != nothing
                say("ERROR: Found multiple files with RKTK ID $id.")
                exit()
            else
                result = filename
            end
        end
    end
    result
end

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

function log_score(x)::Int
    bx = BigFloat(x)
    if     iszero(bx)  ; typemax(Int)
    elseif isfinite(bx); round(Int, -100 * log10(bx))
    else               ; 0
    end
end

function scaled_log_score(x)
    bx = BigFloat(x)
    if     iszero(bx)  ; typemax(Int)
    elseif isfinite(bx); round(Int, 10000 - 3333 * log10(bx) / 2)
    else               ; 0
    end
end

score_str(x)::String = lpad(clamp(log_score(x), 0, 9999), 4, '0')
scaled_score_str(x)::String = lpad(clamp(scaled_log_score(x), 0, 9999), 4, '0')

rktk_score_str(opt::RKOCBFGSOptimizer{T}) where {T <: Real} =
    score_str(opt.current_objective_value[]) * '-' * score_str(norm(opt.current_gradient)) *
        '-' * scaled_score_str(norm(opt.current_point))

rktk_filename(opt, id::RKTKID)::String =
    rktk_score_str(opt) * '-' * string(id) * ".txt"

numstr(x)::String = @sprintf("%#-18.12g", BigFloat(x))
shortstr(x)::String = @sprintf("%#.5g", BigFloat(x))

function print_help()::Nothing
    say("Usage: julia RKTK.jl <command> [parameters...]")
    say()
    say("RKTK provides the following <command> options:")
    say("    search <order> <num-stages> <precision>")
    say("    refine <rktk-id> <precision>")
    say()
end

function print_table_header()::Nothing
    say(" Iteration │  Objective value  │   Gradient norm   │",
                    "  Last step size   │    Point norm     │ Type")
    say("───────────┼───────────────────┼───────────────────┼",
                    "───────────────────┼───────────────────┼──────")
end

function print_table_row(iter, obj_value, grad_norm,
                         step_size, point_norm, type)::Nothing
    say(" ", lpad(iter, 9, ' '), " | ",
        numstr(obj_value), "│ ", numstr(grad_norm),  "│ ",
        numstr(step_size), "│ ", numstr(point_norm), "│ ", type)
end

function rmk_table_row(iter, obj_value, grad_norm,
                       step_size, point_norm, type)::Nothing
    rmk(" ", lpad(iter, 9, ' '), " | ",
        numstr(obj_value), "│ ", numstr(grad_norm),  "│ ",
        numstr(step_size), "│ ", numstr(point_norm), "│ ", type)
end

function print_table_row(opt, type)::Nothing
    print_table_row(opt.iteration_count[], opt.current_objective_value[], norm(opt.current_gradient),
                    opt.last_step_size[1], norm(opt.current_point), type)
end

function rmk_table_row(opt, type)::Nothing
    rmk_table_row(opt.iteration_count[], opt.current_objective_value[], norm(opt.current_gradient),
                  opt.last_step_size[1], norm(opt.current_point), type)
end

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

function rkoc_optimizer(::Type{T}, order::Int, num_stages::Int,
        x_init::Vector{BigFloat}, num_iters::Int) where {T <: Real}
    obj_func, grad_func = rkoc_explicit_backprop_functors(T, order, num_stages)
    num_vars = div(num_stages * (num_stages + 1), 2)
    @assert length(x_init) == num_vars
    opt = BFGSOptimizer(obj_func, grad_func, T.(x_init), inv(T(1_000_000)))
    opt.iteration_count[] = num_iters
    opt
end

function rkoc_optimizer(::Type{T}, id::RKTKID,
                        filename::String) where {T <: Real}
    trajectory = filter(!isempty, strip.(split(read(filename, String), "\n\n")))
    point_data = filter(!isempty, strip.(split(trajectory[end], '\n')))
    header = split(point_data[1])
    x_init = BigFloat.(point_data[2:end])
    num_iters = parse(Int, header[1])
    rkoc_optimizer(T, id.order, id.num_stages, x_init, num_iters), header
end

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

function save_to_file(opt::RKOCBFGSOptimizer{T}, id::RKTKID) where {T <: Real}
    filename = rktk_filename(opt, id)
    rmk("Saving progress to file ", filename, "...")
    old_filename = find_filename_by_id(".", id)
    if old_filename != nothing
        if old_filename != filename
            mv(old_filename, filename)
        end
        trajectory = filter(!isempty,
            strip.(split(read(filename, String), "\n\n")))
        point_data = filter(!isempty, strip.(split(trajectory[end], '\n')))
        header = split(point_data[1])
        num_iters = parse(Int, header[1])
        prec = parse(Int, header[2])
        if (precision(BigFloat) <= prec) && (opt.iteration_count[] <= num_iters)
            rmk("No save necessary.")
            return
        end
    end
    file = open(filename, "a+")
    println(file, opt.iteration_count[], ' ', precision(BigFloat), ' ',
            log_score(opt.current_objective_value[]), ' ', log_score(norm(opt.current_gradient)))
    for x in opt.current_point
        println(file, BigFloat(x))
    end
    println(file)
    close(file)
    rmk("Save complete.")
end

const TERM = isa(stdout, Base.TTY)

function run!(opt::RKOCBFGSOptimizer{T}, id::RKTKID) where {T <: Real}
    print_table_header()
    print_table_row(opt, "NONE")
    save_to_file(opt, id)
    last_print_time = last_save_time = time_ns()
    while true
        step!(opt)
        if opt.has_converged[]
            print_table_row(opt, "DONE")
            save_to_file(opt, id)
            return
        end
        current_time = time_ns()
        if current_time - last_save_time > UInt(60_000_000_000)
            save_to_file(opt, id)
            last_save_time = current_time
        end
        if opt.last_step_type[] == DZOptimization.GradientDescentStep
            print_table_row(opt, "GRAD")
            last_print_time = current_time
        elseif TERM && (current_time - last_print_time > UInt(80_000_000))
            rmk_table_row(opt, "BFGS")
            last_print_time = current_time
        end
    end
end

function run!(opt::RKOCBFGSOptimizer{T}, id::RKTKID,
              duration_ns::UInt) where {T <: Real}
    save_to_file(opt, id)
    start_time = last_save_time = time_ns()
    while true
        _, objective_decreased = step!(opt)
        current_time = time_ns()
        if (!objective_decreased) || (current_time - start_time > duration_ns)
            save_to_file(opt, id)
            return !objective_decreased
        end
        if current_time - last_save_time > UInt(60_000_000_000)
            save_to_file(opt, id)
            last_save_time = current_time
        end
    end
end

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

function search(::Type{T}, id::RKTKID) where {T <: Real}
    setprecision(approx_precision(T))
    num_vars = div(id.num_stages * (id.num_stages + 1), 2)
    optimizer = rkoc_optimizer(T, id.order, id.num_stages,
        rand(BigFloat, num_vars), 0)
    say("Running $T search $id.\n")
    run!(optimizer, id)
    say("\nCompleted $T search $id.\n")
end

# function multisearch(::Type{T}, order::Int, num_stages::Int) where {T <: Real}
#     num_threads = nthreads()
#     num_opts = num_threads - 1
#     num_vars = div(num_stages * (num_stages + 1), 2)
#     say("Constructing optimizers...")
#     optimizers = [rkoc_optimizer(
#             T, order, num_stages, rand(BigFloat, num_vars), 0)
#         for _ = 1 : num_opts]
#     bfgs_used = zeros(Bool, num_opts)
#     objective_decreased = zeros(Bool, num_opts)
#     for i = 1 : num_opts
#         print_table_row(optimizers[i], "NONE")
#         bfgs_used[i], objective_decreased[i] = step!(optimizers[i])
#     end
#     @threads for i = 1 : num_threads
#         if i <= num_opts
#             while true
#                 bfgs_used[i], objective_decreased[i] = step!(optimizers[i])
#             end
#         else
#             while true
#                 for i = 1 : num_opts
#                     print("\033[F")
#                 end
#                 for i = 1 : num_opts
#                     print_table_row(optimizers[i],
#                         ifelse(objective_decreased[i],
#                             ifelse(bfgs_used[i], "BFGS", "GRAD"), "DONE"))
#                 end

#             end
#         end
#     end
# end

function refine(::Type{T}, id::RKTKID, filename::String) where {T <: Real}
    setprecision(approx_precision(T))
    say("Running ", T, " refinement $id.\n")
    optimizer, header = rkoc_optimizer(T, id, filename)
    if precision(BigFloat) < parse(Int, header[2])
        say("WARNING: Refining at lower precision than source file.\n")
    end
    starting_iteration = optimizer.iteration_count[]
    run!(optimizer, id)
    ending_iteration = optimizer.iteration_count[]
    if ending_iteration > starting_iteration
        say("\nRepeating $T refinement $id.\n")
        refine(T, id, find_filename_by_id(".", id))
    else
        say("\nCompleted $T refinement $id.\n")
    end
end

function clean(::Type{T}) where {T <: Real}
    setprecision(approx_precision(T))
    optimizers = Tuple{Int,RKTKID,RKOCBFGSOptimizer{T}}[]
    for filename in readdir()
        if match(RKTKID_REGEX, filename) != nothing
            id = find_rktkid(filename)
            optimizer, header = rkoc_optimizer(T, id, filename)
            if precision(BigFloat) < parse(Int, header[2])
                say("ERROR: Cleaning at lower precision than source file \"",
                    filename, "\".")
                exit()
            end
            push!(optimizers,
                (log_score(optimizer.current_objective_value[]), id, optimizer))
        end
    end
    say("Found ", length(optimizers), " RKTK files.")
    while true
        num_optimizers = length(optimizers)
        sort!(optimizers, by=(t -> t[1]), rev=true)
        completed = zeros(Bool, num_optimizers)
        @threads for i = 1 : num_optimizers
            _, id, optimizer = optimizers[i]
            old_score = rktk_score_str(optimizer)
            start_iter = optimizer.iteration_count[]
            completed[i] = run!(optimizer, id, UInt(10_000_000_000))
            stop_iter = optimizer.iteration_count[]
            new_score = rktk_score_str(optimizer)
            say(ifelse(completed[i], "    Cleaned ", "    Working "),
                id, " (", stop_iter - start_iter, " iterations: ",
                old_score, " => ", new_score, ") on thread ", threadid(), ".")
        end
        next_optimizers = Tuple{Int,RKTKID,RKOCBFGSOptimizer{T}}[]
        for i = 1 : length(optimizers)
            if !completed[i]
                _, id, optimizer = optimizers[i]
                push!(next_optimizers,
                    (log_score(optimizer.current_objective_value[]), id, optimizer))
            end
        end
        if length(next_optimizers) == 0
            say("All RKTK files cleaned!")
            break
        end
        optimizers = next_optimizers
        say(length(optimizers), " RKTK files remaining.")
    end
end

function benchmark(::Type{T}, order::Int, num_stages::Int,
        benchmark_secs) where {T <: Real}
    setprecision(approx_precision(T))
    num_vars = div(num_stages * (num_stages + 1), 2)
    x_init = [BigFloat(i) / num_vars for i = 1 : num_vars]
    construction_ns = time_ns()
    opt = rkoc_optimizer(T, order, num_stages, x_init, 0)
    start_ns = time_ns()
    benchmark_ns = round(typeof(start_ns), benchmark_secs * 1_000_000_000)
    terminated_early = false
    while time_ns() - start_ns < benchmark_ns
        _, made_progress = step!(opt)
        if !made_progress
            terminated_early = true
            break
        end
    end
    Int(start_ns - construction_ns), opt.iteration_count[], terminated_early
end

function benchmark(::Type{T}, order::Int, num_stages::Int,
        benchmark_secs, num_trials::Int) where {T <: Real}
    construction_secs = Float64[]
    iteration_counts = Int[]
    success = true
    for _ = 1 : num_trials
        construction_ns, iteration_count, terminated_early =
            benchmark(T, order, num_stages, benchmark_secs / num_trials)
        push!(construction_secs, construction_ns / 1_000_000_000)
        push!(iteration_counts, iteration_count)
        if terminated_early
            success = false
            break
        end
    end
    if success
        say(rpad("$T: ", 16, ' '),
            shortstr(mean(iteration_counts)), " ± ",
            shortstr(std(iteration_counts)))
    else
        say(rpad("$T: ", 16, ' '), "Search terminated too early")
    end
end

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

function get_order(n::Int)
    result = tryparse(Int, ARGS[n])
    if (result == nothing) || (result < 1) || (result > 20)
        say("ERROR: Parameter $n (\"$(ARGS[n])\") must be ",
            "an integer between 1 and 20.")
        exit()
    end
    result
end

function get_num_stages(n::Int)
    result = tryparse(Int, ARGS[n])
    if (result == nothing) || (result < 1) || (result > 99)
        say("ERROR: Stage parameter $n (\"$(ARGS[n])\") must be ",
            "an integer between 1 and 99.")
        exit()
    end
    result
end

function main()

    if (length(ARGS) == 0) || ("-h" in ARGS) || ("--help" in ARGS) let
        print_help()

    end elseif uppercase(ARGS[1]) == "SEARCH" let
        order, num_stages = get_order(2), get_num_stages(3)
        prec = parse(Int, ARGS[4])
        while true
            search(precision_type(prec), RKTKID(order, num_stages, uuid4()))
        end

    # end elseif uppercase(ARGS[1]) == "MULTISEARCH" let
    #     order, num_stages = get_order(2), get_num_stages(3)
    #     prec = parse(Int, ARGS[4])
    #     multisearch(precision_type(prec), order, num_stages)

    end elseif uppercase(ARGS[1]) == "REFINE" let
        id = find_rktkid(ARGS[2])
        if id == nothing
            say("ERROR: Invalid RKTK ID ", ARGS[2], ".")
            say("RKTK IDs have the form ",
                "RKTK-XXYY-ZZZZZZZZ-ZZZZ-ZZZZ-ZZZZ-ZZZZZZZZZZZZ.\n")
            exit()
        end
        filename = find_filename_by_id(".", id)
        if filename == nothing
            say("ERROR: No file exists with RKTK ID $id.\n")
            exit()
        end
        prec = parse(Int, ARGS[3])
        refine(precision_type(prec), id, filename)

    end elseif uppercase(ARGS[1]) == "CLEAN" let
        prec = parse(Int, ARGS[2])
        clean(precision_type(prec))

    end elseif uppercase(ARGS[1]) == "BENCHMARK" let
        order, num_stages = get_order(2), get_num_stages(3)
        benchmark_secs = parse(Float64, ARGS[4])
        num_trials = parse(Int, ARGS[5])
        for T in (Float32, Float64, Float64x2, Float64x3, Float64x4,
                    Float64x5, Float64x6, Float64x7, Float64x8)
            benchmark(T, order, num_stages, benchmark_secs, num_trials)
        end
        for T in (Float32, Float64, Float64x2, Float64x3, Float64x4,
                    Float64x5, Float64x6, Float64x7, Float64x8)
            setprecision(approx_precision(T))
            benchmark(BigFloat, order, num_stages, benchmark_secs, num_trials)
        end
        say()

    end else let
        say("ERROR: Unrecognized <command> option \"", ARGS[1], "\".\n")
        print_help()

    end end
end

main()

ext: jl


text:
"""
    Utils

Module containing various utility functions for `WaveProp`.
"""
module Utils

using DocStringExtensions
using StaticArrays

using WaveProp

export
    svector,
    matrix_to_blockmatrix,
    blockmatrix_to_matrix,
    blockvector_to_vector,
    vector_to_blockvector,
    notimplemented,
    abstractmethod,
    assert_extension,
    assert_concrete_type

"""
    svector(f,n)

Just like [`Base.ntuple`](https://docs.julialang.org/en/v1/base/base/#Base.ntuple), but convert output to an `SVector`.
"""
svector(f,n) = ntuple(f,n) |> SVector

"""
    blockmatrix_to_matrix(A::Matrix{B}) where {B<:SMatrix}

Convert a `Matrix{B}`, where `B<:SMatrix`, to the equivalent `Matrix{T}`, where `T = eltype(B)`
"""
function blockmatrix_to_matrix(A::Matrix{B})  where B <: SMatrix
    T = eltype(B)
    sblock = size(B)
    ss     = size(A).*sblock # matrix size when viewed as matrix over T
    Afull = Matrix{T}(undef,ss)
    for i=1:ss[1], j=1:ss[2]
        bi, ind_i = divrem(i-1,sblock[1]) .+ (1,1)
        bj, ind_j = divrem(j-1,sblock[2]) .+ (1,1)
        Afull[i,j] = A[bi,bj][ind_i,ind_j]
    end
    return Afull
end

"""
    blockvector_to_vector(A::Vector{B}) where {B<:SVector}

Convert a `Vector{B}`, where `B<:SVector`, to the equivalent `Vector{T}`, where `T = eltype(B)`
"""
function blockvector_to_vector(A::Vector{B})  where B <: SVector
    T = eltype(B)
    reinterpret(T,A) |> collect
end

"""
    matrix_to_blockmatrix(A::Matrix,B)

Convert a `Matrix{T}` to a `Matrix{B}`, where `B<:Type{SMatrix}`. The element
type of `B` must match that of `A`, and the size of `A` must be divisible by the
size of `B` along each dimension.
"""
function matrix_to_blockmatrix(A::Matrix,B::Type{<:SMatrix})
    @assert eltype(A) == eltype(B)
    @assert sum(size(A) .% size(B)) == 0 "block size $(size(B)) not compatible with size of A=$(size(A))"
    sblock = size(B)
    nblock = div.(size(A),sblock)
    Ablock = Matrix{B}(undef,nblock)
    for i in 1:nblock[1]
        istart = (i-1)*sblock[1] + 1
        iend = i*sblock[1]
        for j in 1:nblock[2]
            jstart = (j-1)*sblock[2] + 1
            jend   = j*sblock[2]
            Ablock[i,j] = A[istart:iend,jstart:jend]
        end
    end
    return Ablock
end

"""
    vector_to_blockvector(A::Vector,B)

Convert a `Vector{T}` to a `Vector{B}`, where `B<:Type{SVector}`. The element
type of `B` must match that of `A`, and the size of `A` must be divisible by the
size of `B` along each dimension.
"""
function vector_to_blockvector(A::Vector,B::Type{<:SVector})
    @assert eltype(A) == eltype(B)
    @assert sum(size(A) .% size(B)) == 0 "block size $(size(B)) not compatible with size of A=$(size(A))"
    T = eltype(B)
    reinterpret(B,A) |> collect
end

"""
    notimplemented()

Things which should probably be implemented at some point.
"""
function notimplemented()
    error("not (yet) implemented")
end

"""
    abstractmethod

A method of an `abstract type` for which concrete subtypes are expected
to provide an implementation.
"""
function abstractmethod(T)
    error("this method needs to be implemented by the concrete subtype $T.")
end

"""
    assert_extension(fname,ext,[msg])

Check that `fname` is of extension `ext`. Print the message `msg` as an assertion error otherwise.
"""
function assert_extension(fname::String,ext::String,msg="file extension must be $(ext)")
    r = Regex("$(ext)\$")
    @assert occursin(r,fname) msg
end

function assert_concrete_type(T::DataType)
    isconcretetype(T) || throw(ConcreteInferenceError(T))
end

function debug(mod="WaveProp")
    @eval ENV["JULIA_DEBUG"] = $(mod)
end

struct ConcreteInferenceError <: Exception
    T
end
Base.showerror(io::IO, e::ConcreteInferenceError) = print(io, "unable to infer concrete type from function signature: T=$(e.T)" )

end # module

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
mutable struct TAPGraphIsing{TJ, TH}
    N::Int
    adjlist::Vector{Vector{Int}}
    J::Vector{Vector{TH}}
    H::Vector{TJ}
    m::Vector{Float64}
    m_old::Vector{Float64}
end

mag(g::TAPGraphIsing, i::Integer) = g.m[i]
mags(g::TAPGraphIsing) = g.m
deg(g::TAPGraphIsing, i::Integer) = length(g.adjlist[i])

function TAPGraphIsing(net::Network; T=1)
    @assert has_eprop(net, "J")
    @assert has_vprop(net, "H")

    adjlist = adjacency_list(net)
    N = nv(net)
    J = [[eprop(net, e)["J"] / T for e in edges(net, i)] for i=1:N]
    H = [vprop(net, i)["H"] / T for i=1:N]

    m = zeros(N)
    m_old = zeros(N)
    TAPGraphIsing(N, adjlist, J, H, m, m_old)
end

function initrand!(g::TAPGraphIsing; μ=0, σ=1)
    for i=1:g.N
        g.m[i] = μ + σ * randn()
        g.m_old[i] = g.m[i]
    end
end

function oneTAPiter!(g::TAPGraphIsing; dump=0.1)
    @unpack N, J, H, m, m_old = g
    Δ = 0.
    mnew = zeros(N)
    for i=1:N
        h = H[i]
        for (ki, j) in enumerate(g.adjlist[i])
            h += J[i][ki]*m[j] - J[i][ki]^2 * (1-m[j]^2) * m_old[i]
        end
        mnew[i] = tanh(h)
    end
    for i=1:N
        m_old[i] = m[i]
        m[i] = (1-dump) * mnew[i] + dump * m_old[i]
        Δ = max(Δ, abs(m[i] - m_old[i]))
    end
    return Δ
end

function converge!(g::TAPGraphIsing; maxiters::Int=1000,
                                    ϵ::Float64=1e-6,
                                    verbose::Bool=false,
                                    dump=0)
    for it=1:maxiters
        verbose && print("it=$it ... ")
        Δ = oneTAPiter!(g; dump)
        verbose && @printf("Δ=%.2g \n", Δ)
        if Δ < ϵ
            verbose && println("Converged!")
            break
        end
    end
end

function corr_disc_nn(g::TAPGraphIsing, i::Int, j::Int)
    @unpack N, adjlist, m = g
    ki = findfirst(adjlist[i], j)
    kj = findfirst(adjlist[j], i)
    @assert(ki > 0)
    @assert(kj > 0)

    J = g.J[i][ki]
    mi = mag(g, i)
    mj = mag(g, j)

    c = J * (1-mi^2)*(1-mj^2) / (1+ J * mi * mj) + mi*mj
    return c
end

corr_conn_nn(g::TAPGraphIsing, i::Int, j::Int) = corr_disc_nn(g, i, j) - mag(g, i)*mag(g,j)

function corr_disc_nn(g::TAPGraphIsing)
    corrs = Vector{Vector{Float64}}(g.N)
    for i=1:g.N
        corrs[i] = zeros(length(g.adjlist[i]))
    end

    for i=1:g.N
        for (ki,j) in enumerate(g.adjlist[i])
            corrs[i][ki] = corr_disc_nn(g,i,j)
        end
    end
    corrs
end

function run_tap(net::Network; T=1,
                            maxiters=1000,
                            ϵ=1e-6,
                            verbose=true,
                            dump=0)
    g = TAPGraphIsing(net; T)
    initrand!(g, μ=0, σ=1)
    converge!(g; maxiters, ϵ, verbose, dump)
    return g
end


text:
#!/usr/bin/env julia

using Luxor

function get_png_files(folder)
    cd(folder)
    imagelist = filter(f -> !startswith(f, ".") && endswith(f, "png"), readdir(folder))
    return map(realpath, imagelist)
end

imagelist = get_png_files(dirname(@__FILE__) * "/../docs/figures")

width, height = 2000, 2000
fname = "/tmp/paint-images.pdf"
Drawing(width, height, fname)
origin()
background("grey50")

for imgfile in imagelist
    img = readpng(imgfile)
    rand(Bool) ? placeimage(img, rand(-width/2:width/2), rand(-height/2:height/2), rand(0.25:0.1:0.75)) : placeimage(img, Point(rand(-width/2:width/2), rand(-height/2:height/2)), rand(0.25:0.1:0.75))
end

finish()
println("finished test: output in $(fname)")


text:
using Knet
include("../helper.jl")
include("accuracy.jl")
include("loss.jl")
include("mpii.batcher.jl")
include("../models.jl")
include("../training.jl")
include("../plots.jl")
include("../modular.resnet.jl")

test_image_count = 1000
validation_image_count =  28883 - 25600 - 1000
train_image_count = 25600
read_image_w = 512
read_image_h = 512
pre_full_path = "/kuacc/users/gsoykan20/comp541_term_project/deeper-cut/"

function load_and_save_batches(save_tag; should_use_pmap = true)
    dtrn, dval, dtst, data_items = get_mpii_batches_and_data_items(1; should_shuffle=true, should_use_pmap = should_use_pmap);
trn_data_items = data_items[begin : end - validation_image_count - test_image_count];
val_data_items = data_items[end - validation_image_count - test_image_count + 1 : end - test_image_count]
tst_data_items = data_items[end - test_image_count + 1 : end];

    println("start saving dtrn")
    Knet.save("$(pre_full_path)$(save_tag)-dtrn.jld2", "$(save_tag)-dtrn", dtrn)

    println("start saving dval")
    Knet.save("$(pre_full_path)$(save_tag)-dval.jld2", "$(save_tag)-dval", dval)

    println("start saving dtst")
    Knet.save("$(pre_full_path)$(save_tag)-dtst.jld2", "$(save_tag)-dtst", dtst)

    println("start saving trn data items")
    Knet.save("$(pre_full_path)$(save_tag)-trn_data_items.jld2", "$(save_tag)-trn_data_items", trn_data_items)

    println("start saving val data items")
    Knet.save("$(pre_full_path)$(save_tag)-val_data_items.jld2", "$(save_tag)-val_data_items", val_data_items)

    println("start saving tst data items")
    Knet.save("$(pre_full_path)$(save_tag)-tst_data_items.jld2", "$(save_tag)-tst_data_items", tst_data_items)

    println("save completed!")
    return (dtrn, dval, dtst, trn_data_items, val_data_items, tst_data_items)
end


text:
using CSV, GLM, Plots, Random; pyplot()
Random.seed!(0)

df = CSV.read("../data/weightHeight.csv", copycols = true)
n = size(df)[1]
df[shuffle(1:n),:] = df
df[[10,40,60,130,140,175,190,200],:Sex] .= "O"

conts = Dict(:Sex=>DummyCoding(base="F",levels=["F","M","O"]))
model1 = lm(@formula(Height ~ Weight * Sex), df, contrasts=conts)
model2 = lm(@formula(Height ~ Weight & Sex), df, contrasts=conts)
model3 = lm(@formula(Height ~ Weight + Weight & Sex), df, contrasts=conts)

a0, a1, a2, a3, a4, a5 = coef(model1)
b0, b1, b2, b3 = coef(model2)
c0, c1, c2, c3 = coef(model3)

println(model1)
println(model2)
println(model3)
println("Model2 and Model3 equivalence: ",
    round.((b0 - c0, b1 - c1, b2 - (c1+c2), b3 - (c1+c3)),digits=5))

pred1(weight,sex) =  a0 + a1 * weight +
				(sex=="M")*(a2+a4*weight)+
				(sex=="O")*(a3+a5*weight)
pred2(weight,sex) =  b0 + weight*(b1*(sex=="F") + b2*(sex=="M") + b3*(sex=="O"))

males,females,other=df[df.Sex .=="M",:],df[df.Sex .=="F",:],df[df.Sex .=="O",:]

xlim = [0, maximum(df.Weight)]
scatter(males.Weight, males.Height, c=:blue, msw=0, label="Male")
plot!(xlim, pred1.(xlim,"M"), c=:blue, lw = 1.5, label="Male fit 1")
plot!(xlim, pred2.(xlim,"M"), c=:blue, linestyle=:dash, label="Male fit 2")
scatter!(females.Weight, females.Height, c=:red, msw=0, label="Female")
plot!(xlim, pred1.(xlim,"F"), c=:red, lw=1.5,label="Female fit 1",xlims=(xlim))
plot!(xlim, pred2.(xlim,"F"), c=:red, linestyle=:dash, label="Female fit 2")
scatter!(other.Weight, other.Height, c=:green, msw=0, label="Other",
	xlabel="Weight (kg)", ylabel="Height (cm)", legend=:topleft)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
module KrylovMethods

	using LinearAlgebra
	using SparseArrays
	using Printf

	import LinearAlgebra.BLAS

	include("cg.jl")
	include("blockCG.jl")
	include("cgls.jl")
	include("bicgstb.jl")
	include("blockBiCGSTB.jl")
	include("gmres.jl")
	include("fgmres.jl")
	include("blockFGMRES.jl")
	include("lanczosBidiag.jl")
	include("ssor.jl")
	include("lsqr.jl")
	include("lanczosTridiag.jl")
	include("lanczos.jl")
	include("minres.jl")
	include("gs.jl")
end


text:
#=
This module doesn't work until now for PyCall works ill.
However if you include it in Main, everything is okay.
=#

module Process
export to_dict, process
using PyCall

function to_dict(py_obj :: PyObject) :: Any

    pisa = pybuiltin("isinstance")
    phas = pybuiltin("hasattr")
    ast = pyimport("ast")
    AST = ast[:AST]
    if pisa(py_obj, AST)

        ty     = py_obj[:__class__][:__name__]
        fields = py_obj[:_fields] |> collect
        map(fields) do field
            field = Symbol(field)
            value = py_obj[field]
            field =>
                if value isa Vector || value isa Tuple
                    map(to_dict, value)
                else
                    value |> to_dict
                end
        end |>
        function (iteritems)
            dict = Dict{Symbol, Any}(iteritems)
            if phas(py_obj, "lineno")
                dict[:lineno] = py_obj[:lineno]
            end
            if phas(py_obj, "col_offset")
                dict[:colno] = py_obj[:col_offset]
            end
            dict[:class] = ty
            dict
        end
    else
        @error "Invalid pyobj. Expected $(string(AST)), got $(py_obj[:__class__])."
    end
end

to_dict(num :: Number) = num
to_dict(str :: String) = str
to_dict(::Nothing) = nothing

function process(codes:: String)

    ast = pyimport("ast")
    node = ast[:parse](codes)
    to_dict(node)
end

end


text:
requiredTypes = [ "CareerPath" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( dirname( Base.source_path() ), "..", "types", reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end

function ActivateLevel(careerpath::CareerPath, level::DataType)
    push!(careerpath.ActivatedLevels, LevelsCodes[level])
end

function isActivatedLevel(careerpath::CareerPath, level::DataType)
    return LevelsCodes[level] in careerpath.ActivatedLevels
end

function addStateToCareerPath(cp::CareerPath, st::CareerStatus)
    push!(cp.Path, st)
end

function set_CareerPath_Length(cp::CareerPath, length::Int)
    cp.Length = length
end

function compute_Length(cp::CareerPath)
    if cp.Length == 0
        for i in 1:length(cp.Path)
            cp.Length += cp.Path[i].MinStay
        end
    end
end

function get_status_index_after_duration(cp::CareerPath, duration::Int)

    spentleng = 0
    index = -1
    for i in 1:length(cp.Path)
        spentleng += cp.Path[i].MinStay
        if spentleng >= duration
            index = i
            break
        end
    end
    return index
end

function AttritionAtYear(cp::CareerPath, year::Int)

    currentyear = 0
    attrition = 0
    index = 1
    sumprenode = 0
    while true
        currentyear += 1

        if currentyear >= year
            break
        else
            if currentyear > sumprenode + cp.Path[index].MinStay
                sumprenode += cp.Path[index].MinStay
                index += 1
            end
            attrition += cp.Path[index].Attrition
        end
    end
    return attrition
end

max_stars_repo_name: Omazaria/CareerPathBeDef

max_issues_repo_name: Omazaria/CareerPathBeDef

max_forks_repo_name: Omazaria/CareerPathBeDef


text:
using SparseArrays, LuxurySparse
import YaoBase: diff, autostatic

"""
    cunmat(nbit::Int, cbits::NTuple{C, Int}, cvals::NTuple{C, Int}, U0::AbstractMatrix, locs::NTuple{M, Int}) where {C, M} -> AbstractMatrix

control-unitary matrix
"""
function cunmat end

"""
    u1ij!(target, i, j, a, b, c, d)
single u1 matrix into a target matrix.

!!! note
    For coo, we take an additional parameter
        * ptr: starting position to store new data.
"""
function u1ij! end

"""
    setcol!(csc::SparseMatrixCSC, icol::Int, rowval::AbstractVector, nzval) -> SparseMatrixCSC

set specific col of a CSC matrix
"""
@inline function setcol!(csc::SparseMatrixCSC, icol::Int, rowval::AbstractVector, nzval)
    @inbounds begin
        S = csc.colptr[icol]
        E = csc.colptr[icol+1] - 1
        csc.rowval[S:E] = rowval
        csc.nzval[S:E] = nzval
    end
    csc
end

"""
    num_nonzero(nbits, nctrls, U, [N])

Return number of nonzero entries of the matrix form of control-U gate. `nbits`
is the number of qubits, and `nctrls` is the number of control qubits.
"""
@inline function num_nonzero(nbits::Int, nctrls::Int, U, N::Int = 1 << nbits)
    return N + (1 << (nbits - nctrls - log2dim1(U))) * (length(U) - size(U, 2))
end

@inline function num_nonzero(
    nbits::Int,
    nctrls::Int,
    U::SDSparseMatrixCSC,
    N::Int = 1 << nbits,
)
    return N + (1 << (nbits - nctrls - log2dim1(U))) * (nnz(U) - size(U, 2))
end

"""
    getcol(csc::SDparseMatrixCSC, icol::Int) -> (View, View)

get specific col of a CSC matrix, returns a slice of (rowval, nzval)
"""
@inline function getcol(csc::SDSparseMatrixCSC, icol::Int)
    @inbounds begin
        S = csc.colptr[icol]
        E = csc.colptr[icol+1] - 1
        view(csc.rowval, S:E), view(csc.nzval, S:E)
    end
end

@inline function reorder_unitary(
    nbit::Int,
    cbits::NTuple{C,Int},
    cvals::NTuple{C,Int},
    U0::AbstractMatrix,
    locs::NTuple{M,Int},
) where {C,M}
    # reorder a unirary matrix.
    U = all(diff(locs) .> 0) ? U0 : reorder(U0, collect(locs) |> sortperm)
    locked_bits = [cbits..., locs...]
    locked_vals = [cvals..., zeros(Int, M)...]
    locs_raw =
        [i + 1 for i in itercontrol(nbit, setdiff(1:nbit, locs), zeros(Int, nbit - M))]
    ic = itercontrol(nbit, locked_bits, locked_vals)
    return U |> autostatic, ic, locs_raw |> staticize
end

adaptive_pow2(n::Int) = adaptive_pow2(UInt(n))

function adaptive_pow2(n::UInt)
    n < 62 ? Int64(1) << n : n < 126 ? Int128(1) << n : big(1) << n
end

const LARGE_MATRIX_WARN = 62

function large_mat_check(n::Int)
    if n > LARGE_MATRIX_WARN
        error(
            "matrix is too large, expect n <= $LARGE_MATRIX_WARN, got $n, integer overflows",
        )
    end
    return nothing
end

function cunmat(n::Int, cbits::NTuple, cvals::NTuple, U::IMatrix, locs::NTuple)
    large_mat_check(n)
    IMatrix{1 << n}()
end

"""
    unmat(nbit::Int, U::AbstractMatrix, locs::NTuple) -> AbstractMatrix

Return the matrix representation of putting matrix at locs.
"""
function unmat(nbit::Int, U::AbstractMatrix, locs::NTuple)
    large_mat_check(nbit)
    cunmat(nbit::Int, (), (), U, locs)
end

############################### Dense Matrices ###########################
function u1mat(nbit::Int, U1::AbstractMatrix, ibit::Int)
    large_mat_check(nbit)
    unmat(nbit, U1, (ibit,))
end

u1mat(nbit::Int, U1::Adjoint, ibit::Int) = u1mat(nbit, copy(U1), ibit)

function u1mat(nbits::Int, U1::SDMatrix, ibit::Int)
    large_mat_check(nbits)
    mask = bmask(ibit)
    N = 1 << nbits
    a, c, b, d = U1
    step = 1 << (ibit - 1)
    step_2 = 1 << ibit
    NNZ = num_nonzero(nbits, 0, U1, N)

    colptr = Vector{Int}(1:2:2*N+1)
    rowval = Vector{Int}(undef, NNZ)
    nzval = Vector{eltype(U1)}(undef, NNZ)

    mat = SparseMatrixCSC(N, N, colptr, rowval, nzval)
    for j = 0:step_2:N-step
        @inbounds @simd for i = j+1:j+step
            u1ij!(mat, i, i + step, a, b, c, d)
        end
    end
    return mat
end

@inline function u1ij!(csc::SparseMatrixCSC, i::Int, j::Int, a, b, c, d)
    @inbounds begin
        csc.rowval[2*i-1] = i
        csc.rowval[2*i] = j
        csc.rowval[2*j-1] = i
        csc.rowval[2*j] = j

        csc.nzval[2*i-1] = a
        csc.nzval[2*i] = c
        csc.nzval[2*j-1] = b
        csc.nzval[2*j] = d
    end
    csc
end

@inline function unij!(mat::SparseMatrixCSC, locs, U::SDMatrix)
    @simd for j = 1:size(U, 2)
        @inbounds setcol!(mat, locs[j], locs, view(U, :, j))
    end
    csc
end

function cunmat(
    nbit::Int,
    cbits::NTuple{C,Int},
    cvals::NTuple{C,Int},
    U0::Adjoint,
    locs::NTuple{M,Int},
) where {C,M}
    cunmat(nbit, cbits, cvals, copy(U0), locs)
end

function cunmat(
    nbit::Int,
    cbits::NTuple{C,Int},
    cvals::NTuple{C,Int},
    U0::SDMatrix,
    locs::NTuple{M,Int},
) where {C,M}
    large_mat_check(nbit)
    MM = size(U0, 1)
    U, ic, locs_raw = reorder_unitary(nbit, cbits, cvals, U0, locs)

    N = 1 << nbit
    NNZ = num_nonzero(nbit, C, U0, N)
    colptr = Vector{Int}(undef, N + 1)
    rowval = Vector{Int}(undef, NNZ)
    nzval = ones(eltype(U0), NNZ)

    @inbounds colptr[1] = 1
    ctest = controller(cbits, cvals)
    @inbounds @simd for b in basis(nbit)
        if ctest(b)
            colptr[b+2] = colptr[b+1] + MM
        else
            colptr[b+2] = colptr[b+1] + 1
            rowval[colptr[b+1]] = b + 1
        end
    end

    mat = SparseMatrixCSC(N, N, colptr, rowval, nzval)
    controldo(ic) do i
        unij!(mat, locs_raw .+ i, U)
    end
    return mat
end

############################### SparseMatrix ##############################
function cunmat(
    nbit::Int,
    cbits::NTuple{C,Int},
    cvals::NTuple{C,Int},
    U0::SparseMatrixCSC{Tv},
    locs::NTuple{M,Int},
)::SparseMatrixCSC{Tv} where {C,M,Tv}
    large_mat_check(nbit)
    N = 1 << nbit
    U, ic, locs_raw = reorder_unitary(nbit, cbits, cvals, U0, locs)

    NNZ = num_nonzero(nbit, C, U0, N)
    colptr = Vector{Int}(undef, N + 1)
    rowval = Vector{Int}(undef, NNZ)
    nzval = ones(eltype(U0), NNZ)
    @inbounds colptr[1] = 1

    ns = diff(U.colptr) |> autostatic
    Ns = ones(Int, N)
    controldo(ic) do i
        @inbounds Ns[locs_raw.+i] = ns
    end
    @inbounds @simd for j = 1:N
        colptr[j+1] = colptr[j] + Ns[j]
        if Ns[j] == 1
            rowval[colptr[j]] = j
        end
    end

    mat = SparseMatrixCSC(N, N, colptr, rowval, nzval)
    controldo(ic) do i
        unij!(mat, locs_raw .+ i, U)
    end
    mat
end

@inline function unij!(mat::SparseMatrixCSC, locs, U::SDSparseMatrixCSC)
    @simd for j = 1:size(U, 2)
        rows, vals = getcol(U, j)
        @inbounds setcol!(mat, locs[j], view(locs, rows), vals)
    end
    csc
end

############################# PermMatrix ###############################
@inline function unij!(pm::PermMatrix, locs::AbstractVector, U::SDPermMatrix)
    M = size(U, 1)
    @inbounds pm.perm[locs] = locs[U.perm]
    @inbounds pm.vals[locs] = U.vals
    return pm
end

function _initialize_output(nbit::Int, nctrl::Int, U::SDPermMatrix{T}) where {T}
    N = 1 << nbit
    PermMatrix(collect(1:N), ones(T, N))
end

function cunmat(
    nbit::Int,
    cbits::NTuple{C,Int},
    cvals::NTuple{C,Int},
    U0::SDPermMatrix,
    locs::NTuple{M,Int},
) where {C,M}
    large_mat_check(nbit)
    U, ic, locs_raw = reorder_unitary(nbit, cbits, cvals, U0, locs)
    pm = _initialize_output(nbit, C, U0)
    controldo(ic) do i
        unij!(pm, locs_raw .+ i, U)
    end
    return pm
end

############################ Diagonal ##########################
function _initialize_output(nbit::Int, nctrl::Int, U::SDDiagonal{T}) where {T}
    Diagonal(ones(T, 1 << nbit))
end

function cunmat(
    nbit::Int,
    cbits::NTuple{C,Int},
    cvals::NTuple{C,Int},
    U0::SDDiagonal,
    locs::NTuple{M,Int},
) where {C,M}
    large_mat_check(nbit)
    U, ic, locs_raw = reorder_unitary(nbit, cbits, cvals, U0, locs)
    dg = _initialize_output(nbit, C, U0)
    controldo(ic) do i
        unij!(dg, locs_raw .+ i, U)
    end
    return dg
end

@inline function unij!(dg::SDDiagonal, locs, U)
    @inbounds dg.diag[locs] = U.diag
    return dg
end


text:
#=
NodeWiseTriangleCount:
- Julia version: 1.6.2
- Author: samuel
- Date: 2021-09-12
=#

using SuiteSparseGraphBLAS
using SparseArrays
using LinearAlgebra

include("/Users/samuel/Documents/Jupyter/GraphBlasAlgorithmsFinal/src/NodeWiseTriangleCount.jl")
dim = 10000
matrix = GBMatrix(sparse(rand(0:1, dim, dim)))
@time nodewise_tc(matrix)


text:
function add!(res::AbstractArray{Float64}, x::AbstractArray{Float64}, y::Float64, n::Int64=length(x))
   for i in 1:n
       @inbounds res[i] = x[i] + y
   end
   nothing
end

add!(res::AbstractArray{Float64}, x::Float64, y::AbstractArray{Float64}, n::Int64=length(y)) = add!(res, y, x, n)

add!(x::AbstractArray{Float64}, y::Float64, n::Int64=length(x))=add!(x, x, y, n)
add!(x::Float64, y::AbstractArray{Float64}, n::Int64=length(y))=add!(y, y, x, n)

plusone!(res::AbstractArray{Float64}, x::AbstractArray{Float64}, n::Int64=length(x)) = add!(res, x, 1.0, n)
plusone!(x::AbstractArray{Float64}, n=length(x)) = plusone!(x, x, n)

function divide!(res::AbstractArray{Float64}, x::AbstractArray{Float64}, y::AbstractArray{Float64}, n::Int64=length(x))
    for i in 1:n
        @inbounds res[i] = x[i] / y[i]
    end
    nothing
end

divide!(res::AbstractArray{Float64}, x::AbstractArray{Float64}, y::Float64, n::Int64=length(x)) = multiply!(res, x, 1/y, n)
function divide!(res::AbstractArray{Float64}, x::Float64, y::AbstractArray{Float64}, n::Int64=length(y))
    for i in 1:n
        @inbounds res[i] = x / y[i]
    end
    nothing
end
divide!(x::AbstractArray{Float64}, y::AbstractArray{Float64}, n::Int64=length(x)) = divide!(x, x, y, n)
# divide!(x::Float64, y::AbstractArray{Float64}, n::Int64=length(y)) = divide!(y, x, y, n)
divide!(x::AbstractArray{Float64}, y::Float64, n::Int64=length(x)) = divide!(x, x, y, n)

rcp!(res::AbstractArray{Float64}, x::AbstractArray{Float64}, n::Int64=length(x)) = divide!(res, 1.0, x, n)
rcp!(x::AbstractArray{Float64}, n::Int64=length(x))=rcp!(x, x, n)

function multiply!(res::AbstractArray{Float64}, x::AbstractArray{Float64}, y::Float64, n::Int64=length(x))
    for i in 1:n
        @inbounds res[i] = x[i] * y
    end
    nothing
end

function multiply!(res::AbstractArray{Float64}, x::Float64, y::AbstractArray{Float64}, n::Int64=length(y))
    for i in 1:n
        @inbounds res[i] = x * y[i]
    end
    nothing
end

# multiply!(x::Float64, y::AbstractArray{Float64}, n::Int64=length(y)) = multiply!(y, x, y, n)
multiply!(x::AbstractArray{Float64}, y::Float64, n::Int64=length(x)) = multiply!(x, x, y, n)

function sqr!(res::AbstractArray{Float64}, x::AbstractArray{Float64}, n::Int64=length(x))
   for i in 1:n
       @inbounds res[i] = x[i]*x[i]
   end
   nothing
end
sqr!(x::AbstractArray{Float64}, n::Int64=length(x)) = sqr!(x, x, n)

function H1(y, mu, sigmas)
    (y .- mu)./sigmas./sigmas
end
function H2(y, mu, sigmas)
    z = (y .- mu)./sigmas
    (z.^2 .-1)./sigmas^2./2
end

function H3(y, mu, sigmas)
    z = (y .- mu)./sigmas
    (z.^3 .-3.*z) ./ sigmas^3./6
end

function H4(y, mu, sigmas)
    z = (y .- mu)./sigmas
    (z.^4 .-6.*z.^2 .+ 3) ./ sigmas^4 ./ 24
end


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/09/2014

module DynamicModel

export
    AbstractDynamicModel,
    addObserver,

    initialize,
    update,

    initializeDynamicModel,
    simulateDynamicModel,

    SimpleADM,
    SimpleADMInitialState,
    SimpleADMOutputState,
    SimpleADMCommand,

    LLADM,
    LLADMState,
    LLADMOutputState,
    LLADMCommand

import CommonInterfaces: initialize, update

using AbstractDynamicModelImpl

using SimpleADMImpl

using LLADMImpl

end

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
mutable struct Rule
    left::Any
    right::Any
    expr::Expr # original expression
    mode::Symbol # can be :rewrite or :dynamic
    right_fun::Union{Nothing, Tuple{Vector{Symbol}, Function}}
end

# operator symbols for simple term rewriting
const rewrite_syms = [:(=>)]
# operator symbols for regular pattern matching rules, "dynamic rules"
# that eval the right side at reduction time.
# might be used to implement big step semantics
const dynamic_syms = [:(|>)]

const raw_syms = [:(↦)]

"""
Construct a `Rule` from a quoted expression.
You can also use the [`@rule`] macro to
create a `Rule`.
## Symbolic Rules

Rules defined as `left_hand => right_hand` are
called `symbolic` rules. Application of a `symbolic` Rule
is a replacement of the `left_hand` pattern with
the `right_hand` substitution, with the correct instantiation
of pattern variables. Function call symbols are not treated as pattern
variables, all other identifiers are treated as pattern variables.
Literals such as `5, :e, "hello"` are not treated as pattern
variables.

## Dynamic Rules

Rules defined as `left_hand |> right_hand` are
called `dynamic` rules. Dynamic rules behave like anonymous functions.
Instead of a symbolic substitution, the right hand of
a dynamic `|>` rule is evaluated during rewriting:
matched values are bound to pattern variables as in a
regular function call. This allows for dynamic computation
of

## Type Assertions

Type assertions are supported in the left hand of rules
to match and access literal values both when using classic
rewriting and EGraph based rewriting.
To use a type assertion pattern, add `::T` after
a pattern variable in the `left_hand` of a rule.

---

## Examples

Symbolic rule
```julia
Rule(:(a * b => b * a))
```

```julia
Rule(:(a::Number * b::Number |> a*b))
```
"""
function Rule(e::Expr; mod::Module=@__MODULE__)
    e = rmlines(copy(e))
    mode = :undef
    mode = getfunsym(e)
    l, r = e.args[iscall(e) ? (2:3) : (1:2)]

    right_fun = nothing

    if mode ∈ raw_syms
        mode = :raw
        l = mod.eval(l, mod)
    elseif mode ∈ dynamic_syms # right hand execution, dynamic rules in egg
        mode = :dynamic
        l = interpolate_dollar(l, mod)
        l = df_walk(x -> eval_types_in_assertions(x, mod), l; skip_call=true)
        right_fun = genrhsfun(l, r, mod)
    elseif mode ∈ rewrite_syms # right side is quoted, symbolic replacement
        mode = :rewrite
        l = interpolate_dollar(l, mod)
        l = df_walk(x -> eval_types_in_assertions(x, mod), l; skip_call=true)
    else
        error(`rule "$e" is not in valid form.\n`)
    end

    e.args[iscall(e) ? 2 : 1] = l
    return Rule(l, r, e, mode, right_fun)
end

macro rule(e)
    Rule(e; mod=__module__)
end

# string representation of the rule
function Base.show(io::IO, x::Rule)
    println(io, "Rule(:(", x.expr, "))")
end

"""
Generates a tuple containing the list of formal parameters (`Symbol`s)
and the [`RuntimeGeneratedFunction`](@ref) corresponding to the right hand
side of a `:dynamic` [`Rule`](@ref).
"""
function genrhsfun(left, right, mod::Module)
    # remove type assertions in left hand
    lhs = df_walk( x -> (isexpr(x, :(::)) ? x.args[1] : x), left; skip_call=true )

    # collect variable symbols in left hand
    lhs_vars = Set{Symbol}()
    df_walk( x -> (if x isa Symbol; push!(lhs_vars, x); end; x), left; skip_call=true )
    params = Expr(:tuple, :_egraph, lhs_vars...)

    ex = :($params -> $right)
    (collect(lhs_vars), closure_generator(mod, ex))
end

# TODO is there anything better than eval to use here?
"""
When creating a theory, type assertions in the left hand contain symbols.
We want to replace the type symbols with the real type values, to fully support
the subtyping mechanism during pattern matching.
"""
function eval_types_in_assertions(x, mod::Module)
    if isexpr(x, :(::))
        !(x.args[1] isa Symbol) && error("Type assertion is not on metavariable")
        Expr(:(::), x.args[1], mod.eval(x.args[2]))
    else x
    end
end


text:
module TypesTests

using Gridap.TensorValues
using Test
using Gridap.Arrays: get_array
using StaticArrays

# Constructors (TensorValue)

a = SMatrix{2,2}(1,2,3,4)
t = TensorValue(a)
@test isa(t,TensorValue{2,2,Int})
@test convert(SMatrix{2,2,Int},t) == [1 3;2 4]

a = MMatrix{2,2}(1,2,3,4)
t = TensorValue(a)
@test isa(t,TensorValue{2,2,Int})
@test convert(SMatrix{2,2,Int},t) == [1 3;2 4]

t = TensorValue{2}((1,2,3,4))
@test isa(t,TensorValue{2,2,Int})
@test convert(SMatrix{2,2,Int},t) == [1 3;2 4]

t = TensorValue{2}(1,2,3,4)
@test isa(t,TensorValue{2,2,Int})
@test convert(SMatrix{2,2,Int},t) == [1 3;2 4]

t = TensorValue(1,2,3,4)
@test isa(t,TensorValue{2,2,Int})
@test convert(SMatrix{2,2,Int},t) == [1 3;2 4]

t = TensorValue((1,2,3,4))
@test isa(t,TensorValue{2,2,Int})
@test convert(SMatrix{2,2,Int},t) == [1 3;2 4]

t = TensorValue{1}(10)
@test isa(t,TensorValue{1,1,Int})
@test convert(SMatrix{1,1,Int},t) == 10*ones(1,1)

t = TensorValue{1}((10,))
@test isa(t,TensorValue{1,1,Int})
@test convert(SMatrix{1,1,Int},t) == 10*ones(1,1)

t = TensorValue(1,2.0,3,4)
@test isa(t,TensorValue{2,2,Float64})
@test convert(SMatrix{2,2,Float64},t) == [1 3;2 4]

t = TensorValue{2}(1,2.0,3,4)
@test isa(t,TensorValue{2,2,Float64})
@test convert(SMatrix{2,2,Float64},t) == [1 3;2 4]

t = TensorValue{2,2}(1,2.0,3,4)
@test isa(t,TensorValue{2,2,Float64})
@test convert(SMatrix{2,2,Float64},t) == [1 3;2 4]

t = TensorValue{2,2,Int}(1,2.0,3,4)
@test isa(t,TensorValue{2,2,Int})
@test convert(SMatrix{2,2,Int},t) == [1 3;2 4]

# Constructors (SymTensorValue)

s = SymTensorValue( (11,21,22) )
@test isa(s,SymTensorValue{2,Int})
@test convert(SMatrix{2,2,Int},s) == [11 21;21 22]

s = SymTensorValue(11,21,22)
@test isa(s,SymTensorValue{2,Int})
@test convert(SMatrix{2,2,Float64},s) == [11.0 21.0;21.0 22.0]

s = SymTensorValue{2}( (11,21,22) )
@test isa(s,SymTensorValue{2,Int})
@test convert(SMatrix{2,2,Int},s) == [11 21;21 22]

s = SymTensorValue{2}(11,21,22)
@test isa(s,SymTensorValue{2,Int})
@test convert(SMatrix{2,2,Float64},s) == [11.0 21.0;21.0 22.0]

s = SymTensorValue{2,Int}( (11,21,22) )
@test isa(s,SymTensorValue{2,Int})
@test convert(SMatrix{2,2,Int},s) == [11 21;21 22]

s = SymTensorValue{2,Float64}(11,21,22)
@test isa(s,SymTensorValue{2,Float64})
@test convert(SMatrix{2,2,Float64},s) == [11.0 21.0;21.0 22.0]

s = SymTensorValue{0,Int}( () )
@test isa(s,SymTensorValue{0,Int})
@test convert(SMatrix{0,0,Int},s) == Array{Any,2}(undef,0,0)

s = SymTensorValue{0,Int}()
@test isa(s,SymTensorValue{0,Int})
@test convert(SMatrix{0,0,Int},s) == Array{Any,2}(undef,0,0)

s = SymTensorValue(11,21.0,22)
@test isa(s,SymTensorValue{2,Float64})
@test convert(SMatrix{2,2,Float64},s) == [11.0 21.0;21.0 22.0]

s = SymTensorValue{2}(11,21.0,22)
@test isa(s,SymTensorValue{2,Float64})
@test convert(SMatrix{2,2,Float64},s) == [11.0 21.0;21.0 22.0]

s = SymTensorValue{2,Int}(11,21.0,22)
@test isa(s,SymTensorValue{2,Int})
@test convert(SMatrix{2,2,Int},s) == [11.0 21.0;21.0 22.0]

# Constructors (SymFourthOrderTensorValue)

s = SymFourthOrderTensorValue( (1111,1121,1122, 2111,2121,2122, 2211,2221,2222) )
@test isa(s,SymFourthOrderTensorValue{2,Int})
@test Tuple(s) == (1111,1121,1122, 2111,2121,2122, 2211,2221,2222)

s = SymFourthOrderTensorValue(1111,2111,2211, 1121,2121,2221, 1122,2122,2222)
@test isa(s,SymFourthOrderTensorValue{2,Int})
@test Tuple(s) == (1111,2111,2211, 1121,2121,2221, 1122,2122,2222  )

s = SymFourthOrderTensorValue{2}( (1111,2111,2211, 1121,2121,2221, 1122,2122,2222) )
@test isa(s,SymFourthOrderTensorValue{2,Int})
@test Tuple(s) == (1111,2111,2211, 1121,2121,2221, 1122,2122,2222  )

s = SymFourthOrderTensorValue{2}(1111,2111,2211, 1121,2121,2221, 1122,2122,2222)
@test isa(s,SymFourthOrderTensorValue{2,Int})
@test Tuple(s) == (1111,2111,2211, 1121,2121,2221, 1122,2122,2222  )

s = SymFourthOrderTensorValue{2,Int}( (1111,2111,2211, 1121,2121,2221, 1122,2122,2222) )
@test isa(s,SymFourthOrderTensorValue{2,Int})
@test Tuple(s) == (1111,2111,2211, 1121,2121,2221, 1122,2122,2222  )

s = SymFourthOrderTensorValue{2,Float64}(1111,2111,2211, 1121,2121,2221, 1122,2122,2222)
@test isa(s,SymFourthOrderTensorValue{2,Float64})
@test Tuple(s) == (1111.0,2111.0,2211.0, 1121.0,2121.0,2221.0, 1122.0,2122.0,2222.0)

s = SymFourthOrderTensorValue{0,Int}( () )
@test isa(s,SymFourthOrderTensorValue{0,Int})
@test Tuple(s) == ()

s = SymFourthOrderTensorValue{0,Int}()
@test isa(s,SymFourthOrderTensorValue{0,Int})
@test Tuple(s) == ()

s = SymFourthOrderTensorValue(1111,2111,2211.0, 1121,2121.0,2221, 1122,2122,2222)
@test isa(s,SymFourthOrderTensorValue{2,Float64})
@test Tuple(s) == (1111,2111,2211, 1121,2121,2221, 1122,2122,2222  )

s = SymFourthOrderTensorValue{2}(1111,2111,2211.0, 1121,2121.0,2221, 1122,2122,2222)
@test isa(s,SymFourthOrderTensorValue{2,Float64})
@test Tuple(s) == (1111,2111,2211, 1121,2121,2221, 1122,2122,2222  )

s = SymFourthOrderTensorValue{2,Int}(1111,2111,2211.0, 1121,2121.0,2221, 1122,2122,2222)
@test isa(s,SymFourthOrderTensorValue{2,Int})
@test Tuple(s) == (1111,2111,2211, 1121,2121,2221, 1122,2122,2222  )

# Constructors (VectorValue)

a = SVector(1)
g = VectorValue(a)
@test isa(g,VectorValue{1,Int})
@test convert(SVector{1,Int},g) == [1,]

a = SVector(1,2,3,4)
g = VectorValue(a)
@test isa(g,VectorValue{4,Int})
@test convert(SVector{4,Int},g) == [1,2,3,4]

a = MVector(1,2,3,4)
g = VectorValue(a)
@test isa(g,VectorValue{4,Int})
@test convert(MVector{4,Int},g) == [1,2,3,4]

g = VectorValue{4}((1,2,3,4))
@test isa(g,VectorValue{4,Int})
@test convert(SVector{4,Int},g) == [1,2,3,4]

g = VectorValue{1}((1,))
@test isa(g,VectorValue{1,Int})
@test convert(SVector{1,Int},g) == [1,]

g = VectorValue{0,Int}(())
@test isa(g,VectorValue{0,Int})
@test convert(SVector{0,Int},g) == []

g = VectorValue{4}(1,2,3,4)
@test isa(g,VectorValue{4,Int})
@test convert(SVector{4,Int},g) == [1,2,3,4]

g = VectorValue{1}(1)
@test isa(g,VectorValue{1,Int})
@test convert(SVector{1,Int},g) == [1,]

g = VectorValue{1,Float64}(1)
@test isa(g,VectorValue{1,Float64})
@test convert(SVector{1,Float64},g) == [1,]

g = VectorValue{0,Int}()
@test isa(g,VectorValue{0,Int})
@test convert(SVector{0,Int},g) == []

g = VectorValue{4,Float64}((1,2,3,4))
@test isa(g,VectorValue{4,Float64})
@test convert(SVector{4,Float64},g) == [1,2,3,4]

g = VectorValue{4,Float64}(1,2,3,4)
@test isa(g,VectorValue{4,Float64})
@test convert(SVector{4,Float64},g) == [1,2,3,4]

g = VectorValue(1,2,3,4)
@test isa(g,VectorValue{4,Int})
@test convert(SVector{4,Int},g) == [1,2,3,4]

g = VectorValue((1,2,3,4))
@test isa(g,VectorValue{4,Int})
@test convert(SVector{4,Int},g) == [1,2,3,4]

g = VectorValue(1)
@test isa(g,VectorValue{1,Int})
@test convert(SVector{1,Int},g) == [1,]

g = VectorValue(1.0,2,3.0,4)
@test isa(g,VectorValue{4,Float64})
@test convert(SVector{4,Float64},g) == [1,2,3,4]

g = VectorValue{4}(1.0,2,3.0,4)
@test isa(g,VectorValue{4,Float64})
@test convert(SVector{4,Float64},g) == [1,2,3,4]

g = VectorValue{4,Int}(1.0,2,3.0,4)
@test isa(g,VectorValue{4,Int})
@test convert(SVector{4,Int},g) == [1,2,3,4]

g = VectorValue((1.0,2,3.0,4))
@test isa(g,VectorValue{4,Float64})
@test convert(SVector{4,Float64},g) == [1,2,3,4]

# Initializers

z = zero(TensorValue{3,3,Int,9})
@test isa(z,TensorValue{3,3,Int,9})
@test convert(SMatrix{3,3,Int},z) == zeros(Int,(3,3))

z = zero(SymTensorValue{3,Int})
@test isa(z,SymTensorValue{3,Int,6})
@test convert(SMatrix{3,3,Int},z) == zeros(Int,(3,3))

z = zero(ThirdOrderTensorValue{3,3,3,Int,27})
@test isa(z,ThirdOrderTensorValue{3,3,3,Int,27})
@test Tuple(z) == Tuple(zeros(Int,(27)))

z = zero(SymFourthOrderTensorValue{2,Int})
@test isa(z,SymFourthOrderTensorValue{2,Int,9})
@test Tuple(z) == Tuple(zeros(Int,(9)))

z = zero(VectorValue{3,Int})
@test isa(z,VectorValue{3,Int})
@test convert(SVector{3,Int},z) == zeros(Int,3)

z = one(TensorValue{3,3,Int,9})
@test isa(z,TensorValue{3,3,Int,9})
@test convert(SMatrix{3,3,Int},z) == [1 0 0; 0 1 0; 0 0 1]
s = one(z)
@test convert(SMatrix{3,3,Int},s) == [1 0 0; 0 1 0; 0 0 1]

z = one(SymFourthOrderTensorValue{2,Int})
@test isa(z,SymFourthOrderTensorValue{2})
@test Tuple(z) == (1.,0.,0., 0.,0.5,0., 0.,0.,1.)

r = rand(VectorValue{3,Int})
@test isa(r,VectorValue{3,Int})
@test r ≠ rand(typeof(r))

r = rand(TensorValue{3,3,Int})
@test isa(r,TensorValue{3,3,Int,9})
@test r ≠ rand(typeof(r))

r = rand(ThirdOrderTensorValue{3,3,3,Int})
@test isa(r,ThirdOrderTensorValue{3,3,3,Int,27})
@test r ≠ rand(typeof(r))

r = rand(SymTensorValue{3,Int})
@test isa(r,SymTensorValue{3,Int,6})
@test r ≠ rand(typeof(r))

r = rand(SymFourthOrderTensorValue{3,Int})
@test isa(r,SymFourthOrderTensorValue{3,Int,36})
@test r ≠ rand(typeof(r))

# Conversions

a = @SVector ones(Int,3)
b = convert(VectorValue{3,Int},a)
@test isa(b,VectorValue{3,Int})

a = ones(Int,3)
b = convert(VectorValue{3,Int},a)
@test isa(b,VectorValue{3,Int})

a = ones(Int,1)
b = convert(VectorValue{1,Int},a)
@test isa(b,VectorValue{1,Int})

a = (1,2,2,1,3,2)
V = TensorValue{3,2,Int,6}
b = convert(V,a)
@test isa(b,V)
b = V[a,a,a,]
@test isa(b,Vector{V})

a = (11,21,22)
V = SymTensorValue{2,Int,3}
b = convert(V,a)
@test isa(b,V)
b = V[a,a,a,]
@test isa(b,Vector{V})

a = (1111,1121,1122, 2111,2121,2122, 2211,2221,2222)
V = SymFourthOrderTensorValue{2,Int,9}
b = convert(V,a)
@test isa(b,V)
@test b[2,2,2,1] == 2221
@test b[2,2,1,2] == 2221
@test b[2,1,2,1] == 2121
@test b[1,2,2,1] == 2121
b = V[a,a,a,]
@test isa(b,Vector{V})

# Misc operations on the type itself

V = VectorValue{3,Int}
@test length(V) == 3
@test size(V) == (3,)

V = VectorValue{3}
@test length(V) == 3
@test size(V) == (3,)

# Custom type printing

v = TensorValue{3,2,Float64}(1,2,3,4,5,6)
s = "(1.0, 2.0, 3.0, 4.0, 5.0, 6.0)"
@test string(v) == s

v = SymTensorValue{3,Int}(1, 0, 0, 1, 0, 1)
s = "(1, 0, 0, 1, 0, 1)"
@test string(v) == s

v = SymFourthOrderTensorValue{2,Int}(1111,1121,1122, 2111,2121,2122, 2211,2221,2222)
s = "(1111, 1121, 1122, 2111, 2121, 2122, 2211, 2221, 2222)"
@test string(v) == s

# Third order tensors

a = SArray{Tuple{2,2,2}}(1,2,3,4,5,6,7,8)
t = ThirdOrderTensorValue(a)
@test isa(t,ThirdOrderTensorValue{2,2,2,Int})

t = ThirdOrderTensorValue(1,2,3,4,5,6,7,8)
@test isa(t,ThirdOrderTensorValue{2,2,2,Int})

t = ThirdOrderTensorValue{2}(1,2,3,4,5,6,7,8)
@test isa(t,ThirdOrderTensorValue{2,2,2,Int})

t = ThirdOrderTensorValue{2,2,2}(1,2,3,4,5,6,7,8)
@test isa(t,ThirdOrderTensorValue{2,2,2,Int})

t = ThirdOrderTensorValue{2,2,2,Float64}(1,2,3,4,5,6,7,8)
@test isa(t,ThirdOrderTensorValue{2,2,2,Float64})

t = ThirdOrderTensorValue(1,2.0,3,4,5,6,7,8)
@test isa(t,ThirdOrderTensorValue{2,2,2,Float64})

t = ThirdOrderTensorValue{2}(1,2.0,3,4,5,6,7,8)
@test isa(t,ThirdOrderTensorValue{2,2,2,Float64})

t = ThirdOrderTensorValue{2,2,2}(1,2.0,3,4,5,6,7,8)
@test isa(t,ThirdOrderTensorValue{2,2,2,Float64})

t = ThirdOrderTensorValue{2,2,2,Int}(1,2.0,3,4,5,6,7,8)
@test isa(t,ThirdOrderTensorValue{2,2,2,Int})

# Misc

v = VectorValue(3,2,1)
m = mutable(v)
@test m == get_array(v)
@test isa(m,MVector)

v = TensorValue{2,3}(1,2,3,4,5,6)
m = mutable(v)
@test m == get_array(v)
@test isa(m,MMatrix)

M = Mutable(VectorValue{3,Int})
@test M == MVector{3,Int}
m = zero(M)
v = VectorValue(m)
@test isa(v,VectorValue{3,Int})

@test num_components(Int) == 1
@test num_components(Float64) == 1
@test num_components(1.0) == 1
@test num_components(1) == 1
@test num_components(VectorValue{3,Float64}) == 3
@test num_components(VectorValue(1,2,3)) == 3
@test num_components(TensorValue(1,2,3,4)) == 4
@test num_components(SymTensorValue(1,2,3)) == 4
@test num_components(SymFourthOrderTensorValue(1111,1121,1122, 2111,2121,2122, 2211,2221,2222)) == 16

a = VectorValue(1,2,3,4)
@test change_eltype(a,Float64) == VectorValue{4,Float64}

a = TensorValue(1,2,3,4)
@test change_eltype(a,Float64) == TensorValue{2,2,Float64,4}

@test change_eltype(1,Float64) == Float64

a = TensorValue(1,2,3,4)
@test isa(Tuple(a),Tuple)
@test Tuple(a) == a.data

p = VectorValue(2,3)
t = diagonal_tensor(p)
@test t == TensorValue(2,0,0,3)

p = VectorValue(1,2,3)
t = diagonal_tensor(p)
@test t == TensorValue(1,0,0,0,2,0,0,0,3)

a = SymTensorValue(11,21,22)
@test change_eltype(a,Float64) == SymTensorValue{2,Float64,3}
@test isa(Tuple(a),Tuple)
@test Tuple(a) == a.data

a = SymFourthOrderTensorValue(1111,1121,1122, 2111,2121,2122, 2211,2221,2222)
@test change_eltype(a,Float64) == SymFourthOrderTensorValue{2,Float64,9}
@test isa(Tuple(a),Tuple)
@test Tuple(a) == a.data

end # module TypesTests


text:
# Miscellaneous and alternative algorithms to compare with the primary
# implementations. Not exported.

struct Poehn1985 <: NeXLAlgorithm end

"""
    jumpratio(ashell::AtomicSubShell, ::Type{Poehn1985})

An implement of jump ratios attributed to
  * Poehn, Wernisch, Hanke (1985) X-ray Spectrom 14(3):120, 1985
"""
function jumpratio(ashell::AtomicSubShell, ::Type{Poehn1985})
    zz = ashell.z
    if (zz > 11) && (zz < 50) && (ashell.subshell.index == 1)
        return 17.54 - 0.6608 * zz + 0.01427 * zz^2 - 1.1e-4 * zz^3
    elseif (zz > 39) && (zz < 83)
        if ashell.subshell.index == 4
            return 20.03 - 0.7732zz + 0.01159 * zz^2 - 5.835e-5 * zz^3
        elseif ashell.subshell.index == 3
            return 1.41
        elseif ashell.shel.index == 2
            return 1.16
        end
    end
    error("Unsupported edge: $(ashell)")
end

"""
    klinewidths(elm::Element)

Linewidth of the K shell according to Bambynek'1974 errata to Bambynek 1972.
Shown to be reliable for Z>36 or so.
"""
klinewidths(elm::Element) = 1.73e-6 * z(elm)^3.93

struct Burhop1965 <: NeXLAlgorithm end

"""
    fluorescenceyield(z::Int, ::Type{Burhop1965})

An approximate expression for the K-shell fluorescence yield due to
E.H.S Burhop, J. Phys. Radium, 16, 625 (1965)
"""
function fluorescenceyield(z::Int, ::Type{Burhop1965})
    d = -0.044 + 0.0346z - 1.35e-6z^3
    return d^4 / (1.0 + d^4)
end

fluorescenceyield(z::Int) = fluorescenceyield(z, Burhop1965)


text:
function fib_seeder{T<:Integer}(n::T)
    a = zeros(BigInt, n)
    a[1] = one(BigInt)
    return a
end

function fib{T<:Integer}(n::T, k::T)
    NFib(n, k, fib_seeder)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using OptimizationLH
using Test

@testset "Model" begin
    mdl = Model();

    # Helpers for computing stats
    choiceCl_kV, Nc = OptimizationLH.choice_classes(mdl);
    @test all(choiceCl_kV .>= 1)
    @test all(choiceCl_kV .<= Nc)
    @test isa(choiceCl_kV, Vector{Int})

    mp = OptimizationLH.make_random_parameters(mdl, 215);
    @test length(mp.alphaV) == mdl.K

    # Handling guesses
    guessV = OptimizationLH.mp_to_guess(mdl, mp);
    @test length(guessV) == OptimizationLH.n_params(mdl)
    @test isa(guessV, Vector{Float64})
    mp2 = OptimizationLH.guess_to_mp(mdl, guessV);
    @test mp2.alphaV == mp.alphaV
    guess2V = OptimizationLH.mp_to_guess(mdl, mp2);
    @test guessV == guess2V

    # Simulation
    prob_mjM = rand(mdl.M, mdl.N);
    choice_sjM = OptimizationLH.simulate_choices(mdl, prob_mjM ./ sum(prob_mjM, dims = 1), 242);
    @test isa(choice_sjM, Matrix{Int})
    @test size(choice_sjM) == (mdl.nSim, mdl.N)
    @test all(choice_sjM .>= 1)
    @test all(choice_sjM .<= mdl.M)

    # Solving model
    mStats = OptimizationLH.solve(mdl, mp);

    # Check that perturbing each guess changes objective
    @test OptimizationLH.check_changing_guesses(mdl);

    @test OptimizationLH.calc_many_deviations();

    # Calibration
    calibrate(mdl, :LN_SBPLX)
end

# --------------


text:
export Trajectory

"""
    Trajectory{names,types,Tbs}(trajectories::Tbs)

A container of different `trajectories`.
Usually you won't use it directly.
"""
struct Trajectory{names,types,Tbs} <: AbstractTrajectory{names,types}
    trajectories::Tbs
end

"A helper function to access inner fields"
Base.getindex(t::Trajectory, s::Symbol) = getproperty(t.trajectories, s)

get_trace(t::Trajectory, s::Symbol) = t[s]

function Base.push!(t::Trajectory, kv::Pair{Symbol})
    k, v = kv
    push!(t[k], v)
end

Base.pop!(t::Trajectory, s::Symbol) = pop!(t[s])


text:
module StructuralEstimation

using Reexport
@reexport using StructuralBase
@reexport using MDPTools
@reexport using NestedFixedPoint
@reexport using NestedPseudoLikelihood

end # module


text:
using Documenter
using Fluxion

makedocs(
    modules = [Fluxion],
    format = :html,
    sitename = "Fluxion.jl",
    authors = "Bogdan Opanchuk",
    pages = [
        "Home" => "index.md",
    ],
)

deploydocs(
    repo = "github.com/fluxion/Fluxion.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
const DEFAULT_ATTRIBS = Dict(
    "time" => Dict(
        "units" => "days since 0000-00-01 00:00:00",
        "standard_name" => "time"
    ),
    "lon" => Dict(
        "units" => "degrees_east",
        "standard_name" => "longitude",
        "valid_range" => Float32[-180.0, 360.0]
    ),
    "lat" => Dict(
        "units" => "degrees_north",
        "standard_name" => "latitude",
        "valid_range" => Float32[-90.0, 90.0]
    ),
    "level" => Dict(
        "units" => "millibars",
        "long_name" => "pressure_level",
    ),
)

"""
    ncwrite(file::String, Xs; globalattr = Dict())
Write the given `ClimArray` instances (any iterable of `ClimArray`s or a single `ClimArray`)
to a `.nc` file following CF standard conventions using NCDatasets.jl.
Optionally specify global attributes for the `.nc` file.

The metadata of the arrays in `Xs`, as well as their dimensions, are properly written
in the `.nc` file and any necessary type convertions happen automatically.

**WARNING**: We assume that any dimensions shared between the `Xs` are identical.

See also [`ncread`](@ref).
"""
function ncwrite(file::String, X::ClimArray; globalattr = Dict())
    ncwrite(file, (X,); globalattr)
end
function ncwrite(file::String, Xs; globalattr = Dict())
    NCDataset(file, "c"; attrib = globalattr) do ds
        for (i, X) in enumerate(Xs)
            n = string(X.name)
            if n == ""
                n = "x$i"
                @warn "$i-th ClimArray has no name, naming it $(n) instead."
            end
            println("processing variable $(n)...")
            add_dims_to_ncfile!(ds, dims(X))
            attrib = X.attrib
            if (isnothing(attrib) || attrib == DimensionalData.NoMetadata())
                attrib = Dict()
            end
            dnames = dim_to_commonname.(dims(X))
            data = Array(X)
            NCDatasets.defVar(ds, n, data, (dnames...,); attrib)
        end
    end
end

function add_dims_to_ncfile!(ds::NCDatasets.AbstractDataset, dimensions::Tuple)
    dnames = dim_to_commonname.(dimensions)
    dims_in_ds = [x[1] for x in ds.dim]
    for (d, dname) ∈ zip(dimensions, dnames)
        dname ∈ dims_in_ds && continue
        println("writing dimension $dname...")
        v = gnv(d); l = length(v)
        NCDatasets.defDim(ds, dname, l) # add dimension entry
        if d isa Coord
            # Define clon/clat variables with this dimension
            lons = getindex.(v, 1); lats = getindex.(v, 2)
            NCDatasets.defVar(ds, "clon", lons, (dname, ); attrib = DEFAULT_ATTRIBS["lon"])
            NCDatasets.defVar(ds, "clat", lats, (dname, ); attrib = DEFAULT_ATTRIBS["lat"])
        else
            # this conversion to DateTime is necessary because CFTime.jl doesn't support Date
            eltype(v) == Date && (v = DateTime.(v))
            attrib = DimensionalData.metadata(d)
            if (isnothing(attrib) || attrib == DimensionalData.NoMetadata()) &&
                    haskey(DEFAULT_ATTRIBS, dname)
                @warn "Dimension $dname has no attributes, adding default attributes."
                attrib = DEFAULT_ATTRIBS[dname]
            end
            # write dimension values as a variable as well (mandatory)
            NCDatasets.defVar(ds, dname, v, (dname, ); attrib = attrib)
        end
    end
end


text:
module ChipSort
import Base.pop!

export
    chipsort_small!, chipsort_medium!, chipsort_large, chipsort!,
    sort_net,
    transpose_vecs, transpose_blocks!, transpose_vecs_tall, transpose_vecs_wide, transpose!,
    bitonic_merge, merge_vecs, build_multi_merger, bitonic_merge_interleaved, chipsort_merge_medium,
    DataBuffer, MergeNode, pop!,
    combsort!, insertion_sort!,
    sort_blocks!, sort_vecs!

include("utils.jl")
include("sorting-networks.jl")
include("transpose-vecs.jl")
include("bitonic-merge-network.jl")
include("sort-array.jl")
include("comb-sort.jl")
include("k-way-merge.jl")
include("sort-array-medium.jl")
include("data-buffers.jl")

end # module


text:
module Clang

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("Clang was not build properly. Please run Pkg.build(\"Clang\").")
end
include(depsjl_path)
# Module initialization function
function __init__()
    check_deps()
end

include("cindex.jl")
include("wrap_c.jl")
include("wrap_cpp.jl")

export cindex
export wrap_c
export wrap_cpp

end


text:
"""
    clean_name(obj)::Symbol
Normalize the value such that it qualifies as a lowercase valid identifier.
# Example
```jldoctest
julia> EconJobMarket.clean_name(Symbol("3This  is a Cat:egory!4cΩ."))
:_3this_is_a_cat_egory!4cΩ_
```
"""
clean_name(obj)::Symbol =
    join(isascii(elem) ? lowercase(elem) : elem for elem ∈ string(obj)) |>
    normalizename


text:
function findOccurrences(text, first, second)
    str_list = split(text, " ");
    ind = 1
    result = []

    while true
        third = try
            str_list[ind+2]
        catch e
            e isa BoundsError && return result
        end
        _first = str_list[ind]
        _second = str_list[ind+1]

        if _first == first && _second == second
            push!(result, third)
        end
        ind += 1
    end
end

# test below
using Test
@test findOccurrences("we will we will rock you", "we", "will") == ["we", "rock"]
@test findOccurrences("alice is a good girl she is a good student", "a", "good") == ["girl","student"]

max_stars_repo_name: ccqpein/Arithmetic-Exercises

max_issues_repo_name: ccqpein/Arithmetic-Exercises

max_forks_repo_name: ccqpein/Arithmetic-Exercises


text:
# deprecations
@deprecate measure_collapseto!(args...; collapseto = 0, kwargs...) measure!(
    ResetTo(collapseto),
    args...;
    kwargs...,
)
@deprecate measure_resetto!(args...; resetto = 0, kwargs...) measure!(
    ResetTo(resetto),
    args...;
    kwargs...,
)
@deprecate measure_remove!(args...; kwargs...) measure!(
    RemoveMeasured(),
    args...;
    kwargs...,
)


text:
using EnvelopedArrays
using Base.Test

@testset "EnvelopedArrays.jl" begin
   # Write your own tests here.
   a = [1.0 2.0 3.0; 4.0 5.0 6.0]
   x = EnvelopedArray(1.0,a)
   x_vcat = vcat(1.0,a)
   @test all( x .== x_vcat)
   @test x[1] == x_vcat[1]
   @test x[2] == x_vcat[2]
   @test x[end] == x_vcat[end]
   @test x[2,2] == x_vcat[2,2]
   @test all(x[1:3,1] .== x_vcat[1:3,1])

   y = EnvelopedArray(1,a)
   @test x == y

   @test envelope(x) == 1.0
   @test parent(x) == a
end


text:
# test stacks and queues
@testset "Stacks" begin

    # Stack

    s = Stack(Int, 5)
    n = 100

    @test isa(s, Stack{Int})
    @test length(s) == 0
    @test isempty(s)
    @test_throws ArgumentError top(s)
    @test_throws ArgumentError pop!(s)

    for i = 1 : n
        push!(s, i)
        @test top(s) == i
        @test !isempty(s)
        @test length(s) == i
    end

    for i = 1 : n
        x = pop!(s)
        @test x == n - i + 1
        if i < n
            @test top(s) == n - i
        else
            @test_throws ArgumentError top(s)
        end
        @test isempty(s) == (i == n)
        @test length(s) == n - i
    end

    s = Stack(Int,1)
    push!(s, 10)
    @test length(empty!(s)) == 0

    #test that iter returns a LIFO collection

    stk = Stack(Int, 10)
    #an array to check iteration sequence against
    arr = Int64[]

    for i = 1:n
        push!(stk,i)
        push!(arr,i)
    end

    #test iterator

    index = length(arr)
    for i in stk
        @test(arr[index] == i)
        index -= 1
    end

    index = 1
    for i in reverse_iter(stk)
        @test(arr[index] == i)
        index += 1
    end

    @test arr == [i for i in reverse_iter(stk)]
    @test reverse(arr) == [i for i in stk]

end # @testset Stack

# Queue
@testset "Queue" begin

    s = Queue(Int, 5)
    n = 100

    @test isa(s, Queue{Int})
    @test length(s) == 0
    @test isempty(s)
    @test_throws ArgumentError front(s)
    @test_throws ArgumentError back(s)
    @test_throws ArgumentError dequeue!(s)

    for i = 1 : n
        enqueue!(s, i)
        @test front(s) == 1
        @test back(s) == i
        @test !isempty(s)
        @test length(s) == i
    end

    for i = 1 : n
        x = dequeue!(s)
        @test x == i
        if i < n
            @test front(s) == i + 1
            @test back(s) == n
        else
            @test_throws ArgumentError front(s)
            @test_throws ArgumentError back(s)
        end
        @test isempty(s) == (i == n)
        @test length(s) == n - i
    end

    #test that iter returns a FIFO collection

    q = Queue(Int, 10)
    #an array to check iteration sequence against
    arr = Int64[]

    for i = 1:n
        enqueue!(q,i)
        push!(arr,i)
    end

    #test iterator

    index = 1
    for i in q
        @test(arr[index] == i)
        index += 1
    end

    index = length(arr)
    for i in reverse_iter(q)
        @test(arr[index] == i)
        index -= 1
    end

    @test arr == [i for i in q]
    @test reverse(arr) == [i for i in reverse_iter(q)]

end # @testset Queue


text:
struct NodeSet <: FEMContainer
    node_set_id::NodeSetID
    num_nodes::Int64
    nodes::Vector{Int64}
    function NodeSet(exo_id::ExoID, node_set_id::NodeSetID)
        num_nodes, _ = read_node_set_parameters(exo_id, node_set_id)
        node_set_nodes = read_node_set_nodes(exo_id, node_set_id)
        return new(node_set_id, num_nodes, node_set_nodes)
    end
end

Base.show(io::IO, node_set::NodeSet) =
print(io, "NodeSet:\n",
          "\tNode set ID     = ", node_set.node_set_id, "\n",
          "\tNumber of nodes = ", node_set.num_nodes, "\n")

# more verbose types
#
NodeSets = Vector{NodeSet}
NodeSetIDs = Vector{NodeSetID}

function read_node_set_ids(exo_id::ExoID, num_node_sets::Int64)
    node_set_ids = Array{NodeSetID}(undef, num_node_sets)
    error = ccall((:ex_get_ids, libexodus), ExodusError,
                  (ExoID, ExodusConstant, Ref{NodeSetID}),
                  exo_id, EX_NODE_SET, node_set_ids)
    exodus_error_check(error, "read_node_set_ids")
    return node_set_ids
end

function read_node_set_parameters(exo_id::ExoID, node_set_id::NodeSetID)
    num_nodes = Ref{Int64}(0)
    num_df = Ref{Int64}(0)
    error = ccall((:ex_get_set_param, libexodus), ExodusError,
                  (ExoID, ExodusConstant, NodeSetID, Ref{Int64}, Ref{Int64}),
                  exo_id, EX_NODE_SET, node_set_id, num_nodes, num_df)
    exodus_error_check(error, "read_node_set_parameters")
    return num_nodes[], num_df[]
end

function read_node_set_nodes(exo_id::ExoID, node_set_id::NodeSetID)
    num_nodes, _ = read_node_set_parameters(exo_id, node_set_id)
    node_set_nodes = Array{Int32}(undef, num_nodes)
    error = ccall((:ex_get_node_set, libexodus), ExodusError,
                  (ExoID, NodeSetID, Ref{Int32}),
                  exo_id, node_set_id, node_set_nodes)
    exodus_error_check(error, "read_node_set_nodes")
    return node_set_nodes
end

function read_node_sets!(exo_id::ExoID, node_set_ids::NodeSetIDs, node_sets::NodeSets)
    for (n, node_set_id) in enumerate(node_set_ids)
        node_sets[n] = NodeSet(exo_id, node_set_id)
    end
end

function read_node_sets(exo_id::ExoID, node_set_ids::Array{NodeSetID})
    # node_sets = Array{NodeSet}(undef, size(node_set_ids, 1))
    node_sets = NodeSets(undef, size(node_set_ids, 1))
    # for (n, node_set_id) in enumerate(node_set_ids)
    #     node_sets[n] = NodeSet(exo_id, node_set_id)
    # end
    read_node_sets!(exo_id, node_set_ids, node_sets)
    return node_sets
end