text:
function [ ] = integ(tbv,pib,p0,lv,rd,ub, g, cp, c_sound,rhow,rhou,tb, zu, zw,xu,xs,x_dist_in, ...
    latdmpcoef,raydmpcoef,raydmpz,trigpi,qb,um,u,up,wm,w,wp,thm,th,thp,pim,pic,pip,pprt,...
    qvtot,qvm,qv,qvp,qcm,qc,qcp,qrainm,qrain,qrainp,dx,dz,nt,nz,nx,dt,asscoef,cmixh,cmixv,...
    qc0,k1,k2,thvm,thv,thvp,modeloutputpathandfilename,dtoutput)
#UNTITLED7 Summary of this function goes here
#   Detailed explanation goes here

VT=zeros(nx,nz);
ql=zeros(nx,nz);
tv=zeros(nx,nz);
##
# Set Variable name for a given timestep
thisdump = dtoutput;
# Set var_0 equal to var_init
th_0=th;
u_0=u;
w_0=w;
pi_0=pic;
pprt_0=pprt;
qv_0=qv;
qc_0=qc;
qrain_0=qrain;

th_str = 'th_0';
u_str = 'u_0';
w_str = 'w_0';
pi_str = 'pi_0';
pprt_str = 'pprt_0';
qv_str = 'qv_0';
qc_str = 'qc_0';
qrain_str = 'qrain_0';
dt_str = 'dt';
dx_str = 'dx';
dz_str = 'dz';
nz_str = 'nz';
nx_str = 'nx';
nt_str = 'nt';
qb_str = 'qb';
zu_str = 'zu';
zw_str = 'zw';
xu_str = 'xu';
xs_str = 'xs';
tb_str = 'tb';

# Append New Timesteps to .mat file.
full_str = ['save ' modeloutputpathandfilename ' ' th_str ' ' u_str ' ' w_str ' ' pi_str ' ' pprt_str ' ' qv_str ' ' qc_str ' ' qrain_str ' ' dt_str ' ' dx_str ' ' dz_str ' ' nt_str ' ' nx_str  ' ' nz_str ' ' qb_str ' ' zu_str ' ' zw_str ' ' xu_str ' ' xs_str ' ' tb_str];
eval(full_str);

##
# um(:,:) = ub;
# u(:,:) = ub;
wm[:,:] .= 0.0;
w[:,:] .= 0.0;
# pic(:,:) = 0.0         # Uncomment to make the initial step non-hydrostatic
# pim = pic;
thm = th;
moist=2;
# mtncenter = 25000;

u[260:340,1:3] .= 0;
w[260:340,1:3] .= 0;
u[270:330,4:7] .= 0;
w[270:330,4:7] .= 0;
u[280:320,8:11] .= 0;
w[280:320,8:11] .= 0;
u[290:310,12:14] .= 0;
w[290:310,12:14] .= 0;

# u(560:640,1:3) = 0;
# w(560:640,1:3) = 0;
# u(570:630,4:7) = 0;
# w(570:630,4:7) = 0;
# u(580:620,8:11) = 0;
# w(580:620,8:11) = 0;
# u(590:610,12:14) = 0;
# w(590:610,12:14) = 0;

um = u;
wm = w;

# Forecast Loop

for n=1:nt
    #display(n)

    if n == 1
        d2t = dt;
        dtx = dt/dx;
        dtz = dt/dz;
    else
        d2t = 2*dt;
        dtx = 2.0*dt/dx;
        dtz = 2.0*dt/dz;
    end
    ## Leap-Frog Scheme

    #  do loop for u
    for k=2:nz-1
        for i=2:nx-1
            up[i,k]=um[i,k] - 0.25*dtx*((u[i+1,k]+u[i,k])^2             ...
                -(u[i-1,k]+u[i,k])^2)            ...
                - 0.25*dtz*(rhow[k+1]*(w[i,k+1]+w[i-1,k+1])  ...
                *(u[i,k+1]+u[i  , k ])  ...
                -rhow[ k ]*(w[i, k ]+w[i-1, k ])  ...
                *(u[i,k-1]+u[i  , k ])) ...
                /rhou[k] ...
                - dtx*cp*tbv[k]*(pic[i,k]-pic[i-1,k]);
        end
    end

    #  do loop for w with buoyancy in Fovell's form in terms of theta' and qv'
    #  first prepare for w loop
    for k=1:nz
        for i=1:nx
            if moist == 2
                tv[i,k]=(tb[k]+th[i,k])*(1+0.61*(qb[k]+qv[i,k]))-tbv[k];
                ql[i,k]=qc[i,k]+qrain[i,k];
            else
                tv[i,k]=th[i,k];
                ql[i,k]=0.0;
            end
        end
    end
    #  now do w loop
    for k=3:nz-1
        for i=2:nx-1
            wp[i,k]=wm[i,k] - 0.25*dtx*( (u[i+1,k]+u[i+1,k-1])           ...
                *(w[i+1,k]+w[ i , k ])           ...
                - (u[ i ,k]+u[ i ,k-1])           ...
                *(w[ i ,k]+w[i-1, k ]))          ...
                - 0.25*dtz*(rhou[ k ]*(w[i,k+1]+w[i,k])^2   ...
                -rhou[k-1]*(w[i,k]+w[i,k-1])^2)  ...
                /rhow[k] ...
                - dtz*cp*0.5*(tbv[k]+tbv[k-1])               ...
                *(pic[i,k]-pic[i,k-1])             ...
                + d2t*g*(                                    ...
                0.5*(tv[i,k]/tbv[k]+tv[i,k-1]/tbv[k-1])...
                -0.5*(ql[i,k]+ql[i,k-1]));

        end
    end

    #
    #  do loop for theta using flux form
    #     for k=2:nz-1
    #         for i=2:nx-1
    #             thp[i,k]=thm[i,k] - 0.5*dtx*(u[i+1,k]*(th[i+1,k]+th[i ,k])  ...
    #                 -u[i ,k]*(th[i ,k]+th[i-1,k])) ...
    #                 - 0.5*dtz*( rhow[k+1]*w[i,k+1]             ...
    #                 *(th[i,k+1]+th[i, k ])          ...
    #                 - rhow[ k ]*w[i, k ]             ...
    #                 *(th[i, k ]+th[i,k-1]))         ...
    #                 /rhou[k] ...
    #                 - 0.5*dtz*( rhow[k+1]*w[i,k+1]             ...
    #                 *(tb[k+1] - tb[ k ])            ...
    #                 + rhow[ k ]*w[i, k ]             ...
    #                 *(tb[ k ] - tb[k-1]))           ...
    #                 /rhou[k];
    #
    #         end
    #     end

    #  do loop for theta using flux form except advective form on third term
    for k=2:nz-1
        for i=2:nx-1
            thp[i,k]=thm[i,k] - 0.5*dtx*(u[i+1,k]*(th[i+1,k]+th[ i ,k])  ...
                -u[ i ,k]*(th[ i ,k]+th[i-1,k])) ...
                - 0.5*dtz*( rhow[k+1]*w[i,k+1]             ...
                *(th[i,k+1]+th[i, k ])          ...
                - rhow[ k ]*w[i, k ]             ...
                *(th[i, k ]+th[i,k-1]))         ...
                /rhou[k] ...
                - 0.5*dtz*( w[i,k+1]*(tb[k+1] - tb[ k ])   ...
                +w[i, k ]*(tb[ k ] - tb[k-1]));
        end
    end

    #  do loop for pi
    for k=2:nz-1
        for i=2:nx-1
            pip[i,k]=pim[i,k] - ((c_sound^2)/(rhou[k]*cp*(tbv[k]^2)))   ...
                *( dtx*rhou[k]*tbv[k]*(u[i+1,k]-u[i,k])   ...
                +dtz*0.5*                               ...
                (rhow[k+1]*w[i,k+1]*(tb[k+1]+tb[k]) ...
                -rhow[ k ]*w[i, k ]*(tb[k]+tb[k-1])));
        end
    end

    # Dimensional Pressure
    for k=2:nz-1
        for i=2:nx-1
            pprt[i,k] = pip[i,k]*cp*rhou[k]*tb[k]*(1.0+0.61*qb[k]);
        end
    end

    ## Call Diffusion Scheme

    [up,um] = diffusion(ub,up,um,cmixh,cmixv,dx,dz,dt,nx,nz);
    [wp,wm] = diffusion(zeros(1,nz),wp,wm,cmixh,cmixv,dx,dz,dt,nx,nz);
    [thp,thm] = diffusion(zeros(1,nz),thp,thm,cmixh,cmixv,dx,dz,dt,nx,nz);
    [pip,pim] = diffusion(zeros(1,nz),pip,pim,cmixh,cmixv,dx,dz,dt,nx,nz);

    ## Apply Rayleigh Damping Top to Dry Variables

    for k=2: nz-1
        for i=2: nx-1
            if zu[k] > raydmpz
                [up[i,k]] = rayleighTop(zu[k],zu[nz-1],up[i,k],ub[k],raydmpcoef,raydmpz,trigpi,i,k);
                [wp[i,k]] = rayleighTop(zw[k],zw[nz-1],wp[i,k],0,raydmpcoef,raydmpz,trigpi,i,k);
                [thp[i,k]] = rayleighTop(zu[k],zu[nz-1],thp[i,k],0,raydmpcoef,raydmpz,trigpi,i,k);
                [pip[i,k]] = rayleighTop(zu[k],zu[nz-1],pip[i,k],0,raydmpcoef,raydmpz,trigpi,i,k);
            end

        end
    end

    ## Leapfrog Advection of Moist Variables

    #  do loop for qv using flux form except advective form on third term
    for k=2:nz-1
        for i=2: nx-1
            qvp[i,k]=qvm[i,k] - 0.5*dtx*(u[i+1,k]*(qv[i+1,k]+qv[i ,k])  ...
                -u[i ,k]*(qv[i ,k]+qv[i-1,k])) ...
                - 0.5*dtz*( rhow[k+1]*w[i,k+1]             ...
                *(qv[i,k+1]+qv[i, k ])          ...
                - rhow[ k ]*w[i, k ]             ...
                *(qv[i, k ]+qv[i,k-1]))         ...
                /rhou[k] ...
                - 0.5*dtz*( w[i,k+1]*(qb[k+1] - qb[ k ])   ...
                +w[i, k ]*(qb[ k ] - qb[k-1]));
        end
    end

    #  do loop for qc using flux form except advective form on third term
    for k=2:nz-1
        for i=2:nx-1
            qcp[i,k]=qcm[i,k] - 0.5*dtx*(u[i+1,k]*(qc[i+1,k]+qc[i ,k])  ...
                -u[i ,k]*(qc[i ,k]+qc[i-1,k])) ...
                - 0.5*dtz*( rhow[k+1]*w[i,k+1]             ...
                *(qc[i,k+1]+qc[i, k ])          ...
                - rhow[ k ]*w[i, k ]             ...
                *(qc[i, k ]+qc[i,k-1]))         ...
                /rhou[k];
        end
    end

    #  do loop for qrain using flux form except advective form on third term
    for k=2:nz-1
        for i=2:nx-1

            qrainp[i,k]=qrainm[i,k] - 0.5*dtx*(u[i+1,k]*(qrain[i+1,k]+qrain[i ,k])  ...
                -u[i ,k]*(qrain[i ,k]+qrain[i-1,k])) ...
                - 0.5*dtz*( rhow[k+1]*(w[i,k+1]-           ...
                0.5*(VT[i,k+1]+VT[i,k]))       ...
                *(qrain[i,k+1]+qrain[i, k ])          ...
                - rhow[ k ]*(w[i, k ]-           ...
                0.5*(VT[i,k]+VT[i,k-1]))       ...
                *(qrain[i, k ]+qrain[i,k-1]))         ...
                /rhou[k];
        end
    end

    for k=2: nz-1
        for i= 2: nx-1

            # Compute qvtotal

            qvtot[i,k] = qb[k]+qvp[i,k];

        end
    end

    ## "Fill" any negative water produced by advection of qc, qrain

    qvtot[qvtot<0]=0;

    #  REMOVE ALL NEGATIVE WATER PRODUCED BY ADVECTION
    for k=2:nz-1
        for i=2:nx-1
            qvp[i,k]=max(qb[k]+qvp[i,k],0.0)-qb[k];
            qcp[i,k]=max(qcp[i,k],0.0);
            qrainp[i,k]=max(qrainp[i,k],0.0);
        end
    end

    ## Call Kessler Mircophysics Scheme

    [ qvp,qcp,qrainp,thp,VT ] = kessler(n, qvp,qvtot,qcp,qrainp,dt,qc0,k1,k2,cp,rd,lv,p0,pib,thp,pip,rhou,rhow,tb,qb,VT,nx,nz );

    ##  Apply computational diffusion to moist variables (qv,qc,qrain)

    [qvp,qvm] = diffusion(zeros(1,nz),qvp,qvm,cmixh,cmixv,dx,dz,dt,nx,nz);
    [qcp,qcm] = diffusion(zeros(1,nz),qcp,qcm,cmixh,cmixv,dx,dz,dt,nx,nz);
    [qrainp,qrainm] = diffusion(zeros(1,nz),qrainp,qrainm,cmixh,cmixv,dx,dz,dt,nx,nz);

    ## Apply Raleigh Top to Moist Variables
    for k=2: nz-1
        for i=2: nx - 1
            if zu[k] > raydmpz
                [qvp[i,k]] = rayleighTop(zu[k],zu[nz-1],qvp[i,k],0,raydmpcoef,raydmpz,trigpi,i,k);
                [qcp[i,k]] = rayleighTop(zu[k],zu[nz-1],qcp[i,k],0,raydmpcoef,raydmpz,trigpi,i,k);
                [qrainp[i,k]] = rayleighTop(zu[k],zu[nz-1],qrainp[i,k],0,raydmpcoef,raydmpz,trigpi,i,k);
            end
        end
    end

    ## Call Asselin Filter

    [um,u,up] = asselin(um,u,up,asscoef,nx,nz );
    [wm,w,wp] = asselin(wm,w,wp,asscoef,nx,nz );
    [thm,th,thp] = asselin(thm,th,thp,asscoef,nx,nz );
    [pim,pic,pip] = asselin(pip,pic,pip,asscoef,nx,nz );
    [qvm,qv,qvp] = asselin(qvm,qv,qvp,asscoef,nx,nz );
    [qcm,qc,qcp] = asselin(qcm,qc,qcp,asscoef,nx,nz );
    [qrainm,qrain,qrainp] = asselin(qrainm,qrain,qrainp,asscoef,nx,nz );

    ## Force top and bottom boundary

    # For w-wind

    #     for i = 1 : nx
    wp[:,1] .= 0.0;
    wp[:,2] .= 0.0;
    wp[:,nz] .= 0.0;

    #     end

    ## Use zero gradient at boundaries for u, theta, p, qv,qc, and qrain at
    # top and bottom
    for i=2 : nx-1
        up[i,1] = up[i,2];
        up[i,nz] = up[i,nz-1];
        thp[i,1] = thp[i,2];
        thp[i,nz] = thp[i,nz-1];
        pip[i,1] = pip[i,2];
        pip[i,nz] = pip[i,nz-1];
        qvp[i,1] = qvp[i,2];
        qvp[i,nz] = qvp[i,nz-1];
        qcp[i,1] = qcp[i,2];
        qcp[i,nz] = qcp[i,nz-1];
        qrainp[i,1] = qrainp[i,2];
        qrainp[i,nz] = qrainp[i,nz-1];
    end
    #
