text:
"""
gwasvctest(plinkFile::String = "", covFile::String = "", traitFile::String = "",
annotationFile::String = "", outFile::String = "", test::String = "eRLRT",
kinship::String = "GRM", kernel::String = "GRM", pvalueComputing::String = "chi2",
windowSize::Int = 50, infLambda::Float64 = 0.0, MemoryLimit::Int = 200000000,
nNullSimPts::Int = 1000000, nMMmax::Int = 0, nBlockAscent::Int = 1000,
nNullSimNewtonIter::Int = 15, tolX::Float64 = 1e-4, snpWtType::String = "")

Test the SNP set effect with presence of additive genetic effect and
enviromental effect.

Optional input name-value pairs:

"plinkFile" - Input Plink files names without extension.

"covFile" - Input file to form the covariate matrix X.

"traitFile" - Input file to form the response vector y.

"annotationFile" - Input file for the annotation of SNPs.

"test"- "eLRT"|"eRLRT"|"eScore"|"none", the requested test, it also
dictates the estimation method, default is "eRLRT".

"kinship" - "GRM"|"MoM"|"theoretical"|"none", indicates method used to
compute the kinship matrix, default is "GRM". You can also use your
own kinship matrix, you can do that by replacing "none" with the file
name in which store the pre-calculated kinship matrix.

"kernel" - "GRM"|"IBS1"|"IBS2"|"IBS3", indicates method used to compute
the kernel matrix S, default is "GRM".

"pvalueComputing" - "MonteCarlo"|"chi2", indicates the method (Monte
Carlo simulaion or Chi square approximation) used to compute the
p-value for eLRT or eRLRT, default is "chi2".

"windowSize" - # SNPs in a SNP set, default is 50.

"infLambda" - Any rational number with float precision. Control the
trade-off for approximating kinship matrix, default is 0.0.

"MemoryLimit" - Max memory used to read and store the raw data, default
is 200000000 (equivalent to 200 megabytes).

"nNullSimPts"- # simulation samples, default is 1,000,000.

"nMMmax" - Max MM iterations, default is 10 (eLRT) or 1000 (eRLRT).

"nBlockAscent" - Max block ascent iterations, default is 1000.

"nNullSimNewtonIter" - Max Newton iteration, default is 15.

"tolX" - Tolerance in change of parameters, default is 1e-4.

"snpWtType" - "beta"|"invvar", add weights to the kernel matrix S,
default is "" (add no weights by default).

Output name-value pair:

"outFile" - Output file which store the P-values for SNP sets effects.
"""
function gwasvctest(args...; covFile::String = "", device::String = "CPU",
                    infLambda::Float64 = 0.0, kernel::String = "GRM",
                    kinship::String = "GRM", nMMmax::Int = 0,
                    nBlockAscent::Int = 1000, nNullSimPts::Int = 1000000,
                    nNullSimNewtonIter::Int = 15, outFile::String = "",
                    plinkFile::String = "", snpWtType::String = "",
                    test::String = "eRLRT", tolX::Float64 = 1e-4,
                    traitFile::String = "", MemoryLimit::Int = 200000000,
                    pvalueComputing::String = "chi2", windowSize::Int = 50,
                    annotationFile::String = "")

  ## parse data from files
  plinkBedfile = string(plinkFile, ".bed");
  plinkBimfile = string(plinkFile, ".bim");
  plinkFamfile = string(plinkFile, ".fam");

  # BIM file: chr, rs#, morgan, bp position, allele1, allele2
  bimdata = readdlm(plinkBimfile, String);
  chrID = bimdata[:, 1];
  snpID = bimdata[:, 2];
  nSNP = length(snpID);
  snpPos = bimdata[:, 4];

  # FAM file: fam ID, ind ID, father ID, mother ID, sex, phenotype
  famdata = readdlm(plinkFamfile);
  nPer = length(famdata[:, 1]);

  # prepare for reading binary data
  map2geno = [2; NaN; 1; 0];
  bin2geno = Array(Float64, 4, 2 ^ 8);
  countbin = 1;
  for iB4 = 0 : 3
    tmpper4 = map2geno[iB4 + 1];
    for iB3 = 0 : 3
      tmpper3 = map2geno[iB3 + 1];
      for iB2 = 0 : 3
        tmpper2 = map2geno[iB2 + 1];
        for iB1 = 0 : 3
          tmpper1 = map2geno[iB1 + 1];
          bin2geno[1, countbin] = tmpper1;
          bin2geno[2, countbin] = tmpper2;
          bin2geno[3, countbin] = tmpper3;
          bin2geno[4, countbin] = tmpper4;
          countbin += 1;
        end
      end
    end
  end

  # BED file:
  fid = open(plinkBedfile);
  bedheader = mmap_array(Uint8, (3, 1), fid);
  # decide BED file version and snp/individual-major
  # PLINK coding (bits->genotype): 00->1/1, 01->1/2, 10->Missing, 11->2/2
  # We use minor allele counts (0,1,2,nan) coding here
  # geno = 2 => NaN
  # geno = 1 => 1
  # geno = 0 => 2
  # geno = 3 => 0
  if bits(bedheader[1, 1]) == "01101100" && bits(bedheader[2, 1]) == "00011011"
    # v1.0 BED file
    if bits(bedheader[3, 1]) == "00000001"
      # SNP-major
      #rawdata = mmap_array(Uint8, (iceil(nPer / 4), nSNP), fid, 3);
      rawdata = mmap_array(Int8, (iceil(nPer / 4), nSNP), fid, 3);
      SNPSize = ifloor(MemoryLimit / (32 * iceil(nPer / 4)));
    else
      # individual-major
      error("gwasvctest:bedwrongn\n",
            "# Individual-major BED file found!",
            "Please transform to SNP-major BED file using PLINK");
    end
  else
    # v0.99 BED file: individual-major
    error("gwasvctest:bedwrongn\n",
          "# Individual-major BED file found!",
          "Please transform to SNP-major BED file using PLINK");
  end
  close(fid);
  # number of read times
  nReads = iceil(nSNP / SNPSize);

  # check if provide SNP annotation file
  if isempty(annotationFile)
    flagAnnotate = false;
  else
    flagAnnotate = true;
    grpInfo = zeros(Int64, nSNP);
    offsetSize = zeros(Int64, nSNP);
    geneName = Array(String, nSNP);
    readAnnotate!(annotationFile, snpID, nSNP, grpInfo, offsetSize, geneName);
    grpInfo = grpInfo[grpInfo .> 0];
    nGrp = length(grpInfo);
    offsetSize = offsetSize[1 : nGrp];
    geneName = geneName[1 : nGrp];
  end

  ## determine the additive genetic component

  if kinship == "GRM" || kinship == "MoM"

    genoOri = Array(Float64, nPer, SNPSize);
    geno = Array(Float64, nPer, SNPSize);
    kinMat = zeros(nPer, nPer);
    tmpscalar = 0.0;
    nSNPred = 0;
    offset = 0;
    sumIsnan = Array(Float64, 1, SNPSize);
    nChrObs = Array(Float64, 1, SNPSize);
    mafOri = Array(Float64, SNPSize);
    maf = Array(Float64, SNPSize);
    oneConst = ones(nPer);
    tmpvecKernel = Array(Float64, SNPSize);

    for idxRead = 1 : nReads

      # read the binary data
      if idxRead == nReads
        curSNPSize = nSNP - (idxRead - 1) * SNPSize;
      else
        curSNPSize = SNPSize;
      end

      readgeno!(genoOri, curSNPSize, nPer, SNPSize, idxRead, bin2geno, rawdata,
                offset, false);

      tmpMatIsnan = isnan(genoOri);
      sum!(sumIsnan, tmpMatIsnan);
      for i = 1 : curSNPSize
        nChrObs[1, i] = 2 * (nPer - sumIsnan[1, i]);
      end

      # minor allele frequencies and get rid of mono-allelic SNPs
      counter = 0;
      offset = (idxRead - 1) * SNPSize;
      for i = 1 : curSNPSize
        mafOri[i] = 0.0;
        for j = 1 : nPer
          if !tmpMatIsnan[j, i]
            mafOri[i] += genoOri[j, i];
          end
        end
        mafOri[i] = mafOri[i] / nChrObs[1, i];
        if mafOri[i] != 0
          counter += 1;
          maf[counter] = mafOri[i];
          pGenoOri = pointer(genoOri) + (i - 1) * nPer * sizeof(Float64);
          pGeno = pointer(geno) + (counter - 1) * nPer * sizeof(Float64);
          BLAS.blascopy!(nPer, pGenoOri, 1, pGeno, 1);
        end
      end
      nSNPred += (curSNPSize - counter);

      # impute missing genotype by expected minor allele counts
      for j = 1 : counter
        for i = 1 : nPer
          if isnan(geno[i, j])
            geno[i, j] = 2 * maf[j];
          end
        end
      end

      if kinship == "GRM"
        # center and scale genotype matrix by minor allele frequencies
        gMAF = maf[1 : counter];
        if counter < SNPSize
          geno = geno[:, 1 : counter];
          tmpvecKernel = tmpvecKernel[1 : counter];
        end
        BLAS.ger!(-2.0, oneConst, gMAF, geno);
        for i = 1 : counter
          tmpvecKernel[i] = 1.0 / sqrt(2 * gMAF[i] * (1 - gMAF[i]));
        end
        scale!(geno, tmpvecKernel);
        # estimate kinship by genetic relation matrix (GRM) from dense markers
        # TODO: Check Yang et al paper the precise formula
        BLAS.gemm!('N', 'T', 1.0, geno, geno, 1.0, kinMat);
        if counter != SNPSize
          geno = [geno Array(Float64, nPer, SNPSize - counter)];
          tmpvecKernel = [tmpvecKernel; Array(Float64, SNPSize - counter)];
        end
      elseif kinship == "MoM"
        if counter < SNPSize
          geno = geno[:, 1 : counter];
        end
        # shift to {-1,0,1} encoding
        geno = geno - 1.0;
        # estimate kinship by method of moment from dense markers
        gMAF = maf[1 : counter];
        tmpscalar += sumabs2(gMAF) + sumabs2(1 - gMAF);
        BLAS.gemm!('N', 'T', 0.5, geno, geno, 1.0, kinMat);
        if counter != SNPSize
          geno = [geno Array(Float64, nPer, SNPSize - counter)];
        end
      end

    end

    if kinship == "GRM"
      scale!(kinMat, 1.0 / (nSNP - nSNPred));
    elseif kinship == "MoM"
      kinMat += (nSNP - nSNPred) / 2.0 - tmpscalar;
      scale!(kinMat, 1.0 / (nSNP - nSNPred - tmpscalar));
    end

  elseif kinship == "theoretical"

    # use theoretical kinship computed from pedigree information
    # use the kinship function provided by Lange
    IndID = famdata[:, 2];
    AugIndID = unique([famdata[:, 2]; famdata[:, 3]; famdata[:, 4]]);
    AugIndID = AugIndID[AugIndID .!= 0];
    AugIndID = sort(AugIndID);
    keepIdx = findin(AugIndID, IndID);
    vals = [1 : length(AugIndID)];
    IndDict = Dict(AugIndID, vals);
    faID = famdata[:, 3];
    moID = famdata[:, 4];
    father = zeros(Int64, length(AugIndID));
    mother = zeros(Int64, length(AugIndID));
    for i = 1 : nPer
      if faID[i] != 0
        father[keepIdx[i]] = IndDict[faID[i]];
        mother[keepIdx[i]] = IndDict[moID[i]];
      end
    end
    kinMat = kinshipcoef(father, mother, 1, length(AugIndID));
    kinMat = kinMat[keepIdx, keepIdx];

  elseif kinship != "none"

    # user provided kinship
    kinMat = readdlm(kinship, ',');
    if size(kinMat, 1) != nPer || size(kinMat, 2) != nPer
      error("gwasvctest:kinshipwrongn\n",
            "# individuals in kinship file does not match plink files");
    end

  end

  ## get responses
  if isempty(traitFile)
    # no trait file provided: retrieve from plink files
    y = famdata[:, 6];
    y = convert(Array{Float64, 1}, y);
  else
    # user supplied trait values (only take the first column) in Plink
    # phenotype file (no headerline) format: FamID IndID Trait1 Trait2 ...
    # Only Trait1 will be read
    y = readdlm(traitFile);
    if size(y, 1) != nPer
      error("gwasvctest:ywrongn\n",
            "# individuals in trait file does not match plink files");
    end
    y = y[:, 3];
    y = convert(Array{Float64, 1}, y);
  end

  ## get covariates
  if isempty(covFile)
    # no covariates provided
    X = ones(nPer, 1);
  else
    # user provided covariates (intercept if exists has to be provided)
    X = readdlm(covFile);
    if size(X, 1) != nPer
      error("gwasvctest:covwrongn\n",
            "# individuals in covariate file does not match plink files");
    end
    X = [ones(nPer, 1) X[:, 3:end]];
    X = convert(Array{Float64, 2}, X);
  end

  ## individuals with misisng phenotype are removed from analysis
  ## missing covariate values are imputed by avarages

  keepIdx = !isnan(y);
  nPerKeep = countnz(keepIdx);
  y = y[keepIdx];
  X = X[keepIdx, :];
  meanX = zeros(size(X, 2));
  for i = 1:length(meanX)
    meanX[i] = mean(X[:, i][!isnan(X[:, i])]);
  end
  ridx, cidx = ind2sub(size(X), find(isnan(X)));
  X[sub2ind(size(X), find(isnan(X)))] = meanX[cidx];
  if kinship != "none"
    kinMat = kinMat[keepIdx, keepIdx];
  end

  ## pre-processing for testing

  # if provide annotation file, avoid large gene set
  if flagAnnotate
    numExceed = countnz(grpInfo .> ifloor(nPerKeep / 4));
    for i = 1 : numExceed
      tmpidx = find(grpInfo .> ifloor(nPerKeep / 4))[1];
      innerItr = iceil(grpInfo[tmpidx] / windowSize);
      tmpgrpInfo = zeros(Int64, innerItr);
      tmpoffsetSize = zeros(Int64, innerItr);
      for j = 1 : innerItr
        if j < innerItr
          tmpgrpInfo[j] = windowSize;
          tmpoffsetSize[j] = offsetSize[tmpidx] + (j - 1) * windowSize;
        else
          tmpgrpInfo[j] = grpInfo[tmpidx] - (j - 1) * windowSize;
          tmpoffsetSize[j] = offsetSize[tmpidx] + (j - 1) * windowSize;
        end
      end
      grpInfo = [grpInfo[1:tmpidx-1]; tmpgrpInfo; grpInfo[tmpidx+1:end]];
      offsetSize = [offsetSize[1:tmpidx-1]; tmpoffsetSize; offsetSize[tmpidx+1:end]];
      geneName = vcat(geneName[1:tmpidx-1], repmat([geneName[tmpidx]], innerItr, 1),
                      geneName[tmpidx+1:end]);
    end

    windowSize = maximum(grpInfo);
    nReads = length(grpInfo);
  end

  # simulate standard normals and square it
  WPreSim = randn(windowSize, nNullSimPts);
  for idxWc = 1 : nNullSimPts
    for idxWr = 1 : windowSize
      WPreSim[idxWr, idxWc] = WPreSim[idxWr, idxWc] * WPreSim[idxWr, idxWc];
    end
  end

  if kinship == "none"
    # only one non-trivial variance component: SNP set

    # pre-compute basis of N(X') and projected y
    Xsvd = svdfact(X, thin = false);
    rankX = countnz(Xsvd[:S] .> nPerKeep * eps(Xsvd[:S][1]));
    XtNullBasis = Xsvd[:U][:, rankX + 1 : end]';
    ynew = zeros(nPerKeep - rankX);
    BLAS.gemv!('N', 1.0, XtNullBasis, y, 0.0, ynew);

    rankQPhi = 0;
    evalPhiAdj = Float64[];
    XPhitNullBasis = [Float64[] Float64[]];
    yShift = Float64[];
    KPhiAdj = [Float64[] Float64[]];
    weightedW = [Float64[] Float64[]];
    QPhi = [Float64[] Float64[]];
    PrePartialSumW = [Float64[] Float64[]];
    PreTotalSumW = [Float64[] Float64[]];
    nPreRank = min(nPerKeep, windowSize);
    tmpn = nPerKeep - rankX;

  else
    # >=2 non-trivial variance components

    ynew = Float64[];

    # obtain a basis of X and (possibly) extra variance components
    # Xsvd[:U] -- QX
    Xsvd = svdfact(X, thin = false);
    rankX = countnz(Xsvd[:S] .> nPerKeep * eps(Xsvd[:S][1]));
    XtNullBasis = Xsvd[:U][:, rankX + 1 : end];
    QX = Xsvd[:U][:, 1 : rankX];

    # eigen-decomposition of the additive component
    # Phieig[:vectors] -- evecPhi
    # Phieig[:values] -- evalPhi
    Phieig = eigfact(kinMat);
    rankPhi = countnz(Phieig[:values] .>
                      max(infLambda, nPerKeep * eps(maximum(Phieig[:values]))));
    #evecPhi = Phieig[:vectors];
    #evalPhi = Phieig[:values];
    # enforce first entry of each eigenvector to be >=0
    # for cross-platform reproducibility
    idxEvecPhi = vec(Phieig[:vectors][1, :]) .< 0;
    #Phieig[:vectors][:, idxEvecPhi] = - Phieig[:vectors][:, idxEvecPhi];
    scale!(Phieig[:vectors][:, idxEvecPhi], -1.0);
    #evecPhi[:, idxEvecPhi] = - evecPhi[:, idxEvecPhi];

    # determine a rank for additive component and a rank for snp set S
    if kernel == "GRM"
      rankS = windowSize;
    elseif kernel == "IBS1"
      rankS = windowSize + 1;
    elseif kernel == "IBS2"
      rankS = 2 * windowSize;
    elseif kernel == "IBS3"
      rankS = 2 * windowSize + 1;
    end
    rankPhi = min(rankPhi, ifloor((nPerKeep - rankX - rankS) / 2));

    # low rank approximation
    sortIdx = sortperm(Phieig[:values], rev = true);
    evalPhi = Phieig[:values][sortIdx[1 : rankPhi]];
    evecPhi = Phieig[:vectors][:, sortIdx[1 : rankPhi]];

    # obtain orthonormal basis of (additive component - covariates)
    # See Golub and Van Loan (1996) Algorithm 12.4.3 & Theorem 12.4.2 on p604
    QPhisvd = svdfact(BLAS.gemm('T', 'N', XtNullBasis, evecPhi));
    PhiKeepIdx = (1 - QPhisvd[:S]) .< 1e-6;
    rankQPhi = countnz(PhiKeepIdx);
    #QPhi = evecPhi * QPhisvd[:V][:, 1 : rankQPhi];
    QPhi = Array(Float64, nPerKeep, rankQPhi);
    BLAS.gemm!('N', 'N', 1.0, evecPhi, QPhisvd[:V][:, 1 : rankQPhi], 0.0, QPhi);
    XPhitNullBasis = null([QX evecPhi]');

    # obtain the square root of the rotated Phi
    # PhiAdjsvd[:U] -- W
    # PhiAdjsvd[:S] -- evalPhiAdj
    tmpMat = Array(Float64, rankQPhi, rankPhi);
    #println("rankPhi: ", rankPhi);
    #println("col of evecPhi: ", size(evecPhi, 2));
    BLAS.gemm!('T', 'N', 1.0, QPhi, evecPhi, 0.0, tmpMat);
    scale!(tmpMat, sqrt(evalPhi));
    PhiAdjsvd = svdfact(tmpMat, thin = false);
    #PhiAdjsvd = svdfact(BLAS.gemm('T', 'N', QPhi, evecPhi) .*
    #                    sqrt(evalPhi)', thin = false);
    #W = PhiAdjsvd[:U];
    evalPhiAdj = PhiAdjsvd[:S] .^ 2;
    # enforce first entry of each eigenvector to be >=0
    idxW = vec(PhiAdjsvd[:U][1, :]) .< 0;
    #PhiAdjsvd[:U][:, idxW] = - PhiAdjsvd[:U][:, idxW];
    scale!(PhiAdjsvd[:U][:, idxW], -1.0);
    #KPhiAdj = PhiAdjsvd[:U] .* sqrt(evalPhiAdj' / minimum(evalPhiAdj) - 1);
    KPhiAdj = PhiAdjsvd[:U][:, :];
    scale!(KPhiAdj, sqrt(evalPhiAdj / minimum(evalPhiAdj) - 1));
    InvSqrtevalPhiAdj = similar(evalPhiAdj);
    for i = 1 : length(evalPhiAdj)
      InvSqrtevalPhiAdj[i] = 1 / sqrt(evalPhiAdj[i]);
    end
    weightedW = PhiAdjsvd[:U][:, :];
    scale!(weightedW, InvSqrtevalPhiAdj);

    # precompute shift in Y
    yShift = QPhi' * y;

    # prepare for simulating Chi Squares and sums of Chi Squares
    nPreRank = min(length(evalPhiAdj), windowSize);
    tmpn = length(evalPhiAdj);

  end

  # simulate Chi Squares and sums of Chi Squares
  PrePartialSumW = Array(Float64, nNullSimPts, nPreRank+1);
  PreTotalSumW = Array(Float64, nNullSimPts, nPreRank+1);
  tmpSumVec = Array(Float64, nNullSimPts);
  p1 = pointer(PrePartialSumW);
  BLAS.blascopy!(nNullSimPts, rand(Distributions.Chisq(tmpn), nNullSimPts), 1, p1, 1);
  for i = 1 : nNullSimPts
    #pW = pointer(WPreSim) + (i - 1) * windowSize * sizeof(Float64);
    #tmpSumVec[i] = BLAS.asum(1, pW, 1);
    tmpSumVec[i] = 0.0;
    PreTotalSumW[i, 1] = tmpSumVec[i] + PrePartialSumW[i, 1];
  end
  for j = 1 : nPreRank
    p1 = pointer(PrePartialSumW) + j * nNullSimPts * sizeof(Float64);
    BLAS.blascopy!(nNullSimPts, rand(Distributions.Chisq(tmpn - j), nNullSimPts), 1, p1, 1);
    for i = 1 : nNullSimPts
      tmpSumVec[i] += WPreSim[j, i];
      PreTotalSumW[i, j+1] = tmpSumVec[i] + PrePartialSumW[i, j+1];
    end
  end

  # prepare output
  if isempty(outFile)
    outFile = string(plinkFile, "-julia.out");
  end
  fid = open(outFile, "w");
  if !flagAnnotate
    println(fid, "StartSNP,EndSNP,pvalue");
  else
    println(fid, "StartSNP,EndSNP,Chr,StartPos,EndPos,GeneName,#SNPs,pvalue");
  end
  close(fid);

  #### test group by group

  ## no annotation file
  if !flagAnnotate

    # preparation for parallel
    if nprocs() == 1
      ncores = 1;
    else
      ncores = nprocs() - 1;
    end
    nSNPParallel = Array(Int, ncores, 1);
    nReadsParallel = Array(Int, ncores, 1);
    rawdataParallel = Array(Matrix{Int8}, ncores, 1);
    snpIDParallel = Array(Vector{String}, ncores, 1);
    fill!(nSNPParallel, int(nSNP / ncores));
    nSNPParallel[ncores] = nSNP - sum(nSNPParallel[1 : ncores-1]);
    totaloffset = 0;
    for i = 1 : ncores
      if i > 1
        totaloffset += nSNPParallel[i-1];
      end
      nReadsParallel[i] = iceil(nSNPParallel[i] / SNPSize);
      rawdataParallel[i] =
        rawdata[:, totaloffset+1 : totaloffset+nSNPParallel[i]];
      snpIDParallel[i] = snpID[totaloffset+1 : totaloffset+nSNPParallel[i]];
    end
    # loop over windows
    outresults = pmap(loopwinFixsize, nSNPParallel, nReadsParallel,
                      rawdataParallel, snpIDParallel, {bin2geno for i = 1:ncores},
                      {nPer for i = 1:ncores}, {SNPSize for i = 1:ncores},
                      {flagAnnotate for i = 1:ncores},
                      {snpWtType for i = 1:ncores}, {keepIdx for i = 1:ncores},
                      {nPerKeep for i = 1:ncores}, {kernel for i = 1:ncores},
                      {kinship for i = 1:ncores}, {test for i = 1:ncores},
                      {XtNullBasis for i = 1:ncores}, {WPreSim for i = 1:ncores},
                      {nBlockAscent for i = 1:ncores}, {nMMmax for i = 1:ncores},
                      {nNullSimPts for i = 1:ncores}, {tolX for i = 1:ncores},
                      {pvalueComputing for i = 1:ncores},
                      {device for i = 1:ncores},
                      {nNullSimNewtonIter for i = 1:ncores},
                      {rankQPhi for i = 1:ncores}, {evalPhiAdj for i = 1:ncores},
                      {XPhitNullBasis for i = 1:ncores},
                      {yShift for i = 1:ncores}, {y for i = 1:ncores},
                      {KPhiAdj for i = 1:ncores}, {weightedW for i = 1:ncores},
                      {QPhi for i = 1:ncores}, {PrePartialSumW for i = 1:ncores},
                      {PreTotalSumW for i = 1:ncores},
                      {windowSize for i = 1:ncores},
                      {nPreRank for i = 1:ncores}, {X for i = 1:ncores},
                      {ynew for i = 1:ncores});

    # write output
    fid = open(outFile, "a");
    for i = 1 : ncores
      for j = 1 : size(outresults[i], 1)
        println(fid, outresults[i][j, 1], ",", outresults[i][j, 2], ",",
                outresults[i][j, 3]);
      end
    end
    close(fid);

  ## provide annotation file
  else

    # preparation for parallel
    if nprocs() == 1
      ncores = 1;
    else
      ncores = nprocs() - 1;
    end
    nReadsParallel = Array(Int, ncores, 1);
    rawdataParallel = Array(Matrix{Int8}, ncores, 1);
    grpInfoParallel = Array(Vector{Int}, ncores, 1);
    offsetSizeParallel = Array(Vector{Int}, ncores, 1);
    snpIDParallel = Array(Vector{String}, ncores, 1);
    chrIDParallel = Array(Vector{String}, ncores, 1);
    snpPosParallel = Array(Vector{String}, ncores, 1);
    geneNameParallel = Array(Vector{String}, ncores, 1);
    fill!(nReadsParallel, int(nReads / ncores));
    nReadsParallel[ncores] = nReads - sum(nReadsParallel[1 : ncores-1]);
    totaloffset = 0;
    for i = 1 : ncores
      if i == 1
        tmpoffset = 0;
      else
        tmpoffset = sum(nReadsParallel[1:i-1]);
        totaloffset += sum(grpInfoParallel[i-1]);
      end
      grpInfoParallel[i] = grpInfo[tmpoffset+1 : tmpoffset+nReadsParallel[i]];
      geneNameParallel[i] = geneName[tmpoffset+1 : tmpoffset+nReadsParallel[i]];
      offsetSizeParallel[i] =
        offsetSize[tmpoffset+1 : tmpoffset+nReadsParallel[i]] -
        offsetSize[tmpoffset+1];
      rawdataParallel[i] =
        rawdata[:, totaloffset+1 : totaloffset+sum(grpInfoParallel[i])];
      snpIDParallel[i] =
        snpID[totaloffset+1 : totaloffset+sum(grpInfoParallel[i])];
      chrIDParallel[i] =
        chrID[totaloffset+1 : totaloffset+sum(grpInfoParallel[i])];
      snpPosParallel[i] =
        snpPos[totaloffset+1 : totaloffset+sum(grpInfoParallel[i])];
    end

    # loop over genes
    outresults = pmap(loopwinAnnot, nReadsParallel, rawdataParallel,
                      grpInfoParallel, offsetSizeParallel, snpIDParallel,
                      {bin2geno for i = 1:ncores}, {nPer for i = 1:ncores},
                      {SNPSize for i = 1:ncores}, {flagAnnotate for i = 1:ncores},
                      {snpWtType for i = 1:ncores}, {keepIdx for i = 1:ncores},
                      {nPerKeep for i = 1:ncores}, {kernel for i = 1:ncores},
                      {kinship for i = 1:ncores}, {test for i = 1:ncores},
                      {XtNullBasis for i = 1:ncores}, {WPreSim for i = 1:ncores},
                      {nBlockAscent for i = 1:ncores}, {nMMmax for i = 1:ncores},
                      {nNullSimPts for i = 1:ncores}, {tolX for i = 1:ncores},
                      {pvalueComputing for i = 1:ncores},
                      {device for i = 1:ncores},
                      {nNullSimNewtonIter for i = 1:ncores},
                      {rankQPhi for i = 1:ncores}, {evalPhiAdj for i = 1:ncores},
                      {XPhitNullBasis for i = 1:ncores},
                      {yShift for i = 1:ncores}, {y for i = 1:ncores},
                      {KPhiAdj for i = 1:ncores}, {weightedW for i = 1:ncores},
                      {QPhi for i = 1:ncores}, {PrePartialSumW for i = 1:ncores},
                      {PreTotalSumW for i = 1:ncores},
                      {windowSize for i = 1:ncores},
                      {nPreRank for i = 1:ncores}, {X for i = 1:ncores},
                      snpPosParallel, chrIDParallel, geneNameParallel,
                      {ynew for i = 1:ncores});

    # write output
    fid = open(outFile, "a");
    for i = 1 : ncores
      for j = 1 : size(outresults[i], 1)
        println(fid, outresults[i][j, 1], ",", outresults[i][j, 2], ",",
                outresults[i][j, 3], ",", outresults[i][j, 4], ",",
                outresults[i][j, 5], ",", outresults[i][j, 6], ",",
                outresults[i][j, 7], ",", outresults[i][j, 8]);
      end
    end
    close(fid);

  end

end


text:
function parse_cmdline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "mcmc"
            help = "Run MCMC to optimize hyperparameters"
            action = :command
        "predict"
            help = "Calculate the posterior of a GP at tdata given data"
            action = :command
        "sample"
            help = "Sample the posterior of a GP"
            action = :command
        "fitplot"
            help = "Diagnostic plots showing the posteriors of different components of the GP"
            action = :command
        "select"
            help = "Output model selection parameters; requires --mcmc and --mcmc2"
            action = :command
        "--log"
            help = "Direct logging output to a file"
        "--debug"
            help = "Enable debug logging"
            action = :store_true
    end

    @add_arg_table! s["mcmc"] begin
        "formula"
            help = "Gaussian Process formula"
            required = true
        "--data", "-x"
            help = "Input data files, accepts \"stdin\". ;-separated, use : to provide additional flags which can be combined: \"#:\" transposes the table, \",:\" reads as CSV, \"~:\" reads as TSV (default). Other characters before : give the column/row to join the tables with, e.g. id:data.tsv;#subjectid:subjects.tsv will use the id column of data.tsv and subjectid row of subjects.tsv."
            required = true
        "--rmv_outliers"
            help = "Outlier removal method for training data (none|fence)"
            default = "none"
        "--outlier_fields"
            help = ";-separated list of additional fields to include in outlier removal"
            default = ""
        "--outlier_ignore"
            help = ";-separated list of fields to ignore for outlier removal"
            default = ""
        "--mcmc", "-m"
            help = "MCMC samples for hyperparameters; if provided, the chain will be extended. Does NOT support --data format flags"
        "--burnin", "-r"
            help = "Burn-in samples"
            arg_type = Int
            default = 0
        "--thin", "-t"
            help = "Thinning"
            arg_type = Int
            default = 1
        "--samples", "-n"
            help = "Number of MCMC samples to generate (after anneal, burn-in, and thinning)"
            arg_type = Int
            default = 100
        "--output", "-o"
            help = "Output filename, accepts \"stdout\", does NOT support --data format flags"
            default = "stdout"
    end

    @add_arg_table! s["predict"] begin
        "formula"
            help = "Gaussian Process formula"
            required = true
        "--data", "-x"
            help = "Input data files, accepts \"stdin\". ;-separated, use : to provide additional flags which can be combined: \"#:\" transposes the table, \",:\" reads as CSV, \"~:\" reads as TSV (default). Other characters before : give the column/row to join the tables with, e.g. id:data.tsv;#subjectid:subjects.tsv will use the id column of data.tsv and subjectid row of subjects.tsv."
            required = true
        "--rmv_outliers"
            help = "Outlier removal method for training data (none|fence)"
            default = "none"
        "--outlier_fields"
            help = ";-separated list of additional fields to include in outlier removal"
            default = ""
        "--outlier_ignore"
            help = ";-separated list of fields to ignore for outlier removal"
            default = ""
        "--mcmc", "-m"
            help = "MCMC samples for hyperparameters, does NOT support --data format flags"
        "--atdata", "-t"
            help = "Data files providing variable values at which to make predictions"
        "--at"
            help = "Predict at these points (alternative to tdata); format: ;-separated variable=start:step:end OR variable=value OR variable=Julia distribution OR variable/group=value"
        "--output", "-o"
            help = "Output filename, accepts \"stdout\", supports --data format flags"
            default = "stdout"
    end

    @add_arg_table! s["sample"] begin
        "formula"
            help = "Gaussian Process formula"
            required = true
        "--data", "-x"
            help = "Input data files, accepts \"stdin\". ;-separated, use : to provide additional flags which can be combined: \"#:\" transposes the table, \",:\" reads as CSV, \"~:\" reads as TSV (default). Other characters before : give the column/row to join the tables with, e.g. id:data.tsv;#subjectid:subjects.tsv will use the id column of data.tsv and subjectid row of subjects.tsv."
        "--rmv_outliers"
            help = "Outlier removal method for training data (none|fence)"
            default = "none"
        "--outlier_fields"
            help = ";-separated list of additional fields to include in outlier removal"
            default = ""
        "--outlier_ignore"
            help = ";-separated list of fields to ignore for outlier removal"
            default = ""
        "--mcmc", "-m"
            help = "MCMC samples for hyperparameters, does NOT support --data format flags"
        "--atdata", "-t"
            help = "Data files providing variable values at which to sample the GP"
        "--at"
            help = "Predict at these points (alternative to tdata); format: ;-separated variable=start:step:end OR variable=value OR variable=Julia distribution OR variable/group=value"
        "--output", "-o"
            help = "Output filename, accepts \"stdout\", supports --data format flags"
            default = "stdout"
        "--plot"
            help = "Output filename for the sample plots"
        "--plotx"
            help = "X axis of the sample plot; can also be \"x:group\" to provide a grouping variable as well"
    end

    @add_arg_table! s["fitplot"] begin
        "formula"
            help = "Gaussian Process formula"
            required = true
        "--data", "-x"
            help = "Input data files, accepts \"stdin\". ;-separated, use : to provide additional flags which can be combined: \"#:\" transposes the table, \",:\" reads as CSV, \"~:\" reads as TSV (default). Other characters before : give the column/row to join the tables with, e.g. id:data.tsv;#subjectid:subjects.tsv will use the id column of data.tsv and subjectid row of subjects.tsv."
            required = true
        "--rmv_outliers"
            help = "Outlier removal method for training data (none|fence)"
            default = "none"
        "--outlier_fields"
            help = ";-separated list of additional fields to include in outlier removal"
            default = ""
        "--outlier_ignore"
            help = ";-separated list of fields to ignore for outlier removal"
            default = ""
        "--mcmc", "-m"
            help = "MCMC samples for hyperparameters, does NOT support --data format flags"
        "--component", "--comp"
            help = "Plot fit plots for components with the given variables"
        "--output", "-o"
            help = "Output filename"
            default = "fitplots.pdf"
    end

    @add_arg_table! s["select"] begin
        "--mcmc", "-m"
            help = "MCMC samples for hyperparameters, does NOT support --data format flags"
        "--mcmc2", "-c"
            help = "MCMC samples for hyperparameters of the second model (for model comparison), does NOT support --data format flags"
    end

    return parse_args(s)
end

ext: jl

max_stars_repo_name: biobakery/GaPLAC

max_issues_repo_name: biobakery/GaPLAC

max_forks_repo_name: biobakery/GaPLAC


text:
"""
    evaluate_results(;[simname=""])

Convert simulation results for data exchange with other frameworks.

# Arguments
- `simname`: Path to a specific .jld file to be converted

See also: [Config](@ref), [cfg_sim](@ref), [cfg_ot](@ref), [cfg_ag](@ref), [cfg_feed](@ref), [cfg_mech](@ref)
"""
function evaluate_results(
    ;
    simname::String=""
)

    if !in("dataexchange", readdir())
        mkdir("dataexchange")
    end

    evalqueue = Array{Simulation, 1}[]

    if simname != ""

        raw_data = load(joinpath("results", simname))
        data = raw_data[first(keys(raw_data))]

        push!(evalqueue, data)

    else
        for simname in readdir("results")
            raw_data = load(joinpath("results", simname))
            data = raw_data[first(keys(raw_data))]

            push!(evalqueue, data)
        end
    end

    results = DataFrame(
        RunName = String[],
        RepNr = Int64[],
        Densities = Float64[],
        OutdegreeSD = Float64[],
        IndegreeSD = Float64[],
        OutdegreeIndegreeRatioMean = Float64[],
        ClosenessCentralityMean = Float64[],
        BetweennessCentralityMean = Float64[],
        EigenCentralityMean = Float64[],
        ClustCoeff = Float64[],
        CommunityCount = Int64[],
        CommunityCountNontrivial = Int64[],
        ConnectedComponents = Int64[],
        OpinionSD = Float64[],
        OpChangeDeltaMean = Float64[],
        PublOwnOpinionDiff = Float64[],
        SupernodeOutdegree = Int64[],
        SupernodeCloseness = Float64[],
        SupernodeBetweenness = Float64[],
        SupernodeEigen = Float64[],
        SupernodeOpinion = Float64[],
        Supernode1st2ndOpdiff = Float64[],
        CommunityOpMeanSDs = Float64[]
    )

    for simarray in evalqueue

        repcount = length(simarray)
        nodes_ranked = [first.(sort([(i, outdegree(simarray[j].final_state[1], i)) for i in 1:nv(simarray[j].final_state[1])], by=last, rev=true)) for j in 1:repcount]

        opinionsd = [std([agent.opinion for agent in simarray[i].final_state[2]]) for i in 1:repcount]
        opchange_delta_mean = [mean([abs(simarray[j].init_state[2][i].opinion - simarray[j].final_state[2][i].opinion) for i in 1:simarray[j].config.network.agent_count]) for j in 1:repcount]
        densities = [density(simarray[i].final_state[1]) for i in 1:repcount]
        outdegree_sd = [std(outdegree(simarray[i].final_state[1])) for i in 1:repcount]
        outdegree_mean = [mean(outdegree(simarray[i].final_state[1])) for i in 1:repcount]
        indegree_sd = [std(indegree(simarray[i].final_state[1])) for i in 1:repcount]
        outinratio = [mean([outdegree(simarray[i].final_state[1], j) / degree(simarray[i].final_state[1], j) for j in 1:nv(simarray[i].final_state[1]) if degree(simarray[i].final_state[1], j) > 0]) for i in 1:repcount]
        publownopiniondiff = [mean([abs(agent.perceiv_publ_opinion - agent.opinion) for agent in simarray[i].final_state[2]]) for i in 1:repcount]
        closeness_centrality_mean = [mean(closeness_centrality(simarray[i].final_state[1])) for i in 1:repcount]
        betweenness_centrality_mean = [mean(betweenness_centrality(simarray[i].final_state[1])) for i in 1:repcount]
        eigen_centrality_mean = [mean(eigenvector_centrality(simarray[i].final_state[1])) for i in 1:repcount]
        supernode_outdegree = [outdegree(simarray[i].final_state[1], nodes_ranked[i][1]) for i in 1:repcount]
        supernode_closeness = [closeness_centrality(simarray[i].final_state[1])[nodes_ranked[i][1]] for i in 1:repcount]
        supernode_betweenness = [betweenness_centrality(simarray[i].final_state[1])[nodes_ranked[i][1]] for i in 1:repcount]
        supernode_eigen = [eigenvector_centrality(simarray[i].final_state[1])[nodes_ranked[i][1]] for i in 1:repcount]
        supernode_opinion = [simarray[i].final_state[2][nodes_ranked[i][1]].opinion for i in 1:repcount]
        supernode1st2nd_opdiff = [abs(simarray[i].final_state[2][nodes_ranked[i][1]].opinion - simarray[i].final_state[2][nodes_ranked[i][2]].opinion) for i in 1:repcount]
        clust_coeff = [global_clustering_coefficient(simarray[i].final_state[1]) for i in 1:repcount]
        conn_components = [length(connected_components(simarray[i].final_state[1])) for i in 1:repcount]

        # Calc the SD of the opinion means of the identified clusters
        n_communities = Int64[]
        community_opinion_mean_sds = Float64[]
        n_communities_nontrivial = Int64[]
        for i in 1:repcount
            label_prop = label_propagation(simarray[i].final_state[1])[1]
            community_sizes = last.(collect(countmap(label_prop)))
            push!(n_communities_nontrivial, length([i for i in community_sizes if i > 1]))
            push!(n_communities, maximum(label_prop))
            if maximum(label_prop) == 1
                push!(community_opinion_mean_sds, 0)
                continue
            end
            push!(community_opinion_mean_sds, std([mean([agent.opinion for agent in simarray[i].final_state[2] if agent.id in findall(x->x==j, label_prop)]) for j in 1:maximum(label_prop)]))
        end

        if simarray[1].batch_name != ""
            runname = ["Run$i" for i in 1:length(simarray)]
        else
            runname = simarray[1].name
        end

        append!(
            results,
            DataFrame(
                RunName = runname,
                RepNr = collect(1:repcount),
                Densities = densities,
                OutdegreeSD = outdegree_sd,
                IndegreeSD = indegree_sd,
                OutdegreeIndegreeRatioMean = outinratio,
                ClosenessCentralityMean = closeness_centrality_mean,
                BetweennessCentralityMean = betweenness_centrality_mean,
                EigenCentralityMean = eigen_centrality_mean,
                ClustCoeff = clust_coeff,
                CommunityCount = n_communities,
                CommunityCountNontrivial = n_communities_nontrivial,
                ConnectedComponents = conn_components,
                OpinionSD = opinionsd,
                OpChangeDeltaMean = opchange_delta_mean,
                PublOwnOpinionDiff = publownopiniondiff,
                SupernodeOutdegree = supernode_outdegree,
                SupernodeCloseness = supernode_closeness,
                SupernodeBetweenness = supernode_betweenness,
                SupernodeEigen = supernode_eigen,
                SupernodeOpinion = supernode_opinion,
                Supernode1st2ndOpdiff = supernode1st2nd_opdiff,
                CommunityOpMeanSDs = community_opinion_mean_sds
            )
        )

    end

    if simarray[1].batch_name == ""
        if !in(simarray[1].name, readdir("dataexchange"))
            mkdir(joinpath("dataexchange", simarray[1].name))
        end

        CSV.write(
        joinpath(
        "dataexchange",
        simarray[1].name,
        "results_" * simarray[1].name * ".csv"
        ),
        results
        )

    else
        if !in(simarray[1].batch_name, readdir("dataexchange"))
            mkdir(joinpath("dataexchange", simarray[1].batch_name))
        end

        CSV.write(
        joinpath(
        "dataexchange",
        simarray[1].batch_name,
        "results.csv"
        ),
        results
        )
    end
end


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

mutable struct Bus
    name::String
    offset::Int
    load::Vector{Float64}
    units::Vector
    price_sensitive_loads::Vector
end

mutable struct CostSegment
    mw::Vector{Float64}
    cost::Vector{Float64}
end

mutable struct StartupCategory
    delay::Int
    cost::Float64
end

mutable struct Unit
    name::String
    bus::Bus
    max_power::Vector{Float64}
    min_power::Vector{Float64}
    must_run::Vector{Bool}
    min_power_cost::Vector{Float64}
    cost_segments::Vector{CostSegment}
    min_uptime::Int
    min_downtime::Int
    ramp_up_limit::Float64
    ramp_down_limit::Float64
    startup_limit::Float64
    shutdown_limit::Float64
    initial_status::Union{Int,Nothing}
    initial_power::Union{Float64,Nothing}
    provides_spinning_reserves::Vector{Bool}
    startup_categories::Vector{StartupCategory}
end

mutable struct TransmissionLine
    name::String
    offset::Int
    source::Bus
    target::Bus
    reactance::Float64
    susceptance::Float64
    normal_flow_limit::Vector{Float64}
    emergency_flow_limit::Vector{Float64}
    flow_limit_penalty::Vector{Float64}
end

mutable struct Reserves
    spinning::Vector{Float64}
end

mutable struct Contingency
    name::String
    lines::Vector{TransmissionLine}
    units::Vector{Unit}
end

mutable struct PriceSensitiveLoad
    name::String
    bus::Bus
    demand::Vector{Float64}
    revenue::Vector{Float64}
end

Base.@kwdef mutable struct UnitCommitmentInstance
    buses_by_name::Dict{AbstractString,Bus}
    buses::Vector{Bus}
    contingencies_by_name::Dict{AbstractString,Contingency}
    contingencies::Vector{Contingency}
    lines_by_name::Dict{AbstractString,TransmissionLine}
    lines::Vector{TransmissionLine}
    power_balance_penalty::Vector{Float64}
    price_sensitive_loads_by_name::Dict{AbstractString,PriceSensitiveLoad}
    price_sensitive_loads::Vector{PriceSensitiveLoad}
    reserves::Reserves
    shortfall_penalty::Vector{Float64}
    time::Int
    units_by_name::Dict{AbstractString,Unit}
    units::Vector{Unit}
end

function Base.show(io::IO, instance::UnitCommitmentInstance)
    print(io, "UnitCommitmentInstance(")
    print(io, "$(length(instance.units)) units, ")
    print(io, "$(length(instance.buses)) buses, ")
    print(io, "$(length(instance.lines)) lines, ")
    print(io, "$(length(instance.contingencies)) contingencies, ")
    print(
        io,
        "$(length(instance.price_sensitive_loads)) price sensitive loads, ",
    )
    print(io, "$(instance.time) time steps")
    print(io, ")")
    return
end

export UnitCommitmentInstance


text:
## Default algorithm

function SciMLBase.solve(cache::LinearCache, alg::Nothing,
                         args...; kwargs...)
    @unpack A = cache
    if A isa DiffEqArrayOperator
        A = A.A
    end

    # Special case on Arrays: avoid BLAS for RecursiveFactorization.jl when
    # it makes sense according to the benchmarks, which is dependent on
    # whether MKL or OpenBLAS is being used
    if A isa Matrix
        if eltype(A) <: Union{Float32,Float64,ComplexF32,ComplexF64} &&
                    ArrayInterface.can_setindex(cache.b) && (size(A,1) <= 100 ||
                                              (isopenblas() && size(A,1) <= 500)
                                             )
            alg = RFLUFactorization()
            SciMLBase.solve(cache, alg, args...; kwargs...)
        else
            alg = LUFactorization()
            SciMLBase.solve(cache, alg, args...; kwargs...)
        end

    # These few cases ensure the choice is optimal without the
    # dynamic dispatching of factorize
    elseif A isa Tridiagonal
        alg = GenericFactorization(;fact_alg=lu!)
        SciMLBase.solve(cache, alg, args...; kwargs...)
    elseif A isa SymTridiagonal
        alg = GenericFactorization(;fact_alg=ldlt!)
        SciMLBase.solve(cache, alg, args...; kwargs...)
    elseif A isa SparseMatrixCSC
        alg = UMFPACKFactorization()
        SciMLBase.solve(cache, alg, args...; kwargs...)

    # This catches the cases where a factorization overload could exist
    # For example, BlockBandedMatrix
    elseif ArrayInterface.isstructured(A)
        alg = GenericFactorization()
        SciMLBase.solve(cache, alg, args...; kwargs...)

    # This catches the case where A is a CuMatrix
    # Which does not have LU fully defined
    elseif !(A isa AbstractDiffEqOperator)
        alg = QRFactorization()
        SciMLBase.solve(cache, alg, args...; kwargs...)

    # Not factorizable operator, default to only using A*x
    # IterativeSolvers is faster on CPU but not GPU-compatible
    elseif cache.u isa Array
        alg = IterativeSolversJL_GMRES()
        SciMLBase.solve(cache, alg, args...; kwargs...)
    else
        alg = KrylovJL_GMRES()
        SciMLBase.solve(cache, alg, args...; kwargs...)
    end
end

function init_cacheval(alg::Nothing, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)
    if A isa DiffEqArrayOperator
        A = A.A
    end

    # Special case on Arrays: avoid BLAS for RecursiveFactorization.jl when
    # it makes sense according to the benchmarks, which is dependent on
    # whether MKL or OpenBLAS is being used
    if A isa Matrix
        if eltype(A) <: Union{Float32,Float64,ComplexF32,ComplexF64} &&
                    ArrayInterface.can_setindex(b) && (size(A,1) <= 100 ||
                                              (isopenblas() && size(A,1) <= 500)
                                             )
            alg = RFLUFactorization()
            init_cacheval(alg, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)
        else
            alg = LUFactorization()
            init_cacheval(alg, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)
        end

    # These few cases ensure the choice is optimal without the
    # dynamic dispatching of factorize
    elseif A isa Tridiagonal
        alg = GenericFactorization(;fact_alg=lu!)
        init_cacheval(alg, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)
    elseif A isa SymTridiagonal
        alg = GenericFactorization(;fact_alg=ldlt!)
        init_cacheval(alg, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)
    elseif A isa SparseMatrixCSC
        alg = UMFPACKFactorization()
        init_cacheval(alg, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)

    # This catches the cases where a factorization overload could exist
    # For example, BlockBandedMatrix
    elseif ArrayInterface.isstructured(A)
        alg = GenericFactorization()
        init_cacheval(alg, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)

    # This catches the case where A is a CuMatrix
    # Which does not have LU fully defined
    elseif !(A isa AbstractDiffEqOperator)
        alg = QRFactorization()
        init_cacheval(alg, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)

    # Not factorizable operator, default to only using A*x
    # IterativeSolvers is faster on CPU but not GPU-compatible
    elseif u isa Array
        alg = IterativeSolversJL_GMRES()
        init_cacheval(alg, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)
    else
        alg = KrylovJL_GMRES()
        init_cacheval(alg, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)
    end
end


text:
include("GBMV.jl")
include("GEMV.jl")
include("GER.jl")
include("GERC.jl")
include("GERU.jl")
include("HBMV.jl")
include("HEMV.jl")
include("HER.jl")
include("HER2.jl")
include("HPMV.jl")
include("HPR.jl")
include("HPR2.jl")
include("SBMV.jl")
include("SPMV.jl")
include("SPR.jl")
include("SPR2.jl")
include("SYMV.jl")
include("SYR.jl")
include("SYR2.jl")
include("TBMV.jl")
include("TBSV.jl")
include("TPMV.jl")
include("TPSV.jl")
include("TRMV.jl")
include("TRSV.jl")


text:
testdir = dirname(@__FILE__)

run(`thrift -gen jl srvcctrl.thrift`)
run(`thrift -gen jl proto_tests.thrift`)
println("Compiled IDLs...")
cp(joinpath(testdir, "srvcctrl_impl.jl"), joinpath(testdir, "gen-jl", "srvcctrl", "srvcctrl_impl.jl"); remove_destination=true)
cp(joinpath(testdir, "proto_tests_impl.jl"), joinpath(testdir, "gen-jl", "proto_tests", "proto_tests_impl.jl"); remove_destination=true)
println("Added service implementations...")


text:
# stuff done once per exucution--not depending on r
function InitWRE!(o::StrBootTest{T}) where T
	o.LIML && o.Repl.kZ==1 && o.Nw==1 && (o.As = o.β̂s = zeros(1, o.B+1))
	o.S✻UZperpinvZperpZperp = Vector{Matrix{T}}(undef, o.Repl.kZ+1)
	o.S✻UZperp              = Vector{Matrix{T}}(undef, o.Repl.kZ+1)
	o.S✻UY                  = Matrix{Vector{T}}(undef, o.Repl.kZ+1, o.Repl.kZ+1)
	o.S✻UXinvXX             = Vector{Matrix{T}}(undef, o.Repl.kZ+1)
	o.S✻UX                  = Vector{Matrix{T}}(undef, o.Repl.kZ+1)
	o.S✻UU                  = Matrix{Matrix{T}}(undef, o.Repl.kZ+1, o.Repl.kZ+1)
	o.T1L = isone(o.Nw) ? [Matrix{T}(undef, o.Repl.kX, ncols(o.v))] :
												[Matrix{T}(undef, o.Repl.kX, length(o.WeightGrp[1])), Matrix{T}(undef, o.Repl.kX, length(o.WeightGrp[end]))]
	o.T1R = deepcopy(o.T1L)

	S✻⋂XZperp = panelcross21(o, o.Repl._X₁, o.X₂, o.DGP.Zperp, o.info✻⋂)
	S✻⋂ZperpZperp = panelcross11(o, o.Repl.Zperp, o.DGP.Zperp, o.info✻⋂)
	S✻⋂XY₂ = panelcross21(o, o.Repl._X₁, o.X₂, o.Y₂, o.info✻⋂)
	S✻⋂ZperpY₂ = panelcross11(o, o.Repl.Zperp, o.Y₂, o.info✻⋂)
	S✻⋂XX = panelcross22(o, o.Repl._X₁, o.X₂, o.DGP._X₁, o.X₂, o.info✻⋂)  # XXX X₂-X₂ part same for DGP & Repl
	S✻⋂ZperpX = panelcross12(o, o.Repl.Zperp, o.DGP._X₁, o.X₂, o.info✻⋂)
	S✻⋂Zperpy₁ = panelcross11(o, o.Repl.Zperp, o.y₁, o.info✻⋂)
	S✻⋂ZperpZR₁ = panelcross11(o, o.Repl.Zperp, o.DGP._ZR₁, o.info✻⋂)
	S✻⋂ZperpZ = panelcross11(o, o.Repl.Zperp, o.DGP._Z, o.info✻⋂)
	S✻⋂Xy₁ = panelcross21(o, o.Repl._X₁, o.X₂, o.y₁, o.info✻⋂)
	S✻⋂XZ = panelcross21(o, o.Repl._X₁, o.X₂, o.DGP._Z, o.info✻⋂)
	S✻⋂XZR₁ = panelcross21(o, o.Repl._X₁, o.X₂, o.DGP._ZR₁, o.info✻⋂)

	if o.robust && o.bootstrapt && !o.granular
		o.S✻⋂u₁XinvXX0    = o.Repl.invXX * (S✻⋂Xy₁  - S✻⋂XZperp * o.DGP.invZperpZperpZperpy₁  - o.Repl.invZperpZperpZperpX' * (S✻⋂Zperpy₁  - S✻⋂ZperpZperp * o.DGP.invZperpZperpZperpy₁))
		o.∂S✻⋂u₁XinvXX∂β̂  = o.Repl.invXX * (S✻⋂XZ   - S✻⋂XZperp * o.DGP.invZperpZperpZperpZ   - o.Repl.invZperpZperpZperpX' * (S✻⋂ZperpZ   - S✻⋂ZperpZperp * o.DGP.invZperpZperpZperpZ))
		o.∂S✻⋂u₁XinvXX∂γ̈  = o.Repl.invXX * (S✻⋂XY₂  - S✻⋂XZperp * o.DGP.invZperpZperpZperpY₂  - o.Repl.invZperpZperpZperpX' * (S✻⋂ZperpY₂  - S✻⋂ZperpZperp * o.DGP.invZperpZperpZperpY₂))
		o.∂S✻⋂u₁XinvXX∂γ̈Π̂ = o.Repl.invXX * (S✻⋂XX   - S✻⋂XZperp * o.DGP.invZperpZperpZperpX   - o.Repl.invZperpZperpZperpX' * (S✻⋂ZperpX   - S✻⋂ZperpZperp * o.DGP.invZperpZperpZperpX))
		o.∂S✻⋂u₁XinvXX∂r  = o.Repl.invXX * (S✻⋂XZR₁ - S✻⋂XZperp * o.DGP.invZperpZperpZperpZR₁ - o.Repl.invZperpZperpZperpX' * (S✻⋂ZperpZR₁ - S✻⋂ZperpZperp * o.DGP.invZperpZperpZperpZR₁))
		o.S✻⋂U₂XinvXX0 = o.∂S✻⋂u₁XinvXX∂γ̈
		o.∂S✻⋂U₂XinvXX0∂Π̂ = o.∂S✻⋂u₁XinvXX∂γ̈Π̂
	end

	info✻_✻⋂, _ = panelsetupID(o.ID✻⋂, 1:o.nbootclustvar)

	S✻XZperp = @panelsum(o, S✻⋂XZperp, info✻_✻⋂)
	S✻ZperpZperp = @panelsum(o, S✻⋂ZperpZperp, info✻_✻⋂)
	S✻XY₂ = @panelsum(o, S✻⋂XY₂, info✻_✻⋂)
	S✻ZperpY₂ = @panelsum(o, S✻⋂ZperpY₂, info✻_✻⋂)
	S✻XX = @panelsum(o, S✻⋂XX, info✻_✻⋂)
	S✻ZperpX = @panelsum(o, S✻⋂ZperpX, info✻_✻⋂)
	S✻Zperpy₁ = @panelsum(o, S✻⋂Zperpy₁, info✻_✻⋂)
	S✻ZperpZR₁ = @panelsum(o, S✻⋂ZperpZR₁, info✻_✻⋂)
	S✻ZperpZ = @panelsum(o, S✻⋂ZperpZ, info✻_✻⋂)
	S✻Xy₁ = @panelsum(o, S✻⋂Xy₁, info✻_✻⋂)
	S✻XZ = @panelsum(o, S✻⋂XZ, info✻_✻⋂)
	S✻XZR₁ = @panelsum(o, S✻⋂XZR₁, info✻_✻⋂)

	o.S✻Zperp_DGPy₁ = S✻Zperpy₁  - S✻ZperpZperp * o.DGP.invZperpZperpZperpy₁
	o.S✻Zperp_DGPZ  = S✻ZperpZ   - S✻ZperpZperp * o.DGP.invZperpZperpZperpZ
	o.S✻Zperp_DGPY₂ = S✻ZperpY₂  - S✻ZperpZperp * o.DGP.invZperpZperpZperpY₂
	o.S✻Zperp_DGPX = S✻ZperpX   - S✻ZperpZperp * o.DGP.invZperpZperpZperpX
	o.S✻Zperp_DGPZR₁ = S✻ZperpZR₁ - S✻ZperpZperp * o.DGP.invZperpZperpZperpZR₁

	o.S✻u₁X0    = S✻Xy₁  - S✻XZperp * o.DGP.invZperpZperpZperpy₁ - o.Repl.invZperpZperpZperpX' * o.S✻Zperp_DGPy₁
	o.∂S✻u₁X∂β̂  = S✻XZ   - S✻XZperp * o.DGP.invZperpZperpZperpZ  - o.Repl.invZperpZperpZperpX' * o.S✻Zperp_DGPZ
	o.∂S✻u₁X∂γ̈  = S✻XY₂  - S✻XZperp * o.DGP.invZperpZperpZperpY₂ - o.Repl.invZperpZperpZperpX' * o.S✻Zperp_DGPY₂
	o.∂S✻u₁X∂γ̈Π̂ = S✻XX   - S✻XZperp * o.DGP.invZperpZperpZperpX  - o.Repl.invZperpZperpZperpX' * o.S✻Zperp_DGPX
	o.∂S✻u₁X∂r  = S✻XZR₁ - S✻XZperp * o.DGP.invZperpZperp * o.DGP.Zperp_ZR₁ - o.Repl.invZperpZperpZperpX' * o.S✻Zperp_DGPZR₁
	o.S✻U₂X0 = o.∂S✻u₁X∂γ̈ * o.Repl.RparY
	o.∂S✻U₂X∂Π̂ = o.∂S✻u₁X∂γ̈Π̂

	if o.LIML || !o.robust || !isone(o.κ)
		# S✻y₁y₁ = panelcross11(o, o.y₁, o.y₁, o.info✻)
		# S✻y₁Y₂ = panelcross11(o, o.y₁, o.Y₂, o.info✻)
		# S✻y₁Zperp = panelcross11(o, o.y₁, o.DGP.Zperp, o.info✻)
		# S✻ZZ = panelcross11(o, o.Repl._Z, o.DGP._Z, o.info✻)
		# S✻ZZperp = panelcross11(o, o.Repl._Z, o.DGP.Zperp, o.info✻)
		# S✻y₁Z = panelcross11(o, o.y₁, o.DGP._Z, o.info✻)
		# S✻ZY₂ = panelcross11(o, o.Repl._Z, o.Y₂, o.info✻)
		# S✻ZX = panelcross12(o, o.Repl._Z, o.DGP._X₁, o.X₂, o.info✻)
		# S✻ZZR₁ = panelcross11(o, o.Repl._Z, o.DGP._ZR₁, o.info✻)
		# S✻y₁X = panelcross12(o, o.y₁, o.DGP._X₁, o.X₂, o.info✻)
		# S✻y₁ZR₁ = panelcross11(o, o.y₁, o.DGP._ZR₁, o.info✻)

		o.Reply₁_DGPy₁  = panelcross11(o, o.Repl.y₁, o.DGP.y₁, o.info✻)
		o.Reply₁_DGPZR₁ = panelcross11(o, o.Repl.y₁, o.DGP.ZR₁, o.info✻)
		o.Reply₁_DGPZ   = panelcross11(o, o.Repl.y₁, o.DGP.Z, o.info✻)
		o.Reply₁_DGPY₂  = panelcross11(o, o.Repl.y₁, o.DGP.Y₂, o.info✻)
		o.Reply₁_DGPX   = panelcross11(o, o.Repl.y₁, [o.DGP.X₁ o.DGP.X₂], o.info✻)

		o.ReplZ_DGPy₁  = panelcross11(o, o.Repl.Z, o.DGP.y₁, o.info✻)
		o.ReplZ_DGPZR₁ = panelcross11(o, o.Repl.Z, o.DGP.ZR₁, o.info✻)
		o.ReplZ_DGPZ   = panelcross11(o, o.Repl.Z, o.DGP.Z, o.info✻)
		o.ReplZ_DGPY₂  = panelcross11(o, o.Repl.Z, o.DGP.Y₂, o.info✻)
		o.ReplZ_DGPX   = panelcross11(o, o.Repl.Z, [o.DGP.X₁ o.DGP.X₂], o.info✻)

		if o.REst
			ZR₁r₁ = o.Repl.ZR₁ * o.r₁
			o.Reply₁_DGPy₁ .-= panelcross11(o, ZR₁r₁, o.DGP.y₁, o.info✻)
			o.Reply₁_DGPZR₁ .-= panelcross11(o, ZR₁r₁, o.DGP.ZR₁, o.info✻)
			o.Reply₁_DGPZ .-= panelcross11(o, ZR₁r₁, o.DGP.Z, o.info✻)
			o.Reply₁_DGPY₂ .-= panelcross11(o, ZR₁r₁, o.DGP.Y₂, o.info✻)
			o.Reply₁_DGPX .-= panelcross11(o, ZR₁r₁, [o.DGP.X₁ o.DGP.X₂], o.info✻)
		end

		# o.∂S✻y₁U₂∂γ̈ = S✻y₁Y₂ - o.Repl.invZperpZperpZperpy₁' * S✻ZperpY₂ - (S✻y₁Zperp - o.Repl.invZperpZperpZperpy₁' * S✻ZperpZperp) * o.DGP.invZperpZperpZperpY₂
		# o.S✻y₁u₁0 = -o.∂S✻y₁U₂∂γ̈ * o.Repl.RparY
		# o.∂S✻y₁u₁∂β̂ = S✻ZZ - S✻ZZperp * o.DGP.invZperpZperpZperpZ - (S✻y₁Z - S✻y₁Zperp * o.DGP.invZperpZperpZperpZ)
		# o.∂S✻y₁u₁∂γ̈ = S✻ZY₂ - o.Repl.invZperpZperpZperpy₁' * S✻ZperpY₂ - (S✻ZZperp - o.Repl.invZperpZperpZperpy₁' * S✻ZperpZperp) * o.DGP.invZperpZperpZperpY₂
		# o.∂S✻y₁u₁∂γ̈Π̂ = S✻ZX   - S✻ZZperp * o.DGP.invZperpZperpZperpX - (S✻y₁X - S✻y₁Zperp * o.DGP.invZperpZperpZperpX)
		# o.∂S✻y₁u₁∂r = S✻ZZR₁  - S✻ZZperp * o.DGP.invZperpZperpZperpZR₁ - (S✻y₁ZR₁ - S✻y₁Zperp * o.DGP.invZperpZperpZperpZR₁)
		# o.S✻y₁U₂0 = S✻y₁y₁ - o.Repl.invZperpZperpZperpy₁' * S✻Zperpy₁ - o.∂S✻y₁U₂∂γ̈
		# o.∂S✻y₁U₂∂β̂ = S✻y₁Z - S✻y₁Zperp * o.DGP.invZperpZperpZperpZ - o.Repl.invZperpZperpZperpy₁' * S✻ZperpZ
		# o.∂S✻y₁U₂∂γ̈Π̂ = S✻y₁X - o.Repl.invZperpZperpZperpy₁' * S✻ZperpX - (S✻y₁Zperp - o.Repl.invZperpZperpZperpy₁' * S✻ZperpZperp) * o.DGP.invZperpZperpZperpX
		# o.∂S✻y₁U₂∂Π̂ =  o.∂S✻y₁U₂∂γ̈Π̂
		# o.∂S✻y₁U₂∂r = S✻y₁ZR₁ - o.Repl.invZperpZperpZperpy₁' * S✻ZperpZR₁ - (S✻y₁Zperp - o.Repl.invZperpZperpZperpy₁' * S✻ZperpZperp) * o.DGP.invZperpZperpZperpZR₁

		# if o.REst
		# 	o.∂S✻y₁u₁∂β̂ .-= o.r₁' * (panelcross11(o, o.Repl._ZR₁, o.DGP._Z, o.info✻) - panelcross11(o, o.Repl._ZR₁, o.DGP.Zperp, o.info✻) * o.DGP.invZperpZperpZperpZ -
		# 											o.Repl.invZperpZperpZperpZR₁' * (S✻ZperpZ -
		# 																												S✻ZperpZperp * o.DGP.invZperpZperpZperpZ) -
		# 																												(panelcross11(o, o.Repl._ZR₁, o.DGP._Z, o.info✻) -
		# 																												panelcross11(o, o.Repl._ZR₁, o.DGP.Zperp, o.info✻) * o.DGP.invZperpZperpZperpZ -
		# 																												o.Repl.invZperpZperpZperpZR₁' * (S✻ZperpZ -
		# 																																												 S✻ZperpZperp * o.DGP.invZperpZperpZperpZ)))
		# 	o.∂S✻y₁u₁∂γ̈ .-= o.r₁' * ((panelcross11(o, o.Repl._ZR₁, o.Y₂, o.info✻) -
		# 										o.Repl.invZperpZperpZperpZR₁' * panelcross11(o, o.Repl.Zperp, o.Y₂, o.info✻)) -
		# 									(panelcross11(o, o.Repl._ZR₁, o.DGP.Zperp, o.info✻) -
		# 										o.Repl.invZperpZperpZperpZR₁' * S✻ZperpZperp) * o.DGP.invZperpZperpZperpY₂ -
		# 										((panelcross11(o, o.Repl._ZR₁, o.Y₂, o.info✻) -
		# 										 o.Repl.invZperpZperpZperpZR₁' * panelcross11(o, o.Repl.Zperp, o.Y₂, o.info✻)) -
		# 											(panelcross11(o, o.Repl._ZR₁, o.DGP.Zperp, o.info✻) -
		# 										 o.Repl.invZperpZperpZperpZR₁' * S✻ZperpZperp) * o.DGP.invZperpZperpZperpY₂))
		# 	o.∂S✻y₁u₁∂γ̈Π̂ .-= o.r₁' * ((panelcross12(o, o.Repl._ZR₁, o.DGP._X₁, o.X₂, o.info✻) -
		# 										panelcross12(o, o.Repl.Zperp * o.Repl.invZperpZperpZperpZR₁, o.DGP._X₁, o.X₂, o.info✻)) -
		# 										(panelcross11(o, o.Repl._ZR₁, o.DGP.Zperp, o.info✻) -
		# 										o.Repl.invZperpZperpZperpZR₁' * S✻ZperpZperp) * o.DGP.invZperpZperpZperpX -
		# 										((panelcross12(o, o.Repl._ZR₁, o.DGP._X₁, o.X₂, o.info✻) -
		# 									 panelcross12(o, o.Repl.Zperp * o.Repl.invZperpZperpZperpZR₁, o.DGP._X₁, o.X₂, o.info✻)) -
		# 										 (panelcross11(o, o.Repl._ZR₁, o.DGP.Zperp, o.info✻) -
		# 									 o.Repl.invZperpZperpZperpZR₁' * S✻ZperpZperp) * o.DGP.invZperpZperpZperpX))
		# 	o.∂S✻y₁u₁∂r .-= o.r₁' * ((panelcross11(o, o.Repl._ZR₁, o.DGP._ZR₁, o.info✻) -
		# 									o.Repl.invZperpZperpZperpZR₁' * S✻ZperpZR₁) -
		# 									(panelcross11(o, o.Repl._ZR₁, o.DGP.Zperp, o.info✻) -
		# 									o.Repl.invZperpZperpZperpZR₁' * S✻ZperpZperp) * o.DGP.invZperpZperpZperpZR₁ -
		# 									((panelcross11(o, o.Repl._ZR₁, o.DGP._ZR₁, o.info✻) -
		# 									o.Repl.invZperpZperpZperpZR₁' * S✻ZperpZR₁) -
		# 									 (panelcross11(o, o.Repl._ZR₁, o.DGP.Zperp, o.info✻) -
		# 									o.Repl.invZperpZperpZperpZR₁' * S✻ZperpZperp) * o.DGP.invZperpZperpZperpZR₁))
		# 	o.S✻y₁U₂0 .-= o.r₁' * (panelcross11(o, o.Repl._ZR₁, o.y₁, o.info✻) -
		# 	                        o.Repl.invZperpZperpZperpZR₁' * panelcross11(o, o.Repl.Zperp, o.y₁, o.info✻) +
		# 													(panelcross11(o, o.Repl._ZR₁, o.Y₂, o.info✻) -
		# 													o.Repl.invZperpZperpZperpZR₁' * panelcross11(o, o.Repl.Zperp, o.Y₂, o.info✻)) -
		# 													(panelcross11(o, o.Repl._ZR₁, o.DGP.Zperp, o.info✻) -
		# 													o.Repl.invZperpZperpZperpZR₁' * S✻ZperpZperp) * o.DGP.invZperpZperpZperpY₂) * o.Repl.RparY
		# 	o.∂S✻y₁U₂∂β̂ .-= o.r₁' * (panelcross11(o, o.Repl._ZR₁, o.DGP._Z, o.info✻) - panelcross11(o, o.Repl._ZR₁, o.DGP.Zperp, o.info✻) * o.DGP.invZperpZperpZperpZ -
		# 														o.Repl.invZperpZperpZperpZR₁' * S✻ZperpZ)
		# 	o.∂S✻y₁U₂∂γ̈ .-= o.r₁' * ((panelcross11(o, o.Repl._ZR₁, o.Y₂, o.info✻) -
    #                       o.Repl.invZperpZperpZperpZR₁' * panelcross11(o, o.Repl.Zperp, o.Y₂, o.info✻)) -
    #                      (panelcross11(o, o.Repl._ZR₁, o.DGP.Zperp, o.info✻) -
    #                       o.Repl.invZperpZperpZperpZR₁' * S✻ZperpZperp) * o.DGP.invZperpZperpZperpY₂)
		# 	o.∂S✻y₁U₂∂γ̈Π̂ .-= o.r₁' * ((panelcross12(o, o.Repl._ZR₁, o.DGP._X₁, o.X₂, o.info✻) -
		# 										 panelcross12(o, o.Repl.Zperp * o.Repl.invZperpZperpZperpZR₁, o.DGP._X₁, o.X₂, o.info✻)) -
    #                       (panelcross11(o, o.Repl._ZR₁, o.DGP.Zperp, o.info✻) -
    #                       o.Repl.invZperpZperpZperpZR₁' * S✻ZperpZperp) * o.DGP.invZperpZperpZperpX)
		# 	o.∂S✻y₁U₂∂Π̂ .+= o.r₁' * ((panelcross11(o, o.Repl._ZR₁, o.Y₂, o.info✻) -
		# 														o.Repl.invZperpZperpZperpZR₁' * panelcross11(o, o.Repl.Zperp, o.Y₂, o.info✻)) -
		# 														(panelcross11(o, o.Repl._ZR₁, o.DGP.Zperp, o.info✻) -
		# 														o.Repl.invZperpZperpZperpZR₁' * S✻ZperpZperp) * o.DGP.invZperpZperpZperpY₂)
		# 	o.∂S✻y₁U₂∂r .-= o.r₁' * ((panelcross11(o, o.Repl._ZR₁, o.DGP._ZR₁, o.info✻) -
		# 														o.Repl.invZperpZperpZperpZR₁' * S✻ZperpZR₁) -
		# 													(panelcross11(o, o.Repl._ZR₁, o.DGP.Zperp, o.info✻) -
		# 														o.Repl.invZperpZperpZperpZR₁' * S✻ZperpZperp) * o.DGP.invZperpZperpZperpZR₁)
		# end
	end
end

function PrepWRE!(o::StrBootTest{T}) where T
  r₁ = o.null ? [o.r₁ ; o.r] : o.r₁
	EstimateIV!(o.DGP, o, r₁)
  MakeResidualsIV!(o.DGP, o)
  Ü₂par = view(o.DGP.Ü₂ * o.Repl.RparY,:,:)

	_S✻U₂X = o.S✻U₂X0 - o.∂S✻U₂X∂Π̂ * o.DGP.Π̂ * o.Repl.RparY  # panelsum2(o, o.Repl.X₁, o.Repl.X₂, uwt, o.info✻)'
	(o.LIML || o.bootstrapt || !isone(o.κ)) &&
		(_S✻U₂Zperp = (o.S✻Zperp_DGPY₂ - o.S✻Zperp_DGPX * o.DGP.Π̂) * o.Repl.RparY)
	o.robust && o.bootstrapt && !o.granular &&
		(_S✻⋂U₂XinvXX = o.S✻⋂U₂XinvXX0 - o.∂S✻⋂U₂XinvXX0∂Π̂ * o.DGP.Π̂ * o.Repl.RparY)
	if o.LIML || !o.robust || !isone(o.κ)
		Reply₁_DGPU₂ = o.Reply₁_DGPY₂ - o.Reply₁_DGPX * o.DGP.Π̂
		ReplZ_DGPU₂  = o.ReplZ_DGPY₂  - o.ReplZ_DGPX  * o.DGP.Π̂
		Reply₁_DGPU₂par = Reply₁_DGPU₂ * o.Repl.RparY
		ReplZ_DGPU₂par = ReplZ_DGPU₂ * o.Repl.RparY
	end

  @inbounds for i ∈ 0:o.Repl.kZ  # precompute various clusterwise sums
		uwt = i>0 ? view(Ü₂par,:,i) : view(o.DGP.u⃛₁,:)::Union{Vector{T}, SubArray{T, 1}}

		# S_✻(u .* X), S_✻(u .* Zperp) for residuals u for each endog var; store transposed
		if iszero(i)  # panelsum2(o, o.Repl.X₁, o.Repl.X₂, uwt, o.info✻)'
			o.S✻UX[1] = panelsum2(o, o.Repl.X₁, o.Repl.X₂, uwt, o.info✻)'  #dropdims(o.S✻u₁X0 - o.∂S✻u₁X∂β̂ * o.DGP.β̈ + (o.∂S✻u₁X∂γ̈ - o.∂S✻u₁X∂γ̈Π̂ * o.DGP.Π̂) * o.DGP.γ̈; dims=3)
			# o.DGP.restricted &&
			# 	(o.S✻UX[1] .-= dropdims(o.∂S✻u₁X∂r * r₁; dims=3))
		else
			o.S✻UX[i+1] = panelsum2(o, o.Repl.X₁, o.Repl.X₂, uwt, o.info✻)'  # view(_S✻U₂X,:,:,i)
		end

		o.S✻UXinvXX[i+1] = o.Repl.invXX * o.S✻UX[i+1]

		if o.LIML || !isone(o.κ) || o.bootstrapt
			if iszero(i)  # panelsum(o, o.Repl.Zperp, uwt, o.infoBootData)'
				o.S✻UZperp[i+1] = view(o.S✻Zperp_DGPy₁ - o.S✻Zperp_DGPZ * o.DGP.β̈ + (o.S✻Zperp_DGPY₂ - o.S✻Zperp_DGPX * o.DGP.Π̂) * o.DGP.γ̈,:,:,1)
				o.DGP.restricted &&
					(o.S✻UZperp[i+1] .-= view(o.S✻Zperp_DGPZR₁ * r₁,:,:,1))
			else
				o.S✻UZperp[i+1] = view(_S✻U₂Zperp,:,:,i)
			end
			o.S✻UZperpinvZperpZperp[i+1] = o.Repl.invZperpZperp * o.S✻UZperp[i+1]
			o.NFE>0 && (o.CTFEU[i+1] = crosstabFE(o, uwt, o.info✻))
		end

		if o.LIML || !isone(o.κ) || !o.robust
			if iszero(i)  # panelsum2(o, o.Repl.y₁par, o.Repl.Z, uwt, o.info✻)
				o.S✻UY[1,1] = dropdims(o.Reply₁_DGPy₁ - o.Reply₁_DGPZR₁ * r₁ - o.Reply₁_DGPZ * o.DGP.β̈ - Reply₁_DGPU₂ * o.DGP.γ̈; dims=(1,3))
				for j ∈ 1:o.Repl.kZ
					o.S✻UY[1,j+1] = @view (o.ReplZ_DGPy₁ - o.ReplZ_DGPZR₁ * r₁ - o.ReplZ_DGPZ * o.DGP.β̈ - ReplZ_DGPU₂ * o.DGP.γ̈)[j,:,i]
				end
      else
				o.S✻UY[i+1,1] = @view Reply₁_DGPU₂par[1,:,i]
				for j ∈ 1:o.Repl.kZ
					o.S✻UY[i+1,j+1] = @view ReplZ_DGPU₂par[j,:,i]
				end
			end
			for j ∈ 0:i
				o.S✻UU[i+1,j+1] = vec(@panelsum(o, j>0 ? view(Ü₂par,:,j) : view(o.DGP.u⃛₁,:), uwt, o.info✻))
			end
		end

		if o.robust && o.bootstrapt
			if !o.granular  # Within each bootstrap cluster, groupwise sum by all-error-cluster-intersections of u.*X and u.Zperp (and times invXX or invZperpZperp)
				_S✻⋂U₂XinvXX = view(o.S✻⋂U₂XinvXX0 - o.∂S✻⋂U₂XinvXX0∂Π̂ * o.DGP.Π̂ * o.Repl.RparY,:,:,i)
				if iszero(i)
					_S✻⋂UXinvXX = dropdims(o.S✻⋂u₁XinvXX0 - o.∂S✻⋂u₁XinvXX∂β̂ * o.DGP.β̈ + (o.∂S✻⋂u₁XinvXX∂γ̈  - o.∂S✻⋂u₁XinvXX∂γ̈Π̂ * o.DGP.Π̂) * o.DGP.γ̈; dims=3)
					o.DGP.restricted
						(_S✻⋂UXinvXX .-= dropdims(o.∂S✻⋂u₁XinvXX∂r * r₁; dims=3))
				else
					_S✻⋂UXinvXX = view(_S✻⋂U₂XinvXX,:,:,i)
				end
				for g ∈ 1:o.N✻
					o.S✻⋂u₁XinvXX[i+1,g] = view(_S✻⋂UXinvXX,:,o._ID✻⋂[o.info✻[g]],1)'  # panelsum(o, o.Repl.XinvXX, u, o.infoCT⋂✻[g])
				end
			end

			i>0 && (o.S✻UPX[i+1] = o.Repl.XinvXX * o.S✻UX[i+1])
			o.S✻UMZperp[i+1] = o.Repl.Zperp * o.S✻UZperpinvZperpZperp[i+1]  # S_* diag(U ̈_(∥j) ) Z_⊥ (Z_⊥^' Z_⊥ )^(-1) Z_(⊥g)^'  over all g

			if iszero(i)  # subtract crosstab of observation by ∩-group of u
				o.S✻UMZperp[  1][o.crosstabBootind] .-= o.DGP.u⃛₁
			else
				o.S✻UMZperp[i+1][o.crosstabBootind] .-= view(Ü₂par,:,i)
			end

			o.NFE>0 &&
				(o.S✻UMZperp[i+1] .+= view(o.invFEwt .* o.CTFEU[i+1], o._FEID, :))  # CT_(*,FE) (U ̈_(parj) ) (S_FE S_FE^' )^(-1) S_FE
		end
  end
	nothing
end

# For WRE, and with reference to Y = [y₁ Z], given 0-based columns indexes within it, ind1, ind2, return all bootstrap realizations of
# Y[:,ind1]'((1-κ)*M_Zperp-κ*M_Xpar)*Y[:,ind2] for κ constant across replications
# ind1 can be a rowvector
# (only really the Hessian when we narrow Y to Z)
function HessianFixedkappa(o::StrBootTest{T}, ind1::Vector{S} where S<:Integer, ind2::Integer, κ::Number, w::Integer) where T
  dest = Matrix{T}(undef, length(ind1), ncols(o.v))
  @inbounds for i ∈ eachindex(ind1, axes(dest,1))
		_HessianFixedkappa!(o, view(dest,i:i,:), ind1[i], ind2, κ, w)
  end
  dest
end

function _HessianFixedkappa!(o::StrBootTest{T}, dest::AbstractMatrix, ind1::Integer, ind2::Integer, κ::Number, w::Integer) where T
  if !(o.Repl.Yendog[ind1+1] || o.Repl.Yendog[ind2+1])  # if both vars exog, result = order-0 term only, same for all draws
		!iszero(κ) &&
			coldot!(o, dest, view(o.Repl.XZ,:,ind1), view(o.Repl.invXXXZ,:,ind2))
		if !isone(κ)
			if iszero(κ)
				fill!(dest, o.Repl.YY[ind1+1,ind2+1])
			else
				dest .= κ .* dest .+ (1 - κ) .* o.Repl.YY[ind1+1,ind2+1]
			end
		end
	else
		if !iszero(κ)  # repititiveness in this section to maintain type stability
			if o.Repl.Yendog[ind1+1]
				T1L = o.T1L[isone(o.Nw) || w<o.Nw ? 1 : 2]  # use preallocated destinations
				mul!(T1L, o.S✻UX[ind1+1], o.v)
				if iszero(ind1)
					T1L .+= o.Repl.Xy₁par
				else
					T1L .+= view(o.Repl.XZ,:,ind1)
				end
				if o.Repl.Yendog[ind2+1]
					T1R = o.T1R[isone(o.Nw) || w<o.Nw ? 1 : 2]  # use preallocated destinations
					mul!(T1R, o.S✻UXinvXX[ind2+1], o.v)
					if iszero(ind2)
						T1R .+=  o.Repl.invXXXy₁par
					else
						T1R .+= view(o.Repl.invXXXZ,:,ind2)
					end
					coldot!(o, dest, T1L, T1R)
				else
					coldot!(o, dest, T1L, view(o.Repl.invXXXZ,:,ind2))
				end
			else
				if o.Repl.Yendog[ind2+1]
					T1R = o.T1R[isone(o.Nw) || w<o.Nw ? 1 : 2]  # use preallocated destinations
					mul!(T1R, o.S✻UXinvXX[ind2+1], o.v)
					if iszero(ind2)
						T1R .+=  o.Repl.invXXXy₁par
					else
						T1R .+= view(o.Repl.invXXXZ,:,ind2)
					end
					coldot!(o, dest, view(o.Repl.XZ,:,ind1), T1R)
				else
					dest .= coldot(o, view(o.Repl.XZ,:,ind1), view(o.Repl.invXXXZ,:,ind2))
				end
			end
		end
		if !isone(κ)
			if o.Repl.Yendog[ind1+1]
				T2 = o.S✻UZperpinvZperpZperp[ind1+1]'o.S✻UZperp[ind2+1]  # quadratic term
				T2[diagind(T2)] .-= ind1 ≤ ind2 ? o.S✻UU[ind2+1, ind1+1] : o.S✻UU[ind1+1, ind2+1]  # minus diagonal crosstab
				o.NFE>0 &&
					(T2 .+= o.CTFEU[ind1+1]'(o.invFEwt .* o.CTFEU[ind2+1]))
					if iszero(κ)
						dest .= o.Repl.YY[ind1+1,ind2+1] .+ colquadformminus!(o, (                            o.S✻UY[ind2+1,ind1+1] .+ o.S✻UY[ind1+1,ind2+1])'o.v, T2, o.v)
					else
						dest .=   κ .* dest .+ (1 - κ)   .* colquadformminus!(o, (o.Repl.YY[ind1+1,ind2+1] .+ o.S✻UY[ind2+1,ind1+1] .+ o.S✻UY[ind1+1,ind2+1])'o.v, T2, o.v)
					end
					# if iszero(ind1) || iszero(ind2)
				# 	tmp = iszero(ind1) ?
				# 					iszero(ind2) ? o.S✻UY[1] * 2 : o.S✻UY[ind2+1] :
				# 					o.S✻UY[ind1+1]
				# 	if iszero(κ)
				# 		dest .= o.Repl.YY[ind1+1,ind2+1] .+ colquadformminus!(o, (                            tmp)'o.v, T2, o.v)
				# 	else
				# 		dest .=  κ .* dest .+ (1 - κ)    .* colquadformminus!(o, (o.Repl.YY[ind1+1,ind2+1] .+ tmp)'o.v, T2, o.v)
				# 	end
				# else
				# 	if iszero(κ)
				# 		dest .= o.Repl.YY[ind1+1,ind2+1]
				# 	else
				# 		dest .= κ .* dest .+ (1 - κ) .* colquadformminus!(o, o.Repl.YY[ind1+1,ind2+1]'o.v, T2, o.v)
				# 	end
				# end
			elseif iszero(κ)
				dest .= o.Repl.YY[ind1+1,ind2+1]
			else
				dest .= κ .* dest .+ (1 - κ) .* o.Repl.YY[ind1+1,ind2+1]
			end
		end
  end
	nothing
end

# put threaded loops in functions to prevent type instability https://discourse.julialang.org/t/type-inference-with-threads/2004/3
function FillingLoop1!(o::StrBootTest{T}, dest::Matrix{T}, ind1::Integer, ind2::Integer, _β̂::AbstractMatrix{T}) where T
	Threads.@threads for i ∈ 1:o.clust[1].N
		PXY✻ = reshape(o.Repl.PXZ[i,ind1], :, 1)
		o.Repl.Yendog[ind1+1] && (PXY✻ = PXY✻ .+ view(o.S✻UPX[ind1+1],i,:)'o.v)

		if iszero(ind2)
			dest[i,:]   = colsum(PXY✻ .* (o.Repl.y₁[i] .- view(o.S✻UMZperp[1],i,:))'o.v)
		elseif o.Repl.Yendog[ind2+1]
			dest[i,:] .-= colsum(PXY✻ .* (o.Repl.Z[i,ind2] * _β̂ .- view(o.S✻UMZperp[ind2+1],i,:)'β̂v))
		else
			dest[i,:] .-= colsum(PXY✻ .* (o.Repl.Z[i,ind2] * _β̂))
		end
	end
	nothing
end
function FillingLoop2!(o::StrBootTest{T}, dest::Matrix{T}, ind1::Integer, ind2::Integer, _β̂::AbstractMatrix{T}) where T
	Threads.@threads for i ∈ 1:o.clust[1].N
		S = o.info⋂[i]
		PXY✻ = o.Repl.Yendog[ind1+1] ? view(o.Repl.PXZ,S,ind1) .+ view(o.S✻UPX[ind1+1],S,:) * o.v :
																	reshape(o.Repl.PXZ[S,ind1], :, 1)

		if iszero(ind2)
			dest[i,:]   = colsum(PXY✻ .* (o.Repl.y₁[S] .- view(o.S✻UMZperp[1],S,:) * o.v))
		else
			dest[i,:] .-= colsum(PXY✻ .* (o.Repl.Yendog[ind2+1] ? o.Repl.Z[S,ind2] * _β̂ .- view(o.S✻UMZperp[ind2+1],S,:) * β̂v :
																																o.Repl.Z[S,ind2] * _β̂                                       ))
		end
	end
	nothing
end
function FillingLoop3!(o::StrBootTest{T}, T₁::Matrix{T}, ind1::Integer, ind2::Integer) where T
	Threads.@threads for i ∈ 1:o.N✻
		T₁[o.IDCT⋂✻[i], i] .+= o.S✻⋂u₁XinvXX[ind2+1,i] * view(o.Repl.XZ,:,ind1)
	end
	nothing
end

# Workhorse for WRE CRVE sandwich filling
# Given a column index within it, ind1, and a matrix β̂s of all the bootstrap estimates,
# return all bootstrap realizations of P_X * Z[:,ind1]_g ' û₁g^*b
# for all groups in the intersection of all error clusterings
# return value has one row per ⋂ cluster, one col per bootstrap replication
function Filling(o::StrBootTest{T}, ind1::Integer, β̂s::AbstractMatrix) where T
	if o.granular
   	if o.Nw == 1  # create or avoid NxB matrix?
			PXY✻ = reshape(o.Repl.PXZ[:,ind1], :, 1)  # store as matrix to reduce compiler confusion
			o.Repl.Yendog[ind1+1] && (PXY✻ = PXY✻ .+ o.S✻UPX[ind1+1] * o.v)

			dest = @panelsum(o, PXY✻ .* (o.Repl.y₁ .- o.S✻UMZperp[1] * o.v), o.info⋂)

			@inbounds for ind2 ∈ 1:o.Repl.kZ
				_β̂ = view(β̂s,ind2,:)'
				dest .-= @panelsum(o, PXY✻ .* (o.Repl.Yendog[ind2+1] ?  view(o.Repl.Z,:,ind2) * _β̂ .- o.S✻UMZperp[ind2+1] * (o.v .* _β̂) :
															                              (view(o.Repl.Z,:,ind2) * _β̂)                                      ), o.info⋂)
			end
		else  # create pieces of each N x B matrix one at a time rather than whole thing at once
			dest = Matrix{T}(undef, o.clust[1].N, ncols(o.v))  # XXX preallocate this & turn Filling into Filling! ?
			@inbounds for ind2 ∈ 0:o.Repl.kZ
				ind2>0 && (β̂v = o.v .* (_β̂ = view(β̂s,ind2,:)'))

				if o.purerobust
					FillingLoop1!(o, dest, ind1, ind2, _β̂)
				else
					FillingLoop2!(o, dest, ind1, ind2, _β̂)
				end
			end
    end
  else  # coarse error clustering
		@inbounds for ind2 ∈ 0:o.Repl.kZ
			β̂v = iszero(ind2) ? o.v : o.v .* (_β̂ = -view(β̂s,ind2,:)')

			# T1 * o.v will be 1st-order terms
			T₁ = o.Repl.Yendog[ind1+1] ? o.Repl.S⋂YX[ind2+1]'o.S✻UXinvXX[ind1+1] : Matrix{T}(undef,0,0)  #  S_∩ (Y_(∥j):*X_∥) (X_∥^' X_∥)^(-1) [S_* (U ̈_(∥i):*X_∥)]^'

			if o.Repl.Yendog[ind2+1]  # add CT_(⋂,*) (P_(X_par ) Y_(pari).*U ̈_(parj) )
				if o.NClustVar == o.nbootclustvar && iszero(o.subcluster)  # simple case of one clustering: full crosstab is diagonal
					tmp = view(o.Repl.XZ,:,ind1)
					if length(T₁)>0
						T₁[diagind(T₁)] .+=         o.S✻UXinvXX[ind2+1]'tmp
					else
						T₁                = reshape(o.S✻UXinvXX[ind2+1]'tmp, :, 1)  # keep T₁ as vector rather than Diagonal matrix; probably better for fusion loop
					end
				else
					!o.Repl.Yendog[ind1+1] && (T₁ = o.JN⋂N✻)
					FillingLoop3(o, T₁, ind1, ind2)
				end
				ncols(o.Repl.Zperp) > 0 && (T₁ = T₁ .- o.Repl.S⋂PXYZperp[ind1+1] * o.S✻UZperpinvZperpZperp[ind2+1])  # subtract S_∩ (P_(X_∥ ) Y_(∥i):*Z_⊥ ) (Z_⊥^' Z_⊥ )^(-1) [S_* (U ̈_(∥j):*Z_⊥ )]^'
				o.NFE               > 0 && (T₁ = T₁ .- o.Repl.CT_FE⋂PY[ind1+1] * o.CTFEU[ind2+1])
			end

			if iszero(ind2)  # order-0 and -1 terms
				if iszero(ncols(T₁))
					dest = o.Repl.FillingT₀[ind1+1,1]
				elseif isone(ncols(T₁))
					dest = o.Repl.FillingT₀[ind1+1,1] .+ T₁ .* β̂v
				else
					dest = T₁ * o.v; dest .+= o.Repl.FillingT₀[ind1+1,1]
				end
			else  # y component
				if iszero(ncols(T₁))  # - x*β̈ components
					dest .+= o.Repl.FillingT₀[ind1+1,ind2+1] * _β̂
				elseif isone(ncols(T₁))
					dest .+= o.Repl.FillingT₀[ind1+1,ind2+1] * _β̂ .+ T₁ .* β̂v
				else
					dest .+= o.Repl.FillingT₀[ind1+1,ind2+1] * _β̂
					o.matmulplus!(dest, T₁, β̂v)
				end
			end

			if o.Repl.Yendog[ind1+1] && o.Repl.Yendog[ind2+1]
				for i ∈ 1:o.clust[1].N
					S = o.info⋂[i]
					o.colquadformminus!(dest, i, o.v, view(o.S✻UPX[ind1+1],S,:)'view(o.S✻UMZperp[ind2+1],S,:), β̂v)
				end
			end
		end
  end
  dest
end
function MakeWREStats!(o::StrBootTest{T}, w::Integer) where T
  if isone(o.Repl.kZ)  # optimized code for 1 coefficient in bootstrap regression
		if o.LIML
			YY₁₁   = HessianFixedkappa(o, [0], 0, zero(T), w)  # κ=0 => Y*MZperp*Y
			YY₁₂   = HessianFixedkappa(o, [0], 1, zero(T), w)
			YY₂₂   = HessianFixedkappa(o, [1], 1, zero(T), w)
			YPXY₁₁ = HessianFixedkappa(o, [0], 0, one(T) , w)  # κ=1 => Y*PXpar*Y
			YPXY₁₂ = HessianFixedkappa(o, [0], 1, one(T) , w)
			YPXY₂₂ = HessianFixedkappa(o, [1], 1, one(T) , w)
			YY₁₂YPXY₁₂ = YY₁₂ .* YPXY₁₂
			x₁₁ = YY₂₂ .* YPXY₁₁ .- YY₁₂YPXY₁₂      # elements of YY✻^-1 * YPXY✻ up to factor of det(YY✻)
			x₁₂ = YY₂₂ .* YPXY₁₂ .- YY₁₂ .* YPXY₂₂
			x₂₁ = YY₁₁ .* YPXY₁₂ .- YY₁₂ .* YPXY₁₁
			x₂₂ = YY₁₁ .* YPXY₂₂ .- YY₁₂YPXY₁₂
			κs = (x₁₁ .+ x₂₂)./2; κs = 1 ./ (1 .- (κs .- sqrt.(κs.^2 .- x₁₁ .* x₂₂ .+ x₁₂ .* x₂₁)) ./ (YY₁₁ .* YY₂₂ .- YY₁₂ .* YY₁₂))  # solve quadratic equation for smaller eignenvalue; last term is det(YY✻)
			!iszero(o.Fuller) && (κs .-= o.Fuller / (o._Nobs - o.kX))
			o.As = κs .* (YPXY₂₂ .- YY₂₂) .+ YY₂₂
			_β̂s = (κs .* (YPXY₁₂ .- YY₁₂) .+ YY₁₂) ./ o.As
		else
			o.As = HessianFixedkappa(o, [1], 1, o.κ, w)
			_β̂s  = HessianFixedkappa(o, [1], 0, o.κ, w) ./ o.As
		end

		if o.null
			o.numerw = _β̂s .+ (o.Repl.Rt₁ - o.r) / o.Repl.RRpar
		else
			o.numerw = _β̂s .- o.DGP.β̂₀
			isone(w) && (o.numerw[1] = _β̂s[1] + (o.Repl.Rt₁ - o.r) / o.Repl.RRpar)
		end

		@storeWtGrpResults!(o.numer, o.numerw)

		if o.bootstrapt
			if o.robust
				J⋂s = Filling(o, 1, _β̂s) ./ o.As
				@inbounds for c ∈ 1:o.NErrClustCombs  # sum sandwich over error clusterings
					nrows(o.clust[c].order)>0 &&
						(J⋂s = J⋂s[o.clust[c].order,:])
					@clustAccum!(denom, c, coldot(o, @panelsum(o, J⋂s, o.clust[c].info)))
				end
			else
				denom = (HessianFixedkappa(o, [0], 0, zero(T), w) .- 2 .* _β̂s .* HessianFixedkappa(o, [0], 1, zero(T), w) .+ _β̂s.^2 .* HessianFixedkappa(o, [1], 1, zero(T), w)) ./ o._Nobs ./ o.As  # classical error variance
			end
			@storeWtGrpResults!(o.dist, o.sqrt ? o.numerw ./ sqrt.(denom) : o.numerw .^ 2 ./ denom)
			denom *= o.Repl.RRpar[1]^2
		end

  else  # WRE bootstrap for more than 1 retained coefficient in bootstrap regression

		β̂s = zeros(T, o.Repl.kZ, ncols(o.v))
		A = Vector{Matrix{T}}(undef, ncols(o.v))

		if o.LIML
			YY✻   = [HessianFixedkappa(o, collect(0:i), i, zero(T), w) for i ∈ 0:o.Repl.kZ] # κ=0 => Y*MZperp*Y
			o.YPXY✻ = [HessianFixedkappa(o, collect(0:i), i,  one(T), w) for i ∈ 0:o.Repl.kZ] # κ=1 => Y*PXpar*Y

			@inbounds for b ∈ axes(o.v,2)
				for i ∈ 0:o.Repl.kZ
					o.YY✻_b[1:i+1,i+1]   = YY✻[i+1][:,b]  # fill uppper triangles, which is all that invsym() looks at
					o.YPXY✻_b[1:i+1,i+1] = o.YPXY✻[i+1][:,b]
				end
				o.κ = 1/(1 - eigvals(invsym(o.YY✻_b) * Symmetric(o.YPXY✻_b))[1])
				!iszero(o.Fuller) && (o.κ -= o.Fuller / (o._Nobs - o.kX))
				β̂s[:,b] = (A[b] = invsym(o.κ*o.YPXY✻_b[2:end,2:end] + (1-o.κ)*o.YY✻_b[2:end,2:end])) * (o.κ*o.YPXY✻_b[1,2:end]' + (1-o.κ)*YY✻_b[1,2:end]')
			end
		else
			δnumer =  HessianFixedkappa(o, collect(1:o.Repl.kZ), 0, o.κ, w)
			δdenom = [HessianFixedkappa(o, collect(1:i), i, o.κ, w) for i ∈ 1:o.Repl.kZ]

			@inbounds Threads.@threads for b ∈ axes(o.v,2)
				for i ∈ 1:o.Repl.kZ
					o.δdenom_b[1:i,i] = view(δdenom[i],:,b)  # fill uppper triangle
				end
				β̂s[:,b] = (A[b] = invsym(o.δdenom_b)) * view(δnumer,:,b)
			end
		end

		if o.bootstrapt
			if o.robust
				@inbounds for i ∈ 1:o.Repl.kZ  # avoid list comprehension construction because of compiler type detection issue
					o.Zyg[i] = Filling(o, i, β̂s)
				end
			else
				YY✻ = [HessianFixedkappa(o, collect(i:o.Repl.kZ), i, zero(T), w) for i ∈ 0:o.Repl.kZ]  # κ=0 => Y*MZperp*Y
			end
		end

		@inbounds for b ∈ reverse(axes(o.v,2))
			if o.null || w==1 && b==1
				o.numer_b .= o.Repl.RRpar * view(β̂s,:,b) + o.Repl.Rt₁ - o.r
			else
				o.numer_b .= o.Repl.RRpar * (view(β̂s,:,b) - o.DGP.β̂₀)
			end

			if o.bootstrapt
				if o.robust  # Compute denominator for this WRE test stat
					for i ∈ 1:o.Repl.kZ  # XXX replace with 3-D array?
						o._J⋂[:,i] = view(o.Zyg[i],:,b)
					end
					J⋂ = o._J⋂ * (A[b] * o.Repl.RRpar')

					for c ∈ 1:o.NErrClustCombs
						(!isone(o.NClustVar) && nrows(o.clust[c].order)>0) &&
							(J⋂ = J⋂[o.clust[c].order,:])
						J_b = @panelsum(o, J⋂, o.clust[c].info)
						@clustAccum!(denom, c, J_b'J_b)
					end
				else  # non-robust
					for i ∈ 0:o.Repl.kZ
						o.YY✻_b[i+1,i+1:o.Repl.kZ+1] = view(YY✻[i+1],:,b)  # fill upper triangle
					end
					denom = (o.Repl.RRpar * A[b] * o.Repl.RRpar') * [-one(T) ; β̂s[:,b]]'Symmetric(o.YY✻_b) * [-one(T) ; β̂s[:,b]] / o._Nobs  # 2nd half is sig2 of errors
				end
				if o.sqrt
					o.dist[b+first(o.WeightGrp[w])-1] = o.numer_b[1] / sqrt(denom[1])
				else
					o.dist[b+first(o.WeightGrp[w])-1] = o.numer_b'invsym(denom)*o.numer_b  # hand-code for 2-dimensional?
				end
			end
			o.numer[:,b+first(o.WeightGrp[w])-1] = o.numer_b  # slight inefficiency: in usual bootstrap-t case, only need to save numerators in numer if getdist("numer") is coming because of svmat(numer)
		end
	end
	w==1 && o.bootstrapt && (o.statDenom = denom)  # original-sample denominator
	nothing
end


text:
"""
Configurable options with defaults (values 0 and NaN indicate unlimited):
```
x_abstol::Real = 1e-50,
x_reltol::Real = 2eps(),
f_abstol::Real = 1e-50,
f_reltol::Real = 2eps(),
f_mindec::Real = 1e-4,
g_restol::Real = sqrt(eps()),
x_abstol_break::Real = Inf,
x_reltol_break::Real = Inf,
f_abstol_break::Real = Inf,
f_reltol_break::Real = Inf,
g_restol_break::Real = Inf,
f_calls_limit::Int = 0,
g_calls_limit::Int = 0,
h_calls_limit::Int = 0,
allow_f_increases::Bool = true,
min_iterations::Int = 0,
max_iterations::Int = 1_000,
warn_iterations::Int = max_iterations,
show_trace::Bool = false,
store_trace::Bool = false,
extended_trace::Bool = false,
show_every::Int = 1,
verbosity::Int = 1
```
"""
struct Options{T}
    x_abstol::T
    x_reltol::T
    f_abstol::T
    f_reltol::T
    f_mindec::T
    g_restol::T
    x_abstol_break::T
    x_reltol_break::T
    f_abstol_break::T
    f_reltol_break::T
    g_restol_break::T
    f_calls_limit::Int
    g_calls_limit::Int
    h_calls_limit::Int
    allow_f_increases::Bool
    min_iterations::Int
    max_iterations::Int
    warn_iterations::Int
    show_trace::Bool
    store_trace::Bool
    extended_trace::Bool
    show_every::Int
    verbosity::Int
end

function Options(;
        x_tol = nothing,
        f_tol = nothing,
        g_tol = nothing,
        x_abstol::Real = 1e-50,
        x_reltol::Real = 2eps(),
        f_abstol::Real = 1e-50,
        f_reltol::Real = 2eps(),
        f_mindec::Real = 1e-4,
        g_restol::Real = sqrt(eps()),
        x_abstol_break::Real = Inf,
        x_reltol_break::Real = Inf,
        f_abstol_break::Real = Inf,
        f_reltol_break::Real = Inf,
        g_restol_break::Real = Inf,
        f_calls_limit::Int = 0,
        g_calls_limit::Int = 0,
        h_calls_limit::Int = 0,
        allow_f_increases::Bool = true,
        min_iterations::Int = 0,
        max_iterations::Int = 1_000,
        warn_iterations::Int = max_iterations,
        show_trace::Bool = false,
        store_trace::Bool = false,
        extended_trace::Bool = false,
        show_every::Int = 1,
        verbosity::Int = 1)

    show_every = show_every > 0 ? show_every : 1

    if extended_trace
       show_trace = true
    end
    if !(x_tol === nothing)
        x_abstol = x_tol
    end
    if !(g_tol === nothing)
        g_restol = g_tol
    end
    if !(f_tol === nothing)
        f_reltol = f_tol
    end

    Options(promote(x_abstol, x_reltol, f_abstol, f_reltol, f_mindec, g_restol,
                    x_abstol_break, x_reltol_break, f_abstol_break, f_reltol_break, g_restol_break)...,
        f_calls_limit, g_calls_limit, h_calls_limit, allow_f_increases, min_iterations, max_iterations, warn_iterations,
        show_trace, store_trace, extended_trace, show_every, verbosity)
end

function Base.show(io::IO, o::SimpleSolvers.Options)
    for k in fieldnames(typeof(o))
        v = getfield(o, k)
        if v isa Nothing
            @printf io "%24s = %s\n" k "nothing"
        else
            @printf io "%24s = %s\n" k v
        end
    end
end

x_abstol(o::Options) = o.x_abstol
x_reltol(o::Options) = o.x_reltol
f_abstol(o::Options) = o.f_abstol
f_reltol(o::Options) = o.f_reltol
f_mindec(o::Options) = o.f_mindec
g_restol(o::Options) = o.g_restol

verbosity(o::Options) = o.verbosity


text:
include("defects.jl")
"""
    δG_R(z::ComplexF64,
     pairs::Vector{Tuple{GrapheneCoord,GrapheneCoord}},
     s::GrapheneSystem)

The correction to the real-space graphene Green's function in the presence of
defects as a function of complex energy `z`.

The function returns a vector of `ComplexF64` for each [`GrapheneCoord`](@ref)
    in `pairs`.

# Arguments
* `z`: complex energy
* `pairs`: pairs of [`GrapheneCoord`](@ref)'s for which `δG_R` is calculated
* `s`: [`GrapheneSystem`](@ref) for which `δG_R` is calculated
"""
function δG_R(
    z::ComplexF64,
    pairs::Vector{Tuple{GrapheneCoord,GrapheneCoord}},
    s::GrapheneSystem,
)
    # Calculate the scattering matrix D which is the same
    # for every pair of coordinates
    prop_mat = propagator_matrix(z, s.scattering_atoms)
    if length(s.imps) == 0
        D =
            s.Δ *
            inv(Diagonal(ones(length(s.scattering_atoms))) .- prop_mat * s.Δ)
    else
        Γ0 = 1 ./ (z .- s.imps) |> Diagonal
        D =
            (s.Δ .+ s.V * Γ0 * adjoint(s.V)) * inv(
                Diagonal(ones(length(s.scattering_atoms))) .-
                prop_mat * (s.Δ .+ s.V * Γ0 * adjoint(s.V)),
            )
    end

    PropVectorRs = [
        [graphene_propagator(x, p[2], z) for x in s.scattering_atoms] for
        p in pairs
    ]

    PropVectorLs = [
        pairs[idx][1] == pairs[idx][2] ? permutedims(PropVectorRs[idx]) :
        [
            graphene_propagator(pairs[idx][1], x, z) for
            x in s.scattering_atoms
        ] |> permutedims for idx = 1:length(pairs)
    ]

    res = [(PropVectorLs[ii]*D*PropVectorRs[ii])[1] for ii = 1:length(pairs)]
    return res
end

"""
    G_R(z::ComplexF64,
    pairs::Vector{Tuple{GrapheneCoord,GrapheneCoord}},
    s::GrapheneSystem)

The full real-space graphene Green's function in the presence of
defects as a function of complex energy `z`.

# Arguments
* `z`: complex energy
* `pairs`: pairs of [`GrapheneCoord`](@ref)'s for which `G_R` is calculated
* `s`: [`GrapheneSystem`](@ref) for which `G_R` is calculated
"""
function G_R(
    z::ComplexF64,
    pairs::Vector{Tuple{GrapheneCoord,GrapheneCoord}},
    s::GrapheneSystem,
)
    res =
        [graphene_propagator(p[1], p[2], z) for p in pairs] + δG_R(z, pairs, s)
    return res
end

"""
    δΓ(z::ComplexF64, s::GrapheneSystem)

The correction to the impurity Green's function due to the impurities'
interaction with graphene.

# Arguments
* `z`: complex energy
* `s`: [`GrapheneSystem`](@ref) for which `δΓ` is calculated
"""
function δΓ(z::ComplexF64, s::GrapheneSystem)
    if isempty(s.imps)
        error("No impurity states in the system")
    else
        Γ0 = 1 ./ (z .- s.imps) |> Diagonal
        prop_mat = propagator_matrix(z, s.scattering_atoms)
        Λ =
            prop_mat +
            prop_mat *
            s.Δ *
            inv(Diagonal(ones(length(s.scattering_atoms))) - prop_mat * s.Δ) *
            prop_mat
        res =
            Γ0 *
            adjoint(s.V) *
            Λ *
            inv(
                Diagonal(ones(length(s.scattering_atoms))) -
                s.V * Γ0 * adjoint(s.V) * Λ,
            ) *
            s.V *
            Γ0
        return res
    end
end

"""
    Γ(z::ComplexF64, s::GrapheneSystem)

The full impurity Green's function with the correction due to the impurities'
interaction with graphene.

# Arguments
* `z`: complex energy
* `s`: [`GrapheneSystem`](@ref) for which `Γ` is calculated
"""
function Γ(z::ComplexF64, s::GrapheneSystem)
    Γ0 = 1 ./ (z .- s.imps) |> Diagonal
    res = Γ0 + δΓ(z, s)
    return res
end


text:
"""
    distFunction(states1::Vector{Float64}, states2::Vector{Float64})::Array{Float64,2}

Calculates the distance abs(x-y) between all points of `states1` and `states2`.
"""
function distFunction(states1::Vector{Float64}, states2::Vector{Float64})::Array{Float64,2}
    n1= length(states1); n2= length(states2)
    dMatrix= Array{Float64}(undef, (n1,n2))
    for i= 1:n1, j= 1:n2
		dMatrix[i,j]= abs(states2[j]- states1[i])
    end
    return dMatrix
end

"""
    Wasserstein(p1::Vector{Float64}, p2::Vector{Float64}, distMatrix::Array{Float64,2}, r::Float64=1.)

Calculates the Wasserstein with marginal distributions `p1` and `p2` for the distance matrix `distMatrix.

Optionally the Hölder norm can be specified. The distance together with the optimal transport density is outputted.
"""
function Wasserstein(p1::Vector{Float64}, p2::Vector{Float64}, distMatrix::Array{Float64,2}, r::Float64=1.)
	ontoSimplex!(p1); ontoSimplex!(p2)
    n1= length(p1);   n2= length(p2)

    A= kron(ones(n2)', Matrix{Float64}(I, n1, n1))
    B= kron(Matrix{Float64}(I, n2, n2), ones(n1)')

	model = Model(Clp.Optimizer)
	set_silent(model)
	@variable(model, x[i=1:n1*n2] >= 0)
	@objective(model, Min, vec(distMatrix.^r)' * x)
	@constraint(model, [A;B] * x .== [p1;p2])
	JuMP.optimize!(model)
    return (distance= (objective_value(model))^(1/ r), π= reshape(value.(x), (n1, n2)))
end

#	Sinkhorn-Knopp iteration algorithm
function Sinkhorn(p1::Vector{Float64}, p2::Vector{Float64}, distMatrix::Array{Float64,2}, r::Float64= 1., λ::Float64= 1.)
	ontoSimplex!(p1); ontoSimplex!(p2); count= 0
	βr= Array{Float64}(undef, length(p1));
	γc= ones(size(p2))		# guess a starting value
	distMatrix.^= r; K= exp.(-λ * distMatrix)
	while count < 1000		# Sinkhorn iteration
		γc= γc./ (p2'*γc)		# rescale
		βr= p1./ (K* γc)		# vector operation
		γc= p2./ (K'* βr)		# vector operation
		count += 1				# iteration count
	end
#	println("r=", βr, p1'*βr); println("c=", γc, p2'*γc)
	π= Diagonal(βr)*K*Diagonal(γc)
	return (distance= (sum(π.* distMatrix)) ^(1/r), π= π)
end

"""
    nestedWasserstein(trr1,trr2,r=2)

Calculates the nested distance based on the Wasserstein distance for two trees `trr1` and `trr2` Optionally the Hölder norm `r` can be specified.
"""
function nestedWasserstein(trr1,trr2,r=1.0)
    T = height(trr1)
    #cumulProb!(trr1); cumulProb!(trr2);

    #d_new = Array{Float64}(undef,length(nodes(trr1,T)),length(nodes(trr2,T))).*0.0

    # initialize with all combinations of leave nodes of trr1 and trr2
    l1 = leaves(trr1)[1]; l2 = leaves(trr2)[1];
    println("distFunctionStart")
    d_new = distFunction(trr1.state[l1],trr2.state[l2])
    println("distFunctionEnd")

    for t= T-1:-1:0

        for (index1,i) in enumerate(nodes(trr1,t))
            #probability of transitioning from i to somewhere
            p1 = trr1.probability[trr1.children[i+1]]
            for (index2,j) in enumerate(nodes(trr2,t))
                p2 = trr2.probability[trr2.children[j+1]]
                #subset d_old for relevant index1/i and index2/j values
                ind1 = [findall(x->x==trr1.children[i+1][k], l1)[1] for k=1:length(trr1.children[i+1]) ]
                ind2 = [findall(x->x==trr2.children[j+1][k], l2)[1] for k=1:length(trr2.children[j+1]) ]
                d_old = d_new[ind1,ind2]
                # p1 is conditional (marginal) probability of transition for trr1
                d_new[index1,index2] = Wasserstein(p1,p2,d_old,r)[1]
            end

        end
        l1 = nodes(trr1,t); l2 = nodes(trr2,t)
        println(t)
    end
    return(d_new[1,1])
end


text:
using Documenter, TestProgress

makedocs(;
    modules=[TestProgress],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/invenia/TestProgress.jl/blob/{commit}{path}#L{line}",
    sitename="TestProgress.jl",
    authors="Invenia Technical Computing Corporation",
    assets=[
        "assets/invenia.css",
        "assets/logo.png",
    ],
)

deploydocs(;
    repo="github.com/invenia/TestProgress.jl",
)


text:
using Pkg
cd(@__DIR__)
Pkg.activate(".")
@show pwd()
include("../src/EMpht.jl")
]