text:
importall LegHemodynamics

function main()

rstflag = "yes" # restarting from scratch or previous simulation
hemoflag = "no" # 10% hemorrhage from left femoral artery
saveflag = "yes" # save solution to .mat file
coupleflag = "no" # coupling to 3D liver tissue model
timeflag = "yes" # solver timing
assimflag = "yes" # patient data assimilation via EnKF

# build solution struct or generate ensemble
if assimflag == "no"
    if rstflag == "no"
        loadfile = "arterylist_2.txt"; # default artery data file for new sim
    elseif rstflag == "yes"
        loadfile = "test.mat"; # restart file
    end
    system = LegHemodynamics.buildall(loadfile;numbeatstotal=1,restart=rstflag);
    savefile = "test.mat" # filename for saving (only used if saveflag == "yes")
elseif assimflag == "yes"
    ensemblesize = 10;
    if rstflag == "no"
        loadfiles = ["arterylist_2_thick_branches.txt" for i=1:ensemblesize];
    elseif rstflag == "yes"
        loadfiles = ["small_dev_7_$i.mat" for i=1:ensemblesize];
    end
    systems = pmap((a1)->LegHemodynamics.buildall(a1;numbeatstotal=1,restart=rstflag),loadfiles);
    savefiles = ["converged_small_dev_1_$i.mat" for i=1:ensemblesize];
    if rstflag == "no" # randomize ICs on new ensemble
        soln = pmap((a1)->LegHemodynamics.applycustomics!(a1),systems);
        systems = [soln[i] for i=1:ensemblesize];
    elseif rstflag == "yes" # ensure uniform time step
        h = [];
        for i = 1:length(systems)
            push!(h,systems[i].solverparams.h)
        end
        for i = 1:length(systems)
            systems[i].solverparams.h = minimum(h);
        end
        systems = pmap((a1,a2)->LegHemodynamics.rediscretizet!(a1,a2),systems,zeros(Int64,length(systems)));
    end
end

# create interpolation object for upstream flow velocity data (needed for proximal BC)
patientdata = MAT.matread("flowrates_mL.mat");
tdata = reshape(patientdata["t"],101);
udata = reshape(patientdata["u_pop_cms"],101)*LegHemodynamics.cmTom;
if assimflag == "no"
    spl = Dierckx.Spline1D(tdata,udata;k=1);
elseif assimflag == "yes"
    spl = [Dierckx.Spline1D(tdata,udata;k=1) for i=1:ensemblesize];
    udata = reshape(patientdata["u_ant_cms"],101)*LegHemodynamics.cmTom;
    splant = Dierckx.Spline1D(tdata,udata;k=1);
    udata = reshape(patientdata["u_post_cms"],101)*LegHemodynamics.cmTom;
    splpost = Dierckx.Spline1D(tdata,udata;k=1);
    udata = reshape(patientdata["u_per_cms"],101)*LegHemodynamics.cmTom;
    splper = Dierckx.Spline1D(tdata,udata;k=1);
end

# timers
if assimflag == "no"
    times = LegHemodynamics.CVTimer();
elseif assimflag == "yes"
    times = [LegHemodynamics.CVTimer() for i=1:ensemblesize];
end

# collect all IDs of terminal/non-terminal branches
terms = Int64[];
splits = Int64[];
if assimflag == "no"
    for i = 1:length(system.branches.ID)
        if isempty(system.branches.children[i])
            push!(terms,i);
        else
            push!(splits,i);
        end
    end
elseif assimflag == "yes"
    for i = 1:length(systems[1].branches.ID)
        if isempty(systems[1].branches.children[i])
            push!(terms,i);
        else
            push!(splits,i);
        end
    end
    term_itr = [terms for i=1:ensemblesize];
    split_itr = [splits for i=1:ensemblesize];
end

# ensemble distributions, allocators, and scalings
if assimflag == "yes"
    nparams = 15;
    nstates = 3;
    nmeas = 3;
    errors = LegHemodynamics.Errors(nparams);

    X = [zeros(nstates) for i=1:ensemblesize];
    θ = [zeros(nparams) for i=1:ensemblesize];
    Y = [zeros(nmeas) for i=1:ensemblesize];

    xhat = zeros(nstates);
    θhat = zeros(nparams);
    yhat = zeros(nmeas);

    yi = [zeros(nmeas) for i=1:ensemblesize];

    # parameter scalings
    θs = zeros(nparams);
    for i=1:ensemblesize
        θs[1] += systems[i].branches.term[14].R[2]; # anterior
        θs[2] += systems[i].branches.term[24].R[2]; # peroneal
        θs[3] += systems[i].branches.term[32].R[2]; # posterior
        θs[4] += systems[i].branches.term[14].C[1];
        θs[5] += systems[i].branches.term[24].C[1];
        θs[6] += systems[i].branches.term[32].C[1];
        θs[7] += systems[i].branches.termscalings[1];
        θs[8] += systems[i].branches.termscalings[2];
        θs[9] += systems[i].branches.termscalings[3];
        θs[10] += systems[i].branches.termscalings[4];
        θs[11] += systems[i].branches.termscalings[5];
        θs[12] += systems[i].branches.termscalings[6];
        # θs[9] += systems[i].branches.A0[14];
        # θs[10] += systems[i].branches.A0[24];
        # θs[11] += systems[i].branches.A0[31];
        θs[13] += systems[i].branches.beta[14];
        θs[14] += systems[i].branches.beta[24];
        θs[15] += systems[i].branches.beta[32];
        # θs[9] += systems[i].branches.YModinPa[14];
        # θs[10] += systems[i].branches.YModinPa[24];
        # θs[11] += systems[i].branches.YModinPa[31];
    end
    θs /= ensemblesize;

    # RTPS allocators
    p = 0.5; # RTPS relaxation amount, ∃ [0.5,0.95]
    c = zeros(ensemblesize);
    σxb = zeros(nstates);
    σxa = zeros(nstates);
    σtb = zeros(nparams);
    σta = ones(nparams);

    tout = Float64[];
    xout = Vector{Float64}[];
    xoutv = Vector{Float64}[];
    yout = Vector{Float64}[];
    youtv = Vector{Float64}[];
    innov = Vector{Float64}[];
    ioutv = Vector{Float64}[];
    θout = Vector{Float64}[];
    θoutv = Vector{Float64}[];
    Pθout = Vector{Float64}[];
    Pθoutv = Vector{Float64}[];
    Pxout = Vector{Float64}[];
    Pxoutv = Vector{Float64}[];
    lbx = Vector{Float64}[];
    ubx = Vector{Float64}[];
    lbxv = Vector{Float64}[];
    ubxv = Vector{Float64}[];

    # default parameter values for non-measured terminal constraints
    Rdefault = 0.75*[77.6909,141.4980,114.7901,102.6618,176.7149,257.4914,18.6383,
        107.4268,291.8191,191.9240,149.3493,35.9008,53.0083,150.1083,259.9696,
        142.5643,95.1046]*LegHemodynamics.gTokg/(LegHemodynamics.mmTom^4);
    Cdefault = [0.08972235,0.04926296,0.06072485,0.06789874,0.03944552,0.02707124,
        0.39707822,0.06488707,0.02388675,0.03631963,0.04667319,0.19569034,0.13150028,
        0.04643721,0.02681318,0.04830696,0.00723676]*(LegHemodynamics.mmTom^4)/LegHemodynamics.gTokg;
end

# time step counter
if assimflag == "no"
    n = system.solverparams.nstart;
elseif assimflag == "yes"
    n = [systems[1].solverparams.nstart for i=1:ensemblesize];
end

# solver loop
tic();
if assimflag == "no"
    while system.solverparams.numbeats < system.solverparams.numbeatstotal
        # time within heart cycle
        tp = system.t[n+1] - sum(system.solverparams.th*system.solverparams.numbeats);
        if mod(n,1000) == 0
            println("Reached time step $n. t = $tp s.")
        end
        # interpolate proximal flow velocity from patient data
        uprox = spl(tp);
        # TVD RK3 time integration
        LegHemodynamics.tvdrk3!(system,times,n,splits,terms,-uprox);
        tic();
        # arterial pressure update
        LegHemodynamics.arterialpressure!(system,n);
        # check for start of new cardiac cycle
        oldnumbeats = system.solverparams.numbeats;
        LegHemodynamics.setnumbeats!(system,n);
        if oldnumbeats < system.solverparams.numbeats
            # update time discretization
            if system.solverparams.numbeats < system.solverparams.numbeatstotal
                LegHemodynamics.updatediscretization!(system);
            end
        end
        n+=1;
        times.tr += toq();
    end
    times.tt += toc();
