text:
__precompile__
include("quadrature.jl")
include("Basis.jl")
include("TimeSolver.jl")

using ProgressMeter
using LinearAlgebra
using Einsum

struct DLRSolver
    # spatial grid of cell interfaces
    x;

    # quadrature
    q::Quadrature;

    # DLRSolver settings
    settings::Settings;

    # spatial basis functions
    basis::Basis;

    # low-rank solution matrices
    X::Array{Float64,2}
    W::Array{Float64,2}
    S::Array{Float64,2}

    # preallocated matrices for Rhs
    A::Array{Float64,3}
    B::Array{Float64,3}
    Y::Array{Float64,3}
    Y1::Array{Float64,2}
    ACons::Array{Float64,3}

    fluxS::Array{Float64,2}
    fluxL::Array{Float64,2}

    yL::Array{Float64,2}
    yS::Array{Float64,2}
    yK::Array{Float64,2}

    # Dirichlet BCs
    uL::Array{Float64,1}
    uR::Array{Float64,1}

    # time solver
    rkUpdate::TimeSolver

    # tridiagonal stencil matrices
    L1I::Tridiagonal{Float64, Vector{Float64}};
    L2::Tridiagonal{Float64, Vector{Float64}};

    # constructor
    function DLRSolver(settings)
        x = settings.x;
        r = settings.r;
        q = Quadrature(settings.Nq,"Gauss");
        basis = Basis(q,settings);

        # note that these are actually the hat variables
        X = zeros(settings.Nx,r)
        S = zeros(r,r)
        W = zeros(settings.N,r)

        A = zeros(r,r,r);
        B = zeros(r,r,settings.Nq^2);
        Y = zeros(r,r,r);
        Y1 = zeros(r,r);

        fluxS = zeros(r,r);
        fluxL = zeros(r,settings.N^2);
        yL = zeros(r,settings.N^2);
        yS = zeros(r,r);
        yK = zeros(settings.Nx,r);

        uL = zeros(settings.Nq^2);
        uR = zeros(settings.Nq^2);

        ACons = zeros(settings.NCons,settings.N^2,settings.N^2);

        rkUpdate = TimeSolver(settings);

        onesBC = ones(settings.Nx);
        onesBC[1] = 0.0; onesBC[end] = 0.0;
        onesBCL = ones(settings.Nx-1);
        onesBCL[end] = 0.0;
        onesBCR = ones(settings.Nx-1);
        onesBCR[1] = 0.0;

        L1I = 0.5*Tridiagonal(onesBCL,-2*onesBC,onesBCR)/settings.dt;
        L2 = 0.25*Tridiagonal(-onesBCL,zeros(settings.Nx),onesBCR)/settings.dx;

        new(x,q,settings,basis,X,S,W,A,B,Y,Y1,ACons,fluxS,fluxL,yL,yS,yK,uL,uR,rkUpdate,L1I,L2);
    end
end

function projector(u::Array{Float64,1},a::Array{Float64,1})
    factor = u'a/(u'u);
    return  factor .* u;
end

