text:
function get_eigen(rate, nvar, min_value = 1e-4)
    """Compute eigen values using exponential decay function.

    .. math::
        \\lambda_i = \\text{exp}^{-\\gamma(i-1)}

    :param rate: rate of exponential decay factor
    :param nvar: Number of variables (number of eigenvalues to compute)
    :param min_value: Lower limit for smallest eigenvalue
    :return: A list of eigenvalues

    """

    if rate < 0
        throw(DomainError(rate, "Need to be greater or equals to zero."))
    else
        exp.(-rate * collect(0:(nvar - 1)))
    end # if
end


text:
using ESDL
using Test
using Dates
#Open a remote cube
@testset "Cube Access" begin

c=Cube()

@testset "Access single variable" begin
d = subsetcube(c,variable="air_temperature_2m",lon=(10,11),lat=(51,50),
                time=(Date("2002-01-01"),Date("2008-12-31")))

@test d.subset==(18:21, 17:20, 93:414)
@test d.axes[1].values==10.125:0.25:10.875
@test d.axes[2].values==50.875:-0.25:50.125
@test isa(d.axes[1],LonAxis)
@test isa(d.axes[2],LatAxis)
end

@testset "Access multiple variables" begin
d2 = subsetcube(c,variable=["air_temperature_2m","gross_primary_productivity"],lon=(10,11),lat=(50,51),
                time=(Date("2002-01-01"),Date("2008-12-31")))

@test d2.cataxis.values == ["air_temperature_2m", "gross_primary_productivity"]
foreach(d2.cubelist) do cc
  @test cc.subset ==(18:21, 17:20, 93:414)
  @test cc.axes[1].values==10.125:0.25:10.875
  @test cc.axes[2].values==50.875:-0.25:50.125
  @test first(cc.axes[3].values) == Dates.Date(2002,1,5)
  @test last(cc.axes[3].values) == Dates.Date(2008, 12, 30)
end
end

@testset "Test values in MemCube" begin
d = subsetcube(c,variable="air_temperature_2m",lon=(10,11),lat=(51,50),
                time=(Date("2002-01-01"),Date("2008-12-31")))
d2 = subsetcube(c,variable=["air_temperature_2m","gross_primary_productivity"],lon=(10,11),lat=(50,51),
                time=(Date("2002-01-01"),Date("2008-12-31")))
data1=readcubedata(d)
data2=readcubedata(d2)

@test size(data1.data)==(4,4,322)
@test size(data2.data)==(4,4,322,2)

@test isapprox(data1.data[1,1,1:10],Float32[267.9917, 269.9631, 276.71036, 280.88998,
  280.90665, 277.02243, 274.5466, 276.919, 279.96243, 279.42276])

@test isapprox(data2.data[1,1,1:10,1],Float32[267.9917, 269.9631, 276.71036, 280.88998,
  280.90665, 277.02243, 274.5466, 276.919, 279.96243, 279.42276])

@test caxes(data1)[1:2]==CubeAxis[LonAxis(10.125:0.25:10.875),LatAxis(50.875:-0.25:50.125)]

tax = caxes(data1)[3]
@test isa(tax, TimeAxis)
@test tax.values[1] == Date(2002,01,05)
@test tax.values[end] == Date(2008,12,30)
@test length(tax.values) == 7*46

@test caxes(data2)[[1,2,4]]==CubeAxis[LonAxis(10.125:0.25:10.875),LatAxis(50.875:-0.25:50.125),VariableAxis(["air_temperature_2m","gross_primary_productivity"])]

tax = caxes(data2)[3]
@test isa(tax, TimeAxis)
@test tax.values[1] == Date(2002,01,05)
@test tax.values[end] == Date(2008,12,30)
@test length(tax.values) == 7*46

end

@testset "Coordinate extraction" begin
d = subsetcube(c,variable="air_temperature_2m",lon=(10,11),lat=(51,50),
                time=(Date("2002-01-01"),Date("2008-12-31")))
data1=readcubedata(d)
# Test reading of coordinate list
ll=[10.1 50.2;10.5 51.1;10.8 51.1]
llcube = readcubedata(extractLonLats(data1,ll))
@test llcube.data[1,:]==data1.data[1,4,:]
@test llcube.data[2,:]==data1.data[3,1,:]
@test llcube.data[3,:]==data1.data[4,1,:]
end

@testset "Accessing regions" begin
#Test access datacube by region
d3 = subsetcube(c,variable="gross_primary_productivity",region="Austria",time=Date("2005-01-01"))
@test d3.axes==[LonAxis(9.625:0.25:14.875),LatAxis(48.875:-0.25:47.375)]
end

@testset "Saving and loading cubes" begin
  d = subsetcube(c,variable="air_temperature_2m",lon=(10,31),lat=(51,50),
                  time=(Date("2002-01-01"),Date("2008-12-31")))
  data1=readcubedata(d)
  #Test saving cubes
  dire=tempname()
  ESDLdir(dire)
  saveCube(data1,"mySavedCube")

  data3=readcubedata(loadCube("mySavedCube"))
  @test data1.axes==data3.axes
  @test data1.data==data3.data

  # Test loadOrGenerate macro
  d=subsetcube(c,time=Date(2001)..Date(2005),lon=(10,11),lat=(50,51),variable=["gross_primary_productivity","net_ecosystem_exchange"])

  rmCube("Anomalies")
  @loadOrGenerate danom=>"Anomalies" begin
      danom = removeMSC(d)
  end

  @test danom isa ESDL.CubeMem

  @loadOrGenerate danom=>"Anomalies" begin
      error("This should never execute")
  end;
  @test danom isa ESDL.Cubes.ConcatCube

  saveCube(danom, "mySavedZArrayCube")

  @test danom isa ESDL.Cubes.ConcatCube

  danom=readcubedata(danom)
  danom2=readcubedata(loadCube("mySavedZArrayCube"))

  @test danom.axes==danom2.axes
  @test all(map(isequal,danom.data,danom2.data))

  ncf = tempname()
  exportcube(danom,ncf)

  using NetCDF, Dates
  #Test exportcube
  @test ncread(ncf,"Lon") == 10.125:0.25:10.875
  @test ncread(ncf,"Lat") == 50.875:-0.25:50.125
  @test ncgetatt(ncf,"Time","units") == "days since 2001-01-01"
  @test getAxis("Time",danom).values .- DateTime(2001) == Millisecond.(Day.(ncread(ncf,"Time")))

  anc = replace(ncread(ncf,"gross_primary_productivity")[:,:,:],-9999.0=>missing)
  @test all(isequal.(anc, permutedims(danom[:,:,:,1],(2,3,1))))
  neear = replace(danom[:,:,:,2],missing=>-9999.0)
  @test all(isequal.(ncread(ncf,"net_ecosystem_exchange")[:,:,:],permutedims(neear,(2,3,1))))
end
end

ext: jl


text:
function is_palidrome(n::Int64)
	a=Int64[]
	while n>=1
		push!(a,n%10);
		n=n-n%10;
		n/=10
		n=Int64(n)

	end
	len=length(a);
	if( len==1)  return true
	end
	for i=1:len
		if a[i]!=a[len+1-i]
			return false
		end
	end
	return true
end

function compute()
	N=1000
	ans=0
	for i=1:N
		for j=1:N
			if is_palidrome(i*j)
				ans=max(ans,i*j)
			end
		end
	end
	return ans
end

print(compute())


text:
using QMTK
using QMTK.Consts.Pauli
using Compat.Test

@testset "local hamiltonian" begin

    mat = σ₁⊗σ₂
    h = LocalHamiltonian(mat)
    rhs = SubSites(Bit, 1, 0)
    rhs_idx = Int(rhs) + 1
    itr = LHIterator(h, rhs)

    for (val, lhs) in itr
        lhs_idx = Int(lhs) + 1
        @test val == mat[lhs_idx, rhs_idx]
    end

end # testset


text:
# julia functions for bootstrap
# Claudia October 2015
# Cecile April 2016

"""
    readBootstrapTrees(listfile; relative2listfile=true)

Read the list of file names in `listfile`, then read all the trees in each of
these files. Output: vector of vectors of trees (networks with h>0 allowed).

`listfile` should be the name of a file containing the path/name to multiple
bootstrap files, one on each line (no header). Each named bootstrap file should
contain multiple trees, one per line (such as bootstrap trees from a single gene).

The path/name to each bootstrap file should be relative to `listfile`.
Otherwise, use option `relative2listfile=false`, in which case the file names
are interpreted as usual: relative to the user's current directory
if not given as absolute paths.
"""
function readBootstrapTrees(filelist::AbstractString; relative2listfile=true::Bool)
    filelistdir = dirname(filelist)
    bootfiles = CSV.read(filelist, header=false, types=Dict(1=>String))
    size(bootfiles)[2] > 0 ||
        error("there should be a column in file $filelist: with a single bootstrap file name on each row (no header)")
    ngenes = size(bootfiles)[1]
    bf = (relative2listfile ? joinpath.(filelistdir, bootfiles[!,1]) : bootfiles[!,1])
    treelists = Array{Vector{HybridNetwork}}(undef, ngenes)
    for igene in 1:ngenes
        treelists[igene] = readMultiTopology(bf[igene])
        print("read $igene/$ngenes bootstrap tree files\r") # using \r for better progress display
    end
    return treelists
end

"""
    sampleBootstrapTrees(vector of tree lists; seed=0::Integer, generesampling=false, row=0)
    sampleBootstrapTrees!(tree list, vector of tree lists; seed=0::Integer, generesampling=false, row=0)

Sample bootstrap gene trees, 1 tree per gene.
Set the seed with keyword argument `seed`, which is 0 by default.
When `seed=0`, the actual seed is set using the clock.
Assumes a vector of vectors of networks (see `readBootstrapTrees`),
each one of length 1 or more (error if one vector is empty, tested in `bootsnaq`).

- site resampling: always, from sampling one bootstrap tree from each given list.
  This tree is sampled at **random** unless `row>0` (see below).
- gene resampling: if `generesampling=true` (default is false),
  genes (i.e. lists) are sampled with replacement.
- `row=i`: samples the ith bootstrap tree for each gene.
  `row` is turned back to 0 if gene resampling is true.

output: one vector of trees. the modifying function (!) modifies the input tree list and returns it.
"""
function sampleBootstrapTrees(trees::Vector{Vector{HybridNetwork}};
                              seed=0::Integer, generesampling=false::Bool, row=0::Integer)
    bootTrees = Array{HybridNetwork}(undef, length(trees))
    sampleBootstrapTrees!(bootTrees, trees, seed=seed, generesampling=generesampling, row=row)
end

function sampleBootstrapTrees!(bootTrees::Vector{HybridNetwork}, trees::Vector{Vector{HybridNetwork}};
                              seed=0::Integer, generesampling=false::Bool, row=0::Integer)
    numgen = length(trees) ## number of genes
    numgen>0 || error("needs at least 1 array of trees")
    numgen <= length(bootTrees) || error("the input tree list needs to be of length $numgen at least")
    if (generesampling) row=0; end
    if row==0
      if seed == 0
        t = time()/1e9
        a = split(string(t),".")
        seed = parse(Int,a[2][end-4:end]) #better seed based on clock
        println("using seed $(seed) for bootstrap trees")
      end
      Random.seed!(seed)
      if generesampling
        indxg = sample(1:numgen, numgen) # default is with replacement. good!
      end
    end
    for g in 1:numgen
        ig = (generesampling ? indxg[g] : g )
        if row==0
          indxt = sample(1:length(trees[ig]),1)[1]
        else
          indxt = row
          length(trees[ig]) >= row || error("gene $g has fewer than $row bootstrap trees.")
        end
        bootTrees[g] = trees[ig][indxt]
    end
    return bootTrees
end

"""
    sampleCFfromCI(data frame, seed=0)
    sampleCFfromCI!(data frame, seed=0)

Read a data frame containing CFs and their credibility intervals, and
sample new obsCF uniformly within the CIs.
These CFs are then rescaled to sum up to 1 for each 4-taxon sets.
Return a data frame with taxon labels in first 4 columns, sampled obs CFs in columns 5-7
and credibility intervals in columns 8-13.

- The non-modifying function creates a new data frame (with re-ordered columns) and returns it.
  If `seed=-1`, the new df is a deep copy of the input df, with no call to the random
  number generator. Otherwise, `seed` is passed to the modifying function.
- The modifying function overwrites the input data frame with the sampled CFs and returns it.
  If `seed=0`, the random generator is seeded from the clock. Otherwise the random generator
  is seeded using `seed`.

Warning: the modifying version does *not* check the data frame: assumes correct columns.

optional argument: `delim=','` by default: how columns are delimited.
"""
function sampleCFfromCI(df::DataFrame, seed=0::Integer)
    @debug "order of columns should be: t1,t2,t3,t4,cf1234,cf1324,cf1423,cf1234LO,cf1234HI,..."
    size(df,2) == 13 || size(df,2) == 14 || @warn "sampleCFfromCI function assumes table from TICR: CF, CFlo, CFhi"
    obsCFcol = [findfirst(isequal(:CF12_34), DataFrames.names(df)),
                findfirst(isequal(:CF13_24), DataFrames.names(df)),
                findfirst(isequal(:CF14_23), DataFrames.names(df))]
    nothing ∉ obsCFcol || error("""CF columns were not found: should be named like 'CF12_34'""")
    obsCFcol == [5,8,11] ||
        @warn """CF columns were found, but not in the expected columns.
                Lower/upper bounds of credibility intervals assumed in columns 6,7, 9,10 and 12,13."""
    colsTa = [1,2,3,4]        # column numbers for taxon names
    colsCI = [6,7,9,10,12,13] # for lower/upper CI bounds
    length(findall(in(obsCFcol), colsTa)) ==0 ||
        error("CFs found in columns 1-4 where taxon labels are expected")
    length(findall(in(obsCFcol), colsCI)) ==0 ||
        error("CFs found in columns where credibility intervals are expected")
    newdf = df[:, [colsTa; obsCFcol; colsCI] ]
    if seed==-1
      return newdf
    else
      return sampleCFfromCI!(newdf::DataFrame, seed)
    end