elseif assimflag == "yes"
    nsamp = [300 for i=1:ensemblesize];
    while systems[1].solverparams.numbeats < systems[1].solverparams.numbeatstotal
        # normal time integration between measurements
        rflag = "predict";
        soln = pmap((a1,a2,a3,a4,a5,a6,a7)->LegHemodynamics.advancetime!(a1,a2,a3,a4,a5,a6,a7;
            runtype=rflag),systems,spl,times,n,nsamp,term_itr,split_itr);
        systems = [soln[i][1] for i=1:ensemblesize];
        n = [soln[i][2] for i=1:ensemblesize];
        times = [soln[i][3] for i=1:ensemblesize];
        # data assimilation
        if assimflag == "yes" && mod((n[1]-1)/nsamp[1],1) == 0
            println("Current time step: $(n[1])")
            println("Current time: $(systems[1].t[n[1]])")

            # time within heart cycle
            tp = systems[1].t[n[1]] - sum(systems[1].solverparams.th*systems[1].solverparams.numbeats);
            # non-dimensional measurement (velocity)
            y = -[splant(tp),splpost(tp),splper(tp)]/LegHemodynamics.us;
            # y = -[splant(tp)]/LegHemodynamics.us;
            println("Scaled velocity measurements at t = $tp s: $y")

            # # generate measurement replicates
            # for i = 1:length(systems)
            #     yi[i][1] = y[1] + rand(Distributions.Normal(0,errors.odev[1]/LegHemodynamics.us));
            #     yi[i][2] = y[2] + rand(Distributions.Normal(0,errors.odev[2]/LegHemodynamics.us));
            #     yi[i][3] = y[3] + rand(Distributions.Normal(0,errors.odev[3]/LegHemodynamics.us));
            # end
            #
            # # forecast parameters and means
            # for i = 1:length(systems)
            #     θ[i][1] = systems[i].branches.term[14].R[2];
            #     θ[i][2] = systems[i].branches.term[24].R[2];
            #     θ[i][3] = systems[i].branches.term[32].R[2];
            #     θ[i][4] = systems[i].branches.term[14].C[1];
            #     θ[i][5] = systems[i].branches.term[24].C[1];
            #     θ[i][6] = systems[i].branches.term[32].C[1];
            #     θ[i][7] = systems[i].branches.termscalings[1];
            #     θ[i][8] = systems[i].branches.termscalings[2];
            #     θ[i][9] = systems[i].branches.termscalings[3];
            #     θ[i][10] = systems[i].branches.termscalings[4];
            #     θ[i][11] = systems[i].branches.termscalings[5];
            #     θ[i][12] = systems[i].branches.termscalings[6];
            #     # θ[i][9] = systems[i].branches.A0[14];
            #     # θ[i][10] = systems[i].branches.A0[24];
            #     # θ[i][11] = systems[i].branches.A0[31];
            #     θ[i][13] = systems[i].branches.beta[14];
            #     θ[i][14] = systems[i].branches.beta[24];
            #     θ[i][15] = systems[i].branches.beta[32];
            #     # θ[i][9] = systems[i].branches.YModinPa[14];
            #     # θ[i][10] = systems[i].branches.YModinPa[24];
            #     # θ[i][11] = systems[i].branches.YModinPa[31];
            # end
            #
            # # forecast parameters (dimensional)
            # θ = [LegHemodynamics.paramwalk!(errors,θs,θ[i]) for i in 1:length(systems)];
            #
            # # non-dimensionalize parameters
            # for i = 1:length(systems)
            #     θ[i] = θ[i]./θs;
            # end
            #
            # # RTPS prior standard deviation
            # for i in 1:length(θ[1])
            #     for j in 1:length(systems)
            #         c[j] = θ[j][i];
            #     end
            #     σtb[i] = std(c;corrected=true);
            # end
            #
            # # parameters back into model
            # for i = 1:length(soln)
            #     systems[i].branches.term[14].R[2] = θ[i][1]*θs[1];
            #     systems[i].branches.term[24].R[2] = θ[i][2]*θs[2];
            #     systems[i].branches.term[32].R[2] = θ[i][3]*θs[3];
            #     systems[i].branches.term[14].C[1] = θ[i][4]*θs[4];
            #     systems[i].branches.term[24].C[1] = θ[i][5]*θs[5];
            #     # systems[i].branches.term[24].C[1] = 2.5e-11;
            #     systems[i].branches.term[32].C[1] = θ[i][6]*θs[6];
            #     systems[i].branches.termscalings[1] = θ[i][7]*θs[7];
            #     systems[i].branches.termscalings[2] = θ[i][8]*θs[8];
            #     systems[i].branches.termscalings[3] = θ[i][9]*θs[9];
            #     systems[i].branches.termscalings[4] = θ[i][10]*θs[10];
            #     systems[i].branches.termscalings[5] = θ[i][11]*θs[11];
            #     systems[i].branches.termscalings[6] = θ[i][12]*θs[12];
            #     # k = 1;
            #     for j in [2,4;] # popliteal
            #         systems[i].branches.term[j].C[1] *= 1/3*(systems[i].branches.termscalings[1]+
            #             systems[i].branches.termscalings[3]+systems[i].branches.termscalings[5]);
            #         systems[i].branches.term[j].R[2] *= 1/3*(systems[i].branches.termscalings[2]+
            #             systems[i].branches.termscalings[4]+systems[i].branches.termscalings[6]);
            #         # systems[i].branches.term[j].C[1] = 1/3*(systems[i].branches.termscalings[1]+
            #             # systems[i].branches.termscalings[3]+systems[i].branches.termscalings[5])*Cdefault[k];
            #         # systems[i].branches.term[j].R[2] = 1/3*(systems[i].branches.termscalings[2]+
            #             # systems[i].branches.termscalings[4]+systems[i].branches.termscalings[6])*Rdefault[k];
            #         # k+=1;
            #     end
            #     for j in [7,9,11,13;] # anterior
            #         systems[i].branches.term[j].C[1] *= systems[i].branches.termscalings[1];
            #         systems[i].branches.term[j].R[2] *= systems[i].branches.termscalings[2];
            #         # systems[i].branches.term[j].C[1] = systems[i].branches.termscalings[1]*Cdefault[k];
            #         # systems[i].branches.term[j].R[2] = systems[i].branches.termscalings[2]*Rdefault[k];
            #         # k+=1;
            #     end
            #     for j in [17,19,21,23;] # posterior
            #         systems[i].branches.term[j].C[1] *= systems[i].branches.termscalings[3];
            #         systems[i].branches.term[j].R[2] *= systems[i].branches.termscalings[4];
            #         # systems[i].branches.term[j].C[1] = systems[i].branches.termscalings[3]*Cdefault[k];
            #         # systems[i].branches.term[j].R[2] = systems[i].branches.termscalings[4]*Rdefault[k];
            #         # k+=1;
            #     end
            #     for j in [25,27,29,31;] # peroneal
            #         systems[i].branches.term[j].C[1] *= systems[i].branches.termscalings[5];
            #         systems[i].branches.term[j].R[2] *= systems[i].branches.termscalings[6];
            #         # systems[i].branches.term[j].C[1] = systems[i].branches.termscalings[5]*Cdefault[k];
            #         # systems[i].branches.term[j].R[2] = systems[i].branches.termscalings[6]*Rdefault[k];
            #         # k+=1;
            #     end
            #     for j in [1:4;] # popliteal
            #         systems[i].branches.beta[j] *= 1/3*(θ[i][13]*θs[13]/systems[i].branches.beta[14]+
            #             θ[i][14]*θs[14]/systems[i].branches.beta[24]+θ[i][15]*θs[15]/systems[i].branches.beta[32]);
            #     end
            #     for j in [5,7:13;] # anteriortibial
            #         systems[i].branches.beta[j] *= θ[i][13]*θs[13]/systems[i].branches.beta[14];
            #     end
            #     for j in [6,15,17:23;] # posteriortibial
            #         systems[i].branches.beta[j] *= θ[i][14]*θs[14]/systems[i].branches.beta[24];
            #     end
            #     for j in [16,25:31;] # peroneal
            #         systems[i].branches.beta[j] *= θ[i][15]*θs[15]/systems[i].branches.beta[32];
            #     end
            #     systems[i].branches.beta[14] = θ[i][13]*θs[13];
            #     systems[i].branches.beta[24] = θ[i][14]*θs[14];
            #     systems[i].branches.beta[32] = θ[i][15]*θs[15];
            #     # systems[i].branches.A0[14] = θ[i][9]*θs[9];
            #     # systems[i].branches.A0[24] = θ[i][10]*θs[10];
            #     # systems[i].branches.A0[31] = θ[i][11]*θs[11];
            #     # for j in [14,24,32]
            #     #     systems[i].branches.term[j].R[1] = systems[i].solverparams.rho*
            #     #         systems[i].branches.c0[j][end]/systems[i].branches.A0[j];
            #     # end
            #     for j = 1:length(systems[i].branches.ID)
            #         systems[i].branches.c0[j][end] = sqrt(0.5*systems[i].branches.beta[j]/
            #             systems[i].solverparams.rho)*systems[i].branches.A0[j]^0.25
            #         if isempty(systems[i].branches.children[j])
            #             systems[i].branches.term[j].R[1] = systems[i].solverparams.rho*
            #                 systems[i].branches.c0[j][end]/systems[i].branches.A0[j];
            #         end
            #         # systems[i].branches.beta[j] = sqrt(pi)*systems[i].branches.thicknessinmm[j]*LegHemodynamics.mmTom*
            #             # systems[i].branches.YModinPa[j]/((1-0.5^2)*systems[i].branches.A0[j]);
            #     end
            # end
            #
            # # find new time step to satisfy CFL condition for all members
            # h = [];
            # for i = 1:length(soln)
            #     for j in 1:length(systems[i].branches.ID)
            #         push!(h,systems[i].solverparams.CFL*systems[i].branches.k[j]/
            #             systems[i].branches.c0[j][end]);
            #     end
            # end
            # if minimum(h) != systems[1].solverparams.h
            #     # println("Old time step size: $(systems[1].solverparams.h) s")
            #     for i = 1:length(soln)
            #         systems[i].solverparams.h = minimum(h);
            #     end
            #     # println("New time step size: $(systems[1].solverparams.h) s")
            #     # re-discretize
            #     if systems[1].solverparams.numbeats < systems[1].solverparams.numbeatstotal
            #         systems = pmap((a1,a2)->LegHemodynamics.rediscretizet!(a1,a2),systems,n);
            #     end
            # end
            #
            # # forecast step
            # rflag = "assim";
            # soln = pmap((a1,a2,a3,a4,a5,a6,a7)->LegHemodynamics.advancetime!(a1,a2,a3,a4,a5,a6,a7;
            #     runtype=rflag),systems,spl,times,n,nsamp,term_itr,split_itr);
            # systems = [soln[i][1] for i=1:ensemblesize];
            # n = [soln[i][2] for i=1:ensemblesize];
            # times = [soln[i][3] for i=1:ensemblesize];
            #
            # # forecast state, measurements
            # X = [zeros(nstates) for i=1:ensemblesize];
            # Y = [zeros(nmeas) for i=1:ensemblesize];
            #
            # for i = 1:ensemblesize
            #     X[i][1] = systems[i].branches.Q[14][1,n[1]]/
            #         (systems[i].branches.A[14][1,n[1]]*LegHemodynamics.us);
            #     X[i][2] = systems[i].branches.Q[24][1,n[1]]/
            #         (systems[i].branches.A[24][1,n[1]]*LegHemodynamics.us);
            #     X[i][3] = systems[i].branches.Q[32][1,n[1]]/
            #         (systems[i].branches.A[32][1,n[1]]*LegHemodynamics.us);
            #     Y[i][1] = systems[i].branches.Q[14][1,n[1]]/
            #         (systems[i].branches.A[14][1,n[1]]*LegHemodynamics.us);
            #     Y[i][2] = systems[i].branches.Q[24][1,n[1]]/
            #         (systems[i].branches.A[24][1,n[1]]*LegHemodynamics.us);
            #     Y[i][3] = systems[i].branches.Q[32][1,n[1]]/
            #         (systems[i].branches.A[32][1,n[1]]*LegHemodynamics.us);
            # end
            #
            # # forecast mean state, parameters, measurement
            # xhat = mean(X);
            # yhat = mean(Y);
            # θhat = mean(θ);
            # # println("Normalized ̂x, first forecast: $xhat")
            # println("Normalized ̂y, first forecast: $yhat")
            # println("Normalized ̂θ, first forecast: $θhat")
            #
            # # forecast params./meas. cross covariance, meas. covariance
            # Pty = zeros(length(θhat),length(yhat))
            # Pyy = zeros(length(yhat),length(yhat))
            # for i = 1:length(soln)
            #     Pty += *((θ[i] .- θhat),(Y[i] .- yhat)');
            #     Pyy += *((Y[i] .- yhat),(Y[i] .- yhat)');
            # end
            # Pty ./= (ensemblesize);
            # Pyy ./= (ensemblesize);
            #
            # # println("Parameter-measurement cross-covariance:")
            # # display(Pty)
            #
            # # add meas. noise to meas. covariance (allows invertibility)
            # Pyy += diagm([(errors.odev[1]/LegHemodynamics.us)^2,
            #     (errors.odev[2]/LegHemodynamics.us)^2,
            #     (errors.odev[3]/LegHemodynamics.us)^2],0);
            # # Pyy += diagm([(errors.odev[1]/LegHemodynamics.us)^2],0);
            #
            # # println("Measurement covariance:")
            # # display(Pyy)
            #
            # # parameter Kalman gain
            # K = Pty*inv(Pyy);
            # # println("Parameter Kalman gain:")
            # # display(K)
            #
            # # parameter analysis step
            # for i = 1:length(soln)
            #     # println("ith measurement replicate: $(yi[i])")
            #     # println("ith forecast measurement: $(Y[i])")
            #     θ[i][:] += K*(yi[i] .- Y[i]);
            # end
            #
            # # RTPS parameter covariance inflation
            # for i in 1:length(θ[1])
            #     for j in 1:length(systems)
            #         c[j] = θ[j][i];
            #     end
            #     σta[i] = std(c;corrected=true);
            # end
            # θhat = mean(θ);
            # for i = 1:length(soln)
            #     θ[i] .= θ[i] .+ p.*((σtb.-σta)./σta).*(θ[i].-θhat);
            # end
            #
            # # analysis parameters back into ensemble members
            # for i = 1:length(soln)
            #     for j = 1:nparams
            #         if θ[i][j]*θs[j] <= errors.lb[j]
            #             # println("Warning: analysis θ_$j below lower bound for member $i.
            #                 # Setting to lower bound of $(errors.lb[j]).")
            #             θ[i][j] = errors.lb[j]/θs[j];
            #         end
            #         if θ[i][j]*θs[j] >= errors.ub[j]
            #             # println("Warning: analysis θ_$j above upper bound for member $i.
            #                 # Setting to upper bound of $(errors.ub[j]).")
            #             θ[i][j] = errors.ub[j]/θs[j];
            #         end
            #     end
            #     systems[i].branches.term[14].R[2] = θ[i][1]*θs[1];
            #     systems[i].branches.term[24].R[2] = θ[i][2]*θs[2];
            #     systems[i].branches.term[32].R[2] = θ[i][3]*θs[3];
            #     systems[i].branches.term[14].C[1] = θ[i][4]*θs[4];
            #     systems[i].branches.term[24].C[1] = θ[i][5]*θs[5];
            #     # systems[i].branches.term[24].C[1] = 2.5e-11;
            #     systems[i].branches.term[32].C[1] = θ[i][6]*θs[6];
            #     systems[i].branches.termscalings[1] = θ[i][7]*θs[7];
            #     systems[i].branches.termscalings[2] = θ[i][8]*θs[8];
            #     systems[i].branches.termscalings[3] = θ[i][9]*θs[9];
            #     systems[i].branches.termscalings[4] = θ[i][10]*θs[10];
            #     systems[i].branches.termscalings[5] = θ[i][11]*θs[11];
            #     systems[i].branches.termscalings[6] = θ[i][12]*θs[12];
            #     # k = 1;
            #     for j in [2,4;] # popliteal
            #         systems[i].branches.term[j].C[1] *= 1/3*(systems[i].branches.termscalings[1]+
            #             systems[i].branches.termscalings[3]+systems[i].branches.termscalings[5]);
            #         systems[i].branches.term[j].R[2] *= 1/3*(systems[i].branches.termscalings[2]+
            #             systems[i].branches.termscalings[4]+systems[i].branches.termscalings[6]);
            #         # systems[i].branches.term[j].C[1] = 1/3*(systems[i].branches.termscalings[1]+
            #             # systems[i].branches.termscalings[3]+systems[i].branches.termscalings[5])*Cdefault[k];
            #         # systems[i].branches.term[j].R[2] = 1/3*(systems[i].branches.termscalings[2]+
            #             # systems[i].branches.termscalings[4]+systems[i].branches.termscalings[6])*Rdefault[k];
            #         # k+=1;
            #     end
            #     for j in [7,9,11,13;] # anterior
            #         systems[i].branches.term[j].C[1] *= systems[i].branches.termscalings[1];
            #         systems[i].branches.term[j].R[2] *= systems[i].branches.termscalings[2];
            #         # systems[i].branches.term[j].C[1] = systems[i].branches.termscalings[1]*Cdefault[k];
            #         # systems[i].branches.term[j].R[2] = systems[i].branches.termscalings[2]*Rdefault[k];
            #         # k+=1;
            #     end
            #     for j in [17,19,21,23;] # posterior
            #         systems[i].branches.term[j].C[1] *= systems[i].branches.termscalings[3];
            #         systems[i].branches.term[j].R[2] *= systems[i].branches.termscalings[4];
            #         # systems[i].branches.term[j].C[1] = systems[i].branches.termscalings[3]*Cdefault[k];
            #         # systems[i].branches.term[j].R[2] = systems[i].branches.termscalings[4]*Rdefault[k];
            #         # k+=1;
            #     end
            #     for j in [25,27,29,31;] # peroneal
            #         systems[i].branches.term[j].C[1] *= systems[i].branches.termscalings[5];
            #         systems[i].branches.term[j].R[2] *= systems[i].branches.termscalings[6];
            #         # systems[i].branches.term[j].C[1] = systems[i].branches.termscalings[5]*Cdefault[k];
            #         # systems[i].branches.term[j].R[2] = systems[i].branches.termscalings[6]*Rdefault[k];
            #         # k+=1;
            #     end
            #     for j = 1:length(term_itr[1])
            #         if term_itr[i][j] != 14 && term_itr[i][j] != 24 && term_itr[i][j] != 32
            #             # systems[i].branches.term[term_itr[i][j]].C[1] *= systems[i].branches.termscalings[1];
            #             # if systems[i].branches.term[term_itr[i][j]].C[1] < 1e-11
            #             #     systems[i].branches.term[term_itr[i][j]].C[1] = 1e-11;
            #             # end
            #             if systems[i].branches.term[term_itr[i][j]].C[1] < 0.5*Cdefault[j]
            #                 # println("Warning: Analysis terminal C < lower bound for member $i,
            #                     # artery $j. Setting to 1e-11.")
            #                 systems[i].branches.term[term_itr[i][j]].C[1] = 0.5*Cdefault[j];
            #             elseif systems[i].branches.term[term_itr[i][j]].C[1] < 1e-12
            #                 systems[i].branches.term[term_itr[i][j]].C[1] = 1e-12;
            #             elseif systems[i].branches.term[term_itr[i][j]].C[1] > 2*Cdefault[j];
            #                 systems[i].branches.term[term_itr[i][j]].C[1] = 2*Cdefault[j];
            #             end
            #             # systems[i].branches.term[term_itr[i][j]].R[2] *= systems[i].branches.termscalings[2];
            #             if systems[i].branches.term[term_itr[i][j]].R[2] < 0.5*Rdefault[j]
            #                 systems[i].branches.term[term_itr[i][j]].R[2] = 0.5*Rdefault[j];
            #             elseif systems[i].branches.term[term_itr[i][j]].R[2] > 2*Rdefault[j]
            #                 systems[i].branches.term[term_itr[i][j]].R[2] = 2*Rdefault[j];
            #             end
            #         end
            #     end
            #     for j in [1:4;] # popliteal
            #         systems[i].branches.beta[j] *= 1/3*(θ[i][13]*θs[13]/systems[i].branches.beta[14]+
            #             θ[i][14]*θs[14]/systems[i].branches.beta[24]+θ[i][15]*θs[15]/systems[i].branches.beta[32]);
            #     end
            #     for j in [5,7:13;] # anteriortibial
            #         systems[i].branches.beta[j] *= θ[i][13]*θs[13]/systems[i].branches.beta[14];
            #     end
            #     for j in [6,15,17:23;] # posteriortibial
            #         systems[i].branches.beta[j] *= θ[i][14]*θs[14]/systems[i].branches.beta[24];
            #     end
            #     for j in [16,25:31;] # peroneal
            #         systems[i].branches.beta[j] *= θ[i][15]*θs[15]/systems[i].branches.beta[32];
            #     end
            #     systems[i].branches.beta[14] = θ[i][13]*θs[13];
            #     systems[i].branches.beta[24] = θ[i][14]*θs[14];
            #     systems[i].branches.beta[32] = θ[i][15]*θs[15];
            #     # systems[i].branches.A0[14] = θ[i][9]*θs[9];
            #     # systems[i].branches.A0[24] = θ[i][10]*θs[10];
            #     # systems[i].branches.A0[31] = θ[i][11]*θs[11];
            #     for j = 1:length(systems[i].branches.ID)
            #         # systems[i].branches.beta[j] = sqrt(pi)*systems[i].branches.thicknessinmm[j]*LegHemodynamics.mmTom*
            #             # systems[i].branches.YModinPa[j]/((1-0.5^2)*systems[i].branches.A0[j]);
            #         systems[i].branches.c0[j][end] = sqrt(0.5*systems[i].branches.beta[j]/
            #             systems[i].solverparams.rho)*systems[i].branches.A0[j]^0.25
            #         if isempty(systems[i].branches.children[j])
            #             systems[i].branches.term[j].R[1] = systems[i].solverparams.rho*
            #                     systems[i].branches.c0[j][end]/systems[i].branches.A0[j];
            #         end
            #     end
            #     # for j in [14,24,32]
            #     #     systems[i].branches.term[j].R[1] = systems[i].solverparams.rho*
            #     #             systems[i].branches.c0[j][end]/systems[i].branches.A0[j];
            #     # end
            # end

            # recalculate and output parameter averages (non-dimensional)
            θhat = mean(θ);
            append!(θout,[θhat])
            # println("Normalized analysis parameter averages: $θhat")

            # analysis parameter variance (for post-processing)
            Ptt = zeros(length(θhat))
            for i = 1:length(soln)
                Ptt += (θ[i] .- θhat).^2;
            end
            Ptt ./= ensemblesize;
            append!(Pθout,[Ptt])

            # # find new time step to satisfy CFL condition for all members
            # h = [];
            # for i = 1:length(soln)
            #     for j in 1:length(systems[i].branches.ID)
            #         push!(h,systems[i].solverparams.CFL*systems[i].branches.k[j]/
            #             systems[i].branches.c0[j][end]);
            #     end
            # end
            # if minimum(h) != systems[1].solverparams.h
            #     # println("Old time step size: $(systems[1].solverparams.h) s")
            #     for i = 1:length(soln)
            #         systems[i].solverparams.h = minimum(h);
            #     end
            #     println("New time step size: $(systems[1].solverparams.h) s")
            #     # re-discretize
            #     if systems[1].solverparams.numbeats < systems[1].solverparams.numbeatstotal
            #         systems = pmap((a1,a2)->LegHemodynamics.rediscretizet!(a1,a2),systems,n);
            #     end
            # end
            #
            # # corrected forecast w/ analysis parameters
            # rflag = "assim";
            # soln = pmap((a1,a2,a3,a4,a5,a6,a7)->LegHemodynamics.advancetime!(a1,a2,a3,a4,a5,a6,a7;
            #     runtype=rflag),systems,spl,times,n,nsamp,term_itr,split_itr);
            # systems = [soln[i][1] for i=1:ensemblesize];
            # n = [soln[i][2] for i=1:ensemblesize];
            # times = [soln[i][3] for i=1:ensemblesize];

            # second forecast state, measurements
            X = [zeros(nstates) for i=1:ensemblesize];
            Y = [zeros(nmeas) for i=1:ensemblesize];

            for i = 1:ensemblesize
                X[i][1] = systems[i].branches.Q[14][1,n[1]]/
                    (systems[i].branches.A[14][1,n[1]]*LegHemodynamics.us);
                X[i][2] = systems[i].branches.Q[24][1,n[1]]/
                    (systems[i].branches.A[24][1,n[1]]*LegHemodynamics.us);
                X[i][3] = systems[i].branches.Q[32][1,n[1]]/
                    (systems[i].branches.A[32][1,n[1]]*LegHemodynamics.us);
                Y[i][1] = systems[i].branches.Q[14][1,n[1]]/
                    (systems[i].branches.A[14][1,n[1]]*LegHemodynamics.us);
                Y[i][2] = systems[i].branches.Q[24][1,n[1]]/
                    (systems[i].branches.A[24][1,n[1]]*LegHemodynamics.us);
                Y[i][3] = systems[i].branches.Q[32][1,n[1]]/
                    (systems[i].branches.A[32][1,n[1]]*LegHemodynamics.us);
            end

            # second forecast mean state, measurement
            xhat = mean(X);
            yhat = mean(Y);
            # println("Normalized ̂x, second forecast: $xhat")
            # println("Normalized ̂y, second forecast: $yhat")

            # output second forecast measurement
            append!(yout,[yhat])
            append!(innov,[yhat-y])

            # output ensemble average state
            append!(xout,[xhat])

            # state variance (for post-processing)
            Pxx = zeros(length(xhat))
            for i = 1:length(soln)
                Pxx += (X[i] .- xhat).^2;
            end
            Pxx ./= ensemblesize;
            append!(Pxout,[Pxx])

            # state 2-sigma quantiles (for post-processing)
            for i = 1:length(xhat)
                xq = Float64[];
                for j = 1:length(soln)
                    push!(xq,X[j][i])
                end
                q = quantile(xq,[0.025,0.975]);
                if i == 1
                    append!(lbx,[ones(1)*q[1]])
                    append!(ubx,[ones(1)*q[2]])
                else
                    push!(lbx[end],q[1])
                    push!(ubx[end],q[2])
                end
            end

            # output measurement time
            push!(tout,systems[1].t[n[1]])
        end
    end
    tt = toc();
    for i = 1:ensemblesize
        times[i].tt += tt;
    end
