text:
#include("./src/01_species_list_creation.jl")
#include("./src/04_eqn_component_functions.jl")
#include("./src/05_network_metrics.jl")
#include("./src/08_perturbation_sims_functions.jl")
using DifferentialEquations
using DelimitedFiles

# set all positional values
organism_id_pos = 1
organism_type_pos = 2
body_size_pos = 3
lower_prey_limit_pos = 4
upper_prey_limit_pos = 5
habitat_midpoint_pos = 6
habitat_lower_limit_pos = 7
habitat_upper_limit_pos = 8
nutritional_value_pos = 9
biomass_pos = 10
network_id_pos = 11
alive_pos = 12
# pulse perturbation
#
# only save some things

"""
    perturbation_eq_sims!(species_list, network,
                            sim_length,
                            network_id,
                            perturbation_type
                            )

Function runs through simulations from end of initial sims to 15000, emplying
perturbations called separately.
"""
function perturbation_eq_sims!(
    species_list::Array,
    sim_length::Int,
    network_id::Int,
    perturbation_type::String
    )

    # allocate x0
    x0 = Float64[]
    for i in 1:size(species_list, 1)
        append!(x0, species_list[i, biomass_pos])
    end

    # keep track of how many extinctions there were
    extinct_sp = Float64[]

    # calculate perturbation values --- c = z*u*(1/f)
    c = 0.5
    if perturbation_type == "pulse"
        u = 100
        f = 1
        z = c/(u*(1/f))
        perturbation_times = 1001:1100
        transient_times = 1101:1251
        normal_times = 1252:15000

        # loop through perturbation times
        time_start = 1000.0
        while time_start =< perturbation_times[end]

            # put species list and network in same object
            p = [species_list, network]

            # set time and defie ODE problem

            #### THIS IS SO INEFFICIENT - figure out a way to take it to zero at each step

            ### should be able to do all the things i want via callbacks including
                # - changing values to zero
                # - changing the parameters

            tspan = (time_start, time_start+9.9)
            prob = ODEProblem(biomass!, x0, tspan, p)
            sol = solve(prob)

            # keep only last timestep and check for extinction
            x0 = sol.u[end]
            check_extinctions!(x0)
            species_list[:,biomass_pos] = x0

            # get new network
            network = create_interaction_matrix(species_list)

            # change alive to extinct if species is extinct
            for i in 1:size(species_list, 1)
                if species_list[i, biomass_pos] == 0
                    species_list[i, alive_pos] = 0
                    append!(extinct_sp, species_list[i, organism_id_pos])
                end
            end

            time_start += 10.0
        end

    elseif perturbation_type == "periodic pulse"
        u = 100
        f = 14
        z = c/(u*(1/f))
        perturbation_times = []
        transient_times = []
        for i in range(1000,14000, step = 1000)
            append!(perturbation_times, i:i+99)
            append!(transient_times, i+100:i+250)
        end
        normal_times = []
        for i in 1001:15000
            if i ∉ perturbation_times
                append!(normal_times, i)
            end
        end
    elseif perturbation_type == "press"
        u = 14000
        f = 1
        z = c/(u*(1/f))
        perturbation_times = 1001:15000
        normal_times = []
    end

    # bind back the network_id and write the species list
    species_list = hcat(species_list, fill(network_id, 50, 1))
    writedlm("./data/intermediate-files/species-lists-pre-perturbation/species_list_$network_id.csv", species_list, ",")

    # create and write network
    network = create_interaction_matrix(species_list)
    network = hcat(fill(network_id, 50, 1), network)
    writedlm("./data/intermediate-files/networks-pre-perturbation/network_$network_id.csv", network, ",")

end

# prep network
# run, save random run number 5 <= x <= 15
# based on the multiple types, figure out what time points to change the perturbation amount
# do pulse one first

# save all initial parameters (maybe a seaprate object)

# sims

function perturbation_simulation(
    species_list:::Array,
    network::Array,
    type_of_perturbation::String,
    x::Array
    )

# pre-allocate memory for the objects we want to save from the time series
# save initial parameters
# sim from the starting biomass out to 5000 - save that point
# ^ maybe every 10 points stop it and check for extinction in any of them
# ^ also, add a component to the species_list to keep track of alive (1) or
#    extinct (0)
# check the shape of the web at that point & save
# sim 10 at a time, check for extinctions, check for perturbation ,
#    get the parameters again

# read in species list
# make network

# determine what kind of perturbation (aka set parameters)
# figure out what time points are going to be what
# do perturbation for x time
# save transients
# after transient time, figure out how long until next perturbation or till the end
# pick 5 sort of equal points in that time and save those

# set all positional values
organism_id_pos = 1
organism_type_pos = 2
body_size_pos = 3
lower_prey_limit_pos = 4
upper_prey_limit_pos = 5
habitat_midpoint_pos = 6
habitat_lower_limit_pos = 7
habitat_upper_limit_pos = 8
nutritional_value_pos = 9
biomass_pos = 10
network_id_pos = 11
alive_pos = 12

# set parameters - pulse
c = 0.5
x = 100
f = 1
z = c/(x*(1/f))

# put species list and network in same object
# extract all values to be used w/in function from the p value
# species_list, network, base_mortality, amount = p
p = [species_list, network, 0.3, z]

# define times
tspan = (0.0, 50.0)

# define the time stops
const tstop = [5.0;8.0]

#  allocate x0
x0 = Float64[]
for i in 1:size(species_list, 1)
    append!(x0, species_list[i, biomass_pos])
end

prob = ODEProblem(biomass_perturb, x0, tspan, p)
sol = solve(prob,Tsit5(),callback = cbs, tstops=tstop)
plot(sol)

# set parameters - pulse
c = 0.5
u = 100
f = 1
z = c/(u*(1/f))

# set constant stop values
const stop1 = [5.0]
const stop2 = [8.0]

# callback to check for extinctions

# callback to do pulse perturbation
# increase mortality (start pulse)
function inc_mort_cond(u,t,integrator)
    t in stop1
end

# decrease mortality (stop pulse)
function dec_mort_cond(u,t,integrator)
    t in stop2
end

# apply an affect when condition is reached
function affect_inc!(integrator)
    integrator.p[4] = 1.0 + z
end

function affect_dec!(integrator)
    integrator.p[4] = 1.0
end

save_positions = (true, true)
cb_inc = DiscreteCallback(inc_mort_cond, affect_inc!, save_positions=save_positions)
save_positions = (false, true)
cb_dec = DiscreteCallback(dec_mort_cond, affect_inc!, save_positions=save_positions)

# put all callbacks into the callback set
cb = TerminateSteadyState(1e-12, 1e-9, DiffEqCallbacks.allDerivPass)
cbs = CallbackSet(cb_inc, cb_dec, cb)

# put species list and network in same object
p = [species_list, network, 0.3, z]

# define times
tspan = (0.0, 100.0)

# define the time stops
const tstop = [5.0;8.0]

#  allocate x0
x0 = Float64[]
for i in 1:size(species_list, 1)
    append!(x0, species_list[i, biomass_pos])
end

prob = ODEProblem(biomass_perturb, u0, tspan, p)
sol = solve(prob,Tsit5(),callback = cbs, tstops=tstop)
plot(sol)

sol.u

perturb_sp = rand(1:50, 17)
for i ∈ perturb_sp
    print(3*i)
end
for i ∉ perturb_sp & ∈ 1:50
    print(i)
end

all =Array(1:50)

intersect(perturb_sp, all)
setdiff(all, perturb_sp) #order matters here

"""
    biomass_perturb(du, u, p, t)

Calling biomass functions, create the actual equation for each species, to be
passed to the simulation. Do this by randomly picking ~1/3
"""
function biomass_perturb(du, u, p, t)
    # pick set of random species to perturb
    all = Array(1:size(species_list, 1))
    perturb_sp = rand(Array(1:size(species_list, 1)),
                        Int(round(size(species_list, 1)*(1/3))))
    # extract all values to be used w/in function from the p value
    species_list, network, base_mortality, amount = p
    # population du
    for i in 1:size(species_list, 1)
        # if species to be perturbed, pass amount, else pass 1 as amount
        if i ∈ perturb_sp
            du[i] = growth_rate(i, species_list, u) +
                    biomass_gain(species_list, network, u, 2, i) -
                    biomass_loss(species_list, network, u, 2, i) -
                    mortality_change(i, u, base_mortality, amount)
        else
            du[i] = growth_rate(i, species_list, u) +
                    biomass_gain(species_list, network, u, 2, i) -
                    biomass_loss(species_list, network, u, 2, i) -
                    mortality_change(i, u, base_mortality, 1)
        end
    end
end

# testing getting any array into a Matrix

vector = [[1.2, 1.3, 1.4],[1.2, 1.5, 1.7]]
matrix = zeros(2,length(vector[1]))
for i in 1:size(matrix,1)
    matrix[i,:] = vector[i]
end
length(vector[1])

max_stars_repo_name: colebrookson/trait-based-rewiring

max_issues_repo_name: colebrookson/trait-based-rewiring

max_forks_repo_name: colebrookson/trait-based-rewiring


text:
using Test

using TimeSeries: TimeArray, timestamp
using TimeSeriesResampler: resample, TimeFrame, ohlc, mean, sum, std
using TimeFrames
using Dates

function variation(a; n=1)
    a[1+n:end] - a[1:end-n]
end

@testset "resample" begin
    # Define a sample timeseries (prices for example)
    idx = DateTime(2010,1,1):Dates.Minute(1):DateTime(2010,4,1)
    idx = idx[1:end-1]
    N = length(idx)
    y1 = rand(-1.0:0.01:1.0, N)
    y1 = 1000 .+ cumsum(y1)
    y2 = rand(-1.0:0.01:1.0, N)
    y2 = 500 .+ cumsum(y2)

    ta1 = TimeArray(collect(idx), y1, [:y1])

    #df = DataFrame(Date=idx, y=y)
    ta2 = TimeArray(collect(idx), hcat(y1, y2), [:y1, :y2])
    #println("ta=")
    #println(ta)

    a_ta = [ta1, ta2]

    # Define how datetime should be grouped (timeframe)
    a_tf = [
        TimeFrame(dt -> floor(dt, Dates.Minute(15))),  # using a lambda function
        TimeFrame(Minute(15)),  # using a TimeFrame object (from TimeFrames.jl)
        TimeFrame("15T"),  # using a string TimeFrame shortcut to create a TimeFrame
        "15T",  # using a string TimeFrame shortcut
    ]

    for ta in a_ta
        #println("ta=")
        #println(ta)

        for tf in a_tf
            println(tf)

            # resample using OHLC values
            ta_ohlc = ohlc(resample(ta, tf))
            #println("ta_ohlc=")
            #println(ta_ohlc)
            @test mean(variation(timestamp(ta_ohlc))) == Dates.Minute(15)

            ## group-by by 1 column
            ta_ohlc = ohlc(resample(ta2, tf)[:y1])
            #println("ta_ohlc=")
            #println(ta_ohlc)
            @test mean(variation(timestamp(ta_ohlc))) == Dates.Minute(15)

            ## group-by by 2 columns
            ta_ohlc = ohlc(resample(ta2, tf)[:y1, :y2])
            #println("ta_ohlc=")
            #println(ta_ohlc)
            @test mean(variation(timestamp(ta_ohlc))) == Dates.Minute(15)

            # resample using mean values
            ta_mean = mean(resample(ta, tf))
            #println("ta_mean=")
            #println(ta_mean)
            @test mean(variation(timestamp(ta_mean))) == Dates.Minute(15)

            # Define an other sample timeseries (volume for example)
            vol = rand(0:0.01:1.0, N)
            ta_vol = TimeArray(collect(idx), vol, [:vol])
            #println("ta_vol=")
            #println(ta_vol)
            @test mean(variation(timestamp(ta_vol))) == Dates.Minute(1)

            # resample using sum values
            ta_vol_sum = sum(resample(ta_vol, tf))
            #println("ta_vol_sum=")
            #println(ta_vol_sum)
            @test mean(variation(timestamp(ta_vol_sum))) == Dates.Minute(15)

            # resample using std values
            ta_vol_std = std(resample(ta_vol, tf))
            #println("ta_vol_std=")
            #println(ta_vol_std)
            @test mean(variation(timestamp(ta_vol_std))) == Dates.Minute(15)
        end
    end
end

ext: jl


text:
using Pda
using Test

@testset "Pda.jl" begin
    # Write your tests here.
end


text:
module GenerateFiles

include("../src/CustomTypes.jl")
include("QNMFileUtils.jl")
using  .CustomTypes
import .QNMFileUtils as QF

export generate

function generate(
      prename::String,
      nr::myI,
      nl::myI,
      s::myI,
      m::myI,
      l::myI,
      n::myI,
      avals::Vector{myF},
      qnm_tolerance::myF,
      coef_tolerance::myF,
      epsilon::myF)
   #println("bash -c 'rm -f $(pwd)/HigherPrecData/*.h5'")
   #run(`bash -c 'rm -f $(pwd)/HigherPrecData/*.h5'`)

   lmin = max(abs(s),abs(m))
   for l=[lmin,lmin+1,lmin+2]
      try
         println("s=$s,m=$m,n=$n,l=$l")
            QF.generate_data(
               "$(pwd())/qnmfiles/$(prename)s$(s)_m$(m)_n$(n)",
               tomyI(nr), tomyI(nl),
               tomyI(s), tomyI(n), tomyI(l), tomyI(m),
               [tomyF(a) for a in avals],
               qnm_tolerance=tomyF(qnm_tolerance),
               coef_tolerance=tomyF(coef_tolerance),
               epsilon=tomyF(epsilon)
            )
      catch err
         println(err)
         continue
      end
   end
   return nothing
end

end


text:
module BenchParallelTriangularSum