end

function sampleCFfromCI!(df::DataFrame, seed=0::Integer)
    if seed == 0
        t = time()/1e9
        a = split(string(t),".")
        seed = parse(Int,a[2][end-4:end]) #better seed based on clock
        println("using seed $(seed) for bootstrap table")
    end
    Random.seed!(seed)
    for i in 1:size(df,1)
        c1 = (df[i, 9]-df[i, 8])*rand()+df[i, 8]
        c2 = (df[i,11]-df[i,10])*rand()+df[i,10]
        c3 = (df[i,13]-df[i,12])*rand()+df[i,12]
        suma = c1+c2+c3
        df[i,5] = c1/suma
        df[i,6] = c2/suma
        df[i,7] = c3/suma
    end
    return df
end

sampleCFfromCI(file::AbstractString; delim=','::Char,seed=0::Integer) = sampleCFfromCI(CSV.read(file, delim=delim),seed)

# function that will do bootstrap of snaq estimation in series
# it repeats optTopRuns nrep times
# it has the same arguments as optTopRuns except for:
# - need data table of CF with conf intervals (instead of d DataCF),
#   or vector of vector of HybridNetworks
# - new argument nrep: number of bootstrap replicates (default 10)
# - new argument runs2: percentage of bootstrap replicates to start in the best network, by default 0.25
# - new argument bestNet: to start the optimization. if prcnet>0.0 and bestNet is not input as argument from a previous run, it will estimate it inside
# - quartetfile if it was used in original data ("none" if all quartets used)
# recall: optTopRuns! does *not* modify its input starting network
function optTopRunsBoot(currT0::HybridNetwork, data::Union{DataFrame,Vector{Vector{HybridNetwork}}},
                        hmax::Integer, liktolAbs::Float64, Nfail::Integer, ftolRel::Float64,ftolAbs::Float64,xtolRel::Float64,xtolAbs::Float64,
                        verbose::Bool, closeN::Bool, Nmov0::Vector{Int},
                        runs1::Integer, outgroup::AbstractString, filename::AbstractString, seed::Integer, probST::Float64,
                        nrep::Integer, runs2::Integer, bestNet::HybridNetwork, quartetfile::AbstractString)
    println("BOOTSTRAP OF SNAQ ESTIMATION")
    writelog = true
    if filename != ""
        logfile = open(string(filename,".log"),"w")
        write(logfile, "BOOTSTRAP OF SNAQ ESTIMATION \n")
    else
        writelog = false
        logfile = stdout
    end

    inputastrees = isa(data, Vector{Vector{HybridNetwork}})
    inputastrees || isa(data, DataFrame) ||
        error("Input data not recognized: $(typeof(data))")

    if runs1>0 && runs2>0
        str = """Will use this network as starting topology for $runs1 run(s) for each bootstrap replicate:
                 $(writeTopologyLevel1(currT0))
                 and this other network for $runs2 run(s):
                 $(writeTopologyLevel1(bestNet))
                 """
        writelog && write(logfile, str)
        print(str)
    end

    if inputastrees # allocate memory, to be re-used later
        newtrees = sampleBootstrapTrees(data, row=1)
        newd = readTrees2CF(newtrees, quartetfile=quartetfile, writeTab=false, writeSummary=false)
        taxa = unionTaxa(newtrees)
    else
        newdf = sampleCFfromCI(data, -1) # column names check, newdf has obsCF in columns 5-7
        # seed=-1: deep copy only, no rand()
        newd = readTableCF!(newdf, collect(1:7)) # allocate memory for DataCF object
    end
    if runs1>0 && isTree(currT0) # get rough first estimate of branch lengths in startnet
        updateBL!(currT0, newd)
    end

    if seed == 0
        t = time()/1e9
        a = split(string(t),".")
        seed = parse(Int,a[2][end-4:end]) #better seed based on clock
    end
    println("main seed $(seed)")
    writelog && write(logfile,"\nmain seed $(seed)\n")
    writelog && flush(logfile)
    Random.seed!(seed)
    seedsData = round.(Int,floor.(rand(nrep)*100000)) # seeds to sample bootstrap data
    if runs1>0
        seeds = round.(Int,floor.(rand(nrep)*100000)) # seeds for all optimizations from currT0
    end
    if runs2>0
      seedsOtherNet = round.(Int,floor.(rand(nrep)*100000)) # for runs starting from other net
    end

    bootNet = HybridNetwork[]

    writelog && write(logfile,"\nBEGIN: $(nrep) replicates\n$(Libc.strftime(time()))\n")
    writelog && flush(logfile)
    for i in 1:nrep
        str = "\nbegin replicate $(i)\nbootstrap data simulation: seed $(seedsData[i])\n"
        writelog && write(logfile, str)
        print(str)
        if !inputastrees
            sampleCFfromCI!(newdf, seedsData[i])
            readTableCF!(newd, newdf, [5,6,7])
        else
            sampleBootstrapTrees!(newtrees, data, seed=seedsData[i])
            calculateObsCFAll!(newd,taxa) # do not use readTrees2CF: to save memory and gc time
        end
        if runs1>0
            str = "estimation, $runs1 run" * (runs1>1 ? "s" : "") * ": seed $(seeds[i])\n"
            writelog && write(logfile, str)
            print(str)
            rootname = ""
            @debug begin rootname = string(filename,"_",i);
                         "rootname set to $rootname"; end
            net1 = optTopRuns!(currT0, liktolAbs, Nfail, newd, hmax,ftolRel, ftolAbs, xtolRel, xtolAbs, verbose, closeN, Nmov0, runs1, outgroup,
                               rootname,seeds[i],probST)
            if runs2==0
                net = net1
            end
        end
        if runs2>0
            str = "estimation, $runs2 run" * (runs2>1 ? "s" : "") * " starting from other net: seed $(seedsOtherNet[i])\n"
            writelog && write(logfile, str)
            print(str)
            rootname = ""
            @debug begin rootname = string(filename,"_",i,"_startNet2");
                         "rootname set to $rootname"; end
            net2 = optTopRuns!(bestNet, liktolAbs, Nfail, newd, hmax,ftolRel, ftolAbs, xtolRel, xtolAbs, verbose, closeN, Nmov0, runs2, outgroup,
                               rootname,seedsOtherNet[i],probST)
            if runs1==0
                net = net2
            end
        end
        if runs1>0 && runs2>0
            net = (net1.loglik < net2.loglik ? net1 : net2)
        end
        writelog && flush(logfile)

        push!(bootNet, net)
        str = (outgroup=="none" ? writeTopologyLevel1(net) : writeTopologyLevel1(net,outgroup))
        if writelog
            write(logfile, str)
            write(logfile,"\n")
            flush(logfile)
        end
        println(str) # net also printed by each optTopRuns! but separately from the 2 starting points
    end # of the nrep bootstrap replicates
    writelog && close(logfile)

    if writelog
      s = open(string(filename,".out"),"w")
      for n in bootNet
        if outgroup == "none"
            write(s,"$(writeTopologyLevel1(n))\n")
        else
            write(s,"$(writeTopologyLevel1(n,outgroup))\n")
        end
        # "with -loglik $(n.loglik)" not printed: not comparable across bootstrap networks
      end
      close(s)
    end
    return bootNet
end

# like snaq, only calls optTopRunsBoot
# undocumented arguments: closeN, Nmov0
"""
    bootsnaq(T::HybridNetwork, df::DataFrame)
    bootsnaq(T::HybridNetwork, vector of tree lists)

Bootstrap analysis for SNaQ.
Bootstrap data can be quartet concordance factors (CF),
drawn from sampling uniformly in their credibility intervals,
as given in the data frame `df`.
Alternatively, bootstrap data can be gene trees sampled from
a vector of tree lists: one list of bootstrap trees per locus
(see `readBootstrapTrees` to generate this,
from a file containing a list of bootstrap files: one per locus).

From each bootstrap replicate, a network
is estimated with snaq!, with a search starting from topology `T`.
Optional arguments include the following, with default values in parentheses:

- `hmax` (1): max number of reticulations in the estimated networks
- `nrep` (10): number of bootstrap replicates.
- `runs` (10): number of independent optimization runs for each replicate
- `filename` ("bootsnaq"): root name for output files. No output files if "".
- `seed` (0 to get a random seed from the clock): seed for random number generator
- `otherNet` (empty): another starting topology so that each replicate will start prcnet% runs on otherNet and (1-prcnet)% runs on `T`
- `prcnet` (0): percentage of runs starting on `otherNet`; error if different than 0.0, and otherNet not specified.
- `ftolRel`, `ftolAbs`, `xtolRel`, `xtolAbs`, `liktolAbs`, `Nfail`,
  `probST`, `verbose`, `outgroup`: see `snaq!`, same defaults.

If `T` is a tree, its branch lengths are first optimized roughly with [`updateBL!`](@ref)
(by using the average CF of all quartets defining each branch and calculating the coalescent units
corresponding to this quartet CF).
If `T` has one or more reticulations, its branch lengths are taken as is to start the search.
The branch lengths of `otherNet` are always taken as is to start the search.
"""
function bootsnaq(startnet::HybridNetwork, data::Union{DataFrame,Vector{Vector{HybridNetwork}}};
                  hmax=1::Integer, liktolAbs=likAbs::Float64, Nfail=numFails::Integer,
                  ftolRel=fRel::Float64, ftolAbs=fAbs::Float64, xtolRel=xRel::Float64, xtolAbs=xAbs::Float64,
                  verbose=false::Bool, closeN=true::Bool, Nmov0=numMoves::Vector{Int},
                  runs=10::Integer, outgroup="none"::AbstractString, filename="bootsnaq"::AbstractString,
                  seed=0::Integer, probST=0.3::Float64, nrep=10::Integer, prcnet=0.0::Float64,
                  otherNet=HybridNetwork()::HybridNetwork, quartetfile="none"::AbstractString)

    inputastrees = isa(data, Vector{Vector{HybridNetwork}})
    inputastrees || isa(data, DataFrame) ||
        error("Input data not recognized: $(typeof(data))")

    if !inputastrees
    (DataFrames.names(data)[[6,7,9,10,12,13]] == [:CF12_34_lo,:CF12_34_hi,:CF13_24_lo,:CF13_24_hi,:CF14_23_lo,:CF14_23_hi]) ||
      @warn """assume table with CI from TICR: CFlo, CFhi in columns 6,7; 9,10; and 12,13.
              Found different column names: $(DataFrames.names(data)[[6,7,9,10,12,13]])"""
    else # check 1+ genes, each with 1+ trees, all with h=0.
        ngenes = length(data)
        ngenes > 0 || error("empty list of bootstrap trees (0 genes)")
        for igene in 1:ngenes
            btr = data[igene]
            length(btr) > 0 || error("no bootstrap trees for $(igene)th gene")
            for itree in 1:length(btr)
                btr[itree].numHybrids == 0 || error("network $itree is not a tree for $(igene)th gene")
            end
        end
    end
    prcnet >= 0 || error("percentage of times to use the best network as starting topology should be positive: $(prcnet)")
    prcnet = (prcnet <= 1.0) ? prcnet : prcnet/100
    runs2 = round(Int, runs*prcnet)            # runs starting from otherNet
    runs1 = runs - runs2                       # runs starting from startnet

    if runs1>0
        startnet=readTopologyLevel1(writeTopologyLevel1(startnet)) # does not modify startnet outside
        flag = checkNet(startnet,true) # light checking only
        flag && error("starting topology suspected not level-1")
        try
            checkNet(startnet)
        catch err
            println("starting topology is not of level 1:")
            rethrow(err)
        end
    end
    runs2 == 0 || otherNet.numTaxa > 0 ||
        error("""otherNet not given and prcnet>0. Please set prcnet to 0 to start optimizations
                from the same network always, or else provide an other network "otherNet"
                to start the optimization from this other network in pcrnet % of runs.""")
    if runs2 > 0
        otherNet=readTopologyLevel1(writeTopologyLevel1(otherNet))
        flag = checkNet(otherNet,true) # light checking only
        flag && error("starting topology 'otherNet' suspected not level-1")
        try
            checkNet(otherNet)
        catch err
            println("starting topology 'otherNet' not a level 1 network:")
            rethrow(err)
        end
    end

    # for multiple alleles: expand into two leaves quartets like sp1 sp1 sp2 sp3.
    if (@isdefined originald) && !isempty(originald.repSpecies) ## not defined if treefile empty, but not needed
        expandLeaves!(originald.repSpecies,startnet)
    end

    optTopRunsBoot(startnet,data,hmax, liktolAbs, Nfail,ftolRel, ftolAbs, xtolRel, xtolAbs,
                   verbose, closeN, Nmov0, runs1, outgroup, filename,
                   seed, probST, nrep, runs2, otherNet, quartetfile)