end

# reshape ensemble averages into vectors of individual time series
if assimflag == "yes"
    for i in 1:length(xout[1]) # i indexes state variables
        xo = Float64[];
        Pxo = Float64[];
        lb = Float64[];
        ub = Float64[];
        for j in 1:length(xout) # j indexes time steps
            push!(xo,xout[j][i])
            push!(Pxo,Pxout[j][i])
            push!(lb,lbx[j][i])
            push!(ub,ubx[j][i])
        end
        append!(xoutv,[xo])
        append!(Pxoutv,[Pxo])
        append!(lbxv,[lb])
        append!(ubxv,[ub])
    end

    for i in 1:length(θout[1]) # i indexes state variables
        θo = Float64[];
        Pθo = Float64[];
        for j in 1:length(θout) # j indexes time steps
            push!(θo,θout[j][i])
            push!(Pθo,Pθout[j][i])
        end
        append!(θoutv,[θo])
        append!(Pθoutv,[Pθo])
    end

    for i in 1:length(yout[1])
        yo = Float64[];
        for j in 1:length(yout)
            push!(yo,yout[j][i])
        end
        append!(youtv,[yo])
    end
end

if saveflag == "yes"
    if assimflag == "no"
        file = MAT.matopen(savefile, "w")
        write(file,"system",system)
        close(file)
    elseif assimflag == "yes"
        for i=1:ensemblesize
            file = MAT.matopen(savefiles[i],"w");
            write(file,"system",systems[i])
            close(file)
        end
        vnames = ["t" "x" "Px" "theta" "thetascale" "Pt" "lb" "ub"];
        for i in 1:length(vnames)
            file = MAT.matopen("$(vnames[i])_converged_small_dev_1.mat","w");
            if vnames[i] == "t"
                write(file,"t",tout)
            elseif vnames[i] == "x"
                write(file,"x",xoutv)
            elseif vnames[i] == "Px"
                write(file,"Px",Pxoutv)
            elseif vnames[i] == "theta"
                write(file,"theta",θoutv)
            elseif vnames[i] == "thetascale"
                write(file,"thetascale",θs)
            elseif vnames[i] == "Pt"
                write(file,"Pt",Pθoutv)
            elseif vnames[i] == "lb"
                write(file,"lb",lbxv)
            elseif vnames[i] == "ub"
                write(file,"ub",ubxv)
            end
            close(file)
        end
    end
end

if assimflag == "no"
    if timeflag == "yes"
        return system, n, times
    elseif timeflag == "no"
        return system, n
    end
elseif assimflag == "yes"
    return n,tout,xoutv,youtv,ioutv,θoutv,Pθoutv,Pxoutv,lbxv,ubxv,θs
end

end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
# https://discourse.julialang.org/t/maintaining-a-fixed-size-top-n-values-list/78868/7
struct BoundedBinaryHeap{T, O <: Base.Ordering} <: DataStructures.AbstractHeap{T}
    ordering::O
    valtree::Vector{T}
    n::Int # maximum length

    function BoundedBinaryHeap{T}(n::Integer, ordering::Base.Ordering) where T
        n ≥ 1 || throw(ArgumentError("max heap size $n must be ≥ 1"))
        new{T, typeof(ordering)}(ordering, sizehint!(Vector{T}(), n), n)
    end

    function BoundedBinaryHeap{T}(n::Integer, ordering::Base.Ordering, xs::AbstractVector) where T
        n ≥ length(xs) || throw(ArgumentError("initial array is larger than max heap size $n"))
        valtree = sizehint!(DataStructures.heapify(xs, ordering), n)
        new{T, typeof(ordering)}(ordering, valtree, n)
    end
end

BoundedBinaryHeap(n::Integer, ordering::Base.Ordering, xs::AbstractVector{T}) where T = BoundedBinaryHeap{T}(n, ordering, xs)

BoundedBinaryHeap{T, O}(n::Integer) where {T, O<:Base.Ordering} = BoundedBinaryHeap{T}(n, O())
BoundedBinaryHeap{T, O}(n::Integer, xs::AbstractVector) where {T, O<:Base.Ordering} = BoundedBinaryHeap{T}(n, O(), xs)

Base.length(h::BoundedBinaryHeap) = length(h.valtree)
Base.isempty(h::BoundedBinaryHeap) = isempty(h.valtree)
@inline Base.first(h::BoundedBinaryHeap) = h.valtree[1]

Base.pop!(h::BoundedBinaryHeap) = DataStructures.heappop!(h.valtree, h.ordering)

function Base.push!(h::BoundedBinaryHeap, v)
    if length(h) < h.n
        DataStructures.heappush!(h.valtree, v, h.ordering)
    elseif Base.Order.lt(h.ordering, @inbounds(h.valtree[1]), v)
        DataStructures.percolate_down!(h.valtree, 1, v, h.ordering)
    end
    return h
end

ext: jl


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

# %%
t2j(s) = s
function t2j(s::Tuple)
    a = (:block, :if, :tuple, :parameters, :kw, :quote, :(=), :(.=), :.)
    s[1] in a && return Expr(s[1], t2j.(s[2:end])...)
    s[1] === :lambda && return Expr(:(->), Expr(:tuple, s[2]...), t2j(s[3]))
    s[1] === :for && return Expr(s[1], Expr(:(=), t2j(s[2]), t2j(s[3])), t2j(s[4]))
    s[1] === :using && return Expr(s[1], Expr.(:., s[2:end])...)
    Expr(:call, t2j.(s)...)
end
macro tupp(x) t2j(Core.eval(__module__, x)) end

# %%
(:for, :k, (:(:), 1, 5),
    (:if, (:iseven, :k),
        (:println, "even: ", :k),
        (:println, "odd:  ", :k))) |> t2j

# %%
@tupp (:for, :k, (:(:), 1, 5),
    (:if, (:iseven, :k),
        (:println, "even: ", :k),
        (:println, "odd:  ", :k)))

# %%
(:block,
    (:(=), :U, (:lambda, (:u,), (:u, :u))),
    (:(=), :F, (:lambda, (:u,), (:lambda, (:n, :a, :b),
                (:if, (:(==), :n, 0), :a,
                    ((:u, :u), (:-, :n, 1), :b, (:+, :a, :b)))))),
    (:(=), :Fib, (:U, :F)),
    (:for, :k, (:(:), 1, 10),
        (:println, "Fib(", :k, ") = ", (:Fib, :k, 0, 1)))) |> t2j

# %%
@tupp (:block,
    (:(=), :U, (:lambda, (:u,), (:u, :u))),
    (:(=), :F, (:lambda, (:u,), (:lambda, (:n, :a, :b),
                (:if, (:(==), :n, 0), :a,
                    ((:u, :u), (:-, :n, 1), :b, (:+, :a, :b)))))),
    (:(=), :Fib, (:U, :F)),
    (:for, :k, (:(:), 1, 10),
        (:println, "Fib(", :k, ") = ", (:Fib, :k, 0, 1))))