#     for k = 1 : nz
#         for i = 1 : nx
#             if (zu[k] < (0.3*xs[i]-4500)) && (xs[i] < 25000) && (xs[i] > 15000)
#                 th[i,k] = th[i,k+1];
#                 pic[i,k] = pic[i,k+1];
#             end
#             if (zu[k] < (-0.3*xs[i]+10500)) && (xs[i] > 25000) && (xs[i] < 35000)
#                 th[i,k] = th[i,k+1];
#                 pic[i,k] = pic[i,k+1];
#             end
#
#             if (zu[k] < (0.3*xu[i]-4500)) && (xu[i] < 25000) && (xu[i] > 15000)
#                 u[i,k] = 0;
#             end
#             if (zu[k] < (-0.3*xu[i]+10500)) && (xu[i] > 25000) && (xu[i] < 35000)
#                 u[i,k] = 0;
#             end
#
#             if (zw[k] < (0.3*xs[i]-4500)) && (xs[i] < 25000) && (xs[i] > 15000)
#                 w[i,k] = 0;
#             end
#             if (zw[k] < (-0.3*xs[i]+10500)) && (xs[i] > 25000) && (xs[i] < 35000)
#                 w[i,k] = 0;
#             end
#         end
#     end

up[260:340,1:3] .= 0;
wp[260:340,1:3] .= 0;
up[270:330,4:7] .= 0;
wp[270:330,4:7] .= 0;
up[280:320,8:11] .= 0;
wp[280:320,8:11] .= 0;
up[290:310,12:14] .= 0;
wp[290:310,12:14] .= 0;

# up(560:640,1:3) = 0;
# wp(560:640,1:3) = 0;
# up(570:630,4:7) = 0;
# wp(570:630,4:7) = 0;
# up(580:620,8:11) = 0;
# wp(580:620,8:11) = 0;
# up(590:610,12:14) = 0;
# wp(590:610,12:14) = 0;

        ## Call Lateral Boundary Sponge

        for k=1: nz
            for i=1: nx

                [up[i,k]] = lateralSponge(up[i,k],ub[k],xu[i],x_dist_in,dx,nx,trigpi,latdmpcoef,1);
                [wp[i,k]] = lateralSponge(wp[i,k],0,xs[i],x_dist_in,dx,nx,trigpi,latdmpcoef,0);
                [thp[i,k]] = lateralSponge(thp[i,k],0,xs[i],x_dist_in,dx,nx,trigpi,latdmpcoef,0);
                [pip[i,k]] = lateralSponge(pip[i,k],0,xs[i],x_dist_in,dx,nx,trigpi,latdmpcoef,0);
                [qvp[i,k]] = lateralSponge(qvp[i,k],0,xs[i],x_dist_in,dx,nx,trigpi,latdmpcoef,0);
                [qcp[i,k]] = lateralSponge(qcp[i,k],0,xs[i],x_dist_in,dx,nx,trigpi,latdmpcoef,0);
                [qrainp[i,k]] = lateralSponge(qrainp[i,k],0,xs[i],x_dist_in,dx,nx,trigpi,latdmpcoef,0);
            end
        end

        ## Apply periodic lateral boundaries for all k

        up[1,:] = up[nx-1,:];
        up[nx,k] = up[2,k];
        wp[1,:] = wp[nx-1,:];
        wp[nx,k] = wp[2,k];
        thp[1,:] = thp[nx-1,:];
        thp[nx,k] = thp[2,k];
        pip[1,:] = pip[nx-1,:];
        pip[nx,k] = pip[2,k] ;
        qvp[1,:] = qvp[nx-1,:];
        qvp[nx,k] = qvp[2,k];
        qcp[1,:] = qcp[nx-1,:];
        qcp[nx,k] = qcp[2,k];
        qrainp[1,:] = qrainp[nx-1,:];
        qrainp(nx,:) = qrainp(2,:);

        ## Create Output Matrix with all times

        timenow = dt*n;

        if timenow >= thisdump

            # Set Variable name for a given timestep
            th_str = ['th_' num2str(n*dt)];
            u_str = ['u_' num2str(n*dt)];
            w_str = ['w_' num2str(n*dt)];
            pi_str = ['pi_' num2str(n*dt)];
            pprt_str = ['pprt_' num2str(n*dt)];
            qv_str = ['qv_' num2str(n*dt)];
            qc_str = ['qc_' num2str(n*dt)];
            qrain_str = ['qrain_' num2str(n*dt)];

            # Set var_n equal to varp

            th_str2 = [th_str '=thp;'];
            u_str2 = [u_str '=up;'];
            w_str2 = [w_str '=wp;'];
            pi_str2 = [pi_str '=pip;'];
            pprt_str2 = [pprt_str '=pprt;'];
            qv_str2 = [qv_str '=qvp;'];
            qc_str2 = [qc_str '=qcp;'];
            qrain_str2 = [qrain_str '=qrainp;'];

            eval(th_str2)
            eval(u_str2)
            eval(w_str2)
            eval(pi_str2)
            eval(pprt_str2)
            eval(qv_str2)
            eval(qc_str2)
            eval(qrain_str2)

            # Append New Timesteps to .mat file.
            full_str = ['save ' modeloutputpathandfilename ' ' th_str ' ' u_str ' ' w_str ' ' pi_str ' ' pprt_str ' ' qv_str ' ' qc_str ' ' qrain_str ' ' '-append'];
            eval(full_str);

            thisdump = thisdump + dtoutput;

        end

        ## Update Variables
        for i = 1: nx
            for k = 1 : nz

                um[i,k] = u[i,k];
                u[i,k] = up[i,k];
                wm[i,k] = w[i,k];
                w[i,k] = wp[i,k];
                pim[i,k] = pic[i,k];
                pic[i,k] = pip[i,k];
                thm[i,k] = th[i,k];
                th[i,k] = thp[i,k];
                #             thvm[i,k] = thv[i,k];
                #             thv[i,k] = thvp[i,k];
                qvm[i,k] = qv[i,k];
                qv[i,k] = qvp[i,k];
                qcm[i,k] = qc[i,k];
                qc[i,k] = qcp[i,k];
                qrainm[i,k] = qrain[i,k];
                qrain[i,k] = qrainp[i,k];

            end
        end

    end # End Forecast Loop

end # End Function

max_stars_repo_name: jeffreycunn/toy-atmospheric-model

max_issues_repo_name: jeffreycunn/toy-atmospheric-model

max_forks_repo_name: jeffreycunn/toy-atmospheric-model


text:
using Pkg
Pkg.activate(@__DIR__)
using SigmaRidgeRegression
using Plots
using StatsBase
using Statistics
using LaTeXStrings
using Random

# To add to tests.jl
# SigmaRidgeRegression.fixed_point_function(hs, γs, [1.0; Inf])
# SigmaRidgeRegression.risk_formula(hs, γs, αs, [1.0; 10000])
_linestyles = [:dot :dashdot :dash]
_main_cols = [:grey :purple :green]

id_design = BlockCovarianceDesign([IdentityCovarianceDesign(), IdentityCovarianceDesign()])

function theoretical_and_realized_mse(
    γs, αs, design::BlockCovarianceDesign; n=400, nreps=50, ntest=20_000
)
    grp = GroupedFeatures(round.(Int, γs .* n))
    design = set_groups(design, grp)
    hs = spectrum.(design.blocks)

    λs1 = SigmaRidgeRegression.optimal_ignore_second_group_λs(γs, αs)
    λs2 = SigmaRidgeRegression.optimal_single_λ(γs, αs)
    λs3 = SigmaRidgeRegression.optimal_λs(γs, αs)

    all_λs = (λs1, λs2, λs3)
    opt_risk_theory = Matrix{Float64}(undef, 1, length(all_λs))
    risk_empirical = Matrix{Float64}(undef, nreps, length(all_λs))

    for (i, λs) in enumerate(all_λs)
        opt_risk_theory[1, i] = SigmaRidgeRegression.risk_formula(hs, γs, αs, λs)
    end

    for j in 1:nreps
        ridge_sim = GroupRidgeSimulationSettings(;
            groups=grp,
            ntrain=n,
            ntest=ntest,
            Σ=design,
            response_model=RandomLinearResponseModel(; αs=αs, grp=grp),
        )
        sim_res = simulate(ridge_sim)

        X_train = sim_res.X[sim_res.resampling_idx[1][1],:]
        Y_train = sim_res.Y[sim_res.resampling_idx[1][1]]

        X_test = sim_res.X[sim_res.resampling_idx[1][2],:]
        Y_test = sim_res.Y[sim_res.resampling_idx[1][2]]

        for (i, λs) in enumerate(all_λs)
            risk_empirical[j, i] = mse_ridge(
                StatsBase.fit(
                    MultiGroupRidgeRegressor(; groups=grp, λs=λs, center=false, scale=false),
                    X_train,
                    Y_train,
                    grp,
                ),
                X_test,
                Y_test,
            )
        end
    end
    risk_empirical = mean(risk_empirical; dims=1)
    return (theoretical=opt_risk_theory, empirical=risk_empirical, all_λs=all_λs)
end

function oracle_risk_plot(
    γs,
    sum_alpha_squared;
    design=id_design,
    ylim=(0, 2.5),
    n=1000,
    title=nothing,
    legend=nothing,
    kwargs...,
)
    ratio_squared = range(0.0, 1.0; length=30)

    αs_squared = ratio_squared .* sum_alpha_squared
    bs_squared = reverse(ratio_squared) .* sum_alpha_squared

    risks = [
        theoretical_and_realized_mse(
            γs, sqrt.([αs_squared[i]; bs_squared[i]]), design; n=n, kwargs...
        ) for i in 1:length(ratio_squared)
    ]
    theoretical_risks = vcat(map(r -> r.theoretical, risks)...) .- 1
    empirical_risks = vcat(map(r -> r.empirical, risks)...) .- 1

    labels =
        [L"$\;$Optimal $\blambda = (\lambda, \infty)$" L"$\;$Optimal $\blambda = (\lambda, \lambda)$" L"$\;$Optimal $\blambda = (\lambda_1, \lambda_2)$"]

    ylabel = L"$\risk{\blambda}- \sigma^2$"
    xlabel = L"\alpha_1^2/(\alpha_1^2 + \alpha_2^2)"
    pl = plot(
        ratio_squared,
        theoretical_risks;
        color=_main_cols,
        linestyle=_linestyles,
        ylim=ylim,
        xguide=xlabel,
        yguide=ylabel,
        legend=legend,
        label=labels,
        background_color_legend=:transparent,
        foreground_color_legend=:transparent,
        grid=false,
        title=title,
        frame=:box,
        plot_titlefontsize=0.5,
        thickness_scaling=2.2,
        legendfontsize=12,
        size=(650, 500),
    )
    plot!(
        pl,
        ratio_squared,
        empirical_risks;
        seriestype=:scatter,
        color=_main_cols,
        markershape=:utriangle,
        markerstrokealpha=0.0,
        markersize=4,
        label=nothing,
    )
    return pl
end

Random.seed!(10)
pgfplotsx()

