text:
"""
    heterozygosity(data::PopObj, mode = "locus")
Calculate observed and expected heterozygosity in a `PopObj`
## Example
heterozygosity(nancycats(), "population" )
### Modes
- `"locus"` : heterozygosity per locus (default)
- `"sample"` or `"ind"` or `"individual"` : heterozygosity per individual/sample
- `"population"` or `"pop"` : heterozygosity per population (PopObj.samples.population)
"""
function heterozygosity(data::PopObj, mode = "locus")
    if mode == "locus"
        obs = het_observed(data)
        exp = het_expected(data)
        locinames = loci(data)
        return DataFrame(locus = locinames, het_obs = obs, het_exp = exp)

    elseif lowercase(mode) == "sample" || lowercase(mode) == "ind" || lowercase(mode) == "individual"
        return DataFrame(name = samples(data), het_obs = het_sample(data))

    elseif lowercase(mode) == "pop" || lowercase(mode) == "population"
        obs = het_population_obs(data)
        exp = het_population_exp(data)
        merged_het = merge(obs, exp)
        locinames = loci(data)

        # remake obs/exp distinction
        obs_pops = obs["pops_obs"] .* "_obs"
        exp_pops = exp["pops_exp"] .* "_exp"

        # interleave pop_obs/exp for column names
        out_df = []
        for (i,j) in zip(obs_pops, exp_pops)
            tmp = DataFrame(locus = locinames, het_obs = merged_het[i], het_exp = merged_het[j])
            push!(out_df, tmp)
        end
    end
    return Tuple(out_df)
end

const het = heterozygosity
const He = heterozygosity

"""
    het_expected(data::PopObj)
Calculate the expected heterozygosity for each locus in a `PopObj`. Returns an
array of heterozygosity values.
"""
function het_expected(data::PopObj)
    het_vals = Vector{Float64}()
    @inbounds for locus in eachcol(data.loci, false)
        # get allele freqs at locus and isolate the freqs
        a = allele_freq(locus)
        het = 1 - (values(a) .^ 2 |> sum)
        push!(het_vals, het)  # push the sum of hetz to the output array
    end
    return het_vals
end

"""
    het_observed(data::PopObj)
Calculate the observed heterozygosity for each locus in a `PopObj`. Returns a
vector of heterozygosity values.
"""
function het_observed(data::PopObj)
    # faster and cheaper matrix-to-matrix comparison
    if all(data.samples.ploidy .== 2) == true
        allele_1 = allele_matrix(data.loci, 1)
        allele_2 = allele_matrix(data.loci, 2)
        het_matrix = allele_1 .!= allele_2
        return map(eachcol(het_matrix)) do i
            mean(i |> skipmissing)
        end
    else
        het_vals = Vector{Float64}()
        for locus in eachcol(data.loci, false)
            a = geno_freq(locus)  # get genotype freqs at locus
            tmp = 0
            genos = keys(a) |> collect
            for geno in genos
                if length(unique(geno)) != 1
                    tmp += a[geno]     # if true, add freq to total in tmp
                end
            end
            push!(het_vals, tmp)
        end
        return het_vals
    end
end

"""
    het_population_exp(data::PopObj)
Return a `Dict` of the expected heterozygosity per population for each locus in
a `PopObj`
"""
function het_population_exp(data::PopObj)
    # get population order and store in its own dict key
    popid = unique(data.samples.population |> collect)
    fixed_popid = [replace(i, " " => "_") for i in popid]
    d = Dict()
    d["pops_exp"] = fixed_popid

    y = deepcopy(data.loci)
    insertcols!(y, 1, :population => data.samples.population)
    y_subdf = groupby(y[!, :], :population)
    for pop in y_subdf
        pop_het_vals = Vector{Union{Missing, Float64}}()
        for locus in eachcol(pop[!, :2:end], false)
            a = allele_freq(locus)
            n = length(locus |> skipmissing |> collect)
            het = 1 - (values(a) .^ 2 |> sum)
            het_adjusted = (n/(n-1)) * het
            push!(pop_het_vals, het_adjusted)  # push the sum of hetz to the output array
        end

        # add "_obs" for distinction vs "_exp"
        popname = replace(pop.population[1], " " => "_")*"_exp"
        d[popname] = pop_het_vals
    end
        return d
end

"""
    het_population_obs(data::PopObj)
Return a Dict of the observed heterozygosity per population for each locus in
a `PopObj`
"""
function het_population_obs(data::PopObj)
    # get population order and store in its own dict key
    popid = unique(data.samples.population |> collect)
    fixed_popid = [replace(i, " " => "_") for i in popid]
    d = Dict()
    d["pops_obs"] = fixed_popid

    y = deepcopy(data.loci)
    insertcols!(y, 1, :population => data.samples.population)
    y_subdf = groupby(y[!, :], :population)

    # faster and cheaper diploid version
    if all(data.samples.ploidy .== 2) == true
        for pop in y_subdf
            ploidy = unique(data.samples.ploidy)
            allele_1 = allele_matrix(select(pop, Not(:population)), 1)
            allele_2 = allele_matrix(select(pop, Not(:population)), 2)
            het_matrix = allele_1 .!= allele_2
            pop_het_vals = map(eachcol(het_matrix)) do i
                mean(i |> skipmissing)
            end

            # add "_obs" for distinction vs "_exp"
            popname = replace(pop.population[1], " " => "_")*"_obs"
            d[popname] = pop_het_vals
        end
            return d
    else
        for pop in y_subdf
            pop_het_vals = Vector{Union{Missing,Float64}}()
            for locus in eachcol(pop[!, :2:end], false)
                # get genotype freqs at locus
                a = geno_freq(locus)
                n = length(locus |> skipmissing |> collect)

                # add condition if the locus is missing from that subpop
                if a === missing
                    push!(pop_het_vals, missing)
                    continue
                end
                tmp = 0
                genos = keys(a) |> collect
                for geno in genos
                    if length(unique(geno)) != 1
                        tmp += a[geno]     # if true, add freq to total in tmp
                    end
                end
                het_adjusted = (n/(n-1)) * tmp
                push!(pop_het_vals, het_adjusted)
            end

            # add "_obs" for distinction vs "_exp"
            popname = replace(pop.population[1], " " => "_")*"_obs"
            d[popname] = pop_het_vals
        end
        return d
    end
end

"""
    het_sample(data::PopObj)
Calculate the observed heterozygosity for each individual in a `PopObj`. Returns
a vector of heterozygosity values corresponding to the order of samples in the PopObj.
"""
function het_sample(data::PopObj)
    # efficient diploid method
    if all(data.samples.ploidy .== 2) == true
        allele_1 = allele_matrix(data.loci, 1)
        allele_2 = allele_matrix(data.loci, 2)
        het_matrix = allele_1 .!= allele_2
        return map(eachrow(het_matrix)) do i
            mean(i |> skipmissing)
        end
    else
        # psuedo-transpose to format genotypes as vectors per individual
        by_ind = map(i -> get_sample_genotypes(data, i), samples(data))

        # calculate observed heterozygosity like for loci
        map(het_sample, by_ind)
    end
end

"""
    het_sample(individual::Vector{<:Union{Missing, Tuple{Vararg}}}))
Calculate the observed heterozygosity for an individual in a `PopObj`. Returns
an array of heterozygosity values.
"""
function het_sample(individual::Vector{<:Union{Missing, NTuple{N,<:Integer}}}) where N
    # calculate observed heterozygosity for an individual
    het = 0
    genos = skipmissing(individual) |> collect
    for geno in genos
        if length(unique(geno)) != 1
            het += 1
        end
    end
    return het/length(genos)
end

function het_sample(individual::Vector{<:Union{Missing, Tuple{Vararg}}})
    # calculate observed heterozygosity for an individual
    het = 0
    genos = skipmissing(individual) |> collect
    for geno in genos
        if length(unique(geno)) != 1
            het += 1
        end
    end
    return het/length(genos)
end