# %%
(:block,
    (:using, :Plots),
    (:(=), :n, 20),
    (:(=), :x, (:range, 0, 2, (:kw, :length, :n))),
    (:(=), :y, (:+, (:., :sinpi, (:tuple, :x)), (:*, 0.2, (:randn, :n)))),
    (:(=), :xs, (:range, 0, 2, (:kw, :length, 200))),
    (:(=), :X, (:.^, :x, (:transpose, (:(:), 0, 3)))),
    (:(=), :b, (:\, :X, :y)),
    (:scatter, :x, :y, (:kw, :label, "sample")),
    (:plot!, :xs, (:., :sinpi, (:tuple, :xs)),
        (:kw, :label, "sinpi(x)"),
        (:kw, :color, (:quote, :black)),
        (:kw, :ls, (:quote, :dash))),
    (:plot!, :xs, (:., :evalpoly, (:tuple, :xs, (:Ref, :b))),
        (:kw, :label, "degree-3 polynomial"),
        (:kw, :color, 2), (:kw, :lw, 2))) |> t2j

# %%
@tupp (:block,
    (:using, :Plots),
    (:(=), :n, 20),
    (:(=), :x, (:range, 0, 2, (:kw, :length, :n))),
    (:(=), :y, (:+, (:., :sinpi, (:tuple, :x)), (:*, 0.2, (:randn, :n)))),
    (:(=), :xs, (:range, 0, 2, (:kw, :length, 200))),
    (:(=), :X, (:.^, :x, (:transpose, (:(:), 0, 3)))),
    (:(=), :b, (:\, :X, :y)),
    (:scatter, :x, :y, (:kw, :label, "sample")),
    (:plot!, :xs, (:., :sinpi, (:tuple, :xs)),
        (:kw, :label, "sinpi(x)"),
        (:kw, :color, (:quote, :black)),
        (:kw, :ls, (:quote, :dash))),
    (:plot!, :xs, (:., :evalpoly, (:tuple, :xs, (:Ref, :b))),
        (:kw, :label, "degree-3 polynomial"),
        (:kw, :color, 2), (:kw, :lw, 2)))

# %%

max_stars_repo_path: 0017/tupp - tiny TUPle Processor.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0017/tupp - tiny TUPle Processor.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0017/tupp - tiny TUPle Processor.jl

max_forks_repo_name: genkuroki/public


text:
using CategoricalArrays
using DataFrames
using PooledArrays
using Random

fullgc() = (GC.gc(true); GC.gc(true); GC.gc(true); GC.gc(true))

@assert length(ARGS) == 6
@assert ARGS[3] in ["int", "pool", "cat", "str"]
@assert ARGS[4] in ["uniq", "dup", "manydup"]
@assert ARGS[5] in ["sort", "rand"]
@assert ARGS[6] in ["1", "2"]

@info ARGS

llen = parse(Int, ARGS[1])
rlen = parse(Int, ARGS[2])
@assert llen > 1000
@assert rlen > 2000

pad = maximum(length.(string.((llen, rlen))))

if ARGS[3] == "int"
    if ARGS[4] == "uniq"
        col1 = [1:llen;]
        col2 = [1:rlen;]
    elseif ARGS[4] == "dup"
        col1 = repeat(1:llen ÷ 2, inner=2)
        col2 = repeat(1:rlen ÷ 2, inner=2)
    else
        @assert ARGS[4] == "manydup"
        col1 = repeat(1:llen ÷ 20, inner=20)
        col2 = repeat(1:rlen ÷ 20, inner=20)
    end
elseif ARGS[3] == "pool"
    if ARGS[4] == "dup"
        col1 = PooledArray(repeat(string.(1:llen ÷ 2, pad=pad), inner=2))
        col2 = PooledArray(repeat(string.(1:rlen ÷ 2, pad=pad), inner=2))
    else
        @assert ARGS[4] == "manydup"
        col1 = PooledArray(repeat(string.(1:llen ÷ 20, pad=pad), inner=20))
        col2 = PooledArray(repeat(string.(1:rlen ÷ 20, pad=pad), inner=20))
    end
elseif ARGS[3] == "cat"
    if ARGS[4] == "dup"
        col1 = categorical(repeat(string.(1:llen ÷ 2, pad=pad), inner=2))
        col2 = categorical(repeat(string.(1:rlen ÷ 2, pad=pad), inner=2))
    else
        @assert ARGS[4] == "manydup"
        col1 = categorical(repeat(string.(1:llen ÷ 20, pad=pad), inner=20))
        col2 = categorical(repeat(string.(1:rlen ÷ 20, pad=pad), inner=20))
    end
else
    @assert ARGS[3] == "str"
    if ARGS[4] == "uniq"
        col1 = string.(1:llen, pad=pad)
        col2 = string.(1:rlen, pad=pad)
    elseif ARGS[4] == "dup"
        col1 = repeat(string.(1:llen ÷ 2, pad=pad), inner=2)
        col2 = repeat(string.(1:rlen ÷ 2, pad=pad), inner=2)
    else
        @assert ARGS[4] == "manydup"
        col1 = repeat(string.(1:llen ÷ 20, pad=pad), inner=20)
        col2 = repeat(string.(1:rlen ÷ 20, pad=pad), inner=20)
    end
end

Random.seed!(1234)

if ARGS[5] == "rand"
    shuffle!(col1)
    shuffle!(col2)
else
    @assert ARGS[5] == "sort"
end

if ARGS[6] == "1"
    df1 = DataFrame(id1 = col1)
    df2 = DataFrame(id1 = col2)
    innerjoin(df1[1:1000, :], df2[1:2000, :], on=:id1)
    innerjoin(df2[1:2000, :], df1[1:1000, :], on=:id1)
    fullgc()
    @time innerjoin(df1, df2, on=:id1)
    fullgc()
    @time innerjoin(df2, df1, on=:id1)
else
    @assert ARGS[6] == "2"
    df1 = DataFrame(id1 = col1, id2 = col1)
    df2 = DataFrame(id1 = col1, id2 = col1)
    innerjoin(df1[1:1000, :], df2[1:2000, :], on=[:id1, :id2])
    innerjoin(df2[1:2000, :], df1[1:1000, :], on=[:id1, :id2])
    fullgc()
    @time innerjoin(df1, df2, on=[:id1, :id2])
    fullgc()
    @time innerjoin(df2, df1, on=[:id1, :id2])
end


text:
using Documenter, DataSets

makedocs(;
    modules=[DataSets],
    format=Documenter.HTML(),
    pages=[
        "Introduction" => "index.md",
        "Tutorial" => "tutorial.md",
        "API Reference" => "reference.md",
        "Design Discussion" => "design.md",
    ],
    repo="https://github.com/JuliaComputing/DataSets.jl/blob/{commit}{path}#L{line}",
    sitename="DataSets.jl",
    authors = "Chris Foster and contributors: https://github.com/JuliaComputing/DataSets.jl/graphs/contributors"
)

deploydocs(;
    repo="github.com/JuliaComputing/DataSets.jl",
    push_preview=true
)


text:
module Helper

export initialize!, readData, setModel, sampleParam!
export allocateSim, allocateEquil, allocEqArray, setupLab
using TypeDef
using DelimitedFiles, Statistics, Printf, Random

function readData(
  model ::Int64,       # ID for estimated model to import from input folder
  mcid  ::Int64,       # ID for Monte Carlo estimate to import
  b1    ::Int64,       # 1st bootstrap sample ID
  b2    ::Int64,       # last boostrap sample ID
  vInt  ::Bool=false,  # Vertical integrated market structure
  opt   ::Bool=true,   # Reset output folders
  ms15  ::Bool=false,  # MS expanded by 50%
  mvno  ::Int64=0;     # Retail marginal cost of MVNO
  ver   ::String="824",# Data version (final version by default)
  nSim  ::Int64=200)   # Number of simulation draws for random coefficients

  tol = 1e-14;
  if vInt==false
    sig  = "b"
  else
    sig  = "v"
  end
  pathRef   = "./"
  if ms15==false
    tagms15 = ""
  elseif ms15==true
    tagms15 = "ms15"
  end
  pathFile  = string("dat",ver,"/m",model,tagms15,"/mc",mvno,"/",sig,b1,"/")
  pathImpf  = string("dat",ver,"/m",model,tagms15,"/")
  pathOutp  = string(pathRef,"output/",pathFile)
  fileOutp  = string(pathOutp,"sim",ver,".jld")
  fileMcost = string(pathOutp,"mcost",ver,".jld")
  fileBase  = string(pathOutp,"base",ver,".jld")
  fileLog   = string(pathOutp,"log",ver,".txt")
  pathImp   = string(pathRef,"input/",pathImpf)
  env = envType(model,mcid,nSim,ver,pathOutp,pathImp,fileOutp,fileMcost,
                fileBase,fileLog,tol,b1,b2,mvno)

  if opt==true
    rm(env.pathOut,force=true,recursive=true)
    mkpath(env.pathOut)
  end
  workPath = "$(homedir())/work/kantar/brand/workfiles/"
  impData  = readdlm(workPath*string("demand",ver,tagms15,".csv"),',';header=true)
  impX     = readdlm(workPath*string("Xinput",ver,".csv"),',';header=true)
  impY     = readdlm(workPath*string("income.csv"),',';header=true)
  impParam = readdlm(string(pathImp,"gmmStage2.csv"),',')
  expDelta = readdlm(string(pathImp,"expDelta.csv"))[:,1]
  θ        = convert(Vector{Float64},impParam[mcid,2:(end-2)])

  qtr     = convert(Vector{Int64},impData[1][:,1])
  region  = convert(Vector{Int64},impData[1][:,2])
  network = convert(Vector{Int64},impData[1][:,3])
  brand   = convert(Vector{Int64},impData[1][:,4])
  pop     = convert(Vector{Float64},impData[1][:,5])
  ms      = convert(Vector{Float64},impData[1][:,6])
  ms0     = convert(Vector{Float64},impData[1][:,7])
  price   = convert(Vector{Float64},impData[1][:,8])
  prepost = convert(Vector{Int64},impData[1][:,13])
  idObs   = convert(Vector{Int64},impData[1][:,14])
  idLag   = convert(Vector{Int64},impData[1][:,15])

  Xinput  = impX[1]
  product = sortslices(unique([brand prepost]; dims=1); dims=1)
  set     = setType(unique(region),unique(qtr),product)

  nObs     = size(region,1)
  nProduct = size(product,1)
  prodindx = zeros(Int64,nObs)
  for i=1:nObs
    for j=1:nProduct
      if (brand[i]== set.product[j,1]) .& (prepost[i]==set.product[j,2])
        prodindx[i] = j
      end
    end
    if prodindx[i]==0
      @printf("Cannot find product index for brand %d prepost %d at obs %d\n",
              brand[i],prepost[i],i)
      error("Program stopped.\n")
    end
  end
  dat,dim,θ₂,θ₁,indxeθ = setModel(model,Xinput,region,qtr,network,brand,
    prepost,product,prodindx,ms,price,log.(expDelta),θ,env,set)

  # Set up indx
  mkt = BitArray(undef,dim.obs,0)
  for g in eachindex(set.region)
    for t in eachindex(set.qtr)
      mask = (dat.region.==set.region[g]).&(dat.qtr.==set.qtr[t])
      mkt = [mkt mask]
      push!(dim.prodt, sum(mask))
    end
  end
  if size(mkt,2)!=dim.market
    error("Mask dimension not matching unique markets.")
  end
  market    = zeros(Int64,dim.region)
  active    = falses(dim.product,dim.region)
  anyActive = falses(dim.product)
  mno       = falses(dim.product)
  mvno      = falses(dim.product)
  mvnp      = falses(dim.mvno,3)
  indx = indxType(market,active,anyActive,mno,mvno,mvnp)

  # mat
  mat,mask = allocateMat(impY,model,pop,dat,dim,set,vInt)

  # par0
  αᵢ = -θ₂[1] ./ mat.y
  if (env.model==104) | (env.model==105)
    αᵢ *= 5
  end
  par0 = parType(θ₂,θ₁,αᵢ,copy(dat.δ),indxeθ)

  # sim
  sim = allocateSim(dim)

  return dat, dim, env, mat, par0, set, sim, indx, mask
end

function setModel(
  model   ::Int64,
  Xinput  ::Matrix{Float64},
  region  ::Vector{Int64},
  qtr     ::Vector{Int64},
  network ::Vector{Int64},
  brand   ::Vector{Int64},
  prepost ::Vector{Int64},
  product ::Matrix{Int64},
  prodindx::Vector{Int64},
  ms      ::Vector{Float64},
  price   ::Vector{Float64},
  δ       ::Vector{Float64},
  θ       ::Vector{Float64},
  env     ::envType,
  set     ::setType)

  X       = [Xinput ones(size(Xinput,1))]
  prepaid = convert(Vector{Float64},prepost.==1)
  fbloque = convert(Vector{Float64},(prepost.==3))
  free    = convert(Vector{Float64},brand.==4)
  fightbr = convert(Vector{Float64},(brand.>=5) .& (brand.<=7))
  lowcost = convert(Vector{Float64},(free.==1) .| (fightbr.==1))
  ant4gs  = sum(Xinput[:,5:6],dims=2)
  o1      = ones(length(price))

  nDᵢ = 0
  nνᵢ = 0
  nRᵢ = 0
  indxeθ = Vector{Int64}(undef,0)

  if model==0      # RC logit I in the manuscript
    Xᵢ = [-price fbloque*10 -ant4gs]
    labΣᵢᵢ = ["Price/y";"Fbloque/y";"Ant4G/y"]
    nDᵢ = 3; nRᵢ = 1
  elseif model==15 # RC logit II in the manuscript
    Xᵢ = [-price -ant4gs o1 fbloque*10 prepaid*10]
    labΣᵢᵢ = ["Price/y";"Ant4G/y";"Intercept/y";"Fbloque/y";"Prepaid/y"]
    nDᵢ = 5; nRᵢ = 1
  else
    error("No model number matched.")
  end

  nXᵢ = size(Xᵢ,2)
  θ₂  = θ[1:nXᵢ]
  θ₁  = θ[(nXᵢ+1):end]
  neθ = length(indxeθ)
  θ₂[indxeθ] = exp.(θ₂[indxeθ])
  ξ = δ - X * θ₁
  nθ₂ = length(θ₂)
  nθ₁ = size(X,2)
  nθ  = length(θ)

  nPerson  = env.nSim
  nObs     = size(X,1)
  nTime    = length(set.qtr)
  nRegion  = length(set.region)
  nMarket  = size(unique([region qtr]; dims=1),1)
  nNetwork = length(unique(network))
  nBrand   = length(unique(brand))
  nPrepost = length(unique(prepost))
  nProduct = size(product,1)
  nHost    = 3
  nMno     = 13 # num of MNO products
  nMvno    = 9  # num of MVNO products
  nProdt = Array{Int64,1}(undef,0)

  dat = datType(region,qtr,network,brand,prepost,prodindx,ms,price,δ,X,Xᵢ,ξ)
  dim = dimType(nRegion,nTime,nNetwork,nBrand,nPrepost,nProduct,nPerson,
                nMarket,nObs,nProdt,nMno,nMvno,nθ₂,nθ₁,nθ,nXᵢ,nDᵢ,nνᵢ,nRᵢ)

  return dat,dim,θ₂,θ₁,indxeθ