using BenchmarkTools
using FLoops

# This is a rather contrived example for exploring the problems with
# which `nestlevel > 1` is useful.
function sum_triangular(exc, xs, f = f)
    init = zero(eltype(xs))
    @floop exc for a0 in xs,
        a1 in f(a0),
        a2 in f(a1),
        a3 in f(a2),
        a4 in f(a3),
        a5 in f(a4),
        a6 in f(a4),
        a7 in f(a4),
        a8 in f(a4),
        a9 in f(a4)

        m = a0 * a1 * a2 * a3 * a4 * a5 * a6 * a7 * a8 * a9
        @reduce(s = init + m)
    end
    return s
end

const SUITE = BenchmarkGroup()

f(x) = 1:x
xs = 1:8
SUITE["threaded"] = @benchmarkable sum_triangular(ThreadedEx(), $xs)
let s1 = SUITE["threaded_nest"] = BenchmarkGroup()
    for nestlevel in 1:3
        s1[:nestlevel=>nestlevel] = @benchmarkable sum_triangular(
            ThreadedEx(basesize = 1, nestlevel = $(Val(nestlevel))),
            $xs,
        )
    end
end
SUITE["sequential"] = @benchmarkable sum_triangular(SequentialEx(), $xs)

end  # module
BenchParallelTriangularSum.SUITE


text:
# This type is for testing the escaping of quotes in HTML output. It needs to be
# outside of the module so that Julia 1.0 would not complain about the type
# being defined in local scope.
struct QuoteTestType{T} end

module TestIO

using Test, DataFrames, CategoricalArrays, Dates, Markdown
import Main: QuoteTestType

# Test LaTeX export
@testset "LaTeX export" begin
    df = DataFrame(A = Int64.( 1:4 ),
                B = ["\$10.0", "M&F", "A~B", "\\alpha"],
                C = ["A", "B", "C", "S"],
                D = [1.0, 2.0, missing, 3.0],
                E = CategoricalArray(["a", missing, "c", "d"]),
                F = Vector{String}(undef, 4),
                G = [ md"[DataFrames.jl](http://juliadata.github.io/DataFrames.jl)", md"###A", md"``\frac{A}{B}``", md"*A*b**A**"]
                )
    str = """
        \\begin{tabular}{r|ccccccc}
        \t& A & B & C & D & E & F & G\\\\
        \t\\hline
        \t& Int64 & String & String & Float64? & Cat…? & String & MD\\\\
        \t\\hline
        \t1 & 1 & \\\$10.0 & A & 1.0 & a & \\emph{\\#undef} & \\href{http://juliadata.github.io/DataFrames.jl}{DataFrames.jl} \\\\
        \t2 & 2 & M\\&F & B & 2.0 & \\emph{missing} & \\emph{\\#undef} & \\#\\#\\#A \\\\
        \t3 & 3 & A\\textasciitilde{}B & C & \\emph{missing} & c & \\emph{\\#undef} & \$\\frac{A}{B}\$ \\\\
        \t4 & 4 & \\textbackslash{}\\textbackslash{}alpha & S & 3.0 & d & \\emph{\\#undef} & \\emph{A}b\\textbf{A} \\\\
        \\end{tabular}
        """

    @test repr(MIME("text/latex"), df) == str
    @test repr(MIME("text/latex"), eachcol(df)) == str
    @test repr(MIME("text/latex"), eachrow(df)) == str

    @test_throws ArgumentError DataFrames._show(stdout, MIME("text/latex"),
                                                DataFrame(ones(2,2), :auto), rowid=10)
end

@testset "Huge LaTeX export" begin
    df = DataFrame(a=1:1000)
    ioc = IOContext(IOBuffer(), :displaysize => (10, 10), :limit => false)
    show(ioc, "text/latex", df)
    @test length(String(take!(ioc.io))) > 10000

    ioc = IOContext(IOBuffer(), :displaysize => (10, 10), :limit => true)
    show(ioc, "text/latex", df)
    @test length(String(take!(ioc.io))) < 10000
end

#Test HTML output for IJulia and similar
@testset "HTML output" begin
    df = DataFrame(Fish = ["Suzy", "Amir"], Mass = [1.5, missing])
    io = IOBuffer()
    show(io, "text/html", df)
    str = String(take!(io))
    @test str == "
" *
                 "
2 rows × 2 columns
" *
                 "
" *
                 "

	Fish	Mass
" *
                 "
	String	" *
                 "Float64?

" *
                 "

1	Suzy	1.5
" *
                 "
2	Amir	missing
" *
                 "

"

    df = DataFrame(Fish = Vector{String}(undef, 2), Mass = [1.5, missing])
    io = IOBuffer()
    show(io, "text/html", df)
    str = String(take!(io))
    @test str == "
" *
                 "
2 rows × 2 columns
" *
                 "
" *
                 "

	Fish	Mass
" *
                 "
	String	" *
                 "Float64?

" *
                 "

1	#undef	1.5
" *
                 "
2	#undef	missing
" *
                 "

"

    io = IOBuffer()
    show(io, "text/html", eachrow(df))
    str = String(take!(io))
    @test str == "
2×2 DataFrameRows
" *
                 "

" *
                 "	Fish	Mass
" *
                 "
	String	Float64?

" *
                 "
1	#undef	1.5
" *
                 "
2	#undef	missing

"

    io = IOBuffer()
    show(io, "text/html", eachcol(df))
    str = String(take!(io))
    @test str == "
2×2 DataFrameColumns
" *
                 "

" *
                 "	Fish	Mass
" *
                 "
	String	Float64?

" *
                 "
1	#undef	1.5
" *
                 "
2	#undef	missing

"

    io = IOBuffer()
    show(io, "text/html", df[1, :])
    str = String(take!(io))
    @test str == "
DataFrameRow (2 columns)
" *
                 "

" *
                 "

	Fish	Mass

	" *
                 "String	Float64?

1	" *
                 "#undef	1.5

"

    io = IOBuffer()
    show(io, MIME"text/html"(), df, summary=false)
    str = String(take!(io))
    @test str == "

" *
                 "	Fish	Mass
" *
                 "
	String	Float64?

" *
                 "
1	#undef	1.5
" *
                 "
2	#undef	missing

"

    io = IOBuffer()
    show(io, MIME"text/html"(), eachrow(df), summary=false)
    str = String(take!(io))
    @test str == "

" *
                 "	Fish	Mass
" *
                 "
	String	Float64?

" *
                 "
1	#undef	1.5
" *
                 "
2	#undef	missing

"

    io = IOBuffer()
    show(io, MIME"text/html"(), eachcol(df), summary=false)
    str = String(take!(io))
    @test str == "

" *
                 "	Fish	Mass
" *
                 "
	String	Float64?

" *
                 "
1	#undef	1.5
" *
                 "
2	#undef	missing

"

    io = IOBuffer()
    show(io, MIME"text/html"(), df[1, :], summary=false)
    str = String(take!(io))
    @test str == "

	Fish	" *
                 "Mass

	String	Float64?

" *
                 "

1	#undef	1.5

"

    io = IOBuffer()
    show(IOContext(io, :limit => true, :displaysize => (10, 10)), MIME"text/html"(),
         DataFrame(Int64[1 2 3 4 5 6 7 8 9], :auto))
    str = String(take!(io))
    @test str == "

1 rows × 9 columns (omitted printing of 7 columns)
" *
                 "

	x1	x2

	" *
                 "Int64	Int64

" *
                 "

1	1	2

"

    @test_throws ArgumentError DataFrames._show(stdout, MIME("text/html"),
                                                DataFrame(ones(2,2), :auto), rowid=10)

    df = DataFrame(
        A=Int64[1,4,9,16],
        B = [
            md"[DataFrames.jl](http://juliadata.github.io/DataFrames.jl)",
            md"###A",
            md"``\frac{A}{B}``",
            md"*A*b**A**" ]
    )

    @test repr(MIME("text/html"), df) ==
    "

4 rows × 2 columns
" *
    "

	A	B

	" *
    "Int64	MD

" *
    "

1	1	" *
    "

DataFrames.jl
\n

" *
    "
2	4

###A
\n

" *
    "
3	9

$\\frac{A}{B}$
\n

" *
    "
4	16

AbA
\n

" *
    "

"

    # Test that single and double quotes get escaped properly
    df = DataFrame(
        xs = ["'", "\"", "<foo>'</bar>"],
        ys = [QuoteTestType{'\''}(), QuoteTestType{'"'}, QuoteTestType{Symbol("\"'")}()],
        zs = QuoteTestType{'"'}[QuoteTestType{'"'}(), QuoteTestType{'"'}(), QuoteTestType{'"'}()],
    )
    io = IOBuffer()
    show(io, "text/html", df)
    str = String(take!(io))
    @test str ==
        "

3 rows × 3 columns
" *
        "

" *
            "
" *
                "	" *
                "xs	" *
                "ys	" *
                "zs	" *
            "

" *
                "	" *
                "String	" *
                "Any	" *
                "QuoteTes…	" *
            "
" *
        "

" *
            "
" *
                "1	" *
                "'	" *
                "QuoteTestType{'\\\\''}()	" *
                "QuoteTestType{'"'}()	" *
            "

" *
                "2	" *
                ""	" *
                "QuoteTestType{'"'}	" *
                "QuoteTestType{'"'}()	" *
            "

" *
                "3	" *
                "<foo>'</bar>	" *
                "QuoteTestType{Symbol("\\\\"'")}()	" *
                "QuoteTestType{'"'}()	" *
            "
" *
        "

"
end

# test limit attribute of IOContext is used
@testset "limit attribute" begin
    df = DataFrame(a=1:1000)
    ioc = IOContext(IOBuffer(), :displaysize => (10, 10), :limit => false)
    show(ioc, "text/html", df)
    @test length(String(take!(ioc.io))) > 10000

    ioc = IOContext(IOBuffer(), :displaysize => (10, 10), :limit => true)
    show(ioc, "text/html", df)
    @test length(String(take!(ioc.io))) < 10000
end

@testset "printtable" begin
    df = DataFrame(A = 1:3,
                B = 'a':'c',
                C = ["A", "B", "C"],
                D = CategoricalArray(string.('a':'c')),
                E = CategoricalArray(["A", "B", missing]),
                F = Vector{Union{Int, Missing}}(1:3),
                G = missings(3),
                H = fill(missing, 3))

    @test sprint(DataFrames.printtable, df) ==
        """
        "A","B","C","D","E","F","G","H"
        1,"a","A","a","A","1",missing,missing
        2,"b","B","b","B","2",missing,missing
        3,"c","C","c",missing,"3",missing,missing
        """
end

@testset "csv/tsv output" begin
    df = DataFrame(a = [1,2], b = [1.0, 2.0])

    for x in [df, eachcol(df), eachrow(df)]
        @test sprint(show, "text/csv", x) == """
        "a","b"
        1,1.0
        2,2.0
        """
        @test sprint(show, "text/tab-separated-values", x) == """
        "a"\t"b"
        1\t1.0
        2\t2.0
        """
    end
end

@testset "Markdown as text/plain and as text/csv" begin
    df = DataFrame(
        A=Int64[1,4,9,16,25,36,49,64],
        B = [
            md"[DataFrames.jl](http://juliadata.github.io/DataFrames.jl)",
            md"``\frac{x^2}{x^2+y^2}``",
            md"# Header",
            md"This is *very*, **very**, very, very, very, very, very, very, very long line" ,
            md"",
            Markdown.parse("∫αγ∞1∫αγ∞2∫αγ∞3∫αγ∞4∫αγ∞5∫αγ∞6∫αγ∞7∫αγ∞8∫αγ∞9∫αγ∞0∫αγ∞1∫αγ∞2∫αγ∞3"),
            Markdown.parse("∫αγ∞1∫αγ∞\n"*
                           "  * 2∫αγ∞3∫αγ∞4\n"*
                           "  * ∫αγ∞5∫αγ\n"*
                           "  * ∞6∫αγ∞7∫αγ∞8∫αγ∞9∫αγ∞0"),
            Markdown.parse("∫αγ∞1∫αγ∞2∫αγ∞3∫αγ∞4∫αγ∞5∫αγ∞6∫αγ∞7∫αγ∞8∫αγ∞9∫αγ∞0∫α\n"*
                           "  * γ∞1∫α\n"*
                           "  * γ∞2∫αγ∞3∫αγ∞4∫αγ∞5∫αγ∞6∫αγ∞7∫αγ∞8∫αγ∞9∫αγ∞0"),
        ]
    )
    if VERSION < v"1.6.0-DEV"
        @test sprint(show, "text/plain", df) ==
            """
            8×2 DataFrame
             Row │ A      B
                 │ Int64  MD
            ─────┼──────────────────────────────────────────
               1 │     1    DataFrames.jl (http://juliadat…
               2 │     4    \\frac{x^2}{x^2+y^2}
               3 │     9    Header\\n  ≡≡≡≡≡≡≡≡
               4 │    16    This is very, very, very, very…
               5 │    25
               6 │    36    ∫αγ∞1∫αγ∞2∫αγ∞3∫αγ∞4∫αγ∞5∫αγ∞6…
               7 │    49    ∫αγ∞1∫αγ∞\\n\\n    •    2∫αγ∞3∫α…
               8 │    64    ∫αγ∞1∫αγ∞2∫αγ∞3∫αγ∞4∫αγ∞5∫αγ∞6…"""
    else
        @test sprint(show, "text/plain", df) ==
            """
            8×2 DataFrame
             Row │ A      B
                 │ Int64  MD
            ─────┼──────────────────────────────────────────
               1 │     1    DataFrames.jl (http://juliadat…
               2 │     4    \\frac{x^2}{x^2+y^2}
               3 │     9    Header\\n  ≡≡≡≡≡≡≡≡
               4 │    16    This is very, very, very, very…
               5 │    25
               6 │    36    ∫αγ∞1∫αγ∞2∫αγ∞3∫αγ∞4∫αγ∞5∫αγ∞6…
               7 │    49    ∫αγ∞1∫αγ∞\\n\\n    •  2∫αγ∞3∫αγ∞…
               8 │    64    ∫αγ∞1∫αγ∞2∫αγ∞3∫αγ∞4∫αγ∞5∫αγ∞6…"""
    end

    @test sprint(show, "text/csv", df) ==
        """
        \"A\",\"B\"
        1,\"[DataFrames.jl](http://juliadata.github.io/DataFrames.jl)\"
        4,\"\$\\\\frac{x^2}{x^2+y^2}\$\"
        9,\"# Header\"
        16,\"This is *very*, **very**, very, very, very, very, very, very, very long line\"
        25,\"\"
        36,\"∫αγ∞1∫αγ∞2∫αγ∞3∫αγ∞4∫αγ∞5∫αγ∞6∫αγ∞7∫αγ∞8∫αγ∞9∫αγ∞0∫αγ∞1∫αγ∞2∫αγ∞3\"
        49,\"∫αγ∞1∫αγ∞\\n\\n  * 2∫αγ∞3∫αγ∞4\\n  * ∫αγ∞5∫αγ\\n  * ∞6∫αγ∞7∫αγ∞8∫αγ∞9∫αγ∞0\"
        64,\"∫αγ∞1∫αγ∞2∫αγ∞3∫αγ∞4∫αγ∞5∫αγ∞6∫αγ∞7∫αγ∞8∫αγ∞9∫αγ∞0∫α\\n\\n  * γ∞1∫α\\n  * γ∞2∫αγ∞3∫αγ∞4∫αγ∞5∫αγ∞6∫αγ∞7∫αγ∞8∫αγ∞9∫αγ∞0\"
        """