"""
    hwe_test(data::PopObj; by_pop::Bool = false; correction = "none")
Calculate chi-squared test of HWE for each locus and returns observed and
expected heterozygosity with chi-squared, degrees of freedom and p-values for
each locus. Use `by_pop = true` to perform this separately for each population
(default: by_pop = false) and return a NamedTuple of DataFrames with the names
corresponding to the population names. Use `correction =` to specify a P-value
correction method for multiple testing.

#### example
`hwe_test(gulfsharks(), correction = "bh")` \n
`hwe_test(gulfsharks(), by_pop = true, correction = "bh")` \n

### `correction` methods (case insensitive)
- `"bonferroni"` : Bonferroni adjustment
- `"holm"` : Holm adjustment
- `"hochberg"` : Hochberg adjustment
- `"bh"` or `"b-h"` : Benjamini-Hochberg adjustment
- `"by"` or `"b-y"`: Benjamini-Yekutieli adjustment
- `"bl"` or `"b-l"` : Benjamini-Liu adjustment
- `"hommel"` : Hommel adjustment
- `"sidak"` : Šidák adjustment
- `"forward stop"` or `"fs"` : Forward-Stop adjustment
- `"bc"` or `"b-c"` : Barber-Candès adjustment
"""
function hwe_test(data::PopObj; by_pop::Bool = false, correction::String = "none")
    if !by_pop
        output = [
            Vector{Union{Missing,Float64}}(),
            Vector{Union{Missing,Int32}}(),
            Vector{Union{Missing,Float64}}()
            ]
        for locus in eachcol(data.loci, false)
            chisq, df, pval = locus_chi_sq(locus)
            push!.(output, [chisq, df, pval])
        end

        het = heterozygosity(data)
        insertcols!(het, size(het,2)+1, χ² = output[1])
        insertcols!(het, size(het,2)+1, DF = output[2])
        insertcols!(het, size(het,2)+1, P = output[3])

        # corrections
        @info "Χ² test for conformation to Hardy-Weinberg Equilibrium"
        if correction == "none"
            return het
        else
            corrected_P = multitest_missing(het.P, correction)

            # add the adjusted p-vals back to the dataframe and return
            insertcols!(het, 7, Pcorr = corrected_P )
        end
    else
        out_array = []
        # get population order and store in its own dict key
        popid = unique(data.samples.population |> collect)

        y = deepcopy(data.loci)
        insertcols!(y, 1, :population => data.samples.population)
        y_subdf = groupby(y[!, :], :population) |> collect
        het = heterozygosity(data, "pop")
        for (eachpop, het_pop) in zip(y_subdf, het)
            output = [[], [], []]
            for locus in eachcol(eachpop, false)[2:end]
                chisq, df, pval = locus_chi_sq(locus)
                push!.(output, [chisq, df, pval])
            end
            insertcols!(het_pop, size(het_pop,2)+1, χ² = output[1] |> Array{Union{Missing,Float64},1})
            insertcols!(het_pop, size(het_pop,2)+1, DF = output[2] |> Array{Union{Missing,Float64},1})
            insertcols!(het_pop, size(het_pop,2)+1, P = output[3] |> Array{Union{Missing,Float64},1})
            push!(out_array, het_pop)
        end
        @info "Χ² test for conformation to Hardy-Weinberg Equilibrium by Population"

        ## corrections
        if correction == "none"
            d = Dict([Symbol(i) => j  for (i,j) in zip(popid, out_array)])
            return NamedTuple{Tuple(Symbol.(popid))}(values(d))
        else
            p_array = []

            #for each in out_array
            for i in out_array
                append!(p_array, i.P)
            end

            # convert p_array to suitible Type
            p_vals = p_array |> Vector{Union{Missing,Float64}}

            # do actual correction
            corrected_P = multitest_missing(p_vals, correction)

            # split corrected pvals by #loci
            p_corrs = Iterators.partition(corrected_P, size(out_array[1],1)) |> collect

            # add each segment of split-pvals to the appropriate df
            for i in 1:length(out_array)
                insertcols!(out_array[i], size(out_array[i],2)+1, Pcorr = p_corrs[i])
            end
            d = Dict([Symbol(i) => j  for (i,j) in zip(popid, out_array)])
            return NamedTuple{Tuple(Symbol.(popid))}(values(d))
        end
    end
end

const hwe = hwe_test

"""
    locus_chi_sq(locus::Vector{<:Union{Missing, NTuple{N,<:Integer}}}) where N
Calculate the chi square statistic and p-value for a locus
Returns a tuple with chi-square statistic, degrees of freedom, and p-value
"""
function locus_chi_sq(locus::Vector{<:Union{Missing, NTuple{N,<:Integer}}}) where N
    # Give the function a locus from a genpop object and it will perform the ChiSquared test for HWE
    number_ind = count(i -> i !== missing, locus)

    # Get expected number of genotypes in a locus
    the_allele_dict = allele_freq(locus)

    # split the appropriate pairs into their own vectors
    alleles = Vector{String}()
    p = Vector{Float64}()
    for (i,j) in pairs(the_allele_dict)
        push!(alleles, "$i")
        push!(p, j)
    end

    #Calculate Expected Genotype numbers
    expected_genotype_freq = p * transpose(p) .* number_ind
    expected_genotype_freq = vec(expected_genotype_freq)

    alleles = (alleles .* ",") .* permutedims(alleles)
    alleles = Vector{String}.(sort.(split.(alleles, ",")))
    alleles = [parse.(Int16, i) |> Tuple for i in alleles]
    alleles = vec(alleles)

    expected = Dict()
    for (geno, freq) in zip(alleles, expected_genotype_freq)
        expected[geno] = get!(expected, geno, 0) + freq
    end

    ## Get observed number of genotypes in a locus
    observed = geno_freq(locus)
    for j in keys(observed)
        observed[j] = observed[j] * number_ind
    end

    chisq_stat = expected
    for genotype in keys(expected)
        o = get(observed, genotype, 0)
        e = get(expected, genotype, 0)

        chisq_stat[genotype] = (o - e)^2 / e
    end
    chisq_stat = values(chisq_stat) |> sum
    df = (length(alleles) - length(the_allele_dict)) / 2

    if df > 0
        chi_sq_dist = Distributions.Chisq(df)
        p_val = 1 - cdf(chi_sq_dist, chisq_stat)
    else
        p_val = missing
    end
    return (chisq_stat, df, p_val)
end

"""
    locus_chi_sq(locus::SubArray{<:Union{Missing, NTuple{N,<:Integer}},1}) where N
Calculate the chi square statistic and p-value for a locus of a `PopObj` split by
population using `groupby()`. Returns a tuple with chi-square statistic, degrees
of freedom, and p-value.
"""
function locus_chi_sq(locus::SubArray{<:Union{Missing, NTuple{N,<:Integer}},1}) where N
    #Give the function a locus from a genpop object and it will perform the ChiSquared test for HWE
    number_ind = count(i -> i !== missing, locus)

    # split the appropriate pairs into their own vectors
    alleles = Vector{String}()
    p = Vector{Float64}()
    for (i,j) in pairs(the_allele_dict)
        push!(alleles, "$i")
        push!(p, j)
    end

    # create a matrix of allele pair combinations
    alleles = (alleles .* ",") .* permutedims(alleles)
    alleles = Vector{String}.(sort.(split.(alleles, ",")))
    alleles = [parse.(Int16, i) |> Tuple for i in alleles]
    alleles = vec(alleles)

    #Calculate Expected Genotype numbers
    expected_genotype_freq = p * transpose(p) .* number_ind
    expected_genotype_freq = vec(expected_genotype_freq)

    expected = Dict()
    for (geno, freq) in zip(alleles, expected_genotype_freq)
        expected[geno] = get!(expected, geno, 0) + freq
    end

    ## Get observed number of genotypes in a locus
    observed = geno_freq(locus)
    for j in keys(observed)
        observed[j] = observed[j] * number_ind
    end

    chisq_stat = expected
    for genotype in keys(expected)
        o = get(observed, genotype, 0)
        e = get(expected, genotype, 0)

        chisq_stat[genotype] = (o - e)^2 / e
    end
    chisq_stat = values(chisq_stat) |> sum
    df = (length(alleles) - length(the_allele_dict)) / 2

    if df > 0
        chi_sq_dist = Distributions.Chisq(df)
        p_val = 1 - cdf(chi_sq_dist, chisq_stat)
    else
        p_val = missing
    end
    return (chisq_stat, df, p_val)
end

"""
    multitest_missing(pvals::Array{Float64,1}, correction::String)
Modification to `MultipleTesting.adjust` to include `missing` values in the
returned array. Missing values are first removed from the array, the appropriate
correction made, then missing values are re-added to the array at their original
positions. See MultipleTesting.jl docs for full more detailed information.
#### example
`multitest_missing([0.1, 0.01, 0.005, 0.3], "bh")`

### `correction` methods (case insensitive)
- `"bonferroni"` : Bonferroni adjustment
- `"holm"` : Holm adjustment
- `"hochberg"` : Hochberg adjustment
- `"bh"` or `"b-h"` : Benjamini-Hochberg adjustment
- `"by"` or `"b-y"`: Benjamini-Yekutieli adjustment
- `"bl"` or `"b-l"` : Benjamini-Liu adjustment
- `"hommel"` : Hommel adjustment
- `"sidak"` : Šidák adjustment
- `"forward stop"` or `"fs"` : Forward-Stop adjustment
- `"bc"` or `"b-c"` : Barber-Candès adjustment
"""
function multitest_missing(pvals::Array{<:Union{Missing, <:AbstractFloat},1}, correction::String)
    # make seperate array for non-missing P vals
    p_no_miss = skipmissing(pvals) |> collect
    # get indices of where original missing are
    miss_idx = findall(i -> i === missing, pvals)

    # make a dict of all possible tests and their respective functions
    d = Dict(
        "bonferroni" => Bonferroni(),
        "holm" => Holm(),
        "hochberg" => Hochberg(),
        "bh" => BenjaminiHochberg(),
        "b-h" => BenjaminiHochberg(),
        "by" => BenjaminiYekutieli(),
        "b-y" => BenjaminiYekutieli(),
        "bl" => BenjaminiLiu(),
        "b-l" => BenjaminiLiu(),
        "hommel" => Hommel(),
        "sidak" => Sidak(),
        "forward stop" => ForwardStop(),
        "fs" => ForwardStop(),
        "bc" => BarberCandes(),
        "b-c" => BarberCandes(),
    )

    correct = adjust(p_no_miss, d[lowercase(correction)]) |> Array{Union{Missing, Float64},1}

    # re-add missing to original positions
    for i in miss_idx
        insert!(correct, i, missing)
    end
    return correct
end

## JASON LOOK AT ME PLEEEEEASE