end

function allocateSim(
  dim::dimType)

  # Sim allocate
  s   = zeros(dim.product)
  p   = zeros(dim.product)
  num = zeros(dim.product)
  D   = zeros(dim.product)
  ms  = zeros(dim.product,dim.region)
  ∂s∂p   = zeros(dim.product,dim.product)
  ∂S∂p   = zeros(dim.product,dim.product)
  ∂D∂p   = zeros(dim.product,dim.product)
  μ̃      = zeros(dim.product,dim.person,dim.region)
  ∂²S∂p² = zeros(dim.product,dim.product,dim.product)
  ∂²D∂p² = zeros(dim.product,dim.product,dim.product)

  sim = simType(s,p,num,D,ms,∂s∂p,∂S∂p,∂D∂p,μ̃,∂²S∂p²,∂²D∂p²)

  return sim
end

function allocateMat(
  impY,
  model::Int64,
  pop  ::Vector{Float64},
  dat  ::datType,
  dim  ::dimType,
  set  ::setType,
  vInt ::Bool=false)

  workPath = "$(homedir())/work/kantar/brand/workfiles/simdraws/"
  csvname  = string(workPath,"lhsnorm",dim.Rᵢ,"d",dim.person,".csv")
  if isfile(csvname)
    impRand = readdlm(csvname,',';header=false)
  else
    # Random.seed!(1234)
    # impRand = randn(dim.person,dim.Rᵢ)
    csvname  = string(workPath,"lhsnorm",dim.Rᵢ,"d200.csv")
    impRand = readdlm(csvname,',';header=false)[1:dim.person,:]
    @printf("Warning: Using user-specified random numbers instead.\n")
  end
  randDraw = impRand[1:dim.person,:]'
  z = randDraw[1,:]
  ν = randDraw[2:end,:]
  y = zeros(dim.person,dim.region)
  ȳ = impY[1]
  if model<100
    for i=1:dim.person
      z[i] = max(z[i], -1.5615)
    end
    for r=1:dim.region
      for i=1:dim.person
        y[i,r] = max(ȳ[r,2] + ȳ[r,3]*z[i], 300)/100
      end
    end
  elseif model>=100
    for r=1:dim.region
      for i=1:dim.person
        y[i,r] = max(exp(ȳ[r,4] + ȳ[r,5]*z[i]), 300)/100
      end
    end
  end

  mat_pop = zeros(dim.region);
  for r=1:dim.region
    g = set.region[r]
    mat_pop[r] = pop[findfirst(dat.region.==g)]
  end
  ms = zeros(dim.product,dim.market)
  p  = zeros(dim.product,dim.qtr)
  μ̃  = zeros(dim.product,dim.person,dim.market)

  active = falses(dim.product,dim.market)
  for n=1:dim.obs
    t = dat.qtr[n]-3
    g = dat.region[n]
    i = dat.product[n]
    r = findfirst(set.region.==g)
    j = dim.region*(t-1) + r
    # @printf("n:%d t:%d r:%d j:%d\n",n,t,r,j)
    active[i,j] = true
  end

  host = [1,1,1,2,2,2,3,3,3,4,1,3,2,1,1,1,2,2,2,3,3,3]
  firm = [sortslices(unique([dat.product dat.network]; dims=1); dims=1) host]

  mno  = (firm[:,2].<=4)
  mvno = .!(mno)
  free = falses(dim.product)
  @. free[firm[:,2]==4] = true
  ∂D∂p = falses(dim.product,dim.product)
  own  = falses(dim.product,dim.product)
  host = falses(dim.product,dim.product)
  if vInt==false # Wholesale-retail
    for j in firm[mno,1]
      @. ∂D∂p[:,j] = (firm[:,2]==firm[j,2])
      @. host[:,j] = (firm[:,3]==firm[j,3])
    end
    for j in firm[mvno,1]
      @. ∂D∂p[:,j] = (firm[:,2]==firm[j,2])
      @. host[:,j] = (firm[:,2]==firm[j,2])
    end
  else # Vertical integration
    for j in firm[mno,1]
      @. ∂D∂p[:,j] = (firm[:,3]==firm[j,3])
      @. host[:,j] = (firm[:,3]==firm[j,3])
    end
    for j in firm[mvno,1]
      @. host[:,j] = (firm[:,3]==firm[j,3])
    end
  end
  user = trues(dim.product)
  mask = maskType(mno,mvno,free,∂D∂p)

  for j=1:dim.product
    @. own[:,j] = (firm[:,2]==firm[j,2])
  end
  mvno = falses(size(firm,1),3)
  mvnp = falses(dim.mvno,3)
  for i=1:3
    @. mvno[:,i] = (firm[:,3]==i) & mask.mvno
    @. mvnp[:,i] = (firm[mask.mvno,3]==i)
  end
  mat = matType(y,ν,mat_pop,ms,p,μ̃,mvno,mvnp,own,host,active,firm)

  return mat, mask
end

function initialize!(
  mat::matType,
  dat::datType,
  dim::dimType,
  par::parType,
  set::setType)

  for n=1:dim.obs
    t = dat.qtr[n]-3
    r = findfirst(set.region.==dat.region[n])
    i = dat.product[n]
    j = dim.region*(t-1)+r
    mat.ms[i,j] = dat.ms[n]
    mutil = par.δ[n]
    for k=1:dim.person
      μ = 0.0
      for l=2:dim.Dᵢ # Heterogeneity net of price
        μ += dat.Xᵢ[n,l] * par.θ₂[l] / mat.y[k,r]
      end
      for l=(dim.Dᵢ+1):(dim.Dᵢ+dim.νᵢ)
        μ += dat.Xᵢ[n,l] * par.θ₂[l] * mat.ν[l-dim.Dᵢ,k]
      end
      mat.μ̃[i,k,j] = mutil + μ
    end
    mat.p[i,t] = dat.price[n]
  end

end

function allocateEquil(
  dim::dimType)

  pdline = trues(dim.product)
  price  = zeros(dim.product,dim.qtr)
  mcost  = [zeros(dim.product,dim.qtr) for t=1:2]
  wprice = zeros(dim.product,dim.qtr)
  markup = zeros(dim.product,dim.qtr,2)
  profit = zeros(dim.product,dim.qtr,2)
  csplus = zeros(dim.qtr)
  demand = zeros(dim.product,dim.qtr)
  fail   = falses(dim.qtr)

  equil = equilType(pdline,price,mcost,wprice,markup,profit,csplus,demand,fail)
  eq    = Array{equilType,1}(undef,16)
  for i=1:16
    eq[i] = deepcopy(equil)
  end
  # Set up pdline corresponding to each scenario
  # Free:  10
  # Sosh:  11
  # Red:   13
  # B&You: 12
  eq[2].pdline[12] = false;
  eq[3].pdline[10] = false;
  @. eq[4].pdline[[10,12]] = false;

  eq[5].pdline[11] = false;
  @. eq[6].pdline[[11,12]] = false;
  @. eq[7].pdline[[10,11]] = false;
  @. eq[8].pdline[[10,11,12]] = false;

  eq[9].pdline[13] = false;
  @. eq[10].pdline[[12,13]] = false;
  @. eq[11].pdline[[10,13]] = false;
  @. eq[12].pdline[[10,12,13]] = false;

  @. eq[13].pdline[[11,13]] = false;
  @. eq[14].pdline[[11,12,13]] = false;
  @. eq[15].pdline[[10,11,13]] = false;
  @. eq[16].pdline[[10,11,12,13]] = false;

  return eq
end

function allocEqArray(
  b1 ::Int64,
  b2 ::Int64,
  dim::dimType)

  nBoot = b2-b1+1
  eqArray = Array{Array{equilType,1},1}(undef,nBoot)
  for i=1:(b2-b1+1)
    eqArray[i] = allocateEquil(dim)
  end
  return eqArray

end

function sampleParam!(
  b   ::Int64,
  dim ::dimType,
  env ::envType,
  mat ::matType,
  par0::parType)

  imp = readdlm(string(env.pathImp,"paramDraws.csv"),',')
  θ = imp[b,1:dim.θ]
  δ = log.(imp[b,(dim.θ+1):end])
  if length(δ)!=dim.obs
    @printf("Length of δ ≂̸ dim.obs at paramDraws.csv.\n")
  end
  α = θ[1]
  θ[par0.indxeθ] = exp.(θ[par0.indxeθ])

  θ₂ = θ[1:dim.θ₂]
  θ₁ = θ[(dim.θ₂+1):end]
  αᵢ = -θ₂[1] ./ mat.y
  if (env.model==104) | (env.model==105)
    αᵢ *= 5
  end
  par = parType(θ₂,θ₁,αᵢ,δ,par0.indxeθ)

  return par
end

function sampleParam!(
  b   ::Int64,
  imp ::Array{Float64,2},
  dim ::dimType,
  env ::envType,
  mat ::matType,
  par0::parType)

  θ = imp[b,1:dim.θ]
  δ = log.(imp[b,(dim.θ+1):end])
  if length(δ)!=dim.obs
    @printf("Length of δ ≂̸ dim.obs at paramDraws.csv.\n")
  end
  α = θ[1]
  θ[par0.indxeθ] = exp.(θ[par0.indxeθ])

  θ₂ = θ[1:dim.θ₂]
  θ₁ = θ[(dim.θ₂+1):end]
  αᵢ = -θ₂[1] ./ mat.y

  par = parType(θ₂,θ₁,αᵢ,δ,par0.indxeθ)

  return par
end

function setupLab()

  network = ["Orange","SFR","Bouygues","Free","MVNO:O","MVNO:S","MVNO:B"]
  plab    = ["Prepaid","Postpaid","F.bloqué"]
  prepost = [repeat(plab,3);repeat(["Postpaid"],4);repeat(plab,3)]
  product = ["Orange","Orange","Orange","SFR","SFR","SFR",
             "Bouygues","Bouygues","Bouygues","Free","Sosh","B&You","Red",
             "MVNO:O","MVNO:O","MVNO:O","MVNO:S","MVNO:S","MVNO:S",
             "MVNO:B","MVNO:B","MVNO:B"]
  prodIdx = [1,2,3,11,4,5,6,13,7,8,9,12,10,14,15,16,17,18,19,20,21,22]
  lab = labType(network,prepost,product,prodIdx)

  return lab

end

end # end of module

max_stars_repo_name: yutec/FightBrand

max_issues_repo_name: yutec/FightBrand

max_forks_repo_name: yutec/FightBrand


text:
"Abstract type representing differences between states."
abstract type Diff end

"""
    combine!(diff::Diff, ds::Diff...)

Combine state differences, modifying the first `Diff` in place.
"""
combine!(d::Diff) = d
combine!(d1::Diff, d2::Diff) = error("Not implemented.")
combine!(d1::Diff, d2::Diff, diffs::Diff...) =
    combine!(combine!(d1, d2), diffs...)

"""
    combine(diff::Diff, ds::Diff...)

Combine state differences, returning a fresh `Diff`.
"""
combine(diffs::Diff...) =
    combine!(empty(diffs[1]), diffs...)

"Convert a `Diff` to an equivalent `Term`."
as_term(diff::Diff) = error("Not implemented.")

"Return whether a `Diff` is redundant (i.e. does nothing)."
is_redundant(diff::Diff) = error("Not implemented.")

"Return an empty `Diff` of the same type as the input `Diff`."
Base.empty(diff::Diff) = error("Not implemented.")


text:
using Random
function initqtree!(wc, i::Integer; backgroundcolor=(0, 0, 0, 0), spacing=getparameter(wc, :spacing))
    img = wc.imgs[i]
    mimg = wordmask(img, backgroundcolor, spacing)
    t = qtree(mimg, wc.params[:groundsize])
    c = isassigned(wc.qtrees, i) ? getcenter(wc.qtrees[i]) : wc.params[:groundsize] ÷ 2
    setcenter!(t, c)
    wc.qtrees[i] = t
end
initqtree!(wc, i; kargs...) = initqtree!.(wc, index(wc, i); kargs...)
"Initialize word's images and other resources with specified style"
function initwords!(wc, i::Integer; backgroundcolor=(0, 0, 0, 0), spacing=getparameter(wc, :spacing),
                    fontsize=getfontsizes(wc, i), color=wc.params[:colors][i],
                    angle=wc.params[:angles][i], font=getfonts(wc, i))
    img, svg = prepareword(wc.words[i], fontsize, color, angle,
        font=font, backgroundcolor=backgroundcolor, border=spacing)
    wc.imgs[i] = img
    wc.svgs[i] = svg
    initqtree!(wc, i, backgroundcolor=backgroundcolor, spacing=spacing)
    nothing