end

@testset "Markdown as HTML" begin
    df = DataFrame(
        A=Int64[1,4,9,16,25,36,49,64],
        B = [
            md"[DataFrames.jl](http://juliadata.github.io/DataFrames.jl)",
            md"``\frac{x^2}{x^2+y^2}``",
            md"# Header",
            md"This is *very*, **very**, very, very, very, very, very, very, very long line" ,
            md"",
            Markdown.parse("∫αγ∞1∫αγ∞2∫αγ∞3∫αγ∞4∫αγ∞5∫αγ∞6∫αγ∞7∫αγ∞8∫αγ∞9∫αγ∞0" *
                           "∫αγ∞1∫αγ∞2∫αγ∞3∫αγ∞4∫αγ∞5∫αγ∞6∫αγ∞7∫αγ∞8∫αγ∞9∫αγ∞0"),
            Markdown.parse("∫αγ∞1∫αγ∞2∫αγ∞3∫αγ∞4∫αγ∞5∫αγ∞6∫αγ\n"*
                           "  * ∞7∫αγ\n"*
                           "  * ∞8∫αγ\n"*
                           "  * ∞9∫αγ∞0∫α\nγ∞1∫αγ∞2∫αγ∞3∫αγ∞4∫αγ∞5∫αγ∞6∫αγ∞7∫αγ∞8∫αγ∞9∫αγ∞0"),
            Markdown.parse("∫αγ∞1∫αγ∞2∫αγ∞3∫αγ∞4∫αγ∞5∫αγ∞6∫αγ∞7∫αγ∞8∫αγ∞9∫αγ∞0∫α\n"*
                           "  * γ∞1∫α\n"*
                           "  * γ∞2∫αγ∞3∫αγ∞4∫αγ∞5∫αγ∞6∫αγ∞7∫αγ∞8∫αγ∞9∫αγ∞0"),
        ]
    )
    @test sprint(show,"text/html",df) ==
        "

8 rows × 2 columns
" *
        "

" *
            "
	A	B
" *
            "
	Int64	MD
" *
        "
" *
        "
" *
        "
1	1
" *
            "
DataFrames.jl
\n

" *
        "
2	4

$\\frac{x^2}{x^2+y^2}$
\n

" *
        "
3	9

Header
\n

" *
        "
4	16
" *
            "
This is very, very, very, very, very, very, very, very, very long line
\n" *
        "

" *
        "
5	25

" *
        "
6	36
" *
            "
∫αγ∞1∫αγ∞2∫αγ∞3∫αγ∞4∫αγ∞5∫αγ∞6∫αγ∞7∫αγ∞8∫αγ∞9∫αγ∞0∫αγ∞1∫αγ∞2∫αγ∞3∫αγ∞4∫αγ∞5∫αγ∞6∫αγ∞7∫αγ∞8∫αγ∞9∫αγ∞0
\n" *
        "

" *
        "
7	49
" *
            "
∫αγ∞1∫αγ∞2∫αγ∞3∫αγ∞4∫αγ∞5∫αγ∞6∫αγ
\n
\n

∞7∫αγ
\n
\n

∞8∫αγ
\n
\n

∞9∫αγ∞0∫α
\n
\n
\n
γ∞1∫αγ∞2∫αγ∞3∫αγ∞4∫αγ∞5∫αγ∞6∫αγ∞7∫αγ∞8∫αγ∞9∫αγ∞0
\n" *
        "

" *
        "
8	64
" *
            "
∫αγ∞1∫αγ∞2∫αγ∞3∫αγ∞4∫αγ∞5∫αγ∞6∫αγ∞7∫αγ∞8∫αγ∞9∫αγ∞0∫α
" *
            "\n
\n" *
                "

γ∞1∫α
\n
\n" *
                "

γ∞2∫αγ∞3∫αγ∞4∫αγ∞5∫αγ∞6∫αγ∞7∫αγ∞8∫αγ∞9∫αγ∞0
\n
\n" *
        "
\n" * "

"
end

@testset "empty data frame and DataFrameRow" begin
    df = DataFrame(a = [1,2], b = [1.0, 2.0])

    @test sprint(show, "text/csv", df[:, 2:1]) == ""
    @test sprint(show, "text/tab-separated-values", df[:, 2:1]) == ""
    @test sprint(show, "text/html", df[:, 2:1]) ==
          "

0 rows × 0 columns
" *
          "

" *
          "

"
    @test sprint(show, "text/latex", df[:, 2:1]) ==
          "\\begin{tabular}{r|}\n\t& \\\\\n\t\\hline\n\t& \\\\\n\t\\hline\n\\end{tabular}\n"

    @test sprint(show, "text/csv", @view df[:, 2:1]) == ""
    @test sprint(show, "text/tab-separated-values", @view df[:, 2:1]) == ""
    @test sprint(show, "text/html", @view df[:, 2:1]) ==
          "

0 rows × 0 columns
" *
          "

" *
          "

"
    @test sprint(show, "text/latex", @view df[:, 2:1]) ==
          "\\begin{tabular}{r|}\n\t& \\\\\n\t\\hline\n\t& \\\\\n\t\\hline\n\\end{tabular}\n"

    @test sprint(show, "text/csv", df[1, 2:1]) == ""
    @test sprint(show, "text/tab-separated-values", df[1, 2:1]) == ""
    @test sprint(show, "text/html", df[1, 2:1]) ==
          "
DataFrameRow (0 columns)

" *
          "

"
    @test sprint(show, "text/latex", df[1, 2:1]) ==
          "\\begin{tabular}{r|}\n\t& \\\\\n\t\\hline\n\t& \\\\\n\t\\hline\n\\end{tabular}\n"
end

@testset "consistency" begin
    df = DataFrame(a = [1, 1, 2, 2], b = [5, 6, 7, 8], c = 1:4)
    push!(df.c, 5)
    @test_throws AssertionError sprint(show, "text/html", df)
    @test_throws AssertionError sprint(show, "text/latex", df)
    @test_throws AssertionError sprint(show, "text/csv", df)
    @test_throws AssertionError sprint(show, "text/tab-separated-values", df)

    df = DataFrame(a = [1, 1, 2, 2], b = [5, 6, 7, 8], c = 1:4)
    push!(DataFrames._columns(df), df[:, :a])
    @test_throws AssertionError sprint(show, "text/html", df)
    @test_throws AssertionError sprint(show, "text/latex", df)
    @test_throws AssertionError sprint(show, "text/csv", df)
    @test_throws AssertionError sprint(show, "text/tab-separated-values", df)
end

@testset "summary tests" begin
    df = DataFrame(ones(2,3), :auto)

    for (v, s) in [(df, "2×3 DataFrame"),
                   (view(df, :, :), "2×3 SubDataFrame"),
                   (df[1, :], "3-element DataFrameRow"),
                   (DataFrames.index(df), "data frame with 3 columns"),
                   (DataFrames.index(df[1:1, 1:1]), "data frame with 1 column"),
                   (DataFrames.index(view(df, 1:1, 1:1)), "data frame with 1 column"),
                   (eachrow(df), "2-element DataFrameRows"),
                   (eachcol(df), "3-element DataFrameColumns")]
        @test summary(v) == s
        io = IOBuffer()
        summary(io, v)
        @test String(take!(io)) == s
    end
end

@testset "eltypes tests" begin
    df = DataFrame(A = Int32.(1:3), B = ["x", "y", "z"])

    io = IOBuffer()
    show(io, MIME("text/plain"), df, eltypes=true)
    str = String(take!(io))
    @test str == """
        3×2 DataFrame
         Row │ A      B
             │ Int32  String
        ─────┼───────────────
           1 │     1  x
           2 │     2  y
           3 │     3  z"""

    io = IOBuffer()
    show(io, MIME("text/plain"), eachcol(df), eltypes=true)
    str = String(take!(io))
    @test str == """
        3×2 DataFrameColumns
         Row │ A      B
             │ Int32  String
        ─────┼───────────────
           1 │     1  x
           2 │     2  y
           3 │     3  z"""

    io = IOBuffer()
    show(io, MIME("text/plain"), eachrow(df), eltypes=true)
    str = String(take!(io))
    @test str == """
        3×2 DataFrameRows
         Row │ A      B
             │ Int32  String
        ─────┼───────────────
           1 │     1  x
           2 │     2  y
           3 │     3  z"""

    io = IOBuffer()
    show(io, MIME("text/plain"), df, eltypes=false)
    str = String(take!(io))
    @test str == """
        3×2 DataFrame
         Row │ A  B
        ─────┼──────
           1 │ 1  x
           2 │ 2  y
           3 │ 3  z"""

    io = IOBuffer()
    show(io, MIME("text/plain"), eachcol(df), eltypes=false)
    str = String(take!(io))
    @test str == """
        3×2 DataFrameColumns
         Row │ A  B
        ─────┼──────
           1 │ 1  x
           2 │ 2  y
           3 │ 3  z"""

    io = IOBuffer()
    show(io, MIME("text/plain"), eachrow(df), eltypes=false)
    str = String(take!(io))
    @test str == """
        3×2 DataFrameRows
         Row │ A  B
        ─────┼──────
           1 │ 1  x
           2 │ 2  y
           3 │ 3  z"""

    io = IOBuffer()
    show(io, MIME("text/html"), df, eltypes=true)
    str = String(take!(io))
    @test str == "

3 rows × 2 columns
" *
                 "

" *
                 "	A	B
" *
                 "
	Int32	String

" *
                 "
1	1	x
" *
                 "
2	2	y

3	3	z

"

    io = IOBuffer()
    show(io, MIME("text/html"), eachcol(df), eltypes=true)
    str = String(take!(io))
    @test str == "
3×2 DataFrameColumns

" *
                 "	A	B
" *
                 "
	Int32	String

" *
                 "
1	1	x
" *
                 "
2	2	y

3	3	z

"

    io = IOBuffer()
    show(io, MIME("text/html"), eachrow(df), eltypes=true)
    str = String(take!(io))
    @test str == "
3×2 DataFrameRows

" *
                 "	A	B
" *
                 "
	Int32	String

" *
                 "
1	1	x
" *
                 "
2	2	y

3	3	z

"

    io = IOBuffer()
    show(io, MIME("text/html"), df, eltypes=false)
    str = String(take!(io))
    @test str == "

3 rows × 2 columns
" *
                 "

" *
                 "	A	B

" *
                 "
1	1	x
" *
                 "
2	2	y

3	3	z

"

    io = IOBuffer()
    show(io, MIME("text/html"), eachcol(df), eltypes=false)
    str = String(take!(io))
    @test str == "
3×2 DataFrameColumns

" *
                 "	A	B

" *
                 "
1	1	x
" *
                 "
2	2	y

3	3	z

"

    io = IOBuffer()
    show(io, MIME("text/html"), eachrow(df), eltypes=false)
    str = String(take!(io))
    @test str == "
3×2 DataFrameRows

" *
                 "	A	B

" *
                 "
1	1	x
" *
                 "
2	2	y

3	3	z

"

    for x in [df, eachcol(df), eachrow(df)]
        io = IOBuffer()
        show(io, MIME("text/latex"), x, eltypes=true)
        str = String(take!(io))
        @test str == """
        \\begin{tabular}{r|cc}
        \t& A & B\\\\
        \t\\hline
        \t& Int32 & String\\\\
        \t\\hline
        \t1 & 1 & x \\\\
        \t2 & 2 & y \\\\
        \t3 & 3 & z \\\\
        \\end{tabular}
        """

        io = IOBuffer()
        show(io, MIME("text/latex"), x, eltypes=false)
        str = String(take!(io))
        @test str == """
        \\begin{tabular}{r|cc}
        \t& A & B\\\\
        \t\\hline
        \t1 & 1 & x \\\\
        \t2 & 2 & y \\\\
        \t3 & 3 & z \\\\
        \\end{tabular}
        """
    end
end