function het_obstest(data::PopObj)
    het_vals = Vector{Float64}()
    for locus in eachcol(data.loci, false)
        # get genotype freqs at locus
        a = geno_freq(locus)
        tmp = 0
        # total number of non-missing samples per locus
        n = skipmissing(locus) |> collect |> length
        genos = keys(a) |> collect
        for geno in genos
            if length(unique(geno)) != 1
                # if true, add freq to total in tmp
                tmp += a[geno]
            end
        end
        # include het correction
        het_adjust = tmp * (n/ (n-1))
        push!(het_vals, het_adjust)
    end

    return het_vals
end


text:
function random_tree(num_features::Int64, depth::Int64)
    if depth == 0
        return Leaf(rand())
    else
        return Inner(rand(1:num_features),
                     rand(),
                     random_tree(num_features, depth - 1),
                     random_tree(num_features, depth -1))
    end
end

random_feature_vector(num_features::Int64) = rand(num_features)
random_feature_vectors(num_examples::Int64, num_features::Int64) =
    map(_ -> random_feature_vector(num_features), 1:num_examples)

random_forest(num_trees::Int64, num_features::Int64, depth::Int64) =
    Forest([random_tree(num_features, depth) for _ in 1:num_trees])

construct_evaluators(f::Forest) =
    [:naive =>  f,
     :flattened => FlatForest(f),
     :compiled => CompiledForest(f)]

ext: jl


text:
n, p = 100, 5
((X, y, θ), (X1, y1, θ1)) = generate_continuous(n, p; seed=52)

@testset "linreg" begin
    lr  = LinearRegression(fit_intercept=false)
    lr1 = LinearRegression()
    β_ref = X \ y
    @test β_ref == fit(lr, X, y)

    # fit_intercept
    β_ref = X1 \ y1
    @test β_ref == fit(lr1, X, y1)

    # == iterative solvers
    β_cg = fit(lr1, X, y1; solver=CG())
    @test norm(β_cg - β_ref) / norm(β_ref) ≤ 1e-12
end

@testset "ridgereg" begin
    λ   = 1.0
    rr  = RidgeRegression(lambda=λ, fit_intercept=false)
    rr1 = RidgeRegression(λ, penalize_intercept=true)

    β_ref  = (X'X + λ*I) \ (X'y)
    β_ref1 = (X1'X1 + λ*I) \ (X1'y1)
    @test β_ref ≈ fit(rr, X, y)
    @test β_ref1 ≈ fit(rr1, X, y1)

    β_cg  = fit(rr, X, y; solver=CG())
    β_cg1 = fit(rr1, X, y1; solver=CG())

    @test norm(β_cg - β_ref) / norm(β_ref)   ≤ 1e-12
    @test norm(β_cg1 - β_ref1) / norm(β_ref) ≤ 1e-12
end

n, p = 500, 100
((X, y, θ), (X1, y1, θ1)) = generate_continuous(n, p;  seed=512, sparse=0.1)

@testset "lasso" begin
    # no intercept
    λ  = 50
    lr = LassoRegression(λ; fit_intercept=false)
    J  = objective(lr, X, y)
    θ_ls    = X \ y
    θ_ista  = fit(lr, X, y, solver=ISTA())
    θ_fista = fit(lr, X, y)

    @test isapprox(J(θ_ls),    586.25012, rtol=1e-5)
    @test isapprox(J(θ_fista), 539.11397, rtol=1e-5) # <- ref value
    @test isapprox(J(θ_ista),  539.11397, rtol=1e-5)
    # sparsity
    @test nnz(θ_ls)    == 100  # not sparse
    @test nnz(θ_fista) == 12 # sparse
    @test nnz(θ_ista)  == 12

    # with intercept
    lr1 = LassoRegression(λ, penalize_intercept=true)
    J   = objective(lr1, X, y1)
    θ_ls    = X1 \ y1
    θ_fista = fit(lr1, X, y1)
    θ_ista  = fit(lr1, X, y1, solver=ISTA())
    @test isapprox(J(θ_ls),    214.35485, rtol=1e-5)
    @test isapprox(J(θ_fista), 178.56434, rtol=1e-5)
    @test isapprox(J(θ_ista),  178.56433, rtol=1e-5)
    # sparsity
    @test nnz(θ_ls)    == 101 # not sparse
    @test nnz(θ_fista) == 5   # sparse
    @test nnz(θ_ista)  == 5

    # with intercept and not penalizing intercept
    lr1 = LassoRegression(λ)
    J   = objective(lr1, X, y1)
    θ_ls    = X1 \ y1
    θ_fista = fit(lr1, X, y1)
    θ_ista  = fit(lr1, X, y1, solver=ISTA())
    @test isapprox(J(θ_ls),    214.31393, rtol=1e-5)
    @test isapprox(J(θ_fista), 178.54691, rtol=1e-5)
    @test isapprox(J(θ_ista),  178.54690, rtol=1e-5)
    # sparsity
    @test nnz(θ_ls)    == 101 # not sparse
    @test nnz(θ_fista) == 6   # sparse
    @test nnz(θ_ista)  == 6

    if DO_COMPARISONS
        lr_sk = SKLEARN_LM.Lasso(alpha=λ/n, random_state=156123)
        lr_sk.fit(X, y1)
        θ1_sk = vcat(lr_sk.coef_[:], lr_sk.intercept_)
        @test isapprox(J(θ1_sk), 178.5, rtol=1e-3)
        @test nnz(θ1_sk) == 6
    end
end

@testset "elnet" begin
    ρ = 0.3
    α = 0.1
    λ = α * (1 - ρ) * n
    γ = α * ρ * n
    enr = ElasticNetRegression(λ, γ; penalize_intercept=true)
    J   = objective(enr, X, y1)
    θ_ls    = X1 \ y1
    θ_fista = fit(enr, X, y1)
    θ_ista  = fit(enr, X, y1, solver=ISTA())
    @test isapprox(J(θ_ls),    140.07105, rtol=1e-5)
    @test isapprox(J(θ_fista), 123.91270, rtol=1e-5)
    @test isapprox(J(θ_ista),  123.91268, rtol=1e-5)
    # sparsity
    @test nnz(θ_ls)    == 101 # not sparse
    @test nnz(θ_fista) == 7   # sparse
    @test nnz(θ_ista)  == 7

    # not penalizing intercept
    enr = ElasticNetRegression(λ, γ)
    J   = objective(enr, X, y1)
    θ_ls    = X1 \ y1
    θ_fista = fit(enr, X, y1)
    θ_ista  = fit(enr, X, y1, solver=ISTA())
    @test isapprox(J(θ_ls),    140.05876, rtol=1e-5)
    @test isapprox(J(θ_fista), 123.89384, rtol=1e-5)
    @test isapprox(J(θ_ista),  123.89383, rtol=1e-5)
    # sparsity
    @test nnz(θ_ls)    == 101 # not sparse
    @test nnz(θ_fista) == 8   # sparse
    @test nnz(θ_ista)  == 8

    if DO_COMPARISONS
        enr_sk = SKLEARN_LM.ElasticNet(alpha=α, l1_ratio=ρ, random_state=5515)
        enr_sk.fit(X, y1)
        θ_sk = vcat(enr_sk.coef_[:], enr_sk.intercept_)
        @test isapprox(J(θ_sk), 123.89, rtol=1e-3)
        @test nnz(θ_sk) == 8
    end
end


text:
# The file contains example from readme

# -----------------------------------------------------------------------------

import Pkg;
Pkg.add(url="https://github.com/sumiya11/RationalFunctionFields.jl")

# -----------------------------------------------------------------------------

using Nemo
using RationalFunctionFields

# -----------------------------------------------------------------------------

R, (x, y) = PolynomialRing(QQ, ["x", "y"])

generators = [
            (x^2 + y^2) // R(1),
            1 // (x^2)
]

# -----------------------------------------------------------------------------

rff = RationalFunctionField(generators)

print( contains(rff, y^2, p=0.999) )

print( contains(rff, y^2, p=1) )


text:
# # [Cantilever with a Tip Force](@id tipforce)
#
# This example shows how to predict the behavior of a cantilever beam that is subjected
# to a constant shear force at the tip.
#
# ![](../assets/tipforce-drawing.svg)
#
#-
#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`tipforce.ipynb`](@__NBVIEWER_ROOT_URL__/examples/tipforce.ipynb).
#-

using GXBeam, LinearAlgebra

L = 1
EI = 1e6

## shear force (applied at end)
λ = 0:0.5:16
p = EI/L^2
P = λ*p

## create points
nelem = 16
x = range(0, L, length=nelem+1)
y = zero(x)
z = zero(x)
points = [[x[i],y[i],z[i]] for i = 1:length(x)]

## index of endpoints of each beam element
start = 1:nelem
stop = 2:nelem+1

## compliance matrix for each beam element
compliance = fill(Diagonal([0, 0, 0, 0, 1/EI, 0]), nelem)

## create assembly of interconnected nonlinear beams
assembly = Assembly(points, start, stop, compliance=compliance)

## pre-initialize system storage
system = System(assembly, true)

## run an analysis for each prescribed tip load
states = Vector{AssemblyState{Float64}}(undef, length(P))
for i = 1:length(P)

    ## create dictionary of prescribed conditions
    prescribed_conditions = Dict(
        ## fixed left side
        1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0,
            theta_z=0),
        ## shear force on right tip
        nelem+1 => PrescribedConditions(Fz = P[i])
    )

    ## perform a static analysis
    static_analysis!(system, assembly;
        prescribed_conditions=prescribed_conditions)

    ## post-process the results
    states[i] = AssemblyState(system, assembly;
        prescribed_conditions=prescribed_conditions)