end
initwords!(wc, i; kargs...) = initword!.(wc, index(wc, i); kargs...)
function initwords!(wc::WC; maxiter=5, tolerance=0.02)
    params = wc.params
    weights = wc.weights
    wc.weights .= weights ./ √(sum(weights.^2 .* length.(wc.words)) / length(weights))
    si = sortperm(wc.weights, rev=true)
    wc.words = @view wc.words[si]
    wc.weights = @view wc.weights[si]
    wc.imgs = @view wc.imgs[si]
    wc.svgs = @view wc.svgs[si]
    wc.qtrees = @view wc.qtrees[si]
    wc.params[:colors] = @view wc.params[:colors][si]
    wc.params[:angles] = @view wc.params[:angles][si]
    wc.params[:wordids] = @view wc.params[:wordids][si]
    wc.params[:indsmap] = nothing
    println("set density = $(params[:density])")
    find_weight_scale!(wc, density=params[:density], maxiter=maxiter, tolerance=tolerance)
    printfontsizes(wc)
    initword!.(wc, 1:length(wc.words))
    setstate!(wc, nameof(initwords!))
    wc
end
initword! = initwords!
function printfontsizes(wc)
    nsmall = findlast(i->getfontsizes(wc, i)<=wc.params[:minfontsize], length(wc):-1:1)
    nsmall === nothing && (nsmall = 0)
    println("fontsize ∈ [$(getfontsizes(wc, length(wc))), $(getfontsizes(wc, 1))]")
    if nsmall > 0
        perc = round(Int, nsmall/length(wc)*100)
        println("$nsmall words($perc%) are limited to the minimum font size.")
        if perc > 70
            msg = "It seems too crowded. Word size may be seriously distorted. You need to reduce the number of words or set a larger mask."
            ratio = contentsize_proposal(wc.words, wc.weights) / √getparameter(wc, :contentarea)
            if ratio > 1.1
                msg = msg * " Recommended mask scaling: ratio=$(round(ratio, digits=3))."
            end
            @warn msg
        end

    end
    if getfontsizes(wc, 1) == wc.params[:maxfontsize]
        @warn "Some words are limited to the maximum font size. Please set a `maxfontsize` in `wordcloud` or set a `maxweight` in `processtext`."
    end
end
"""
* placewords!(wc)
* placewords!(wc, style=:uniform)
* placewords!(wc, style=:gathering)
* placewords!(wc, style=:gathering, level=5) #`level` controls the intensity of gathering, typically between 4 and 6, defaults to 5.
* placewords!(wc, style=:gathering, level=6, rt=0) #rt=0, rectangle; rt=1, ellipse; rt=2, rhombus. defaults to 1.
There is also a bool keyword argument `centerlargestword`, which can be set to center the largest word.
When you have set `style=:gathering`, you should disable repositioning in `generate!` at the same time, especially for big words. e.g. `generate!(wc, reposition=0.7)`.
The keyword argument `reorder` is a function to reorder the words, which affects the order of placement. Like `reverse`, `WordCloud.shuffle`.
"""
function placewords!(wc::WC; style=:uniform, rt=:auto, centerlargestword=:auto, reorder=identity, callback=x->x, kargs...)
    if STATEIDS[getstate(wc)] < STATEIDS[:initwords!]
        initwords!(wc)
    end
    @assert style in [:uniform, :gathering]
    if centerlargestword == :auto
        c = wc.params[:groundsize] ÷ 2 # can't wc.maskqtree[1][end÷2]. 1D index goes wrong.
        kernelsize = Stuffing.kernelsize
        centerlargestword = wc.maskqtree[1][c, c] == QTrees.EMPTY && (
            length(wc.qtrees) < 2
            || (length(wc.qtrees) >= 2
                && wc.weights[2] / wc.weights[1] < 0.5
                && prod(kernelsize(wc.qtrees[2])) / prod(kernelsize(wc.qtrees[1])) < 0.5))
        if centerlargestword
            println("center the largest word $(repr(getwords(wc, 1)))")
        end
    end
    arg = ()
    if centerlargestword
        setcenter!(wc.qtrees[1],  wc.params[:groundsize] .÷ 2)
        arg = (2:length(wc.qtrees) |> collect,)
        callback(1)
    end
    qtrees = reorder(wc.qtrees)
    if length(wc.qtrees) > 0 + centerlargestword
        if style == :gathering
            if rt == :auto
                if hasparameter(wc, :rt)
                    rt = getparameter(wc, :rt)
                    println("use wordcloud's parameter `:rt=>$rt` to place words")
                else
                    rt = 1
                    println("rt = 1, ellipse")
                end
            end
            p = min(50, 2 / rt)
            ind = Stuffing.place!(deepcopy(wc.maskqtree), qtrees, arg...;
                    roomfinder=findroom_gathering, p=p, callback=callback, kargs...)
        else
            ind = Stuffing.place!(deepcopy(wc.maskqtree), qtrees, arg...;
                    roomfinder=findroom_uniform, callback=callback, kargs...)
        end
        if ind === nothing error("no room for placement") end
    end
    setstate!(wc, nameof(placewords!))
    setparameter!(wc, 0, :epoch)
    wc
end

"rescale!(wc::WC, ratio::Real)\nRescale all words's size. set `ratio`<1 to shrink, set `ratio`>1 to expand."
function rescale!(wc::WC, ratio::Real)
    qts = wc.qtrees
    centers = getcenter.(qts)
    wc.params[:scale] *= ratio
    initword!.(wc, 1:length(wc))
    setcenter!.(wc.qtrees, centers)
    wc
end

recolor_reset!(wc, i::Integer) = initword!(wc, i)
recolor_reset!(wc, w=:; kargs...) = recolor_reset!.(wc, index(wc, w); kargs...)
function counter(iter; C=Dict{eltype(iter),Int}())
    for e in iter
        C[e] = get(C, e, 0) + 1
    end
    C
end
function mostfrequent(iter; C=Dict{eltype(iter),Int}())
    C = counter(iter; C=C)
    argmax(C)
end
function recolor_main!(wc, i::Integer; background=getmask(wc))
    bg = ARGB.(background)
    img = getimages(wc, i)
    x, y = getpositions(wc, i)
    bg, img = Render.overlappingarea(bg, img, x, y)
    m = wordmask(img, (0, 0, 0, 0), 0)
    bkv = @view bg[m]
    c = mostfrequent(bkv)
    initword!(wc, i, color=c)
end
recolor_main!(wc, w=:; kargs...) = recolor_main!.(wc, index(wc, w); kargs...)
function recolor_average!(wc, i::Integer; background=getmask(wc))
    bg = ARGB.(background)
    img = getimages(wc, i)
    x, y = getpositions(wc, i)
    bg, img = Render.overlappingarea(bg, img, x, y)
    m = wordmask(img, (0, 0, 0, 0), 0)
    bkv = @view bg[m]
    c = sum(bkv) / length(bkv)
    initword!(wc, i, color=c)
end
recolor_average!(wc, w=:; kargs...) = recolor_average!.(wc, index(wc, w); kargs...)

function recolor_blending!(wc, i::Integer; alpha=0.5, background=getmask(wc))
    bg = background
    img = getimages(wc, i)
    x, y = getpositions(wc, i)
    bg, img = Render.overlappingarea(bg, img, x, y)
    m = wordmask(img, (0, 0, 0, 0), 0)
    bg = @view bg[m]
    img = @view img[m]
    alphas = Colors.alpha.(img)
    img .= Render.overlay.(img, convert.(eltype(img), Colors.alphacolor.(bg, alpha)))
    img .= Colors.alphacolor.(img, alphas)
    nothing
end
recolor_blending!(wc, w=:; kargs...) = recolor_blending!.(wc, index(wc, w); kargs...)

function recolor_clipping!(wc, i::Integer; background=getmask(wc))
    bg = background
    img = getimages(wc, i)
    x, y = getpositions(wc, i)
    bg, img = Render.overlappingarea(bg, img, x, y)
    m = wordmask(img, (0, 0, 0, 0), 0)
    bg = @view bg[m]
    img = @view img[m]
    img .= convert.(eltype(img), Colors.alphacolor.(bg, Colors.alpha.(img)))
    nothing
end
recolor_clipping!(wc, w=:; kargs...) = recolor_clipping!.(wc, index(wc, w); kargs...)
"""
recolor the words in `wc` in different styles with the background picture.
The styles supported are `:average`, `:main`, `:clipping`, `:blending`, and :reset (to undo all effects of others).
e.g.
* recolor!(wc, style=:average)
* recolor!(wc, style=:mian)
* recolor!(wc, style=:clipping, background=blur(getmask(wc))) # `background` is optional
* recolor!(wc, style=:blending, alpha=0.3) # `background` and `alpha` are optional
* recolor!(wc, style=:reset)

The effects of `:average`, `:main` and `:clipping` are only determined by the `background`. But the effect of `:blending` is also affected by the previous word color. Therefore, `:blending` can also be used in combination with others
The results of `clipping` and `blending` can not be exported as SVG files, use PNG instead.
"""
function recolor!(wc, args...; style=:average, kargs...)
    if style == :average
        recolor_average!(wc, args...; kargs...)
    elseif style == :main
        recolor_main!(wc, args...; kargs...)
    elseif style == :blending
        recolor_blending!(wc, args...; kargs...)
    elseif style == :clipping
        recolor_clipping!(wc, args...; kargs...)
    elseif style == :reset
        recolor_reset!(wc, args...; kargs...)
    else
        error("unknown style $style")
    end
    nothing
end
"""
# Positional Args
* wc: the wordcloud to fit
* nepoch: training epoch nums
# Keyword Args
* patient: number of epochs before repositioning
* reposition: a Bool value to turn on/off teleport, a Float number `p` between 0~1 indicating the repositioning ratio (Minimum `p`), a Int number `n` equivalent to `i -> i >= n`, a Function index::Int -> doteleport::Boll, or a white list collision.
* trainer: appoint a training engine
"""
function fit!(wc, args...; reposition=true, optimiser=SGD(), krags...)
    reposition isa Union{Function,Number} || (reposition = index(wc, reposition)) # Bool <: Number
    if STATEIDS[getstate(wc)] < STATEIDS[:placewords!]
        placewords!(wc)
    end
    qtrees = [wc.maskqtree, wc.qtrees...]
    ep, nc = train!(qtrees, args...; reposition=reposition, optimiser=optimiser, krags...)
    wc.params[:epoch] += ep
    if nc == 0
        setstate!(wc, nameof(fit!))
        # @assert isempty(outofkernelbounds(wc.maskqtree, wc.qtrees))
        # colllist = first.(totalcollisions(qtrees))
        # @assert length(colllist) == 0
    else
        setstate!(wc, nameof(placewords!))
    end
    wc
end
function printcollisions(wc)
    qtrees = [wc.maskqtree, wc.qtrees...]
    colllist = first.(totalcollisions(qtrees))
    get_text(i) = i > 1 ? wc.words[i - 1] : "#MASK#"
    collwords = [(get_text(i), get_text(j)) for (i, j) in colllist]
    if length(colllist) > 0
        @warn "Have $(length(colllist)) collisions. Try setting a larger `nepoch` and `retry`, or lower `density` in `wordcloud` to fix that"
        println("$collwords")
    end
end

"""
# Positional Args
* wc: the wordcloud to fit
* nepoch: training epoch nums
# Keyword Args
* retry: shrink & retrain times, defaults to 3, set to `1` to disable shrinking
* patient: number of epochs before repositioning
* reposition: a Bool value to turn on/off teleport, a Float number `p` between 0~1 indicating the repositioning ratio (Minimum `p`), a Int number `n` equivalent to `i -> i >= n`, a Function index::Int -> doteleport::Boll, or a white list collision.
* trainer: appoint a training engine
"""
function generate!(wc::WC, args...; retry=3, krags...)
    if STATEIDS[getstate(wc)] < STATEIDS[:placewords!]
        placewords!(wc)
    end
    for r in 1:retry
        if r != 1
            rescale!(wc, 0.97)
            dens = textoccupying(getwords(wc), getfontsizes(wc), getfonts(wc)) / wc.params[:contentarea]
            println("▸$r. try scale = $(wc.params[:scale]). The density is reduced to $dens")
            printfontsizes(wc)
        else
            println("▸$r. scale = $(wc.params[:scale])")
        end
        fit!(wc, args...; krags...)
        if getstate(wc) == :fit!
            break
        end
    end
    if STATEIDS[getstate(wc)] >= STATEIDS[:fit!]
        println("$(wc.params[:epoch]) epochs")
        setstate!(wc, nameof(generate!))
    else # check
        printcollisions(wc)
    end
    wc
end

STATES = nameof.([wordcloud, initwords!, placewords!, fit!, generate!])
STATEIDS = Dict([s => i for (i, s) in enumerate(STATES)])

"""
keep some words and ignore the others, then execute the function. It's the opposite of `ignore`.
* keep(fun, wc, ws::String) #keep a word
* keep(fun, wc, ws::Set{String}) #kepp all words in ws
* keep(fun, wc, ws::Vector{String}) #keep all words in ws
* keep(fun, wc, inds::Union{Integer, Vector{Integer}})
* keep(fun, wc::WC, mask::AbstractArray{Bool}) #keep words. length(mask)==length(wc)
"""
function keep(fun, wc::WC, mask::AbstractArray{Bool})
    mem = [wc.words, wc.weights, wc.imgs, wc.svgs, wc.qtrees,
            wc.params[:colors], wc.params[:angles], wc.params[:wordids], wc.params[:indsmap]]
    wc.words = @view wc.words[mask]
    wc.weights = @view wc.weights[mask]
    wc.imgs = @view wc.imgs[mask]
    wc.svgs = @view wc.svgs[mask]
    wc.qtrees = @view wc.qtrees[mask]
    wc.params[:colors] = @view wc.params[:colors][mask]
    wc.params[:angles] = @view wc.params[:angles][mask]
    wc.params[:wordids] = @view wc.params[:wordids][mask]
    wc.params[:indsmap] = nothing
    r = nothing
    try
        r = fun()
    finally
        wc.words = mem[1]
        wc.weights = mem[2]
        wc.imgs = mem[3]
        wc.svgs = mem[4]
        wc.qtrees = mem[5]
        wc.params[:colors] = mem[6]
        wc.params[:angles] = mem[7]
        wc.params[:wordids] = mem[8]
        wc.params[:indsmap] = mem[9]
    end
    r