@testset "improved printing of special types" begin
    df = DataFrame(A=Int64.(1:9), B = Vector{Any}(undef, 9))
    df.B[1:8] = [df, # DataFrame
                 df[1,:], # DataFrameRow
                 view(df,1:1, :), # SubDataFrame
                 eachrow(df), # DataFrameColumns
                 eachcol(df), # DataFrameRows
                 groupby(df, :A), # GroupedDataFrame
                 missing,
                 nothing]

    io = IOBuffer()
    show(io, df)
    str = String(take!(io))

    @test str == """
    9×2 DataFrame
     Row │ A      B
         │ Int64  Any
    ─────┼──────────────────────────────────────────
       1 │     1  9×2 DataFrame
       2 │     2  2-element DataFrameRow
       3 │     3  1×2 SubDataFrame
       4 │     4  9-element DataFrameRows
       5 │     5  2-element DataFrameColumns
       6 │     6  GroupedDataFrame with 9 groups b…
       7 │     7  missing
       8 │     8
       9 │     9  #undef"""

    # TODO: update when https://github.com/KristofferC/Crayons.jl/issues/47 is resolved
    if VERSION >= v"1.6" && Base.get_have_color()
        io = IOBuffer()
        show(IOContext(io, :color => true), df)
        str = String(take!(io))
        @test str == """
        \e[1m9×2 DataFrame\e[0m
        \e[1m Row \e[0m│\e[1m A     \e[0m\e[1m B                                 \e[0m
        \e[1m     \e[0m│\e[90m Int64 \e[0m\e[90m Any                               \e[0m
        ─────┼──────────────────────────────────────────
           1 │     1 \e[90m 9×2 DataFrame                     \e[0m
           2 │     2 \e[90m 2-element DataFrameRow            \e[0m
           3 │     3 \e[90m 1×2 SubDataFrame                  \e[0m
           4 │     4 \e[90m 9-element DataFrameRows           \e[0m
           5 │     5 \e[90m 2-element DataFrameColumns        \e[0m
           6 │     6 \e[90m GroupedDataFrame with 9 groups b… \e[0m
           7 │     7 \e[90m missing                           \e[0m
           8 │     8 \e[90m                                   \e[0m
           9 │     9 \e[90m #undef                            \e[0m"""
    end

    io = IOBuffer()
    show(io, MIME("text/html"), df)
    str = String(take!(io))
    @test str == "

9 rows × 2 columns
" *
                 "

	A	B
" *
                 "
	Int64	Any

" *
                 "
" *
                 "
1	1	9×2 DataFrame
" *
                 "
2	2	2-element DataFrameRow
" *
                 "
3	3	1×2 SubDataFrame
" *
                 "
4	4	9-element DataFrameRows
" *
                 "
5	5	2-element DataFrameColumns
" *
                 "
6	6	GroupedDataFrame with 9 groups based on key: A
" *
                 "
7	7	missing
" *
                 "
8	8
" *
                 "
9	9	#undef

"

    io = IOBuffer()
    show(io, MIME("text/latex"), df)
    str = String(take!(io))
    @test str == """
    \\begin{tabular}{r|cc}
    \t& A & B\\\\
    \t\\hline
    \t& Int64 & Any\\\\
    \t\\hline
    \t1 & 1 & \\emph{9×2 DataFrame} \\\\
    \t2 & 2 & \\emph{2-element DataFrameRow} \\\\
    \t3 & 3 & \\emph{1×2 SubDataFrame} \\\\
    \t4 & 4 & \\emph{9-element DataFrameRows} \\\\
    \t5 & 5 & \\emph{2-element DataFrameColumns} \\\\
    \t6 & 6 & \\emph{GroupedDataFrame with 9 groups based on key: A} \\\\
    \t7 & 7 & \\emph{missing} \\\\
    \t8 & 8 &  \\\\
    \t9 & 9 & \\emph{\\#undef} \\\\
    \\end{tabular}
    """

    @test_throws UndefRefError show(io, MIME("text/csv"), df)
    @test_throws UndefRefError show(io, MIME("text/tab-separated-values"), df)

    df[end, 2] = "\""
    push!(df, (10, Symbol("\"")))
    push!(df, (11, '"'))
    df.B[6] = groupby(df, :A)
    io = IOBuffer()
    show(io, MIME("text/csv"), df)
    str = String(take!(io))
    @test str == """
    "A","B"
    1,"11×2 DataFrame"
    2,"2-element DataFrameRow"
    3,"1×2 SubDataFrame"
    4,"11-element DataFrameRows"
    5,"2-element DataFrameColumns"
    6,"GroupedDataFrame with 11 groups based on key: A"
    7,missing
    8,nothing
    9,"\\""
    10,"\\""
    11,"\\""
    """

    io = IOBuffer()
    show(io, MIME("text/tab-separated-values"), df)
    str = String(take!(io))
    @test str == """
    "A"\t"B"
    1\t"11×2 DataFrame"
    2\t"2-element DataFrameRow"
    3\t"1×2 SubDataFrame"
    4\t"11-element DataFrameRows"
    5\t"2-element DataFrameColumns"
    6\t"GroupedDataFrame with 11 groups based on key: A"
    7\tmissing
    8\tnothing
    9\t"\\""
    10\t"\\""
    11\t"\\""
    """
end

@testset "check truncate keyword argument" begin
    df = DataFrame(x = "0123456789"^10)

    # no truncation
    io = IOBuffer()
    show(io, MIME("text/html"), df)
    str = String(take!(io))
    @test str == "

1 rows × 1 columns
" *
                 "

	x
" *
                 "
	String

" *
                 "

1	" *
                 "01234567890123456789012345678901234567890123456789" *
                 "01234567890123456789012345678901234567890123456789	"*
                 "

"

    # no truncation
    io = IOBuffer()
    show(io, MIME("text/latex"), df)
    str = String(take!(io))
    @test str == """
    \\begin{tabular}{r|c}
    \t& x\\\\
    \t\\hline
    \t& String\\\\
    \t\\hline
    \t1 & 0123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789 \\\\
    \\end{tabular}
    """

    # no truncation
    io = IOBuffer()
    show(io, MIME("text/csv"), df)
    str = String(take!(io))
    @test str == "\"x\"\n\"01234567890123456789012345678901234567890123456789" *
                 "01234567890123456789012345678901234567890123456789\"\n"

    # no truncation
    io = IOBuffer()
    show(io, MIME("text/tab-separated-values"), df)
    str = String(take!(io))
    @test str == "\"x\"\n\"01234567890123456789012345678901234567890123456789" *
                 "01234567890123456789012345678901234567890123456789\"\n"

    # default truncation
    io = IOBuffer()
    show(io, MIME("text/plain"), df)
    str = String(take!(io))
    @test str == """
        1×1 DataFrame
         Row │ x
             │ String
        ─────┼───────────────────────────────────
           1 │ 01234567890123456789012345678901…"""

    io = IOBuffer()
    show(io, df)
    str = String(take!(io))
    @test str == """
        1×1 DataFrame
         Row │ x
             │ String
        ─────┼───────────────────────────────────
           1 │ 01234567890123456789012345678901…"""

    # no truncation
    io = IOBuffer()
    show(io, df, truncate=0)
    str = String(take!(io))
    @test str == """
        1×1 DataFrame
         Row │ x
             │ String
        ─────┼──────────────────────────────────────────────────────────────────────────────────────────────────────
           1 │ 0123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789"""

    # custom truncation
    io = IOBuffer()
    show(io, df, truncate=1)
    str = String(take!(io))
    @test str == """
        1×1 DataFrame
         Row │ x
             │ String
        ─────┼────────
           1 │ 01234…"""

    df = DataFrame(x12345678901234567890 = "0123456789"^10)
    io = IOBuffer()
    show(io, df, truncate=1, rowlabel=:r12345678901234567890)
    str = String(take!(io))
    @test str == """
        1×1 DataFrame
         r12345678901234567890 │ x12345678901234567890
                               │ String
        ───────────────────────┼───────────────────────
                             1 │ 01234567890123456789…"""

end

end # module


text:
@testset "AlgAssRelOrd" begin
  include("AlgAssRelOrd/Order.jl")
  include("AlgAssRelOrd/Ideal.jl")
end


text:
using Test
using DynamicPolynomials
# using Revise
using PolynomialBasis
PB = PolynomialBasis

function allapprox(v1,v2)
    return all(v1 .≈ v2)
end

basis = PB.HermiteTensorProductBasis(2)

testvals = zeros(16)
testvals[1] = 1.0
@test allapprox(basis(-1.,-1.),testvals)
@test allapprox(basis([-1.,-1.]),testvals)

fill!(testvals,0.0)
testvals[5] = 1.0
@test allapprox(basis(-1.,1.),testvals)
@test allapprox(basis([-1.,1.]),testvals)

fill!(testvals,0.0)
testvals[9] = 1.0
@test allapprox(basis(1,-1),testvals)
@test allapprox(basis([1,-1]),testvals)

fill!(testvals,0.0)
testvals[13] = 1.0
@test allapprox(basis(1,1),testvals)
@test allapprox(basis([1,1]),testvals)

testvals = zeros(16,2)
testvals[2,2] = 1.0
testvals[3,1] = 1.0
@test allapprox(gradient(basis,-1,-1),testvals)
@test allapprox(gradient(basis,[-1,-1]),testvals)

fill!(testvals,0.0)
testvals[6,2] = 1.0
testvals[7,1] = 1.0
@test allapprox(gradient(basis,-1,1),testvals)
@test allapprox(gradient(basis,[-1,1]),testvals)

fill!(testvals,0.0)
testvals[10,2] = 1.0
testvals[11,1] = 1.0
@test allapprox(gradient(basis,1,-1),testvals)
@test allapprox(gradient(basis,[1,-1]),testvals)

fill!(testvals,0.0)
testvals[14,2] = 1.0
testvals[15,1] = 1.0
@test allapprox(gradient(basis,1,1),testvals)
@test allapprox(gradient(basis,[1,1]),testvals)


text:
@reexport module Utils

using Random: randperm
export mergedict, getiris
export skipmean,skipmedian,skipstd
export aggregatorclskipmissing

export holdout,
       kfold,
       score,
       infer_eltype,
       nested_dict_to_tuples,
       nested_dict_set!,
       nested_dict_merge,
       create_transformer

using TSML: Transformer
using Statistics
using DataFrames
using CSV

import MLBase: Kfold

"""
    holdout(n, right_prop)

Holdout method that partitions a collection
into two partitions.

- `n`: Size of collection to partition
- `right_prop`: Percentage of collection placed in right partition

Returns: two partitions of indices, left and right
"""
function holdout(n, right_prop)
  shuffled_indices = randperm(n)
  partition_pivot = round(Int,right_prop * n)
  right = shuffled_indices[1:partition_pivot]
  left = shuffled_indices[partition_pivot+1:end]
  return (left, right)
end

"""
    kfold(num_instances, num_partitions)

Returns k-fold partitions.

- `num_instances`: total number of instances
- `num_partitions`: number of partitions required

Returns: training set partition.
"""
function kfold(num_instances, num_partitions)
  return collect(Kfold(num_instances, num_partitions))
end

"""
    score(metric::Symbol, actual, predicted)

Score learner predictions against ground truth values.

Available metrics:
- :accuracy

- `metric`: metric to assess with
- `actual`: ground truth values
- `predicted`: predicted values

Returns: score of learner
"""
function score(metric::Symbol, actual, predicted)
  if metric == :accuracy
    mean(actual .== predicted) * 100.0
  else
    error("Metric $metric not implemented for score.")
  end
end

"""
    infer_eltype(vector::Vector)

Returns element type of vector unless it is Any.
If Any, returns the most specific type that can be
inferred from the vector elements.

- `vector`: vector to infer element type on

Returns: inferred element type
"""
function infer_eltype(vector::Vector)
  # Obtain element type of vector
  vec_eltype = eltype(vector)

  # If element type of Vector is Any and not empty,
  # and all elements are of the same type,
  # then return the vector elements' type.
  if vec_eltype == Any && !isempty(vector)
    all_elements_same_type = all(x -> typeof(x) == typeof(first(vector)), vector)
    if all_elements_same_type
      vec_eltype = typeof(first(vector))
    end
  end

  # Return inferred element type
  return vec_eltype
end

"""
    nested_dict_to_tuples(dict::Dict)

Converts nested dictionary to list of tuples

- `dict`: dictionary that can have other dictionaries as values

Returns: list where elements are ([outer-key, inner-key, ...], value)
"""
function nested_dict_to_tuples(dict::Dict)
  set = Set()
  for (entry_id, entry_val) in dict
    if typeof(entry_val) <: Dict
      inner_set = nested_dict_to_tuples(entry_val)
      for (inner_entry_id, inner_entry_val) in inner_set
        new_entry = (vcat([entry_id], inner_entry_id), inner_entry_val)
        push!(set, new_entry)
      end
    else
      new_entry = ([entry_id], entry_val)
      push!(set, new_entry)
    end
  end
  return set
end

"""
    nested_dict_set!(dict::Dict, keys::Array{T, 1}, value) where {T}

Set value in a nested dictionary.

- `dict`: nested dictionary to assign value
- `keys`: keys to access nested dictionaries in sequence
- `value`: value to assign
"""
function nested_dict_set!(dict::Dict, keys::Array{T, 1}, value) where {T}
  inner_dict = dict
  for key in keys[1:end-1]
    inner_dict = inner_dict[key]
  end
  inner_dict[keys[end]] = value
end