end

#!jl nothing #hide

#
# The analytical solution to this problem has been presented by several authors.  Here
# we follow the solution by H. J. Barten in "On the Deflection of a Cantilever Beam",
# after incorporating the corrections they submitted for finding the tip angle.
#

import Elliptic

δ = range(pi/4, pi/2, length=10^5)[2:end-1]

k = @. cos(pi/4)/sin(δ)
λ_a = @. (Elliptic.F(pi/2, k^2) - Elliptic.F(δ,  k^2))^2

θ_a = @. 2*(pi/4 - acos(k))

ξ_a = @. sqrt(2*sin(θ_a)/λ_a) .- 1

η_a = @. 1-2/sqrt(λ_a)*(Elliptic.E(pi/2, k^2) - Elliptic.E(δ, k^2))

#!jl nothing #hide

#
# Plotting the results reveals that the analytical and computational solutions show
# excellent agreement.
#

using Plots
#md using Suppressor #hide
pyplot()
#!jl nothing #hide

#-

#md @suppress_err begin #hide

u = [states[i].points[end].u[1] for i = 1:length(P)]
θ = [states[i].points[end].theta[2] for i = 1:length(P)]
w = [states[i].points[end].u[3] for i = 1:length(P)]

## set up the plot
plot(
    xlim = (0, 16),
    xticks = 0:1:16,
    xlabel = "Nondimensional Force \$\\left(\\frac{PL^2}{EI}\\right)\$",
    ylim = (0, 1.2),
    yticks = 0.0:0.2:1.2,
    ylabel = "Nondimensional Tip Displacements",
    grid = false,
    overwrite_figure=false
    )

plot!([0], [0], color=:black, label="Analytical")
scatter!([0], [0], color=:black, label="GXBeam")
plot!([0], [0], color=1, label="Vertical \$\\left(w/L\\right)\$")
plot!([0], [0], color=2, label="Horizontal \$\\left(-u/L\\right)\$")
plot!([0], [0], color=3, label="\$ \\theta/(\\pi/2) \$")

plot!(λ_a, η_a, color=1, label="")
scatter!(λ, w/L, color=1, label="")

plot!(λ_a, -ξ_a, color=2, label="")
scatter!(λ, -u/L, color=2, label="")

plot!(λ_a, θ_a*2/pi, color=3, label="")
scatter!(λ, -4*atan.(θ/4)*2/pi, color=3, label="")

#!nb plot!(show=true)
#md savefig("../assets/tipforce-displacement.svg") #hide
#md closeall() #hide
#md end #hide
#md nothing #hide

#md # ![](../assets/tipforce-displacement.svg)


text:
#include("./src/SeisBeamforming.jl")
include("/Users/kurama/Documents/kurama/research/SCECInternship2019/Beamforming/SeisBeamforming/src/SeisBeamforming.jl")
using .SeisBeamforming
using SeisIO, Noise, JLD2, FileIO, PlotlyJS, Geodesy, Statistics, LinearAlgebra, Dates

#==================================================#
# Input Parameters

channel = ["BP1"] # channel name to do beam forming; ["BP1", "HHZ" is available]

bp_len = 60 * 60 # time window to conduct beam forming
bp_step = 60 * 30 # overlap of time window
freq_min = 0.1 # low freq at band pass filter before FFT
freq_max = 5.0 # high freq at band pass filter before FFT
sampling_freq = 20.0 # sampling frequency used to downsampling at FFT
# examined parameters
freq_range = [0.2, 3.0] # frequency range for evaluating beam power: used to choose omega from FFTfreq [1/s]
fexamin_num = 10# num of examined frequency

vel_range  = [1000, 5000]      # evaluating velocity range [cmin cmax] [m/s]
Δvel = 200 # grid space of velocity [m/s]

Δθ   = 10 # grid space of azimuth of trial slowness vector [deg]

finame      = "./dataset/BPnetwork.jld2"
fodir   = "./dataset"
foname  = "BP_beamform" # data is saved at ./dataset/$foname.jld2

IsPlotFigure = false # plot figures to debug
#==================================================#

mkpath(fodir)
fopath=joinpath(fodir, foname*".jld2")

InputDictionary = Dict([
        "channel"  	=> channel,
	"bp_len"	=> bp_len,
	"bp_step"	=> bp_step,
	"freq_min"	=> freq_min,
	"freq_max"	=> freq_max,
	"sampling_freq"	=> sampling_freq,
        "freq_range"	=> freq_range,
	"fexamin_num"	=> fexamin_num,
        "vel_range" 	=> vel_range,
        "Δvel"   	=> Δvel,
        "Δθ"     	=> Δθ,
        "finame"	=> finame,
        "fopath"        => fopath,
	"IsPlotFigure"  => IsPlotFigure
    ])

# mass request with input Dictionary
seisnoisebeamform(InputDictionary)

max_stars_repo_name: kura-okubo/SeisBeamforming

max_issues_repo_name: kura-okubo/SeisBeamforming

max_forks_repo_name: kura-okubo/SeisBeamforming


text:
"""
    function subs_(exp,ob,tar)

Substitutes all encounters of a symbolic expression or symbol **ob** by a target expression **tar** in **exp**.
Is update true, only substitutes those symbols **ob** that are being assigned.

Example
```@Julia
> subs_(:(x += 3*x),:x,:y,update=false)
:(y+=3*y)
> subs_(:(x += 3*x),:x,:y,update=true)
:(y+=3*x)
```
"""
function subs_(exp::Expr,ob::Union{Expr,Symbol},tar::Union{Expr,Symbol,<:Number};update=false)
    for (pos,a) in enumerate(exp.args)
        if a == ob && !update
            exp.args[pos] = tar
        elseif a == ob && exp.head in [:(=),:(+=),:(-=),:(%=),:(/=)] && pos == 1
            exp.args[pos] = tar
        elseif typeof(a) == Expr
            a = subs_(a,ob,tar,update=update)
        end
    end
    return exp
end

"""
    function subsArguments_(exp,ob,tar)

Substitutes all encounters of a symbol **ob** by an expanded list of symbols **tar** in **exp**.
Mainly used to put a list of arguments into a function argument list.
"""
function subsArguments_(exp::Expr,ob::Symbol,tar::Array{Symbol,1})

    if typeof(tar) == Symbol
        return subs_(exp,ob,tar)
    elseif typeof(tar) <: Array
        add = ""
        for i in tar
            add = string(add,i,",")
        end

        return Meta.parse(replace(string(exp),string(ob)=>add[1:end-1]))
    else
        error("Arguments sent are wrong.")
    end
end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
module _Pi

include("pi_chudnovsky.jl")
include("pi_machin.jl")
include("pi_iterative.jl")

end


text:
# Zadatak 8.1
include("labyrinth.jl")
include("zadatak81DFS.jl")

# definisanje lavirinta
nbRows = 8
nbColumns = 8
G = emptyLabyrinth(nbRows, nbColumns)
# definisanje vrsta cvorova
defineNodesType!(G, 1:nbRows*nbColumns, 'P'); #PUT
defineNodesType!(G, [5 9 10 11 13 15 21 23 24 26 28 29 31 34 44 45 46 47 50 55 58 60 61 63], 'Z'); #ZID

# poziv za trazenje puta od dvora 57 do cvora 1
path = DFS!(G, 57, 62);

println("Ispis putanje i izgled lavirinta posle obilaska")
println("Path = $path")
printLabyrinthPath(G, path)


text:
module MaxwellBase

# @reexport makes all exported symbols of the exported packages available in module using
# MaxwellBase.  Basically, it is equivalent to calling `using ...` prepended by @reexport
# automatically upon calling `using MaxwellBase`.
using Reexport
@reexport using LinearAlgebra, SparseArrays, StaggeredGridCalculus, GeometryPrimitives
using AbbreviatedTypes
using SimpleConstants
using Base.Threads: @threads  # used in param.jl

export S²ComplexF1, S²ComplexF2, S²ComplexF3, MatParam, ParamInd, ObjInd

## Type aliases
const S²ComplexF1 = S²ComplexF{1,1}
const S²ComplexF2 = S²ComplexF{2,4}
const S²ComplexF3 = S²ComplexF{3,9}

const MatParam = Union{Number,AbsVecNumber,AbsMatNumber}
const ParamInd = UInt8  # change this to handle more than 2⁸ = 256 materials
const ObjInd = UInt16  # change this to handle more than 2¹⁶ = 65536 objects

# The order of inclusion matters: if types or functions in file A are used in file B, file A
# must be included first.
include("enumtype.jl")
include("util.jl")
include("unit.jl")
include("material.jl")
include("object.jl")
include("field.jl")
include("assignment.jl")
include("smoothing.jl")
include("param.jl")

end # module MaxwellBase


text:
using PackageTest
using Test

@testset "PackageTest.jl" begin
    # 2x + y
    @test my_f(2,1) == 5
end


text:
using Sandbox
using Base.Test

# write your own tests here
@test 1 == 1

@test @stub :read


text:
require("DataFrames")
using DataFrames

# The current DataFrame type hierarchy relatively transparently implements
# the notion of a set of named, equal-length, heterogeneous-typed vectors that may include
# missing data.