function qrGramSchmidt(A::Array{Float64,2})
    XNew,S = qr(A);
    XNew = Matrix(XNew);
    S = Matrix(S);
    r = size(S,2);
    X = XNew[:,1:r]
    S = S[1:r,1:r]
    return X,S;
    N = size(A,1);
    r = size(A,2);

    Q = zeros(N,r);
    R = zeros(r,r);

    for k = 1:r
        Q[:,k] .= A[:,k];
        for j = 1:k-1
            Q[:,k] .-= projector(Q[:,j],A[:,k]);
        end
    end

    for k = 1:r
        Q[:,k] ./= sqrt(Q[:,k]'Q[:,k]); # normalize
    end

    for k = 1:r
        for j = r:-1:k
            R[k,j] = Q[:,k]'A[:,j];
        end
    end

    return Q,R;

end

function RhsK(obj::DLRSolver,K::Array{Float64,2},W::Array{Float64,2})
    Nx = obj.settings.Nx;
    r = obj.settings.r; # DLR rank
    Nq = obj.settings.Nq;
    N = obj.settings.N;
    fXi = 0.25;
    dt = obj.settings.dt;
    flux = zeros(r);
    dx = obj.settings.dx;

    WQuad = EvalAtQuad(obj.basis,W)';

    # Compute A_{i,j,m} = E[W_i W_j W_m]
    WQuad = EvalAtQuad(obj.basis,W)';
    for i = 1:r
        for j = 1:r
            for m = 1:r
                obj.A[i,j,m] = Integral(obj.q,WQuad[i,:].*WQuad[j,:].*WQuad[m,:].*fXi);
            end
        end
    end

    for j = 2:(Nx-1)
        for p = 1:r
            flux[p] = 0.0;
            for l = 1:r
                for m = 1:r
                    flux[p] += 1/4/dx * (K[j+1,l]*K[j+1,m]-K[j-1,l]*K[j-1,m])*obj.A[l,m,p];
                end
            end
        end
        for p = 1:r
            obj.yK[j,p] = 1/2/dt .* (K[j+1,p]-2*K[j,p]+K[j-1,p]) - flux[p];
        end
    end
    return obj.yK;
end

function RhsS(obj::DLRSolver,X::Array{Float64,2},W::Array{Float64,2},L::Array{Float64,2})
    Nx = obj.settings.Nx;
    r = obj.settings.r; # DLR rank
    Nq = obj.settings.Nq^2;
    N = obj.settings.N^2;
    fXi = 0.25;
    dt = obj.settings.dt;
    # Compute A_{i,j,m} = E[L_i L_j phi_m]
    LQuad = EvalAtQuad(obj.basis,L)';
    WQuad = EvalAtQuad(obj.basis,W)';
    for i = 1:r
        for j = 1:r
            for m = 1:r
                obj.A[i,j,m] = Integral(obj.q,LQuad[i,:].*LQuad[j,:].*WQuad[m,:].*fXi);
            end
        end
    end

    for p = 1:r
        for l = 1:r
            obj.Y1[p,l] = 0.0;
            for j = 2:(Nx-1)
                if obj.settings.stabilization != 1
                    if obj.settings.stabilization == 2
                        obj.Y1[p,l] += -1/2/dt * X[j,p]*(X[j+1,l]-2*X[j,l]+X[j-1,l]);# DLR first, discretize second; stable version for projector splitting
                    else
                        obj.Y1[p,l] += 1/2/dt * X[j,p]*(X[j+1,l]-2*X[j,l]+X[j-1,l]); # discretize first, DLR second
                    end
                end
            end
        end
    end

    for p = 1:r
        for m = 1:r
            for l = 1:r
                obj.Y[m,l,p] = 0;
                for j = 2:(Nx-1)
                    obj.Y[m,l,p] += 1/4/obj.settings.dx * X[j,p]*(X[j+1,l]*X[j+1,m]-X[j-1,l]*X[j-1,m]);
                end
            end
        end
    end

    for q = 1:r
        for l = 1:r
            obj.fluxS[q,l] = 0.0;
            for p = 1:r
                for m = 1:r
                    obj.fluxS[q,l] += obj.Y[m,p,q]*obj.A[m,p,l];
                end
            end
        end
        for l = 1:r
            obj.yS[q,l] = 0.0
            for m = 1:r
                for i = 1:N
                    obj.yS[q,l] += obj.Y1[q,m]*L[i,m]*W[i,l];
                end
            end
            obj.yS[q,l] -= obj.fluxS[q,l];
        end
    end
    return obj.yS;

end

function RhsL(obj::DLRSolver,X::Array{Float64,2},L::Array{Float64,2},recompute::Bool=true)
    Nx = obj.settings.Nx;
    r = obj.settings.r; # DLR rank
    Nq = obj.settings.Nq^2;
    N = obj.settings.N^2;
    fXi = 0.25;
    # Compute A_{i,j,m} = E[L_i L_j phi_m]
    LQuad = EvalAtQuad(obj.basis,L)';
    for i = 1:r
        for j = 1:r
            for m = 1:N
                obj.B[i,j,m] = Integral(obj.q,LQuad[i,:].*LQuad[j,:].*obj.basis.PhiQuad[:,m].*fXi);
            end
        end
    end

    if recompute
        for p = 1:r
            for m = 1:r
                for l = 1:r
                    obj.Y[m,l,p] = 0;
                    for j = 2:(Nx-1)
                        obj.Y[m,l,p] += 1/4/obj.settings.dx * X[j,p]*(X[j+1,l]*X[j+1,m]-X[j-1,l]*X[j-1,m]);
                    end
                end
            end
        end

        for p = 1:r
            for l = 1:r
                obj.Y1[p,l] = 0.0;
                for j = 2:(Nx-1)
                    if obj.settings.stabilization != 1
                        obj.Y1[p,l] += 1/2/obj.settings.dt * X[j,p]*(X[j+1,l]-2*X[j,l]+X[j-1,l]);
                    end
                end
            end
        end
    end

    for p = 1:r
        for i = 1:N
            obj.fluxL[p,i] = 0.0;
            for l = 1:r
                for m = 1:r
                    obj.fluxL[p,i] += obj.Y[m,l,p]*obj.B[m,l,i];
                end
            end
        end
    end

    for p = 1:r
        for i = 1:N
            obj.yL[p,i] = 0.0;
            for l = 1:r
                obj.yL[p,i] += obj.Y1[p,l]*L[i,l];
            end
            obj.yL[p,i] -= obj.fluxL[p,i];
        end
    end

    return obj.yL;
end

function SetupIC(obj::DLRSolver)
    u = zeros(obj.settings.N*obj.settings.N,obj.settings.Nx);
    uCons = zeros(obj.settings.Nx,obj.settings.NCons);
    uVals = zeros(obj.settings.Nq^2)
    for j = 1:obj.settings.Nx
        for q = 1:obj.settings.Nq
            for k = 1:obj.settings.Nq
                uVals[(k-1)*obj.settings.Nq+q] = obj.settings.IC(obj.x[j],obj.q.xi[k],obj.q.xi[q])[1];
            end
        end
        u[:,j] = ComputeMomentsDLR(obj.basis,uVals*0.25);
        uCons[j,:] = ComputeMomentsCons(obj.basis,uVals*0.25);
    end
    return u,uCons;
end

function FactorIC(obj::DLRSolver)
    r = obj.settings.r;
    # Set up initial condition
    u, uCons = SetupIC(obj);
    u = Array(u');

    # Low-rank approx of init data:
    X,S,W = svd(u);

    # rank-r truncation:
    X = X[:,1:r];
    W = W[:,1:r];
    S = Array(Diagonal(S));
    S = S[1:r, 1:r];
    return X,S,W,uCons;
end

# matrix projector-splitting integrator
function SolveBackward(obj::DLRSolver)
    t = 0.0;
    dt = obj.settings.dt;
    dx = obj.x[2]-obj.x[1];
    tEnd = obj.settings.tEnd;
    Nx = obj.settings.Nx;
    r = obj.settings.r; # DLR rank
    Nq = obj.settings.Nq^2;
    N = obj.settings.N^2; # here, N is the number of quadrature points

    fXi = 0.25;

    # Set up initial condition
    u,uCons = SetupIC(obj);

    # Low-rank approx of init data:
    X,S,W = svd(u');

    # rank-r truncation:
    X = X[:,1:r];
    W = W[:,1:r];
    S = Diagonal(S);
    S = S[1:r, 1:r];

    K = zeros(Nx,r);
    K1 = zeros(Nx,r);
    KNew = X*S;
    L = zeros(N,r);
    L1 = zeros(N,r);
    S1 = zeros(r,r);

    Nt = Integer(round(tEnd/dt));

    # compute Dirichlet values if they are independent of time
    for k = 1:obj.settings.Nq
        for q = 1:obj.settings.Nq
            obj.uL[(k-1)*obj.settings.Nq+q] = obj.settings.IC(obj.x[1],obj.q.xi[k],obj.q.xi[q])[1];
            obj.uR[(k-1)*obj.settings.Nq+q] = obj.settings.IC(obj.x[end],obj.q.xi[k],obj.q.xi[q])[1];
        end
    end

    prog = Progress(Nt,1)

    # time loop
    #@gif
    for n = 1:Nt

        ###### K-step ######
        K .= X*S;

        WQuad = EvalAtQuad(obj.basis,W)';

        # impose BCs
        for i = 1:r
            K[1,i] = Integral(obj.q,obj.uL.*WQuad[i,:].*fXi);
            K[end,i] = Integral(obj.q,obj.uR.*WQuad[i,:].*fXi);
        end

        if obj.settings.rkType == "Heun"
            K1 .= K .+ dt*RhsK(obj,K,W);
            K1 .= K1 .+ dt*RhsK(obj,K1,W);
            K .= 0.5.*(K.+K1);
        elseif obj.settings.rkType == "Euler"
            K .= K .+ dt*RhsK(obj,K,W);
        elseif obj.settings.rkType == "SSP"
            K .= UpdateK(obj.rkUpdate,obj,K,W);
        end

        X,S = qr(K); # optimize by choosing XFull, SFull
        X = X[:, 1:obj.settings.r];
        S = S[1:obj.settings.r, 1:obj.settings.r];

        ###### S-step ######

        L .= W*S';
        if obj.settings.rkType == "Heun"
            S1 .= S .- dt.*RhsS(obj,X,W,L);
            L .= W*S';
            S1 .= S1 .- dt.*RhsS(obj,X,W,L);
            S .= 0.5.*(S.+S1);
        elseif obj.settings.rkType == "Euler"
            S .= S .- dt.*RhsS(obj,X,W,L);
        elseif obj.settings.rkType == "SSP"
            S .= UpdateS(obj.rkUpdate,obj,X,S,W);
        end

        ###### L-step ######
        L .= W*S';

        if obj.settings.rkType == "Heun"
            L1 .= L .+ dt*RhsL(obj,X,L,false)';
            L1 .= L1 .+ dt*RhsL(obj,X,L1,false)';
            L .= 0.5.*(L.+L1);
        elseif obj.settings.rkType == "Euler"
            L .= L .+ dt*RhsL(obj,X,L,false)';
        elseif obj.settings.rkType == "SSP"
            L .= UpdateL(obj.rkUpdate,obj,X,L,false);
        end

        W,S = qr(L);
        #W,S = qr(L);
        W = W[:, 1:obj.settings.r];
        S = S[1:obj.settings.r, 1:obj.settings.r];

        S .= S';

        # apply filter step
        Filter(obj,W);

        next!(prog) # update progress bar

        t = t+dt;
    end

    # return end time and solution
    return t, X,S,W;
end

# unconventional integrator
function SolveForward(obj::DLRSolver)
    t = 0.0;
    dt = obj.settings.dt;
    dx = obj.x[2]-obj.x[1];
    tEnd = obj.settings.tEnd;
    Nx = obj.settings.Nx;
    r = obj.settings.r; # DLR rank
    Nq = obj.settings.Nq^2;
    N = obj.settings.N^2; # here, N is the number of quadrature points

    fXi = 0.25;

    # Set up initial condition
    u,uCons = SetupIC(obj);

    # Low-rank approx of init data:
    X,S,W = svd(u');

    # rank-r truncation:
    X = X[:,1:r];
    W = W[:,1:r];
    S = Array(Diagonal(S));
    S = S[1:r, 1:r];

    K = zeros(Nx,r);
    K1 = zeros(Nx,r);
    KNew = X*S;
    L = zeros(N,r);
    L1 = zeros(N,r);
    S1 = zeros(r,r);

    PhiQuad = obj.basis.PhiQuad;

    Nt = Integer(round(tEnd/dt));

    # compute Dirichlet values if they are independent of time
    for k = 1:obj.settings.Nq
        for q = 1:obj.settings.Nq
            obj.uL[(k-1)*obj.settings.Nq+q] = obj.settings.IC(obj.x[1],obj.q.xi[k],obj.q.xi[q])[1];
            obj.uR[(k-1)*obj.settings.Nq+q] = obj.settings.IC(obj.x[end],obj.q.xi[k],obj.q.xi[q])[1];
        end
    end

    prog = Progress(Nt,1)

    # time loop
    #@gif
    for n = 1:Nt

        ################## K-step ##################
        K .= X*S;

        # impose BCs
        WQuad = EvalAtQuad(obj.basis,W)';
        for i = 1:r
            K[1,i] = Integral(obj.q,obj.uL.*WQuad[i,:].*fXi);
            K[end,i] = Integral(obj.q,obj.uR.*WQuad[i,:].*fXi);
        end

        if obj.settings.rkType == "Heun"
            K1 .= K .+ dt*RhsK(obj,K,W);
            K1 .= K1 .+ dt*RhsK(obj,K1,W);
            K .= 0.5.*(K.+K1);
        elseif obj.settings.rkType == "Euler"
            K .= K .+ dt*RhsK(obj,K,W);
        elseif obj.settings.rkType == "SSP"
            K .= UpdateK(obj.rkUpdate,obj,K,W);
        end

        XNew,STmp = qrGramSchmidt(K);

        MUp = XNew' * X;

        ################## L-step ##################
        L .= W*S';

        if obj.settings.rkType == "Heun"
            L1 .= L .+ dt*RhsL(obj,X,L)';
            L1 .= L1 .+ dt*RhsL(obj,X,L1)';
            L .= 0.5.*(L.+L1);
        elseif obj.settings.rkType == "Euler"
            L .= L .+ dt*RhsL(obj,X,L)';
        elseif obj.settings.rkType == "SSP"
            L .= UpdateL(obj.rkUpdate,obj,X,L);
        end

        WNew,STmp = qrGramSchmidt(L);

        NUp = WNew' * W;

        W .= WNew;
        X .= XNew;

        ################## S-step ##################
        S .= MUp*S*(NUp')

        L = W*S';

        if obj.settings.rkType == "Heun"
            S1 .= S .+ dt.*RhsS(obj,X,W,L);
            L .= W*S';
            S1 .= S1 .+ dt.*RhsS(obj,X,W,L);
            S .= 0.5.*(S.+S1);
        elseif obj.settings.rkType == "Euler"
            S .= S .+ dt.*RhsS(obj,X,W,L);
        elseif obj.settings.rkType == "SSP"
            S .= UpdateS(obj.rkUpdate,obj,X,S,W,false);
        end

        # apply filter step
        Filter(obj,W);

        next!(prog) # update progress bar

        t = t+dt;
    end

    #println(norm(uHat))

    # return end time and solution
    return t, X,S,W;
end

function RhsNodal(obj::DLRSolver,uQ::Array{Float64,2})
    rhs = zeros(size(uQ))
    Nx = obj.settings.Nx;
    for j = 2:(Nx-1)
        for k = 1:obj.settings.Nq^2
            rhs[j,k] = 0.5*(uQ[j+1,k]-2*uQ[j,k]+uQ[j-1,k])/obj.settings.dt - 0.25/obj.settings.dx*(uQ[j+1,k]^2-uQ[j-1,k]^2);
        end
    end
    return rhs;
end

function RhsNodalSplit(obj::DLRSolver,uQ::Array{Float64,2})
    rhs = zeros(size(uQ))
    Nx = obj.settings.Nx;

    for j = 2:(Nx-1)
        for k = 1:obj.settings.Nq^2
            rhs[j,k] = 0.5*(uQ[j+1,k]-2*uQ[j,k]+uQ[j-1,k])/obj.settings.dt - 0.25/obj.settings.dx*(uQ[j+1,k]^2-uQ[j-1,k]^2);
        end
    end
    return rhs;
end

function SolveNaiveSplit(obj::DLRSolver)
    t = 0.0;
    dt = obj.settings.dt;
    dx = obj.x[2]-obj.x[1];
    tEnd = obj.settings.tEnd;
    Nx = obj.settings.Nx;
    Nq = obj.settings.Nq^2;
    N = obj.settings.N^2; # here, N is the number of quadrature points

    fXi = 0.25;

    # save basis
    PhiQuad = Array(obj.basis.PhiQuad');
    PhiQuadCons = Array(obj.basis.PhiQuadCons');
    PhiQuadW = Array(obj.basis.PhiQuadW');
    PhiQuadWCons = Array(obj.basis.PhiQuadWCons');

    # Set up initial condition
    u, uCons = SetupIC(obj);
    u = Array(u');
    uQ = u*PhiQuad;
    #uNew = u;

    Nt = Integer(round(tEnd/dt));

    # compute Dirichlet values if they are independent of time
    for k = 1:obj.settings.Nq
        for q = 1:obj.settings.Nq
            obj.uL[(k-1)*obj.settings.Nq+q] = obj.settings.IC(obj.x[1],obj.q.xi[k],obj.q.xi[q])[1];
            obj.uR[(k-1)*obj.settings.Nq+q] = obj.settings.IC(obj.x[end],obj.q.xi[k],obj.q.xi[q])[1];
        end
    end

    prog = Progress(Nt,1)

    # time loop
    for n = 1:Nt

        uQ = u*PhiQuad + uCons*PhiQuadCons;

        println(maximum(dt*RhsNodal(obj,uQ)))

        #uQ = uQ .+ dt*RhsNodal(obj,uQ);

        uCons = uCons .+ dt*RhsNodalSplit(obj,uQ)*fXi*PhiQuadWCons;

        u = u .+ dt*RhsNodalSplit(obj,uQ)*fXi*PhiQuadW;

        next!(prog) # update progress bar

        t = t+dt;
    end

    # compute moments
    uQ = u*PhiQuad + uCons*PhiQuadCons;
    u = uQ*Array(obj.basis.PhiQuadWFull')*fXi;

    # return end time and solution
    return t, u;
    #return t, uQ*fXi*PhiQuadW;
end

function SolveNaiveSplitUnconventionalIntegrator(obj::DLRSolver)
    t = 0.0;
    dt = obj.settings.dt;
    dx = obj.x[2]-obj.x[1];
    tEnd = obj.settings.tEnd;
    Nx = obj.settings.Nx;
    Nq = obj.settings.Nq^2;
    N = obj.settings.N^2; # here, N is the number of quadrature points
    r = obj.settings.r;

    fXi = 0.25;

    # save basis
    PhiQuad = Array(obj.basis.PhiQuad');
    PhiQuadCons = Array(obj.basis.PhiQuadCons');
    PhiQuadW = Array(obj.basis.PhiQuadW');
    PhiQuadWCons = Array(obj.basis.PhiQuadWCons');

    # Set up initial condition
    u, uCons = SetupIC(obj);
    u = Array(u');
    uQ = u*PhiQuad;

    # Low-rank approx of init data:
    X,S,W = svd(u);

    # rank-r truncation:
    X = X[:,1:r];
    W = W[:,1:r];
    S = Array(Diagonal(S));
    S = S[1:r, 1:r];

    K = zeros(Nx,r);
    K1 = zeros(Nx,r);
    KNew = X*S;
    L = zeros(N,r);
    L1 = zeros(N,r);
    S1 = zeros(r,r);

    Nt = Integer(round(tEnd/dt));

    # compute Dirichlet values if they are independent of time
    for k = 1:obj.settings.Nq
        for q = 1:obj.settings.Nq
            obj.uL[(k-1)*obj.settings.Nq+q] = obj.settings.IC(obj.x[1],obj.q.xi[k],obj.q.xi[q])[1];
            obj.uR[(k-1)*obj.settings.Nq+q] = obj.settings.IC(obj.x[end],obj.q.xi[k],obj.q.xi[q])[1];
        end
    end

    prog = Progress(Nt,1)

    # time loop
    for n = 1:Nt

        # reconstruct solution (inefficient)
        uQ = X*S*W'*PhiQuad + uCons*PhiQuadCons;

        # update conservative part (inefficient)
        uCons = uCons .+ dt*RhsNodalSplit(obj,uQ)*fXi*PhiQuadWCons;

        ###### K-step ######
        K .= X*S;
        WQuadW = PhiQuadW*W; # compute W_{ki}*w_k , where W \in R^{Nq x r}

        #K .= K .+ dt*Integral(obj.q,F(obj,K*EvalAtQuadDLR(obj.basis,W))*EvalAtQuadDLR(obj.basis,W).*fXi);
        K .= K .+ dt*RhsNodalSplit(obj,uQ)*WQuadW.*fXi;

        XNew,STmp = qrGramSchmidt(K); # optimize bei choosing XFull, SFull

        MUp = XNew' * X;

        ###### L-step ######
        L = W*S';

        #L .= L .+ dt*(X'*F(obj,X*EvalAtQuadDLR(obj.basis,L)))';
        L .= L .+ dt*(X'*RhsNodalSplit(obj,uQ)*PhiQuadW.*fXi)';

        WNew,STmp = qrGramSchmidt(L);

        NUp = WNew' * W;
        W .= WNew;
        X .= XNew;

        ################## S-step ##################
        S .= MUp*S*(NUp')
        uQ = X*S*W'*PhiQuad + uCons*PhiQuadCons;
        WQuadW = PhiQuadW*W;

        S .= S .+ dt.*X'*RhsNodalSplit(obj,uQ)*WQuadW.*fXi;

        next!(prog) # update progress bar

        t = t+dt;
    end

    # compute moments
    uQ = X*S*W'*PhiQuad + uCons*PhiQuadCons;
    u = uQ*Array(obj.basis.PhiQuadWFull')*fXi;

    # return end time and solution
    return t, u;
end

function SolveSplitUnconventionalIntegrator(obj::DLRSolver,X::Array{Float64,2},S::Array{Float64,2},W::Array{Float64,2},uCons::Array{Float64,2})
    t = 0.0;
    dt = obj.settings.dt;
    dx = obj.x[2]-obj.x[1];
    tEnd = obj.settings.tEnd;
    Nx = obj.settings.Nx;
    N = obj.settings.N^2; # here, N is the number of quadrature points
    r = obj.settings.r;
    NCons = obj.settings.NCons
    w = obj.q.wTens;

    fXi = 0.25;

    # save basis
    PhiQuad = Array(obj.basis.PhiQuad');
    PhiQuadCons = Array(obj.basis.PhiQuadCons');
    PhiQuadW = Array(obj.basis.PhiQuadW');

    K = zeros(Nx,r);
    L = zeros(N,r);

    Nt = Integer(round(tEnd/dt));

    # compute Dirichlet values if they are independent of time
    for k = 1:obj.settings.Nq
        for q = 1:obj.settings.Nq
            obj.uL[(k-1)*obj.settings.Nq+q] = obj.settings.IC(obj.x[1],obj.q.xi[k],obj.q.xi[q])[1];
            obj.uR[(k-1)*obj.settings.Nq+q] = obj.settings.IC(obj.x[end],obj.q.xi[k],obj.q.xi[q])[1];
        end
    end

    prog = Progress(Nt,1)

    Y = zeros(r,r,r);
    YCons = zeros(NCons,NCons,r);
    YCross = zeros(NCons,r,r);

    Y1 = zeros(r,r,r);
    Y1Cons = zeros(NCons,NCons,r);
    Y1Cross = zeros(NCons,r,r);
    Y1ConsFull = zeros(NCons,NCons,NCons);

    fluxK1 = zeros(Nx,r)
    fluxK2 = zeros(Nx,r)
    fluxK3 = zeros(Nx,r)

    fluxL1 = zeros(r,N^2)
    fluxL2 = zeros(r,N^2)
    fluxL3 = zeros(r,N^2)

    fluxS1 = zeros(r,N^2)
    fluxS2 = zeros(r,N^2)
    fluxS3 = zeros(r,N^2)

    fluxCons1 = zeros(Nx,NCons)
    fluxCons2 = zeros(Nx,NCons)
    fluxCons3 = zeros(Nx,NCons)

    # time loop
    for n = 1:Nt

        # reconstruct solution (inefficient)
        WQuad = PhiQuad'*W; # compute W_{ki}*w_k , where W \in R^{Nq x r}
        K .= X*S;

        # update conservative part (inefficient)
        @einsum Y1Cross[l,m,i] := WQuad[k,m]*WQuad[k,i]*PhiQuadCons[l,k]*w[k].*fXi
        @einsum Y1Cons[l,m,i] := PhiQuadCons[m,k]*WQuad[k,i]*PhiQuadCons[l,k]*w[k].*fXi
        @einsum Y1ConsFull[l,m,i] := PhiQuadCons[m,k]*PhiQuadCons[i,k]*PhiQuadCons[l,k]*w[k].*fXi

        @einsum fluxCons1[j,l] := K[j,m]*Y1Cross[l,m,i]*K[j,i];
        @einsum fluxCons2[j,l] := K[j,i]*Y1Cons[l,m,i]*uCons[j,m];
        @einsum fluxCons3[j,l] := uCons[j,m]*Y1ConsFull[l,m,i]*uCons[j,i];

        uConsNew = uCons .+ dt*(obj.L1I*uCons .- obj.L2*(fluxCons1+2.0*fluxCons2+fluxCons3));

        ###### K-step ######
        @einsum Y[l,m,i] := WQuad[k,l]*WQuad[k,m]*WQuad[k,i]*w[k].*fXi
        @einsum YCross[l,m,i] := PhiQuadCons[l,k]*WQuad[k,m]*WQuad[k,i]*w[k].*fXi
        @einsum YCons[l,m,i] := PhiQuadCons[l,k]*PhiQuadCons[m,k]*WQuad[k,i]*w[k].*fXi

        # value uQ^2*WQuadW.*fXi
        @einsum fluxK1[j,i] := K[j,l]*Y[l,m,i]*K[j,m]
        @einsum fluxK2[j,i] := uCons[j,l]*YCross[l,m,i]*K[j,m]
        @einsum fluxK3[j,i] := uCons[j,l]*YCons[l,m,i]*uCons[j,m]

        K .= K .+ dt*(obj.L1I*K .- obj.L2*(fluxK1.+2*fluxK2.+fluxK3));

        XNew,STmp = qrGramSchmidt(K); # optimize bei choosing XFull, SFull

        MUp = XNew' * X;

        ###### L-step ######
        L = W*S';
        LQuad = PhiQuad'*L;

        XL1X = X'*obj.L1I*X;
        XL2 = X'*obj.L2

        @einsum Y[l,m,i] := XL2[i,j]*X[j,l]*X[j,m];
        @einsum YCross[l,m,i] := XL2[i,j]*uCons[j,l]*X[j,m];
        @einsum YCons[l,m,i] := XL2[i,j]*uCons[j,l]*uCons[j,m];

        @einsum fluxL1[i,k] := Y[l,m,i]*LQuad[k,l]*LQuad[k,m]
        @einsum fluxL2[i,k] := YCross[l,m,i]*PhiQuadCons[l,k]*LQuad[k,m]
        @einsum fluxL3[i,k] := YCons[l,m,i]*PhiQuadCons[l,k]*PhiQuadCons[m,k]

        L .= L .+ dt*(XL1X*L' .- (fluxL1 .+ 2*fluxL2 .+ fluxL3)*PhiQuadW.*fXi)';

        WNew,STmp = qrGramSchmidt(L);

        NUp = WNew' * W;
        W .= WNew;
        X .= XNew;

        ################## S-step ##################
        S .= MUp*S*(NUp')

        WQuad = PhiQuad'*W;
        K .= X*S;

        @einsum Y[l,m,i] := WQuad[k,l]*WQuad[k,m]*WQuad[k,i]*w[k].*fXi
        @einsum YCross[l,m,i] := PhiQuadCons[l,k]*WQuad[k,m]*WQuad[k,i]*w[k].*fXi
        @einsum YCons[l,m,i] := PhiQuadCons[l,k]*PhiQuadCons[m,k]*WQuad[k,i]*w[k].*fXi

        XL2 = X'*obj.L2
        @einsum Y1[l,m,i] := XL2[i,j]*K[j,l]*K[j,m];
        @einsum Y1Cross[l,m,i] := XL2[i,j]*uConsNew[j,l]*K[j,m];
        @einsum Y1Cons[l,m,i] := XL2[i,j]*uConsNew[j,l]*uConsNew[j,m];

        @einsum fluxS1[j,i] := Y[l,m,i]*Y1[l,m,j];
        @einsum fluxS2[j,i] := YCross[l,m,i]*Y1Cross[l,m,j];
        @einsum fluxS3[j,i] := YCons[l,m,i]*Y1Cons[l,m,j];

        S .= S .+ dt.*(XL1X*S .- (fluxS1.+2*fluxS2.+fluxS3));

        uCons .= uConsNew

        next!(prog) # update progress bar

        t = t+dt;
    end

    # return end time and solution
    return t, X,S,W,uCons;
end

function SolveSplitProjectorSplittingIntegrator(obj::DLRSolver,X::Array{Float64,2},S::Array{Float64,2},W::Array{Float64,2},uCons::Array{Float64,2})
    t = 0.0;
    dt = obj.settings.dt;
    dx = obj.x[2]-obj.x[1];
    tEnd = obj.settings.tEnd;
    Nx = obj.settings.Nx;
    N = obj.settings.N^2; # here, N is the number of quadrature points
    r = obj.settings.r;
    NCons = obj.settings.NCons
    w = obj.q.wTens;

    fXi = 0.25;

    # save basis
    PhiQuad = Array(obj.basis.PhiQuad');
    PhiQuadCons = Array(obj.basis.PhiQuadCons');
    PhiQuadW = Array(obj.basis.PhiQuadW');

    K = zeros(Nx,r);
    L = zeros(N,r);

    Nt = Integer(round(tEnd/dt));

    # compute Dirichlet values if they are independent of time
    for k = 1:obj.settings.Nq
        for q = 1:obj.settings.Nq
            obj.uL[(k-1)*obj.settings.Nq+q] = obj.settings.IC(obj.x[1],obj.q.xi[k],obj.q.xi[q])[1];
            obj.uR[(k-1)*obj.settings.Nq+q] = obj.settings.IC(obj.x[end],obj.q.xi[k],obj.q.xi[q])[1];
        end
    end

    prog = Progress(Nt,1)

    Y = zeros(r,r,r);
    YCons = zeros(NCons,NCons,r);
    YCross = zeros(NCons,r,r);

    Y1 = zeros(r,r,r);
    Y1Cons = zeros(NCons,NCons,r);
    Y1Cross = zeros(NCons,r,r);
    Y1ConsFull = zeros(NCons,NCons,NCons);

    fluxK1 = zeros(Nx,r)
    fluxK2 = zeros(Nx,r)
    fluxK3 = zeros(Nx,r)

    fluxL1 = zeros(r,N^2)
    fluxL2 = zeros(r,N^2)
    fluxL3 = zeros(r,N^2)

    fluxS1 = zeros(r,N^2)
    fluxS2 = zeros(r,N^2)
    fluxS3 = zeros(r,N^2)

    fluxCons1 = zeros(Nx,NCons)
    fluxCons2 = zeros(Nx,NCons)
    fluxCons3 = zeros(Nx,NCons)

    # time loop
    for n = 1:Nt

        # reconstruct solution (inefficient)
        WQuad = PhiQuad'*W; # compute W_{ki}*w_k , where W \in R^{Nq x r}
        K .= X*S;

        # update conservative part (inefficient)
        @einsum Y1Cross[l,m,i] := WQuad[k,m]*WQuad[k,i]*PhiQuadCons[l,k]*w[k].*fXi
        @einsum Y1Cons[l,m,i] := PhiQuadCons[m,k]*WQuad[k,i]*PhiQuadCons[l,k]*w[k].*fXi
        @einsum Y1ConsFull[l,m,i] := PhiQuadCons[m,k]*PhiQuadCons[i,k]*PhiQuadCons[l,k]*w[k].*fXi

        @einsum fluxCons1[j,l] := K[j,m]*Y1Cross[l,m,i]*K[j,i];
        @einsum fluxCons2[j,l] := K[j,i]*Y1Cons[l,m,i]*uCons[j,m];
        @einsum fluxCons3[j,l] := uCons[j,m]*Y1ConsFull[l,m,i]*uCons[j,i];

        uConsNew = uCons .+ dt*(obj.L1I*uCons .- obj.L2*(fluxCons1+2.0*fluxCons2+fluxCons3));

        ###### K-step ######
        @einsum Y[l,m,i] := WQuad[k,l]*WQuad[k,m]*WQuad[k,i]*w[k].*fXi
        @einsum YCross[l,m,i] := PhiQuadCons[l,k]*WQuad[k,m]*WQuad[k,i]*w[k].*fXi
        @einsum YCons[l,m,i] := PhiQuadCons[l,k]*PhiQuadCons[m,k]*WQuad[k,i]*w[k].*fXi

        # value uQ^2*WQuadW.*fXi
        @einsum fluxK1[j,i] := K[j,l]*Y[l,m,i]*K[j,m]
        @einsum fluxK2[j,i] := uCons[j,l]*YCross[l,m,i]*K[j,m]
        @einsum fluxK3[j,i] := uCons[j,l]*YCons[l,m,i]*uCons[j,m]

        K .= K .+ dt*(obj.L1I*K .- obj.L2*(fluxK1.+2*fluxK2.+fluxK3));

        X,S = qrGramSchmidt(K); # optimize bei choosing XFull, SFull

        ################## S-step ##################

        @einsum Y[l,m,i] := WQuad[k,l]*WQuad[k,m]*WQuad[k,i]*w[k].*fXi
        @einsum YCross[l,m,i] := PhiQuadCons[l,k]*WQuad[k,m]*WQuad[k,i]*w[k].*fXi
        @einsum YCons[l,m,i] := PhiQuadCons[l,k]*PhiQuadCons[m,k]*WQuad[k,i]*w[k].*fXi

        XL1X = X'*obj.L1I*X;
        XL2 = X'*obj.L2
        @einsum Y1[l,m,i] := XL2[i,j]*K[j,l]*K[j,m];
        @einsum Y1Cross[l,m,i] := XL2[i,j]*uConsNew[j,l]*K[j,m];
        @einsum Y1Cons[l,m,i] := XL2[i,j]*uConsNew[j,l]*uConsNew[j,m];

        @einsum fluxS1[j,i] := Y[l,m,i]*Y1[l,m,j];
        @einsum fluxS2[j,i] := YCross[l,m,i]*Y1Cross[l,m,j];
        @einsum fluxS3[j,i] := YCons[l,m,i]*Y1Cons[l,m,j];

        S .= S .- dt.*(XL1X*S .- (fluxS1.+2*fluxS2.+fluxS3));

        ###### L-step ######
        L = W*S';
        LQuad = PhiQuad'*L;

        @einsum Y[l,m,i] := XL2[i,j]*X[j,l]*X[j,m];
        @einsum YCross[l,m,i] := XL2[i,j]*uCons[j,l]*X[j,m];
        @einsum YCons[l,m,i] := XL2[i,j]*uCons[j,l]*uCons[j,m];

        @einsum fluxL1[i,k] := Y[l,m,i]*LQuad[k,l]*LQuad[k,m]
        @einsum fluxL2[i,k] := YCross[l,m,i]*PhiQuadCons[l,k]*LQuad[k,m]
        @einsum fluxL3[i,k] := YCons[l,m,i]*PhiQuadCons[l,k]*PhiQuadCons[m,k]

        L .= L .+ dt*(XL1X*L' .- (fluxL1 .+ 2*fluxL2 .+ fluxL3)*PhiQuadW.*fXi)';

        W,S = qrGramSchmidt(L);
        S .= S';

        uCons .= uConsNew

        next!(prog) # update progress bar

        t = t+dt;
    end

    # return end time and solution
    return t, X,S,W,uCons;
end

function SolveUnconventionalIntegrator(obj::DLRSolver,X::Array{Float64,2},S::Array{Float64,2},W::Array{Float64,2})
    t = 0.0;
    dt = obj.settings.dt;
    tEnd = obj.settings.tEnd;
    Nx = obj.settings.Nx;
    N = obj.settings.N^2; # here, N is the number of quadrature points
    r = obj.settings.r;
    w = obj.q.wTens;

    fXi = 0.25;

    K = zeros(Nx,r);
    L = zeros(N,r);

    Nt = Integer(round(tEnd/dt));

    # compute Dirichlet values if they are independent of time
    for k = 1:obj.settings.Nq
        for q = 1:obj.settings.Nq
            obj.uL[(k-1)*obj.settings.Nq+q] = obj.settings.IC(obj.x[1],obj.q.xi[k],obj.q.xi[q])[1];
            obj.uR[(k-1)*obj.settings.Nq+q] = obj.settings.IC(obj.x[end],obj.q.xi[k],obj.q.xi[q])[1];
        end
    end

    prog = Progress(Nt,1)

    Y = zeros(r,r,r);

    Y1 = zeros(r,r,r);

    fluxK1 = zeros(Nx,r)

    fluxL1 = zeros(r,N^2)

    fluxS1 = zeros(r,N^2)

    # time loop
    for n = 1:Nt

        WQuad = obj.basis.PhiQuad*W; # compute W_{ki}*w_k , where W \in R^{Nq x r}
        K .= X*S;

        ###### K-step ######
        @einsum Y[l,m,i] := WQuad[k,l]*WQuad[k,m]*WQuad[k,i]*w[k].*fXi

        # value uQ^2*WQuadW.*fXi
        @einsum fluxK1[j,i] := K[j,l]*Y[l,m,i]*K[j,m]

        K .= K .+ dt*(obj.L1I*K .- obj.L2*fluxK1);

        XNew,STmp = qrGramSchmidt(K); # optimize bei choosing XFull, SFull

        MUp = XNew' * X;

        ###### L-step ######
        L = W*S';
        LQuad = obj.basis.PhiQuad*L;

        XL1X = X'*obj.L1I*X;
        XL2 = X'*obj.L2

        @einsum Y[l,m,i] := XL2[i,j]*X[j,l]*X[j,m];

        @einsum fluxL1[i,k] := Y[l,m,i]*LQuad[k,l]*LQuad[k,m]

        L .= L .+ dt*(XL1X*L' .- fluxL1*obj.basis.PhiQuadW'.*fXi)';

        WNew,STmp = qrGramSchmidt(L);

        NUp = WNew' * W;
        W .= WNew;
        X .= XNew;

        ################## S-step ##################
        S .= MUp*S*(NUp')

        WQuad = obj.basis.PhiQuad*W;
        K .= X*S;

        @einsum Y[l,m,i] := WQuad[k,l]*WQuad[k,m]*WQuad[k,i]*w[k].*fXi

        XL2 = X'*obj.L2
        @einsum Y1[l,m,i] := XL2[i,j]*K[j,l]*K[j,m];

        @einsum fluxS1[j,i] := Y[l,m,i]*Y1[l,m,j];

        S .= S .+ dt.*(XL1X*S .- fluxS1);

        next!(prog) # update progress bar

        t = t+dt;
    end

    # return end time and solution
    return t, X,S,W;
end

function SolveProjectorSplittingIntegrator(obj::DLRSolver,X::Array{Float64,2},S::Array{Float64,2},W::Array{Float64,2})
    t = 0.0;
    dt = obj.settings.dt;
    tEnd = obj.settings.tEnd;
    Nx = obj.settings.Nx;
    N = obj.settings.N^2; # here, N is the number of quadrature points
    r = obj.settings.r;
    w = obj.q.wTens;

    fXi = 0.25;

    K = zeros(Nx,r);
    L = zeros(N,r);

    Nt = Integer(round(tEnd/dt));

    # compute Dirichlet values if they are independent of time
    for k = 1:obj.settings.Nq
        for q = 1:obj.settings.Nq
            obj.uL[(k-1)*obj.settings.Nq+q] = obj.settings.IC(obj.x[1],obj.q.xi[k],obj.q.xi[q])[1];
            obj.uR[(k-1)*obj.settings.Nq+q] = obj.settings.IC(obj.x[end],obj.q.xi[k],obj.q.xi[q])[1];
        end
    end

    prog = Progress(Nt,1)

    Y = zeros(r,r,r);

    Y1 = zeros(r,r,r);

    fluxK1 = zeros(Nx,r)

    fluxL1 = zeros(r,N^2)

    fluxS1 = zeros(r,N^2)

    # time loop
    for n = 1:Nt

        WQuad = obj.basis.PhiQuad*W; # compute W_{ki}*w_k , where W \in R^{Nq x r}
        K .= X*S;

        # impose BCs
        for i = 1:r
            K[1,i] = Integral(obj.q,obj.uL.*WQuad[:,i].*fXi);
            K[end,i] = Integral(obj.q,obj.uR.*WQuad[:,i].*fXi);
        end

        ###### K-step ######
        @einsum Y[l,m,i] := WQuad[k,l]*WQuad[k,m]*WQuad[k,i]*w[k].*fXi

        # value uQ^2*WQuadW.*fXi
        @einsum fluxK1[j,i] := K[j,l]*Y[l,m,i]*K[j,m]

        K .= K .+ dt*(obj.L1I*K .- obj.L2*fluxK1);

        X,S = qrGramSchmidt(K); # optimize bei choosing XFull, SFull

        ################## S-step ##################

        @einsum Y[l,m,i] := WQuad[k,l]*WQuad[k,m]*WQuad[k,i]*w[k].*fXi

        XL1X = X'*obj.L1I*X;
        XL2 = X'*obj.L2
        @einsum Y1[l,m,i] := XL2[i,j]*K[j,l]*K[j,m];

        @einsum fluxS1[j,i] := Y[l,m,i]*Y1[l,m,j];

        S .= S .- dt.*(XL1X*S .- fluxS1);

        ###### L-step ######
        L = W*S';
        LQuad = obj.basis.PhiQuad*L;

        @einsum Y[l,m,i] := XL2[i,j]*X[j,l]*X[j,m];

        @einsum fluxL1[i,k] := Y[l,m,i]*LQuad[k,l]*LQuad[k,m]

        L .= L .+ dt*(XL1X*L' .- fluxL1*obj.basis.PhiQuadW'.*fXi)';

        W,S = qrGramSchmidt(L);
        S .= S';

        next!(prog) # update progress bar

        t = t+dt;
    end

    # return end time and solution
    return t, X,S,W;
end

function Filter(obj::DLRSolver,W::Array{Float64,2})
    lambda = obj.settings.lambda
    N = obj.settings.N;
    if obj.settings.filterType == "L2"
        for i = 1:N
            for l = 1:N
                W[(l-1)*N+i,:] .= W[(l-1)*N+i,:]/(1+lambda*(i-1)^2*i^2+lambda*(l-1)^2*l^2);
            end
        end
    elseif obj.settings.filterType == "EXP"
        epsilonM = eps(Float64);
        c = log( epsilonM );
        for i = 1:size(W,1)
            eta = i/(obj.settings.N+1)
            W[i,:] .= W[i,:]*exp( c * eta^obj.settings.filterOrder )^(lambda*obj.settings.dt);
        end
    end
end

# SSP Update function for K-step
function UpdateK(obj::TimeSolver,solver::DLRSolver,K::Array{Float64,2},W::Array{Float64,2})
    NCells = solver.settings.Nx;
    obj.KRK[1,:,:] .= K;
    for s = 1:obj.rkStages
        obj.Krhs[s,:,:] .= RhsK(solver,obj.KRK[s,:,:],W);;
        obj.KRK[s+1,:,:] .= zeros(NCells,solver.settings.r);
        for j = 1:s
            obj.KRK[s+1,:,:] = obj.KRK[s+1,:,:]+obj.alpha[s,j].*obj.KRK[j,:,:]+obj.dt*obj.beta[s,j].*obj.Krhs[j,:,:];
        end
    end

    return obj.KRK[obj.rkStages+1,:,:];
end

# SSP Update function for S-step
function UpdateS(obj::TimeSolver,solver::DLRSolver,X::Array{Float64,2},S::Array{Float64,2},W::Array{Float64,2},backward::Bool=true)
    obj.SRK[1,:,:] .= S;
    for s = 1:obj.rkStages
        obj.Srhs[s,:,:] .= RhsS(solver,X,W,W*obj.SRK[s,:,:]');
        if backward # if projector splitting integrator is used, sign in rhsS must be changed
            obj.Srhs[s,:,:] .= -obj.Srhs[s,:,:];
        end
        obj.SRK[s+1,:,:] .= zeros(solver.settings.r,solver.settings.r);
        for j = 1:s
            obj.SRK[s+1,:,:] = obj.SRK[s+1,:,:]+obj.alpha[s,j].*obj.SRK[j,:,:]+obj.dt*obj.beta[s,j].*obj.Srhs[j,:,:];
        end
    end

    return obj.SRK[obj.rkStages+1,:,:];
end

# SSP Update function for L-step
function UpdateL(obj::TimeSolver,solver::DLRSolver,X::Array{Float64,2},L::Array{Float64,2},recompute::Bool=true)
    N = solver.settings.N^2;
    obj.LRK[1,:,:] .= L;
    for s = 1:obj.rkStages
        obj.Lrhs[s,:,:] .= RhsL(solver,X,obj.LRK[s,:,:],recompute)';
        obj.LRK[s+1,:,:] .= zeros(N,solver.settings.r);
        for j = 1:s
            obj.LRK[s+1,:,:] = obj.LRK[s+1,:,:]+obj.alpha[s,j].*obj.LRK[j,:,:]+obj.dt*obj.beta[s,j].*obj.Lrhs[j,:,:];
        end
    end

    return obj.LRK[obj.rkStages+1,:,:];
end


text:
# load required modules
using DataFrames
using GeoJSON
using LeafletJS
using Random
using RDatasets

# set seed for reproducible sequence of numbers
Random.seed!(42);

# load geojson data
nbhd = GeoJSON.parsefile("../../data/raw/neighborhoods.geojson")
rodents = GeoJSON.parsefile("../../data/raw/rodents.geojson");

# create layers
layers = Vector{Layer}(undef,0)
push!(layers, Layer(nbhd))
push!(layers, Layer(rodents));

# render map with stament default tile
LeafletMap(layers, [42.35, -71.08])

# render map with stamen watercolor tile
LeafletMap(layers, [42.35, -71.08], provider = Provider.Stamen(:watercolor))

# render map with stamen terrain tile
LeafletMap(layers, [42.35, -71.08], provider = Provider.Stamen(:terrain))

# render map with stamen toposm-features tile
LeafletMap(layers, [42.35, -71.08], provider = Provider.Stamen(Symbol("toposm-features")))

# render map with OSM mapnok tile
LeafletMap(layers, [42.35, -71.08], provider = Provider.OSM())

# render map with OSM black and white tile
LeafletMap(layers, [42.35, -71.08], provider = Provider.OSMBlackAndWhite())

# render map with OSM HOT tile
LeafletMap(layers, [42.35, -71.08], provider = Provider.OSMHOT())

# render map with CARTO light tile
LeafletMap(layers, [42.35, -71.08], provider = Provider.CARTO())

# render map with CARTO dark tile
LeafletMap(layers, [42.35, -71.08], provider = Provider.CARTO(:dark_all))

# render map with NASA GIBS light MODIS Terra tile
LeafletMap(layers, [42.35, -71.08], provider=Provider.NASAGIBS())

ext: jl

max_stars_repo_name: lykmapipo/data-science-learning

max_issues_repo_name: lykmapipo/data-science-learning

max_forks_repo_name: lykmapipo/data-science-learning


text:
function time_derivative(f::Function)
  function time_derivative_f(x,t)
    fxt = zero(return_type(f,x,t))
    _time_derivative_f(f,x,t,fxt)
  end
  time_derivative_f(x::VectorValue) = t -> time_derivative_f(x,t)
  time_derivative_f(t) = x -> time_derivative_f(x,t)
end

const ∂t = time_derivative

function _time_derivative_f(f,x,t,fxt)
  ForwardDiff.derivative(t->f(x,t),t)
end

function _time_derivative_f(f,x,t,fxt::VectorValue)
  VectorValue(ForwardDiff.derivative(t->get_array(f(x,t)),t))
  # VectorValue(ForwardDiff.derivative(t->f(x,t),t))
end

function _time_derivative_f(f,x,t,fxt::TensorValue)
  TensorValue(ForwardDiff.derivative(t->get_array(f(x,t)),t))
end


text:
using LinearAlgebra.BLAS

@testset "BLAS" begin
    let rng = MersenneTwister(123456)
        for _ in 1:10
            x, y, vx, vy = randn.(Ref(rng), [10, 6, 10, 6])
            _dot = (x, y)->BLAS.dot(5, x, 2, y, 1)
            @test check_errs(_dot, _dot(x, y), (x, y), (vx, vy))
        end
    end

    let rng = MersenneTwister(123456)
        for _ in 1:10
            x, vx = randn(rng, 100), randn(rng, 100)
            @test check_errs(nrm2, randn(rng), x, vx)
        end
    end
    let rng = MersenneTwister(123456)
        λ = x->nrm2(50, x, 2)
        for _ in 1:10
            x, vx = randn(rng, 100), randn(rng, 100)
            @test check_errs(λ, randn(rng), x, vx)
        end
    end

    let rng = MersenneTwister(123456)
        λ = x->asum(50, x, 2)
        for _ in 1:10
            x, vx = randn(rng, 100), randn(rng, 100)
            x .+= sign.(x)
            @test check_errs(asum, randn(rng), x, vx)
            @test check_errs(λ, randn(rng), x, vx)
        end
    end

    # Test each of the four permutations of `gemm`.
    let rng = MersenneTwister(123456), N = 100
        for tA in ['T', 'N'], tB in ['T', 'N']
            λ, γ = (α, A, B)->gemm(tA, tB, α, A, B), (A, B)->gemm(tA, tB, A, B)
            for _ in 1:10
                α, vα = randn.([rng, rng])
                A, B, VA, VB = randn.(Ref(rng), [N, N, N, N], [N, N, N, N])
                @test check_errs(λ, λ(α, A, B), (α, A, B), (vα, VA, VB))
                @test check_errs(γ, γ(A, B), (A, B), (VA, VB))
            end
        end
    end

    # Test both permutations of `gemv`.
    let rng = MersenneTwister(123456), N = 100
        for tA in ['T', 'N']
            λ, γ = (α, A, x)->gemv('T', α, A, x), (A, x)->gemv('T', A, x)
            for _ in 1:10
                α, vα = randn.([rng, rng])
                A, VA = randn.(Ref(rng), [N, N], [N, N])
                x, vx = randn.(Ref(rng), [N, N])
                @test check_errs(λ, λ(α, A, x), (α, A, x), (vα, VA, vx))
                @test check_errs(γ, γ(A, x), (A, x), (VA, vx))
            end
        end
    end

    # Test all four permutations of `symm`.
    let rng = MersenneTwister(123456), N = 100
        lmask, umask = Matrix(LowerTriangular(ones(N, N))), Matrix(UpperTriangular(ones(N, N)))
        for side in ['L', 'R'], ul in ['L', 'U']
            λ, γ = (α, A, B)->symm(side, ul, α, A, B), (A, B)->symm(side, ul, A, B)
            for _ in 1:10
                α, vα = randn.([rng, rng])
                A, B, VA, VB = randn.(Ref(rng), [N, N, N, N], [N, N, N, N])
                @test check_errs(λ, λ(α, A, B), (α, A, B), (vα, VA, VB))
                @test check_errs(γ, γ(A, B), (A, B), (VA, VB))
            end
        end
    end

    let rng = MersenneTwister(123456), N = 100
        for ul in ['L', 'U']
            λ, γ = (α, A, x)->symv(ul, α, A, x), (A, x)->symv(ul, A, x)
            for _ in 1:10
                α, vα = randn.([rng, rng])
                A, VA = randn.(Ref(rng), [N, N], [N, N])
                x, vx = randn.(Ref(rng), [N, N])
                @test check_errs(λ, λ(α, A, x), (α, A, x), (vα, VA, vx))
                @test check_errs(γ, γ(A, x), (A, x), (VA, vx))
            end
        end
    end

    let rng = MersenneTwister(123456), N = 10
        for side in ['L', 'R'], ul in ['L', 'U'], tA in ['N', 'T'], dA in ['U', 'N']
            λ = (α, A, B)->trmm(side, ul, tA, dA, α, A, B)
            for _ in 1:10
                α, vα = randn.([rng, rng])
                A, B, VA, VB = randn.(Ref(rng), [N, N, N, N], [N, N, N, N])
                @test check_errs(λ, λ(α, A, B), (α, A, B), (vα, VA, VB))
            end
        end
    end

    let rng = MersenneTwister(123456), N = 10
        for ul in ['L', 'U'], tA in ['N', 'T'], dA in ['U', 'N']
            λ = (A, b)->trmv(ul, tA, dA, A, b)
            for _ in 1:10
                A, VA = randn.(Ref(rng), [N, N], [N, N])
                b, vb = randn.(Ref(rng), [N, N])
                @test check_errs(λ, λ(A, b), (A, b), (VA, vb))
            end
        end
    end

    let rng = MersenneTwister(123456), N = 10
        for side in ['L', 'R'], ul in ['L', 'U'], tA in ['N', 'T'], dA in ['U', 'N']
            λ = (α, A, X)->trsm(side, ul, tA, dA, α, A, X)
            for _ in 1:10
                α, vα = randn.([rng, rng])
                A, X, VA, VX = randn.(Ref(rng), [N, N, N, N], [N, N, N, N])
                A = randn(rng, N, N) + UniformScaling(3)
                @test check_errs(λ, λ(α, A, X), (α, A, X), (vα, VA, VX))
            end
        end
    end

    let rng = MersenneTwister(123456), N = 10
        for ul in ['L', 'U'], tA in ['N', 'T'], dA in ['U', 'N']
            λ = (A, x)->trsv(ul, tA, dA, A, x)
            for _ in 1:10
                A = randn(rng, N, N) + UniformScaling(1)
                A = A'A
                VA = randn(rng, N, N)
                x, vx = randn.(Ref(rng), [N, N])
                @test check_errs(λ, λ(A, x), (A, x), (VA, vx))
            end
        end
    end
end


text:
const category = Dict(:A => :Arts,
    :B => :Education,
    :C => :Environmental,
    :D => :AnimalWelfare,
    :E => :Health,
    :F => :MentalHealth,
    :G => :MedicalDiseases,
    :H => :MedicalResearch,
    :I => :CrimeAndLegal,
    :J => :Employment,
    :K => :FoodAndAgriculture,
    :L => :Housing,
    :M => :PublicSafety,
    :N => :Recreation,
    :O => :YouthDevelopment,
    :P => :HumanServices,
    :Q => :ForeignAffairs,
    :R => :SocialAction,
    :S => :CommunityImprovement,
    :T => :GrantmakingFoundations,
    :U => :ScienceandTechResearch,
    :V => :SocialScienceResearch,
    :W => :SocietyBenefitAndMultipurpose,
    :X => :ReligionRelated,
    :Y => :MutualBenefitOrg,
    :Z => :Unknown,
    :ZZ => :AllOthers
)

const charitytype = Dict(:pc => :privatecharity,
    :pf => :privatefoundation,
    :co => :othercharity
)

const affiliationcodes = Dict(
    Symbol(0) => missing,
    Symbol(1) => :CentralNoGroup,
    Symbol(2) => :Intermediate,
    Symbol(3) => :Independent,
    Symbol(6) => :CentralOfGroupNoChurch,
    Symbol(7) => :Intermediate,
    Symbol(8) => :CentralOfGroupChurch,
    Symbol(9) => :Subordinate,
    missing => missing,
    Symbol(0.0) => missing,
    :None => missing,
    Symbol(1.0) => :CentralNoGroup,
    Symbol(2.0) => :Intermediate,
    Symbol(3.0) => :Independent,
    Symbol(6.0) => :CentralOfGroupNoChurch,
    Symbol(7.0) => :Intermediate,
    Symbol(8.0) => :CentralOfGroupChurch,
    Symbol(9.0) => :Subordinate
)

const conformedaffiliationcodes = Dict(
    missing => :noncentral, #Unknown
    :Unknown => :noncentral, #Unknown
    :CentralNoGroup => :central, #CentralNoGroup
    :Intermediate => :noncentral, #Intermediate
    :Independent => :central, #Independent
    :CentralOfGroupNoChurch => :central, #CentralOfGroupNoChurch
    :Intermediate => :noncentral, #Intermediate
    :CentralOfGroupChurch => :central, #CentralOfGroupChurch
    :Subordinate => :noncentral #Subordinate
)

const naicsCodes = Dict(
    Symbol("11") => :Agriculture,
    Symbol("21") => :Mining,
    Symbol("22") => :Utilities,
    Symbol("23") => :Construction,
    Symbol("31") => :Manufacturing,
    Symbol("32") => :Manufacturing,
    Symbol("33") => :Manufacturing,
    Symbol("42") => :WholesaleTrade,
    Symbol("44") => :RetailTrade,
    Symbol("45") => :RetailTrade,
    Symbol("48") => :Transportation,
    Symbol("49") => :Transportation,
    Symbol("51") => :Information,
    Symbol("52") => :Finance,
    Symbol("53") => :RealEstate,
    Symbol("54") => :ProfServices,
    Symbol("55") => :Management,
    Symbol("56") => :AdminAndWaste,
    Symbol("61") => :Education,
    Symbol("62") => :HealthCare,
    Symbol("71") => :Entertainment,
    Symbol("72") => :Hospitality,
    Symbol("81") => :OtherServices,
    Symbol("92") => :PublicAdmin,
    Symbol("99") => :Unknown,
    missing => :Unknown,
    :No => :Unknown)

    #=NAICS unshortened (2017, from US Census)
    11	Agriculture, Forestry, Fishing and Hunting
    21	Mining, Quarrying, and Oil and Gas Extraction
    22	Utilities
    23	Construction
    31-33	Manufacturing
    42	Wholesale Trade
    44-45	Retail Trade
    48-49	Transportation and Warehousing
    51	Information
    52	Finance and Insurance
    53	Real Estate and Rental and Leasing
    54	Professional, Scientific, and Technical Services
    55	Management of Companies and Enterprises
    56	Administrative and Support and Waste Management and Remediation Services
    61	Educational Services
    62	Health Care and Social Assistance
    71	Arts, Entertainment, and Recreation
    72	Accommodation and Food Services
    81	Other Services (except Public Administration)
    92	Public Administration
    missing => :other,
    :No => :other
    =#

const ntee2naics2names = Dict(
    :AnimalWelfare => :OtherServices,
    :Arts => :Entertainment,
    :CommunityImprovement => :OtherServices,
    :CrimeAndLegal => :ProfServices,
    :Education => :Education,
    :Employment => :OtherServices,
    :Environmental => :OtherServices,
    :FoodAndAgriculture => :Agriculture,
    :ForeignAffairs => :OtherServices,
    :GrantmakingFoundations => :OtherServices,
    :Health => :HealthCare,
    :Housing => :OtherServices,
    :HumanServices => :HealthCare,
    :MedicalDiseases => :OtherServices,
    :MedicalResearch => :ProfServices,
    :MentalHealth => :HealthCare,
    :MutualBenefitOrg => :OtherServices,
    :PublicSafety => :HealthCare,
    :Recreation => :Entertainment,
    :ReligionRelated => :OtherServices,
    :ScienceandTechResearch => :ProfServices,
    :SocialAction => :OtherServices,
    :SocialScienceResearch => :ProfServices,
    :SocietyBenefitAndMultipurpose => :Utilities,
    :YouthDevelopment => :OtherServices,
    :Unknown=>:Unknown,
    missing=>:Unknown,
    :No => :Unknown,
    :AllOthers=>:Unknown
)

const LBM_XFIELDS = Dict(
  :bmsp500=>[:lsp500t30],
  :bmff3=>[:lmkt3t30, :lSMB3, :lHML3],
  :bmff5=>[:lmkt5t30, :lSMB5, :lHML5, :lRMW5, :lCMA5])

const FF_SUFFIX_FIELDS = Dict(
"3"=>[:Mkt_RF, :SMB, :HML, :RF],
"5"=>[:Mkt_RF, :SMB, :HML, :RMW, :CMA, :RF])

#misc dictionary for shortening field names
const SHORT_NAME_INDEX = Dict(
  :qlreturn3yr=>:ret1y,
  :qlreturn3yr=>:ret2y,
  :qlreturn3yr=>:ret3y,
  :qlreturn4yr=>:ret4y,
  :qlreturn5yr=>:ret5y,
  :qpprogramexpense=>:Exp,
  :qpcontributions=>:Contrib,
  :qDpcontributions=>:DifContrib
)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using BinaryBuilder

version = v"2.9.9"

sources = [
    "https://github.com/GNOME/libxml2/archive/v$(version).tar.gz" =>
        "d673f0284cec867ee00872a8152e0c3c09852f17fd9aa93f07579a37534f0bfe",
]

script = raw"""
cd ${WORKSPACE}/srcdir/libxml2-*
./autogen.sh
./configure --prefix=${prefix} --host=${target} --without-python --with-zlib=${prefix}
make -j${nproc} install
"""

products(prefix) = [
    LibraryProduct(prefix, "libxml2", :libxml2)
]

platforms = supported_platforms()

dependencies = [
    "https://github.com/bicycle1885/ZlibBuilder/releases/download/v1.0.3/build_Zlib.v1.2.11.jl",
]

build_tarballs(ARGS, "XML2Builder", version, sources, script, platforms, products, dependencies)


text:
@testset "dump" begin

    n = 1000
    chain = make_chain(n)

    for i in eachindex(chain)
        b1 = chain[i] |> Block |> BTCParser.dump_block_data
        b2 = chain[i] |> BTCParser.dump_block_data
        @test b1 == b2
    end

end

@testset "dump segwit" begin

    # Points to the very first segwit transaction, probably only works on my
    # personal copy of the blockchain:
    #
    # The block is at height 481825 (one based, if you ask bitcoin-cli, height is
    # 481824)
    link = Link(
        BTCParser.to_unsigned((0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
                               0x00, 0x1c, 0x80, 0x18, 0xd9, 0xcb, 0x3b, 0x74,
                               0x2e, 0xf2, 0x51, 0x14, 0xf2, 0x75, 0x63, 0xe3,
                               0xfc, 0x4a, 0x19, 0x02, 0x16, 0x7f, 0x98, 0x93)[end:-1:1]),
        UInt64(976),
        UInt64(0x577479)
    )

    b1 = Block(link) |> BTCParser.dump_block_data
    b2 = link |> BTCParser.dump_block_data

    @test b1 == b2
end


text:
## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-information-about-the-users-current-playback

"""
    player_get_state(;additional_types::String="track", market::String="US")

**Summary**: Get information about the user’s current playback state, including track or episode,
             progress, and active device.

# Optional keywords
- `additional_types::String` : "track" (default) or "episode"
- `market::String` : An ISO 3166-1 alpha-2 country code. If a country code is specified,
                     only episodes that are available in that market will be returned.
                     Default is set to "US".

# Example
```julia-repl
julia> Spotify.player_get_state()[1]
[ Info: We try the request without checking if current grant includes scope user-read-playback-state.
JSON3.Object{Base.CodeUnits{UInt8, String}, Vector{UInt64}} with 10 entries:
  :device                 => {…
  :shuffle_state          => false
  :repeat_state           => "off"
  :timestamp              => 1636493367689
  :context                => {…
  :progress_ms            => 66454
  :item                   => {…
  :currently_playing_type => "track"
  :actions                => {…
  :is_playing             => true
```
"""
function player_get_state(;additional_types::String="track", market::String="US")

    url = "me/player?additional_types=$additional_types&market=$market"

    return Spotify.spotify_request(url; scope = "user-read-playback-state")

end

## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-a-users-available-devices

"""
    player_get_devices()

**Summary**: Get information about a user’s available devices.

# Example
```julia-repl
julia> Spotify.player_get_devices()[1]["devices"]
2-element JSON3.Array{JSON3.Object, Base.CodeUnits{UInt8, String}, SubArray{UInt64, 1, Vector{UInt64}, Tuple{UnitRange{Int64}}, true}}:
{
                   "id": "your_device_id",
            "is_active": false,
   "is_private_session": false,
        "is_restricted": false,
                 "name": "Web Player (Chrome)",
                 "type": "Computer",
       "volume_percent": 100
}
```
"""
function player_get_devices()

    return Spotify.spotify_request("me/player/devices"; scope = "user-read-playback-state")

end

## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-the-users-currently-playing-track

"""
    player_get_current_track(;additional_types::String="track", market::String="US")

**Summary**: Get the object currently being played on the user's Spotify account.

# Optional keywords
- `additional_types::String` : "track" (default) or "episode"
- `market::String` : An ISO 3166-1 alpha-2 country code. If a country code is specified,
                     only episodes that are available in that market will be returned.
                     Default is set to "US".

# Example
```julia-repl
julia> Spotify.player_get_current_track()[1]
[ Info: We try the request without checking if current grant includes scope user-read-playback-state.
JSON3.Object{Base.CodeUnits{UInt8, String}, Vector{UInt64}} with 7 entries:
  :timestamp              => 1636491068506
  :context                => {…
  :progress_ms            => 5265
  :item                   => {…
  :currently_playing_type => "track"
  :actions                => {…
  :is_playing             => true
```
"""
function player_get_current_track(;additional_types::String="track", market::String="US")

    url = "me/player/currently-playing?additional_types=$additional_types&market=$market"

    return Spotify.spotify_request(url; scope = "user-read-playback-state")

end

## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-recently-played

"""
    player_get_recent_tracks(;duration::Int64=1, limit::Int64=20)

**Summary**: Get current user's recently played tracks.

# Optional keywords
- `duration::Int64` : Number of days to look in the past, default is set to 1
- `limit::Int64` : Maximum number of items to return, default is set to 20

# Example
```julia-repl
julia> Spotify.player_get_recent_tracks()[1]
[ Info: We try the request without checking if current grant includes scope user-read-recently-played.
JSON3.Object{Base.CodeUnits{UInt8, String}, Vector{UInt64}} with 5 entries:
  :items   => JSON3.Object[{…
  :next    => "https://api.spotify.com/v1/me/player/recently-played?after=1636123644988&limit=20"
  :cursors => {…
  :limit   => 20
  :href    => "https://api.spotify.com/v1/me/player/recently-played?after=1636410050&limit=20"
```
"""
function player_get_recent_tracks(;duration::Int64=1, limit::Int64=20)

    # Subtract duration from current date and convert to Int64
    starting_date = Dates.datetime2unix(Dates.now() - Dates.Day(duration))
    after = round(Int64, starting_date)

    url = "me/player/recently-played?after=$after&limit=$limit"

    return Spotify.spotify_request(url; scope = "user-read-recently-played")

end


text:
using Documenter, RayTracer

makedocs(
    modules = [RayTracer],
    doctest = true,
    sitename = "RayTracer",
    authors = "Avik Pal",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        assets     = ["assets/raytracer.css"],
    ),
    pages = [
        "Home"              => "index.md",
        "Getting Started"   => [
            "Introduction to Rendering" => "getting_started/teapot_rendering.md",
            "Inverse Lighting"          => "getting_started/inverse_lighting.md",
            "Optimizing using Optim.jl" => "getting_started/optim_compatibility.md",
        ],
        "API Documentation" => [
            "General Utilities"         => "api/utilities.md",
            "Differentiation"           => "api/differentiation.md",
            "Scene Configuration"       => "api/scene.md",
            "Renderers"                 => "api/renderers.md",
            "Optimization"              => "api/optimization.md",
            "Acceleration Structures"   => "api/accelerators.md",
        ],
    ],
)

deploydocs(
    repo = "github.com/avik-pal/RayTracer.jl.git"
)


text:
using LAP_julia
using ImageFiltering: Fill, KernelFactors.gaussian, centered, kernelfactors, imfilter!, padarray, Pad
using LinearAlgebra: qr
using PyPlot: Figure
using ComputationalResources, Images
import SpecialFunctions: erfcinv

# function filt_onebyone!(imgfilt, img, kernel, filter_half_size, points)
#     padded = padarray(img, Pad(:symmetric, filter_half_size, filter_half_size))
#
#     @simd for k in 1:size(points, 2)
#         ind_filt = (points[1, k]:points[1, k], points[2, k]:points[2, k])
#         imfilter!(CPU1(Algorithm.FIRTiled()), imgfilt, padded, kernel, NoPad(), ind_filt)
#     end
# end

"""
    highpass_image(img, window_half_size)

High-pass `img` with a hight-pass filter. This high-pass filter is `1-gaussian`.
"""
function highpass_image(img, window_half_size)

    σ = 2 * window_half_size[1]

    gaus = gaussian(σ, 2 * σ + 1)

    filt_img = imfilter(img, kernelfactors((gaus, gaus)), "symmetric")
    return img .- filt_img
end

"""
    mapped_out(flow)

Return a mask that has `true` at the locations where `flow` would take the pixel out bounds of the image size.
"""
function mapped_out(flow::Flow)::BitArray{2}

    flow_size = size(flow)

    X = ones(flow_size[1]) * collect(1:flow_size[2])'
    Y = collect(1:flow_size[1]) * ones(flow_size[2])'

    return ((X .- real.(flow)) .> flow_size[2]) .|
           ((X .- real.(flow)) .< 1) .|
           ((Y .- imag.(flow)) .> flow_size[1]) .|
           ((Y .- imag.(flow)) .< 1)
end

function mapped_out(inds, img)
    a = indices_spatial(img)
    rng = extrema.(a)
    is_in(ind) = is_in_bounds(ind[1], rng[1]...) && is_in_bounds(ind[2], rng[2]...)::Bool
    return filter(x -> is_in(x), inds)
end

"""
    estimation_noise_variance(img::Image)

Estimate the noise variance of an image y, by evaluating the Median of its Absolute Difference.

Note: Implementation in matlab by Thierry Blu, the Chinese University of Hong kong, Shatin, Hong Kong
"""
function estimation_noise_variance(img::Image)
    kern = centered([0 1 0;1 -4 1;0 1 0]/sqrt(20))
    dy = abs.(imfilter(img.*255, kern, Fill(0,kern)))
    return (median(dy[:])/(erfcinv(0.5)*sqrt(2)))^2
end

"""
    filter_border_inds(inds, img_size, border_width)

Filter out the CartesianIndices `inds` that are within the `border_width` of an image of size `img_size`.
"""
function filter_border_inds(inds, img_size, border_width)
    ind_mask = trues(length(inds))

    @simd for k in 1:length(inds)
        if (inds[k][1] < border_width) ||
            (inds[k][2] < border_width) ||
            (inds[k][1] >= (img_size[1]-border_width)) ||
            (inds[k][2] >= (img_size[2]-border_width))
           ind_mask[k] = false
       end
    end
    return inds[ind_mask]
end

function filt_onebyone!(imgfilt, img, kernel, fhs, points)
    padded = padarray(img, Pad(:symmetric, fhs, fhs))

    @simd for k in 1:size(points, 2)
        ind_filt = (points[1, k]-fhs:points[1, k]+fhs, points[2, k]-fhs:points[2, k]+fhs)
        imfilter!(CPU1(Algorithm.FIRTiled()), imgfilt, padded, kernel, NoPad(), ind_filt)
    end
end

# NOTE: can parallelize
"""
    multi_mat_div_qr(A, b)

Return `E`, where `E[i, :]` is the solutution of the least squares problem ``\\min\\|D_ix - c_i\\|^2`` for each ``D_i``
and ``c_i``, where ``D_i`` is `A[:, :, i]`, ``c_i`` is `b[:, i]` and ``i`` is the size of the second dimension
of `b` and third dimension of `A`. In other words, each row of `E` is the solution to one matrix from `A` and it's
corresponding vector from `b`.
"""
function multi_mat_div_qr(A, b)
    res = Array{Float64,2}(undef, size(A)[1], size(A)[3])
    # Threads.@threads for j in axes(A, 4) check: (https://stackoverflow.com/questions/57678890/batch-matrix-multiplication-in-julia)
    @simd for k in axes(A)[3]
        # try
            @views res[:, k] = qr(A[:, :, k], Val(true)) \ b[:, k]
        # catch
            # @info A[:, :, k], b[:, k]
        # end
        # @views res[:, k] = A[:, :, k] \ b[:, k]
    end
    return transpose(res)
end

# TODO change to inds
function multi_mat_div_at_points(A, b, inds, image_size)
    points = inds_to_points(inds)
    # [2, pix_count]
    res = similar(A, size(A, 1) + 1, size(points, 2)) # example size = (filter_count, point_count)
    A = reshape(A, size(A)[1], size(A)[1], image_size...)
    b = reshape(b, size(b)[1], image_size...)

    for k in axes(points)[2]
        # ind = points[:, k]...
        @views res[:, points[:, k]...] .= cat(qr(A[:, :, points[:, k]...], Val(true)) \ b[:, points[:, k]...], dims=2)
    end
    res = reshape(res, size(res)[1], :)
    return transpose(res)
end

# NOTE: here maybe an average could be better so the size of the window
# doest effect the number range of the results
# also might improve the speed if the pixels are saved into a tmp var
"""
    function window_sum!(result, pixels, img_size, window)

Get a sum of values (using a "symmetric" padding on the borders) in a window `window` around each point and saving the
sum into `result`.

Note: Uses a filter algorithm with a 2D ones kernel.
Note: For small images (around `256x256`) and small filters (window of size 9 and below) is faster than [`window_sum3!`](@ref). <- factcheck this again
"""
function window_sum!(result, pixels, img_size, window)
    # prepare a kernel of ones of the window size
    ones_arr = centered(ones(window[1]))
    ones_kernel = kernelfactors((ones_arr, ones_arr))

    # filtering gets a sum of pixels of window size in each coord
    imfilter!(reshape(result, img_size), reshape(pixels, img_size), ones_kernel, "symmetric")
    return result
end

"""
    function window_sum3!(result, pixels, img_size, window)

Get a sum of values (using a "symmetric" padding on the borders) in a window `window` around each point and saving the
sum into `result`.

Note: Uses a cumsum algorithm.
Note: For small images (around `256x256`) and large filters (window of size 11 and above) is faster than [`window_sum!`](@ref). <- factcheck this again
"""
function window_sum3!(result, pixels, img_size, window)
    w = Int64.((window[1]-1)/2)
    summed = padarray(reshape(pixels, img_size), Pad(:symmetric, (w+1, w+1), (w, w)))
    @views summed[-w, :] .= 0
    @views summed[:, -w] .= 0

    summed = cumsum2d!(summed, summed)
    (c, d) = img_size

    reshape(result, img_size) .= view(summed, (1+w):(w+c), (1+w):(w+d)) .-
                                 view(summed, (1+w):(w+c), (-w):(d-w-1)) .-
                                 view(summed, (-w):(c-w-1), (1+w):(w+d)) .+
                                 view(summed, (-w):(c-w-1), (-w):(d-w-1))
    return result
end

"""
    function window_sum3!(result, pixels, img_size, window, inds)

Get a sum of values (using a "symmetric" padding on the borders) in a window `window` around `inds` and saving the
sum into `result`.

Note: Uses a cumsum algorithm.
"""
function window_sum3_at_inds!(result, pixels, img_size, window, inds)

    w = Int64.((window[1]-1)/2)
    summed = padarray(reshape(pixels, img_size), Pad(:symmetric, (w+1, w+1), (w, w)))
    @views summed[-w, :] .= 0
    @views summed[:, -w] .= 0
    summed = cumsum2d!(summed, summed)

    a = CartesianIndex(w, w)
    b = CartesianIndex(w, -w-1)
    c = CartesianIndex(-w-1, w)
    d = CartesianIndex(-w-1, -w-1)
    for (matrix_ind, result_ind)  in zip(inds, CartesianIndices(result))
        result[result_ind] = summed[matrix_ind+a] -
                          summed[matrix_ind+b] -
                          summed[matrix_ind+c] +
                          summed[matrix_ind+d]
    end
    return result
end

"""
    to_lin_index(cart_ind, size)

Get linear index from CartesianIndex `cart_ind` of 2D matrix of size `size`.
"""
@inline function to_lin_index(cart_ind, size)
    return size[2]*(cart_ind[2]-1)+cart_ind[1]
end

"""
    cumsum2d!(result, A)

Cummulative sum over 2D matrix `A`, storing the result into `result`.
"""
function cumsum2d!(result, A)
    cumsum!(result, A, dims = 1)
    cumsum!(result, result, dims = 2)
    return result
end

"""
    prepare_gaussian_filters(filter_half_size)

Get KernelFactors for the first 3 forward and backward gaussian filters from the LAP paper.
"""
function prepare_gaussian_filters(filter_half_size)

    filter_size = 2 * filter_half_size + 1

    # Calculate separable filters from basis:
    sigma = (filter_half_size + 2) / 4
    centered_inds = centered(-filter_half_size:filter_half_size)
    gaus = gaussian(sigma, filter_size)
    gaus_der = gaus .* centered_inds .* (-1)/sigma^2
    gaus_der_flip = reverse(gaus_der, dims=1)

    forward_kernels = [kernelfactors((gaus, gaus)),
                       kernelfactors((gaus, gaus_der_flip)),
                       kernelfactors((gaus_der_flip, gaus))]

    backward_kernels = [kernelfactors((gaus, gaus)),
                        kernelfactors((gaus, gaus_der)),
                        kernelfactors((gaus_der, gaus))]

    return forward_kernels, backward_kernels
end

function add_at_points(A, B, inds)
    for ind in inds
        if isnan(B[ind])
            continue
        end
        @assert !isnan(A[ind])
        println(B[ind], ind)
        A[ind] = A[ind] + B[ind]
    end
    return A
end

function sum_at_points(A, inds)
    return sum([A[ind] for ind in inds])
end

"""
    gem3d!(A, b)

Reduce the linear systems of equations ``C_n x_n = d_n`` to the row echelon form. The matrixes ``C``
are slices from the first 2 dimensions of `A`, the vectors ``d`` are the slices from the first dimension
of `b` and ``n`` is the size of the third dimension of `A`.

See also: [`gem3d!`](@ref), [`back_substitution3d`](@ref)
"""
function gem3d!(A, b)
    ratios = similar(A, size(A, 3))
    for k in axes(A,1)
        for l in k+1:size(A, 2)
            ratios .= A[l, k, :] ./ A[k, k, :]
            A[l, :, :] .-= ratios' .* A[k, :, :]
            b[l, :] .-= ratios .* b[k, :]
        end
    end
end

"""
    back_substitution3d(A, b)

Solve the linear systems of equations ``C_n x_n = d_n``, where the matrixes ``C`` (in row echelon form)
are slices from the first 2 dimensions of `A`, the vectors ``d`` are the slices from the first dimension
of `b` and ``n`` is the size of the third dimension of `A`.
Returns the coefficients ``x_n`` as a 2D matrix `coeffs`, where the vectors ``x_n`` are the rows of this matrix.

See also: [`gem3d!`](@ref), [`multi_mat_div_gem`](@ref)
"""
function back_substitution3d(A, b)
    row_count = size(A, 1)
    coeffs = similar(A, (size(A, 3), size(A, 1)))
    for k in row_count:-1:1
        coeffs[:,k] = (b[k,:])'
        for m in (k+1):row_count
            coeffs[:,k] = coeffs[:,k]-A[k,m,:] .* coeffs[:,m]
        end
        coeffs[:,k]=coeffs[:,k]./A[k,k,:];
    end
    return coeffs
end

"""
    multi_mat_div_gem(A, b)

Solve the linear systems of equations ``C_n x_n = d_n``, where the matrixes ``C``
are slices from the first 2 dimensions of `A`, the vectors ``d`` are the slices from the first dimension
of `b` and ``n`` is the size of the third dimension of `A`.
Returns the coefficients ``x_n`` as a 2D matrix `coeffs`, where the vectors ``x_n`` are the rows of this matrix.

Note: Uses Gaussian elimination and back substitution.
Note: For small matrices ``C``, this is an order of magnitude faster than

See also: [`gem3d!`](@ref), [`back_substitution3d`](@ref)
"""
function multi_mat_div_gem(A, b)
    gem3d!(A, b)
    return back_substitution3d(A, b)
end


text:
using PastaQ
using Random
using ITensors

Random.seed!(1234)
N = 10  # Number of spins
B = 1.0 # Transverse magnetic field

# Build Ising Hamiltonian
sites = siteinds("Qubit", N)
ampo = AutoMPO()
for j in 1:(N - 1)
  # Ising ZZ interactions
  ampo .+= -1, "Z", j, "Z", j + 1
end
for j in 1:N
  # Transverse field X
  ampo .+= -B, "X", j
end
H = MPO(ampo, sites)

# Run DMRG
dmrg_iter = 5      # DMRG steps
dmrg_cutoff = 1E-10   # Cutoff
Ψ0 = randomMPS(sites) # Initial state
sweeps = Sweeps(dmrg_iter)
maxdim!(sweeps, 10, 20, 30, 40, 50, 100)
cutoff!(sweeps, dmrg_cutoff)
# Run DMRG
println("Running DMRG to get ground state of transverse field Ising model:")
E, Ψ = dmrg(H, Ψ0, sweeps)
@show maxlinkdim(Ψ)
println()

# Generate data
nshots = 10_000
data = getsamples(Ψ, nshots; local_basis=["X", "Y", "Z"])

# Quantum state tomography
# Initialize variational state
χ = maxlinkdim(Ψ)
ψ0 = randomstate(Ψ0; χ=χ)

# Measurements
Energy(ψ::MPS) = inner(ψ, H, ψ)
F(ψ::MPS) = fidelity(ψ, Ψ)

# Initialize observer
obs = Observer([F, Energy, ("Z", 1, "Z", N ÷ 2)])

ZZ = PastaQ.measure(Ψ, ("Z", 1, "Z", N ÷ 2))
@printf("⟨Ψ|Ĥ|Ψ⟩ =  %.5f   ", E)
@printf("⟨Ψ|Ŝᶻ(1)Ŝˣ(%d)|Ψ⟩ = %.5f   (DMRG)\n", N ÷ 2, ZZ)
# Run tomography
println("Running tomography to learn the Ising model ground state from sample data")
ψ = tomography(
  data,
  ψ0;
  optimizer=SGD(; η=0.01),
  batchsize=500,
  epochs=10,
  (observer!)=obs,
  print_metrics=["Energy", "F", "Z(1)Z(5)"],
)

@show maxlinkdim(ψ)
println()


text:
module YaoEGraph

using Metatheory
using Metatheory.EGraphs
@metatheory_init

# Custom type APIs for the GATExpr
using YaoHIR
using Metatheory.TermInterface

TermInterface.gethead(t::Chain) = :call
TermInterface.getargs(t::Chain) = [:chain, t.args...]
TermInterface.gethead(t::Gate) = :call
TermInterface.getargs(t::Gate) = [:gate, t.locations]
TermInterface.gethead(t::Ctrl) = :call
TermInterface.getargs(t::Ctrl) = [:ctrl, t.gate, t.ctrl]

function EGraphs.extractnode(::EGraph, n::ENode{Chain}, extractor::Function)
    return Chain(extractor.(n.args[2:end])...)
end

end


text:
##### Generic methods #####

# sampling

rand!{M<:Matrix}(d::MatrixDistribution, A::AbstractArray{M}) = _rand!(sampler(d), A)
rand(d::MatrixDistribution, n::Int) = _rand!(sampler(d), Array(Matrix{eltype(d)}, n))

# pdf & logpdf

_pdf{T<:Real}(d::MatrixDistribution, x::AbstractMatrix{T}) = exp(_logpdf(d, x))

function logpdf{T<:Real}(d::MatrixDistribution, x::AbstractMatrix{T})
    size(x) == size(d) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    _logpdf(d, x)
end

function pdf{T<:Real}(d::MatrixDistribution, x::AbstractMatrix{T})
    size(x) == size(d) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    _pdf(d, x)
end

function _logpdf!{M<:Matrix}(r::AbstractArray, d::MatrixDistribution, X::AbstractArray{M})
    for i = 1:length(X)
        r[i] = logpdf(r, X[i])
    end
    return r
end

function _pdf!{M<:Matrix}(r::AbstractArray, d::MatrixDistribution, X::AbstractArray{M})
    for i = 1:length(X)
        r[i] = pdf(r, X[i])
    end
    return r
end

function logpdf!{M<:Matrix}(r::AbstractArray, d::MatrixDistribution, X::AbstractArray{M})
    length(X) == length(r) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    _logpdf!(r, d, X)
end

function pdf!{M<:Matrix}(r::AbstractArray, d::MatrixDistribution, X::AbstractArray{M})
    length(X) == length(r) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    _pdf!(r, d, X)
end

function logpdf{M<:Matrix}(d::MatrixDistribution, X::AbstractArray{M})
    T = promote_type(partype(d), eltype(M))
    _logpdf!(Array(T, size(X)), d, X)
end

function pdf{M<:Matrix}(d::MatrixDistribution, X::AbstractArray{M})
    T = promote_type(partype(d), eltype(M))
    _pdf!(Array(T, size(X)), d, X)
end

##### Specific distributions #####

for fname in ["wishart.jl", "inversewishart.jl"]
    include(joinpath("matrix", fname))
end


text:
mutable struct Bar
           baz
           qux::Float64
       end


text:
import Interpolations:
    interpolate,
    BSpline,
    Quadratic,
    Line,
    OnGrid,
    scale,
    gradient1,
    extrapolate,
    Linear,
    Gridded,
    NoInterp,
    Cubic,
    Constant,
    Flat

"""
$(SIGNATURES)

Construct interpolation of N-D array `y` along its last dimension on grid `x`. `method` specifies the interpolation algorithm, which can be either "BSpline" or "Gridded". If `x` is an `AbstractRange`, the default method is "BSpline", otherwise the default `method` is "Gridded". `order` is the interpolation order. If `x` is `AbstractRange`, the default order is 2, otherwise the default order is 1. `extrapolation` specifies the extrapolation methods, which defines what happens if you try to index into the interpolation object with coordinates outside of its bounds in any dimension. The implemented extrapolation methods are "line", "flat", or you can pass a constant to be used as a "fill" value returned for any out-of-bounds evaluation. The default value is 0.
"""
function construct_interpolations(
    x::AbstractRange{S},
    y::AbstractArray{T,N};
    method="BSpline",
    order=2,
    extrapolation=0,
) where {S <: Real,T <: Number,N}
    method = interp_tuple(N, method, order)
    itp = interpolate(y, method)
    itp = extrapolate(itp, extrap_method(extrapolation))
    index = interp_index(x, size(y))
    scale(itp, index...)
end

function construct_interpolations(
    x::AbstractArray{S},
    y::AbstractArray{T,N};
    method="Gridded",
    order=1,
    extrapolation=0,
) where {S <: Real,T <: Number,N}
    if lowercase(method) == "bspline"
        Δ = diff(x)
        if allequal(Δ)
            x = range(x[1], x[end], length=length(x))
            return construct_interpolations(
                x,
                y,
                method=method,
                order=order,
                extrapolation=extrapolation,
            )
        else
            @warn "The grid is not uniform. Using grided linear interpolation."
            method = "Gridded"
        end
    end
    method = interp_tuple(N, method, order)
    index = interp_index(x, size(y))
    itp = interpolate(index, y, method)
    extrapolate(itp, extrap_method(extrapolation))
end

function construct_interpolations(
    x::AbstractRange{S},
    y::AbstractArray{W,1};
    method="BSpline",
    order=2,
    extrapolation=0,
) where {S <: Real,W <: AbstractArray}
    y = cat(y...; dims=ndims(y[1]) + 1)
    construct_interpolations(
        x,
        y,
        method=method,
        order=order,
        extrapolation=extrapolation,
    )
end

"""
$(SIGNATURES)

Construct interpolation of a nested array. If `x` is an `AbstractRange`, the nested array will be converted into multi-dimensional array for interpolation. Otherwise only "Gridded" method and `order` 1 is supported.
"""
function construct_interpolations(
    x::AbstractArray{S,1},
    y::AbstractArray{W,1};
    method="Gridded",
    order=1,
    extrapolation=0,
) where {S <: Real,W <: AbstractArray}
    itp = interpolate((x,), y, interp_method(method, order))
    extrapolate(itp, extrap_method(extrapolation))
end

function interp_index(x, s)
    res = []
    for i = 1:(length(s) - 1)
        push!(res, 1:s[i])
    end
    push!(res, x)
    (res...,)
end

function interp_tuple(N, method, order)
    res = []
    for i = 1:(N - 1)
        push!(res, NoInterp())
    end
    interp_term = interp_method(method, order)
    push!(res, interp_term)
    (res...,)
end

function interp_method(method, order::Integer; boundary=Line(OnGrid()))
    if lowercase(method) == "bspline"
        if order == 0
            res = BSpline(Constant())
        elseif order == 1
            res = BSpline(Linear())
        elseif order == 2
            res = BSpline(Quadratic(boundary))
        elseif order == 3
            res = BSpline(Cubic(boundary))
        else
            throw(ArgumentError("Order $order for method $method is not implemented."))
        end
    elseif lowercase(method) == "gridded"
        if order == 0
            res = Gridded(Constant())
        elseif order == 1
            res = Gridded(Linear())
        else
            throw(ArgumentError("Oder $order for method $method is not implemented."))
        end
    elseif lowercase(method) == "none"
        res = NoInterp()
    else
        throw(ArgumentError("Method $method is invalid."))
    end
    res
end

function extrap_method(method::String)
    if lowercase(method) == "line"
        Line()
    elseif lowercase(method) == "flat"
        Flat()
    else
        throw(ArgumentError("Extrapolation method: $method is not supported."))
    end
end

extrap_method(method::Number) = method

"""
$(SIGNATURES)

Calculate the gradient of `itp` at `s`.
"""
function gradient(itp, s::Number)
    gradient1(itp, s)
end

function gradient(itp, s::AbstractArray)
    [gradient1(itp, x) for x in s]
end


text:
#using LinearAlgebraicRepresentation
#using SparseArrays, DataStructures, NearestNeighbors
#Lar = LinearAlgebraicRepresentation
#using ViewerGL
#GL = ViewerGL

function chaincongruence(W, Delta_0, Delta_1; epsilon=1e-6)
	V, vclasses = vcongruence(W, epsilon)
	EV = cellcongruence(Delta_0, vclasses, dim=1)
	EV = sort(union(map(sort, EV)))
	FV = cellcongruence(Delta_1 * Delta_0, vclasses, dim=2)
	FV = sort(union(map(sort, FV)))
	copFE = coboundary_1(FV::Lar.Cells, EV::Lar.Cells)
	FE = cop2lar(copFE)
println("\n>>>>> final output from chaincongruence\n")
@show V;
@show EV;
@show FV;
@show FE;
	return V, EV, FV, FE
end

function vcongruence(V::Matrix,epsilon)
    vclasses, visited = [], []
    kdtree = NearestNeighbors.KDTree(V);
    for vidx = 1 : size(V, 2) if !(vidx in visited)
       nearvs = NearestNeighbors.inrange(kdtree, V[:,vidx], epsilon)
       push!(vclasses, nearvs)
			 append!(visited, nearvs) end
    end
    W = hcat([sum(V[:,class], dims=2)/length(class) for class in vclasses]...)
    return W, vclasses
end

function cellcongruence(Delta, vclasses; dim)  # ==> OK !!!
  # cells  with old vertices
	cellarray = Lar.cop2lar(Delta)
	# conversion array old -> new vertex
	newvert = Vector(undef, size(Delta,2))
  for (k, class) in enumerate(vclasses)
  	for v in class
  	  newvert[v] = k
  	end
  end
	# conversione delle celle da vecchi a nuovi indici di vertice
	newcells = Vector{Int64}[]
	for cell in cellarray
		newcell = Int64[]
		for v in cell
			push!(newcell, newvert[v])
		end
		length(newcell) > dim ? push!(newcells, newcell) : break
	end
	# remove duplicate cells
	newcells = union(map(sort,newcells))
		# remove duplicate vertices
		newcells = map(union,newcells)
	# remove empty cells
	outcells = filter(x -> !(length(x) <= dim || all(y->y==x[1], x)), newcells)
end

# test elementare con cubo unitario

#using LinearAlgebraicRepresentation
#using SparseArrays, DataStructures, NearestNeighbors
#Lar = LinearAlgebraicRepresentation
#using ViewerGL
#GL = ViewerGL
#
#W = convert(Matrix,[0.5310492999999998 0.8659989999999999 0.14191280000000003; 1.0146684 0.6827212999999999 0.2169682; 0.3477716 0.5268921 0.4947971000000001; 0.8313907882395298 0.3436144447063971 0.5698524407571428; 0.6061046999999998 1.2188832999999994 0.5200012; 1.0897237999999998 1.0356056999999999 0.5950565999999999; 0.42282699999999984 0.8797763999999998 0.8728855; 0.9064461979373597 0.6964987903021808 0.9479408896095312; 0.5310493 0.8659989999999999 0.14191279999999987; 1.0146684000000001 0.6827213 0.21696819999999994; 0.6061047 1.2188833 0.5200011999999999; 1.0897237772434623 1.035605657895053 0.5950566438156151; 0.3477716 0.5268921 0.4947971; 0.8313908 0.3436145 0.5698525000000001; 0.422827 0.8797764000000001 0.8728855; 0.9064462 0.6964988 0.9479409; 0.5310493 0.8659989999999999 0.14191280000000006; 0.34777160000000007 0.5268920999999999 0.4947971000000001; 0.6061047 1.2188833 0.5200012000000002; 0.4228270000000002 0.8797764 0.8728855000000001; 1.0146684 0.6827213 0.21696819999999994; 0.8313908 0.3436145 0.5698525000000001; 1.0897238 1.0356057 0.5950565999999999; 0.9064461675456482 0.6964988122992563 0.9479408949632379]');
#Delta_0 = SparseArrays.sparse([1, 3, 1, 4, 2, 3, 2, 4, 5, 7, 5, 8, 6, 7, 6, 8, 9, 11, 9, 12, 10, 11, 10, 12, 13, 15, 13, 16, 14, 15, 14, 16, 17, 19, 17, 20, 18, 19, 18, 20, 21, 23, 21, 24, 22, 23, 22, 24], [1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15, 16, 16, 17, 17, 18, 18, 19, 19, 20, 20, 21, 21, 22, 22, 23, 23, 24, 24], Int8[-1, -1, 1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1]);
#Delta_1 = SparseArrays.sparse([1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24], Int8[1, -1, -1, 1, 1, -1, -1, 1, 1, -1, -1, 1, 1, -1, -1, 1, 1, -1, -1, 1, 1, -1, -1, 1]);
#
#
#V,EV,FE,FV = chaincongruence(W,Delta_0,Delta_1; epsilon=1e-4);
#EV = Lar.removeDups(map(sort!,EV))
#FV = Lar.removeDups(map(sort!,FV))
#
#
#@show V; # centroids of the 8  classes of W points
#@show EV; # edges-by-vertices
#@show FV; # faces-by-vertices
#@show FE; # faces-by-edges
#VV = [[v] for v=1:size(V,2)]
##
#GL.VIEW(push!(GL.numbering(.5)((V,Lar.Cells[VV,EV,FV]), GL.COLORS[1], 0.5),GL.GLFrame2));


text:
cd("/Users/adrienpajon/Documents/all.nosync/git_repo/pulse_analysis/2020-11-12 xp pulse all seuil/Julia_analysis")

using Plots

using StatsPlots

using JLD

# using Pyplot
pyplot()

include("../functions/analysisPULSE.jl")

####################################################################
gel=load("results/PULSE_analysis_gel.jld")
pleine=load("results/PULSE_analysis_pleine.jld")
pneu=load("results/PULSE_analysis_pneu.jld")

filenameMTiXsens=[gel["filenameMTiXsens"];pleine["filenameMTiXsens"];pneu["filenameMTiXsens"]]
sensorXsensNames=gel["sensorXsensNames"]
Energysignal_MTiXsens=merge(gel["Energysignal_MTiXsens"],pleine["Energysignal_MTiXsens"],pneu["Energysignal_MTiXsens"])
rms_MTiXsens=merge(gel["rms_MTiXsens"],pleine["rms_MTiXsens"],pneu["rms_MTiXsens"])
vdv_MTiXsens=merge(gel["vdv_MTiXsens"],pleine["vdv_MTiXsens"],pneu["vdv_MTiXsens"])
crest_factor_MTiXsens=merge(gel["crest_factor_MTiXsens"],pleine["crest_factor_MTiXsens"],pneu["crest_factor_MTiXsens"])
Energysignal_MTiXsens_polar=merge(gel["Energysignal_MTiXsens_polar"],pleine["Energysignal_MTiXsens_polar"],pneu["Energysignal_MTiXsens_polar"])
rms_MTiXsens_polar=merge(gel["rms_MTiXsens_polar"],pleine["rms_MTiXsens_polar"],pneu["rms_MTiXsens_polar"])
vdv_MTiXsens_polar=merge(gel["vdv_MTiXsens_polar"],pleine["vdv_MTiXsens_polar"],pneu["vdv_MTiXsens_polar"])

splitName,
plot_order,
perturbLabel,
perturbLabel_unique,
perturbIndices,
wheelType,
wheelLabel_unique,
wheelIndices=splitFileName(filenameMTiXsens)

bool_perturb_choice="gel"
bool_perturb=wheelType.==bool_perturb_choice

filenameMTiXsens_=filenameMTiXsens[bool_perturb]

splitName,
plot_order,
perturbLabel,
perturbLabel_unique,
perturbIndices,
wheelType,
wheelLabel_unique,
wheelIndices=splitFileName(filenameMTiXsens_,0.05,0)

include("../functions/analysisPULSE.jl")
plotAnalysisPULSE_MTiXsens(
        Energysignal_MTiXsens,
        rms_MTiXsens,
        vdv_MTiXsens,
        crest_factor_MTiXsens,
        filenameMTiXsens_,
        sensorXsensNames,
        ["x-y";"z"],
        perturbLabel_unique,
        plot_order,
        perturbIndices,
        [:blue :red :green :orange :yellow],
        [:rect :rect :rect :rect :circle],
        string("\n","wheel : ",bool_perturb_choice))

plotAnalysisPULSE_MTiXsens(
        Energysignal_MTiXsens_polar,
        rms_MTiXsens_polar,
        vdv_MTiXsens_polar,
        crest_factor_MTiXsens,
        filenameMTiXsens_,
        sensorXsensNames,
        ["r";"\u03B8"],
        perturbLabel_unique,
        plot_order,
        perturbIndices,
        [:blue :red :green :orange :yellow],
        [:rect :rect :rect :rect :circle],
        string("\n","wheel : ",bool_perturb_choice))

println("PULSE Analysis Finished")

# using JLD
# save("results/PULSE_analysis_gel.jld",
#         "filenameMTiXsens",filenameMTiXsens,
#         "sensorXsensNames",sensorXsensNames,
#         "step_max_perturb_MTiXsens",step_max_perturb_MTiXsens,
#         "Energysignal_MTiXsens",Energysignal_MTiXsens,
#         "rms_MTiXsens",rms_MTiXsens,
#         "vdv_MTiXsens",vdv_MTiXsens,
#         "crest_factor_MTiXsens",crest_factor_MTiXsens,
#         "Energysignal_MTiXsens_polar",Energysignal_MTiXsens_polar,
#         "rms_MTiXsens_polar",rms_MTiXsens_polar,
#         "vdv_MTiXsens_polar",vdv_MTiXsens_polar)


text:
include("dependencies_for_runtests.jl")

using Oceananigans.BoundaryConditions: ImpenetrableBoundaryCondition

""" Take one time step with three forcing functions on u, v, w. """
function time_step_with_forcing_functions(arch)
    @inline Fu(x, y, z, t) = exp(π * z)
    @inline Fv(x, y, z, t) = cos(42 * x)
    @inline Fw(x, y, z, t) = 1.0

    grid = RectilinearGrid(arch, size=(1, 1, 1), extent=(1, 1, 1))
    model = NonhydrostaticModel(grid=grid, forcing=(u=Fu, v=Fv, w=Fw))
    time_step!(model, 1, euler=true)

    return true
end

@inline Fu_discrete_func(i, j, k, grid, clock, model_fields) = @inbounds -model_fields.u[i, j, k]
@inline Fv_discrete_func(i, j, k, grid, clock, model_fields, params) = @inbounds - model_fields.v[i, j, k] / params.τ
@inline Fw_discrete_func(i, j, k, grid, clock, model_fields, params) = @inbounds - model_fields.w[i, j, k]^2 / params.τ

""" Take one time step with a DiscreteForcing function. """
function time_step_with_discrete_forcing(arch)

    Fu = Forcing(Fu_discrete_func, discrete_form=true)

    grid = RectilinearGrid(arch, size=(1, 1, 1), extent=(1, 1, 1))
    model = NonhydrostaticModel(grid=grid, forcing=(u=Fu,))
    time_step!(model, 1, euler=true)

    return true
end

""" Take one time step with ParameterizedForcing forcing functions. """
function time_step_with_parameterized_discrete_forcing(arch)

    Fv = Forcing(Fv_discrete_func, parameters=(τ=60,), discrete_form=true)
    Fw = Forcing(Fw_discrete_func, parameters=(τ=60,), discrete_form=true)

    grid = RectilinearGrid(arch, size=(1, 1, 1), extent=(1, 1, 1))
    model = NonhydrostaticModel(grid=grid, forcing=(v=Fv, w=Fw))
    time_step!(model, 1, euler=true)

    return true
end

""" Take one time step with a Forcing forcing function with parameters. """
function time_step_with_parameterized_continuous_forcing(arch)

    u_forcing = Forcing((x, y, z, t, ω) -> sin(ω * x), parameters=π)

    grid = RectilinearGrid(arch, size=(1, 1, 1), extent=(1, 1, 1))
    model = NonhydrostaticModel(grid=grid, forcing=(u=u_forcing,))
    time_step!(model, 1, euler=true)

    return true
end

""" Take one time step with a Forcing forcing function with parameters. """
function time_step_with_single_field_dependent_forcing(arch, fld)

    forcing = NamedTuple{(fld,)}((Forcing((x, y, z, t, u) -> -u, field_dependencies=:u),))

    grid = RectilinearGrid(arch, size=(1, 1, 1), extent=(1, 1, 1))
    model = NonhydrostaticModel(grid=grid, forcing=forcing,
                                buoyancy=SeawaterBuoyancy(), tracers=(:T, :S))
    time_step!(model, 1, euler=true)

    return true
end

""" Take one time step with a Forcing forcing function with parameters. """
function time_step_with_multiple_field_dependent_forcing(arch)

    u_forcing = Forcing((x, y, z, t, v, w, T) -> sin(v) * exp(w) * T, field_dependencies=(:v, :w, :T))

    grid = RectilinearGrid(arch, size=(1, 1, 1), extent=(1, 1, 1))
    model = NonhydrostaticModel(grid=grid, forcing=(u=u_forcing,),
                                buoyancy=SeawaterBuoyancy(), tracers=(:T, :S))
    time_step!(model, 1, euler=true)

    return true
end

""" Take one time step with a Forcing forcing function with parameters. """
function time_step_with_parameterized_field_dependent_forcing(arch)

    u_forcing = Forcing((x, y, z, t, u, p) -> sin(p.ω * x) * u, parameters=(ω=π,), field_dependencies=:u)

    grid = RectilinearGrid(arch, size=(1, 1, 1), extent=(1, 1, 1))
    model = NonhydrostaticModel(grid=grid, forcing=(u=u_forcing,))
    time_step!(model, 1, euler=true)

    return true
end

function relaxed_time_stepping(arch)
    x_relax = Relaxation(rate = 1/60,   mask = GaussianMask{:x}(center=0.5, width=0.1),
                                      target = LinearTarget{:x}(intercept=π, gradient=ℯ))

    y_relax = Relaxation(rate = 1/60,   mask = GaussianMask{:y}(center=0.5, width=0.1),
                                      target = LinearTarget{:y}(intercept=π, gradient=ℯ))

    z_relax = Relaxation(rate = 1/60,   mask = GaussianMask{:z}(center=0.5, width=0.1),
                                      target = π)

    grid = RectilinearGrid(arch, size=(1, 1, 1), extent=(1, 1, 1))
    model = NonhydrostaticModel(grid=grid, forcing=(u=x_relax, v=y_relax, w=z_relax))
    time_step!(model, 1, euler=true)

    return true
end

function advective_and_multiple_forcing(arch)
    grid = RectilinearGrid(arch, size=(2, 2, 3), extent=(1, 1, 1), halo=(4, 4, 4))

    constant_slip = AdvectiveForcing(UpwindBiasedFifthOrder(), w=1)

    no_penetration = ImpenetrableBoundaryCondition()
    slip_bcs = FieldBoundaryConditions(grid, (Center, Center, Face), top=no_penetration, bottom=no_penetration)
    slip_velocity = ZFaceField(grid, boundary_conditions=slip_bcs)
    velocity_field_slip = AdvectiveForcing(CenteredSecondOrder(), w=slip_velocity)
    simple_forcing(x, y, z, t) = 1

    model = NonhydrostaticModel(; grid, tracers=(:a, :b), forcing=(a=constant_slip, b=(simple_forcing, velocity_field_slip)))
    time_step!(model, 1, euler=true)

    return true
end

@testset "Forcings" begin
    @info "Testing forcings..."

    for arch in archs
        A = typeof(arch)
        @testset "Forcing function time stepping [$A]" begin
            @info "  Testing forcing function time stepping [$A]..."

            @testset "Non-parameterized forcing functions [$A]" begin
                @info "      Testing non-parameterized forcing functions [$A]..."
                @test time_step_with_forcing_functions(arch)
                @test time_step_with_discrete_forcing(arch)
            end

            @testset "Parameterized forcing functions [$A]" begin
                @info "      Testing parameterized forcing functions [$A]..."
                @test time_step_with_parameterized_continuous_forcing(arch)
                @test time_step_with_parameterized_discrete_forcing(arch)
            end

            @testset "Field-dependent forcing functions [$A]" begin
                @info "      Testing field-dependent forcing functions [$A]..."

                for fld in (:u, :v, :w, :T)
                    @test time_step_with_single_field_dependent_forcing(arch, fld)
                end

                @test time_step_with_multiple_field_dependent_forcing(arch)
                @test time_step_with_parameterized_field_dependent_forcing(arch)
            end

            @testset "Relaxation forcing functions [$A]" begin
                @info "      Testing relaxation forcing functions [$A]..."
                @test relaxed_time_stepping(arch)
            end

            @testset "Advective and multiple forcing [$A]" begin
                @info "      Testing advective and multiple forcing [$A]..."
                @test advective_and_multiple_forcing(arch)
            end
        end
    end
end


text:
#=
# Tests vector unknown capability with a simple Poisson-like problem.
=#

### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("vector");

# Try making an optional log
useLog("vectorlog", level=3)

# Set up the configuration
domain(2)
functionSpace(order = 3)

# Specify the problem
mesh(QUADMESH, elsperdim=20, bids=3)

u = variable("u", type=VECTOR)
p = variable("p")

testSymbol("v", type=VECTOR)
testSymbol("w")

# boundary(u, 2, DIRICHLET, [0.1, 0.1])
# boundary(p, 2, DIRICHLET, 0.1)
boundary(u, 1, NEUMANN, ["-pi*sin(2*pi*y)", "-2*pi*sin(3*pi*y)"])
boundary(p, 1, NEUMANN, "-pi*sin(pi*y)")
boundary(u, 2, DIRICHLET, [0, 0])
boundary(p, 2, DIRICHLET, 0)
boundary(u, 3, DIRICHLET, [0, 0])
boundary(p, 3, DIRICHLET, 0)

# Write the weak form
coefficient("f", ["-25*pi*pi*sin(pi*x)*sin(2*pi*y)", "-65*pi*pi*sin(2*pi*x)*sin(3*pi*y)"], type=VECTOR)
coefficient("g", "-2*pi*pi*sin(pi*x)*sin(pi*y)")
coefficient("a", 5)

weakForm([u,p], ["-a*inner(grad(u), grad(v)) - dot(f,v)", "-dot(grad(p), grad(w)) - g*w"])

solve([u,p]);

# # exact solution is [sin(pi*x)*sin(2*pi*y), sin(2*pi*x)*sin(3*pi*y)]
# # check error
erroru = zeros(size(u.values));
maxerru = 0
maxerrp = 0
exactu1(x,y) = sin(pi*x)*sin(2*pi*y);
exactu2(x,y) = sin(2*pi*x)*sin(3*pi*y);
exactp(x,y) = sin(pi*x)*sin(pi*y);

for i=1:size(Finch.grid_data.allnodes,2)
    x = Finch.grid_data.allnodes[1,i];
    y = Finch.grid_data.allnodes[2,i];
    exac = [exactu1(x,y), exactu2(x,y)];
    for j=1:size(Finch.grid_data.allnodes,1)
        erroru[j,i] = u.values[j,i] - exac[j];
        global maxerru;
        maxerru = max(abs(erroru[j,i]),maxerru);
    end
    global maxerrp;
    maxerrp = max(abs(p.values[i] - exactp(x,y)),maxerrp);
end
println("u max error = "*string(maxerru));
println("p max error = "*string(maxerrp));

using Plots
pyplot();
display(plot(Finch.grid_data.allnodes[1,:], Finch.grid_data.allnodes[2,:], erroru[1,:], st=:surface))

# check
# log_dump_config();
# log_dump_prob();

output_values([u,p], "vector2d", format="vtk");

finalize_finch();

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.1
#     language: julia
#     name: julia-1.6
# ---

# %%
struct F{p, T<:Integer} <: Integer
    a::T
    F{p, T}(a::Integer) where {p, T<:Integer} = new{p, T}(mod(T(a), p))
end
F{p}(a::Integer) where p = F{p, typeof(a)}(mod(a, p))
F{p, S}(x::F{p, T}) where {S<:Integer, p, T<:Integer} = F{p, S}(x.a)

Base.promote_rule(::Type{F{p, T}}, ::Type{S}) where {p, T<:Integer, S<:Integer} =
    F{p, promote_type(T, S)}

Base.zero(::Type{F{p, T}}) where {p, T<:Integer} = F{p}(mod(zero(T), p))
Base.one(::Type{F{p, T}}) where {p, T<:Integer} = F{p}(mod(one(T), p))
for op in (:-, :+)
    @eval Base.$op(x::F{p}) where p = F{p}(mod($op(x.a), p))
end
for op in (:-, :+, :*)
    @eval Base.$op(x::F{p}, y::F{p}) where p = F{p}(mod($op(x.a, y.a), p))
end
Base.inv(x::F{p}) where p = F{p}(invmod(x.a, p))
Base.:/(x::F{p}, y::F{p}) where p = x * inv(y)
Base.:\(x::F{p}, y::F{p}) where p = inv(x) * y
Base.:(==)(x::F{p}, y::F{p}) where p = x.a == y.a
Base.:<(x::F{p}, y::F{p}) where p = x.a < y.a

Base.show(io::IO, x::F{p}) where p = print(io, "F", p, '(', x.a, ')')

F7 = F{7, BigInt}
x, y = F7(10), F7(-2)
@show(x, y, zero(x), one(x), +x, -x, x + y, x - y, x * y, x / y, x \ y, x^3, x^-5, x == F7(3), x == 3)
println()
A = F7[1 2; 3 4]
@show(A, 4A, A/4)
println()
using LinearAlgebra
L, U = lu(A)
@show(L, U, L * U, det(A), inv(A), inv(A) * A, A * inv(A));

# %%
"""See https://docs.julialang.org/en/v1/manual/interfaces/"""
Base.iterate(Fp::Type{F{p, T}}) where {p, T<:Integer} = (zero(Fp), zero(T))
function Base.iterate(Fp::Type{F{p, T}}, state) where {p, T<:Integer}
    nextstate = state + 1
    nextstate < p ? (Fp(nextstate), nextstate) : nothing
end
Base.IteratorSize(Fp::Type{F{p, T}}) where {p, T<:Integer} = Base.HasLength()
Base.length(Fp::Type{F{p, T}}) where {p, T<:Integer} = p
Base.eltype(Fp::Type{F{p, T}}) where {p, T<:Integer} = Fp

squares(Fp) = Fp[x^2 for x in Fp]
squareroots(k, Fp) = Fp[x for x in Fp if x^2 == k]
@show(collect(F7), squares(F7), squareroots.(0:6, Ref(F7)));

# %%
using GaloisFields, LinearAlgebra
GF7 = @GaloisField 7
B = GF7[1 2; 3 4]
det(B)

# %%
inv(B)

# %%
lu(B)

# %%
using AbstractAlgebra

@show GF7 = GF(7)
@show P, x = PolynomialRing(GF7, "x")
@show C = GF7[1 2; 3 4]
squares(Fp) = [x^2 for x in Fp]
squareroots(k, Fp) = [x for x in Fp if x^2 == k]
println()
@show(det(C), inv(C), lu(C), charpoly(P, C))
println()
@show(collect(GF7), squares(GF7), squareroots.(0:6, Ref(GF7)));

# %%
@which iterate(GF7)

# %%
@code_warntype iterate(F7, big(1))

# %%
[x for x in F7 if x^2 == 4]

# %%
[x for x in F7 if x^2 == 5]

# %%
@which collect(x for x in F7 if x^2 == 5)

# %%
Base.@default_eltype(x for x in F7 if x^2 == 5)

# %%
Base.@default_eltype(x for x in F7 if x^2 == 5)

# %%
(@macroexpand Base.@default_eltype(x for x in F7 if x^2 == 5)) |> Base.remove_linenums!

# %%
iter = (x for x in F7 if x^2 == 5)

# %%
Base._iterator_upper_bound(iter)

# %%
Core.Compiler.return_type(Base._iterator_upper_bound, Tuple{typeof(iter)})

# %%
@which Base.inferencebarrier(nothing)

# %%
iter2 = (x for x in F7 if x^2 == 4)

# %%
Base._iterator_upper_bound(iter2)

# %%
iter3 = (x for x in GF7 if x^2 == 4)
collect(iter3)

# %%
iter4 = (x for x in GF7 if x^2 == 5)
collect(iter4)

# %%
Core.Compiler.return_type(Base._iterator_upper_bound, Tuple{typeof(iter4)})

# %%
Base._iterator_upper_bound(iter4)

# %%
@code_warntype Base._iterator_upper_bound(iter4)

# %%
@code_warntype Base._iterator_upper_bound(iter2)

# %%
@code_warntype iterate(iter4)

# %%
@code_warntype iterate(iter)

# %%
Core._apply_iterate(Base.iterate, Base.iterate, Core.tuple(iter3.iter), ())

# %%
Core._apply_iterate(Base.iterate, Base.iterate, Core.tuple(iter2.iter), ())

# %%
Core.tuple(iter.iter)

# %%
Core.tuple(iter4.iter)

# %%

max_stars_repo_path: 0010/minimal def of a prime field old.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public