end

"""
`treeEdgesBootstrap(boot_net::Vector{HybridNetwork}, ref_net::HybridNetwork)`

Read a list of bootstrap networks (`boot_net`) and a reference network (`ref_net`),
and calculate the bootstrap support of the tree edges in the reference network.
All minor hybrid edges (γ<0.5) are removed to extract the major tree from
each network. All remaining edges are tree edges, each associated with a bipartition.

output:

- a data frame with one row per tree edge and two columns: edge number, bootstrap support
  (as a percentage)
- the major tree from the reference network, where minor hybrid edges (with γ<0.5)
  have been removed.
"""
function treeEdgesBootstrap(net::Vector{HybridNetwork}, net0::HybridNetwork)
    # estimated network, major tree and matrix
    S = tipLabels(net0)
    tree0 =majorTree(net0, unroot=true)
    M0 = tree2Matrix(tree0,S, rooted=false)

    M = Matrix[]
    tree = HybridNetwork[]
    for n in net
        t = majorTree(n, unroot=true)
        push!(tree,t)
        mm = tree2Matrix(t,S, rooted=false)
        push!(M,mm)
    end

    df = DataFrame(edgeNumber=Int[], proportion=Float64[])

    for i in 1:size(M0,1) #rows in M0: internal edges
        cnt = 0 #count
        for j in 1:length(M) #for every M
            for k in 1:size(M[j],1) #check every row in M
                if M0[i,2:end] == M[j][k,2:end] || M0[i,2:end] == map(x->(x+1)%2,M[j][k,2:end]) #same row
                    #println("found same row: $(M0[i,2:end]) and $(M[j][k,2:end])")
                    cnt += 1
                    break
                end
            end
        end
        push!(df,[M0[i,1] cnt*100/length(M)])
        # fixit later, when edges have an attribute for bootstrap support:
        # set BS to cnt/length(M) for the edge numbered M0[i,1].
    end
    @info """edge numbers in the data frame correspond to the current edge numbers in the network.
       If the network is modified, the edge numbers in the (modified) network might not correspond
       to those in the bootstrap table. Plot the bootstrap values onto the current network with
       plot(network_name, edgeLabel=bootstrap_table_name)"""
    return df, tree0
end

"""
    hybridDetection(net::Vector{HybridNetwork}, net1::HybridNetwork, outgroup::AbstractString)

function can only compare hybrid nodes in networks that have the same underlying major tree
also, need to root all networks in the same place, and the root has to be compatible with the
direction of the hybrid edges

it computes the rooted hardwired distance between networks, the root matters.
input: vector of bootstrap networks (net), estimated network (net1), outgroup

returns

- a matrix with one row per bootstrap network, and 2*number of hybrids in net1,
  column i corresponds to whether hybrid i (`net1.hybrid[i]`) is found in the bootstrap network,
  column 2i+1 corresponds to the estimated gamma on the bootstrap network
  (0.0 if hybrid not found).
  To know the order of hybrids, print `net1.hybrid` or `h.name for h in net1.hybrid`
- list of discrepant trees (trees not matching the main tree in net1)
"""
function hybridDetection(net::Vector{HybridNetwork}, net1::HybridNetwork, outgroup::AbstractString)
    tree1 = majorTree(net1, unroot=true)
    rootnet1 = deepcopy(net1)
    rootatnode!(rootnet1,outgroup)

    # HF for "hybrid found?"
    HFmat = zeros(length(net),net1.numHybrids*2)

    # discrepant trees
    discTrees = HybridNetwork[]
    # major trees
    majorTrees = HybridNetwork[]

    i = 1
    for n in net
        tree = majorTree(n, unroot=true)
        push!(majorTrees,tree)
        RFmajor = hardwiredClusterDistance(tree, tree1, false)
        if RFmajor != 0
            push!(discTrees,tree)
            i+=1
            continue # skip replicate if major tree is *not* correct
        end

        found = zeros(Bool,   net1.numHybrids) # repeats false
        gamma = zeros(Float64,net1.numHybrids) # repeats 0.0
        # re-root estimated network if not rooted correctly
        reroot = true
        if length(n.node[n.root].edge) == 2 # check if root connects to correct outgroup
            for e in n.node[n.root].edge
                for node in e.node
                    if node.name == outgroup
                        reroot = false
                        break
                    end
                end
                if (!reroot) break end
            end
        end
        !reroot || println("Will need to reroot the estimated network...")
        for trueh = 1:net1.numHybrids
            netT = deepcopy(rootnet1)
            displayedNetworkAt!(netT, netT.hybrid[trueh]) # bug: need correct attributes to re-root later...
            for esth = 1:n.numHybrids
                netE = deepcopy(n)
                displayedNetworkAt!(netE, netE.hybrid[esth])
                if reroot
                    rootatnode!(netE, outgroup) # if re-rooting is not possible,
                end                         # then the hybridization doesn't match.
                if hardwiredClusterDistance(netT, netE, true) == 0 # true: rooted
                    found[trueh] = true
                    node = netE.hybrid[1]
                    edges = hybridEdges(node)
                    edges[2]. hybrid || error("edge should be hybrid")
                    !edges[2]. isMajor || error("edge should be minor hybrid")
                    edges[2].gamma <= 0.5 || error("gamma should be less than 0.5")
                    gamma[trueh] = edges[2].gamma
                    break # to exit loop over esth
                end
            end
        end
        HFmat[i,1:net1.numHybrids] = found
        HFmat[i,(net1.numHybrids+1):end] = gamma
        i+=1
    end
    treeMatch=size(HFmat[sum(HFmat,2).>0,:],1) #number of bootstrap trees that match tree1
    println("$(treeMatch) out of $(length(net)) bootstrap major trees match with the major tree in the estimated network")
    println("order of hybrids:")
    for h in net1.hybrid
        println("$(h.name)")
    end
    return HFmat,discTrees
end