# An AbstractDataFrame is something that can be indexed in a [rows,cols] syntax.
# It currently inherits Associative functionality, which may be too broad --
# additional experience will be necessary.

#abstract AbstractDataFrame <: Associative{Any,Any}

# The standard implementation of an ADF is a DataFrame. It has a set of columns
# and a set of column names that refer to the columns. (Note that the column
# names are stored in an Index data structure that allows fast lookups and
# some clever features described below.) Column names are of type ByteString --
# either ASCIIString or UTF8String. The current implementation allows columns
# to be of any type, which may be too broad, as described below.

#type DataFrame <: AbstractDataFrame
#    columns::Vector{Any}
#    colindex::Index
# ...

# a SubDataFrame is a lightweight wrapper around a DataFrame used most frequently in
# split/apply sorts of operations.
#type SubDataFrame <: AbstractDataFrame
#    parent::DataFrame
#    rows::Vector{Int} # maps from subdf row indexes to parent row indexes
# ...

# The typical column type is a DataVector, which inherits from AbstractDataVector. A
# DataVector is a vector of some type, along with a Boolean vector that acts as
# a mask for missing data. The decision to use a mask instead of a particular
# value (e.g., NaN + a particular payload for floats) remains contentious,
# as it is in all other statistical programming languages and frameworks.
#type DataVector{T} <: AbstractDataVector{T}
#    data::Vector{T}
#    na::AbstractVector{Bool}

# DataVectors can be built and used independently of DataFrames. For interactive
# use, the type-referencing constructor is handy, and the NA literal is supported.
dvint = DataVector[1, 2, NA, 4]
# DataVectors can also be contructed from existing arrays. Note that the existing
# array is referenced, not copied. In general, the DataFrames.jl philosophy is
# similar to the overall Julia philosophy -- don't copy data unless requested by
# the user or otherwise necessary.
x = [5:8]
dvint2 = DataArray(x)
dvint2[1] = 10
print(x)
# Note that assignment of NA to a DataVector does not affect the underlying
# vector, just the mask.
dvint2[2] = NA
print(dvint2)
print(x)

# an aspect of DataVectors that needs to be revisited and likely redesigned is that
# they contain meta-data specifying how the NAs should be treated under
# certain operations. To create a new DataVector object (with the same data
# and mask objects) that looks as if the NAs are ignored, use naFilter; to
# define a replacement value, use naReplace.
try sum(dvint) end # fails now -- eventually, this should work.
sum(naFilter(dvint))
sum(naReplace(dvint, 100))

# A PooledDataVector is an implementation of AbstractDataVector that uses indexes into
# a pool of values instead of an array of values, which is space efficient in
# the case of a large number of repeated elements. A common use is repeated
# strings. PDVs don't have a mask -- they can use a 0-index to represent the NA.
pdvstr = PooledDataVector["one", "one", "two", "two", NA, "one", "one"]
idump(pdvstr)  # show the internal structure of a PooledDataVector
# some operations, such as csvDataFrame which loads a CSV file into a DataFrame
# will try to created a PooledDataVector when possible.

# A DataFrame is a mix of DataVectors, PooledDataVectors, and (maybe) other vectors

# this useful constructor from an expression unpacks the assignments, converts the LHSs
# to column names, and converts the RHS to DataVectors.
df1 = DataFrame(quote
    a = shuffle([1:10])
    b = ["A","B"][rand(1:2,10)]
    v1 = randn(10)
end)
typeof(df1["a"])

# Let's look at the structure of the data (`dump` is like R's `str`).
dump(df1)

# The names of the columns in a DF are stored in a new type called an Index (possibly
# needs a less-generic name?)
# an AbstractIndex is a thing that can be used to look up ordered things by name, but that
# will also accept a position or set of positions or range or other things and pass them
# through cleanly.
# This makes referencing columns by name or position very simple and flexible.
ii = Index(["cat", "dog", "bird", "turtle"])
ii["cat"]
ii[["turtle", "dog"]]
ii[2]
ii[2:3]
# Even better, you can define meta-names
set_group(ii, "mammals", ["cat", "dog"])
ii["mammals"]
# Here is the structure of an index, a Dict pointing to a list of names:
dump(ii)
# and this can be very useful shorthand for DFs
set_group(df1, "IVs", ["a", "v1"])
print(df1)
df1["IVs"]
# note that this expansion should work in Formulas

# rbind and cbind should eventually work like their R equivalents; hcat and vcat are
# Julian aliases. They're methods need fleshing out at the moment
cbind(df1, DataFrame(:("cat"=[5:14])))
rbind(df1, DataFrame(:("a"=99; "b"="B"; "v1"=3.1415)))

# a SubDataFrame is a lightweight wrapper around a DataFrame used most frequently in
# split/apply sorts of operations.
df1 = DataFrame(quote
    a = shuffle([1:10])
    b = ["A","B"][rand(1:2,10)]
    v1 = randn(10)
end)
sub(df1, [1:5])
sub(df1, [1:5], "v1")

# Here is the structure of a SubDataFrame:
idump( sub(df1, [1:5]), 2)

# Expressions are used in many places in DataFrames.jl. The functions
# `with`, `within`, and `based_on` evaluate an expression within a
# DataFrame. `with` returns the result of the evaluation. These all
# work by traversing the expression and replacing terms that are in
# the DataFrame and then eval'ing the expression.
with(df1, :( a + sum(v1) ))
# Because `eval` acts globally, it can pull in global variables but
# not local variables. To insert those, they must be escaped.
a2local = df1["a"] + df1["a"]
with(df1, :( a + $a2local))

# `within(df)` returns the a copy of `df` after evaluating the
# expression. This is most useful for creating new columns from
# existing columns. `within!` is the same but changes the original
# DataFrame.
within(df1, :( a2 = 2 + a ))

# `based_on(df)` evaluates to a new DataFrame using the contents of `df`.
df2 = based_on(df1, quote
    a = a + 2
    v1sum = sum(v1)
end)

# Rows of DataFrames may also be selected based on the results of an
# expression, somewhat like R's data.table package.
df1[:( a .> 5 )]
# equivalent to:
df1[with(df1, :(a .> 5)), :]

# TODO: groupedDataFrame

# TODO: grouping operations

# TODO: pipelining
# The pipeline operator (`|`) is useful for many DataFrames.jl operations.
# This is handy at the REPL to evaluate the structure of a result:
df1[:( a .> 5 )] | dump

# TODO: by, colwise

# TODO: stack, unstack


text:
struct SymTrTenField{T,N,N2,L} <: AbstractSymTrTenArray{Complex{T},N}
    rr::SymTrTenArray{T,N,Array{T,N},Array{T,N},Array{T,N},Array{T,N},Array{T,N}}
    c::SymTrTenArray{Complex{T},N,ScalarField{T,N,N2,L},ScalarField{T,N,N2,L},ScalarField{T,N,N2,L},ScalarField{T,N,N2,L},ScalarField{T,N,N2,L}}
    r::SubArray{SymTen{T},N,SymTrTenArray{T,N,Array{T,N},Array{T,N},Array{T,N},Array{T,N},Array{T,N}},Tuple{Base.OneTo{Int},Vararg{Base.Slice{Base.OneTo{Int}},N2}},L}

    function SymTrTenField{T,N,N2,L}(xx::ScalarField{T,N,N2,L},xy::ScalarField{T,N,N2,L},xz::ScalarField{T,N,N2,L},yy::ScalarField{T,N,N2,L},yz::ScalarField{T,N,N2,L}) where {T,N,N2,L}
        c = SymTrTenArray(xx,xy,xz,yy,yz)
        rr = SymTrTenArray(InplaceRealFFT.data(xx),InplaceRealFFT.data(xy),InplaceRealFFT.data(xz),InplaceRealFFT.data(yy),InplaceRealFFT.data(yz))
        r = view(rr, Base.OneTo(size(xx.field.r, 1)), ntuple(i->Colon(), Val(N2))...)
        return new{T,N,N2,L}(rr,c,r)
    end

end

@inline function Base.getproperty(a::S,s::Symbol) where {S<:SymTrTenField}
    if (s === :kx || s === :ky || s === :kz || s === :k)
        return getproperty(getfield(getfield(a,:c),:xx),s)
    else
        return getfield(a,s)
    end
end

@inline FluidTensors.xxvec(v::SymTrTenField) =
    FluidTensors.xxvec(v.c)
@inline FluidTensors.xyvec(v::SymTrTenField) =
    FluidTensors.xyvec(v.c)
@inline FluidTensors.xzvec(v::SymTrTenField) =
    FluidTensors.xzvec(v.c)
@inline FluidTensors.yyvec(v::SymTrTenField) =
    FluidTensors.yyvec(v.c)
@inline FluidTensors.yzvec(v::SymTrTenField) =
    FluidTensors.yzvec(v.c)

SymTrTenField(xx::ScalarField{T,N,N2,L},xy::ScalarField{T,N,N2,L},xz::ScalarField{T,N,N2,L},yy::ScalarField{T,N,N2,L},yz::ScalarField{T,N,N2,L}) where {T,N,N2,L} = SymTrTenField{T,N,N2,L}(xx,xy,xz,yy,yz)

SymTrTenField{T}(dims::NTuple{3,Int},l::NTuple{3,Real}) where T = SymTrTenField(ScalarField{T}(dims,l),ScalarField{T}(dims,l),ScalarField{T}(dims,l),ScalarField{T}(dims,l),ScalarField{T}(dims,l))
SymTrTenField(dims::NTuple{3,Int},l::NTuple{3,Real}) = SymTrTenField(ScalarField(dims,l),ScalarField(dims,l),ScalarField(dims,l),ScalarField(dims,l),ScalarField(dims,l))