end

"""
pin some words as if they were part of the background, then execute the function.
* pin(fun, wc, ws::String) #pin a word
* pin(fun, wc, ws::Set{String}) #pin all words in ws
* pin(fun, wc, ws::Vector{String}) #pin all words in ws
* pin(fun, wc, inds::Union{Integer, Vector{Integer}})
* pin(fun, wc::WC, mask::AbstractArray{Bool}) #pin words. length(mask)==length(wc)
"""
function pin(fun, wc::WC, mask::AbstractArray{Bool})
    maskqtree = wc.maskqtree
    wcmask = wc.mask
    contentarea = wc.params[:contentarea]

    maskqtree2 = deepcopy(maskqtree)
    Stuffing.overlap!(maskqtree2, wc.qtrees[mask])
    wc.maskqtree = maskqtree2
    resultpic = copy(wc.mask)
    wc.mask = overlay!(resultpic, wc.imgs[mask], getpositions(wc, mask))
    wc.params[:contentarea] = occupying(QTrees.kernel(wc.maskqtree[1]), QTrees.FULL)
    r = nothing
    try
        r = ignore(fun, wc, mask)
    finally
        wc.maskqtree = maskqtree
        wc.mask = wcmask
    wc.params[:contentarea] = contentarea
    end
    r
end

function keep(fun, wc, ind::Integer)
    mask = falses(length(wc))
    mask[ind] = 1
    keep(fun, wc, mask)
end
function keep(fun, wc, inds::AbstractArray{<:Integer})
    mask = falses(length(wc))
    mask[inds] .= 1
    keep(fun, wc, mask)
end
function keep(fun, wc, ws)
    keep(fun, wc, index(wc, ws))
end

"""
ignore some words as if they don't exist, then execute the function. It's the opposite of `keep`.
* ignore(fun, wc, ws::String) #ignore a word
* ignore(fun, wc, ws::Set{String}) #ignore all words in ws
* ignore(fun, wc, ws::Vector{String}) #ignore all words in ws
* ignore(fun, wc, inds::Union{Integer, Vector{Integer}})
* ignore(fun, wc::WC, mask::AbstractArray{Bool}) #ignore words. length(mask)==length(wc)
"""
function ignore(fun, wc::WC, mask::AbstractArray{Bool})
    keep(fun, wc, .!mask)
end
function ignore(fun, wc, ind::Integer)
    mask = falses(length(wc))
    mask[ind] = 1
    ignore(fun, wc, mask)
end
function ignore(fun, wc, inds::AbstractArray{<:Integer})
    mask = falses(length(wc))
    mask[inds] .= 1
    ignore(fun, wc, mask)
end
function ignore(fun, wc, ws)
    ignore(fun, wc, index(wc, ws))
end

function pin(fun, wc, ind::Integer)
    mask = falses(length(wc))
    mask[ind] = 1
    pin(fun, wc, mask)
end
function pin(fun, wc, inds::AbstractArray{<:Integer})
    mask = falses(length(wc))
    mask[inds] .= 1
    pin(fun, wc, mask)
end
function pin(fun, wc, ws)
    pin(fun, wc, index(wc, ws))
end


text:
const TEX_TEMPLATE_FOLDER = joinpath(dirname(@__FILE__), "tpl")
const DEFAULT_LATEX_DEST_FOLDER = "./LaTeX"

function latex(data::ModelSelection.ModelSelectionData, originaldata::ModelSelection.ModelSelectionData;
     path::Union{Nothing, String}=DEFAULT_LATEX_DEST_FOLDER)
    tempfolder = tempname()
    mkdir(tempfolder)

    addextras(data, :latex, nothing, tempfolder)
    if size(data.results, 1) > 0
        dict = Dict()
        latex!(dict, data, originaldata)
        for i in 1:size(data.results, 1)
            latex!(dict, data, originaldata, data.results[i])
        end
        create_workspace(tempfolder)
        if data.results[1].ttest
            dict[string(ModelSelection.AllSubsetRegression.ALLSUBSETREGRESSION_EXTRAKEY)]["intelligent_text"] = create_figures(data, tempfolder)
        end
        render_latex(dict, tempfolder)
        zip_folder(tempfolder, path)
    end
    rm(tempfolder, force=true, recursive=true)
end

function get_col_statistics(data)
    nnan = filter(x -> !isnan(x), data)
    if( isempty(nnan) )
        return Dict(
            "avg"      => "",
            "median"   => "",
            "mode"     => "",
            "std"      => "",
            "skew"     => "",
            "kurt"     => ""
        )
    end
    Dict(
        "avg"      => @sprintf("%.3f", mean(nnan)),
        "median"   => @sprintf("%.3f", median(nnan)),
        "mode"     => @sprintf("%.3f", mode(nnan)),
        "std"      => @sprintf("%.3f", std(nnan)),
        "skew"     => @sprintf("%.3f", skewness(nnan)),
        "kurt"     => @sprintf("%.3f", kurtosis(nnan))
    )
end

function latex!(dict::Dict, data::ModelSelection.ModelSelectionData, originaldata::ModelSelection.ModelSelectionData)
    # Preprocessing
    preprocessing_dict = process_dict(data.extras[Preprocessing.PREPROCESSING_EXTRAKEY])
    preprocessing_dict["equation"] = join(map(x -> "$x", filter(x -> x != :_cons, preprocessing_dict["datanames"])), " ")
    preprocessing_dict["datanames"] = string("[:", join(preprocessing_dict["datanames"], ", :"), "]")
    preprocessing_dict["descriptive"] = []

    datanames_index = ModelSelection.create_datanames_index(originaldata.expvars)

    for (i, var) in enumerate(originaldata.expvars)
        orig = originaldata.expvars_data[:, datanames_index[var]]
        obs = collect(skipmissing(orig))
        c_obs = length(obs)
        c_miss = length(orig) - c_obs
        push!(preprocessing_dict["descriptive"], Dict(
            "name" => replace(string(var), "_" => "\\_"),
            "nobs" => length(obs),
            "mean" => @sprintf("%.2f", mean(obs)),
            "std"  => @sprintf("%.2f", std(obs)),
            "max"  => @sprintf("%.2f", maximum(obs)),
            "min"  => @sprintf("%.2f", minimum(obs)),
            "miss" => @sprintf("%.2f", (c_miss/c_obs)*100)
        ))
        if( (i + 1) % 54 == 0 )
            push!(preprocessing_dict["descriptive"], Dict("name" => false))
        end
    end

    if "seasonaladjustment" in keys(preprocessing_dict) && preprocessing_dict["seasonaladjustment"] != nothing
        preprocessing_dict["seasonaladjustment"] = get_array_details(preprocessing_dict["seasonaladjustment"])
    end
    dict[string(Preprocessing.PREPROCESSING_EXTRAKEY)] = preprocessing_dict

    # FeatureExtraction
    if FeatureExtraction.FEATUREEXTRACTION_EXTRAKEY in keys(data.extras)
        data.extras[FeatureExtraction.FEATUREEXTRACTION_EXTRAKEY]
        featureextraction_dict = process_dict(data.extras[FeatureExtraction.FEATUREEXTRACTION_EXTRAKEY])
        if "fe_lag" in keys(featureextraction_dict)
            featureextraction_dict["fe_lag"] =
                get_array_details(featureextraction_dict["fe_lag"])
        end
        if "fe_log" in keys(featureextraction_dict)
            featureextraction_dict["fe_log"] = replace(string(featureextraction_dict["fe_log"]), "Symbol" => "")
        end
        if "fe_inv" in keys(featureextraction_dict)
            featureextraction_dict["fe_inv"] = replace(string(featureextraction_dict["fe_inv"]), "Symbol" => "")
        end
        if "fe_sqr" in keys(featureextraction_dict)
            featureextraction_dict["fe_sqr"] = replace(string(featureextraction_dict["fe_sqr"]), "Symbol" => "")
        end
        if "interaction" in keys(featureextraction_dict)
            featureextraction_dict["interaction"] =
                get_array_simple_details( featureextraction_dict["interaction"])
        end
        if  !("fe_lag" in keys(featureextraction_dict)) &&
            !("fe_log" in keys(featureextraction_dict)) &&
            !("fe_inv" in keys(featureextraction_dict)) &&
            !("fe_sqr" in keys(featureextraction_dict)) &&
            !("interaction" in keys(featureextraction_dict))
            featureextraction_dict = false
        end
        dict[string(FeatureExtraction.FEATUREEXTRACTION_EXTRAKEY)] = featureextraction_dict
    end

    # PreliminarySelection
    if PreliminarySelection.PRELIMINARYSELECTION_EXTRAKEY in keys(data.extras)
        dict[string(PreliminarySelection.PRELIMINARYSELECTION_EXTRAKEY)] = process_dict(data.extras[PreliminarySelection.PRELIMINARYSELECTION_EXTRAKEY])

        betas = dict[string(PreliminarySelection.PRELIMINARYSELECTION_EXTRAKEY)]["lassobetas"]
        betas_dict = []
        for (i, beta) in enumerate(filter(j -> j != 0, betas))
            if( beta != 0)
                push!(betas_dict, Dict("name" => replace(string(data.expvars[i]), "_" => "\\_"), "coef" => @sprintf("%.3f", beta)))
            end
        end

        dict[string(PreliminarySelection.PRELIMINARYSELECTION_EXTRAKEY)]["lassolambda"] =
         @sprintf("%.3f", dict[string(PreliminarySelection.PRELIMINARYSELECTION_EXTRAKEY)]["lassolambda"])

        dict[string(PreliminarySelection.PRELIMINARYSELECTION_EXTRAKEY)]["expvssars"] = betas_dict
    end

    # Output
    if Output.OUTPUT_EXTRAKEY in keys(data.extras)
        dict[string(Output.OUTPUT_EXTRAKEY)] = process_dict(data.extras[Output.OUTPUT_EXTRAKEY])
    end

    return dict
end

function latex!(dict::Dict, data::ModelSelection.ModelSelectionData, originaldata::ModelSelection.ModelSelectionData, result::ModelSelection.AllSubsetRegression.AllSubsetRegressionResult)
    if ModelSelection.AllSubsetRegression.ALLSUBSETREGRESSION_EXTRAKEY in keys(data.extras)
        dict[string(ModelSelection.AllSubsetRegression.ALLSUBSETREGRESSION_EXTRAKEY)] = process_dict(data.extras[ModelSelection.AllSubsetRegression.ALLSUBSETREGRESSION_EXTRAKEY])

        if "fixedvariables" in keys(dict) && size(dict["fixedvariables"], 1) == 0
            delete!(dict["fixedvariables"])
        end

        datanames_index = ModelSelection.create_datanames_index(result.datanames)
        cols = ModelSelection.get_selected_variables(Int64(result.bestresult_data[datanames_index[:index]]), data.expvars, data.intercept)

        modelavg_datanames = ModelSelection.AllSubsetRegression.get_varnames(result.modelavg_datanames)

        d_bestmodel = Dict()
        d_bestmodel["depvar"] = data.depvar
        d_bestmodel["bmexpvars"] = []

        for var in data.expvars
            intercept = if (data.intercept) 1 else 0 end
            d = Dict()
            d["name"] = replace(string(var), "_" => "\\_")
            d["nobs"] = result.nobs
            d["criteria"] = result.criteria

            if (var in data.expvars[cols] && !isnan(result.bestresult_data[datanames_index[Symbol("$(var)_b")]]))
                d["best"] = Dict()
                d["best"]["b"] = @sprintf("%.3f", result.bestresult_data[datanames_index[Symbol("$(var)_b")]])
                if (result.ttest)
                    d["best"]["ttest"] = true
                    d["best"]["bstd"] = @sprintf("%.3f", result.bestresult_data[datanames_index[Symbol("$(var)_bstd")]])
                    t = result.bestresult_data[datanames_index[Symbol("$(var)_t")]]
                    prob_t = pdf(TDist(data.nobs - length(data.expvars) - (if data.intercept 1 else 0 end) ), t)
                    if (prob_t < 0.01)
                        d["best"]["stars"] = "***"
                    else
                        if (prob_t < 0.05)
                            d["best"]["stars"] = "**"
                        else
                            if (prob_t < 0.1)
                                d["best"]["stars"] = "*"
                            else
                                d["best"]["stars"] = ""
                            end
                        end
                    end
                end
            end

            if (result.modelavg && !isnan(result.modelavg_data[datanames_index[Symbol("$(var)_b")]]))
                d_bestmodel["modelavg"] = true
                d["avg"] = Dict()
                d["avg"]["b"] = @sprintf("%.3f", result.modelavg_data[datanames_index[Symbol("$(var)_b")]])
                if (result.ttest)
                    d["avg"]["ttest"] = true
                    d["avg"]["bstd"] = @sprintf("%.3f", result.modelavg_data[datanames_index[Symbol("$(var)_bstd")]])
                    t = result.modelavg_data[datanames_index[Symbol("$(var)_t")]]
                    prob_t = pdf(TDist(data.nobs - length(data.expvars) - (if data.intercept 1 else 0 end) ), t)
                    if (prob_t < 0.01)
                        d["avg"]["stars"] = "***"
                    else
                        if (prob_t < 0.05)
                            d["avg"]["stars"] = "**"
                        else
                            if (prob_t < 0.1)
                                d["avg"]["stars"] = "*"
                            else
                                d["avg"]["stars"] = ""
                            end
                        end
                    end
                end
            end

            if ("best" in keys(d) || "avg" in keys(d))
                push!(d_bestmodel["bmexpvars"], d)
            end
        end

        expvars_dict = map( var -> begin
            t = result.data[:,ModelSelection.get_column_index(Symbol("$(var)_t"), result.datanames)]
            b = result.data[:,ModelSelection.get_column_index(Symbol("$(var)_b"), result.datanames)]

            t_nnan = filter(x -> !isnan(x), t)
            b_nnan = filter(x -> !isnan(x), b)

            t_poscol = filter(x -> x > 0, t_nnan)
            b_poscol = filter(x -> x > 0, b_nnan)

            prob_t = map( it -> pdf(TDist(data.nobs - length(data.expvars) - (if data.intercept 1 else 0 end) ), it), t)
            Dict(
                "name" => var,
                "posshare" => @sprintf("%.3f", size(b_poscol, 1) / size(b_nnan, 1)),
                "sigshare" => @sprintf("%.3f", size(filter(p -> p < 0.1, prob_t), 1) / size(b_nnan, 1)),
                "t" => get_col_statistics(t),
                "b" => get_col_statistics(b)
            )
        end, data.expvars)

        dict[string(ModelSelection.AllSubsetRegression.ALLSUBSETREGRESSION_EXTRAKEY)]["expvars"] = expvars_dict
        dict[string(ModelSelection.AllSubsetRegression.ALLSUBSETREGRESSION_EXTRAKEY)]["expvarswithoutcons"] =
        filter(x -> x["name"] != :_cons, expvars_dict)
        dict[string(ModelSelection.AllSubsetRegression.ALLSUBSETREGRESSION_EXTRAKEY)]["bestmodel"] = d_bestmodel

        if dict[string(ModelSelection.AllSubsetRegression.ALLSUBSETREGRESSION_EXTRAKEY)]["criteria"] != nothing
            dict[string(ModelSelection.AllSubsetRegression.ALLSUBSETREGRESSION_EXTRAKEY)]["criteria"] =
            string("[:", join(dict[string(ModelSelection.AllSubsetRegression.ALLSUBSETREGRESSION_EXTRAKEY)]["criteria"], ", :"), "]")
        end

        if dict[string(ModelSelection.AllSubsetRegression.ALLSUBSETREGRESSION_EXTRAKEY)]["residualtest"] != false
            if "time" in keys(dict[string(ModelSelection.AllSubsetRegression.ALLSUBSETREGRESSION_EXTRAKEY)]) &&
                dict[string(ModelSelection.AllSubsetRegression.ALLSUBSETREGRESSION_EXTRAKEY)]["time"] != nothing
                dict[string(ModelSelection.AllSubsetRegression.ALLSUBSETREGRESSION_EXTRAKEY)]["residualtestfortex2"] = true
            else
                dict[string(ModelSelection.AllSubsetRegression.ALLSUBSETREGRESSION_EXTRAKEY)]["residualtestfortex"] = true
            end
        end
    end
    return dict