"""
`hybridBootstrapSupport(boot_net::Vector{HybridNetwork}, ref_net::HybridNetwork; rooted=false)`

Match hybrid nodes in a reference network with those in an array of networks,
like bootstrap networks.
All networks must be fully resolved, and on the same taxon set.
If `rooted=true`, all networks are assumed to have been properly rooted beforehand.
Otherwise, the origin of each hybrid edge is considered as an unrooted bipartition (default).

Two hybrid edges in two networks are said to match if they share the same "hybrid" clade
(or recipient) and the same "donor clade", which is a sister to the hybrid clade in the network.
Since a hybrid clade has 2 parent edges, it is sister to two clades simultaneously: one is its
major sister (following the major hybrid edge with γ>0.5) and one is its minor sister (following
the major hybrid edge with γ<0.5).

To calculate these hybrid and sister clades at a given hybrid node, all other hybrid edges are
first removed from the network. Then, the hybrid clade is the hardwired cluster (descendants) of
either hybrid edge and major/minor clade is the hardwired cluster of the sibling edge of the
major/minor hybrid parent.
If `rooted=false`, sister clades are considered as bipartitions.

Output:

1. a "node" data frame (see below)
2. an "edge" data frame (see below)
3. a "clade" data frame to describe the make up of all clades found as hybrids or sisters,
   starting with a column `taxa` that lists all taxa. All other columns correspond to a given
   clade and contain true/false values. `true` means that a given taxon belongs in a given clade.
   For a clade named `H1`, for instance, and if the data frame was named `cla`, the
   list of taxa in this clade can be obtained with `cla[:taxa][cla[:H1]]`.
4. an array of gamma values, with one row for each bootstrap network and two columns (major/minor)
   for each hybrid edge in the reference network. If this hybrid edge was found in the bootstrap network
   (i.e. same hybrid and sister clades, after removal of all other hybrid nodes),
   its bootstrap gamma value is recorded here. Otherwise, the gamma entry is 0.0.
5. a vector with the number of each hybrid edge in the reference network,
   in the same order as for the columns in the array of gamma values above.

The "node" data frame has one row per clade and 9 columns giving:

   - `:clade`: the clade's name, like the taxon name (if a hybrid is a single taxon) or
     the hybrid tag (like 'H1') in the reference network
   - `:node`: the node number in the reference network. missing if the clade is not in this network.
   - `:hybridnode`: typically the same node number as above, except for hybrid clades in the
     reference network. For those, the hybrid node number is listed here.
   - `:edge`: number of the parent edge, parent to the node in column 2,
     if found in the ref network. missing otherwise.
   - `:BS_hybrid`: percentage of bootstrap networks in which the clade is found to be a hybrid clade.
   - `:BS_sister`: percentage of bootstrap networks in which the clade is found to be sister to
     some hybrid clade (sum of the next 2 columns)
   - `:BS_major_sister`: percentage of bootstrap networks in which the clade is found to be the
     major sister to some hybrid clade
   - `:BS_minor_sister`: same as previous, but minor
   - `:BS_hybrid_samesisters`: percentage of bootstrap networks in which the clade is found to be
     a hybrid and with the same set of sister clades as in the reference network.
     Applies to hybrid clades found in the reference network only, missing for all other clades.

The "edge" data frame has one row for each pair of clades, and 8 columns:

  - `:edge`: hybrid edge number, if the edge appears in the reference network. missing otherwise.
  - `:hybrid_clade`: name of the clade found to be a hybrid, descendent of 'edge'
  - `:hybrid`: node number of that clade, if it appears in the reference network. missing otherwise.
  - `:sister_clade`: name of the clade that is sister to 'edge', i.e. be sister to a hybrid
  - `:sister`: node number of that clade, if in the ref network.
  - `:BS_hybrid_edge`: percentage of bootstrap networks in which 'edge' is found to be a hybrid
     edge, i.e. when the clade in the 'hybrid' column is found to be a hybrid and the clade in
     the 'sister' column is one of its sisters.
  - `:BS_major`: percentage of bootstrap networks in which 'edge' is found to be a major hybrid
     edge, i.e. when 'hybrid' is found to be a hybrid clade and 'sister' is found to be its
     major sister.
  - `:BS_minor`: same as previous, but minor
"""
function hybridBootstrapSupport(nets::Vector{HybridNetwork}, refnet::HybridNetwork;
         rooted=false::Bool)
    numNets = length(nets)
    numNets>0 || error("there aren't any test (bootstrap) networks")
    numHybs = refnet.numHybrids
    numHybs>0 || error("there aren't any hybrid in reference network")
    try directEdges!(refnet)
    catch err
        if isa(err, RootMismatch)
            err.msg *= "\nPlease change the root in reference network (see rootatnode! or rootonedge!)"
        end
        rethrow(err)
    end
    taxa = tipLabels(refnet)
    ntax = length(taxa)

    # extract hardwired clusters of each tree edge in major tree of reference net,
    # and of each hybrid edge after all other hybrids are removed following the major edge.
    clade = Vector{Bool}[] # list all clades in reference and in bootstrap networks
    treenode = Int[]     # node number for each clade: of tree node if found in reference net
    treeedge = Int[]     # number of tree edge corresponding to the clade (parent of treenode)
    leafname = AbstractString[] # "" if internal node, leaf name if leaf
    # clade, treenode, treeedge: same size. leafname and hybparent: same size initially only
    hybind = Int[]       # indices in 'clade', that appear as hybrid clades in reference
    hybnode = Int[]      # for those hybrid clades: number of hybrid node
    majsisedge = Int[]   # for those hybrid clades: number of major sister edge (in ref net)
    minsisedge = Int[]   #                                    minor
    majsisind = Int[]    # for those hybrid clades: index of major sister clade in 'clade'
    minsisind = Int[]    #                                   minor
    # hybind, hybnode, majsis*, minsis*: same size = number of hybrid nodes in reference network

    if !rooted && length(refnet.node[refnet.root].edge)==2 && any(e -> e.hybrid, refnet.node[refnet.root].edge)
        refnet = deepcopy(refnet) # new binding inside function
        fuseedgesat!(refnet.root, refnet)
    end # issues otherwise: correct tree edge for root bipartitions, find sister clades, ...

    reftre = majorTree(refnet, unroot=true)
    skipone = (!rooted && length(reftre.node[reftre.root].edge)<3) # not count same bipartition twice
    for pe in reftre.edge
        hwc = hardwiredCluster(pe,taxa) # not very efficient, but human readable
        if skipone && refnet.node[refnet.root] ≡ getParent(pe) && sum(hwc)>1
            skipone = false             # wrong algo for trivial 2-taxon rooted tree (A,B);
            println("skip edge $(pe.number)")
        else
            push!(clade, hwc)
            push!(treeedge, pe.number)
            cn = getChild(pe) # child node of pe
            push!(treenode, cn.number)
            push!(leafname, (cn.leaf ? cn.name : ""))
        end
    end
    hybparent = zeros(Int,length(clade)) # 0 if has no hybrid parent node. Index in hybnode otherwise.
    for trueh = 1:numHybs
        net0 = deepcopy(refnet)
        displayedNetworkAt!(net0, net0.hybrid[trueh]) # removes all minor hybrid edges but one
        hn = net0.hybrid[1]
        hemaj, hemin, ce = hybridEdges(hn) # assumes no polytomy at hybrid nodes and correct node.hybrid
        (hemin.hybrid && !hemin.isMajor) || error("edge should be hybrid and minor")
        (hemaj.hybrid &&  hemaj.isMajor) || error("edge should be hybrid and major")
        ic = findfirst(isequal(ce.number), treeedge)
        ic !== nothing || error("hybrid node $(hn.number): child edge not found in major tree")
        hybparent[ic] = trueh
        push!(hybind,ic)
        push!(hybnode, hn.number)
        push!(majsisedge, hemaj.number)
        push!(minsisedge, hemin.number)
        for sis in ["min","maj"]
          he = (sis=="min" ? hemin : hemaj)
          pn = getParent(he) # parent node of sister (origin of gene flow if minor)
          atroot = (!rooted && pn ≡ net0.node[net0.root]) # polytomy at root pn of degree 3: will exclude one child edge
          hwc = zeros(Bool,ntax) # new binding each time. pushed to clade below.
          for ce in pn.edge    # important if polytomy
            if ce ≢ he && pn ≡ getParent(ce)
                hw = hardwiredCluster(ce,taxa)
                if atroot && any(hw & clade[ic]) # sister clade intersects child clade
                    (hw & clade[ic]) == clade[ic] ||
                        @warn "weird clusters at the root in reference, hybrid node $(hn.number)"
                else
                    hwc .|= hw
                end
            end
          end
          i = findfirst(isequal(hwc), clade)
          if (!rooted && i===nothing) i = findfirst(isequal(.!hwc), clade) end
          i !== nothing || error(string("hyb node $(hn.number): ",sis,"or clade not found in main tree"))
          if (sis=="min") push!(minsisind, i)
          else            push!(majsisind, i)
          end
          if sis=="min"  # need to get clade not in main tree: hybrid + minor sister
            pe = nothing # looking for the (or one) parent edge of pn
            for ce in pn.edge
              if pn ≡ getChild(ce)
                  pe=ce
                  break
              end
            end
            # pe == nothing if minor hybrid is at root and rooted=true. Will just miss edge and node number
            # for that clade, but in that case, the minor sister might have been assigned that edge anyway...
            if pe != nothing
              hwc = hardwiredCluster(pe,taxa)
              i = findfirst(isequal(hwc), clade) # i>0: (hybrid + minor sister) can be in main tree if
              # hybrid origin is ancestral, i.e. hybrid clade is nested within minor sister.
              if i===nothing
                push!(clade, hwc)
                push!(treenode, pn.number)
                push!(treeedge, pe.number)
                push!(hybparent, 0)
                push!(leafname, (pn.leaf ? pn.name : ""))
              end
            end
          end
        end
    end
    # for cl in clade @show taxa[cl]; end; @show treenode; @show treeedge; @show hybparent; @show leafname
    # @show hybind; @show hybnode; @show majsisedge; @show minsisedge; @show majsisind; @show minsisind

    # node-associated summaries:
    nclades = length(clade)
    nh = length(hybnode)
    nedges = 2*nh
    BShyb        = zeros(Float64, nclades) # clade = hybrid
    BSmajsis     = zeros(Float64, nclades) # clade = major sister of some hybrid
    BSminsis     = zeros(Float64, nclades) # clade = minor sister of some hybrid
    # edge-associated summaries, i.e. associated to a pair of clades (hyb,sis):
    hybcladei   = repeat(hybind, inner=[2]) # indices in 'clade'
    siscladei   = Array{Int}(undef, nedges) # edge order: (major then minor) for all hybrids
    edgenum     = Array{Int}(undef, nedges)
    for i=1:nh
        siscladei[2*i-1] = majsisind[i]
        siscladei[2*i]   = minsisind[i]
        edgenum[2*i-1] = majsisedge[i]
        edgenum[2*i]   = minsisedge[i]
    end # less desirable order: [majsisind; minsisind], corresponds to outer=[2] for hybcladei
    BShybmajsis = zeros(Float64, nedges) # one clade = hybrid, one clade = major sister
    BShybminsis = zeros(Float64, nedges) # one clade = hybrid, one clade = minor sister
    # network*edge-associated values: only keep those for edges in ref net
    gamma = zeros(Float64,numNets,nedges)
    # node-associate 3-way partition summary:
    BShybsamesis = zeros(Float64, nh) # same hybrid & sister set as in ref net

    nextnum = max(maximum([n.number for n in refnet.edge]),
                  maximum([n.number for n in refnet.node]) ) + 1

    for i = 1:numNets
        net = nets[i]
        length(tipLabels(net))==ntax || error("networks have non-matching taxon sets")
        try directEdges!(net) # make sure the root is admissible
        catch err
          if isa(err, RootMismatch)
            err.msg *= "\nPlease change the root in test network (see rootatnode! or rootatedge!)"
          end
          rethrow(err)
        end
        for esth = 1:net.numHybrids   # try to match estimated hybrid edge
            hwcPar = zeros(Bool,ntax)   # minor sister clade
            hwcChi = zeros(Bool,ntax)   #       child  clade
            hwcSib = zeros(Bool,ntax)   # major sister clade
            net1 = deepcopy(net)
            displayedNetworkAt!(net1, net1.hybrid[esth])
            hn = net1.hybrid[1]
            if !rooted && length(net1.node[net1.root].edge)==2 && any(e -> e.hybrid, net1.node[net1.root].edge)
                fuseedgesat!(net1.root, net1)
            end
            hemaj, hemin, ce = hybridEdges(hn) # assumes no polytomy at hybrid node
            (hemin.hybrid && !hemin.isMajor) || error("edge should be hybrid and minor")
            (hemaj.hybrid &&  hemaj.isMajor) || error("edge should be hybrid and major")
            hardwiredCluster!(hwcChi,hemin,taxa)
            for sis in ["min","maj"]
                he = (sis=="min" ? hemin : hemaj)
                pn = getParent(he) # parent of hybrid edge
                atroot = (!rooted && pn ≡ net1.node[net1.root])
                # if at root: exclude the child edge in the same cycle as he.
                # its cluster includes hwcChi. all other child edges do not interest hwcChi.
                # if (atroot) @show i; @warn "$(sis)or edge is at the root!"; end
                for ce in pn.edge
                  if ce ≢ he && pn ≡ getParent(ce)
                    hwc = hardwiredCluster(ce,taxa)
                    if !atroot || sum(hwc .& hwcChi) == 0 # empty intersection
                      if (sis=="maj") hwcSib .|= hwc;
                      else            hwcPar .|= hwc; end
                    elseif (hwc .& hwcChi) != hwcChi
                        @warn "weird clusters at the root. bootstrap net i=$i, hybrid $(net.hybrid[esth].name)"
                    end
                  end
                end # will use complement too: test network may be rooted differently
            end
            # @show taxa[hwcChi]; @show taxa[hwcPar]
            if all(hwcPar) || all(hwcSib) || all(.!hwcPar) || all(.!hwcSib)
                @warn "parent or sibling cluster is full or empty. bootstrap net i=$i, hybrid $(net.hybrid[esth].name)"
            end

            ihyb = findfirst(isequal(hwcChi), clade)
            newhyb = (ihyb===nothing) # hwcChi not found in clade list
            if newhyb
              push!(clade, hwcChi)
              push!(treenode, nextnum)
              push!(treeedge, nextnum)
              push!(BShyb,    1.0)
              push!(BSmajsis, 0.0)
              push!(BSminsis, 0.0)
              ihyb = length(clade) # was nothing; converted to nex available integer
              nextnum += 1
            else
              BShyb[ihyb] += 1.0
            end

            iSmaj = findfirst(isequal(hwcSib), clade)
            iSmin = findfirst(isequal(hwcPar), clade)
            if (!rooted && iSmaj===nothing) iSmaj = findfirst(isequal(.!hwcSib), clade) end
            if (!rooted && iSmin===nothing) iSmin = findfirst(isequal(.!hwcPar), clade) end
            newmaj = (iSmaj===nothing)
            if newmaj
              push!(clade, hwcSib)
              push!(treenode, nextnum)
              push!(treeedge, nextnum)
              push!(BShyb,    0.0)
              push!(BSmajsis, 1.0)
              push!(BSminsis, 0.0)
              iSmaj = length(clade) # was nothing; now integer
              nextnum += 1
            else
              BSmajsis[iSmaj] += 1.0
            end
            newmin = (iSmin===nothing)
            if newmin
              push!(clade, hwcPar)
              push!(treenode, nextnum)
              push!(treeedge, nextnum)
              push!(BShyb,    0.0)
              push!(BSmajsis, 0.0)
              push!(BSminsis, 1.0)
              iSmin = length(clade) # was nothing; now integer
              nextnum += 1
            else
              BSminsis[iSmin] += 1.0
            end
            samehyb = (newhyb ? Int[] : findall(isequal(ihyb), hybcladei))
            for sis in ["maj","min"]
                newpair = newhyb || (sis=="min" ? newmin : newmaj)
                iSsis = (sis=="min" ? iSmin : iSmaj)
                if !newpair # hyb and sis clades were already found, but not sure if together
                    iish = findfirst(isequal(iSsis), siscladei[samehyb])
                    if iish !== nothing # pair was indeed found
                        ipair = samehyb[iish]
                        if (sis=="min") BShybminsis[ipair] += 1.0
                        else            BShybmajsis[ipair] += 1.0 end
                        if ipair <= nedges # pair = edge in reference network
                            gamma[i,ipair] = (sis=="min" ? hemin.gamma : hemaj.gamma)
                        end
                    else newpair = true; end
                end
                if newpair
                    push!(hybcladei, ihyb)
                    push!(siscladei, iSsis)
                    push!(BShybminsis, (sis=="min" ? 1.0 : 0.0))
                    push!(BShybmajsis, (sis=="min" ? 0.0 : 1.0))
                end
            end
            if ihyb<=nclades && hybparent[ihyb]>0 # hyb clade match in ref net
                th = hybparent[ihyb]
                if ((iSmaj==majsisind[th] && iSmin==minsisind[th]) ||
                    (iSmin==majsisind[th] && iSmaj==minsisind[th]) )
                    BShybsamesis[th] += 1
                end
            end
        end # loop over hybrid nodes
    end # loop over bootstrap nets

    fac = 100/numNets
    BShyb    *= fac # size: length(clade)
    BSmajsis *= fac
    BSminsis *= fac
    BShybmajsis  *= fac # size: length of hybcladei or siscladei
    BShybminsis  *= fac
    BShybsamesis *= fac # size: length(hybnode) = nedges/2

    # combine results into 2 dataframes: node & edge summaries
    # and 1 dataframe to describe clades
    # first: detect nodes with all BS=0
    keepc = ones(Bool,length(clade)) # keep clade h in output tables?
    for h=nclades:-1:1
        if BShyb[h]==0.0 && BSmajsis[h]==0.0 && BSminsis[h]==0.0
            keepc[h] = false
            deleteat!(BShyb,h); deleteat!(BSmajsis,h); deleteat!(BSminsis,h)
        end
    end
    nkeepc = sum(keepc)
    # clade descriptions
    resCluster = DataFrame(taxa=taxa)
    cladestr = Array{String}(undef, length(clade))
    rowh = 1
    for h=1:length(clade)
        nn = treenode[h] # node number in ref net
        cladestr[h] = string("c_", (nn<0 ? "minus" : ""), abs(nn))
        # column name for a clade at node 5: "c_5". At node -5: "c_minus5"
        # because symbol :c_-5 causes an error.
        if h <= nclades &&  leafname[h] != "" # replace "c5" by leaf name, e.g. "taxon1"
            cladestr[h] = leafname[h]
        end
        if h <= nclades &&  hybparent[h]>0 # replace "c5" or leaf name by hybrid name, e.g. "H1"
            na = refnet.hybrid[hybparent[h]].name
            cladestr[h] = (na=="" ? string("H", refnet.hybrid[hybparent[h]].number) : replace(na, r"^#" => ""))
        end
        if keepc[h]
            rowh += 1
            insertcols!(resCluster, rowh, Symbol(cladestr[h]) => clade[h])
        end
    end
    # node summaries
    resNode = DataFrame(clade=cladestr[keepc],
                        node=allowmissing(treenode[keepc]),
                        hybridnode=allowmissing(treenode[keepc]),
                        edge=allowmissing(treeedge[keepc]),
                        BS_hybrid=BShyb, BS_sister = BSmajsis + BSminsis,
                        BS_major_sister=BSmajsis, BS_minor_sister=BSminsis,
                        BS_hybrid_samesisters=Vector{Union{Missing,Float64}}(undef, nkeepc))
    rowh = 1
    for h=1:length(clade)
        if h <= nclades && keepc[h] && hybparent[h]>0
            resNode[rowh,:hybridnode]            =      hybnode[hybparent[h]]
            resNode[rowh,:BS_hybrid_samesisters] = BShybsamesis[hybparent[h]]
        elseif keepc[h]
            resNode[rowh,:BS_hybrid_samesisters] = missing
        end
        if h>nclades # clade *not* in the reference network
            resNode[rowh,:node] = missing
            resNode[rowh,:hybridnode] = missing
            resNode[rowh,:edge] = missing
        end
        if keepc[h]  rowh += 1; end
    end
    insertcols!(resNode, 10, :BS_all => resNode[!,:BS_hybrid]+resNode[!,:BS_sister])
    sort!(resNode, [:BS_all,:BS_hybrid]; rev=true)
    select!(resNode, Not(:BS_all))
    # edge summaries
    resEdge = DataFrame(edge = Vector{Union{Int, Missing}}(undef, length(hybcladei)),
                        hybrid_clade=cladestr[hybcladei],
                        hybrid=Vector{Union{Int, Missing}}(treenode[hybcladei]),
                        sister_clade=cladestr[siscladei],
                        sister=Vector{Union{Int, Missing}}(treenode[siscladei]),
                        BS_hybrid_edge = BShybmajsis+BShybminsis,
                        BS_major=BShybmajsis, BS_minor=BShybminsis)
    for i=1:length(hybcladei)
        h = hybcladei[i]
        if h <= nclades && hybparent[h]>0
            resEdge[i,:hybrid] = hybnode[hybparent[h]]
        end
        if h>nclades            resEdge[i,:hybrid]=missing; end
        if siscladei[i]>nclades resEdge[i,:sister]=missing; end
        if i <= nedges
             resEdge[i,:edge] = edgenum[i]
        else resEdge[i,:edge] = missing
        end
    end
    o = [1:nedges; sortperm(resEdge[nedges+1:length(hybcladei),:BS_hybrid_edge],rev=true) .+ nedges]
    return resNode, resEdge[o,:], resCluster, gamma, edgenum