Base.similar(a::SymTrTenField) = SymTrTenField(similar(a.c.xx),similar(a.c.xy),similar(a.c.xz),similar(a.c.yy),similar(a.c.yz))

function InplaceRealFFT.rfft!(v::SymTrTenField)
    rfft!(v.c.xx)
    rfft!(v.c.xy)
    rfft!(v.c.xz)
    rfft!(v.c.yy)
    rfft!(v.c.yz)
    return v
end

function InplaceRealFFT.brfft!(v::SymTrTenField)
    brfft!(v.c.xx)
    brfft!(v.c.xy)
    brfft!(v.c.xz)
    brfft!(v.c.yy)
    brfft!(v.c.yz)
    return v.r
end

function InplaceRealFFT.irfft!(v::SymTrTenField)
    irfft!(v.c.xx)
    irfft!(v.c.xy)
    irfft!(v.c.xz)
    irfft!(v.c.yy)
    irfft!(v.c.yz)
    return v.r
end

isrealspace(a::SymTrTenField) =
    isrealspace(a.c.xx) && isrealspace(a.c.xy) && isrealspace(a.c.xz) && isrealspace(a.c.yy) && isrealspace(a.c.yz)

Base.show(io::IO,m::MIME"text/plain",a::SymTrTenField) =
    isrealspace(a) ? show(io,m,a.r) : show(io,m,a.c)

setreal!(a::SymTrTenField) = (setreal!(a.c.xx); setreal!(a.c.xy); setreal!(a.c.xz); setreal!(a.c.yy); setreal!(a.c.yz))

setfourier!(a::SymTrTenField) = (setfourier!(a.c.xx); setfourier!(a.c.xy); setfourier!(a.c.xz); setfourier!(a.c.yy); setfourier!(a.c.yz))

fourier!(a::SymTrTenField) = (fourier!(a.c.xx); fourier!(a.c.xy); fourier!(a.c.xz); fourier!(a.c.yy); fourier!(a.c.yz); a)

real!(a::SymTrTenField) = (real!(a.c.xx); real!(a.c.xy); real!(a.c.xz); real!(a.c.yy); real!(a.c.yz); a.r)

################################ SymTenField ########################################

struct SymTenField{T,N,N2,L} <: AbstractSymTenArray{Complex{T},N}
    rr::SymTenArray{T,N,Array{T,N},Array{T,N},Array{T,N},Array{T,N},Array{T,N},Array{T,N}}
    c::SymTenArray{Complex{T},N,ScalarField{T,N,N2,L},ScalarField{T,N,N2,L},ScalarField{T,N,N2,L},ScalarField{T,N,N2,L},ScalarField{T,N,N2,L},ScalarField{T,N,N2,L}}
    r::SubArray{SymTen{T},N,SymTenArray{T,N,Array{T,N},Array{T,N},Array{T,N},Array{T,N},Array{T,N},Array{T,N}},Tuple{Base.OneTo{Int},Vararg{Base.Slice{Base.OneTo{Int}},N2}},L}

    function SymTenField{T,N,N2,L}(xx::ScalarField{T,N,N2,L},xy::ScalarField{T,N,N2,L},xz::ScalarField{T,N,N2,L},yy::ScalarField{T,N,N2,L},yz::ScalarField{T,N,N2,L},zz::ScalarField{T,N,N2,L}) where {T,N,N2,L}
        c = SymTenArray(xx,xy,xz,yy,yz,zz)
        rr = SymTenArray(InplaceRealFFT.data(xx),InplaceRealFFT.data(xy),InplaceRealFFT.data(xz),InplaceRealFFT.data(yy),InplaceRealFFT.data(yz),InplaceRealFFT.data(zz))
        r = view(rr, Base.OneTo(size(xx.field.r, 1)), ntuple(i->Colon(), Val{N2}())...)
        return new{T,N,N2,L}(rr,c,r)
    end

end

@inline function Base.getproperty(a::S,s::Symbol) where {S<:SymTenField}
    if (s === :kx || s === :ky || s === :kz || s === :k)
        return getproperty(getfield(getfield(a,:c),:xx),s)
    else
        return getfield(a,s)
    end
end

@inline FluidTensors.xxvec(v::SymTenField) =
    FluidTensors.xxvec(v.c)
@inline FluidTensors.xyvec(v::SymTenField) =
    FluidTensors.xyvec(v.c)
@inline FluidTensors.xzvec(v::SymTenField) =
    FluidTensors.xzvec(v.c)
@inline FluidTensors.yyvec(v::SymTenField) =
    FluidTensors.yyvec(v.c)
@inline FluidTensors.yzvec(v::SymTenField) =
    FluidTensors.yzvec(v.c)
@inline FluidTensors.zzvec(v::SymTenField) =
    FluidTensors.zzvec(v.c)

SymTenField(xx::ScalarField{T,N,N2,L},xy::ScalarField{T,N,N2,L},xz::ScalarField{T,N,N2,L},yy::ScalarField{T,N,N2,L},yz::ScalarField{T,N,N2,L},zz::ScalarField{T,N,N2,L}) where {T,N,N2,L} = SymTenField{T,N,N2,L}(xx,xy,xz,yy,yz,zz)

SymTenField{T}(dims::NTuple{3,Int},l::NTuple{3,Real}) where T = SymTenField(ScalarField{T}(dims,l),ScalarField{T}(dims,l),ScalarField{T}(dims,l),ScalarField{T}(dims,l),ScalarField{T}(dims,l),ScalarField{T}(dims,l))
SymTenField(dims::NTuple{3,Int},l::NTuple{3,Real}) = SymTenField(ScalarField(dims,l),ScalarField(dims,l),ScalarField(dims,l),ScalarField(dims,l),ScalarField(dims,l),ScalarField(dims,l))

Base.similar(a::SymTenField) = SymTenField(similar(a.c.xx),similar(a.c.xy),similar(a.c.xz),similar(a.c.yy),similar(a.c.yz),similar(a.c.zz))

function InplaceRealFFT.rfft!(v::SymTenField)
    rfft!(v.c.xx)
    rfft!(v.c.xy)
    rfft!(v.c.xz)
    rfft!(v.c.yy)
    rfft!(v.c.yz)
    rfft!(v.c.zz)
    return v
end

function InplaceRealFFT.brfft!(v::SymTenField)
    brfft!(v.c.xx)
    brfft!(v.c.xy)
    brfft!(v.c.xz)
    brfft!(v.c.yy)
    brfft!(v.c.yz)
    brfft!(v.c.zz)
    return v.r
end

function InplaceRealFFT.irfft!(v::SymTenField)
    irfft!(v.c.xx)
    irfft!(v.c.xy)
    irfft!(v.c.xz)
    irfft!(v.c.yy)
    irfft!(v.c.yz)
    irfft!(v.c.zz)
    return v.r
end

isrealspace(a::SymTenField) =
    isrealspace(a.c.xx) && isrealspace(a.c.xy) && isrealspace(a.c.xz) && isrealspace(a.c.yy) && isrealspace(a.c.yz) && isrealspace(a.c.zz)

Base.show(io::IO,m::MIME"text/plain",a::SymTenField) =
    isrealspace(a) ? show(io,m,a.r) : show(io,m,a.c)

setreal!(a::SymTenField) = (setreal!(a.c.xx); setreal!(a.c.xy); setreal!(a.c.xz); setreal!(a.c.yy); setreal!(a.c.yz); setreal!(a.c.zz))

setfourier!(a::SymTenField) = (setfourier!(a.c.xx); setfourier!(a.c.xy); setfourier!(a.c.xz); setfourier!(a.c.yy); setfourier!(a.c.yz); setfourier!(a.c.zz))

fourier!(a::SymTenField) = (fourier!(a.c.xx); fourier!(a.c.xy); fourier!(a.c.xz); fourier!(a.c.yy); fourier!(a.c.yz); fourier!(a.c.zz); a)

real!(a::SymTenField) = (real!(a.c.xx); real!(a.c.xy); real!(a.c.xz); real!(a.c.yy); real!(a.c.yz); real!(a.c.zz); a.r)

################################################################# AntiSymTenField
struct AntiSymTenField{T,N,N2,L} <: AbstractAntiSymTenArray{Complex{T},N}
    rr::AntiSymTenArray{T,N,Array{T,N},Array{T,N},Array{T,N}}
    c::AntiSymTenArray{Complex{T},N,ScalarField{T,N,N2,L},ScalarField{T,N,N2,L},ScalarField{T,N,N2,L}}
    r::SubArray{AntiSymTen{T},N,AntiSymTenArray{T,N,Array{T,N},Array{T,N},Array{T,N}},Tuple{Base.OneTo{Int},Vararg{Base.Slice{Base.OneTo{Int}},N2}},L}

    function AntiSymTenField{T,N,N2,L}(xy::ScalarField{T,N,N2,L},xz::ScalarField{T,N,N2,L},yz::ScalarField{T,N,N2,L}) where {T,N,N2,L}
        c = AntiSymTenArray(xy,xz,yz)
        rr = AntiSymTenArray(InplaceRealFFT.data(xy),InplaceRealFFT.data(xz),InplaceRealFFT.data(yz))
        r = view(rr, Base.OneTo(size(xy.field.r, 1)), ntuple(i->Colon(), Val(N2))...)
        return new{T,N,N2,L}(rr,c,r)
    end