"""
    nested_dict_merge(first::Dict, second::Dict)

Second nested dictionary is merged into first.

If a second dictionary's value as well as the first
are both dictionaries, then a merge is conducted between
the two inner dictionaries.
Otherwise the second's value overrides the first.

- `first`: first nested dictionary
- `second`: second nested dictionary

Returns: merged nested dictionary
"""
function nested_dict_merge(first::Dict, second::Dict)
  target = copy(first)
  for (second_key, second_value) in second
    values_both_dict =
      typeof(second_value) <: Dict &&
      typeof(get(target, second_key, nothing)) <: Dict
    if values_both_dict
      target[second_key] = nested_dict_merge(target[second_key], second_value)
    else
      target[second_key] = second_value
    end
  end
  return target
end

"""
    create_transformer(prototype::Transformer, options=nothing)

Create transformer

- `prototype`: prototype transformer to base new transformer on
- `options`: additional options to override prototype's options

Returns: new transformer
"""
function create_transformer(prototype::Transformer, options=nothing)
  new_options = copy(prototype.options)
  if options != nothing
    new_options = nested_dict_merge(new_options, options)
  end

  prototype_type = typeof(prototype)
  return prototype_type(new_options)
end

"""
    aggregatorclskipmissing(fn::Function)

Function to create aggregator closure with skipmissing features
"""
function aggregatorclskipmissing(fn::Function)
  function skipagg(x::Union{AbstractArray,DataFrame})
    if length(collect(skipmissing(x))) == 0
      return missing
    else
      return fn(skipmissing(x))
    end
  end
  return skipagg
end

function skipmean(x::T) where {T<:Union{AbstractArray,DataFrame}}
  if length(collect(skipmissing(x))) == 0
    missing
  else
    mean(skipmissing(x))
  end
end

function skipmedian(x::T) where {T<:Union{AbstractArray,DataFrame}}
  if length(collect(skipmissing(x))) == 0
    missing
  else
    median(skipmissing(x))
  end
end

function skipstd(x::T) where {T<:Union{AbstractArray,DataFrame}}
  if length(collect(skipmissing(x))) == 0
    missing
  else
    std(skipmissing(x))
  end
end

"""
    mergedict(first::Dict, second::Dict)

Second nested dictionary is merged into first.

If a second dictionary's value as well as the first
are both dictionaries, then a merge is conducted between
the two inner dictionaries.
Otherwise the second's value overrides the first.

- `first`: first nested dictionary
- `second`: second nested dictionary

Returns: merged nested dictionary
"""
function mergedict(first::Dict, second::Dict)
  target = copy(first)
  for (second_key, second_value) in second
    values_both_dict =
      typeof(second_value) <: Dict &&
      typeof(get(target, second_key, nothing)) <: Dict
    if values_both_dict
      target[second_key] = mergedict(target[second_key], second_value)
    else
      target[second_key] = second_value
    end
  end
  return target
end

function getiris()
  iris = CSV.read(joinpath(Base.@__DIR__,"../data","iris.csv"))
  return iris
end

end


text:
function backward{N}(backend::GPUBackend, state::SplitLayerState{N}, inputs::Vector{Blob}, diffs::Vector{Blob})
  if !isa(diffs[1], NullBlob)
    diff = diffs[1]
    len = length(diff)
    one = convert(eltype(diff), 1)
    for i = 2:N
      CuBLAS.axpy(backend.cublas_ctx, len, one, state.blobs_diff[i].ptr, 1, diff.ptr, 1)
    end
  end
end


text:
using lipm
using Test

@testset "lipm.jl" begin
    # Write your tests here.
end


text:
function get_sizes(obj::Objective)
    n = size(obj.Q,1)
    m = size(obj.R,1)
    return n,m
end

function get_sizes(X::Matrix,U::Matrix)
    n,N = size(X)
    m = size(U,1)
    return n,m,N
end

function convertInf!(A::Matrix,infbnd=1.1e20)
    infs = isinf.(A)
    A[infs] = sign.(A[infs])*infbnd
    return nothing
end

function get_weights(method,N::Int)
    if method == :trapezoid
        weights = ones(N)
        weights[[1,end]] = 0.5
    elseif method == :hermite_simpson_separated
        weights = ones(N)*2/3
        weights[2:2:end] = 4/3
        weights[[1,end]] = 1/3
    end
    weights *= dt
    return weights
end

function get_bounds(obj::Objective,N::Int)
    n,m = get_sizes(obj)
    lb = zeros((n+m),N)
    ub = zeros((n+m),N)

    lb[1:n,:] .= obj.x_min
    ub[1:n,:] .= obj.x_max
    lb[n+(1:m),:] .= obj.u_min
    ub[n+(1:m),:] .= obj.u_max

    # Initial Condition
    lb[1:n,1] .= obj.x0
    ub[1:n,1] .= obj.x0

    # Terminal Constraint
    lb[1:n,N] .= obj.xf
    ub[1:n,N] .= obj.xf

    # Convert Infinite bounds
    convertInf!(lb)
    convertInf!(ub)
    return vec(lb), vec(ub)
end

"""
Stack state and controls for all time steps into a single vector of variables
Z = [X1,U1,X2,U2,..,] (all transposed)
"""
function packZ(X,U)
    n, N = size(X)
    m = size(U,1)
    Z = zeros(n+m,N)
    Z[1:n,:] .= X
    Z[n+1:end,1:end] .= U
    Z = vec(Z)
end

function unpackZ(Z, sze)
    n,m,N = sze
    Z = reshape(Z,n+m,N)
    X = Z[1:n,:]
    U = Z[n+1:end,:]
    return X,U
end