end

"""
    summarizeHFdf(HFmat::Matrix)

Summarize data frame output from [`hybridDetection`](@ref).
Output: dataframe with one row per hybrid, and 5 columns:

- hybrid index (order from estimated network, see [`hybridDetection`](@ref),
- number of bootstrap trees that match the
  underlying tree of estimated network
- number of bootstrap networks that have the hybrid
- mean estimated gamma in the bootstrap networks that have the hybrid
- sd estimated gamma in the bootstrap networks that have the hybrid also

last row has index -1, and the third column has the number of networks
that have all hybrids (hybrid index, mean gamma, sd gamma are
meaningless in this last row)
"""
function summarizeHFdf(HFmat::Matrix)
    HFmat2 = HFmat[sum(HFmat,2) .>0,:]
    gt = size(HFmat2,1)
    total = size(HFmat,1)
    numH = round(Int,size(HFmat,2)/2)
    df = DataFrame(hybrid=Int[],goodTrees=Float64[],netWithHybrid=Float64[],meanGamma=Float64[], sdGamma=Float64[])
    for i in 1:numH
        mat = HFmat2[HFmat2[:,i] .> 0, :]
        n = size(mat,1)
        g = mean(mat[:,round(Int,numH+i)])
        s = std(mat[:,round(Int,numH+i)])
        push!(df,[i gt n g s])
    end
    which = Bool[]
    for i in 1:size(HFmat2,1)
        push!(which,sum(HFmat2[i,1:numH]) == numH)
    end
    push!(df, [-1 gt sum(which) -1.0 -1.0])
    return df
end


text:
################################################################################
#
#     Clgrp.jl : Class group computation of maximal orders in number fields
#
# This file is part of hecke.
#
# Copyright (c) 2015: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
# (C) 2015, 2016 Claus Fieker
#
################################################################################
#
# Todo:
#  - make sure the precision for LLL is high enough (by checking that the
#    resulting elements have a reasonable norm/ length? theory?)
#    done
#  - add reasonable verbose printing
#    done
#  - write hnf from upper_triangular
#  - understand/use profiling information (memory, ...)
#
#  - need different norm function: modular resultant? with a large known
#    factor AND without any proof...
#    otherwise, it takes much too long if the ideal is non-trivial
#  DONE (norm_div)
#
#  - move the various factor, issmooth and similar to a sensible
#    spot. This has nothing to do with class groups
#  - the SingleP version:
#      figure out if a union is better than the current version
#      ie have a version for easy primes
#         and a dumb version as fallback
#      make sure stuff works with fractionals as well!
#      just scaling by the den can affect smoothness (incomplete factor base)
#
#
# Clean up:
#  - write show functions
#  - export
#
# Note: enumerating x,0,0,0 is pointless unless x==1
#
################################################################################

export class_group, FactorBase, issmooth, factor, lll_basis,
       unit_group_fac_elem, unit_group, regulator

add_verbose_scope(:ClassGroup)
add_verbose_scope(:ClassGroup_time)
add_verbose_scope(:ClassGroup_gc)

add_assert_scope(:ClassGroup)
set_assert_level(:ClassGroup, 0)
set_assert_level(:LatEnum, 0)

include("Clgp/Ctx.jl")
include("Clgp/FacBase_Euc.jl")
include("Clgp/FacBase_Idl.jl")
include("Clgp/Main_enum.jl")
include("Clgp/Map.jl")
include("Clgp/Proof.jl")
include("Clgp/Rel_add.jl")
include("Clgp/Rel_enum.jl")
include("Clgp/Rel_LLL.jl")
include("Clgp/Main_LLL.jl")
include("Clgp/Rel_Schmettow.jl")
include("Clgp/Saturate.jl")
using .RelSaturate

################################################################################
#
#  Main function
#
################################################################################

function class_group_ctx(O::NfOrd; bound::Int = -1, method::Int = 3, large::Int = 1000, redo::Bool = false, use_aut::Bool = false)

  if !redo
    try
      c = _get_ClassGrpCtx_of_order(O)::ClassGrpCtx{SMat{fmpz}}
      return c
    catch e
    end
  end

  if bound == -1
    bound = Int(ceil(log(abs(discriminant(O)))^2*0.3))
    (bound == 0) && (bound = 1)
  end

  c = class_group_init(O, bound, complete = false, use_aut = use_aut)::ClassGrpCtx{SMat{fmpz}}

  _set_ClassGrpCtx_of_order(O, c)

  c.B2 = bound * large

  if method == 2
    class_group_find_relations2(c)
  else
    d = root(abs(discriminant(O)), 2)
    c.expect = class_group_expected(c, 100)
    class_group_via_lll(c)
  end

  return c
end

################################################################################
#
#  Verification
#
################################################################################

# think of a sensible function name