end

@inline function Base.getproperty(a::S,s::Symbol) where {S<:AntiSymTenField}
    if (s === :kx || s === :ky || s === :kz || s === :k)
        return getproperty(getfield(getfield(a,:c),:xy),s)
    else
        return getfield(a,s)
    end
end

@inline FluidTensors.xyvec(v::AntiSymTenField) =
    FluidTensors.xyvec(v.c)
@inline FluidTensors.xzvec(v::AntiSymTenField) =
    FluidTensors.xzvec(v.c)
@inline FluidTensors.yzvec(v::AntiSymTenField) =
    FluidTensors.yzvec(v.c)

AntiSymTenField(xy::ScalarField{T,N,N2,L},xz::ScalarField{T,N,N2,L},yz::ScalarField{T,N,N2,L}) where {T,N,N2,L} = AntiSymTenField{T,N,N2,L}(xy,xz,yz)

AntiSymTenField{T}(dims::NTuple{3,Int},l::NTuple{3,Real}) where T = AntiSymTenField(ScalarField{T}(dims,l),ScalarField{T}(dims,l),ScalarField{T}(dims,l))
AntiSymTenField(dims::NTuple{3,Int},l::NTuple{3,Real}) = AntiSymTenField(ScalarField(dims,l),ScalarField(dims,l),ScalarField(dims,l))

Base.similar(a::AntiSymTenField) = AntiSymTenField(similar(a.c.xy),similar(a.c.xz),similar(a.c.yz))

function InplaceRealFFT.rfft!(v::AntiSymTenField)
    rfft!(v.c.xy)
    rfft!(v.c.xz)
    rfft!(v.c.yz)
    return v
end

function InplaceRealFFT.brfft!(v::AntiSymTenField)
    brfft!(v.c.xy)
    brfft!(v.c.xz)
    brfft!(v.c.yz)
    return v.r
end

function InplaceRealFFT.irfft!(v::AntiSymTenField)
    irfft!(v.c.xy)
    irfft!(v.c.xz)
    irfft!(v.c.yz)
    return v.r
end

isrealspace(a::AntiSymTenField) =
    isrealspace(a.c.xy) && isrealspace(a.c.xz) && isrealspace(a.c.yz)

Base.show(io::IO,m::MIME"text/plain",a::AntiSymTenField) =
    isrealspace(a) ? show(io,m,a.r) : show(io,m,a.c)

setreal!(a::AntiSymTenField) = (setreal!(a.c.xy); setreal!(a.c.xz); setreal!(a.c.yz))

setfourier!(a::AntiSymTenField) = (setfourier!(a.c.xy); setfourier!(a.c.xz); setfourier!(a.c.yz))

fourier!(a::AntiSymTenField) = (fourier!(a.c.xy); fourier!(a.c.xz); fourier!(a.c.yz); a)

real!(a::AntiSymTenField) = (real!(a.c.xy); real!(a.c.xz); real!(a.c.yz); a.r)

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

const AbstractTensorField{T} = Union{<:SymTenField{T},<:SymTrTenField{T},<:AntiSymTenField{T}}


text:
using Knet,ArgParse,Compat,GZip
using MNIST

function drelu(x)
    return x>0?1:0
end

function predict!(x,w,ao,ho)

    j=1
    for i=1:2:length(w)
        x = w[i]*x .+ w[i+1]
        ao[j].=x
        if i<length(w)-1
            x = relu(x) # max(0,x)
            ho[j].=x
        end
        j+=1
        println(j)
    end

    return x
end

function ffward!(x,ygold,w,ao,ho)
    ypred=predict!(x,w,ao,ho)
    ynorm = logp(ypred,1)
    err = -sum(ygold .* ynorm) / size(ygold,2)
    #err = ypred - ygold
    return err
end