function get_pack(Z, obj)
    n,m = get_sizes(obj)
    return (n,m, length(Z)//(n*m))
end

function get_initial_state(obj::Objective, N::Int)
    n,m = get_sizes(obj)
    X0 = line_trajectory(obj.x0, obj.xf, N)
    U0 = zeros(m,N)
    Z0 = packZ(X0,U0)
end

"""
Evaluate Objective Value
"""
function cost(X,U,weights,obj::Objective)
    N = size(X,2)
    Qf = obj.Qf; xf = obj.xf;
    J = zeros(eltype(X),N)
    for k = 1:N
        J[k] = cost_k(X[:,k],U[:,k],obj)
    end
    J = weights'J
    J += 0.5*(X[:,N] - xf)'*Qf*(X[:,N] - xf)
end

function cost_k(x,u,obj::Objective)
    xf = obj.xf
    0.5*(x-xf)'obj.Q*(x-xf) + 0.5*u'obj.R*u
end

"""
Gradient of Objective
"""
function cost_gradient(X,U,weights,obj::Objective)
    Q = obj.Q; xf = obj.xf; R = obj.R; Qf = obj.Qf;
    n,m,N = get_sizes(X,U)
    grad_f = zeros(eltype(X),(n+m)*N)
    grad_f = reshape(grad_f, n+m, N)
    for k = 1:N
        grad_f[1:n,k] = weights[k]*Q*(X[:,k] - xf) # weight includes dt
        grad_f[n+1:end,k] = weights[k]*R*U[:,k]
    end
    grad_f[1:n,N] += Qf*(X[:,N] - xf)
    return vec(grad_f)
end

"""
Evaluate constraint values
[ dt*(f(x1)+f(x2))/2 - x2 + x1     ] Collocation Constraints
[ dt*(f(x2)+f(x3))/2 - x3 + x2     ]
[           ...                 ]
[ dt*(f(xN-1)+f(xN))/2 - xN + xN-1 ]
"""
function calc_constraints(X,U,method,dt,f::Function)
    n,m,N = get_sizes(X,U)
    g = zeros(eltype(X),(N-1)*n)
    g = reshape(g,n,N-1)
    if method == :trapezoid
        for k = 1:N-1
            # Collocation Constraints
            g[:,k] = dt*( f(X[:,k+1],U[:,k+1]) + f(X[:,k],U[:,k]) )/2 - X[:,k+1] + X[:,k]
        end
    elseif method == :hermite_simpson_separated
        fVal = zeros(eltype(X),n,N)
        for k = 1:N
            fVal[:,k] = f(X[:,k],U[:,k])
        end
        iLow = 1:2:N-1
        iMid = iLow + 1
        iUpp = iMid + 1

        midpoints =  - X[:,iMid] + (X[:,iLow]+X[:,iUpp])/2 + dt*(fVal[:,iLow] - fVal[:,iUpp])/8
        collocation = - X[:,iUpp] + X[:,iLow] + dt*(fVal[:,iLow] + 4*fVal[:,iMid] + fVal[:,iUpp])/6
        g[:,iLow] = collocation
        g[:,iMid] = midpoints
    end
    # reshape(g,n*(N-1),1)
    return vec(g)
end

"""
Constraint Jacobian
"""
function constraint_jacobian(X,U,dt,method,F::Function)
    n,m,N = get_sizes(X,U)
    Z = packZ(X,U)
    jacob_g = zeros((N-1)*n,N*(n+m))
    z = reshape(Z,n+m,N)

    if method == :trapezoid
        nm = n+m
        # First time step
        fz = F(z[:,1])
        jacob_g[1:n,1:n+m] .= dt*fz/2+eye(n,n+m)

        # Loop over time steps
        for k = 2:N-1
            off_1 = (k-1)*(n)
            off_2 = (k-1)*(n+m)
            # Calculate (n,n+m) Jacobian of both states and controls
            fz = F(z[:,k])
            jacob_g[off_1-n+(1:n),off_2+(1:n+m)] .= dt*fz/2 - eye(n,n+m)
            jacob_g[off_1 + (1:n),off_2+(1:n+m)] .= dt*fz/2 + eye(n,n+m)
        end

        # Last time step
        fz = F(z[:,N])
        jacob_g[end-n+1:end,end-n-m+1:end] = dt*fz/2-eye(n,n+m)

    elseif method == :hermite_simpson_separated
        nSeg = Int((N-1)/2)
        jacob_g = zeros(n*(N-1),(n+m)N)
        fz1 = F(z[:,1])
        Inm = Matrix(I,n,n+m)

        function calc_block(k)
            vals = zeros(2n,3(n+m))
            fz = F(z[:,k])
            vals[0+(1:n),(1:n+m)] .= dt*fz/6 + Inm
            vals[n+(1:n),(1:n+m)] .= dt*fz/8 + Inm/2
            fm = F(z[:,k+1])
            vals[0+(1:n),n+m+(1:n+m)] .= 2*dt*fm/3
            vals[n+(1:n),n+m+(1:n+m)] .= -Inm
            fz1 .= F(z[:,k+2])
            vals[0+(1:n),2(n+m)+(1:n+m)] .=  dt*fz1/6 - Inm
            vals[n+(1:n),2(n+m)+(1:n+m)] .= -dt*fz1/8 + Inm/2
            return vals
        end

        for i = 1:nSeg
            off_1 = 2(i-1)*(n)
            off_2 = 2(i-1)*(n+m)
            k = 2i-1

            jacob_g[off_1+(1:2n), off_2+(1:3(n+m))] = calc_block(k)
        end
    end

    return jacob_g
end


text:
#=
This file was used for calculating the offset of Sirius B from Sirius A. We
needed to do this because it is very hard to lock onto Sirius B directly,
in addition to the fact we were using Sirius A as the AO NGS and targeting
Sirius B in an off-axis manner.
=#

using AstroAngles
using AstroLib
using Dates
using Measurements
using SkyCoords
using Unitful, UnitfulAstro

# http://simbad.u-strasbg.fr/simbad/sim-basic?Ident=sirius+a&submit=SIMBAD+search
sirius = ICRSCoords(hms"06:45:08.91728", dms"-16:42:48.0171")
parallax = 379.21e-3 #mas
distance = 1 / parallax # pc

## orbital parameters
## https://www.stelledoppie.it/index2.php?iddoppia=27936
# a = deg2rad((7.4957 ± 0.0025) / distance / 3600) # semi-major axis
a = (7.4957 ± 0.0025)
e = (0.59142 ± 0.00037)          # eccentricity
Ω = deg2rad(45.4 ± 0.071)        # logitude of ascending node
i = deg2rad(136.336 ± 0.04)      # inclination
ω = deg2rad(149.161 ± 0.075)     # argument of periapsis
T = (1994.5715 ± 0.0058)u"yr" # periastron epoch
P = (50.1284 ± 0.0043)u"yr"   # orbital period

## Solve Kepler's equations
# mean anomaly
M(t) = ustrip(u"NoUnits", 2π * (t - T) / P)
# eccentric anomaly
E(t) = kepler_solver(M(t), e)

## equation of a conic section
# true anomaly
ν(t) = 2atan(sqrt((1 + e) / (1 - e)) * tan(E(t)/2))
r(t) = a * (1 - e^2) / (1 + e * cos(ν(t)))

pa(t) = Ω + atan(tan(ν(t) + ω) * cos(i))
sep(t) = r(t) * cos(ν(t) + ω) / cos(pa(t) - Ω)

year_val(epoch) = year_val(Date(epoch))
year_val(epoch::TimeType) = year(epoch) + (dayofyear(epoch) - 1) / daysinyear(epoch)

for epoch in ["2020-02-04", "2020-11-21", "2020-11-28"]
    t = year_val(epoch)u"yr"
    ρ = sep(t)
    θ = pa(t) |> rad2deg

    @info epoch=epoch t ρ θ
end

# t = 2020.8876712328767123u"yr" # time for 11/20/2020
# t = 2020.9094227697923797u"yr" # time for 11/27/2020
# ρ_arcsec = sep(t)
# θ = pa(t)

# ρ = deg2rad(ρ_arcsec / 3600)
# θ_deg = rad2deg(θ)

# println("ρ= $ρ_arcsec\", θ= $(θ_deg)°")
# println("offset: $(ρ_arcsec .* sincosd(θ_deg))")

# sirius_b = offset(sirius, Measurements.value(ρ), Measurements.value(θ))

# dra = sirius_b.ra - sirius.ra
# # ra = sirius_b.ra
# ra_hrs = 24 * dra / 2π
# ra_h, ra_r = divrem(sirius_b.ra * 24 / 2π, 1)
# ra_m, ra_r = divrem(ra_r * 60, 1)
# ra_s = ra_r * 60

# ddec = sirius_b.dec - sirius.dec
# # dec = sirius_b.dec
# dec_degs = rad2deg(ddec)
# dec_d, dec_r = divrem(rad2deg(sirius_b.dec), 1)
# dec_m, dec_r = divrem(dec_r * 60, 1)
# dec_s = dec_r * 60

# println("ρ= $ρ_arcsec\", θ= $(θ_deg)°")
# println("sirius b = $(ra_h)h$(ra_m)m$(ra_s)s $(dec_d):$(dec_m):$(dec_s)")

# println("Δ(ra, dec)= $(ra_hrs * 24 * 3600 / 2π )s $(dec_degs * 3600)\"")
# println("Δ(ra, dec)= $(rad2deg(ra_hrs) * 3600)\" $(dec_degs * 3600)\"")

max_stars_repo_name: mileslucas/sirius-b

max_issues_repo_name: mileslucas/sirius-b

max_forks_repo_name: mileslucas/sirius-b


text:
using Memcache
using Test

@test isa(Memcache, Module)

include("test.jl")


text:
module HTMLElements
  # copied from https://discourse.julialang.org/t/best-tool-for-printing-tables-in-jupyter-notebook/3548/8
  export table

  function table(cur_matrix::Matrix; header_row=[], title=nothing)

    cur_table = ""

    if title != nothing
      cur_table *= "
$title
"
    end

    cur_table *= "
"

    if !isempty(header_row)
      cur_table *= "

"

      for cur_header in header_row
        cur_table *= "$cur_header	"
      end

      cur_table *= "

"
    end

    cur_table *= "
"

    for ii in 1:size(cur_matrix, 1)

      cur_table *= "
"

      for jj in 1:size(cur_matrix, 2)
        cur_table *= ""

        cur_table *= string(cur_matrix[ii, jj])

        cur_table *= "	"
      end

      cur_table *= "
"

    end

    cur_table *= "
"

    cur_table *= "
"

    return HTML(cur_table)

  end

end

max_stars_repo_name: simulkade/peteng

max_issues_repo_name: simulkade/peteng

max_forks_repo_name: simulkade/peteng


text:
module RandomMatrices

using Random

export unitvector, orthogonal, normal, general, symmetric, hermitian

using Random: GLOBAL_RNG
using LinearAlgebra:  reflector!, QR

import Base.randn, Base.rand

"""
unitvector(m, T::Type [, rng::AbstractRNG]) -> Vector{T}

  create random unit vector, uniformly distributed on unit sphere of `m` dimensions.
  `T` may be any real or complex float type.
"""
function unitvector(m::Integer, S::Union{Type{T},Type{Complex{T}}}=Float64, rng::AbstractRNG = GLOBAL_RNG) where T<:AbstractFloat

  unitvector!(zeros(S, m), rng)
end

"""
  unitvector!(v::AbstractVector{T} [, rng::AbstractRNG]) -> v

  Same as unitvector, populating the given vector or 1-dim view `v`.
"""
function unitvector!(v::Union{AbstractVector{T},AbstractVector{Complex{T}}}, rng::AbstractRNG = GLOBAL_RNG) where T<:AbstractFloat

  n = length(v)
  r = zero(T)
  sqnh = sqrt(n) / 2
  while r < sqnh
    randn!(rng, v)
    r = norm(v)
  end
  normalize!(v)
end

"""
  orthogonal(m, n, ::Type{T}=Float64 [, rng::AbstractRNG]) -> LinearAlgebra.QRCompactWYQ

Generate random orthogonal real or unitary complex matrix.
"""
function orthogonal(m::Integer, n::Integer, S::Type=Float64, rng::AbstractRNG=GLOBAL_RNG)
  S = S == Complex ? Complex{Float64} : S
  S <: Union{T,Complex{T}} where T <: AbstractFloat || error("illegal type $S")
  m >= n || error("more independent vectors($n) than dimension($m) requested")
  A = zeros(S, m, n)
  τ = zeros(S, min(m,n))

  n2 = n == m && (S <: Real) && rand(Bool) ? n - 1 : n
  for k = 1:n2
    x = view(A, k:m, k)
    unitvector!(x)
    τk = reflector!(x)
    τ[k] = τk
  end
  QR(A, τ)[:Q]
end

# modify rand and randn to produce BigFloat and complex values
function rand(rng::AbstractRNG, ::Type{BigFloat})
  bigFloat_randfill(rng, rand(rng, Float64))
end
function randn(rng::AbstractRNG, ::Type{BigFloat})
  bigFloat_randfill(rng, randn(rng, Float64))
end
function rand(rng::AbstractRNG, ::Type{Complex{T}}) where T<:AbstractFloat
  a, b = 2rand(T, 2) - 1
  while hypot(a,b) > 1
    a, b = 2rand(T, 2) - 1
  end
  complex(a, b)
end

"""
Append random bits in BigFloat representation of input small float number.
"""
function bigFloat_randfill(rng::AbstractRNG, a::T) where {T<:Union{Float64,Float32,Float16}}
  missing = precision(BigFloat) - precision(a)
  bigi = BigInt(1) << missing
  add = ( rand(rng, BigInt(0):bigi-1) / bigi )
  add *= big"2.0"^(exponent(a) - precision(a))
  add = copysign(add, a)
  a + add
end

"""
  diagonal(dia::AbstractVector{T}) -> Array{T,2}

Create quasi diagonal real matrix with diagonal blocks given by `dia`.
`dia` must contain only real and pairs of conjugate complex numbers.
If complex output matrix is required, use `diagm(dia)`.
"""
function diagonal(dia::Union{AbstractVector{T}, AbstractVector{Complex{T}}}) where T<:AbstractFloat

  m = length(dia)
  mr = m
  j = 1
  while j <= m
    a = dia[j]
    if !isreal(a)
      if j >= m || imag(a) != -imag(dia[j+1])
        throw(InexactError())
      else
        j += 1
      end
    end
    j += 1
  end

  A = zeros(T, mr, mr)

  k = 1
  while k <= m
    dr = real(dia[k])
    di = imag(dia[k])
    if di == 0
      A[k,k] = dr
      k += 1
    else
      A[k,k] = A[k+1,k+1] = dr
      A[k+1,k] = di
      A[k,k+1] = -di
      k += 2
    end
  end
  A
end

"""
  normal(dia::AbstractVector, [,rng]) -> Array{T,2}

Create real or complex normal matrix with eigenvalues given by dia.
"""
function normal(dia::Union{AbstractVector{T}, AbstractVector{Complex{T}}}, S::Union{Type{T},Type{Complex{T}}}=Float64, rng::AbstractRNG = GLOBAL_RNG) where T<:AbstractFloat

  D = S<:Real ? diagonal(dia) : diagm(dia)
  m = size(D, 1)
  Q = orthogonal(m, m, S, rng)
  Q * D * Q'
end

"""
  general(dia::Vector{<:AbstractFloat}, m, n, S::Type=Float64, [,rng]) -> Array{S,2}

Create general real or complex (m,n)-matrix with singular values given.

"""
function general(singular_values::AbstractVector{T}, m::Int, n::Int, S::Union{Type{T},Type{Complex{T}}}=Float64, rng::AbstractRNG=GLOBAL_RNG) where T<:AbstractFloat

  j = length(singular_values)
  j <= min(m, n) || error("too many singular values")
  all(singular_values .> 0) || error("singular values must be positive")

  D = zeros(S, m, n)
  for k = 1:j
    D[k,k] = singular_values[k]
  end
  orthogonal(m, m, S, rng) * D * orthogonal(n, n, S, rng)'
end

end # module

max_stars_repo_name: KlausC/LinAlgBigFloat

max_issues_repo_name: KlausC/LinAlgBigFloat

max_forks_repo_name: KlausC/LinAlgBigFloat


text:
# Types and constructors

struct MonomialBasis{P,V,G,D} <: TensorPolynomialBasis{P,V,G}
  terms::Vector{CartesianIndex{D}}
  orders::NTuple{D,Int}
end

struct MonomialBasisCache{T}
  c::Matrix{T}
  g::Matrix{T}
end

function (::Type{MonomialBasis{P,V}})(filter::Function,order::Int) where {P,V}
  D = _length(P)
  terms = _define_terms(filter,order,D)
  G = _gradient_type(V,P)
  orders = tuple(fill(order,D)...)
  MonomialBasis{P,V,G,D}(terms,orders)
end

function (::Type{MonomialBasis{P,V}})(orders::NTuple{N,Int}) where {P,V,N}
  terms = [ ci for ci in CartesianIndices(orders.+1) if true]
  D = _length(P)
  @assert D == N
  G = _gradient_type(V,P)
  MonomialBasis{P,V,G,D}(terms,orders)
end

# Implementation of the interface

length(b::MonomialBasis{P,V}) where {P,V} = _length(V)*length(b.terms)

ndims(b::MonomialBasis{P}) where P = _length(P)

function ScratchData(b::MonomialBasis{P,V}) where {P,V}
  T = eltype(V)
  dim = _length(P)
  n1d = maximum(b.orders.+1)
  c = zeros(dim,n1d)
  g = zeros(dim,n1d)
  MonomialBasisCache{T}(c,g)
end

function evaluate!(
  v::AbstractVector{V},
  b::MonomialBasis{P,V},
  x::P,
  cache::MonomialBasisCache) where {P,V}
  _evaluate_nd!(v,x,b.orders,b.terms,cache.c)
end

function gradient!(
  v::AbstractVector{G},
  b::MonomialBasis{P,V,G},
  x::P,
  cache::MonomialBasisCache) where {P,V,G}
  _gradient_nd!(v,x,b.orders,b.terms,cache.c,cache.g,V)
end

# Helpers

_q_filter(e,o) = true

function _define_terms(filter,order,dim)
  n1d = order+1
  t = tuple(fill(n1d,dim)...)
  g = tuple(fill(1,dim)...)
  cis = CartesianIndices(t)
  co = CartesianIndex(g)
  [ ci for ci in cis if filter(Tuple(ci-co),order) ]
end

function _evaluate_1d!(v::AbstractMatrix{T},x,order,d) where T
  n = order + 1
  z = one(T)
  @inbounds v[d,1] = z
  @inbounds for i in 2:n
    v[d,i] = x[d]^(i-1)
  end
end

function _gradient_1d!(v::AbstractMatrix{T},x,order,d) where T
  n = order + 1
  z = zero(T)
  @inbounds v[d,1] = z
  @inbounds for i in 2:n
    v[d,i] = (i-1)*x[d]^(i-2)
  end
end

function _evaluate_nd!(
  v::AbstractVector{V},
  x,
  orders,
  terms::AbstractVector{CartesianIndex{D}},
  c::AbstractMatrix{T}) where {V,T,D}

  dim = D
  for d in 1:dim
    _evaluate_1d!(c,x,orders[d],d)
  end

  o = one(T)
  k = 1

  for ci in terms

    s = o
    @inbounds for d in 1:dim
      s *= c[d,ci[d]]
    end

    k = _set_value!(v,s,k)

  end

end

@inline function _set_value!(v::AbstractVector{V},s::T,k) where {V,T}
  m = zero(_mutable(V))
  z = zero(T)
  js = eachindex(m)
  @inbounds for j in js
    @inbounds for i in js
      m[i] = z
    end
    m[j] = s
    v[k] = m
    k += 1
  end
  k
end

@inline function _set_value!(v::AbstractVector{<:Real},s,k)
    v[k] = s
    k+1
end

function _gradient_nd!(
  v::AbstractVector{G},
  x,
  orders,
  terms::AbstractVector{CartesianIndex{D}},
  c::AbstractMatrix{T},
  g::AbstractMatrix{T},
  ::Type{V}) where {G,T,D,V}

  dim = D
  for d in 1:dim
    _evaluate_1d!(c,x,orders[d],d)
    _gradient_1d!(g,x,orders[d],d)
  end

  z = zero(MVector{D,T})
  o = one(T)
  k = 1

  for ci in terms

    s = z
    for i in eachindex(s)
      s[i] = o
    end
    for q in 1:dim
      for d in 1:dim
        if d != q
          s[q] *= c[d,ci[d]]
        else
          s[q] *= g[d,ci[d]]
        end
      end
    end

    k = _set_gradient!(v,s,k,V)

  end

end

@inline function _set_gradient!(
  v::AbstractVector{G},s::MVector{D,T},k,::Type{<:Real}) where {G,D,T}

  v[k] = s
  k+1
end

@inline function _set_gradient!(
  v::AbstractVector{G},s::MVector{D,T},k,::Type{V}) where {V,G,D,T}

  m = zero(_mutable(G))
  w = zero(V)
  z = zero(T)
  for j in eachindex(w)
    @inbounds for i in eachindex(m)
      m[i] = z
    end
    for i in eachindex(s)
      m[i,j] = s[i]
    end
    v[k] = m
    k += 1
  end
  k
end


text:
using JuliaPackageExample
using Test

@testset "JuliaPackageExample.jl" begin
    # 2x + 3y
    @test my_f(2,1) == 7
    @test my_f(2,2) == 10

    # x^2 + y^2
    @test my_g(4,3) == 25
    @test my_g(2,3) == 13
end


text:
### A Pluto.jl notebook ###
# v0.14.8

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 6d960bb1-b4d7-49f3-8549-78786eabfc74
using PlutoUI, Plots, StatPlots, CardioModels

# ╔═╡ d494e27e-7c67-4abc-b491-0dd09e58f169
md"""
# Baselli Model
The model is defined through the following three equations, where $G$ denotes an all-zero and $M$ an all-pole transfer function:

$S_n = G_{S S}(z) S_n \quad + \quad G_{SI}(z)I_n \quad  + \quad G_{S\rho}(z) \rho_n  \quad  + \quad  M_{S}(z) w_{S}$
$I_n = G_{IS}(z)S_n \quad  + \hspace{40mm}  G_{I\rho}(z) \rho_n \quad  + \quad  M_{I}(z) w_{I}$
$\rho_n = \hspace{104mm}M_{\rho}(z) w_{\rho}$

It can be identified from data as an ARXAR-model using the generalized least squares method.
"""

# ╔═╡ e7684c53-fc26-4658-9825-d622c442f394
md"""
Choose one out of four predefined simulations:

Simulation:
$(@bind sim NumberField(1:4, default = 1))
"""

# ╔═╡ e4febc9c-bed3-4a06-9982-49546ff06269
md"""

* $(sim == 1 ? md"**1: respiration enters on S (LF are generated in the S, S loop)**" : md"1: respiration enters on S (LF are generated in the S, S loop)")
* $(sim == 2 ? md"**2: respiration enters on I (LF are generated in the S, S loop)**" : md"2: respiration enters on I (LF are generated in the S, S loop)")
* $(sim == 3 ? md"**3: the LF are generated by S regulating feedback in the S, S loop (respiration enters both on S and on I)**" : md"3 : the LF are generated by S regulating feedback in the S, S loop (respiration enters both on S and on I)")
* $(sim == 4 ? md"**4: LF are generated externally from the closed-loop system (i.e., via a central command-respiration enters both on S and on I)**" : md"4: LF are generated externally from the closed-loop system (i.e., via a central command-respiration enters both on S and on I)")

"""

# ╔═╡ 8ad5bc71-dea2-43cc-b745-9c47ff6d7dff
model = getModel(sim);

# ╔═╡ 8c387e69-c2c6-4e93-897a-bdd47c65764d
plot(model, lab = "", dpi = 500, title = "Model #$(sim)")

# ╔═╡ b4c65d6d-72a8-470d-96c3-3d0c6bd810aa
md"""

N: $(@bind N NumberField(1:10000, default = 512))

Plot Spectral Denstity: $(@bind ispsd CheckBox())
"""

# ╔═╡ 8f948c31-0561-42e6-a0eb-48891126d9aa
S, I, ρ = predict(model,  N);

# ╔═╡ 0993d11c-c05f-48e7-9c1d-c95bba3d7251
begin
	if ispsd
		pS = psdplot(S, lab = "", title = "")
		pI = psdplot(I, lab = "")
		pρ = psdplot(ρ, lab = "")
	else
		pS = plot(S, lab = "", ylab = "S [mmHg]")
		pI = plot(I, lab = "", ylab = "I [s]")
		pρ = plot(ρ, lab = "", ylab = "\\rho [a.u.]", xlab =  "beat #")
	end
	plot(pS, pI, pρ, layout = (3,1), dpi = 500)
end

# ╔═╡ bfe88fac-5336-4d2b-9d03-c2313ccf542c
md"""
# Idendification of model from data

Identify the model from generated data with different model orders.
"""

# ╔═╡ f2963e4f-e032-4fe1-99a8-2fca1c4667c7
md"""
model order L: $(@bind L NumberField(1:20, default = 10))
"""

# ╔═╡ b33cbdd4-0af3-4c3f-833e-334bf22cc964
modelEstimate = fitBaselli(S, I, ρ, na = L, nb = L, nd = L);

# ╔═╡ 86c8b99c-4d2a-43b7-ad8a-bb639791d97c
begin
	plot(model, lab = "", color = :green)
	plot!(modelEstimate, lab = "", color = :red, dpi = 500, size = (1000, 500))
end

# ╔═╡ 86d849c1-5a6b-4d65-9ea1-cfc85cda5072
postprocess(modelEstimate, S, I)

# ╔═╡ Cell order:
# ╠═6d960bb1-b4d7-49f3-8549-78786eabfc74
# ╟─d494e27e-7c67-4abc-b491-0dd09e58f169
# ╟─e7684c53-fc26-4658-9825-d622c442f394
# ╟─e4febc9c-bed3-4a06-9982-49546ff06269
# ╠═8ad5bc71-dea2-43cc-b745-9c47ff6d7dff
# ╠═8c387e69-c2c6-4e93-897a-bdd47c65764d
# ╟─b4c65d6d-72a8-470d-96c3-3d0c6bd810aa
# ╠═8f948c31-0561-42e6-a0eb-48891126d9aa
# ╟─0993d11c-c05f-48e7-9c1d-c95bba3d7251
# ╟─bfe88fac-5336-4d2b-9d03-c2313ccf542c
# ╟─f2963e4f-e032-4fe1-99a8-2fca1c4667c7
# ╠═b33cbdd4-0af3-4c3f-833e-334bf22cc964
# ╠═86c8b99c-4d2a-43b7-ad8a-bb639791d97c
# ╠═86d849c1-5a6b-4d65-9ea1-cfc85cda5072


text:
isdefined(Base, :__precompile__) && __precompile__()

"""
Module for constructing self-contained power system objects.
"""
module PowerSystemExtensions

#################################################################################
# Exports

export CleanEnergyReserve
export InertiaReserve
export ThermalCleanEnergy

#################################################################################
# Imports

import InfrastructureSystems
import InfrastructureSystems:
    InfrastructureSystemsInternal,
    InfrastructureSystemsComponent,
    InfrastructureSystemsType

import PowerSystems
import PowerSystems:
    Device,
    Bus,
    ThermalGen,
    Min_Max,
    OperationalCost,
    ThreePartCost,
    PrimeMovers,
    ThermalFuels,
    Service,
    Reserve,
    ReserveDirection,
    DynamicInjection,
    INFINITE_TIME,
    get_value,
    set_value

const IS = InfrastructureSystems
const PSY = PowerSystems

#################################################################################
# Includes

# Include all auto-generated structs.
include("models/generated/includes.jl")

# Better printing
include("utils/print.jl")

end # module

max_stars_repo_name: sourabhdalvi/PowerSystemExtensions

max_issues_repo_name: sourabhdalvi/PowerSystemExtensions

max_forks_repo_name: sourabhdalvi/PowerSystemExtensions


text:
#252634,#234B57,#207269,#4F9866,#9BB858,#FACF5A
# a blue/green/yellow theme of no specific origin
orange_palette =  reverse([
    colorant"#271924", # dark blue
    colorant"#20545D", # semi dark blue
    colorant"#32856A", # green blue
    colorant"#86B15B", # green,
    colorant"#FACF5A"  # yellow
])
black = orange_palette[5]

_themes[:orange] = PlotTheme(
    bg = black,
    bginside = black,
    fg = orange_palette[1],
    fgtext = orange_palette[2],
    fgguide = orange_palette[2],
    fglegend = orange_palette[2],
    palette = expand_palette(black, orange_palette[1:4]),
    colorgradient = :viridis
)


text:
increment!(n) = n.n += 1
function add_subtract!((n1, n2))
    c1 = isnothing(n1) ? 10 : n1.n
    c2 = isnothing(n2) ? 10 : n2.n
    if !isnothing(n1) n1.n += c2 end
    if !isnothing(n2) n2.n -= c1 end
end

flat_prod(ts::Union{Tuple, NamedTuple, Array, Dict}) = prod([flat_prod(t) for t in ts])
flat_prod(t) = t.n
flat_prod(::Nothing) = 1
function postorder_multiply_subtract(ts, chs)
    c = 0
    if !isempty(chs)
        c = flat_prod(chs)
    end
    c -= flat_prod(ts)
    if isempty(chs)
        return Leaf(c)
    elseif chs isa NamedTuple
        return NTVertex(c, chs)
    else
        return VectorVertex(c, collect(chs))
    end
end

@testset "treemap!" for o in ORDERS
        t = deepcopy(COMPLETE_BINARY_TREE_1)
        treemap!(increment!, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [2,3,4,5,6,7,8]
        t = deepcopy(COMPLETE_BINARY_TREE_2)
        treemap!(increment!, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [2,3,4,5,6,7,8]
        t = deepcopy(T1)
        treemap!(increment!, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [2,3,4,5,6,7,8]
        t = deepcopy(T2)
        treemap!(increment!, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [2,3,4,6,7,8]
        t = deepcopy(T3)
        treemap!(increment!, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [2,3,5,6]
        t = deepcopy(T4)
        treemap!(increment!, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [2,4,7,8]
        t = deepcopy(T5)
        treemap!(increment!, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [2]

        t2, t3 = deepcopy.([T2, T3])
        treemap!(add_subtract!, t2, t3; complete=false, order=o)
        @test  [n.n for n in collect(NodeIterator(t2; order=LevelOrder()))] == [2,4,3,10,6,7]
        @test  [n.n for n in collect(NodeIterator(t3; order=LevelOrder()))] == [0,0,4,0]
        t2, t3 = deepcopy.([T2, T3])
        treemap!(add_subtract!, t2, t3; complete=true, order=o)
        @test  [n.n for n in collect(NodeIterator(t2; order=LevelOrder()))] == [2,4,13,10,16,17]
        @test  [n.n for n in collect(NodeIterator(t3; order=LevelOrder()))] == [0,0,-6,0]
end

@testset "leafmap!" for o in ORDERS
        t = deepcopy(COMPLETE_BINARY_TREE_1)
        leafmap!(increment!, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [1,2,3,5,6,7,8]
        t = deepcopy(COMPLETE_BINARY_TREE_2)
        leafmap!(increment!, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [1,2,3,5,6,7,8]
        t = deepcopy(T1)
        leafmap!(increment!, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [1,2,3,5,6,7,8]
        t = deepcopy(T2)
        leafmap!(increment!, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [1,2,3,6,7,8]
        t = deepcopy(T3)
        leafmap!(increment!, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [1,2,5,6]
        t = deepcopy(T4)
        leafmap!(increment!, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [1,3,7,8]
        t = deepcopy(T5)
        leafmap!(increment!, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [2]

        t2, t3 = deepcopy.([T2, T3])
        leafmap!(add_subtract!, t2, t3; complete=false, order=o)
        @test  [n.n for n in collect(NodeIterator(t2; order=LevelOrder()))] == [1,2,3,10,6,7]
        @test  [n.n for n in collect(NodeIterator(t3; order=LevelOrder()))] == [1,2,4,0]
        t2, t3 = deepcopy.([T2, T3])
        leafmap!(add_subtract!, t2, t3; complete=true, order=o)
        @test  [n.n for n in collect(NodeIterator(t2; order=LevelOrder()))] == [1,2,3,10,16,17]
        @test  [n.n for n in collect(NodeIterator(t3; order=LevelOrder()))] == [1,2,-6,0]
end

@testset "typemap!" for o in ORDERS
        t = deepcopy(COMPLETE_BINARY_TREE_1)
        typemap!(increment!, BinaryVertex, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [2,3,4,4,5,6,7]
        t = deepcopy(COMPLETE_BINARY_TREE_2)
        typemap!(increment!, NTVertex, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [2,3,3,4,6,6,8]
        t = deepcopy(T1)
        typemap!(increment!, Leaf, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [1,2,3,5,5,7,7]
        t = deepcopy(T2)
        typemap!(increment!, Any, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [2,3,4,6,7,8]
        t = deepcopy(T3)
        typemap!(increment!, Union{Leaf, NTVertex}, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [2,3,4,6]
        t = deepcopy(T4)
        typemap!(increment!, BinaryVertex, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [1,4,6,7]
        t = deepcopy(T5)
        typemap!(increment!, NTVertex, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [2]

        t1, t2 = deepcopy.([T1, T2])
        typemap!(add_subtract!, (BinaryVertex, BinaryVertex), t1, t2; complete=false, order=o)
        @test  [n.n for n in collect(NodeIterator(t1; order=LevelOrder()))] == [1,2,6,4,5,6,7]
        @test  [n.n for n in collect(NodeIterator(t2; order=LevelOrder()))] == [1,2,0,5,6,7]
        t1, t2 = deepcopy.([T1, T2])
        typemap!(add_subtract!, (NTVertex, Leaf), t1, t2; complete=true, order=o)
        @test  [n.n for n in collect(NodeIterator(t1; order=LevelOrder()))] == [1,2,3,4,5,6,14]
        @test  [n.n for n in collect(NodeIterator(t2; order=LevelOrder()))] == [1,2,3,5,6,0]
end

@testset "predicatemap!" for o in ORDERS
        t = deepcopy(COMPLETE_BINARY_TREE_1)
        predicatemap!(increment!, HierarchicalUtils._leaf_predicate, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [1,2,3,5,6,7,8]
        t = deepcopy(COMPLETE_BINARY_TREE_2)
        predicatemap!(increment!, HierarchicalUtils._leaf_predicate, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [1,2,3,5,6,7,8]
        t = deepcopy(T1)
        predicatemap!(increment!, HierarchicalUtils._leaf_predicate, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [1,2,3,5,6,7,8]
        t = deepcopy(T2)
        predicatemap!(increment!, HierarchicalUtils._leaf_predicate, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [1,2,3,6,7,8]
        t = deepcopy(T3)
        predicatemap!(increment!, HierarchicalUtils._leaf_predicate, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [1,2,5,6]
        t = deepcopy(T4)
        predicatemap!(increment!, HierarchicalUtils._leaf_predicate, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [1,3,7,8]
        t = deepcopy(T5)
        predicatemap!(increment!, HierarchicalUtils._leaf_predicate, t; order=o)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [2]

        t2, t3 = deepcopy.([T2, T3])
        predicatemap!(add_subtract!, HierarchicalUtils._leaf_predicate, t2, t3; complete=false, order=o)
        @test  [n.n for n in collect(NodeIterator(t2; order=LevelOrder()))] == [1,2,3,10,6,7]
        @test  [n.n for n in collect(NodeIterator(t3; order=LevelOrder()))] == [1,2,4,0]
        t2, t3 = deepcopy.([T2, T3])
        predicatemap!(add_subtract!, HierarchicalUtils._leaf_predicate, t2, t3; complete=true, order=o)
        @test  [n.n for n in collect(NodeIterator(t2; order=LevelOrder()))] == [1,2,3,10,16,17]
        @test  [n.n for n in collect(NodeIterator(t3; order=LevelOrder()))] == [1,2,-6,0]
end

@testset "postorder treemap" begin
        t = treemap(postorder_multiply_subtract, COMPLETE_BINARY_TREE_1)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [701,18,39,-4,-5,-6,-7]
        t = treemap(postorder_multiply_subtract, COMPLETE_BINARY_TREE_1, COMPLETE_BINARY_TREE_1)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [694979,396,1755,-16,-25,-36,-49]

        t = treemap(postorder_multiply_subtract, T2, T3; complete=true)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [3743,96,39,-4,-25,-6,-7]
        t = treemap(postorder_multiply_subtract, T2, T3; complete=false)
        @test  [n.n for n in collect(NodeIterator(t; order=LevelOrder()))] == [-30, -29, -25]
end


text:
function bubble_sort(seq::Array{Int64})::Array{Int64}
    L = length(seq)
    for _ in (0:L-1)
        for n in (1:L-1)
            if seq[n+1] < seq[(n-1)+1]
                seq[(n-1)+1], seq[n+1] = (seq[n+1], seq[(n-1)+1])
            end
        end
    end
    return seq
end

function main()
    unsorted = [14, 11, 19, 5, 16, 10, 19, 12, 5, 12]
    expected = [5, 5, 10, 11, 12, 12, 14, 16, 19, 19]
    @assert(bubble_sort(unsorted) == expected)
    println(join(["OK"], " "))
end

main()


text:
# Copyright (c) 2020-2022, Bank of Canada
# All rights reserved.

"""
    struct Period{FREQ} … end

A FAME period includes the frequency, year and period.
"""
struct Period{FREQ}
    year::Int32
    period::Int32
end

function Period{FREQ}(ind::Integer) where {FREQ}
    if iscase(FREQ)
        return Period{FREQ}(0, convert(Int32, ind))
    elseif FREQ == :undefined
        # in a scalar freq is undefined and periods are set to 0
        return Period{FREQ}(0, 0)
    elseif convert(Int64, ind) == FAME_INDEX_NC
        # in an empty series index is set to INDEX_NC and year,period to -1,-1
        return Period{FREQ}(-1, -1)
    else
        # otherwise we call fame to convert index to year:period
        y = Ref{Int32}(-1)
        p = Ref{Int32}(-1)
        @fame_call_check(fame_index_to_year_period,
            (Cint, Clonglong, Ref{Cint}, Ref{Cint}),
            fame_freq[FREQ], convert(Int64, ind), y, p)
        return Period{FREQ}(y[], p[])
    end
end
Base.show(io::IO, p::Period) = print(io, "$(p.year):$(p.period)")

function FameIndex(p::Period{FREQ})::FameIndex where {FREQ}
    if iscase(FREQ)
        return p.period
    else
        date = Ref{Cint}(-1)
        @fame_call_check(fame_year_period_to_index, (Cint, Ref{Cint}, Cint, Cint),
            fame_freq[FREQ], date, p.year, p.period)
        return date[]
    end
end

"""
    mutable struct FameObject{CL,FT,FR,DT} … end

FAME object of class `CL`, type `FT`, frequency `FR` and data type `DT`.

A `FameObject` is returned by [`quick_info`](@ref) and a vector of
`FameObject`s is returned by [`listdb`](@ref).

Also, a `FameObject` can be constructed directly by calling
`FameObject(name, class, type, frequency)` or
`FameObject(name, class, type, frequency, first_index, last_index)`. The values
of `class`, `type`, and `frequency` can be symbols (like, `:scalar`,
`:precision`, etc.) or integers. Refer to the FAME CHLI documentation for the
code values.
"""
mutable struct FameObject{CL,FT,FR,DT}
    # CL : fame_class
    # FT : fame_type
    # FR : fame_freq
    # DT : type of the `data` field
    name::String
    first_index::Ref{FameIndex}
    last_index::Ref{FameIndex}
    data::DT
end
export FameObject

function FameObject(name, cl, ty, fr, fi = -1, li = -1, data = nothing)
    # When the value is a date the type is a frequency.
    _ty = try
        val_to_symbol(ty, fame_type)
    catch
        val_to_symbol(ty, fame_freq)
    end
    _cl = val_to_symbol(cl, fame_class)
    _fr = val_to_symbol(fr, fame_freq)
    ElType = _ty == :precision ? Float64 :
             _ty == :numeric ? Float32 :
             _ty == :boolean ? Int32 :
             _ty == :string ? String :
             _ty == :date ? FameDate :
             _ty == :namelist ? String :
             # otherwise _ty is a frequency, so the value is a date
             FameDate
    T = _cl == :series ? Vector{ElType} :
        _cl == :scalar ? Ref{ElType} :
        throw(ArgumentError("Can't handle class $cl."))
    if data === nothing
        return FameObject{_cl,_ty,_fr,T}(name, Ref(fi), Ref(li), T())
    else
        return FameObject{_cl,_ty,_fr,T}(name, Ref(fi), Ref(li), data)
    end
end

function Base.show(io::IO, fo::FameObject{CL,TY,FR,DT}) where {CL,TY,FR,DT}
    print(io, fo.name, ": ", join((CL, TY, FR,
            Period{getfreq(fo)}(fo.first_index[]),
            Period{getfreq(fo)}(fo.last_index[])), ","))
end

getfreq(::FameObject{CL,FT,FR,DT}) where {CL,FT,FR,DT} = FR

export quick_info
"""
    quick_info(db, name)

Get information about object named `name` in the given database. The information
inlcudes its class, type, frequency, and range. Return a [`FameObject`](@ref) in
which all these attributes are set correctly, but the data does not contain the
correct values.
"""
function quick_info(db::FameDatabase, name::String)
    cl = Ref{Cint}(-1)
    ty = Ref{Cint}(-1)
    fr = Ref{Cint}(-1)
    findex = Ref{Clonglong}(-1)
    lindex = Ref{Clonglong}(-1)
    @fame_call_check(fame_quick_info,
        (Cint, Cstring, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Clonglong}, Ref{Clonglong}),
        db.key, name, cl, ty, fr, findex, lindex)
    return FameObject(name, cl[], ty[], fr[], findex[], lindex[])
end

### List database

export listdb

"""
    item_option(name, values...)

Set ITEM option. These are used when matching a wildcard to database objects.
Option names include "CLASS", "TYPE", "FREQUENCY", "ALIAS".

When the `values` argument is not given, it defaults to "ON". Otherwise `values`
should be strings and indicate which ITEMs are to be turned on.

"""
function item_option end

function item_option(name::String, values::AbstractString...)
    uname = uppercase(strip(name))
    if isempty(values) || (length(values) == 1 && isempty(strip(values[1])))
        @cfm_call_check(cfmsopt, (Cstring, Cstring), "ITEM $uname", "ON")
        return
    end
    @cfm_call_check(cfmsopt, (Cstring, Cstring), "ITEM $uname", "OFF")
    for val in values
        uval = uppercase(strip(val))
        try
            @cfm_call_check(cfmsopt, (Cstring, Cstring), "ITEM $uname $uval", "ON")
        catch e
            e isa HLIError && @error "Bad $(uname) $(uval)."
            rethrow()
        end
    end
end

"""
    listdb(db::FameDatabase, [wc; filters...])

List objects in the database that match the given wildcard and filters. Return a
[`Vector{FameObject}`](@ref FameObject).

The wildcard `wc` is a string containing wildcard characters. '^' matches any
one character, while '?' matches any zero or more characters. If not given, the
default is '?' which would list the entire database.

The filters:
  * `alias::Bool` - whether or not to match alias names.
  * `class::String` - which class of object to match. Multiple classes can be
    given in a comma-separated string, e.g., `class="SERIES,SCALAR"`.
  * `type::String` - which type of object to match, e.g., `type="NUMERIC,PRECISION"`.
  * `freq::String` - which frequency of object to match, e.g., `freq="QUARTERLY"`.
"""
function listdb end

@inline listdb(dbname::AbstractString, args...; kwargs...) =
    opendb(dbname) do db
        listdb(db, args...; kwargs...)
    end

function listdb(db::FameDatabase, wildcard::String = "?";
    alias::Bool = true, class = "", type = "", freq = "")

    class = string(class)
    type = string(type)
    freq = string(freq)

    @cfm_call_check(cfmsopt, (Cstring, Cstring), "ITEM ALIAS", alias ? "ON" : "OFF")
    item_option("CLASS", split(class, ",")...)
    item_option("TYPE", split(type, ",")...)
    item_option("FREQUENCY", split(freq, ",")...)

    wc_key = Ref{Cint}(-1)
    @fame_call_check(fame_init_wildcard,
        (Cint, Ref{Cint}, Cstring, Cint, Cstring),
        db.key, wc_key, wildcard, 0, C_NULL)
    ret = FameObject[]
    try
        while true
            name = repeat(" ", 101)
            cl = Ref{Cint}()
            ty = Ref{Cint}()
            fr = Ref{Cint}()
            fi = Ref{Clonglong}()
            li = Ref{Clonglong}()
            status = @fame_call(fame_get_next_wildcard,
                (Cint, Cstring, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Clonglong}, Ref{Clonglong}, Cint, Ref{Cint}),
                wc_key[], name, cl, ty, fr, fi, li, length(name) - 1, C_NULL)
            if status == HNOOBJ
                break
            elseif status == HTRUNC
                println("Object name too long and truncated $name")
            else
                check_status(status)
            end
            # FAME pads the string with \0 on the right to the length we gave.
            name = strip(name, '\0')
            # note: fame_get_next_wildcard returns incorrect fi and li for scalars.
            @fame_call_check(fame_quick_info,
                (Cint, Cstring, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Clonglong}, Ref{Clonglong}),
                db.key, name, cl, ty, fr, fi, li)
            fo = FameObject(name, cl[], ty[], fr[], fi[], li[])
            push!(ret, fo)
        end
    finally
        @fame_call(fame_free_wildcard, (Cint,), wc_key[])
    end
    return ret
end
export listdb


text:
# The incredible pressures at this depth are starting to put a strain on your
# submarine. The submarine has polymerization equipment that would produce
# suitable materials to reinforce the submarine, and the nearby
# volcanically-active caves should even have the necessary input elements in
# sufficient quantities.
#
# The submarine manual contains instructions for finding the optimal polymer
# formula; specifically, it offers a polymer template and a list of pair
# insertion rules (your puzzle input). You just need to work out what polymer
# would result after repeating the pair insertion process a few times.
#
# For example:
#
# NNCB
#
# CH -> B
# HH -> N
# CB -> H
# NH -> C
# HB -> C
# HC -> B
# HN -> C
# NN -> C
# BH -> H
# NC -> B
# NB -> B
# BN -> B
# BB -> N
# BC -> B
# CC -> N
# CN -> C
#
# The first line is the polymer template - this is the starting point of the
# process.
#
# The following section defines the pair insertion rules. A rule like AB -> C
# means that when elements A and B are immediately adjacent, element C should
# be inserted between them. These insertions all happen simultaneously.
#
# So, starting with the polymer template NNCB, the first step simultaneously
# considers all three pairs:
#
# - The first pair (NN) matches the rule NN -> C, so element C is inserted
#   between the first N and the second N.
# - The second pair (NC) matches the rule NC -> B, so element B is inserted
#   between the N and the C.
# - The third pair (CB) matches the rule CB -> H, so element H is inserted
#   between the C and the B.
#
# Note that these pairs overlap: the second element of one pair is the first
# element of the next pair. Also, because all pairs are considered
# simultaneously, inserted elements are not considered to be part of a pair
# until the next step.
#
# After the first step of this process, the polymer becomes NCNBCHB.
#
# Here are the results of a few steps using the above rules:
#
# Template:     NNCB
# After step 1: NCNBCHB
# After step 2: NBCCNBBBCBHCB
# After step 3: NBBBCNCCNBBNBNBBCHBHHBCHB
# After step 4: NBBNBNBBCCNBCNCCNBBNBBNBBBNBBNBBCBHCBHHNHCBBCBHCB
#
# This polymer grows quickly. After step 5, it has length 97; After step 10, it
# has length 3073. After step 10, B occurs 1749 times, C occurs 298 times, H
# occurs 161 times, and N occurs 865 times; taking the quantity of the most
# common element (B, 1749) and subtracting the quantity of the least common
# element (H, 161) produces 1749 - 161 = 1588.
#
# Apply 10 steps of pair insertion to the polymer template and find the most
# and least common elements in the result. What do you get if you take the
# quantity of the most common element and subtract the quantity of the least
# common element?

const steps = 10

# read input
iter = eachline("input.txt")
(template, _) = iterate(iter)
iterate(iter) # blank line
rules = Dict([split(l, " -> ") for l in iter])

# Count elements in the starting template and initialize other elements to 0
cnts = Dict([(c, 1) for c in template])
for (_,c) in rules
  c = c[1]
  if !haskey(cnts, c)
    cnts[c] = 0
  end
end

# recursive function to do each step
function step(l, r, depth=1)
  depth > steps && return
  m = rules[l*r][1]
  cnts[m] += 1
  step(l, m, depth+1)
  step(m, r, depth+1)
end

# do eet - iterate over (1,2), (2,3), ..., (len-1, len)
foreach(zip(1:(length(template)-1), 2:length(template))) do (i,j)
  step(template[i], template[j])
end
println("Result: ", abs(-(extrema(last, cnts)...)))