function _validate_class_unit_group(c::ClassGrpCtx, U::UnitGrpCtx)
  @vprint :UnitGroup 1 "Validating unit group and class group ... \n"
  O = U.order

  # The residue of the zeta function cannot be computed for degree 1 (K = Q),
  # so we shortcircuit it.

  h = class_group_current_h(c)
  if degree(O) == 1
    if h == 1 && U.tentative_regulator == 1
      return fmpz(1)
    else
      error("Something odd for K = Q")
    end
  end

  if !isdefined(U, :torsion_units)
    @vprint :UnitGroup 1 "Computing torsion structure ... \n"
    #U.torsion_units = torsion_units(O)
    g, ord = torsion_units_gen_order(O)
    U.torsion_units_order = ord
    U.torsion_units_gen = g
  end

  @vprint :UnitGroup 1 "Torsion structure done!\n"

  w = U.torsion_units_order

  r1, r2 = signature(O)

  if !isdefined(U, :residue)
    @vprint :UnitGroup 1 "Computing residue of Dedekind zeta function ... \n"
    U.residue = zeta_log_residue(O, 0.6931/2)  #log(2)/2
  end
  residue = U.residue

  pre = prec(parent(residue))

  Ar = ArbField(pre, false)

  loghRtrue = Ar(residue) + log(Ar(w)*sqrt(abs(Ar(discriminant(O))))//(Ar(2)^(r1+r2) * const_pi(Ar)^r2))

  # I should assert that radius(loghRtrue) < log(2)

  @assert isfinite(loghRtrue)

  while true
    loghRapprox = log(h* abs(U.tentative_regulator))

    @assert isfinite(loghRapprox)

    if contains(loghRtrue, loghRapprox)
      return fmpz(1)
    elseif !overlaps(loghRtrue, loghRapprox)
      e = exp(loghRapprox - loghRtrue)
      e_fmpz = abs_upper_bound(e, fmpz)
      @vprint :ClassGroup 1 "validate called, index bound is $e_fmpz\n"
      return e_fmpz
    end

    error("Not yet implemented")
  end
end

function class_group_current_h(c::ClassGrpCtx)
  module_trafo_assure(c.M)
  c.h = check_index(c.M)
  return c.h
end

function _class_unit_group(O::NfOrd; bound::Int = -1, method::Int = 3, large::Int = 1000, redo::Bool = false, unit_method::Int = 1, use_aut::Bool = false, GRH::Bool = true)

  @vprint :UnitGroup 1 "Computing tentative class and unit group ... \n"

  @v_do :UnitGroup 1 pushindent()
  c = class_group_ctx(O, bound = bound, method = method, large = large, redo = redo, use_aut = use_aut)
  @v_do :UnitGroup 1 popindent()

  if c.finished
    U = _get_UnitGrpCtx_of_order(O)
    @assert U.finished
    @vprint :UnitGroup 1 "... done (retrieved).\n"
    return c, U, fmpz(1)
  end

  @vprint :UnitGroup 1 "Tentative class number is now $(c.h)\n"

  U = UnitGrpCtx{FacElem{nf_elem, AnticNumberField}}(O)

  need_more = true

  hnftime = 0.0

  r = 0

  do_units = true
  while true
    @v_do :UnitGroup 1 pushindent()
    if do_units
      if unit_method == 1
        @vtime_add_elapsed :UnitGroup 1 c :unit_time r = _unit_group_find_units(U, c)
      else
        @vtime_add_elapsed :UnitGroup 1 c :unit_hnf_time module_trafo_assure(c.M)
        @vtime_add_elapsed :UnitGroup 1 c :unit_time r = _unit_group_find_units_with_transform(U, c)
      end
      @v_do :UnitGroup 1 popindent()
    end
    if isone(r)  # use saturation!!!!
      idx = _validate_class_unit_group(c, U)
      @assert idx == _validate_class_unit_group(c, U)
      stable = 3.5
      if iszero(c.sat_done)
        @vprint :ClassGroup 1 "Finite index, saturating at 2\n"
        while saturate!(c, U, 2, stable)
          @vtime_add_elapsed :UnitGroup 1 c :unit_time r = _unit_group_find_units(U, c)
        end
        idx = _validate_class_unit_group(c, U)
        c.sat_done = 2
      end
      while idx < 20 && idx > 1
        @vprint :ClassGroup 1 "Finishing by saturating up to $idx\n"
        fl = false
        p = 2
        while !fl && p < 2*Int(idx)
          fl = saturate!(c, U, p, stable)
          p = next_prime(p)
        end
        #fl = any(p->saturate!(c, U, p, stable), PrimesSet(1, 2*Int(idx)))
        @assert fl  # so I can switch assertions off...
        c.sat_done = 2*Int(idx)
        @vtime_add_elapsed :UnitGroup 1 c :unit_time r = _unit_group_find_units(U, c)
        n_idx = _validate_class_unit_group(c, U)
        @vprint :ClassGroup 1 "index estimate down to $n_idx from $idx\n"
        @assert idx != n_idx
        idx = n_idx
      end
      if idx == 1
        break
      end
    end
    #TODO: use LLL?
    if need_more
      d = root(abs(discriminant(O)), 2)
      c.expect = class_group_expected(d, degree(O), Int(norm(c.FB.ideals[1])), 100)
      need_more = false
    end
    h_old = class_group_current_h(c)
    class_group_new_relations_via_lll(c, extra = unit_rank(O) - length(U.units) +1)
    if h_old == class_group_current_h(c)
      do_units = true
    else
      do_units = false
    end
  end
  @assert U.full_rank
  _set_UnitGrpCtx_of_order(O, U)

  c.finished = true
  U.finished = true

  @vprint :ClassGroup 1 "hnftime $(c.time[:hnf_time])\n"

  if !GRH
    class_group_proof(c, fmpz(2), factor_base_bound_minkowski(O))
    for (p, _) in factor(c.h)
      while saturate!(c, U, Int(p), 3.5)
      end
    end
  end

  return c, U, _validate_class_unit_group(c, U)
end

function unit_group_ctx(c::ClassGrpCtx; redo::Bool = false)
  O = order(c.FB.ideals[1])
  if !redo
    try
      U = _get_UnitGrpCtx_of_order(O)::UnitGrpCtx
      return U
    catch e
    end
  end

  U = UnitGrpCtx{FacElem{nf_elem, AnticNumberField}}(O)
  need_more = true
  while true
    r = _unit_group_find_units(U, c)
    if r == 0
      if need_more
        d = root(abs(discriminant(O)), 2)
        c.expect = class_group_expected(d, degree(O), Int(norm(c.FB.ideals[1])), 100)
        need_more = false
      end
      class_group_new_relations_via_lll(c, E)
    else
      _set_UnitGrpCtx_of_order(O, U)
      return U
    end
  end
end

function unit_group(c::ClassGrpCtx)
  U = unit_group_ctx(c)
  return unit_group(c, U)
end

function unit_group(c::ClassGrpCtx, U::UnitGrpCtx)
  O = order(c.FB.ideals[1])
  K = nf(O)
  U, mU = unit_group_fac_elem(c, U)

  r = MapUnitGrp{typeof(O)}()
  r.header = Hecke.MapHeader(U, O,
    x->O(evaluate(image(mU, x))),
    x->preimage(mU, FacElem([K(x)], fmpz[1])))
  return U, r
end

@doc Markdown.doc"""
    class_group(O::NfOrd; bound = -1, method = 3, redo = false, large = 1000) -> GrpAbFinGen, Map

Returns a group $A$ and a map $f$ from $A$ to the set of ideals of $O$.
The inverse of the map is the projection onto the group of ideals modulo the
group of principal ideals.
\texttt{redo} allows to trigger a re-computation, thus avoiding the cache.
\texttt{bound}, when given, is the bound for the factor base.
"""
function class_group(O::NfOrd; bound::Int = -1, method::Int = 3,
                     redo::Bool = false, unit_method::Int = 1,
                     large::Int = 1000, use_aut::Bool = false, GRH::Bool = true)
  c, U, b = _class_unit_group(O, bound = bound, method = method, redo = redo, unit_method = unit_method, large = large, use_aut = use_aut, GRH = GRH)

  @assert b == 1
  return class_group(c)
end

@doc Markdown.doc"""
    unit_group(O::NfOrd) -> GrpAbFinGen, Map

Returns a group $U$ and an isomorphism map $f \colon U \to \mathcal O^\times$.
A set of fundamental units of $\mathcal O$ can be
obtained via `[ f(U[1+i]) for i in 1:unit_rank(O) ]`.
`f(U[1])` will give a generator for the torsion subgroup.
"""
function unit_group(O::NfOrd; method::Int = 3, unit_method::Int = 1, use_aut::Bool = false, GRH::Bool = true)
  c, U, b = _class_unit_group(O, method = method, unit_method = unit_method, use_aut = use_aut, GRH = GRH)
  @assert b==1
  return unit_group(c, U)
end

@doc Markdown.doc"""
    unit_group_fac_elem(O::NfOrd) -> GrpAbFinGen, Map

Returns a group $U$ and an isomorphism map $f \colon U \to \mathcal O^\times$.
A set of fundamental units of $\mathcal O$ can be
obtained via `[ f(U[1+i]) for i in 1:unit_rank(O) ]`.
`f(U[1])` will give a generator for the torsion subgroup.
All elements will be returned in factored form.
"""
function unit_group_fac_elem(O::NfOrd; method::Int = 3, unit_method::Int = 1, use_aut::Bool = false, GRH::Bool = true)
  c, U, b = _class_unit_group(O, method = method, unit_method = unit_method, use_aut = use_aut, GRH = GRH)
  @assert b==1
  return unit_group_fac_elem(c, U)
end

@doc Markdown.doc"""
    regulator(O::NfOrd)
Computes the regulator of $O$, ie. the discriminant of the unit lattice.
"""
function regulator(O::NfOrd; method::Int = 3, unit_method::Int = 1, use_aut::Bool = false, GRH::Bool = true)
  c, U, b = _class_unit_group(O, method = method, unit_method = unit_method, use_aut = use_aut, GRH = GRH)
  @assert b==1
  unit_group_fac_elem(c, U)
  return U.tentative_regulator
end

@doc Markdown.doc"""
    regulator(K::AnticNumberField)
Computes the regulator of $K$, ie. the discriminant of the unit lattice
for the maximal order of $K$
"""
function regulator(K::AnticNumberField)
  return regulator(maximal_order(K))
end


text:
module PyPlot2Igor

import PyCall: PyObject
import PyPlot: gcf
import ProgressMeter: Progress, update!

export pyplot2igor

abstract type PlotObject end

struct Scatter <: PlotObject
    o::PyObject
end

struct Line <: PlotObject
    o::PyObject
end

function write_wave(filename::String, array::AbstractVector, name::String)
    open(filename, "a") do f
        print(f, "WAVES/D $(name)\n")
        print(f, "BEGIN\n")
        for element in array
            print(f, "$(element)\n")
        end
        print(f, "END\n")
    end
end

function write_command(filename, command)
    open(filename, "a") do f
        print(f, "X $(command)\n")
    end
end

function itx_file(filename)
    open(filename, "w+") do f
        print(f, "IGOR\n")
    end
end

"""
    pyplot2igor(source_file::String, target_file::String) -> nothing

Scrapes the data from the figures in `source_file` and saves them
as IGOR text format in `target_file`. The macro looks for PyPlot's
`savefig` command. Only figures that get saved by the original script
will get exported.

# Example
```
@pyplot2igor "my_plotting_script.jl" "igor_text.itx"
```
"""
function pyplot2igor(source_file::String, target_file::String)
    fig_counter = 0

    ex_toplevel = Meta.parseall(read(source_file) |> String)

    # these are all the expressions
    exs = ex_toplevel.args

    PyPlot2Igor.itx_file(target_file)

    # initialize progress meter
    p = Progress(length(exs); dt=0.5, desc="Evaluating Expressions... ")

    for (line_number, ex) in enumerate(exs)
        if ex isa LineNumberNode
            continue
        elseif ex.head == :call && has_argument(ex, :savefig)
            fig_counter += 1
            write_plots(target_file, gcf(), fig_counter)
        else
            eval(ex)
        end

        # update progress meter
        update!(p, line_number)
    end
end

function has_argument(ex, search_argument)
    found = false

    for arg in  ex.args
        if arg isa Symbol
            arg == :savefig && (found = true)
        elseif arg isa QuoteNode
            arg.value == :savefig && (found = true)
        elseif arg isa Expr
            has_argument(arg, search_argument) && (found = true)
        end
    end

    found
end

isline(o::PyObject) = hasproperty(o, :get_data)
isscatter(o::PyObject) = hasproperty(o, :get_offsets)

function parse_to_type(obj::PyObject)
    if isline(obj)
        return Line(obj)
    elseif isscatter(obj)
        return Scatter(obj)
    else
        # @info("Could not match type of $obj.")
    end
end

function get_data(l::Line)
    data = l.o.get_data()
    data[1], data[2]
end

function get_data(s::Scatter)
    data = s.o.get_offsets()
    data[:,1], data[:,2]
end

function write_plots(filename, fig, fig_counter)
    for ax in fig.axes
        for (i, child) in ax.get_children() |> enumerate
            typed_child = parse_to_type(child)

            isnothing(typed_child) && continue

            x, y =  get_data(typed_child)

            xwave_name = "F$(fig_counter)X$i"
            ywave_name = "F$(fig_counter)Y$i"

            PyPlot2Igor.write_wave(filename, x, xwave_name)
            PyPlot2Igor.write_wave(filename, y, ywave_name)

            if i == 1
                cmd = "Display $ywave_name vs $xwave_name"
            else
                cmd = "AppendToGraph $ywave_name vs $xwave_name"
            end

            PyPlot2Igor.write_command(filename, cmd)
        end
    end
end

end #module


text:
# For unparametrized destination types
generate_copyto!_signature(dest, dest_type::Symbol, Msig) =
    :(Base.copyto!($(dest)::$(dest_type), applied_obj::$(Msig)))

# For parametrized destination types
function generate_copyto!_signature(dest, dest_type::Expr, Msig)
    dest_type.head == :curly ||
        throw(ArgumentError("Invalid destination specification $(dest)::$(dest_type)"))
    :(Base.copyto!($(dest)::$(dest_type), applied_obj::$(Msig)) where {$(dest_type.args[2:end]...)})
end

function generate_copyto!(body, factor_names, Msig)
    body.head == :(->) ||
        throw(ArgumentError("Invalid copyto! specification"))
    body.args[1].head == :(::) ||
        throw(ArgumentError("Invalid destination specification $(body.args[1])"))
    (dest,dest_type) = body.args[1].args
    copyto!_signature = generate_copyto!_signature(dest, dest_type, Msig)
    f_body = quote
        axes($dest) == axes(applied_obj) || throw(DimensionMismatch("axes must be same"))
        $(factor_names) = applied_obj.args
        $(body.args[2].args...)
        $(dest)
    end
    Expr(:function, copyto!_signature, f_body)
end

"""
    @materialize function op(args...)

This macro simplifies the setup of a few functions necessary for the
materialization of [`Applied`](@ref) objects:

- `ApplyStyle`, used to ensure dispatch of the applied object to the
  routines below

- `copyto!(dest::DestType, applied_obj::Applied{...,op})` performs the
  actual materialization of `applied_obj` into the destination object
  which has been generated by

- `similar` which usually returns a suitable matrix

- `materialize` which makes use of the above functions

# Example

```julia
@materialize function *(Ac::MyAdjointBasis,
                        O::MyOperator,
                        B::MyBasis)
    MyApplyStyle # An instance of this type will be returned by ApplyStyle
    T -> begin # generates similar
        A = parent(Ac)
        parent(A) == parent(B) ||
            throw(ArgumentError("Incompatible bases"))

        # There may be different matrices best representing different
        # situations:
        if ...
            Diagonal(Vector{T}(undef, size(B,1)))
        else
            Tridiagonal(Vector{T}(undef, size(B,1)-1),
                        Vector{T}(undef, size(B,1)),
                        Vector{T}(undef, size(B,1)-1))
        end
    end
    dest::Diagonal{T} -> begin # generate copyto!(dest::Diagonal{T}, ...) where T
        dest.diag .= 1
    end
    dest::Tridiagonal{T} -> begin # generate copyto!(dest::Tridiagonal{T}, ...) where T
        dest.dl .= -2
        dest.ev .= 1
        dest.du .= 3
    end
end
```
"""
macro materialize(expr)
    expr.head == :function || expr.head == :(=) || error("Must start with a function")
    @assert expr.args[1].head == :call
    op = expr.args[1].args[1]

    bodies = filter(e -> !(e isa LineNumberNode), expr.args[2].args)
    length(bodies) < 3 &&
        throw(ArgumentError("At least three blocks required (ApplyStyle, similar, and at least one copyto!)"))

    factor_types = :(<:Tuple{})
    factor_names = :(())
    apply_style = first(bodies)
    apply_style_fun = :(LazyArrays.ApplyStyle(::typeof($op)) = $(apply_style)())

    # Generate Applied signature
    for arg in expr.args[1].args[2:end]
        arg isa Expr && arg.head == :(::) ||
            throw(ArgumentError("Invalid argument specification $(arg)"))
        arg_name, arg_typ = arg.args
        push!(factor_types.args[1].args, :(<:$(arg_typ)))
        push!(factor_names.args, arg_name)
        push!(apply_style_fun.args[1].args, :(::Type{<:$(arg_typ)}))
    end
    Msig = :(LazyArrays.Applied{$(apply_style), typeof($op), $(factor_types)})

    sim_body = bodies[2]
    sim_body.head == :(->) ||
        throw(ArgumentError("Invalid similar specification"))
    T = first(sim_body.args)

    copytos! = map(body -> generate_copyto!(body, factor_names, Msig), bodies[3:end])

    f = quote
        $(apply_style_fun)

        function Base.similar(applied_obj::$Msig, ::Type{$T}=eltype(applied_obj)) where $T
            $(factor_names) = applied_obj.args
            $(sim_body.args[2])
        end

        $(copytos!...)

        LazyArrays.materialize(applied_obj::$Msig) =
            copyto!(similar(applied_obj, eltype(applied_obj)), applied_obj)
    end
    esc(f)
end


text:
# The main parser API.
#
# This is defined separately from parser.jl so that:
# * parser.jl doesn't need to refer to any tree data structures
# * It's clear which parts are the public API
#
# What should the general parsing API look like? Some points to consider:
#
# * After parsing atoms or statements or most other internal rules, it's
#   usual to start in the middle of the input text and end somewhere else in
#   the middle of the input text. So we should taken an index for the start of
#   parsing and supply an index back to the caller after parsing.
#
# * `parseall` is a special case where we expect to consume all the input.
#   Perhaps this is the API which throws an error if we don't consume it all,
#   and doesn't accept an index as input?
#
# * The ParseStream is the fundamental interface which wraps the code string
#   and index up together for input and contains the output events, diagnostics
#   and current stream position after parsing. The user should potentially be
#   able to use this directly. It does, however assume a Julia-compatible token
#   stream.
#
# * It could be useful to support an IO-based interface so that users can parse
#   Julia code intermixed with other DSLs. Documenter.jl and string macros come
#   to mind as examples which could use this. A tricky part is deciding where
#   the input ends: For string macros this is done by the parser, but for
#   Documenter it's probably just done beforehand according to the Markdown
#   code block rules.
#
# * The API should have an interface where a simple string is passed in. How
#   does SourceFile relate to this?
#
# * It's neat for `parse` to be overloadable to produce various output data
#   structures; GreenNode, SyntaxNode, Expr, (etc?) in the same way that
#   Base.parse can be used for non-Julia code.  (Heh... though
#   `Base.parse(Expr, "...")` would also make a certain amount of sense.)
#
# * What's the no-copy API look like? A String can be put into an IOBuffer via
#   unsafe_wrap(Vector{UInt8}, str) ... A SubString likewise. Also there's the
#   `codeunits` function to hold a GC-safe view of string data as an array (but
#   we can't use a Vector{UInt8})

struct ParseError <: Exception
    source::SourceFile
    diagnostics::Vector{Diagnostic}
end

function Base.showerror(io::IO, err::ParseError, bt; backtrace=false)
    println(io, "ParseError:")
    show_diagnostics(io, err.diagnostics, err.source)
end

function Base.showerror(io::IO, err::ParseError)
    println(io, "ParseError:")
    show_diagnostics(io, err.diagnostics, err.source)
end

Base.display_error(io::IO, err::ParseError, bt) = Base.showerror(io, err, bt)

"""
    # Input and output:
    stream = parse(stream::ParseStream; kws...)
    (tree, diagnostics)        = parse(TreeType, io::IOBuffer; kws...)
    (tree, diagnostics, index) = parse(TreeType, str::AbstractString, [index::Integer]; kws...)
    # Keywords
    parse(...; rule=:toplevel, version=VERSION, ignore_trivia=true)

Parse Julia source code from `input`, returning the output in a format
compatible with `input`:

* When `input` is a `ParseStream`, the stream itself is returned and the
  `ParseStream` interface can be used to process the output.
* When `input` is a seekable `IO` subtype, the output is `(tree, diagnostics)`.
  The buffer `position` will be set to the next byte of input.
* When `input` is an `AbstractString, Integer`, or `Vector{UInt8}, Integer` the
  output is `(tree, diagnostics, index)`, where `index` (default 1) is the next
  byte of input.

`rule` may be any of
* `toplevel` (default) — parse a whole "file" of top level statements. In this
  mode, the parser expects to fully consume the input.
* `statement` — parse a single statement, or statements separated by semicolons.
* `atom` — parse a single syntax "atom": a literal, identifier, or
  parenthesized expression.

`version` (default `VERSION`) may be used to set the syntax version to
any Julia version `>= v"1.0"`. We aim to parse all Julia syntax which has been
added after v"1.0", emitting an error if it's not compatible with the requested
`version`.

See also [`parseall`](@ref) for a simpler but less powerful interface.
"""
function parse(stream::ParseStream; rule::Symbol=:toplevel)
    ps = ParseState(stream)
    if rule === :toplevel
        parse_toplevel(ps)
    elseif rule === :statement
        parse_stmts(ps)
    elseif rule === :atom
        parse_atom(ps)
    else
        throw(ArgumentError("Unknown grammar rule $rule"))
    end
    stream
end

function parse(::Type{T}, io::IO;
               rule::Symbol=:toplevel, version=VERSION, kws...) where {T}
    stream = ParseStream(io; version=version)
    parse(stream; rule=rule)
    tree = build_tree(T, stream; kws...)
    seek(io, last_byte(stream))
    tree, stream.diagnostics
end

# Generic version of parse for all other cases where an index must be passed
# back - ie strings and buffers
function parse(::Type{T}, input...;
               rule::Symbol=:toplevel, version=VERSION, kws...) where {T}
    stream = ParseStream(input...; version=version)
    parse(stream; rule=rule)
    tree = build_tree(T, stream; kws...)
    tree, stream.diagnostics, last_byte(stream) + 1
end

"""
    parseall(TreeType, input...;
             rule=:toplevel,
             version=VERSION,
             ignore_trivia=true)

Experimental convenience interface to parse `input` as Julia code, emitting an
error if the entire input is not consumed. `input` can be a string or any other
valid input to the `ParseStream` constructor. By default `parseall` will ignore
whitespace and comments before and after valid code but you can turn this off
by setting `ignore_trivia=false`.

A `ParseError` will be thrown if any errors occurred during parsing.

See [`parse`](@ref) for a more complete and powerful interface to the parser,
as well as a description of the `version` and `rule` keywords.
"""
function parseall(::Type{T}, input...; rule=:toplevel, version=VERSION,
                  ignore_trivia=true, filename=nothing) where {T}
    stream = ParseStream(input...; version=version)
    if ignore_trivia && rule != :toplevel
        bump_trivia(stream, skip_newlines=true)
        empty!(stream)
    end
    parse(stream; rule=rule)
    if (ignore_trivia  && peek(stream, skip_newlines=true) != K"EndMarker") ||
       (!ignore_trivia && (peek(stream, skip_newlines=false, skip_whitespace=false) != K"EndMarker"))
        emit_diagnostic(stream, error="unexpected text after parsing $rule")
    end
    if any_error(stream.diagnostics)
        source = SourceFile(sourcetext(stream, steal_textbuf=true), filename=filename)
        throw(ParseError(source, stream.diagnostics))
    end
    # TODO: Figure out a more satisfying solution to the wrap_toplevel_as_kind
    # mess that we've got here.
    # * It's kind of required for GreenNode, as GreenNode only records spans,
    #   not absolute positions.
    # * Dropping it would be ok for SyntaxNode and Expr...
    tree = build_tree(T, stream; wrap_toplevel_as_kind=K"toplevel", filename=filename)
    if !isempty(stream.diagnostics)
        # Crudely format any warnings to the current logger.
        buf = IOBuffer()
        show_diagnostics(IOContext(buf, stdout), stream,
                         SourceFile(sourcetext(stream, steal_textbuf=true), filename=filename))
        @warn Text(String(take!(buf)))
    end
    tree
end


text:
using SafeTestsets

module GyroKinetics4dTests

    using Test
    using GeometricIntegrators
    using ChargedParticleDynamics.GyroKinetics4d

    const Δt = 1.
    const nt = 1

    const nl = 100
    const nx = 10
    const ny = 10

    export test_gyro_kinetics_4d_erk4, test_gyro_kinetics_4d_glrk, test_gyro_kinetics_4d_strang
    export nl, nx, ny

    function test_gyro_kinetics_4d_erk4(ode, ωabs; Δt = Δt)
        int = Integrator(ode, getTableauERK4(), Δt)
        sol = integrate(ode, int, nt)
        @test true
    end

    function test_gyro_kinetics_4d_glrk(ode, ωabs; Δt = Δt)
        int = Integrator(ode, getTableauGLRK(1), Δt)
        sol = integrate(ode, int, nt)
        @test true
    end

    function test_gyro_kinetics_4d_strang(ode, ωabs; Δt = Δt)
        tab = getTableauStrang()
        mpi = Tuple(IntegratorConstructor(eltype(ode.q₀), ndims(ode), getTableauGLRK(1)) for i in 1:6)
        int = IntegratorComposition(ode, mpi, tab, Δt)
        sol = integrate(ode, int, nt)
        @test true
    end

end

@safetestset "Gyrokinetic GC Model in 4D with ITER-like Solov'ev Equilibrium with X-Point                         " begin

    using GeometricIntegrators.Solvers
    using ChargedParticleDynamics.GyroKinetics4d.GuidingCenter4dSolovevIterXpoint
    # using ChargedParticleDynamics.GyroKinetics4d.GuidingCenter4dSolovevIterXpoint: scaling_factor, transform_q̃_to_q_jacobian!, transform_q̃_to_q_rhs!
    using ..GyroKinetics4dTests

    # test coordinate transformation
    q̃, params = initial_conditions_trapped()
    t  = 0.
    @test transform_q_to_q̃(t, transform_q̃_to_q(t, q̃, params), params) ≈ q̃ atol=1E-6

    # test ODE
    test_gyro_kinetics_4d_erk4(guiding_center_4d_ode(initial_conditions_trapped()...), ωabs; Δt=1.)
    test_gyro_kinetics_4d_glrk(guiding_center_4d_ode(initial_conditions_trapped()...), ωabs; Δt=1.)
    # test_gyro_kinetics_4d_glrk(guiding_center_4d_ode(initial_conditions_barely_passing()...), Δt=1.)
    # test_gyro_kinetics_4d_glrk(guiding_center_4d_ode(initial_conditions_barely_trapped()...), Δt=1.)
    # test_gyro_kinetics_4d_glrk(guiding_center_4d_ode(initial_conditions_deeply_passing()...), Δt=1.)
    # test_gyro_kinetics_4d_glrk(guiding_center_4d_ode(initial_conditions_deeply_trapped()...), Δt=1.)

    # test SODE
    test_gyro_kinetics_4d_strang(guiding_center_4d_sode(initial_conditions_trapped()...), ωabs; Δt=1.)
    # test_gyro_kinetics_4d_strang(guiding_center_4d_sode(initial_conditions_barely_passing()...), Δt=1.)
    # test_gyro_kinetics_4d_strang(guiding_center_4d_sode(initial_conditions_barely_trapped()...), Δt=1.)
    # test_gyro_kinetics_4d_strang(guiding_center_4d_sode(initial_conditions_deeply_passing()...), Δt=1.)
    # test_gyro_kinetics_4d_strang(guiding_center_4d_sode(initial_conditions_deeply_trapped()...), Δt=1.)

end


text:
# Core Ops

function raft_execute(op, func, argv)
    return func(argv)
end

function raft_init_node(nodeId, clusterId)
    db_put(log, "NodeId", (nodeId, clusterId))
    nodeId
end

function raft_get_node()
    db_get(log, "NodeId")
end

"""
    raft_init_cluster()
Cluster bootstrap. Insert an entry in log with clusterid, node id, ipadress and
port of this node. When new nodes join this entry will be replicated to other
nodes.
"""
function raft_init_cluster(args)
    clusterId != nothing && throw(RavanaException("This node is part of cluster " * clusterId))
    thisNode = nothing
    (arg) = args

    global nodeId    = rand(UInt128)
    global clusterId = rand(UInt128)
    global nodes = raftNode[]
    thisNode = raftNode(ipAddress, ipPort, nodeId, clusterId)
    push!(nodes, thisNode)
    arg = nodes

    raft_init_node(nodeId, clusterId)
    db_put(log, 1, (1, OP_INIT_CLUSTER, arg))
    db_put(log, "RavanaCluster", arg)

    global rInfo = raftInfo(1, nodeId, 1, 1, 1)
    raft_save_info(rInfo)
    raft_set_state(FOLLOWER)

    return arg
end

raft_noop() = return

function raft_cluster_add_node(thisNode, clusterNode)

end

function raft_cluster_get_config(args)
    @debug "in raft_cluster_get_config()"
    return nodes
end

# Ravana ops
const OP_INIT_CLUSTER            = Int32(1)
const OP_NOOP                    = Int32(2)
const OP_GET_CLUSTER_CONF        = Int32(3)
const OP_ADD_NODE                = Int32(4)

const OP_UNKNOWN                 = Int32(10000)

# Lookup table for Core Ravana Ops
const op_table = Dict(OP_INIT_CLUSTER      =>  raft_init_cluster,
                      OP_NOOP              =>  raft_noop,
                      OP_GET_CLUSTER_CONF  =>  raft_cluster_get_config,
                      OP_ADD_NODE          =>  raft_cluster_add_node)


text:
F81(πA::F, πC::F, πG::F, πT::F;
    safe::Bool=true) where F <: Float64 =
F81rel(πA, πC, πG, πT)

F81(β::F,
    πA::F, πC::F, πG::F, πT::F;
    safe::Bool=true) where F <: Float64 =
F81abs(β, πA, πC, πG, πT, safe=safe)

function F81(θ_vec::A,
             π_vec::A;
             safe::Bool=true) where A <: AbstractArray
  if safe && length(π_vec) != 4
    error("Incorrect base frequency vector length")
  end
  if length(θ_vec) == 0
    return F81rel(π_vec[DNA_A], π_vec[DNA_C], π_vec[DNA_G], π_vec[DNA_T],
                  safe=safe)
  elseif length(θ_vec) == 1
    return F81abs(θ_vec[1],
                  π_vec[DNA_A], π_vec[DNA_C], π_vec[DNA_G], π_vec[DNA_T],
                  safe = safe)
  else
    error("Parameter vector length incompatiable with absolute or relative rate form of substitution model")
  end
end

function F81rel(θ_vec::A,
                π_vec::A;
                safe::Bool=true) where A <: AbstractArray
  if safe
    if length(θ_vec) != 0
      error("Incorrect parameter vector length")
    elseif length(π_vec) != 4
      error("Incorrect base frequency vector length")
    end
  end
  return F81rel(π_vec[DNA_A], π_vec[DNA_C], π_vec[DNA_G], π_vec[DNA_T], safe=safe)
end

function F81abs(θ_vec::A,
                π_vec::A;
                safe::Bool=true) where A <: AbstractArray
  if safe
    if length(θ_vec) != 1
      error("Incorrect parameter vector length")
    elseif length(π_vec) != 4
      error("Incorrect base frequency vector length")
    end
  end
  return F81abs(θ_vec[1], π_vec[DNA_A], π_vec[DNA_C], π_vec[DNA_G], π_vec[DNA_T], safe=safe)
end


text:
using Test, CompGrids, ParallelStencil, MPI, ImplicitGlobalGrid

@testset "ParallelStencil grids" begin
    size_t = (32,32,32)
    length = (1, 2, 3)

    for idim=1:3
        for mpi in (false, true)

            if mpi
                if idim==1
                    ParallelStencil.@reset_parallel_stencil()
                    @init_backend(ParallelStencil, Threads, true,  1);
                elseif idim==2
                    ParallelStencil.@reset_parallel_stencil()
                    @init_backend(ParallelStencil, Threads, true,  2);
                elseif idim==3
                    ParallelStencil.@reset_parallel_stencil()
                    @init_backend(ParallelStencil, Threads, true,  3);
                end
            else
                if idim==1
                    ParallelStencil.@reset_parallel_stencil()
                    @init_backend(ParallelStencil, Threads, false, 1);
                elseif idim==2
                    ParallelStencil.@reset_parallel_stencil()
                    @init_backend(ParallelStencil, Threads, false, 2);
                elseif idim==3
                    ParallelStencil.@reset_parallel_stencil()
                    @init_backend(ParallelStencil, Threads, false, 3);
                end
            end
            grid = RegularRectilinearCollocatedGrid(size=(size_t[1:idim]...,), extent=length[1:idim])
            @test grid.L[1] == 1.0
            if idim>1
                @test grid.L[2] == 2.0
            end
            if idim>2
                @test grid.L[3] == 3.0
            end
        end
    end

    # Specify with x-coordinates
    grid = RegularRectilinearCollocatedGrid(size=(10,), x=(-5, 12))
    @test grid.L[1] == 17.0

    # Specify as single numbers (instead of Tuple)
    grid = RegularRectilinearCollocatedGrid(size=10, extent=100.0)
    @test grid.L[1] == 100.0

    # Create fields on the grid as well
    ParallelStencil.@reset_parallel_stencil()
    @init_backend(ParallelStencil, Threads, true, 2);
    grid = RegularRectilinearCollocatedGrid(size=(10,12), extent=(100.0,50), fields=(T=0,P=11.22))
    @test grid.fields[:P][1] == 11.22

    # Different methods to set values (using idx is identical to the PETSc implementation)
    grid.fields.P .= 0.0
    for i in grid.corners.lower:grid.corners.upper
        grid.fields.T[i] = i[1]
        grid.fields.P[grid.idx(i)] = i[1]
    end

    # test on one core
    @test sum(abs.(grid.fields.P - grid.fields.T)) == 0.0

end


text:
using SparseArrays
using LinearAlgebra
using Plots
using RecipesBase

include("C:/Users/lucas/OneDrive/Documents/GitHub/GDB.jl/src/domain/grid.jl")
include("C:/Users/lucas/OneDrive/Documents/GitHub/GDB.jl/src/domain/barrier.jl")
include("C:/Users/lucas/OneDrive/Documents/GitHub/GDB.jl/src/domain/operators.jl")
include("C:/Users/lucas/OneDrive/Documents/GitHub/GDB.jl/src/domain/interpolation.jl")
include("C:/Users/lucas/OneDrive/Documents/GitHub/GDB.jl/experimental/LBM.jl")
include("C:/Users/lucas/OneDrive/Documents/GitHub/GDB.jl/experimental/periodic_derivatives.jl")
include("C:/Users/lucas/OneDrive/Documents/Dartmouth/Research/Code/LBMEigen.jl")

function thermal_lattice_boltzmann_method(streamtensor, Dx, Dy, Dxx, Dyy, xres, yres, f, g, tau, tauc, deltax, deltat)
    # vector at each point with nine dimensions corresponding to each direction 0-8
    c = deltax/deltat
    nu = tau*c^2/3
    # xres = grid._Nx
    # yres = grid._Ny

    fnew = streamtensor * f
    gnew = streamtensor * g
    # println(size(ftransform))

    rho = zeros((xres*yres))
    for i=1:xres*yres
        rho[i] = sum(fnew[i,:])
    end

    p = zeros(xres*yres)
    for i=1:xres*yres
        p[i] = c^2/3 * rho[i]
    end

    u = zeros(xres*yres, 2)
    for i=1:xres*yres
        ux = (fnew[i,2] - fnew[i,4] + 1/sqrt(2) * (fnew[i,6] + fnew[i,9] - fnew[i,7] - fnew[i,8])) / rho[i]
        uy = (fnew[i,3] - fnew[i,5] + 1/sqrt(2) * (fnew[i,6] + fnew[i,8] - fnew[i,9] - fnew[i,8])) / rho[i]
        u[i,1] = ux
        u[i,2] = uy
    end

    # derivative matrices are operators that act on other matrices

    q = zeros(xres*yres, 9)

    grad = [Dx, Dy]
    divu = Dx*u[:,1] + Dy*u[:,2]
    gradp = zeros(xres*yres, 2)
    gradp[:,1] = Dx*p[:]
    gradp[:,2] = Dy*p[:]
    laplu = zeros(xres*yres, 2)
    laplu[:,1] = (Dxx+Dyy)*u[:,1]
    laplu[:,2] = (Dxx+Dyy)*u[:,2]

    # unitvecs = [[0,0] [c,0] [0,c] [-c,0] [0,-c] 1/sqrt(2)*[c,c] 1/sqrt(2)*[-c,c] 1/sqrt(2)*[-c,-c] 1/sqrt(2)*[c,-c]]
    unitvecs = zeros(9,2)
    unitvecs[:,1] = [0, c, 0, -c, 0, 1/sqrt(2)*c, -1/sqrt(2)*c, -1/sqrt(2)*c, 1/sqrt(2)*c]
    unitvecs[:,2] = [0, 0, c, 0, -c, 1/sqrt(2)*c, 1/sqrt(2)*c, -1/sqrt(2)*c, -1/sqrt(2)*c]

    vgradu = zeros(xres*yres, 9, 2)
    for k=1:9
        vgradu[:,k,1] = (unitvecs[k,1] .- u[:,1]) .* (Dx*u[:,1]) + (unitvecs[k,2] .- u[:,2]) .* (Dy*u[:,1])
        vgradu[:,k,2] = (unitvecs[k,1] .- u[:,1]) .* (Dx*u[:,2]) + (unitvecs[k,2] .- u[:,2]) .* (Dy*u[:,2])
    end

    # .* for element wise multiplication, similar for division
    for i=1:xres*yres
        for a=1:9
            q[i,a] = - 1/rho[i]*dot(unitvecs[a,:]-u[i,:], gradp[i,:])
            + nu*(dot(unitvecs[a,:]-u[i,:], laplu[i,:]))
            + dot(unitvecs[a,:]-u[i,:], vgradu[i,a,:])
        end
    end

    w = zeros(9)
    w[1] = 4/9
    w[2] = 1/9
    w[3] = 1/9
    w[4] = 1/9
    w[5] = 1/9
    w[6] = 1/36
    w[7] = 1/36
    w[8] = 1/36
    w[9] = 1/36

    s = zeros((xres*yres, 9))
    usquared = zeros((xres*yres))
    for i=1:xres*yres
        usquared[i] = dot(u[i,:],u[i,:])
        s[i,:] = w[:] .* (3*(unitvecs*u[i,:])/c^2 .+ 9/2*((unitvecs*u[i,:])).^2/c^3 .- 3/2*usquared[i]/c^2)
    end

    feq = zeros((xres*yres, 9))
    for i=1:xres*yres
        feq[i,:] = w[:].*rho[i] .+ rho[i].*s[i,:]
    end

    fnext = zeros((xres*yres, 9))
    for i=1:xres*yres
        for k=1:9
            # fnext[i,j,k] = 1/(1 - 1/tau) * (fnew[i,j,k] - 1/tau*feq[i,j,k])
            # fnext[i,j,k] = fnew[i,j,k] - (1/tau)*(fnew[i,j,k]-feq[i,j,k])
            fnext[i,:] = fnew[i,:] .- (.5*deltat/(tau+.5*deltat)).*(fnew[i,:].-feq[i,:])
        end
    end

    rhoepsilon = zeros(xres*yres)
    for i=1:xres*yres
        rhoepsilon[i] = sum(gnew[i,:] .- deltat/2 .*fnext[i,:].*q[i,:])
    end

    geq = zeros((xres*yres, 9))
    for i=1:xres*yres
        geq[i,1] = -2/3*rhoepsilon[i] * usquared[i]/c^2
        geq[i,2] = rhoepsilon[i]/9 * (1.5 + 1.5*dot(unitvecs[2,:],u[i,:])/c^2 + 4.5*(dot(unitvecs[2,:],u[i,:]))^2/c^4 - 1.5*usquared[i]/c^2)
        geq[i,3] = rhoepsilon[i]/9 * (1.5 + 1.5*dot(unitvecs[3,:],u[i,:])/c^2 + 4.5*(dot(unitvecs[3,:],u[i,:]))^2/c^4 - 1.5*usquared[i]/c^2)
        geq[i,4] = rhoepsilon[i]/9 * (1.5 + 1.5*dot(unitvecs[4,:],u[i,:])/c^2 + 4.5*(dot(unitvecs[4,:],u[i,:]))^2/c^4 - 1.5*usquared[i]/c^2)
        geq[i,5] = rhoepsilon[i]/9 * (1.5 + 1.5*dot(unitvecs[5,:],u[i,:])/c^2 + 4.5*(dot(unitvecs[5,:],u[i,:]))^2/c^4 - 1.5*usquared[i]/c^2)
        geq[i,6] = rhoepsilon[i]/36 * (3 + 6*dot(unitvecs[6,:],u[i,:])/c^2 + 4.5*(dot(unitvecs[6,:],u[i,:]))^2/c^4 - 1.5*usquared[i]/c^2)
        geq[i,7] = rhoepsilon[i]/36 * (3 + 6*dot(unitvecs[7,:],u[i,:])/c^2 + 4.5*(dot(unitvecs[7,:],u[i,:]))^2/c^4 - 1.5*usquared[i]/c^2)
        geq[i,8] = rhoepsilon[i]/36 * (3 + 6*dot(unitvecs[8,:],u[i,:])/c^2 + 4.5*(dot(unitvecs[8,:],u[i,:]))^2/c^4 - 1.5*usquared[i]/c^2)
        geq[i,9] = rhoepsilon[i]/36 * (3 + 6*dot(unitvecs[9,:],u[i,:])/c^2 + 4.5*(dot(unitvecs[9,:],u[i,:]))^2/c^4 - 1.5*usquared[i]/c^2)
    end

    gnext = zeros(xres*yres, 9)
    for i=1:xres*yres
        gnext[i,:] = gnew[i,:] .- deltat/(tauc + .5*deltat).*(gnew[i,:] .- geq[i,:]) .- tauc/(tauc + .5*deltat)*deltat.*fnext[i,:].*q[i,:]
    end

    return fnext, gnext, q, rho, usquared
end

N = 100
tau = 1000
tauc = 1000
deltax = 1
deltat = 1

grd = Grid(100, 100, 1; Nbuffer=0) do

    r0 = Float64[-1,-1]
    r1 = Float64[1,1]
    inside(x,y) = true
    inside, r0, r1

end

# xres = grd._Nx
# yres = grd._Ny

xres = 100
yres = 100

stensor = PeriodicStreamTensor(xres, yres)

sx = 5
sy = 5
dx = 1
dy = 1
pDx = periodic_derivative(1, 0, sx, sy, xres, yres, dx, dy)
pDy = periodic_derivative(0, 1, sx, sy, xres, yres, dx, dy)
pDxx = periodic_derivative(2, 0, sx, sy, xres, yres, dx, dy)
pDyy = periodic_derivative(0, 2, sx, sy, xres, yres, dx, dy)

Minvx = stencil1d(sx)
Minvy = stencil1d(sy)
Dx = derivative_matrix(1, 0, Minvx, Minvy, grd)
Dy = derivative_matrix(0, 1, Minvx, Minvy, grd)
Dxx = derivative_matrix(2, 0, Minvx, Minvy, grd)
Dyy = derivative_matrix(0, 2, Minvx, Minvy, grd)

fsine = zeros(xres, yres, 9)
gsine = zeros(xres, yres, 9)
rhoinit = zeros(xres, yres)
uinit = zeros(xres, yres, 9)
for i=1:xres
    for j=1:yres
        for k=1:9
            fsine[i,j,k] = .5 + .4*sin(i*2*pi/(xres))*sin(j*2*pi/(yres))
        end
        rhoinit[i,j] = sum(fsine[i,j,:])
        uinit[i,j,:] = fsine[i,j,:] ./ rhoinit[i,j]
        for k=1:9
            gsine[i,j,k] = (1 - uinit[i,j,k])^2/2 * fsine[i,j,k]
        end
    end
end

fs2 = zeros(xres*yres, 9)
gs2 = zeros(xres*yres, 9)

for j=1:yres
    for i=1:xres
        fs2[(j-1)*xres+i,:] = fsine[i,j,:]
        gs2[(j-1)*xres+i,:] = gsine[i,j,:]
    end
end

rhofull = zeros(xres*yres)
anim = @animate for i=1:N
    fs2, gs2, ex, ex1, ex2 = thermal_lattice_boltzmann_method(stensor, pDx, pDy, pDxx, pDyy, xres, yres, fs2, gs2, tau, tauc, deltax, deltat)
    fplot = zeros(xres, yres)
    gplot = zeros(xres, yres)
    explot = zeros(xres, yres)
    exp1 = zeros(xres, yres)
    exp2 = zeros(xres, yres)
    ucheck = zeros(xres, yres)
    if i==1
        rhofull = ex2
    else
        rhofull = hcat(rhofull, ex2)
    end
    for i=1:xres
        for j=1:yres
            exp1[i,j] = ex1[(j-1)*xres+i]
            exp2[i,j] = ex2[(j-1)*xres+i]
            for k=1:9
                fplot[i,j] = fplot[i,j] + fs2[(j-1)*xres+i,k]
                gplot[i,j] = gplot[i,j] + gs2[(j-1)*xres+i,k]
                explot[i,j] = explot[i,j] + ex[(j-1)*xres+i,k]
                ucheck[i,j] = ucheck[i,j] + ex[(j-1)*xres+i,k]/ex1[(j-1)*xres+i]
            end
        end
    end
    println(fplot[Int32(xres/2),Int32(yres/2)])
    println(gplot[Int32(xres/2),Int32(yres/2)])
    heatmap(gplot, clims=(0,5))
end
# gif(anim, fps=10)

p = size(rhofull)[2]
covarmatrix = 1/(p-1) * rhofull * transpose(rhofull)
evals, evecs = eigs(covarmatrix, nev=20, which=:LM)

evecstransform = zeros(xres, yres)
for j=1:yres
    for i=1:xres
        evecstransform[i,j] = evecs[(j-1)*xres+i,3]
    end
end

heatmap(evecstransform)


text:
module Surprise
# package code goes here

using Persa
include("methods.jl")

end # module