function updateParams!(x,w,ho,δao,err;lr=0.5)
    w[1] -= ((δao[1] * x') .* lr)
    j=2
    for i=3:2:(length(w)-2)
        w[i] -= ((δao[j] * ho[j-1]') .* lr)
        j+=1
    end
    #w[length(w)-1] -= ((fill!(zeros(10,100),err) * ho[j-1]') .* lr)
end

function dfa_train(w,ao,ho,B, dtrn; lr=.5, epochs=10)

    δao = Any[]
    for epoch=1:epochs
        for (x,y) in dtrn
            err = ffward!(x, y,w,ao,ho)
            for i=1:(length(ao)-1)
                push!(δao,(B[i].*err) .* drelu.(ao[i]))
            end
            updateParams!(x,w,ho,δao,err;lr=lr)
            println(sum(w[3]))
        end
    end

    return w
end

function accuracy(w,ao,ho, dtst)

    ncorrect = ninstance = nloss = 0

    for (x, ygold) in dtst

        ypred = predict!(x,w,ao,ho)

        ynorm = ypred .- log(sum(exp(ypred),1))
        nloss += -sum(ygold .* ynorm)
        ncorrect += sum(ygold .* (ypred .== maximum(ypred,1)))
        ninstance += size(ygold,2)
    end
    return (ncorrect/ninstance, nloss/ninstance)
end

function initParams(h,batchsize; atype=Array{Float32}, winit=0.1)
    w = Any[]
    ao = Any[]
    ho = Any[]
    B = Any[]
    x = 28*28
    for y in [h..., 10]
        push!(w, convert(atype, winit*randn(y,x)))
        push!(w, convert(atype, randn(y, 1)))
        push!(ao,convert(atype,zeros(y,batchsize)))
        push!(ho,convert(atype,zeros(y,batchsize)))
        push!(B,convert(atype,rand(y,x)))
        x = y
    end
    return (w,ao,ho,B)
end

function loaddata()
    global xtrn,ytrn,xtst,ytst
    info("Loading MNIST...")
    xtrn = gzload("train-images-idx3-ubyte.gz")[17:end]
    xtst = gzload("t10k-images-idx3-ubyte.gz")[17:end]
    ytrn = gzload("train-labels-idx1-ubyte.gz")[9:end]
    ytst = gzload("t10k-labels-idx1-ubyte.gz")[9:end]
end

function gzload(file; path=Knet.dir("data",file), url="http://yann.lecun.com/exdb/mnist/$file")
    isfile(path) || download(url, path)
    f = gzopen(path)
    a = @compat read(f)
    close(f)
    return(a)
end

function minibatch(x, y, batchsize; atype=Array{Float32}, xrows=784, yrows=10, xscale=255)
    xbatch(a)=convert(atype, reshape(a./xscale, xrows, div(length(a),xrows)))
    ybatch(a)=(a[a.==0]=10; convert(atype, sparse(convert(Vector{Int},a),1:length(a),one(eltype(a)),yrows,length(a))))
    xcols = div(length(x),xrows)
    xcols == length(y) || throw(DimensionMismatch())
    data = Any[]
    for i=1:batchsize:xcols-batchsize+1
        j=i+batchsize-1
        push!(data, (xbatch(x[1+(i-1)*xrows:j*xrows]), ybatch(y[i:j])))
    end
    return data
end

#function main()

    o=Dict{Symbol,Any}()
    o[:seed]=-1; o[:batchsize]=100; o[:hidden]=[100;50;30;20];
    o[:epochs]=10; o[:lr]=0.5; o[:atype]="Array{Float32}";
    o[:gcheck]=0; o[:winit]=0.1; o[:fast]=false;

    if !o[:fast]
        println("opts=",[(k,v) for (k,v) in o]...)
    end

    atype = eval(parse(o[:atype]))

    (w,ao,ho,B) = initParams(o[:hidden],o[:batchsize]; atype=atype, winit=o[:winit])

    if !isdefined(MNIST,:xtrn); loaddata(); end

    global dtrn = minibatch(xtrn, ytrn, o[:batchsize]; atype=atype)
    global dtst = minibatch(xtst, ytst, o[:batchsize]; atype=atype)

    report(epoch)=println((:epoch,epoch,:trn,accuracy(w,ao,ho,dtrn),:tst,accuracy(w,ao,ho,dtst)))

    @time for epoch=1:o[:epochs]
        dfa_train(w,ao,ho,B,dtrn; lr=o[:lr], epochs=1)
        report(epoch)
    end

    return w
#end

#main()


text:
include("imports.jl")

@testset ExtendedTestSet "Server updates" begin
    server = testserver()
    @test_nowarn addsource!(server, p"hello.md", XNode(:body, [XLeaf("Hello")]))
    @test_nowarn addrewrite!(server, p"hello.md")
    @test_nowarn addbuild!(server, p"hello.md")
    @test !haskey(server.project.sources, p"hello.md")
    @test_nowarn applyupdates!(server)
    @test haskey(server.project.sources, p"hello.md")
    @test haskey(server.project.outputs, p"hello.md")
    @test isfile(joinpath(server.builder.dir, p"hello.md.html"))
end

@testset ExtendedTestSet "FileServer" begin
    dir = Path(mktempdir())
    fs = FileServer(dir)
    start(fs)
    stop(fs)
end

@testset ExtendedTestSet "ServeFiles" begin
    server = testserver()
    mode = ServeFiles()
    @test geteventsource(mode, server, Channel()) isa Pollen.FileServer
end

@testset ExtendedTestSet "ServeFilesLazy" begin
    server = testserver()
    mode = ServeFilesLazy()
    @test geteventsource(mode, server, Channel()) isa Pollen.FileServer
end

@testset ExtendedTestSet "servereventsources" begin
    server = testserver()
    mode = ServeFiles()
    @test Pollen.servereventsources(server, mode, Channel()) |> length == 2
end


text:
using FiberNlse
# Simulation dimension
Nₜ, Nₗ = (1000,1000);

# Fiber properties
D = -17*ps/nm/km; # Dispersion factor
α = 0.046/km; # Losses
γ = 1.1/W/km; # Nonlinear factor
L = 5.0*km; # Fiber length

# Signal properties
T = 40*ps; # Signal duration
λ = 1550*nm; # Wavelength
τ = 5*ps; # Pulse duration
N = 2 # Soliton number
sim,t,l = FiberNlse.configure(Nₜ,Nₗ,D, γ, α, L, T, λ);

# Input construction
Pp = abs((sim.β2/γ/τ^2)*N^2) # Soliton power
Ψₒ = @. sqrt(Pp)/cosh(t/τ) # Soliton formula
FiberNlse.inputSignal(sim,Ψₒ);

FiberNlse.simulate(sim, false); # run the simulation

# Visualization

#= using GLMakie

begin
fig = Figure(resolution=(860, 400), fontsize=12)
ax1 = Axis3(fig[1,1]; aspect=(1,1,1),ylabel="Distance [km]", xlabel="Time [ps]", zlabel="Power [W]")
ax2 = Axis(fig[1,2], ylabel="Distance [km]", xlabel="Time [ps]")
hm = surface!(ax1, t/ps, l/km, abs.(sim.Ψ').^2)
heatmap!(ax2, t/ps, l/km, abs.(sim.Ψ').^2)
Colorbar(fig[1, 3], hm, height=Relative(0.5), label="Power [W]")
fig
end =#


text:
@inline function AbstractRegistry(registry::LocalFolderRegistry)
    return registry
end

@inline function AbstractRegistry(registry::GitRegistry)
    return registry
end

@inline function AbstractRegistry(url::AbstractString)
    registry = GitRegistry(url)
    return registry
end

@inline function make_registry_list(registries::AbstractVector;
                                    include_default_registries::Bool)
    registry_list = Vector{AbstractRegistry}(undef, 0)
    for registry in registries
        push!(registry_list, AbstractRegistry(deepcopy(registry)))
    end
    default_registries = [GitRegistry("https://github.com/JuliaRegistries/General.git")]
    if include_default_registries
        for default_registry in default_registries
            push!(registry_list, AbstractRegistry(default_registry))
        end
    end
    unique!(registry_list)
    return registry_list
end


text:
# ---------------------------------------------------------------

cxx"""
  #include <QOpenGLWidget>
  class MyGLCanvas : public QOpenGLWidget
  {

  public:
    // {};

    MyGLCanvas(int idx) : _idx(idx) {}
    int getidx() { return _idx; }

  protected:
    void paintGL() {
      $:(draw( icxx"return this;" )::Void);
    }

  private:
    int _idx;
  };

  class MyCanvas : public QWidget
  {

  public:
    // {};

    MyCanvas(int idx) : _idx(idx) {}
    int getidx() { return _idx; }

  protected:
    // just call the julia method `draw` and pass the MyCanvas pointer
    void paintEvent(QPaintEvent *event)
    {
      $:(draw( icxx"return this;" )::Void);
    }

  private:
    int _idx;
  };
"""

# ---------------------------------------------------------------

# cxx"""
#   // class MyWindow : public QRasterWindow
#   // {
#   // public:
#   //   MyWindow(int idx, const char* title, int w, int h) {
#   //     std::cout << title << w << h << std::endl;
#   //     m_idx = idx;
#   //     setTitle(title);
#   //     resize(w,h);
#   //   }

#   // protected:

#   //   void render(QPainter *p) {
#   //     std::cout << "here!" << std::endl;
#   //     $:(render( icxx"return p;" )::Void);
#   //   }

#   // private:
#   //   int m_idx;
#   // };
# """

# ---------------------------------------------------------------


text:
struct CallableBondConstantVol{C <: BusinessCalendar, DC <: DayCount, D <: Dates.Period} <: CallableBondVolatilityStructure
  settlementDays::Int
  calendar::C
  volatility::Quote
  dc::DC
  maxBondTenor::D
end

CallableBondConstantVol(settlementDays::Int, cal::BusinessCalendar, fwdYieldVol::Quote, dc::DayCount) =
                      CallableBondConstantVol(settlementDays, cal, fwdYieldVol, dc, Dates.Year(100))


text:
####
#### Sort FunctionalTable by rows.
####

function sort_in_memory(ft::FunctionalTable, ordering::TableOrdering)
    FunctionalTable(sort!(collect(ft), lt = (a, b) -> isless_ordering(ordering, a, b)),
                    TrustOrdering(ordering))
end

function Base.sort(ft::FunctionalTable,
                   column_ordering_specifications::Tuple{Vararg{ColumnOrderingSpecification}})
    sort_ordering = table_ordering(column_ordering_specifications)
    # prefix of an existing order is costless
    is_prefix(sort_ordering, ordering(ft)) && return FunctionalTable(ft, TrustOrdering(sort_ordering))
    # TODO columns which are already sorted need not be resorted
    # TODO columns which do not need to be sorted can be skipped and a permutation applied
    # TODO select appropriate sorting method using heuristics for table size & sorting
    sort_in_memory(ft, sort_ordering)
end

# convenience wrapper
function Base.sort(ft::FunctionalTable,
          column_ordering_specifications::ColumnOrderingSpecification...)
    sort(ft, column_ordering_specifications)
end


text:
module IndianGPA

using SPPL

n = @sppl begin
    Nationality ~ SPPL.Choice([:India => 0.5, :USA => 0.5])
    if Nationality == :India
        Perfect ~ SPPL.Bernoulli(0.1)
        if Perfect == 1
            GPA ~ SPPL.Atomic(10)
        else
            GPA ~ SPPL.Uniform(0, 10)
        end
    else
        Perfect ~ SPPL.Bernoulli(0.15)
        if Perfect == 1
            GPA ~ SPPL.Atomic(4)
        else
            GPA ~ SPPL.Uniform(0, 4)
        end
    end
end

println(n.model.prob(n.Perfect << set(1)))
event = (n.Nationality << set(:USA)) & (n.GPA > 3) | ((8 < n.GPA) < 10)
model_condition = n.model.condition(event)

end # module


text:
# include("Projects/NaquadahBrowser/src/NaquadahBrowser.jl")
using Gtk, Gtk.ShortNames, Cairo, Graphics   # Colors
println(pwd())

if Sys.islinux()
    if basename(pwd()) != "NaquadahBrowser" # Make sure we're in the right dir
        cd("Projects/NaquadahBrowser")
    end
    defaultPage = "file:///src/SamplePages/test.sml"
    global PATH = pwd() * "/src/SamplePages/"
elseif is_windows() # Do something windows-y here... needs tested by someone with a MS machine
    if basename(pwd()) != "NaquadahBrowser" # Make sure we're in the right dir
        cd("Projects\\NaquadahBrowser")
    end
    defaultPage = "file:///src\\SamplePages\\test.sml"
    global PATH = pwd() * "src\\SamplePages\\"
end
#using NaquadahCore
include("NaquadahCore.jl")
# ======================================================================================
# Win()  Create a window
# ======================================================================================
struct Win
    window::Any #Gtk.GtkWindowLeafGtk.GtkCanvas
    canvas::Any #Gtk.GtkCanvas
    controls::Element # Array{Element,1}
    documents::Vector{Page}
    document::Page
    function Win()
        win = Window("Naquadah: An experimental web browser", 1000, 600) # later we can make it save last page size
        canvas = Canvas()
        push!(win, canvas)
        document = Page(defaultPage)
        document.canvas = canvas
        node = document.children[1]
        node.DOM = Dict( ">" => "window", "width" => 1000, "display" => "block", "padding" => [0,0,0,0], "nodes"=>[])
        FetchPage(document, defaultPage)

        document.win = win
        return new(win, canvas, node, [document], document)
    end
end

# ======================================================================================
# This is where drawing starts after a page is fetched.
# ======================================================================================
function DrawANode(document::Page)
    c = document.canvas
    node = document.children[1] #.children[2]
    page = node.children[3]

        #  This is efectively the resize event.
           @guarded draw(c) do widget

                ScrollY = 0.0
                ctx = getgc(c)
                document.height::Float64   = height(c)
                document.width::Float64    = width(c)

                if page.scroll.y < 0 # Don't scroll above zero/top
                    ScrollY = page.scroll.y
                    page.scroll.y = 0.0
                    VmoveAllChildren(page, abs(ScrollY), false) # FROM: LayoutBuild.jl
                end

                setUpWindow(document, document.width, document.height)    # FROM: LayoutBuild.jl
                allAtributesToLayout(document, node)
                AttatchEvents(document, c)            # FROM: Events.jl
                CreateLayoutTree(document, node)      # FROM: LayoutBegin.jl
                page.scroll.y = ScrollY
                VmoveAllChildren(page, ScrollY, false)
                DrawViewport(ctx, document, node)     # FROM: GraphDraw.jl
           end
    show(c);
end
# ======================================================================================

function main()
    win = Win(); # Win(controls)
    DrawANode(win.document)
end

main()