end

function latex!(dict::Dict, data::ModelSelection.ModelSelectionData, originaldata::ModelSelection.ModelSelectionData, result::ModelSelection.CrossValidation.CrossValidationResult)
    if ModelSelection.CrossValidation.CROSSVALIDATION_EXTRAKEY in keys(data.extras)
        dict[string(ModelSelection.CrossValidation.CROSSVALIDATION_EXTRAKEY)] =
        process_dict(data.extras[ModelSelection.CrossValidation.CROSSVALIDATION_EXTRAKEY])

        datanames_index = ModelSelection.create_datanames_index(result.datanames)

        d_bestmodel = []

        for var in data.expvars
            intercept = if (data.intercept) 1 else 0 end
            d = Dict()
            d["name"] = replace(string(var), "_" => "\\_")
            if (!isnan(result.average_data[datanames_index[Symbol("$(var)_b")]]))
                d["mean"] = Dict()
                d["mean"]["b"] = @sprintf("%.6f", result.average_data[datanames_index[Symbol("$(var)_b")]])
                if (result.ttest)
                    d["mean"]["ttest"] = true
                    d["mean"]["bstd"] = @sprintf("%.6f", result.average_data[datanames_index[Symbol("$(var)_bstd")]])
                    t = result.average_data[datanames_index[Symbol("$(var)_b")]] / result.average_data[datanames_index[Symbol("$(var)_bstd")]]
                    d["mean"]["t"] = t
                    prob_t = pdf(TDist(data.nobs - length(data.expvars) - (if data.intercept 1 else 0 end) ), t)
                    if (prob_t < 0.01)
                        d["mean"]["stars"] = "***"
                    else
                        if (prob_t < 0.05)
                            d["mean"]["stars"] = "**"
                        else
                            if (prob_t < 0.1)
                                d["mean"]["stars"] = "*"
                            else
                                d["mean"]["stars"] = ""
                            end
                        end
                    end
                end
                d["median"] = Dict()
                d["median"]["b"] = @sprintf("%.6f", result.median_data[datanames_index[Symbol("$(var)_b")]])
                if (result.ttest)
                    d["median"]["ttest"] = true
                    d["median"]["bstd"] = @sprintf("%.6f", result.median_data[datanames_index[Symbol("$(var)_bstd")]])
                    t = result.median_data[datanames_index[Symbol("$(var)_b")]] / result.median_data[datanames_index[Symbol("$(var)_bstd")]]
                    d["median"]["t"] = t
                    prob_t = pdf(TDist(data.nobs - length(data.expvars) - (if data.intercept 1 else 0 end) ), t)
                    if (prob_t < 0.01)
                        d["median"]["stars"] = "***"
                    else
                        if (prob_t < 0.05)
                            d["median"]["stars"] = "**"
                        else
                            if (prob_t < 0.1)
                                d["median"]["stars"] = "*"
                            else
                                d["median"]["stars"] = "-"
                            end
                        end
                    end
                end
            end

            if ("mean" in keys(d) || "median" in keys(d))
                push!(d_bestmodel, d)
            end
        end

        dict[string(ModelSelection.CrossValidation.CROSSVALIDATION_EXTRAKEY)]["kfoldvars"] = d_bestmodel

        dict[string(ModelSelection.CrossValidation.CROSSVALIDATION_EXTRAKEY)]["outsample"] = Dict(
            "median" => @sprintf("%.6f", result.median_data[datanames_index[:rmseout]]),
            "mean" => @sprintf("%.6f", result.average_data[datanames_index[:rmseout]])
        )

    end
    return dict
end

function process_dict(dict)
    dict_str = Dict()
    for key in keys(dict)
        key_str = string(key)
        if dict[key] != nothing
            dict_str[key_str] = dict[key]
        end
    end
    return dict_str
end

"""
Copy the required files from the tpl folder. The existent content will be replaced. The folder is created if not exists.
"""
function create_workspace(destfolder::AbstractString)
    cp(TEX_TEMPLATE_FOLDER, destfolder, force=true)
end

"""
Write template based on dict info
"""
function render_latex(dict, destfolder::AbstractString)
    io = open(joinpath(destfolder, "main.tex"), "w")
    render = render_from_file(joinpath(TEX_TEMPLATE_FOLDER, "../tpl.tex"), dict)
    write(io, render)
    close(io)
end

function zip_folder(sourcefolder::AbstractString, destfolder::AbstractString)
    filenamesforzip = map(x -> joinpath(sourcefolder,x), readdir(sourcefolder))
    create_zip(destfolder, filenamesforzip)
end

dropnans(res, var) = res.results[findall(x -> !isnan(x), res.results[:, res.header[var]]), res.header[var]]

"""
Create required figures by the template. Generate png images into the dest folder
"""
function create_figures(data, destfolder)
    expvars2 = filter(x -> x != :_cons, data.expvars)
    #deleteat!(expvars2, ModelSelection.get_column_index(:_cons, data.expvars))

    criteria_diff = Array{Any}(undef, size(expvars2,1), 2)

    for (i, expvar) in enumerate(expvars2)
        bcol = ModelSelection.get_column_index(Symbol("$(expvar)_b"), data.results[1].datanames)
        tcol = ModelSelection.get_column_index(Symbol("$(expvar)_t"), data.results[1].datanames)
        r2col = ModelSelection.get_column_index(:r2adj, data.results[1].datanames)

        x = data.results[1].data[findall(x -> !isnan(x), data.results[1].data[:, bcol]), bcol]
        y = data.results[1].data[findall(x -> !isnan(x), data.results[1].data[:, tcol]), tcol]

        biden = kde((x, y), npoints=(100,100))

        contour(biden.x, biden.y, biden.density; xlabel="Coef. $expvar", ylabel="t-test $expvar")
        savefig(joinpath(destfolder, "contour_$(expvar)_b_t.png"))

        wireframe(biden.x, biden.y, biden.density; xlabel="Coef. $expvar", ylabel="t-test $expvar", camera=(45,45))
        savefig(joinpath(destfolder, "wireframe_$(expvar)_b_t.png"))

        criteria_with = data.results[1].data[findall(x -> !isnan(x), data.results[1].data[:, bcol]), r2col]
        criteria_without = data.results[1].data[findall(x -> isnan(x), data.results[1].data[:, bcol]), r2col]

        uniden_with = kde(criteria_with)
        uniden_without = kde(criteria_without)

        p1 = plot(range(min(criteria_with...), stop = max(criteria_with...), length = 150), z -> pdf(uniden_with, z))
        p1 = plot!(range(min(criteria_without...), stop = max(criteria_without...), length = 150), z -> pdf(uniden_without,z))
        plot(p1, label=["Including $(expvar)" "Excluding $(expvar)"], ylabel = "Adj. R2")
        savefig(joinpath(destfolder, "Kdensity_criteria_$(expvar).png"))

        p2 = violin(["Including $(expvar)" "Excluding $(expvar)"], [criteria_with, criteria_without], leg = false, marker = (0.1, stroke(0)), alpha = 0.50, color = :blues)
        p2 = boxplot!(["Including $(expvar)" "Excluding $(expvar)"], [criteria_with, criteria_without], leg = false, marker = (0.3, stroke(2)), alpha = 0.6, color = :orange)
        plot(p2, ylabel = "Adj. R2")
        savefig(joinpath(destfolder, "BoxViolinDot_$(expvar).png"))

        criteria_diff[i, 1] = mean(criteria_with) - mean(criteria_without)
        criteria_diff[i, 2] = "$(expvar)"
    end

    a = sortslices(criteria_diff, dims = 1)
    labels = convert(Array{String}, a[:,2])
    bar(labels, a[:,1], legend = false, color = :blues, orientation = :horizontal, xlabel="Average impact of each variable on the Adj. R2")
    savefig(joinpath(destfolder, "cov_relevance.png"))

    positivegainsvariables = findall(x -> x > 0, a[:,1])
    negativegainsvariables = findall(x -> x < 0, a[:,1])

    intelligent_text = Dict()
    intelligent_text["numofpositivegainsvariables"] = size(positivegainsvariables, 1)
    intelligent_text["numofnegativegainsvariables"] = size(negativegainsvariables, 1)

    intelligent_text["multiple_positive_variables"] = intelligent_text["numofpositivegainsvariables"] > 1
    intelligent_text["multiple_negative_variables"] = intelligent_text["numofpositivegainsvariables"] > 1
    intelligent_text["no_negative_variables"] = intelligent_text["numofnegativegainsvariables"] == 0

    intelligent_text["bestvar"] = a[1,2]
    intelligent_text["bestvar_gainsinperc"] = a[1,1]

    intelligent_text["worstvar"] = a[end,2]
    intelligent_text["worstvar_gainsinperc"] = a[end,1]

    intelligent_text
end


text:
"""
Initialize harmonic external potential for electronic system.
"""
function init_pot_harm_3d( LF::LF3dGrid, ω, center )
  Npoints = LF.Nx * LF.Ny * LF.Nz
  Vpot = zeros( Float64,Npoints )
  for ip = 1:Npoints
    Vpot[ip] = 0.5 * ω^2 * norm(LF.lingrid[:,ip] - center[:])^2
  end
  return Vpot
end

"""
Initialize ion-electron potential according to pure Coulombic
potential.
This potential is singular at r = 0.
"""
function init_pot_Hcoul( LF::LF3dGrid, center )
  Npoints = LF.Nx * LF.Ny * LF.Nz
  Vpot = zeros( Float64,Npoints )
  for ip = 1:Npoints
    r = norm( LF.lingrid[:,ip] - center[:] )
    Vpot[ip] = -1.0/r
  end
  return Vpot
end

"""
Initialize ion-electron potential for hydrogen atom using
analytic form and parameters used in Gygi's paper.
(need reference)
"""
function init_pot_Hps( LF::LF3dGrid, center )
  Npoints = LF.Nx * LF.Ny * LF.Nz
  Vpot = zeros( Float64,Npoints )

  const rc1 = 0.25
  const rc2 = 0.284
  const aa  = -1.9287
  const bb  = 0.3374

  for ip=1:Npoints
    r = norm( LF.lingrid[:,ip] - center[:] )
    Vpot[ip] = -1.0/r * erf( r/rc1 ) + (aa + bb*r^2)*exp(-(r/rc2)^2)
  end
  return Vpot
end

"""
Initialize local ion-electron potential for hydrogen atom according
to HGH pseudopotential parameters.

This function is used only for testing purpose.

This function should not be used to initialize ion-electron potential for
periodic system.
"""
function init_pot_Hps_HGH( LF::LF3dGrid, center )
  Npoints = LF.Nx * LF.Ny * LF.Nz
  Vpot = zeros( Float64,Npoints )

  const Zion = 1
  const rloc = 0.2
  const C1 = -4.0663326
  const C2 = 0.6678322

  # TODO Add journal reference
  for ip = 1:Npoints
    r = norm( LF.lingrid[:,ip] - center[:] )
    if r < 1e-7
      println("WARNING: small r:", r);
      println("Using limiting value:")
      Vpot[ip] = -2*Zion/(sqrt(2*pi)*rloc) + C1
    else
      rrloc = r/rloc
      Vpot[ip] = -Zion/r * erf( r/(sqrt(2.0)*rloc) ) +
               (C1 + C2*rrloc^2)*exp(-0.5*(rrloc)^2)
    end
  end
  return Vpot
end