using PGFPlotsX
empty!(PGFPlotsX.CUSTOM_PREAMBLE)
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage{bm}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\newcommand{\blambda}{\bm{\lambda}}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\newcommand{\risk}[1]{\bm{R}(#1)}")

nreps = 1
title_curve_1 = L"\gamma_1 = \gamma_2 = \frac{1}{4},\;\; \alpha_1^2 + \alpha_2^2 = 1"
curve_1 = oracle_risk_plot(
    [0.25, 0.25], 1.0; legend=:topleft, nreps=nreps, title=title_curve_1
)
plot!(curve_1; tex_output_standalone=true)
savefig(curve_1, "oracle_risk1.tikz")

function generate_risk_plots(base_plot_name; nreps=nreps, kwargs...)
    title_curve_1 = L"\gamma_1 = \gamma_2 = \frac{1}{4},\;\; \alpha_1^2 + \alpha_2^2 = 1"
    curve_1 = oracle_risk_plot(
        [0.25, 0.25], 1.0; legend=:topleft, nreps=nreps, title=title_curve_1, kwargs...
    )

    title_curve_2 = L"\gamma_1 = \frac{1}{10},\; \gamma_2 = \frac{4}{10},\;\; \alpha_1^2 + \alpha_2^2 = 1"
    curve_2 = oracle_risk_plot(
        [0.1, 0.4], 1.0; legend=nothing, nreps=nreps, title=title_curve_2, kwargs...
    )

    title_curve_3 = L"\gamma_1 = \gamma_2 = 1,\;\; \alpha_1^2 + \alpha_2^2 = 1"
    curve_3 = oracle_risk_plot(
        [1.0, 1.0], 1.0; legend=nothing, nreps=nreps, title=title_curve_3, kwargs...
    )

    title_curve_4 = L"\gamma_1 = \gamma_2 = \frac{1}{4},\;\; \alpha_1^2 + \alpha_2^2 = 2"
    curve_4 = oracle_risk_plot(
        [0.25, 0.25], 2.0; legend=nothing, nreps=nreps, title=title_curve_4, kwargs...
    )

    title_curve_5 = L"\gamma_1 = \frac{1}{10},\; \gamma_2 = \frac{4}{10},\;\; \alpha_1^2 + \alpha_2^2 = 2"
    curve_5 = oracle_risk_plot(
        [0.1, 0.4], 2.0; legend=nothing, nreps=nreps, title=title_curve_5, kwargs...
    )

    title_curve_6 = L"\gamma_1 = \gamma_2 = 1,\;\; \alpha_1^2 + \alpha_2^2 = 2"
    curve_6 = oracle_risk_plot(
        [1.0, 1.0], 2.0; legend=nothing, nreps=nreps, title=title_curve_6, kwargs...
    )

    for (i, c) in enumerate([curve_1, curve_2, curve_3, curve_4, curve_5, curve_6])
        savefig(c, "$(base_plot_name)$i.tikz")
    end
end

generate_risk_plots("identity_covariance/oracle_risk.jl")

exponential_design = BlockCovarianceDesign([
    ExponentialOrderStatsCovarianceDesign(; rate=0.5),
    ExponentialOrderStatsCovarianceDesign(; rate=0.5),
])

generate_risk_plots(
    "exponential_covariance/oracle_risk"; design=exponential_design, ylim=(0, 6)
)

ext: jl


text:
using DynamicIterators
using DynamicIterators: dub, _lastiterate, Sampled, State, Controlled
using Test
using Trajectories
import DynamicIterators.dyniterate

@testset "Examples" begin
      include("../example/metropolishastings.jl")
      include("../example/cumsum.jl")

end

c = collect
cf = collectfrom

# Arnolds cat map, repeats after n iterations
# https://en.wikipedia.org/wiki/Arnold%27s_cat_map

@testset "Evolve" begin

      A = [ 0 0 0 0 0 0
            0 1 1 1 1 0
            0 1 0 0 1 0
            0 1 1 1 1 0
            0 1 0 0 1 0
            0 0 0 0 0 0
      ]
      arnold_imap(x, (m,n) = size(x)) = CartesianIndex(1 + (2x[1]+ x[2] - 3)%m, 1 + (x[1] + x[2] - 2) % n)
      arnold(A) = [A[arnold_imap(x, size(A))] for x in CartesianIndices(A)]
      A0 = copy(A)
      P = Evolve(arnold)
      let A = A0
            for i in 1:12
                  A = evolve(P, A)
            end
            @test A == A0
      end

      i, A = evolve(P, (0=>A0), 12)
      @test A == A0
      @test i == 12

      (i, A), _ = dyniterate(P, Control(0=>A0), 12)
      @test A == A0
      @test i == 12

      X = trace(P, (0=>A0), endtime(12))
      @test X isa Trajectory{Array{Int64,1},Array{Array{Int64,2},1}}
      @test keys(X) == 0:12
      @test last(values(X)) == A0

      As = collectfrom(P, A0, 13)
      @test As[1] == As[13]

      @test evolve(1:10, 5) == 6
      @test evolve(1:10, 10) == nothing

      # broken?
      @test evolve(1:10, 11) == nothing

end
@test collect(from(1:14, 10)) == [11, 12, 13, 14]

@testset "time" begin
      @test dyniterate(1:10, Steps(5, 3)) == (8, 8)
      @test dyniterate(TimeLift(1:2:10), NewKey(1=>nothing, 5)) == (5=>1, 5=>1)
end

@testset "Mix" begin
      M = mix((x,y) -> (x+y, y), 0:20000, 0:100)
      F = from(M, (0,0))
      @show x, s = iterate(F)
      @show x, s = iterate(F, s)
      @show x, s = iterate(F, s)
      @test _lastiterate(M, (0,0)) == (100*101÷2 + 100, 100)

      m = collectfrom(M, (0,0))
      @test m[end]  == (100*101÷2 + 100, 100)
      @test eltype(m) == Tuple{Int,Int}

end

@testset "random" begin
      @test all([Randn(0.0)] .== collectfrom(WhiteNoise(), Randn(0.0), 9))

      @test collectfrom(Sampled(WhiteNoise()), (0 => 0.1), 10) isa Array{Pair{Int64,Float64},1}

      @test Base.IteratorEltype(from(WhiteNoise(), Sample(0 => 0.0))) == Base.EltypeUnknown()

      @test Base.IteratorEltype(from(WhiteNoise(), Start(0 => DynamicIterators.Randn(0.0)))) == Base.HasEltype()
      @test eltype(from(WhiteNoise(), Start(0 => 0.0))) == typeof(0 => 0.0)

      @show collectfrom(WhiteNoise(), Sample(Start(0 => 0.0)), 10)
      @test collectfrom(WhiteNoise(), Sample(Start(0 => 0.0)), 10) isa Array{Pair{Int64,Float64},1}

      @test eltype(from(Sampled(WhiteNoise()), (0 => 0.1))) == Pair{Int64,Float64}

      @test eltype(Randn(10)) == Int
      @test eltype(Randn{Int}) == Int

      @test collectfrom(InhomogeneousPoisson(x -> sin(x) + 1, 2.0), (0.0=>0), 10) isa  Array{Pair{Float64,Int64},1}
end

collatz(n) = n % 2 == 0 ? n÷2 : 3n + 1

function bare_collatz(k, n)
      for i in 1:n-1
            k = collatz(k)
      end
      k
end

@testset "control" begin
      F = from(Controlled(1:2:20, Evolve(collatz)), 1=>14)
      ϕ = iterate(F)
      @test (1 => 14, (1, Control(1 => 14))) == ϕ
      ϕ = iterate(F, ϕ[2])
      @test (3 => 22, (3, Control(3 => 22))) == ϕ
      @test collectfrom(F, 1=>14) isa Array{Pair{Int64,Int64},1}

      F = from(Controlled(3:2:20, Evolve(collatz)), 1=>14)
      ϕ = iterate(F)
      @test (3 => 22, (3, Control(3 => 22))) == ϕ
      @test collectfrom(F, 1=>14) isa Array{Pair{Int64,Int64},1}

end
struct Squares <: DynamicIterator
end
dyniterate(S::Squares, (state,)::State) = (state*state, State(state+1))
dyniterate(S::Squares, ::Nothing) = (1, State(2))

@testset "TimeLift" begin

      U = TimeLift(Squares())

      @test collect(from(bind(4:2:8, U), NextKeys(State(1))))== [ 4 => 1
                                                                   6 => 9
                                                                   8 => 25]
end

@testset "Synchronize" begin
      P = InhomogeneousPoisson(x -> sin(x) + 1, 2.0)

      PP = synchronize(P, P)
      u = DynamicIterators.state(0.0 => (0,0), PP)
      @show u
      @show collectfrom(PP, u, 10)

end

@inferred _lastiterate(Evolve(collatz), 1=>6171, endtime(10000) )
@inferred lastiterate(Evolve(collatz), 1=>6171, endtime(10000) )
#using BenchmarkTools

let E = Evolve(collatz), st = endtime(10000)

      @test _lastiterate(E, 1=>6171, st )[2] == lastiterate(E, 1=>6171, st )[2] == bare_collatz(6171, 10000)
      @test @allocated(_lastiterate(E, 1=>6171, st ) ) == 0
      @test @allocated(lastiterate(E, 1=>6171, st ) ) == 0

      @time lastiterate(E, 1=>6171, st)
      @time _lastiterate(E, 1=>6171, st)
      @time bare_collatz(6171, 10000)
end

print("Done")


text:
function Base.isapprox(xs::NTuple{N}, ys::NTuple{N}; kwargs...) where {N}
    for i = 1:N
        if !isapprox(xs[i], ys[i]; kwargs...)
            return false
        end
    end
    true
end

@testset "AffinePatches" begin
    @testset "OrthogonalPatch" begin
        x = ProjectiveVectors.embed(rand(ComplexF64, 2))
        patch = state(OrthogonalPatch(), x)
        @test patch isa HomotopyContinuation.OrthogonalPatchState
        @test nequations(patch) == 1

        init!(patch, x)
        @test norm(patch.v̄) ≈ (1.0,) atol = 1e-15
        u = [0.0im]
        evaluate!(u, patch, x)
        @test u[1] ≈ 0.0 atol = 1e-15

        U = zeros(Complex{Float64}, 1, 3)
        jacobian!(U, patch, x)
        @test U ≈ reshape(conj.(LinearAlgebra.normalize(x)), 1, 3)

        y = ProjectiveVectors.embed(rand(ComplexF64, 2))
        onpatch!(y, patch)
        evaluate!(u, patch, y)
        @test u[1] ≈ 0.0 atol = 1e-15
    end

    @testset "EmbeddingPatch" begin
        x = ProjectiveVectors.embed(rand(ComplexF64, 2))
        patch = state(EmbeddingPatch(), x)
        @test patch isa HomotopyContinuation.EmbeddingPatchState
        @test nequations(patch) == 1

        init!(patch, x)
        @test x[3] ≈ 1.0 atol = 1e-15
        u = [0.0im]
        evaluate!(u, patch, x)
        @test u[1] ≈ 0.0 atol = 1e-15

        U = zeros(Complex{Float64}, 1, 3)
        jacobian!(U, patch, x)
        @test U == reshape([0.0, 0, 1], 1, 3)

        y = ProjectiveVectors.embed(rand(ComplexF64, 2))
        onpatch!(y, patch)
        evaluate!(u, patch, y)
        @test u[1] ≈ 0.0 atol = 1e-15
    end

    @testset "RandomPatch" begin
        x = ProjectiveVectors.embed(rand(ComplexF64, 2))
        patch = state(RandomPatch(), x)
        @test patch isa HomotopyContinuation.RandomPatchState
        @test nequations(patch) == 1

        init!(patch, x)
        u = [0.0im]
        evaluate!(u, patch, x)
        @test u[1] ≈ 0.0 atol = 1e-14

        U = zeros(Complex{Float64}, 1, 3)
        jacobian!(U, patch, x)
        @test U == reshape(patch.v, 1, 3)

        y = ProjectiveVectors.embed(rand(ComplexF64, 2))
        onpatch!(y, patch)
        evaluate!(u, patch, y)
        @test u[1] ≈ 0.0 atol = 1e-14
    end

    @testset "FixedPatch" begin
        x = ProjectiveVectors.embed(rand(ComplexF64, 2))
        patch = state(FixedPatch(), x)
        @test patch isa HomotopyContinuation.FixedPatchState
        @test nequations(patch) == 1

        init!(patch, x)
        u = [0.0im]
        evaluate!(u, patch, x)
        u
        @test u[1] ≈ 0.0 atol = 1e-15

        U = zeros(Complex{Float64}, 1, 3)
        jacobian!(U, patch, x)
        @test U ≈ reshape(patch.v̄, 1, 3)

        y = ProjectiveVectors.embed(rand(ComplexF64, 2))
        onpatch!(y, patch)
        evaluate!(u, patch, y)
        @test u[1] ≈ 0.0 atol = 1e-15
    end
end


text:
module Package1

export Type1, perc

include("file1.jl")
include("file2.jl")

# code
mutable struct Type1
    total
end

perc(a::Type1) = a.total * 0.01

end


text:
# Autogenerated wrapper script for gperftools_jll for x86_64-unknown-freebsd
export libprofiler, libtcmalloc, libtcmalloc_and_profiler, libtcmalloc_debug, libtcmalloc_minimal, libtcmalloc_minimal_debug

using LibUnwind_jll
JLLWrappers.@generate_wrapper_header("gperftools")
JLLWrappers.@declare_library_product(libprofiler, "libprofiler.so.0")
JLLWrappers.@declare_library_product(libtcmalloc, "libtcmalloc.so.4")
JLLWrappers.@declare_library_product(libtcmalloc_and_profiler, "libtcmalloc_and_profiler.so.4")
JLLWrappers.@declare_library_product(libtcmalloc_debug, "libtcmalloc_debug.so.4")
JLLWrappers.@declare_library_product(libtcmalloc_minimal, "libtcmalloc_minimal.so.4")
JLLWrappers.@declare_library_product(libtcmalloc_minimal_debug, "libtcmalloc_minimal_debug.so.4")
function __init__()
    JLLWrappers.@generate_init_header(LibUnwind_jll)
    JLLWrappers.@init_library_product(
        libprofiler,
        "lib/libprofiler.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libtcmalloc,
        "lib/libtcmalloc.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libtcmalloc_and_profiler,
        "lib/libtcmalloc_and_profiler.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libtcmalloc_debug,
        "lib/libtcmalloc_debug.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libtcmalloc_minimal,
        "lib/libtcmalloc_minimal.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libtcmalloc_minimal_debug,
        "lib/libtcmalloc_minimal_debug.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Plot helpers

"Axis unit type. Currently supported: `SI`, `RE`."
@enum AxisUnit SI RE
"Color scales type for 2D plots. Currently supported: `Log`, `Linear`, `SymLog`."
@enum ColorScale Log Linear SymLog

"Plotting arguments."
struct PlotArgs
   "Unit of spatial axis"
   axisunit::AxisUnit
   "data array size"
   sizes::Vector{Int}
   "plotting data range"
   plotrange::Vector{Float32}
   "cell IDs in the cut plane"
   idlist::Vector{Int}
   "mapping from original cell order to cut plane"
   indexlist::Vector{Int}
   "title"
   str_title::String
   "xlabel"
   strx::String
   "ylabel"
   stry::String
   "colorbar title"
   cb_title::String
end

"Set plot-related arguments of `var` in `axisunit`."
function set_args(meta::MetaVLSV, var, axisunit::AxisUnit; normal::Symbol=:none, origin=0.0)
   @unpack ncells, coordmin, coordmax = meta

   if normal == :x
      seq = @SVector [2,3]
      dir = 1
   elseif normal == :y || (ncells[2] == 1 && ncells[3] != 1) # polar
      seq = @SVector [1,3]
      dir = 2
   elseif normal == :z || (ncells[3] == 1 && ncells[2] != 1) # ecliptic
      seq = @SVector [1,2]
      dir = 3
   else
      throw(ArgumentError("1D data detected. Please use 1D plot functions."))
   end

   sizes = ncells[seq]
   plotrange = [coordmin[seq[1]], coordmax[seq[1]], coordmin[seq[2]], coordmax[seq[2]]]
   axislabels = ['X', 'Y', 'Z'][seq]

   if normal == :none
      idlist, indexlist = Int[], Int[]
   else
      idlist, indexlist = let sliceoffset = origin - coordmin[dir]
         getslicecell(meta, sliceoffset, dir, coordmin[dir], coordmax[dir])
      end
   end

   # Scale the sizes to the highest refinement level
   sizes *= 2^meta.maxamr # data needs to be refined later

   unitstr = axisunit == RE ? L"$R_E$" : L"$m$"
   strx = axislabels[1]*"["*unitstr*"]"
   stry = axislabels[2]*"["*unitstr*"]"

   str_title = @sprintf "t= %4.1fs" meta.time

   datainfo = readvariablemeta(meta, var)

   cb_title = !isempty(datainfo.variableLaTeX) ?
      datainfo.variableLaTeX * " ["*datainfo.unitLaTeX*"]" : ""

   PlotArgs(axisunit, sizes, plotrange, idlist, indexlist, str_title, strx, stry, cb_title)
end

"Set colormap limits for `data`."
function set_lim(vmin, vmax, data, colorscale::ColorScale=Linear)
   if colorscale in (Linear, SymLog)
      v1 = isinf(vmin) ? minimum(x->isnan(x) ? +Inf : x, data) : vmin
      v2 = isinf(vmax) ? maximum(x->isnan(x) ? -Inf : x, data) : vmax
   else # logarithmic
      datapositive = data[data .> 0.0]
      v1 = isinf(vmin) ? minimum(datapositive) : vmin
      v2 = isinf(vmax) ? maximum(x->isnan(x) ? -Inf : x, data) : vmax
   end

   v1, v2
end

"Return x and y ranges for 2D."
function get_axis(axisunit::AxisUnit, plotrange, sizes)
   if axisunit == RE
      x = LinRange(plotrange[1], plotrange[2], sizes[1]) ./ Vlasiator.Re
      y = LinRange(plotrange[3], plotrange[4], sizes[2]) ./ Vlasiator.Re
   else
      x = LinRange(plotrange[1], plotrange[2], sizes[1])
      y = LinRange(plotrange[3], plotrange[4], sizes[2])
   end
   x, y
end

get_axis(pArgs::PlotArgs) = get_axis(pArgs.axisunit, pArgs.plotrange, pArgs.sizes)

"Obtain data from `meta` of `var` for 2D plotting. Use `op` to select vector components."
function prep2d(meta::MetaVLSV, var, op=:none)
   dataRaw = Vlasiator.getdata2d(meta, var)

   data =
      if ndims(dataRaw) == 3
         if op in (:x, :1)
            @view dataRaw[1,:,:]
         elseif op in (:y, :2)
            @view dataRaw[2,:,:]
         elseif op in (:z, :3)
            @view dataRaw[3,:,:]
         elseif op == :mag
            @views hypot.(dataRaw[1,:,:], dataRaw[2,:,:], dataRaw[3,:,:])
         end
      else
         dataRaw
      end
   data
end

"Return `data` of `var` on a uniform 2D mesh on the finest AMR level."
function prep2dslice(meta::MetaVLSV, var, normal, op, pArgs::PlotArgs)
   @unpack idlist, indexlist = pArgs

   data3D = readvariable(meta, var)

   if startswith(var, "fg_") # field quantities, fsgrid
      throw(ArgumentError("FS grid variable $var plotting in cut currently not supported!"))
   else # moments, dccrg grid
      # vlasov grid, AMR
      if ndims(data3D) == 1
         data2D = data3D[indexlist]

         data = refineslice(meta, idlist, data2D, normal)
      elseif ndims(data3D) == 2
         data2D = data3D[:,indexlist]

         if op in (:x, :y, :z, :1, :2, :3)
            if op in (:x, :1)
               slice = @view data2D[1,:]
            elseif op in (:y, :2)
               slice = @view data2D[2,:]
            elseif op in (:z, :3)
               slice = @view data2D[3,:]
            end
            data = refineslice(meta, idlist, slice, normal)
         elseif op == :mag
            datax = @views refineslice(meta, idlist, data2D[1,:], normal)
            datay = @views refineslice(meta, idlist, data2D[2,:], normal)
            dataz = @views refineslice(meta, idlist, data2D[3,:], normal)
            data = hypot.(datax, datay, dataz)
         end
      end
   end
   data
end

"""
    prep_vdf(meta::MetaVLSV, location; kwargs...)

Return the cell velocities `v1, v2`, bin ranges `r1, r2`, cell VDF values `fweight`,
and strings of labels and titles for VDF plots.

# Optional arguments
- `unit::AxisUnit`: location unit in `SI`, `RE`.
- `unitv::String`: velocity unit in ("km/s", "m/s").
- `limits::Vector{Real}`: velocity space range given in [xmin, xmax, ymin, ymax].
- `slicetype`: symbol for choosing the slice type from :xy, :xz, :yz, :bperp, :bpar, :bpar1.
- `center`: symbol for setting the reference frame from :bulk, :peak.
- `vslicethick`: setting the velocity space slice thickness in the normal direction. If set
to 0, the whole distribution along the normal direction is projected onto a plane. Currently
this is only meaningful when `center` is set such that a range near the bulk/peak normal
velocity is selected!
- `weight::Symbol`: choosing distribution weights from phase space density or particle flux
between `:particle` and `:flux`.
- `flimit`: minimum VDF threshold for plotting.
- `verbose`: display the selection process.
"""
function prep_vdf(meta::MetaVLSV, location;
   species="proton", unit::AxisUnit=SI, unitv="km/s", slicetype=:nothing, vslicethick=0.0,
   center=:nothing, weight=:particle, flimit=-1.0, verbose=false)

   ncells = meta.ncells

   if haskey(meta.meshes, species)
      vmesh = meta.meshes[species]
   else
      throw(ArgumentError("Unable to detect population $species"))
   end

   unit == RE && (location .*= Re)

   # Set unit conversion factor
   unitvfactor = unitv == "km/s" ? 1f3 : 1.0f0

   # Get closest cell ID from input coordinates
   cidReq = getcell(meta, location)
   cidNearest = getnearestcellwithvdf(meta, cidReq)

   # Set normal direction
   if slicetype == :nothing
      slicetype =
         if ncells[2] == 1 && ncells[3] == 1 # 1D, select xz
            :xz
         elseif ncells[2] == 1 # polar
            :xz
         elseif ncells[3] == 1 # ecliptic
            :xy
         end
   end

   if slicetype == :xy
      dir1, dir2, dir3 = 1, 2, 3
      sliceNormal = [0., 0., 1.]
   elseif slicetype == :xz
      dir1, dir2, dir3 = 1, 3, 2
      sliceNormal = [0., 1., 0.]
   elseif slicetype == :yz
      dir1, dir2, dir3 = 2, 3, 1
      sliceNormal = [1., 0., 0.]
   elseif slicetype in (:bperp, :bpar, :bpar1)
      if hasvariable(meta, "B_vol")
         B = readvariable(meta, "B_vol", cidNearest)
      elseif hasvariable(meta, "vg_b_vol")
         B = readvariable(meta, "vg_b_vol", cidNearest)
      end
      BxV = B × Vbulk
      if slicetype == :bperp # slice in b_perp1/b_perp2
         sliceNormal = B ./ norm(B)
      elseif slicetype == :bpar1 # slice in b_parallel/b_perp1 plane
         sliceNormal = B × BxV
         sliceNormal ./= norm(sliceNormal)
      else # slice in b_parallel/b_perp2 plane
         sliceNormal = BxV ./ norm(BxV)
      end
   end

   v1size = vmesh.vblocks[dir1] * vmesh.vblock_size[dir1]
   v2size = vmesh.vblocks[dir2] * vmesh.vblock_size[dir2]

   v1min, v1max = vmesh.vmin[dir1], vmesh.vmax[dir1]
   v2min, v2max = vmesh.vmin[dir2], vmesh.vmax[dir2]

   @assert (v1max - v1min) / v1size ≈ (v2max - v2min) / v2size "Noncubic vgrid detected!"
   cellsize = (v1max - v1min) / v1size

   vcellids, vcellf = readvcells(meta, cidNearest; species)

   V = getvcellcoordinates(meta, vcellids; species)

   if center == :bulk # centered with bulk velocity
      Vbulk =
         if hasvariable(meta, "moments")
            # From a restart file
            readvariable(meta, "restart_V", cidNearest)
         elseif hasvariable(meta, species*"/vg_v")
            # Vlasiator 5
            readvariable(meta, species*"/vg_v", cidNearest)
         elseif hasvariable(meta, species*"/V")
            readvariable(meta, species*"/V", cidNearest)
         else
            readvariable(meta, "V", cidNearest)
         end

      for i in eachindex(V)
         V[i] = V[i] .- Vbulk
      end
   elseif center == :peak # centered on highest VDF-value
      Vpeak = maximum(vcellf)
      for i in eachindex(V)
         V[i] = V[i] .- Vpeak
      end
   end

   # Set sparsity threshold
   if flimit < 0
      flimit =
         if hasvariable(meta, species*"/vg_effectivesparsitythreshold")
            readvariable(meta, species*"/vg_effectivesparsitythreshold", cidNearest)
         elseif hasvariable(meta, species*"/EffectiveSparsityThreshold")
            readvariable(meta, species*"/EffectiveSparsityThreshold", cidNearest)
         else
            1e-16
         end
   end

   # Drop velocity cells which are below the sparsity threshold
   findex_ = vcellf .≥ flimit
   fselect = vcellf[findex_]
   Vselect = V[findex_]

   if slicetype ∈ (:xy, :yz, :xz)
      v1select = [v[dir1] for v in Vselect]
      v2select = [v[dir2] for v in Vselect]
      vnormal = [v[dir3] for v in Vselect]

      strx, stry = getindex(["vx", "vy", "vz"], [dir1, dir2])
   elseif slicetype ∈ (:Bperp, :Bpar, :Bpar1)
      #TODO: NOT working yet!
      if slicetype == :Bperp
         v1select = Vrot2[1,:] # the X axis of the slice is BcrossV=perp1
         v2select = Vrot2[2,:] # the Y axis of the slice is Bcross(BcrossV)=perp2
         vnormal = Vrot2[3,:] # the Z axis of the slice is B
         strx = L"$v_{B \times V}$"
         stry = L"$v_{B \times (B \times V)}$"
      elseif slicetype == :Bpar
         v1select = Vrot2[3,:] # the X axis of the slice is B
         v2select = Vrot2[2,:] # the Y axis of the slice is Bcross(BcrossV)=perp2
         vnormal = Vrot2[1,:] # the Z axis of the slice is -BcrossV=perp1
         strx = L"$v_{B}$"
         stry = L"$v_{B \times V}$"
      elseif slicetype == :Bpar1
         v1select = Vrot2[3,:] # the X axis of the slice is B
         v2select = Vrot2[1,:] # the Y axis of the slice is BcrossV=perp1
         vnormal = Vrot2[2,:] # the Z axis of the slice is Bcross(BcrossV)=perp2
         strx = L"$v_{B}$"
         stry = L"$v_{B \times (B \times V)}$"
      end
   end

   let str = " [$unitv]"
      strx *= str
      stry *= str
   end

   if vslicethick < 0 # Set a proper thickness
      vslicethick =
         if any(sliceNormal .== 1.0) # Assure that the slice cut through at least 1 vcell
            cellsize
         else # Assume cubic vspace grid, add extra space
            cellsize*(√3+0.05)
      end
   end

   # Weights using particle flux or phase-space density
   fweight = weight == :flux ? fselect*norm([v1select, v2select, vnormal]) : fselect

   # Select cells within the slice area
   if vslicethick > 0.0
      v1, v2, fweight = let ind_ = @. (abs(vnormal) ≤ 0.5*vslicethick)
         v1select[ind_], v2select[ind_], fweight[ind_]
      end
   else
      v1, v2 = v1select, v2select
   end

   v1 ./= unitvfactor
   v2 ./= unitvfactor

   r1 = LinRange(v1min, v1max, v1size+1) / unitvfactor
   r2 = LinRange(v2min, v2max, v2size+1) / unitvfactor

   str_title = @sprintf "t = %4.1fs" meta.time

   if verbose
      @info "Original coordinates : $location"
      @info "Original cell        : $(getcellcoordinates(meta, cidReq))"
      @info "Nearest cell with VDF: $(getcellcoordinates(meta, cidNearest))"
      let
         x, y, z = getcellcoordinates(meta, cidNearest)
         @info "cellid $cidNearest, x = $x, y = $y, z = $z"
      end

      if center == :bulk
         @info "Transforming to plasma frame, travelling at speed $Vbulk"
      elseif center == :peak
         @info "Transforming to peak f-value frame, travelling at speed $Vpeak"
      end

      @info "Using VDF threshold value of $flimit."

      if vslicethick > 0
         @info "Performing slice with a counting thickness of $vslicethick"
      else
         @info "Projecting total VDF to a single plane"
      end
   end

   v1, v2, r1, r2, fweight, strx, stry, str_title
end


text:
using BinDeps

@BinDeps.setup

cudautils   = library_dependency("cudautils")
linscan_aqd = library_dependency("linscan_aqd", aliases=["linscan_aqd","linscan_aqd.so"])
linscan_aqd_pairwise_byte = library_dependency("linscan_aqd_pairwise_byte", aliases=["linscan_aqd_pairwise_byte","linscan_aqd_pairwise_byte.so"])
encode_icm_so = library_dependency("encode_icm_so", aliases=["encode_icm_so", "encode_icm_so.so"])

deps = [cudautils, linscan_aqd, linscan_aqd_pairwise_byte, encode_icm_so]

prefix=joinpath(BinDeps.depsdir(linscan_aqd))
linscan_aqdbuilddir = joinpath(BinDeps.depsdir(linscan_aqd),"builds")

# === CPP linscan code ===
# TODO maybe refactor this?
provides(BuildProcess,
    (@build_steps begin
        CreateDirectory(linscan_aqdbuilddir)
        @build_steps begin
            ChangeDirectory(linscan_aqdbuilddir)
            FileRule(joinpath(prefix,"builds","linscan_aqd.so"),@build_steps begin
                `g++ -O3 -shared -fPIC ../src/linscan_aqd.cpp -o linscan_aqd.so -fopenmp`
            end)
        end
    end),linscan_aqd, os = :Unix, installed_libpath=joinpath(prefix,"builds"))

provides(BuildProcess,
    (@build_steps begin
        CreateDirectory(linscan_aqdbuilddir)
        @build_steps begin
            ChangeDirectory(linscan_aqdbuilddir)
            FileRule(joinpath(prefix,"builds","linscan_aqd_pairwise_byte.so"),@build_steps begin
                `g++ -O3 -shared -fPIC ../src/linscan_aqd_pairwise_byte.cpp -o linscan_aqd_pairwise_byte.so -fopenmp`
            end)
        end
    end),linscan_aqd_pairwise_byte, os = :Unix, installed_libpath=joinpath(prefix,"builds"))

provides(BuildProcess,
    (@build_steps begin
        CreateDirectory(linscan_aqdbuilddir)
        @build_steps begin
            ChangeDirectory(linscan_aqdbuilddir)
            FileRule(joinpath(prefix,"builds","encode_icm_so.so"),@build_steps begin
                `rm -f encode_icm_so.so`
                `g++ -O3 -shared -fPIC ../src/encode_icm.cpp -o encode_icm_so.so -fopenmp`
            end)
        end
    end),encode_icm_so, os = :Unix, installed_libpath=joinpath(prefix,"builds"))

# === CUDA code ===
provides(BuildProcess,
    (@build_steps begin
        CreateDirectory(linscan_aqdbuilddir)
        @build_steps begin
            ChangeDirectory(linscan_aqdbuilddir)
            FileRule(joinpath(prefix,"builds","cudautils.so"),@build_steps begin
                `/usr/local/cuda/bin/nvcc -ptx ../src/cudautils.cu -o cudautils.ptx -arch=compute_35`
                `/usr/local/cuda/bin/nvcc --shared -Xcompiler -fPIC -shared ../src/cudautils.cu -o cudautils.so -arch=compute_35`
            end)
        end
    end),cudautils, os = :Unix, installed_libpath=joinpath(prefix,"builds"))

@BinDeps.install Dict([(:linscan_aqd, :linscan_aqd),
                      (:linscan_aqd_pairwise_byte, :linscan_aqd_pairwise_byte),
                      (:encode_icm_so, :encode_icm_so),
                      (:cudautils, :cudautils)])


text:
using MOEA

method = DE()

result = optimize(sphere,
                  method,
		  population = [Individual([2.0, 2.0]) for i in 1:method.populationSize],
		  options = MOEA.Options(store_trace=true))


text:
## Various possible plot interfaces

## Plot interfaces for functions

## plot(x, [y], args..., kwargs...) lineplot
## plot(x::Tuple{Vector}, args...; symboltype::String="o", kwargs...) scatter plot
## plot(f, a::Real, b::Real, args...; kwargs...)  function plot using adaptive point, a, b length two atleast
## plot(fs::Tuple{Function}, a, b; kwargs) parametric plot
## plot(fs::Vector{Function}, a::Real, b::Real, args...; kwargs...) function plot, overlay
## plot(fs::Array{Function, 2}, a::Real, b::Real, args...; kwargs...) table of plots

errs_to_nan(f) = (x) -> try f(x) catch e NaN end

typealias ScatterPlotPoints{T<:Real, S<:Real} (Vector{T}, Vector{S})

function _plot(p::FramedPlot, x::ScatterPlotPoints, args...; symbolkind="circle",kwargs...)
    #XXX: check if args have symbolkind options
    _plot(p, x[1], x[2], args...; symbolkind=symbolkind, kwargs...)
end

function _plot(p::FramedPlot, f::Function, a::Real, b::Real, args...;  kwargs...)
    xs = adaptive_points(f, a, b)
    ys = map(errs_to_nan(f), xs)
    _plot(p, xs, ys, args...; kwargs...)
end

## multiple plots on one
## kwargs vectorized, not recycled
## e.g.:  plot([sin, cos], 0, 2pi, color=["blue", "red"])
function _plot(p::FramedPlot, fs::Vector{Function}, a::Real, b::Real, args...; kwargs...)

    f = fs[1]

    xs = adaptive_points(f, a, b)
    ys = map(errs_to_nan(f), xs)
    kws = [(k, v[1]) for (k,v) in kwargs]
    _plot(p, xs, ys, args...; kws...)

    for i in 2:length(fs)
        xs = adaptive_points(fs[i], a, b)
        ys = map(errs_to_nan(fs[i]), xs)
        kws = [(k, v[i]) for (k,v) in kwargs]
        _plot(p, xs, ys, args...;  kws...)
    end

    p

end

## parametric plot
typealias ParametricFunctionPair (Function, Function)
function _plot(p::FramedPlot, fs::ParametricFunctionPair, a::Real, b::Real, args...; npoints::Int=500, kwargs...)
    us = linspace(a, b, npoints)
    xs = map(errs_to_nan(fs[1]), us)
    ys = map(errs_to_nan(fs[2]), us)
    _plot(p, xs, ys, args...; kwargs...)
end

## Array
## kwargs are vectorized (without recycling)
## e.g.:  plot([sin cos]', 0, 2pi, color=["blue" "red"]')
function plot(fs::Array{Function, 2}, a::Real, b::Real, args...; kwargs...)
    m,n = size(fs)
    tbl = Table(m, n)
    for i in 1:m, j in 1:n
        f = fs[i,j]
        kws = [(k, v[i,j]) for (k, v) in kwargs]
        p = plot(f, a, b, args...; kws...)
        tbl[i,j] = p
    end

    tbl
end

## adaptive plotting
## algorithm from http://yacas.sourceforge.net/Algochapter4.html
## use of heaps follows roughly quadgk.jl

using Base.Collections
import Base.isless, Base.Order.Reverse

immutable Segment
    a::Number
    b::Number
    depth::Int
    E::Real
end

isless(i::Segment, j::Segment) = isless(i.E, j.E)

function evalrule(f, a, b; depth::Int=0)
    xs = linspace(a, b, 7)[2:6] # avoid edges?
    y = [try f(x) catch e NaN end for x in xs]

    if all(map(isnan,y))
        error("Function does not evaluate to a real number at initial set of points in the interval ($a, $b)")
    end

    wiggles(x,y,z) = any(map(u->isinf(u) | isnan(u), [x,y,z])) || (y < min(x,z)) || (y > max(x,z)) ? 1 : 0
    wiggly = [wiggles(y[i:i+2]...) for i in 1:3]

    if depth < 0
        E = 0
    elseif sum(wiggly) > 2
        E = 1
    else
        ## not too wiggly, but may not approximate well enough
        g = y - minimum(y)
        E = (1/3)*abs( (y[1] + 4*y[2] + y[3]) - (y[3] + 4*y[4] + y[5])) # no h = (b-a)/2
    end
    Segment(a, b, depth, E)
end

function adaptive_points(f::Function, a::Real, b::Real;
                             tol::Real=1e-3,
                             max_depth::Int=6)

    n = 100
    s = (a+b)/2 + (b-a)/2* cos((n:-1:0) * pi/n) # non even, to avoid antialiasing. Overkill?
    segs = Segment[]
    for i in 1:n
        heappush!(segs, evalrule(f, s[i], s[i+1], depth=max_depth), Reverse)
    end

    E = segs[1].E

    while E > tol
        s = heappop!(segs, Reverse)
        mid = (s.a + s.b) * 0.5
        s1 = evalrule(f, s.a, mid, depth=s.depth - 1)
        s2 = evalrule(f, mid, s.b, depth=s.depth - 1)
        heappush!(segs, s1, Reverse)
        heappush!(segs, s2, Reverse)
        E = segs[1].E
    end

    x = Float64[]
    [append!(x, [s.a, s.b]) for s in segs]
    x = sort!(x)
    x
end

## Contour plot -- but this is too slow to be usable
## algorithm from MATLAB
## cf. http://www.mathworks.com/help/matlab/creating_plots/contour-plots.html
type Contourc
    contours
end

function contourc(f::Function, x, y; cs::Union(Nothing, Number)=nothing)

    fxy = [f(x,y) for x in x, y in y]

    ## we have edges 1,2,3,4
    ## with 1 on bottom, 2 on right, 3 top, 4 left. So
    ## exiting 1 goes to 3; 2->4, 3->1 and 4->2

    ## helper functions
    prop(c, z0, z1) = (c - z0)/(z1-z0)
    function interp_square(c, i,j)
        ## worry of i,j on boundary of 1:length(x), length(y)
        if (i < 1 || i >= length(x)) || (j < 1 || j >= length(y))
            return -1 * ones(4)
        else
            [prop(c, fxy[i,j], fxy[i+1,j]),
             prop(c, fxy[i+1,j], fxy[i+1, j+1]),
             prop(c, fxy[i,j+1], fxy[i+1, j+1]),
             prop(c, fxy[i,j], fxy[i, j+1])]
        end
    end
    insquare(x) =  any(0 .<= x .<= 1)

    function interp_point(edge, i, j, t)
        if edge == 1
            newx = x[i] + t*(x[i+1] - x[i])
            newy = y[j]
        elseif edge == 2
            newx = x[i+1]
            newy = y[j] + t*(y[j+1] - y[j])
        elseif edge == 3
            newx = x[i] + t*(x[i+1] - x[i])
            newy = y[j+1]
        else
            newx = x[i]
            newy = y[j] + t*(y[j+1] - y[j])
        end
        (newx, newy)
    end

    function which_next(edge, i, j)
        if edge == 1
            (i, j-1)
        elseif edge == 2
            (i+1, j)
        elseif edge == 3
            (i, j+1)
        else
            (i-1, j)
        end
    end

    function next_square(c, i, j, enter_edge, cx, cy, m)
##        println("Chasing $i $j")

        sq = interp_square(c, i, j)
        w = 0 .<= sq .<= 1
        ## check if 2 or more (saddle point) XXX

        ## what is next edge?
        if enter_edge == 1
            next_edge = setdiff((1:4)[w], 3)
        elseif enter_edge == 2
            next_edge = setdiff((1:4)[w], 4)
        elseif enter_edge == 3
            next_edge = setdiff((1:4)[w], 1)
        else
            next_edge = setdiff((1:4)[w], 2)
        end

        if length(next_edge) == 0
            return (cx, cy)
        end

        next_edge = next_edge[1]

        ##
        newx, newy = interp_point(next_edge, i, j, sq[next_edge])
        push!(cx, newx); push!(cy, newy)

        if m[i,j] == 1
            ## already visited
            return (cx, cy)
        end
        m[i,j] = 1
        next_square(c, which_next(next_edge, i,j)..., next_edge, cx, cy, m)
    end

    function chase_square(c, i, j, m)
##        println("chase $i $j")

        sq = interp_square(c, i, j)
        w = 0 .<= sq .<= 1

        ## chase both edges, then put together
        edges = [1:4][w]
        ## should be 2, might be more (saddle point)

        m[i,j] = 1                  # visited

        out = map(edges) do edge
            cx = Float64[]; cy = Float64[]
            t = sq[edge]
            newx, newy = interp_point(edge, i,j,t)
            push!(cx, newx); push!(cy, newy)
            next_square(c, which_next(edge, i,j)..., edge, cx, cy, m)
        end
        ## out is array of tuples
        if !isa(out[1], Nothing)
            ([reverse(out[1][1]), out[2][1]], [reverse(out[1][2]), out[2][2]])
        else
            nothing
        end

    end

## for each level to plot
    if isa(cs, Nothing)
        cs = linspace(min(fxy), max(fxy), 7+2)[2:8]
    else
        cs = [cs]
    end
    contours = {}

    for c in cs
        m = zeros(Int, size(fxy)...)
        c_contours = {}

        for i in 2:length(x)-1, j in 2:length(y)-1
            if m[i,j] == 1
                continue
            end
            sq = interp_square(c, i, j)
            if insquare(sq)
                path = chase_square(c, i,j, m)
##                println("Chased path:", path)
                if !isa(path, Nothing)
                    push!(c_contours, path)
                end
            else
                m[i,j] = 1              # visited
            end
        end
        push!(contours, (c, c_contours))
    end
    Contourc(contours)
end

## contour plot
function plot(f::Contourc; kwargs...)
    p = FramedPlot()
    ## out is array of tuples (c, c_countrous)
    for clevels in f.contours
        for contours in clevels[2]
            add(p, Curve(contours[1], contours[2]))
        end
    end
    for (k, v) in kwargs
        setattr(p, k, v)
    end

    p
end

# ## Some tests
# using Winston
# ## basic
# plot(sin, 0, 2pi) |> Winston.tk
# plot(x -> sin(1/x), 0, 1) |> Winston.tk

# ## attributes
# plot(sin, 0, 2pi, title="Title") |> Winston.tk

# ## parametric
# plot((x -> sin(2x), x -> cos(3x)), 0, 2pi) |> Winston.tk

# ## vector
# plot([sin, cos], 0, 2pi) |> Winston.tk

# ## table with arguments
# plot([sin cos; x->-sin(x) x->-cos(x)], 0, 2pi, title=["f" "f'(x)"; "f''(x)" "f'''(x)"]) |> Winston.tk

# ## contour plot
# f(x,y) = sin(x)*sin(y)
# x = linspace(-pi, pi, 50)
# y = linspace(-pi, pi, 50)
# c = Winston.contourc(f, x, y)
# plot(c) |> Winston.tk


text:
function solve end
function solve! end

function _show(io::IO, object)
    print(io, nameof(typeof(object)), ":")
    pad = get(io, :pad, "")
    lvl = get(io, :lvl, 1)
    names = propertynames(object)
    N = length(names)
    for (n, name) in enumerate(names)
        lvl_ = lvl
        field = getproperty(object, name)
        if field !== nothing
            print(io, "\n", pad, "   ‣ " * string(name) * " ≔ ")
            if lvl < 2
                lvl_ += 1
                show(IOContext(io, :pad => string(pad, "   "), :lvl => lvl_), field)
            else
                summary(io, field)
            end
        end
    end
end

_summary(io::IO, object) = print(io, nameof(typeof(object)))


text:
using Test

include("ArithmeticCoding.jl")

using Main.ArithmeticCoding: encode, BitVector

using Main.ArithmeticCoding.BinaryFractions: binaryfraction, tobitstring, nbits, tofloat, getbit, BinaryFraction, shortestbetween, bitvector

using Main.ArithmeticCoding.SymbolDistributions: SymbolDistribution, getp, getinterval

@testset "tobitstring" begin
    b = "1"
    expected = "1" * repeat("0", nbits - 1)
    actual = tobitstring(b)
    @test length(actual) == nbits
    @test expected == actual

    b = "010010"
    expected = "010010" * repeat("0", nbits - 6)
    actual = tobitstring(b)
    @test length(actual) == nbits
    @test expected == actual
end

@testset "binaryfraction" begin
    d = .5
    expansion = "1"
    expected = tobitstring(expansion)
    actual = bitstring(binaryfraction(d))
    @test expected == actual

    d = 1/2 + 1/4
    expansion = "11"
    expected = tobitstring(expansion)
    actual = bitstring(binaryfraction(d))
    @test expected == actual

    d = 1/32 + 1/64
    expansion = "000011"
    expected = tobitstring(expansion)
    actual = bitstring(binaryfraction(d))
    @test expected == actual

    d = 0.0
    expansion = "0"
    expected = tobitstring(expansion)
    actual = bitstring(binaryfraction(d))
    @test expected == actual

    d = 0
    expansion = "0"
    expected = tobitstring(expansion)
    actual = bitstring(binaryfraction(d))
    @test expected == actual

    d = 1
    expected = repeat("1", nbits)
    actual = bitstring(binaryfraction(d))
    @test expected == actual

    @test_throws ErrorException binaryfraction(-.1)
    @test_throws ErrorException binaryfraction(1.0)
    @test_throws ErrorException binaryfraction(-1)
    @test_throws ErrorException binaryfraction(2)
end

@testset "getbit" begin
    d = .5
    bf = binaryfraction(d)
    bit = 1
    expected = true
    actual = getbit(bf, bit)
    @test expected == actual

    bit = 2
    expected = false
    actual = getbit(bf, bit)
    @test expected == actual
end

@testset "shortestbetween" begin
    upper = BinaryFraction(0b110)
    lower = BinaryFraction(0b010)
    short = BinaryFraction(0b011)
    expected = bitvector(short)
    actual = shortestbetween(lower, upper)
    @test expected == actual

    upper = BinaryFraction(0b1100)
    lower = BinaryFraction(0b1010)
    short = BinaryFraction(0b1011)
    expected = bitvector(short)
    actual = shortestbetween(lower, upper)
    @test expected == actual
end

@testset "tofloat" begin
    d = .5
    bf = binaryfraction(d)
    @test tofloat(bf) ≈ d

    d = 1
    bf = binaryfraction(d)
    @test tofloat(bf) ≈ d

    d = 0.0
    bf = binaryfraction(d)
    @test tofloat(bf) ≈ d

    d = 0.1
    bf = binaryfraction(d)
    @test tofloat(bf) ≈ d

    d = rand()
    bf = binaryfraction(d)
    @test tofloat(bf) ≈ d
end

@testset "bf_arithmetic" begin
    d1 = .1
    d2 = .2
    bf1 = binaryfraction(d1)
    bf2 = binaryfraction(d2)
    expected = d1 + d2
    actual = tofloat(bf1 + bf2)
    @test expected ≈ actual

    d1 = .1
    d2 = 1 - d1
    bf1 = binaryfraction(d1)
    bf2 = binaryfraction(d2)
    expected = d1 + d2
    actual = tofloat(bf1 + bf2)
    @test expected ≈ actual

    d1 = .1
    d2 = .2
    bf1 = binaryfraction(d1)
    bf2 = binaryfraction(d2)
    expected = d2 - d1
    actual = tofloat(bf2 - bf1)
    @test expected ≈ actual

    expected = d1 < d2
    actual = bf1 < bf2
    @test actual == expected
end

@testset "SymbolDistribution" begin
    symbols = ['a', 'b', 'c']
    p = [.1, .2, .7]
    sd = SymbolDistribution(symbols, p)
    expected = [(0.0, .1), (.1, .3), (.3, 1.0)]
    actual = sd.intervals
    @test all([all(expected[i] .≈ actual[i]) for i=1:length(expected)])

    expected = .1
    actual = getp(sd, 'a')
    @test expected == actual

    expected = (0.1, 0.3)
    actual = getinterval(sd, 'b')
    @test all(expected .≈ actual)
end

@testset "encode" begin
    s = "bbb□"
    sd = SymbolDistribution(['a', 'b', '□'], [.4, .5, .1])
    expected = BitVector("1101000" * repeat("0", nbits - 7))
    actual = encode(s, sd)
    @test expected == actual
end

max_stars_repo_name: chipfranzen/compression

max_issues_repo_name: chipfranzen/compression

max_forks_repo_name: chipfranzen/compression


text:
using Test
using Quaternions

test_associative(x, y, z, ⊗) = @test (x ⊗ y) ⊗ z ≈ x ⊗ (y ⊗ z)
test_commutative(x, y, ⊗) = @test x ⊗ y ≈ y ⊗ x
test_inverse(x, eins, ⊗, inv) = (@test x ⊗ inv(x) ≈ eins; @test inv(x) ⊗ x ≈ eins)
test_neutral(x, eins, ⊗) = (@test x ⊗ eins ≈ x; @test eins ⊗ x ≈ x)
test_monoid(x, y, z, ⊗, eins) = (test_associative(x, y, z, ⊗); test_neutral(x, eins, ⊗))
test_group(x, y, z, ⊗, eins, inv) = (test_monoid(x, y, z, ⊗, eins);test_inverse(x, eins, ⊗, inv))
test_multiplicative(x, y, ⊗, f) = @test f(x ⊗ y) ≈ f(x) ⊗ f(y)

include("test_Quaternion.jl")


text:
# ------------------------------------
# advection term by AUSM
# ------------------------------------
function AUSM(E_adv_hat, F_adv_hat, G_adv_hat, QbaseF, QbaseB, QbaseU, QbaseD, QbaseL, QbaseR,
                QconF, QconB, QconU, QconD, QconL, QconR, cellxmax, cellymax, cellzmax, vecAx, vecAy, vecAz,
                specific_heat_ratio, volume, nval, Minf, ad_scheme, icell)
    g         = specific_heat_ratio
    temp_vec = zeros(nval)   # vecAx used for the pressure term
    Lpsi = zeros(nval)
    Rpsi = zeros(nval)

    mdot = 0.0
    ph = 0.0

    tvA = zeros(3)
    vA  = zeros(3)

    for k in 1+icell:cellzmax -icell
        for j in 1+icell:cellymax -icell
            for i in 1+icell:cellxmax+1 -icell
                volume_av = 0.5*(volume[i-1,j,k] + volume[i,j,k])

                # i-1 cell
                # 規格化
                tvA[1] = 0.5*(vecAx[i-1,j,k,1] + vecAx[i,j,k,1])
                tvA[2] = 0.5*(vecAx[i-1,j,k,2] + vecAx[i,j,k,2])
                tvA[3] = 0.5*(vecAx[i-1,j,k,3] + vecAx[i,j,k,3])
                vA[1] = tvA[1] / (tvA[1]^2 + tvA[2]^2 + tvA[3]^2)^0.5
                vA[2] = tvA[2] / (tvA[1]^2 + tvA[2]^2 + tvA[3]^2)^0.5
                vA[3] = tvA[3] / (tvA[1]^2 + tvA[2]^2 + tvA[3]^2)^0.5

                rhoL = QbaseL[i,j,k,1]
                UL = QbaseL[i,j,k,2]*vA[1] + QbaseL[i,j,k,3]*vA[2] + QbaseL[i,j,k,4]*vA[3]
                pL = QbaseL[i,j,k,5]

                # i cell
                # 規格化
                tvA[1] = 0.5*(vecAx[i,j,k,1] + vecAx[i+1,j,k,1])
                tvA[2] = 0.5*(vecAx[i,j,k,2] + vecAx[i+1,j,k,2])
                tvA[3] = 0.5*(vecAx[i,j,k,3] + vecAx[i+1,j,k,3])
                vA[1] = tvA[1] / (tvA[1]^2 + tvA[2]^2 + tvA[3]^2)^0.5
                vA[2] = tvA[2] / (tvA[1]^2 + tvA[2]^2 + tvA[3]^2)^0.5
                vA[3] = tvA[3] / (tvA[1]^2 + tvA[2]^2 + tvA[3]^2)^0.5

                rhoR = QbaseR[i,j,k,1]
                UR = QbaseR[i,j,k,2]*vA[1] + QbaseR[i,j,k,3]*vA[2] + QbaseR[i,j,k,4]*vA[3]
                pR = QbaseR[i,j,k,5]

                # scheme
                if  ad_scheme == 1
                    mdot, ph = AUSM_plus_half(rhoL, rhoR, UL, UR, pL, pR, g, i, j)
                elseif ad_scheme == 2
                    mdot, ph = AUSM_plusup_half(rhoL, rhoR, UL, UR, pL, pR, Minf, g, i, j)
                elseif ad_scheme == 4
                    velocity = (0.5 * (QbaseL[i,j,k,2]^2 + QbaseL[i,j,k,3]^2 + QbaseL[i,j,k,4]^2
                                        + QbaseR[i,j,k,2]^2 + QbaseR[i,j,k,3]^2 + QbaseR[i,j,k,4]^2))^0.5
                    mdot, ph = SLAU_half(rhoL, rhoR, UL, UR, pL, pR, velocity, g)
                end

                # flux half
                sqAx = (vecAx[i,j,k,1]^2 + vecAx[i,j,k,2]^2 + vecAx[i,j,k,3]^2)^0.5
                temp_vec[2] = vecAx[i,j,k,1] / sqAx
                temp_vec[3] = vecAx[i,j,k,2] / sqAx
                temp_vec[4] = vecAx[i,j,k,3] / sqAx

                for l in 1:nval
                    Lpsi[l] = QconL[i,j,k,l] / QconL[i,j,k,1]
                    Rpsi[l] = QconR[i,j,k,l] / QconR[i,j,k,1]
                end
                Lpsi[5] = (QconL[i,j,k,5] + QbaseL[i,j,k,5]) / QconL[i,j,k,1]
                Rpsi[5] = (QconR[i,j,k,5] + QbaseR[i,j,k,5]) / QconR[i,j,k,1]

                if mdot > 0
                    for l in 1:nval
                        E_adv_hat[i,j,k,l] = (mdot * Lpsi[l] + ph * temp_vec[l]) * sqAx
                    end
                else
                    for l in 1:nval
                        E_adv_hat[i,j,k,l] = (mdot * Rpsi[l] + ph * temp_vec[l]) * sqAx
                    end
                end
                #=
                if i==10 && j ==10 && k==5
                    println(" E ")
                    println(mdot)
                    println(rhoL)
                    println(UL)
                    println(pL)
                    println(rhoR)
                    println(UR)
                    println(pR)
                end
                =#
            end
        end
    end

    for k in 1+icell:cellzmax -icell
        for j in 1+icell:cellymax+1 -icell
            for i in 1+icell:cellxmax -icell
                volume_av = 0.5*(volume[i,j-1,k] + volume[i,j,k])

                # j-1 cell
                # 規格化
                tvA[1] = 0.5*(vecAy[i,j-1,k,1] + vecAy[i,j,k,1])
                tvA[2] = 0.5*(vecAy[i,j-1,k,2] + vecAy[i,j,k,2])
                tvA[3] = 0.5*(vecAy[i,j-1,k,3] + vecAy[i,j,k,3])
                vA[1] = tvA[1] / (tvA[1]^2 + tvA[2]^2 + tvA[3]^2)^0.5
                vA[2] = tvA[2] / (tvA[1]^2 + tvA[2]^2 + tvA[3]^2)^0.5
                vA[3] = tvA[3] / (tvA[1]^2 + tvA[2]^2 + tvA[3]^2)^0.5

                rhoD = QbaseD[i,j,k,1]
                VD = QbaseD[i,j,k,2]*vA[1] + QbaseD[i,j,k,3]*vA[2] + QbaseD[i,j,k,4]*vA[3]
                pD = QbaseD[i,j,k,5]

                # j cell
                tvA[1] = 0.5*(vecAy[i,j,k,1] + vecAy[i,j+1,k,1])
                tvA[2] = 0.5*(vecAy[i,j,k,2] + vecAy[i,j+1,k,2])
                tvA[3] = 0.5*(vecAy[i,j,k,3] + vecAy[i,j+1,k,3])
                vA[1] = tvA[1] / (tvA[1]^2 + tvA[2]^2 + tvA[3]^2)^0.5
                vA[2] = tvA[2] / (tvA[1]^2 + tvA[2]^2 + tvA[3]^2)^0.5
                vA[3] = tvA[3] / (tvA[1]^2 + tvA[2]^2 + tvA[3]^2)^0.5

                rhoU = QbaseU[i,j,k,1]
                VU = QbaseU[i,j,k,2]*vA[1] + QbaseU[i,j,k,3]*vA[2] + QbaseU[i,j,k,4]*vA[3]
                pU = QbaseU[i,j,k,5]

                # scheme
                if  ad_scheme == 1
                    mdot, ph = AUSM_plus_half(rhoD, rhoU, VD, VU, pD, pU, g, i, j)
                elseif ad_scheme == 2
                    mdot, ph = AUSM_plusup_half(rhoD, rhoU, VD, VU, pD, pU, Minf, g, i, j)
                elseif ad_scheme == 4
                    velocity = (0.5 * (QbaseD[i,j,k,2]^2 + QbaseD[i,j,k,3]^2 + QbaseD[i,j,k,4]^2 +
                                    QbaseU[i,j,k,2]^2 + QbaseU[i,j,k,3]^2 + QbaseU[i,j,k,4]^2))^0.5
                    mdot, ph = SLAU_half(rhoD, rhoU, VD, VU, pD, pU, velocity, specific_heat_ratio)
                end

                # flux half
                sqAy = (vecAy[i,j,k,1]^2 + vecAy[i,j,k,2]^2 + vecAy[i,j,k,3]^2)^0.5
                temp_vec[2] = vecAy[i,j,k,1] / sqAy
                temp_vec[3] = vecAy[i,j,k,2] / sqAy
                temp_vec[4] = vecAy[i,j,k,3] / sqAy

                for l in 1:nval
                    Lpsi[l] = QconD[i,j,k,l] / QconD[i,j,k,1]
                    Rpsi[l] = QconU[i,j,k,l] / QconU[i,j,k,1]
                end
                Lpsi[5] = (QconD[i,j,k,5] + QbaseD[i,j,k,5]) / QconD[i,j,k,1]
                Rpsi[5] = (QconU[i,j,k,5] + QbaseU[i,j,k,5]) / QconU[i,j,k,1]

                if mdot > 0
                    for l in 1:nval
                        F_adv_hat[i,j,k,l] = (mdot * Lpsi[l] + ph * temp_vec[l]) * sqAy
                    end
                else
                    for l in 1:nval
                        F_adv_hat[i,j,k,l] = (mdot * Rpsi[l] + ph * temp_vec[l]) * sqAy
                    end
                end
                #=
                if i==10 && j ==10 && k==5
                    println(" F ")
                    println(mdot)
                    println(rhoD)
                    println(VD)
                    println(pD)
                    println(rhoU)
                    println(VU)
                    println(pU)
                end
                =#
                #=
                if isequal(vA[3], NaN) == true
                    println("  ")
                    println(i)
                    println(j)
                    println(k)
                    println(tvA)
                    println(vecAy[i,j,k,:])
                    println(vecAy[i,j+1,k,:])
                end
                =#
            end
        end
    end

    for k in 1+icell:cellzmax+1 -icell
        for j in 1+icell:cellymax -icell
            for i in 1+icell:cellxmax -icell
                volume_av = 0.5*(volume[i,j,k-1] + volume[i,j,k])

                # k-1 cell
                # 規格化
                tvA[1] = 0.5*(vecAz[i,j,k-1,1] + vecAz[i,j,k,1])
                tvA[2] = 0.5*(vecAz[i,j,k-1,2] + vecAz[i,j,k,2])
                tvA[3] = 0.5*(vecAz[i,j,k-1,3] + vecAz[i,j,k,3])
                vA[1] = tvA[1] / (tvA[1]^2 + tvA[2]^2 + tvA[3]^2)^0.5
                vA[2] = tvA[2] / (tvA[1]^2 + tvA[2]^2 + tvA[3]^2)^0.5
                vA[3] = tvA[3] / (tvA[1]^2 + tvA[2]^2 + tvA[3]^2)^0.5

                rhoB = QbaseB[i,j,k,1]
                UB = QbaseB[i,j,k,2]*vA[1] + QbaseB[i,j,k,3]*vA[2] + QbaseB[i,j,k,4]*vA[3]
                pB = QbaseB[i,j,k,5]

                # k cell
                # 規格化
                tvA[1] = 0.5*(vecAz[i,j,k,1] + vecAz[i,j,k+1,1])
                tvA[2] = 0.5*(vecAz[i,j,k,2] + vecAz[i,j,k+1,2])
                tvA[3] = 0.5*(vecAz[i,j,k,3] + vecAz[i,j,k+1,3])
                vA[1] = tvA[1] / (tvA[1]^2 + tvA[2]^2 + tvA[3]^2)^0.5
                vA[2] = tvA[2] / (tvA[1]^2 + tvA[2]^2 + tvA[3]^2)^0.5
                vA[3] = tvA[3] / (tvA[1]^2 + tvA[2]^2 + tvA[3]^2)^0.5

                rhoF = QbaseF[i,j,k,1]
                UF = QbaseF[i,j,k,2]*vA[1] + QbaseF[i,j,k,3]*vA[2] + QbaseF[i,j,k,4]*vA[3]
                pF = QbaseF[i,j,k,5]

                # scheme
                if  ad_scheme == 1
                    mdot, ph = AUSM_plus_half(rhoB, rhoF, UB, UF, pB, pF, g, i, j)
                elseif ad_scheme == 2
                    mdot, ph = AUSM_plusup_half(rhoB, rhoF, UB, UF, pB, pF, Minf, g, i, j)
                elseif ad_scheme == 4
                    velocity = (0.5 * (QbaseF[i,j,k,2]^2 + QbaseF[i,j,k,3]^2 + QbaseF[i,j,k,4]^2
                                        + QbaseB[i,j,k,2]^2 + QbaseB[i,j,k,3]^2 + QbaseB[i,j,k,4]^2))^0.5
                    mdot, ph = SLAU_half(rhoB, rhoF, UB, UF, pB, pF, velocity, g)
                end

                # flux half
                sqAz = (vecAz[i,j,k,1]^2 + vecAz[i,j,k,2]^2 + vecAz[i,j,k,3]^2)^0.5
                temp_vec[2] = vecAz[i,j,k,1] / sqAz
                temp_vec[3] = vecAz[i,j,k,2] / sqAz
                temp_vec[4] = vecAz[i,j,k,3] / sqAz

                for l in 1:nval
                    Lpsi[l] = QconB[i,j,k,l] / QconB[i,j,k,1]
                    Rpsi[l] = QconF[i,j,k,l] / QconF[i,j,k,1]
                end
                Lpsi[5] = (QconB[i,j,k,5] + QbaseB[i,j,k,5]) / QconB[i,j,k,1]
                Rpsi[5] = (QconF[i,j,k,5] + QbaseF[i,j,k,5]) / QconF[i,j,k,1]

                if mdot > 0
                    for l in 1:nval
                        G_adv_hat[i,j,k,l] = (mdot * Lpsi[l] + ph * temp_vec[l]) * sqAz
                    end
                else
                    for l in 1:nval
                        G_adv_hat[i,j,k,l] = (mdot * Rpsi[l] + ph * temp_vec[l]) * sqAz
                    end
                end
                #=
                if i==10 && j ==10 && k==5
                    println("  ")
                    println(mdot)
                end
                =#
            end
        end
    end

    return E_adv_hat, F_adv_hat, G_adv_hat
end

function AUSM_plus_half(rhoL, rhoR, UL, UR, pL, pR, g,i,j)
    # param
    beta  = 1/8
    alpha = 3/16

    # L, R
    aL = (g * pL / rhoL)^0.5
    aR = (g * pR / rhoR)^0.5

    # half
    ah   = 0.5*( aL + aR )
    rhoh = 0.5*( rhoL + rhoR )

    ML = UL/ah
    MR = UR/ah

    M_p4 = 0
    M_m4 = 0
    p_p5 = 0
    p_m5 = 0
    if abs(ML) >= 1
        M_p4 = 0.5*(ML + abs(ML))
        p_p5 = 0.5*(ML + abs(ML)) / ML
    else
        M_p4 = 0.25*(ML+1)^2 + beta*(ML^2-1)^2
        p_p5 = 0.25*(ML+1)^2 * (2-ML) + alpha*ML*(ML^2-1)^2
    end

    if abs(MR) >= 1
        M_m4 = 0.5*(MR - abs(MR))
        p_m5 = 0.5*(MR - abs(MR)) / MR
    else
        M_m4 = -0.25*(MR-1)^2 - beta*(MR^2-1)^2
        p_m5 = 0.25*(MR-1)^2 * (2+MR) - alpha*MR*(MR^2-1)^2
    end

    # Mh = M_p4 + M_m4 + Mp/fa
    Mh = M_p4 + M_m4
#=
    if i==100 && j==200
        println(" MMM_AUSM+ ")
        println(M_p4)
        println(M_m4)
    end
  =#
    # mdot half
    mdot = ah * Mh
    if Mh > 0
        mdot = mdot * rhoL
    else
        mdot = mdot * rhoR
    end

    # p half
    ph = p_p5*pL + p_m5*pR

    return mdot, ph
end

function AUSM_plusup_half(rhoL, rhoR, UL, UR, pL, pR, Minf, g, i,j)
    # param
    beta  = 1/8
    kp    = 0.25
    sigma = 1.0
    ku    = 0.75

    # L, R
    aL = (g * pL / rhoL)^0.5
    aR = (g * pR / rhoR)^0.5

    # half
    ah   = 0.5*( aL + aR )
    rhoh = 0.5*( rhoL + rhoR )

    Mbar = (( UL^2 + UR^2 ) / ( 2 * ah^2 ))^0.5
    Mo   = (min(1,max( Mbar^2, Minf^2 )))^0.5
    fa   = Mo * (2-Mo)

    alpha = 3/16 * ( -4 + 5*fa^2 )

    ML = UL/ah
    MR = UR/ah

    M_p4 = 0
    M_m4 = 0
    p_p5 = 0
    p_m5 = 0
    if abs(ML) >= 1
        M_p4 = 0.5*(ML + abs(ML))
        p_p5 = 0.5*(ML + abs(ML)) / ML
    else
        Mtp  = 0.25*(ML + 1)^2
        Mtm  = -0.25*(ML - 1)^2
        M_p4 = Mtp * (1 - 16*beta*Mtm)
        p_p5 = Mtp * ((2-ML) - 16*alpha*ML*Mtm)
    end

    if abs(MR) >= 1
        M_m4 = 0.5*(MR - abs(MR))
        p_m5 = 0.5*(MR - abs(MR)) / MR
    else
        Mtp  = 0.25*(MR + 1)^2
        Mtm  = -0.25*(MR - 1)^2
        M_m4 = Mtm * (1 + 16*beta*Mtp)
        p_m5 = Mtm * ((-2-MR) + 16*alpha*MR*Mtp)
    end

    # M half
    Mp = -kp * max( 1 - sigma*Mbar^2, 0.0) * (pR - pL)/(rhoh*ah^2)

    Mh = M_p4 + M_m4 + Mp/fa

    #=
    if i==100 && j==200
        println(" MMM_AUSM+up ")
        println(M_p4)
        println(M_m4)
        println(Mp/fa)
    end
    =#

    # mdot half
    mdot = ah * Mh
    if Mh > 0
        mdot = mdot * rhoL
    else
        mdot = mdot * rhoR
    end

    # p half
    pu = -ku * p_p5 * p_m5 * (rhoL + rhoR) * ah *(UR - UL)

    ph = p_p5*pL + p_m5*pR + fa * pu
    return mdot, ph
end

function set_Minf(bdcon, specific_heat_ratio, Rd, nval)
    Minf = 0.0
    check = 0
    for i in 1:6
        if Int(bdcon[i][1]) == 0 || Int(bdcon[i][1]) == 5
            rho = bdcon[i][2]
            u = bdcon[i][3]
            v = bdcon[i][4]
            w = bdcon[i][5]
            p = bdcon[i][6]

            a = (specific_heat_ratio * p / rho)^0.5
            Minf = (u^2 + v^2 + w^2)^0.5 / a

            println(" Minf = " * string(Minf))

            check = 1

        elseif Int(bdcon[i][1]) == 6 || Int(bdcon[i][1]) == 88
            rho = bdcon[i][2]
            u = bdcon[i][3]
            v = bdcon[i][4]
            w = bdcon[i][5]
            T = bdcon[i][nval+2]

            p = (rho*Rd) * T

            a = (specific_heat_ratio * p / rho)^0.5
            Minf = (u^2 + v^2 + w^2)^0.5 / a

            println(" Minf = " * string(Minf))

            check = 1
        end
    end

    if check == 0
        println(" Minf error ")
        println("  ")
        println(" if you don't use inlet condition, ")
        println(" don't use AUSM+up ")
        println("  ")
        throw(UndefVarError(:x))
    end

    return Minf
end

function SLAU_half(rhoL, rhoR, UL, UR, pL, pR, velocity, specific_heat_ratio)
    mdot =1.0
    ph =1.0

    # L, R
    aL = (specific_heat_ratio * pL / rhoL)^0.5
    aR = (specific_heat_ratio * pR / rhoR)^0.5

    # half
    ah   = 0.5*( aL + aR )

    ML = UL/ah
    MR = UR/ah

    temp1 = -max( min( ML, 0.0 ), -1.0)
    temp2 =  min( max( MR, 0.0 ), 1.0)
    g = temp1 * temp2

    barUL = (1-g) * (rhoL*abs(UL) + rhoR*abs(UR))/(rhoL + rhoR) + g*abs(UL)
    barUR = (1-g) * (rhoL*abs(UL) + rhoR*abs(UR))/(rhoL + rhoR) + g*abs(UR)

    hatM = min(1.0, 1/ah * velocity)
    chi  = (1-hatM)^2

    beta_a = 0
    beta_b = 0
    if abs(MR) >= 1
        beta_a = 0.5 * (1+sign(-MR))
    else
        beta_a = 0.25 * (2+MR) * (MR-1)^2
    end
    if abs(ML) >= 1
        beta_b = 0.5 * (1+sign(ML))
    else
        beta_b = 0.25 * (2-ML) * (ML+1)^2
    end

    # mdot half
    mdot = 0.5 * (rhoL*(UL+abs(barUL)) + rhoR*(UR-abs(barUR)) - chi/ah*(pR-pL))

    # p half
    ph = 0.5*(pR+pL) + 0.5*(beta_b-beta_a)*(pL-pR) + (1-chi) * (beta_b+beta_a-1) * 0.5*(pR+pL)

    # SLAU2
    #ph = 0.5*(pR+pL) + 0.5*(beta_b-beta_a) + ((UR^2+UL^2)/2)^0.5 * (beta_b+beta_a-1) * 0.5*(rhoR+rhoL)*ah

    return mdot, ph
end


text:
module GraphHierarchy

using LightGraphs
using SimpleWeightedGraphs
using LinearAlgebra
using SparseArrays
using IterativeSolvers
using StatsBase
using Statistics

export forward_hierarchical_structure
export backward_hierarchical_structure
export hierarchical_structure
export forward_hierarchical_coefficients
export backward_hierarchical_coefficients
export hierarchical_coefficients

function weightedMeanStd(v::Array{T,1} where T <: AbstractFloat, w::Array{T,1} where T <: AbstractFloat)
    difftype = typeof(v[1])
    if difftype == typeof(w[1])
        av0 = mean(v, StatsBase.weights(w))
        return ( 1 - av0 ,std( v, StatsBase.weights(w), mean=av0, corrected=false))
    else
        ww = [ convert( difftype, c) for c in w ]
        av = mean(v, StatsBase.weights(ww))
        return ( 1 - av ,std( v, StatsBase.weights(ww), mean=av, corrected=false))
    end
end

function calculatingStructureAssumingTransposedMatrix(A::SparseMatrixCSC{T1,T2} where T1 <: Number where T2 <: Int)
    d = vec(sum(A, dims=2))
    dout = vec(sum(A, dims=1))
    M = spdiagm(0 => d)-A
    dropzeros!(M)
    HLs = lsqr(M,d)
    pos = findall(!iszero,A)
    HDs = sparse([p[2] for p in pos], [p[1] for p in pos], [HLs[p[1]] - HLs[p[2]] for p in pos], length(d), length(d))
    HDsums = sum(HDs, dims=1)
    ICs = [d[i] > 0 ? 1.0 - sum(HDs[:,i])/d[i] : 1.0 for i in 1:length(d)]
    RCs = [dout[i] > 0 ? sum(HDs[i,:])/dout[i] : 0.0 for i in 1:length(dout)]
    return (HLs, ICs, RCs, HDs )
end

function calculatingStructureAssumingTransposedMatrixBig(A::SparseMatrixCSC{T1,T2} where T1 <: Number where T2 <: Int)
    (x, y, v) = findnz(A)
    vb = -map(BigFloat,v)
    M = sparse(x, y, vb, size(A)[1], size(A)[1])
    d = -vec(sum(M, dims=2))
    dout = -vec(sum(M, dims=1))
    for i in 1:length(d)
        if d[i] > 0
            M[i,i] = d[i]
        end
    end
    HLs = lsqr(M,d)
    HDs = sparse(y, x, [ HLs[x[i]] - HLs[y[i]] for i in 1:length(x) ], length(d), length(d) )
    ICs = [d[i] > 0 ? BigFloat(1) - sum(HDs[:,i])/d[i] : BigFloat(1) for i in 1:length(d)]
    RCs = [dout[i] > 0 ? sum(HDs[i,:])/dout[i] : BigFloat(0)  for i in 1:length(dout)]
    return (HLs, ICs, RCs, HDs )
end

function calculatingStructureAssumingTransposedWeightedMatrix(A::SparseMatrixCSC{T1,T2} where T1 <: Number where T2 <: Int)
    (x, y, v) = findnz(A)
    d = vec(sum(A, dims=2))
    dout = vec(sum(A, dims=1))
    M = spdiagm(0 => d)-A
    dropzeros!(M)
    HLs = lsqr(M,d)
    HDs = sparse(y, x, [ HLs[x[i]] - HLs[y[i]] for i in 1:length(x) ], length(d), length(d) )
    HDsums = sum(HDs, dims=1)
    ICs = [d[i] > 0 ? 1.0 - A[i,:]'*HDs[:,i]/d[i] : 1.0 for i in 1:length(d)]
    RCs = [dout[i] > 0 ? A[:,i]'*HDs[i,:]/dout[i] : 0.0 for i in 1:length(dout)]
    return (HLs, ICs, RCs, HDs )
end

function calculatingStructureAssumingTransposedWeightedMatrixBig(A::SparseMatrixCSC{T1,T2} where T1 <: Number where T2 <: Int)
    (x, y, v) = findnz(A)
    vb = -map(BigFloat,v)
    M = sparse(x, y, vb, size(A)[1], size(A)[1])
    d = -vec(sum(M, dims=2))
    dout = -vec(sum(M, dims=1))
    for i in 1:length(d)
        if d[i] > 0
            M[i,i] = d[i]
        end
    end
    HLs = lsqr(M,d)
    HDs = sparse(y, x, [ HLs[x[i]] - HLs[y[i]] for i in 1:length(x) ], length(d), length(d) )
    HDsums = sum(HDs, dims=1)
    ICs = [d[i] > 0 ? BigFloat(1) - A[i,:]'*HDs[:,i]/d[i] : BigFloat(1) for i in 1:length(d)]
    RCs = [dout[i] > 0 ? A[:,i]'*HDs[i,:]/dout[i] : BigFloat(0) for i in 1:length(dout)]
    return (HLs, ICs, RCs, HDs )
end

function forward_hierarchical_structure(g::SimpleDiGraph{T} where T <: Int; big::Bool=false)
    if big
        return calculatingStructureAssumingTransposedMatrixBig(sparse(transpose(adjacency_matrix(g))))
    else
        return calculatingStructureAssumingTransposedMatrix(sparse(transpose(adjacency_matrix(g))))
    end
end

function backward_hierarchical_structure(g::SimpleDiGraph{T} where T <: Int; big::Bool=false)
    if big
        (xb, yb, wb, zb) = calculatingStructureAssumingTransposedMatrixBig(adjacency_matrix(g))
        return (xb, yb, wb, sparse(transpose(zb)))
    else
        (x, y, w, z) = calculatingStructureAssumingTransposedMatrix(adjacency_matrix(g))
        return (x, y, w, sparse(transpose(z)))
    end
end

function hierarchical_structure(g::SimpleDiGraph{T} where T <: Int; big::Bool=false)
    return (forward_hierarchical_structure(g, big = big),backward_hierarchical_structure(g, big = big))
end

function forward_hierarchical_structure(g::SimpleWeightedDiGraph{T1,T2} where T1 <: Int where T2 <: Number; big::Bool=false, weighted::Bool=true)
    if weighted
        if big
            return calculatingStructureAssumingTransposedWeightedMatrixBig(transpose(LightGraphs.weights(g)))
        else
            return calculatingStructureAssumingTransposedWeightedMatrix(transpose(LightGraphs.weights(g)))
        end
    else
        if big
            return calculatingStructureAssumingTransposedMatrixBig(sparse(transpose((!iszero).(adjacency_matrix(g)))))
        else
            return calculatingStructureAssumingTransposedMatrix(sparse(transpose((!iszero).(adjacency_matrix(g)))))
        end
    end
end

function backward_hierarchical_structure(g::SimpleWeightedDiGraph{T1,T2} where T1 <: Int where T2 <: Number; big::Bool=false, weighted::Bool=true)
    if weighted
        if big
            (xb, yb, wb, zb) = calculatingStructureAssumingTransposedWeightedMatrixBig(sparse(LightGraphs.weights(g)))
            return (xb, yb, wb, sparse(transpose(zb)))
        else
            (x, y, w, z) = calculatingStructureAssumingTransposedWeightedMatrix(sparse(LightGraphs.weights(g)))
            return (x, y, w, sparse(transpose(z)))
        end
    else
        if big
            (xb, yb, wb, zb) = calculatingStructureAssumingTransposedMatrixBig(sparse((!iszero).(adjacency_matrix(g))))
            return (xb, yb, wb, sparse(transpose(zb)))
        else
            (x, y, w, z) = calculatingStructureAssumingTransposedMatrix(sparse((!iszero).(adjacency_matrix(g))))
            return (x, y, w, sparse(transpose(z)))
        end
    end
end

function hierarchical_structure(g::SimpleWeightedDiGraph{T1,T2} where T1 <: Int where T2 <: Number; big::Bool=false, weighted::Bool=true)
    return (forward_hierarchical_structure(g, big = big, weighted = weighted),backward_hierarchical_structure(g, big = big, weighted = weighted))
end

function forward_hierarchical_structure(g::SimpleGraph{T} where T <: Int; big::Bool=false)
    if big
        return calculatingStructureAssumingTransposedMatrixBig(adjacency_matrix(g))
    else
        return calculatingStructureAssumingTransposedMatrix(adjacency_matrix(g))
    end
end

function backward_hierarchical_structure(g::SimpleGraph{T} where T <: Int; big::Bool=false)
    forward_hierarchical_structure(g, big = big)
end

function hierarchical_structure(g::SimpleGraph{T} where T <: Int; big::Bool=false)
    return forward_hierarchical_structure(g, big = big)
end

function forward_hierarchical_structure(g::SimpleWeightedGraph{T1,T2} where T1 <: Int where T2 <: Number; big::Bool=false, weighted::Bool=true)
    if weighted
        if big
            return calculatingStructureAssumingTransposedWeightedMatrixBig(sparse(LightGraphs.weights(g)))
        else
            return calculatingStructureAssumingTransposedWeightedMatrix(sparse(LightGraphs.weights(g)))
        end
    else
        if big
            return calculatingStructureAssumingTransposedMatrixBig(adjacency_matrix(g))
        else
            return calculatingStructureAssumingTransposedMatrix(adjacency_matrix(g))
        end
    end
end

function backward_hierarchical_structure(g::SimpleWeightedGraph{T1,T2} where T1 <: Int where T2 <: Number; big::Bool=false, weighted::Bool=true)
    return forward_hierarchical_structure(g, big = big, weighted = weighted)
end

function hierarchical_structure(g::SimpleWeightedGraph{T1,T2} where T1 <: Int where T2 <: Number; big::Bool=false, weighted::Bool=true)
    return forward_hierarchical_structure(g, big = big, weighted = weighted)
end

function forward_hierarchical_structure(A::SparseMatrixCSC{T1,T2} where T1 <: Number where T2 <: Int; big::Bool=false)
    if big
        return calculatingStructureAssumingTransposedMatrixBig(sparse(transpose(A)))
    else
        return calculatingStructureAssumingTransposedMatrix(sparse(transpose(A)))
    end
end

function backward_hierarchical_structure(A::SparseMatrixCSC{T1,T2} where T1 <: Number where T2 <: Int; big::Bool=false)
    if big
        return calculatingStructureAssumingTransposedMatrixBig(A)
    else
        return calculatingStructureAssumingTransposedMatrix(A)
    end
end

function hierarchical_structure(A::SparseMatrixCSC{T1,T2} where T1 <: Number where T2 <: Int; big::Bool=false)
    return ( forward_hierarchical_structure(A, big=big), backward_hierarchical_structure(A, big=big) )
end

function hierarchical_coefficients(g::Union{SimpleGraph{T},SimpleDiGraph{T}} where T <: Int, diffs::SparseMatrixCSC{T1,T2} where T1 <: AbstractFloat where T2 <: Int)
    diffarray = findnz(diffs)[3]
    av = mean(diffarray)
    return (1 - av, std(diffarray, mean=av, corrected=false))
end

function hierarchical_coefficients(g::Union{SimpleGraph{T},SimpleDiGraph{T}} where T <: Int, hstruc::Tuple{Array{T1,1},Array{T1,1},Array{T1,1},SparseMatrixCSC{T1,T2}} where T1 <: AbstractFloat where T2 <: Int )
    return hierarchical_coefficients(g, hstruc[4])
end

function hierarchical_coefficients(g::SimpleDiGraph{T} where T <: Int , hstruc::Tuple{Tuple{Array{T1,1},Array{T1,1},Array{T1,1},SparseMatrixCSC{T1,T2}},Tuple{Array{T1,1},Array{T1,1},Array{T1,1},SparseMatrixCSC{T1,T2}}} where T1 <: AbstractFloat where T2 <: Int )
    return (hierarchical_coefficients(g, hstruc[1]),hierarchical_coefficients(g, hstruc[2]))
end

function hierarchical_coefficients(g::Union{SimpleWeightedGraph{T1,T2},SimpleWeightedDiGraph{T1,T2}} where T1 <: Int where T2 <: AbstractFloat, diffs::SparseMatrixCSC{T1,T2} where T1 <: AbstractFloat where T2 <: Int; weighted::Bool = true )
    if weighted
        W0 = sparse(LightGraphs.weights(g))
        (xw,yw,vw) = findnz(W0)
        (x,y,v) = findnz(diffs)

        if sum(map(abs,x - xw)) + sum(map(abs,y - yw)) > 0
            error("Hierarchical differences do not match the weighted adjacency matrix:\nDifferences:", diffs, "\nWeighted Adjacency Matrix:", W0)
        end

        return weightedMeanStd(v,vw)
    else
        v = findnz(diffs)[3]
        av = mean(v)
        return (1 - av, std(v, mean=av, corrected=false))
    end
end

function hierarchical_coefficients(g::Union{SimpleWeightedGraph{T1,T2},SimpleWeightedDiGraph{T1,T2}} where T1 <: Int  where T2 <: AbstractFloat, hstruc::Tuple{Array{T1,1},Array{T1,1},Array{T1,1},SparseMatrixCSC{T1,T2}} where T1 <: AbstractFloat where T2 <: Int;  weighted::Bool = true  )
    return hierarchical_coefficients(g, hstruc[4], weighted=weighted)
end

function hierarchical_coefficients(g::SimpleWeightedDiGraph{T1,T2} where T1 <: Int  where T2 <: AbstractFloat , hstruc::Tuple{Tuple{Array{T1,1},Array{T1,1},Array{T1,1},SparseMatrixCSC{T1,T2}},Tuple{Array{T1,1},Array{T1,1},Array{T1,1},SparseMatrixCSC{T1,T2}}} where T1 <: AbstractFloat where T2 <: Int; weighted::Bool = true  )
    return (hierarchical_coefficients(g, hstruc[1], weighted=weighted), hierarchical_coefficients(g, hstruc[2], weighted=weighted))
end

end # module


text:
using Distributions, Statistics

"""
    KernelScenarios(data::Union{Array{Int64,2},Array{Float64,2}},kernelDistribution = Logistic)

    This function takes a matrix of data in 2 dimension (NxT) where N is the number of trajectories and T is the number of stages.
    It uses function closures and above is a function.
    To get results a trajectory from above, `KernelScenarios(data,kernelDistribution)()` gives the result.
    User specifies both data and the distribution of the kernel he/she wants to use. The default kernel is the Logistic kernel.
"""

function KernelScenarios(data::Union{Array{Int64,2},Array{Float64,2}},kernelDistribution = Logistic; Markovian::Bool = true)
    # We use the function closure here because we want to use this function in the stochastic approximation process.
    # That is, we wwant to follow the specifications of the TreeApproximation and Latticeapproximation inputs.
    # The way we have created the stochastic approximation functions is that it takes a function which does not take any inputs.
    # And so the outer function takes the data and the distribution of the kernel while the inner function takes no inputs.
    # To obtain the result then, call "KernelScenarios(data,Logistic)()".
    # Notice we have an input of whether you are creating Markovian trajectories or not.
    # This is important for generating either scenario trees or scenario lattices.
    # The default for Markovian is true hence the trajectory created is Markovian which is used for scenario lattices.
    function closure()
        N,T = size(data)                                                  # Dimensions of the data
        d = 1
        w = fill(1.0,N)                                                     # initialize the weights
        ξ = Array{Float64,1}(undef,T)                                     # trajectory to be created
        for t = 1:T
            w = w / sum(w)                                                # normalized weights
            Nt = sum(w)^2 / sum(w.^2)                                     # effective sample size
            σt = sqrt(sum(w .* ((data[:,t] .- sum(w .* data[:,t])).^2)))  # effective standard devistion
            ht = σt * Nt^(-1/(d + 4)) + eps()
            # Composition method comes in here
            u = rand(Uniform(0,1))
            _,jstar = findmin(cumsum(w) .< u * sum(w))
            # The cumulative sum of weights leads to a high probability of picking a data path near an observation (Sequin && Pichler 2017)
            ξ[t] = rand(kernelDistribution(data[jstar,t],ht))
            # The new generated value is according to the distribution of density of the current stage and dependent on the history of all the data paths.
            # All the rows in the column are associated with a certain weight as follows:
            # Each weight is a product of the kernels of the data at that point.
            if t<T
                for j = 1:N
                    # The choice of the kernel does not have any important effect on density estimation (Jones (1990)).
                    if Markovian
                        # Markovian is used for scenario lattices
                        w[j] = pdf(kernelDistribution(data[j,t],ht),ξ[t])
                    else
                        # Non-Markovian is used for scenario lattices
                        w[j] *= pdf(kernelDistribution(data[j,t],ht),ξ[t])
                    end
                end
            end
        end
        return ξ                         # The length of ξ is equal to the number stages = T
    end
end


text:
using HCIToolbox: get_annulus_segments

"""
    Median()

Classic PSF subtraction using the median of entire data cube.

# References
* [Marois et al. 2006](http://adsabs.harvard.edu/abs/2006ApJ...641..556M) Angular Differential Imaging: A Powerful High-Contrast Imaging Technique
"""
struct Median <: ADIAlgorithm end

reconstruct(::Median, cube, angles, cube_ref=cube) = median(cube_ref, dims=1)

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

# # Optimized median subtraction for a given annulus
# function medsub_annular(cube::AbstractArray{T,3}, angles::AbstractVector; fwhm, ann_size=4, radius_int=0, nframes=4, delta_sep=(0.1, 1), delta_rot=1, kwargs...) where T
#     angles_ = normalize_par_angles(angles)
#     med = median(cube, dims=1)
#     S = cube .- med
#     n, y, x = size(cube)

#     n_annuli = round(Int, (y / 2 - radius_int) / ann_size)

#     # _median_subt_ann_adi
#     ann = 0:n_annuli - 1
#     pa_thresh, inner_r, _ = define_annuli(angles, ann, n_annuli, fwhm, radius_int, ann_size, delta_rot)
#     mask = get_annulus_segments(@view(S[1, :, :]), first(inner_r), ann_size)
#     indices = findall(mask)
#     matrix = S[:, indices]
#     out = similar(matrix)

#     for frame in 1:n
#         if !iszero(pa_thresh)
#             inds_left = find_indices_adi.((angles,), frame, pa_thresh, nframes)
#             @show size(pa_thresh)
#             @show size(inds_left)
#             matrix_disc = @view matrix[inds_left, :, :]
#         else
#             matrix_disc = matrix
#         end
#         ref_psf_opt = median(matrix_disc, dims=1)
#         curr_frame = @view matrix[frame, :, :]
#         @. out[frame, :, :] = curr_frame - ref_psf_opt
#     end
#     # return out, yy, xx, pa_thresh
#     # end _median_subt_ann_adi

#     for  ann in 1:n_annuli
#         out[:, yy[ann], xx[ann]] .= out[ann, :, :]
#     end

# end

# function find_indices_adi(angles, frame, pa_thresh, nframes)
#     n = length(angles)

#     index_prev = findfirst(1:frame - 1) do i
#         abs(angles[frame] - angles[i]) < pa_thresh
#     end
#     isnothing(index_prev) && (index_prev = frame - 1)

#     index_foll = findfirst(frame:n) do k
#         abs(angles[k] - angles[frame]) > pa_thresh
#     end
#     isnothing(index_foll) && (index_foll = n)

#     # nframesi s not None (median annular)
#     window = nframes ÷ 2

#     ind1 = max(index_prev - window, 1)
#     ind2 = index_prev - 1
#     ind3 = index_foll
#     ind4 = min(index_foll + window, n) - 1
#     indices = vcat(ind1:ind2, ind3:ind4)
# end

# """
# Function that defines the annuli geometry using the input parameters.
# Returns the parallactic angle threshold, the inner radius and the annulus
# center for each annulus.
# """
# function define_annuli(angles, ann, N, fwhm, radius_int, width, delta_rot)
#     inner_rad = @. ifelse(ann == N - 1, radius_int + (ann * width - 1), radius_int + ann * width)
#     ann_center = @. inner_rad + (width / 2)
#     thresh = @. rad2deg(2 * atan(delta_rot * fwhm / (2 * ann_center)))
#     mid_range = abs(maximum(angles) - minimum(angles)) / 2
#     @. thresh[thresh ≥ 0.9 * mid_range] = 0.9 * mid_range
#     return thresh, inner_rad, ann_center
# end