text:
module CellFieldsBench

using Gridap
using Gridap.FieldValues
using Gridap.CellValues
using Gridap.CellFunctions
using Gridap.Geometry
using Gridap.Geometry.Cartesian

l = 1000000

include("../test/CellFunctionsTestsMocks.jl")

function doloop(x)
  for xi in x
  end
end

println("+++ CellFieldsBench ( length = $l ) +++")

tiff = varinner(tbv,tbv)

print("TensorBasisBasisInner ->"); @time doloop(tiff)
print("TensorBasisBasisInner ->"); @time doloop(tiff)

vexpand = expand(vbv,sfv2)

print("VectorScalarExpand ->"); @time doloop(vexpand)
print("VectorScalarExpand ->"); @time doloop(vexpand)

using Gridap.Geometry
using Gridap.Quadratures
using Gridap.CellQuadratures
using Gridap.CellIntegration

grid = CartesianGrid(partition=(1000,1000))
trian = Triangulation(grid)

meshcoords = cellcoordinates(trian)
quad = CellQuadrature(trian,degree=2)
points = coordinates(quad)
phi = geomap(trian)
basis = CellBasis(trian)
physbasis = attachgeomap(basis,phi)
physbasisgrad = gradient(physbasis)
valsgrad = evaluate(physbasisgrad,points)
xg = evaluate(phi,points)
vals = evaluate(basis,points)

print("DummyMeshCoords2D ->"); @time doloop(meshcoords)
print("DummyMeshCoords2D ->"); @time doloop(meshcoords)

print("DummyCellBasis2D ->"); @time doloop(vals)
print("DummyCellBasis2D ->"); @time doloop(vals)

print("xg ->"); @time doloop(xg)
print("xg ->"); @time doloop(xg)

print("PhysBasisGradVals ->"); @time doloop(valsgrad)
print("PhysBasisGradVals ->"); @time doloop(valsgrad)

kmatg = varinner(valsgrad, valsgrad)

print("DummyStiffnessMatrix2DAtQPoints ->"); @time doloop(kmatg)
print("DummyStiffnessMatrix2DAtQPoints ->"); @time doloop(kmatg)

a(v,u) = inner(∇(v),∇(u)) + inner(v,u)
V = physbasis
U = physbasis
kmat = integrate(a(V,U),trian,quad)

print("DummyStiffnessMatrix2D ->"); @time doloop(kmat)
print("DummyStiffnessMatrix2D ->"); @time doloop(kmat)

end


text:
# Currently unsupported
# Probably need to parse the associated .afm files to make this work

_latex_to_bakoma = Dict(
    raw"\sqrt"                     => ("cmex10", 0x70),
    raw"\bigcap"                   => ("cmex10", 0x5c),
    raw"\bigcup"                   => ("cmex10", 0x5b),
    raw"\bigodot"                  => ("cmex10", 0x4b),
    raw"\bigoplus"                 => ("cmex10", 0x4d),
    raw"\bigotimes"                => ("cmex10", 0x4f),
    raw"\biguplus"                 => ("cmex10", 0x5d),
    raw"\bigvee"                   => ("cmex10", 0x5f),
    raw"\bigwedge"                 => ("cmex10", 0x5e),
    raw"\coprod"                   => ("cmex10", 0x61),
    raw"\int"                      => ("cmex10", 0x5a),
    raw"\langle"                   => ("cmex10", 0xad),
    raw"\leftangle"                => ("cmex10", 0xad),
    raw"\leftbrace"                => ("cmex10", 0xa9),
    raw"\oint"                     => ("cmex10", 0x49),
    raw"\prod"                     => ("cmex10", 0x59),
    raw"\rangle"                   => ("cmex10", 0xae),
    raw"\rightangle"               => ("cmex10", 0xae),
    raw"\rightbrace"               => ("cmex10", 0xaa),
    raw"\sum"                      => ("cmex10", 0x58),
    raw"\widehat"                  => ("cmex10", 0x62),
    raw"\widetilde"                => ("cmex10", 0x65),
    raw"\{"                        => ("cmex10", 0xa9),
    raw"\}"                        => ("cmex10", 0xaa),
    raw"{"                         => ("cmex10", 0xa9),
    raw"}"                         => ("cmex10", 0xaa),

    raw","                         => ("cmmi10", 0x3b),
    raw"."                         => ("cmmi10", 0x3a),
    raw"/"                         => ("cmmi10", 0x3d),
    raw"<"                         => ("cmmi10", 0x3c),
    raw">"                         => ("cmmi10", 0x3e),
    raw"\alpha"                    => ("cmmi10", 0xae),
    raw"\beta"                     => ("cmmi10", 0xaf),
    raw"\chi"                      => ("cmmi10", 0xc2),
    raw"\combiningrightarrowabove" => ("cmmi10", 0x7e),
    raw"\delta"                    => ("cmmi10", 0xb1),
    raw"\ell"                      => ("cmmi10", 0x60),
    raw"\epsilon"                  => ("cmmi10", 0xb2),
    raw"\eta"                      => ("cmmi10", 0xb4),
    raw"\flat"                     => ("cmmi10", 0x5b),
    raw"\frown"                    => ("cmmi10", 0x5f),
    raw"\gamma"                    => ("cmmi10", 0xb0),
    raw"\imath"                    => ("cmmi10", 0x7b),
    raw"\iota"                     => ("cmmi10", 0xb6),
    raw"\jmath"                    => ("cmmi10", 0x7c),
    raw"\kappa"                    => ("cmmi10", 0x2219),
    raw"\lambda"                   => ("cmmi10", 0xb8),
    raw"\leftharpoondown"          => ("cmmi10", 0x29),
    raw"\leftharpoonup"            => ("cmmi10", 0x28),
    raw"\mu"                       => ("cmmi10", 0xb9),
    raw"\natural"                  => ("cmmi10", 0x5c),
    raw"\nu"                       => ("cmmi10", 0xba),
    raw"\omega"                    => ("cmmi10", 0x21),
    raw"\phi"                      => ("cmmi10", 0xc1),
    raw"\pi"                       => ("cmmi10", 0xbc),
    raw"\psi"                      => ("cmmi10", 0xc3),
    raw"\rho"                      => ("cmmi10", 0xbd),
    raw"\rightharpoondown"         => ("cmmi10", 0x2b),
    raw"\rightharpoonup"           => ("cmmi10", 0x2a),
    raw"\sharp"                    => ("cmmi10", 0x5d),
    raw"\sigma"                    => ("cmmi10", 0xbe),
    raw"\smile"                    => ("cmmi10", 0x5e),
    raw"\tau"                      => ("cmmi10", 0xbf),
    raw"\theta"                    => ("cmmi10", 0xb5),
    raw"\triangleleft"             => ("cmmi10", 0x2f),
    raw"\triangleright"            => ("cmmi10", 0x2e),
    raw"\upsilon"                  => ("cmmi10", 0xc0),
    raw"\varepsilon"               => ("cmmi10", 0x22),
    raw"\varphi"                   => ("cmmi10", 0x27),
    raw"\varrho"                   => ("cmmi10", 0x25),
    raw"\varsigma"                 => ("cmmi10", 0x26),
    raw"\vartheta"                 => ("cmmi10", 0x23),
    raw"\wp"                       => ("cmmi10", 0x7d),
    raw"\xi"                       => ("cmmi10", 0xbb),
    raw"\zeta"                     => ("cmmi10", 0xb3),

    raw"!"                         => ("cmr10", 0x21),
    raw"%"                         => ("cmr10", 0x25),
    raw"&"                         => ("cmr10", 0x26),
    raw"("                         => ("cmr10", 0x28),
    raw")"                         => ("cmr10", 0x29),
    raw"+"                         => ("cmr10", 0x2b),
    raw"0"                         => ("cmr10", 0x30),
    raw"1"                         => ("cmr10", 0x31),
    raw"2"                         => ("cmr10", 0x32),
    raw"3"                         => ("cmr10", 0x33),
    raw"4"                         => ("cmr10", 0x34),
    raw"5"                         => ("cmr10", 0x35),
    raw"6"                         => ("cmr10", 0x36),
    raw"7"                         => ("cmr10", 0x37),
    raw"8"                         => ("cmr10", 0x38),
    raw"9"                         => ("cmr10", 0x39),
    raw":"                         => ("cmr10", 0x3a),
    raw";"                         => ("cmr10", 0x3b),
    raw"="                         => ("cmr10", 0x3d),
    raw"?"                         => ("cmr10", 0x3f),
    raw"@"                         => ("cmr10", 0x40),
    raw"["                         => ("cmr10", 0x5b),
    raw"\#"                        => ("cmr10", 0x23),
    raw"\$"                        => ("cmr10", 0x24),
    raw"\%"                        => ("cmr10", 0x25),
    raw"\Delta"                    => ("cmr10", 0xa2),
    raw"\Gamma"                    => ("cmr10", 0xa1),
    raw"\Lambda"                   => ("cmr10", 0xa4),
    raw"\Omega"                    => ("cmr10", 0xad),
    raw"\Phi"                      => ("cmr10", 0xa9),
    raw"\Pi"                       => ("cmr10", 0xa6),
    raw"\Psi"                      => ("cmr10", 0xaa),
    raw"\Sigma"                    => ("cmr10", 0xa7),
    raw"\Theta"                    => ("cmr10", 0xa3),
    raw"\Upsilon"                  => ("cmr10", 0xa8),
    raw"\Xi"                       => ("cmr10", 0xa5),
    raw"\circumflexaccent"         => ("cmr10", 0x5e),
    raw"\combiningacuteaccent"     => ("cmr10", 0xb6),
    raw"\combiningbreve"           => ("cmr10", 0xb8),
    raw"\combiningdiaeresis"       => ("cmr10", 0xc4),
    raw"\combiningdotabove"        => ("cmr10", 0x5f),
    raw"\combininggraveaccent"     => ("cmr10", 0xb5),
    raw"\combiningoverline"        => ("cmr10", 0xb9),
    raw"\combiningtilde"           => ("cmr10", 0x7e),
    raw"\leftbracket"              => ("cmr10", 0x5b),
    raw"\leftparen"                => ("cmr10", 0x28),
    raw"\rightbracket"             => ("cmr10", 0x5d),
    raw"\rightparen"               => ("cmr10", 0x29),
    raw"\widebar"                  => ("cmr10", 0xb9),
    raw"]"                         => ("cmr10", 0x5d),

    raw"*"                         => ("cmsy10", 0xa4),
    raw"-"                         => ("cmsy10", 0xa1),
    raw"\Downarrow"                => ("cmsy10", 0x2b),
    raw"\Im"                       => ("cmsy10", 0x3d),
    raw"\Leftarrow"                => ("cmsy10", 0x28),
    raw"\Leftrightarrow"           => ("cmsy10", 0x2c),
    raw"\P"                        => ("cmsy10", 0x7b),
    raw"\Re"                       => ("cmsy10", 0x3c),
    raw"\Rightarrow"               => ("cmsy10", 0x29),
    raw"\S"                        => ("cmsy10", 0x78),
    raw"\Uparrow"                  => ("cmsy10", 0x2a),
    raw"\Updownarrow"              => ("cmsy10", 0x6d),
    raw"\Vert"                     => ("cmsy10", 0x6b),
    raw"\aleph"                    => ("cmsy10", 0x40),
    raw"\approx"                   => ("cmsy10", 0xbc),
    raw"\ast"                      => ("cmsy10", 0xa4),
    raw"\asymp"                    => ("cmsy10", 0xb3),
    raw"\backslash"                => ("cmsy10", 0x6e),
    raw"\bigcirc"                  => ("cmsy10", 0xb0),
    raw"\bigtriangledown"          => ("cmsy10", 0x35),
    raw"\bigtriangleup"            => ("cmsy10", 0x34),
    raw"\bot"                      => ("cmsy10", 0x3f),
    raw"\bullet"                   => ("cmsy10", 0xb2),
    raw"\cap"                      => ("cmsy10", 0x5c),
    raw"\cdot"                     => ("cmsy10", 0xa2),
    raw"\circ"                     => ("cmsy10", 0xb1),
    raw"\clubsuit"                 => ("cmsy10", 0x7c),
    raw"\cup"                      => ("cmsy10", 0x5b),
    raw"\dag"                      => ("cmsy10", 0x79),
    raw"\dashv"                    => ("cmsy10", 0x61),
    raw"\ddag"                     => ("cmsy10", 0x7a),
    raw"\diamond"                  => ("cmsy10", 0xa6),
    raw"\diamondsuit"              => ("cmsy10", 0x7d),
    raw"\div"                      => ("cmsy10", 0xa5),
    raw"\downarrow"                => ("cmsy10", 0x23),
    raw"\emptyset"                 => ("cmsy10", 0x3b),
    raw"\equiv"                    => ("cmsy10", 0xb4),
    raw"\exists"                   => ("cmsy10", 0x39),
    raw"\forall"                   => ("cmsy10", 0x38),
    raw"\geq"                      => ("cmsy10", 0xb8),
    raw"\gg"                       => ("cmsy10", 0xc0),
    raw"\heartsuit"                => ("cmsy10", 0x7e),
    raw"\in"                       => ("cmsy10", 0x32),
    raw"\infty"                    => ("cmsy10", 0x31),
    raw"\lbrace"                   => ("cmsy10", 0x66),
    raw"\lceil"                    => ("cmsy10", 0x64),
    raw"\leftarrow"                => ("cmsy10", 0xc3),
    raw"\leftrightarrow"           => ("cmsy10", 0x24),
    raw"\leq"                      => ("cmsy10", 0x2219),
    raw"\lfloor"                   => ("cmsy10", 0x62),
    raw"\ll"                       => ("cmsy10", 0xbf),
    raw"\mid"                      => ("cmsy10", 0x6a),
    raw"\mp"                       => ("cmsy10", 0xa8),
    raw"\nabla"                    => ("cmsy10", 0x72),
    raw"\nearrow"                  => ("cmsy10", 0x25),
    raw"\neg"                      => ("cmsy10", 0x3a),
    raw"\ni"                       => ("cmsy10", 0x33),
    raw"\nwarrow"                  => ("cmsy10", 0x2d),
    raw"\odot"                     => ("cmsy10", 0xaf),
    raw"\ominus"                   => ("cmsy10", 0xaa),
    raw"\oplus"                    => ("cmsy10", 0xa9),
    raw"\oslash"                   => ("cmsy10", 0xae),
    raw"\otimes"                   => ("cmsy10", 0xad),
    raw"\pm"                       => ("cmsy10", 0xa7),
    raw"\prec"                     => ("cmsy10", 0xc1),
    raw"\preceq"                   => ("cmsy10", 0xb9),
    raw"\prime"                    => ("cmsy10", 0x30),
    raw"\propto"                   => ("cmsy10", 0x2f),
    raw"\rbrace"                   => ("cmsy10", 0x67),
    raw"\rceil"                    => ("cmsy10", 0x65),
    raw"\rfloor"                   => ("cmsy10", 0x63),
    raw"\rightarrow"               => ("cmsy10", 0x21),
    raw"\searrow"                  => ("cmsy10", 0x26),
    raw"\sim"                      => ("cmsy10", 0xbb),
    raw"\simeq"                    => ("cmsy10", 0x27),
    raw"\slash"                    => ("cmsy10", 0x36),
    raw"\spadesuit"                => ("cmsy10", 0xc4),
    raw"\sqcap"                    => ("cmsy10", 0x75),
    raw"\sqcup"                    => ("cmsy10", 0x74),
    raw"\sqsubseteq"               => ("cmsy10", 0x76),
    raw"\sqsupseteq"               => ("cmsy10", 0x77),
    raw"\subset"                   => ("cmsy10", 0xbd),
    raw"\subseteq"                 => ("cmsy10", 0xb5),
    raw"\succ"                     => ("cmsy10", 0xc2),
    raw"\succeq"                   => ("cmsy10", 0xba),
    raw"\supset"                   => ("cmsy10", 0xbe),
    raw"\supseteq"                 => ("cmsy10", 0xb6),
    raw"\swarrow"                  => ("cmsy10", 0x2e),
    raw"\times"                    => ("cmsy10", 0xa3),
    raw"\to"                       => ("cmsy10", 0x21),
    raw"\top"                      => ("cmsy10", 0x3e),
    raw"\uparrow"                  => ("cmsy10", 0x22),
    raw"\updownarrow"              => ("cmsy10", 0x6c),
    raw"\uplus"                    => ("cmsy10", 0x5d),
    raw"\vdash"                    => ("cmsy10", 0x60),
    raw"\vee"                      => ("cmsy10", 0x5f),
    raw"\vert"                     => ("cmsy10", 0x6a),
    raw"\wedge"                    => ("cmsy10", 0x5e),
    raw"\wr"                       => ("cmsy10", 0x6f),
    raw"\|"                        => ("cmsy10", 0x6b),
    raw"|"                         => ("cmsy10", 0x6a),
    raw"\_"                        => ("cmtt10", 0x5f)
)

const BakomaSymbolSet = SymbolSet(
    "Bakoma (Computer Modern)",
    Dict([sym => TeXChar(Char(code), join("assets/fonts/$font.ttf"))
          for (sym, (font, code)) in _latex_to_bakoma])
)

ext: jl


text:
using stoosc
rc("font", family="arial")

# Load model and set paraemters (except for A that will be scanned)
model = load_model("model01", "ode")
model = set_parameter_value(model, ["τ", "I"], [1.0, 0.0])

# A summary function that performs a numerical simulation and estimates the
# parameters of the solution
summary_function = function (model=nothing)

    # If the model is not passed, return the names of the parameters
    if isnothing(model)
        return ["p2p", "rms"]
    end

    # Set up the solver and simulate the model
    model = set_timespan(model, (0.0, 1000.0))
    model = set_solver(model, saveat=200.0:0.1:1000.0)
    solution = simulate_population(model)

    # If the simulation is successful, calculate the parameters, otherwise
    # return NaNs
    if solution.success
        x = solution.mean[:, 3]
        xmax = maximum(x)
        xmin = minimum(x)
        p2p = xmax - xmin
        rms = sqrt(mean(x.^2))
    else
        p2p = NaN
        rms = NaN
    end

    return [p2p, rms]

end

# Specify the name of the parameter and the values to scan over
parameter_names = ["A"]
parameter_values = [10 .^ range(-3.0, 0.0; length=100)]

# Perform the parameter scan
df = scan(model, parameter_names, parameter_values, summary_function;
    show_progress=true)

# Create a figure and save it
A_log = log10.(df[:, "A"])
p2p = df[:, "p2p"]
ylim = [minimum(p2p) - 0.001, maximum(p2p) + 0.001]

fig, ax = subplots(figsize=(7, 3))
ax.plot(A_log, p2p, color="black", label="Peak-to-peak amplitude")
ax.plot([-1, -1], ylim, label="Used value of A", color="red")
ax.plot([log10(0.125), log10(0.125)], ylim, label="Analytical bifurcation", color="blue")
ax.legend(edgecolor=:black, framealpha=1.0)
ax.set_xlim(-3, 0)
ax.set_ylim(ylim)
ax.set_xlabel("log\$_{10}\$(A)")
ax.set_ylabel("Peak-to-peak amplitude")
ax.set_title("Numerical parameter scan for A")
fig.tight_layout()

save_figure(fig, "scan_A.svg")

max_stars_repo_name: vkumpost/stoosc

max_issues_repo_name: vkumpost/stoosc

max_forks_repo_name: vkumpost/stoosc


text:
using JAnki
using REPL.TerminalMenus

function main()

card = JAnki.anki_find_notes("transfer")["result"][1]

card_dict = Dict()

card = JAnki.anki_get_card(card)["result"][1]

card_id = card["noteId"]

field_list = collect(keys(card["fields"]))

println("Select what fields you wish to query from these cards:")
selections = request(MultiSelectMenu(field_list, pagesize = 5))
selections = field_list[collect(keys(selections.dict))]

card_info = ""

for selection in selections
	card_info = card_info * " " *  replace(card["fields"][selection]["value"], r"
|
|\n|\t|
|
" => "")
end

card_info
end


text:
  # Some utilities to modify nested QuoteNode's

  "takes :A, :B, :C and return :(A.B.C)"
  qn(a, b) = Expr(:(.), a, QuoteNode(b))
  qn(a, b, rest...) = qn(qn(a, b), rest...)

  "parse :(A.B.C) and returns an array of symbols with :A, :B, :C"
  qs(ex::Expr) = qs(ex, Symbol[])
  qs(s::Symbol) = [s]
  function qs(ex::Expr, out::Vector{Symbol})

      if ex.head == :(.)
          qs(ex.args[1], out)
          qs(ex.args[2], out)
      end
      out
  end
  qs(s::Symbol, out::Vector{Symbol}) = push!(out, s)
  qs(s::QuoteNode, out::Vector{Symbol}) = qs(s.value, out)

  #qn(:A, qs(:(B.C))...) == :(A.B.C)
  #qn(:A, :B, :C) == :(A.B.C)

  """
      eval_symbol(s, eval_in::Module)

  eval s in module eval_in, used for data hint.

  """
  function eval_symbol(ex::Union{Expr, Symbol}, eval_in::Module)

      #this prepends eval_in so we can evaluate in Main directly :(eval_in.ex)
      ex = qn(eval_in, qs(ex)...)

      evalout = try Core.eval(Main, ex) catch err; "" end
  end

  function get_doc(s::Symbol, eval_in::Module)
      Base.Docs.doc(
          Base.Docs.Binding(eval_in, s)
      )
  end
  function get_doc(ex::Expr, eval_in::Module)
      try Core.eval(eval_in, :( @doc $ex )) catch; "" end
  end
  function get_doc(ex::String, eval_in::Module)
      try get_doc(Meta.parse(ex), eval_in::Module) catch err; return err end
  end


text:
#===============================================================================

    test_functions.jl

    Author: Tom Price
    Date: Dec 2018

    automated test script for Julia API
    this code should be run at startup in fresh julia REPL

===============================================================================#

const BBB = match(r"beaglebone", readlines(`uname -a`)[1]) != nothing
const RUN_THIS_CODE_INTERACTIVELY_NOT_ON_INCLUDE = false

import DataFrames.DataFrame
import DataStructures.OrderedDict
import JSON: json, parse, parsefile
using LibCURL

@static if !BBB
    import FactCheck: clear_results
    import BSON: bson
end

const td = readdlm("$LOAD_FROM_DIR/../test/data/test_data.csv", ',', header=true)
const TEST_DATA = DataFrame([
    slicedim(td[1], 2, i) for i in 1:size(td[1],2)],
    map(Symbol, td[2][:]))

#===============================================================================
    interactive code (development environment) >>
===============================================================================#

## example code to generate, run, and save tests
## BSON preferred to JLD because it can save functions and closures
if (RUN_THIS_CODE_INTERACTIVELY_NOT_ON_INCLUDE & !BBB)
    # cd("/home/vagrant/chaipcr/bioinformatics/QpcrAnalysis")
    # JULIA_ENV=development julia -e 'push!(LOAD_PATH,"."); include("../juliaserver.jl")' &
    dir = pwd()
    push!(LOAD_PATH, dir)
    using QpcrAnalysis
    test_functions = QpcrAnalysis.generate_tests()
    check = QpcrAnalysis.test_dispatch(test_functions)

    if all(values(check))
        BSON.bson("$(QpcrAnalysis.LOAD_FROM_DIR)/../test/data/dispatch_tests.bson", test_functions)
        println("All test functions checked and saved")
        ## time functions second time around (after compilation)
        timing = QpcrAnalysis.time_dispatch(test_functions)
    else
        println("Test functions failed check:")
        println(check)
    end
end

## alternative test functions
# QpcrAnalysis.generate_test_script("revised_exec_testfns.jl")
# run(`mv exec_testfns.jl ../build`)

## Amplification timings 12 July 2019

# amp1 = test_functions["amplification single channel"]()
# # open("/tmp/amp1-test.json","w") do f
# #     JSON.print(f, amp1[2]["cq"])
# # end
# amp1_local = amp1[2]["cq"]
# amp1_saved = JSON.parsefile("/tmp/amp1-test.json")
# amp1_local == amp1_saved ## should be true
# using Memento
# setlevel!(QpcrAnalysis.logger, "error")
# # using BenchmarkTools
# # @benchmark test_functions["amplification single channel"]()
# @timev for i in 1:100; test_functions["amplification single channel"](); end;

## commit 7d9864b27b7d8fcc777963fc2085ba800e41faa2
#  26.419600 seconds (25.08 M allocations: 1.561 GiB, 1.56% gc time)
# elapsed time (ns): 26419600000
# gc time (ns):      411273000
# bytes allocated:   1676623840
# pool allocs:       25009327
# non-pool GC allocs:59900
# malloc() calls:    9600
# realloc() calls:   900
# GC pauses:         73

## commit dfe7f96d9e00f03b939cc84c898a58156c05fd41
#  26.820811 seconds (25.08 M allocations: 1.562 GiB, 1.66% gc time)
# elapsed time (ns): 26820811000
# gc time (ns):      445567000
# bytes allocated:   1676798240
# pool allocs:       25007527
# non-pool GC allocs:59900
# malloc() calls:    9600
# realloc() calls:   900
# GC pauses:         73

## commit c4d3ac425fbfa1e65ae352f1e78b1483645754e5
#  31.532593 seconds (25.08 M allocations: 1.562 GiB, 1.49% gc time)
# elapsed time (ns): 31532593000
# gc time (ns):      469400000
# bytes allocated:   1676803712
# pool allocs:       25010026
# non-pool GC allocs:59900
# malloc() calls:    9600
# realloc() calls:   900
# GC pauses:         74

## commit 2e7344b877395abe088d59422b5708fd8165f43d (fixup)
#  27.648082 seconds (25.36 M allocations: 1.572 GiB, 1.60% gc time)
# elapsed time (ns): 27648082000
# gc time (ns):      442063000
# bytes allocated:   1688355424
# pool allocs:       25288926
# non-pool GC allocs:60400
# malloc() calls:    9600
# realloc() calls:   900
# GC pauses:         74

## commit 43161cad1e67aa8073efb63ec3be51bd5c21a3fd
# BenchmarkTools.Trial:
#   memory estimate:  15.99 MiB
#   allocs estimate:  250801
#   --------------
#   minimum time:     269.112 ms (0.00% GC)
#   median time:      285.939 ms (1.94% GC)
#   mean time:        288.227 ms (1.55% GC)
#   maximum time:     322.246 ms (2.45% GC)
#   --------------
#   samples:          18
#   evals/sample:     1

## commit 9e6b315e1f56b3310e0488590bde9d0e0725529b
# BenchmarkTools.Trial:
#   memory estimate:  15.99 MiB
#   allocs estimate:  250801
#   --------------
#   minimum time:     306.455 ms (0.00% GC)
#   median time:      324.136 ms (1.74% GC)
#   mean time:        322.756 ms (1.27% GC)
#   maximum time:     336.128 ms (1.79% GC)
#   --------------
#   samples:          16
#   evals/sample:     1

## commit 00f645a6adfdbb99ab2279b6358f0fc50a21e532
# BenchmarkTools.Trial:
#   memory estimate:  15.99 MiB
#   allocs estimate:  250818
#   --------------
#   minimum time:     290.265 ms (1.82% GC)
#   median time:      301.966 ms (1.76% GC)
#   mean time:        303.233 ms (1.43% GC)
#   maximum time:     316.700 ms (1.77% GC)
#   --------------
#   samples:          17
#   evals/sample:     1

## better
#  29.574523 seconds (25.08 M allocations: 1.562 GiB, 1.54% gc time)
# elapsed time (ns): 29574523000
# gc time (ns):      456487000
# bytes allocated:   1676802240
# pool allocs:       25008593
# non-pool GC allocs:59900
# malloc() calls:    9600
# realloc() calls:   900
# GC pauses:         73

## commit ff41db9d3fd782dd24fdcdfe7835c15faca380e9
#   memory estimate:  16.00 MiB
#   allocs estimate:  250872
#   --------------
#   minimum time:     269.513 ms (0.00% GC)
#   median time:      286.489 ms (1.90% GC)
#   mean time:        283.428 ms (1.46% GC)
#   maximum time:     291.448 ms (2.13% GC)
#   --------------
#   samples:          18
#   evals/sample:     1

## commit cc325bb92525153a641c8083c977b2f202aed130
# BenchmarkTools.Trial:
#   memory estimate:  16.00 MiB
#   allocs estimate:  250883
#   --------------
#   minimum time:     257.397 ms (0.00% GC)
#   median time:      272.069 ms (2.05% GC)
#   mean time:        268.909 ms (1.49% GC)
#   maximum time:     276.962 ms (2.02% GC)
#   --------------
#   samples:          19
#   evals/sample:     1

## commit c964da0fe71de22009e5ea3881d1b0629b179e9c on master
#  30.972977 seconds (25.09 M allocations: 1.562 GiB, 1.51% gc time)
# elapsed time (ns): 30972977000
# gc time (ns):      468801000
# bytes allocated:   1677302448
# pool allocs:       25018510
# non-pool GC allocs:59905
# malloc() calls:    9600
# realloc() calls:   900
# GC pauses:         73

## commit dce3b4a1265df1ca5582283b374b118f0c3d3195
# BenchmarkTools.Trial:
#   memory estimate:  16.00 MiB
#   allocs estimate:  250888
#   --------------
#   minimum time:     262.566 ms (0.00% GC)
#   median time:      279.580 ms (2.11% GC)
#   mean time:        281.085 ms (1.59% GC)
#   maximum time:     308.169 ms (2.30% GC)
#   --------------
#   samples:          18
#   evals/sample:     1

## tweaks
# 27.564401 seconds (25.59 M allocations: 1.598 GiB, 1.62% gc time)
# elapsed time (ns): 27564401000
# gc time (ns):      446231000
# bytes allocated:   1715774576
# pool allocs:       25519126
# non-pool GC allocs:61100
# malloc() calls:    9600
# realloc() calls:   900
# GC pauses:         75

# first attempt using static arrays
# 35.524490 seconds (27.04 M allocations: 9.616 GiB, 3.94% gc time)
# elapsed time (ns): 35524490000
# gc time (ns):      1399246000
# bytes allocated:   10325273296
# pool allocs:       25401356
# non-pool GC allocs:1629005
# malloc() calls:    9700
# realloc() calls:   1200
# GC pauses:         450
# full collections:  1

# 27.231486 seconds (27.77 M allocations: 1.724 GiB, 1.75% gc time)
# elapsed time (ns): 27231486000
# gc time (ns):      476178000
# bytes allocated:   1851260240
# pool allocs:       27699785
# non-pool GC allocs:62205
# malloc() calls:    9600
# realloc() calls:   1100
# GC pauses:         80

# commit 803d603cf6663601a97b22c1937d8b682110b4d2
# 29.464264 seconds (27.88 M allocations: 1.674 GiB, 1.66% gc time)
# elapsed time (ns): 29464264000
# gc time (ns):      487881000
# bytes allocated:   1797779888
# pool allocs:       27815683
# non-pool GC allocs:57505
# malloc() calls:    9600
# realloc() calls:   1100
# GC pauses:         78

# commit a01b3b50a6e753f3ce92d4a33c1cc9adbb6256b7:
# 27.899778 seconds (27.90 M allocations: 1.675 GiB, 1.99% gc time)
# elapsed time (ns): 27899778000
# gc time (ns):      554854000
# bytes allocated:   1798466880
# pool allocs:       27832990
# non-pool GC allocs:57504
# malloc() calls:    9600
# realloc() calls:   1100
# GC pauses:         78

# commit 24b9b1128a41d1c1b7cfd504f147a8f9713c5ce1:
# 27.199742 seconds (30.54 M allocations: 1.727 GiB, 1.81% gc time)
# elapsed time (ns): 27199742000
# gc time (ns):      492167000
# bytes allocated:   1854066288
# pool allocs:       30471464
# non-pool GC allocs:58100
# malloc() calls:    9600
# realloc() calls:   400
# GC pauses:         80

## Meltcurve timings 12 July 2019

# mc3 = test_functions["melting curve single channel"]()
# # open("/tmp/mc3-test.json","w") do f
# #     JSON.print(f, mc3[2]["melt_curve_analysis"])
# # end
# mc3_local = mc3[2]["melt_curve_analysis"]
# mc3_saved = JSON.parsefile("/tmp/mc3-test.json")
# mc3_local == mc3_saved ## should be true
# using Memento
# setlevel!(QpcrAnalysis.logger, "warn")
# # using BenchmarkTools
# # @benchmark test_functions["melting curve single channel"]()
# @timev for i in 1:100; test_functions["melting curve single channel"](); end;

## commit "fixed error in selection of smoothing parameter for melting curve analysis"
#  11.660345 seconds (53.87 M allocations: 8.660 GiB, 11.45% gc time)
# elapsed time (ns): 11660345000
# gc time (ns):      1335187000
# bytes allocated:   9298158176
# pool allocs:       53788340
# non-pool GC allocs:67900
# malloc() calls:    8700
# realloc() calls:   1200
# GC pauses:         406
# full collections:  1

## commit dce3b4a1265df1ca5582283b374b118f0c3d3195
# BenchmarkTools.Trial:
#   memory estimate:  88.80 MiB
#   allocs estimate:  538375
#   --------------
#   minimum time:     122.706 ms (10.05% GC)
#   median time:      126.369 ms (11.01% GC)
#   mean time:        127.271 ms (11.23% GC)
#   maximum time:     149.936 ms (8.78% GC)
#   --------------
#   samples:          40
#   evals/sample:     1

## commit a01b3b50a6e753f3ce92d4a33c1cc9adbb6256b7:
# 12.023655 seconds (52.23 M allocations: 8.743 GiB, 11.55% gc time)
# elapsed time (ns): 12023655000
# gc time (ns):      1388750000
# bytes allocated:   9387864768
# pool allocs:       52151308
# non-pool GC allocs:71200
# malloc() calls:    8700
# realloc() calls:   1200
# GC pauses:         410
# full collections:  1

# commit 73883fd51b000e8c01da05814b0e06b89b153653:
# 13.969750 seconds (52.23 M allocations: 8.743 GiB, 11.36% gc time)
# elapsed time (ns): 13969750000
# gc time (ns):      1586836000
# bytes allocated:   9387864768
# pool allocs:       52151308
# non-pool GC allocs:71200
# malloc() calls:    8700
# realloc() calls:   1200
# GC pauses:         410
# full collections:  1

# commit 24b9b1128a41d1c1b7cfd504f147a8f9713c5ce1:
# 15.593274 seconds (169.50 M allocations: 16.070 GiB, 15.57% gc time)
# elapsed time (ns): 15593274000
# gc time (ns):      2428248000
# bytes allocated:   17254517088
# pool allocs:       169416202
# non-pool GC allocs:64902
# malloc() calls:    13700
# realloc() calls:   600
# GC pauses:         745
# full collections:  3

# mc4 = test_functions["melting curve dual channel"]()
# # open("/tmp/mc4-test.json","w") do f
# #     JSON.print(f, mc4[2]["melt_curve_analysis"])
# # end
# mc4_local = mc4[2]["melt_curve_analysis"]
# mc4_saved = JSON.parsefile("/tmp/mc4-test.json")
# mc4_local == mc4_saved ## should be true
# using Memento
# setlevel!(QpcrAnalysis.logger, "warn")
# # using BenchmarkTools
# # @benchmark test_functions["meltcurve dual channel"]()
# @timev for i in 1:100; test_functions["melting curve dual channel"](); end;

## commit "fixed error in selection of smoothing parameter for melting curve analysis"
#  22.297143 seconds (107.45 M allocations: 15.962 GiB, 10.65% gc time)
# elapsed time (ns): 22297143000
# gc time (ns):      2373774000
# bytes allocated:   17138929472
# pool allocs:       107307234
# non-pool GC allocs:125503
# malloc() calls:    17100
# realloc() calls:   1400
# GC pauses:         747
# full collections:  2

## commit c964da0fe71de22009e5ea3881d1b0629b179e9c
# BenchmarkTools.Trial:
#   memory estimate:  163.67 MiB
#   allocs estimate:  1073208
#   --------------
#   minimum time:     233.911 ms (9.25% GC)
#   median time:      245.089 ms (11.52% GC)
#   mean time:        244.477 ms (11.22% GC)
#   maximum time:     252.217 ms (7.67% GC)
#   --------------
#   samples:          21
#   evals/sample:     1

## commit a01b3b50a6e753f3ce92d4a33c1cc9adbb6256b7:
# BenchmarkTools.Trial:
#   memory estimate:  163.67 MiB
#   allocs estimate:  1073239
#   --------------
#   minimum time:     241.945 ms (9.06% GC)
#   median time:      252.435 ms (11.32% GC)
#   mean time:        253.406 ms (11.08% GC)
#   maximum time:     287.109 ms (10.74% GC)
#   --------------
#   samples:          20
#   evals/sample:     1

# commit 30d20ce41a26d7ad304cf7db57a76e8974f6d055:
# 25.629426 seconds (114.48 M allocations: 16.327 GiB, 11.84% gc time)
# elapsed time (ns): 25629426000
# gc time (ns):      3035671000
# bytes allocated:   17531359168
# pool allocs:       114343408
# non-pool GC allocs:118300
# malloc() calls:    20900
# realloc() calls:   1400
# GC pauses:         764
# full collections:  6

# commit 24b9b1128a41d1c1b7cfd504f147a8f9713c5ce1:
# 30.701926 seconds (330.17 M allocations: 29.970 GiB, 14.68% gc time)
# elapsed time (ns): 30701926000
# gc time (ns):      4506477000
# bytes allocated:   32179894688
# pool allocs:       330025002
# non-pool GC allocs:121402
# malloc() calls:    27300
# realloc() calls:   700
# GC pauses:         1360
# full collections:  5

# mc11 = test_functions["thermal consistency dual channel"]()
# # open("/tmp/mc11-master.json","w") do f
# #     JSON.print(f, mc11[2])
# # end
# master11 = JSON.parsefile("/tmp/mc11-master.json")
# mc11[2] == master11 ## should be true
# @timev for i in 1:100; test_functions["thermal consistency dual channel"](); end;

# standard curve results
# response_body:
# {"targets":[{"target_id":1,"slope":9.119527,"offset":35.45153,"efficiency":-0.223136,"r2":0.075808},{"target_id":2,"slope":-2.499102,"offset":12.596607,"efficiency":1.512718,"r2":0.298742}],"groups":[],"valid":true}
# (true, DataStructures.OrderedDict{String,Any}("targets"=>Any[DataStructures.OrderedDict{String,Any}("target_id"=>1,"slope"=>9.11953,"offset"=>35.4515,"efficiency"=>-0.223136,"r2"=>0.075808), DataStructures.OrderedDict{String,Any}("target_id"=>2,"slope"=>-2.4991,"offset"=>12.5966,"efficiency"=>1.51272,"r2"=>0.298742)],"groups"=>Any[],"valid"=>true))

#===============================================================================
    interactive code (BBB) >>
===============================================================================#

## timing tests on BBB
if (RUN_THIS_CODE_INTERACTIVELY_NOT_ON_INCLUDE & BBB)
    ENV["JULIA_ENV"]="production"
    cd("/root/chaipcr/bioinformatics/QpcrAnalysis")
    push!(LOAD_PATH, pwd())
    using QpcrAnalysis
    include("$(QpcrAnalysis.LOAD_FROM_DIR)/../test/test_functions.jl") ## this file
    test_functions = generate_tests()
    check = test_dispatch(test_functions)
    if all(values(check))
        println("All test functions passed check")
        ## time functions second time around (after compilation)
        timing = time_dispatch(test_functions)
    else
        println("Test functions failed check:")
        println(check)
    end
end

#===============================================================================
    test functions >>
===============================================================================#

function do_curl(
    action          ::AbstractString,
    request_body    ::AbstractString;
    verify          ::Bool = false,
    verbose         ::Bool = true,
)
    ## create text buffer in the form of a closure
    local buffer_contents::String
    local buffer_readptr::Csize_t
    local buffer_length::Csize_t
    #
    function reset_buffer()
        buffer_contents = ""
        buffer_readptr = 1
        buffer_length = 0
    end
    #
    function buffer(
        action      ::Symbol = :get,
        chunk_size  ::Csize_t = Csize_t(0),
        chunk       ::String = ""
    )
        if action == :reset
            reset_buffer()
            return nothing
        elseif action == :length
            return buffer_length
        elseif action == :append
            buffer_contents = buffer_contents * chunk
            buffer_length += chunk_size
            return nothing
        elseif action == :read
            const read_length = min(chunk_size, buffer_length - buffer_readptr + 1)
            const read_range = range(buffer_readptr, read_length)
            buffer_readptr += read_length
            if buffer_readptr > buffer_length
                const return_chunk = buffer_contents[read_range]
                reset_buffer()
                return return_chunk
            else
                return buffer_contents[read_range]
            end
        else ## action == :get
            return buffer_contents
        end
    end
    #
    ## init a curl handle
    const curl = curl_easy_init()
    #
    ## First set the URL that is about to receive our POST
    const exp_id = verify ? "0" : "1"
    const url = "http://127.0.0.1:8081/experiments/" * exp_id * "/" * action
    curl_easy_setopt(curl, CURLOPT_URL, url)
    #
    ## Now specify we want to POST data
    curl_easy_setopt(curl, CURLOPT_POST, 1)
    #
    ## Now specify GET method
    curl_easy_setopt(curl, CURLOPT_CUSTOMREQUEST, "GET")
    #
    ## create custom header to read request
    req_headers = Ptr{Void}(0)
    req_headers = curl_slist_append(req_headers, "Content-Type: application/json")
    curl_easy_setopt(curl, CURLOPT_HTTPHEADER, req_headers)
    #
    ## copy request data into buffer
    const req_text = Array{UInt8}(request_body)
    QpcrAnalysis.print_v(println, verbose, "request_body:", request_body)
    #
    ## callback function to read request
    ## after v0.6: use Ptr{Cvoid} instead of Ptr{Void}
    function curl_read_req(dest::Ptr{Void}, s::Csize_t, n::Csize_t, buf_thunk::Ptr{Void})
        const chunk_size::Csize_t = s * n
        buf = unsafe_pointer_to_objref(buf_thunk)::Function
        eof = (s == 0 || n == 0 || chunk_size < 1)
        if !eof
            const buf_size = buf(:length)
            eof |= iszero(buf_size)
        end
        if eof
            buf(:reset) ## buffer must be emptied before writing the response body into it
            return Csize_t(0)
        end
        ## else
        ## read a chunk of text from buffer
        const source = Array{UInt8}(buf(:read, chunk_size))
        const source_length::Csize_t = length(source)
        ## copy chunk to dest
        ccall(:memcpy, Ptr{Void}, (Ptr{Void}, Ptr{Void}, Csize_t), dest, source, source_length)
        return source_length
    end
    #
    ## use C compatible callback function
    ## after v0.6: use @cfunction instead of cfunction
    const c_curl_read_req =
        cfunction(curl_read_req, Csize_t, (Ptr{Void}, Csize_t, Csize_t, Ptr{Void}))
    curl_easy_setopt(curl, CURLOPT_READFUNCTION, c_curl_read_req)
    #
    ## load request body into buffer
    buffer(:reset)
    buffer(:append, request_body |> length |> Csize_t, request_body)
    #
    ## buffer thunk to our read function
    curl_easy_setopt(curl, CURLOPT_READDATA, pointer_from_objref(buffer))
    #
    ## verbose debug output option
    curl_easy_setopt(curl, CURLOPT_VERBOSE, verbose ? 1 : 0)
    #
    ## set request length
    curl_easy_setopt(curl, CURLOPT_POSTFIELDSIZE, Clong(length(request_body)))
    #
    ## callback function to receive response data
    ## after v0.6: use Ptr{Cvoid} instead of Ptr{Void}
    function curl_write_cb(source::Ptr{Void}, s::Csize_t, n::Csize_t, buf_thunk::Ptr{Void})
        const chunk_size::Csize_t = s * n
        dest = Array{UInt8}(chunk_size)
        ccall(:memcpy, Ptr{Void}, (Ptr{Void}, Ptr{Void}, UInt64), dest, source, chunk_size)
        const chunk = dest |> pointer |> Cstring |> unsafe_string
        # QpcrAnalysis.print_v(println, verbose, "received: " * chunk)
        buf = unsafe_pointer_to_objref(buf_thunk)::Function
        buf(:append, chunk_size, chunk)
        return chunk_size
    end
    #
    ## use C compatible callback function
    ## after v0.6: use @cfunction instead of cfunction
    c_curl_write_cb =
        cfunction(curl_write_cb, Csize_t, (Ptr{Void}, Csize_t, Csize_t, Ptr{Void}))
    curl_easy_setopt(curl, CURLOPT_WRITEFUNCTION, c_curl_write_cb)
    #
    ## thunk to pass to write function
    curl_easy_setopt(curl, CURLOPT_WRITEDATA, pointer_from_objref(buffer))
    #
    ## execute the query
    const res = curl_easy_perform(curl)
    QpcrAnalysis.print_v(println, verbose, "curl response: $res")
    #
    ## get response body
    response_text = buffer()
    QpcrAnalysis.print_v(println, verbose, "response body: " * response_text)
    #
    ## retrieve HTTP code
    if verbose
        const http_code = Array{Clong}(1)
        curl_easy_getinfo(curl, CURLINFO_RESPONSE_CODE, http_code)
        QpcrAnalysis.print_v(println, verbose, "http code: ", http_code[1])
    end
    #
    ## retrieve elapsed time
    if verbose
        const elapsed = Array{Cdouble}(1)
        curl_easy_getinfo(curl, CURLINFO_TOTAL_TIME, elapsed)
        QpcrAnalysis.print_v(println, verbose, "elapsed time: ", elapsed[1])
    end
    #
    ## free headers
    curl_slist_free_all(req_headers)
    #
    ## release handle
    curl_easy_cleanup(curl)
    #
    return(true, response_text)
end

#==============================================================================#

function generate_tests(;
    test_API    ::Bool = true,
    verbose     ::Bool = false
)
    test_functions = OrderedDict()
    strip = [" single"," dual"," channel"]
    for i in 1:size(TEST_DATA)[1]
        for channel_num in [:single_channel, :dual_channel]
            datafile = TEST_DATA[i, channel_num]
            if (datafile != "")
                action_key = TEST_DATA[i, :action]
                # action = Val{QpcrAnalysis.ACT[action_key]}()
                request = JSON.parsefile(
                    "$(QpcrAnalysis.LOAD_FROM_DIR)/../test/data/$datafile.json",
                    dicttype=OrderedDict)
                body = String(JSON.json(request))

                function test_function()
                    QpcrAnalysis.print_v(println, verbose, "Testing $testname")
                    @static BBB || FactCheck.clear_results()
                    if test_API
                        (ok, response_body) = QpcrAnalysis.do_curl(
                            action_key, body; verbose=verbose, verify=false)
                    else
                        (ok, response_body) = QpcrAnalysis.dispatch(
                            action_key, body; verify=false)
                    end
                    response_parsed = JSON.parse(response_body, dicttype=OrderedDict)
                    if (ok && haskey(response_parsed, "valid") && response_parsed["valid"])
                        QpcrAnalysis.print_v(println, verbose, "Passed $testname\n")
                    else
                        QpcrAnalysis.print_v(println, verbose, "Failed $testname\n")
                    end
                    return (ok, response_parsed)
                end

                testname = replace(TEST_DATA[i, :action], r"_"=>" ")
                for str in strip
                    testname = replace(testname, str=>"")
                end
                testname = replace("$testname "*string(channel_num), r"_"=>" ")
                test_functions[testname] = test_function
            end ## if datafile
        end ## single/dual channel (channel_num)
    end ## next action (i)
    return test_functions
end

#==============================================================================#

## generate script to call test functions
## to precompile julia routines for BBB
function generate_test_script(outfile ::String)
    open(outfile, "w") do f
        write(f, """
        println("Starting precompile template !!!")
        push!(LOAD_PATH, "/root/chaipcr/bioinformatics/QpcrAnalysis/")

        println("Using time: ")
        @time using QpcrAnalysis
        println("Done Using!")
        for \$iteration in ["First", "Second"]
            println("\\nAbout to test dispatch. \$iteration time dispatch time:")
        """)
        write_dispatch_calls(f)
        write(f, """
            println("\\nDone dispatch time test")
        end # next iteration
        println("\\nDone with test functions!")
        """)
    end ## close file
end ## generate_test_script()

#==============================================================================#

## write dispatch calls for generate_test_script()
function write_dispatch_calls(f)
    strip = [" single", " dual", " channel"]
    for i in 1:size(TEST_DATA)[1]
        for channel_num in [:single_channel, :dual_channel]
            datafile = TEST_DATA[i, channel_num]
            if (datafile != "")
                action = TEST_DATA[i, :action]
                request = JSON.parsefile("$(QpcrAnalysis.LOAD_FROM_DIR)/../test/data/$datafile.json",
                    dicttype=OrderedDict)
                body = String(JSON.json(request))
                testname = replace(TEST_DATA[i,:action], r"_"=>" ")
                for str in strip
                    testname = replace(testname, str=>"")
                end
                testname = replace("$testname "*string(channel_num), r"_"=>" ")
                write(f, """
                    println("\\nTesting \$testname")
                    action=\"\"\"$action\"\"\"
                    body=\"\"\"$body\"\"\"
                    @time (ok, response_body) =
                        QpcrAnalysis.dispatch(action, body; verify=false)
                    println("OK? \$ok")
                """)
            end ## if datafile
        end ## single/dual channel (channel_num)
    end ## next action (i)
end ## write_dispatch_calls()

#==============================================================================#

## run test functions
## returns true for every test that runs without errors
function test_dispatch(test_functions ::Associative)
    OrderedDict(map(keys(test_functions)) do testname
        println("\nMaking dispatch call: $testname")
        result = test_functions[testname]()
        testname => result[1] && result[2]["valid"]
    end)
end

## time performance
function time_dispatch(test_functions ::Associative)
    OrderedDict(map(keys(test_functions)) do testname
        println("\nMaking dispatch call: $testname")
        @timev result = test_functions[testname]()
        testname => result[1] && result[2]["valid"]
    end)
end

#===============================================================================
    results >>
===============================================================================#

## Meltcurve timings
#
# meltcrv commit  932b24a9be5bb148074830b0fd812618234ccfc1 (don't round mc_denser)
#  10.571982 seconds (61.91 M allocations: 8.668 GiB, 11.33% gc time)
# elapsed time (ns): 10571982383
# gc time (ns):      1197314629
# bytes allocated:   9307300800
# pool allocs:       61839800
# non-pool GC allocs:58000
# malloc() calls:    10900
# realloc() calls:   600
# GC pauses:         405
# full collections:  2

# meltcrv commit  932b24a9be5bb148074830b0fd812618234ccfc1 (remove args from nested funcs)
#  14.142563 seconds (63.41 M allocations: 8.969 GiB, 7.87% gc time)
# elapsed time (ns): 14142562549
# gc time (ns):      1113163232
# bytes allocated:   9630603200
# pool allocs:       63321800
# non-pool GC allocs:79700
# malloc() calls:    12500
# realloc() calls:   600
# GC pauses:         419
# full collections:  2

# meltcrv commit f12f5bda9485e307481be0012fc9ec4555aed0a6 (slowest)
# 18.955378 seconds (49.29 M allocations: 8.766 GiB, 7.62% gc time)
# elapsed time (ns): 18955377866
# gc time (ns):      1443815331
# bytes allocated:   9412268800
# pool allocs:       49198400
# non-pool GC allocs:80300
# malloc() calls:    12500
# realloc() calls:   600
# GC pauses:         410
# full collections:  3

# master commit c39573826114c84d3a516d3ec447c83765871368
# 9.707145 seconds (66.11 M allocations: 8.753 GiB, 12.19% gc time)
# elapsed time (ns): 9707144763
# gc time (ns):      1182875365
# bytes allocated:   9398664064
# pool allocs:       66009104
# non-pool GC allocs:94900
# malloc() calls:    9600
# realloc() calls:   600
# GC pauses:         410
# full collections:  3

#==============================================================================#

# BBB results 2019-01-07
# run by Tom Price as root@10.0.100.231

# Making dispatch call: amplification single channel
#   6.335094 seconds (315.82 k allocations: 13.261 MiB, 3.85% gc time)
# elapsed time (ns): 6335093972
# gc time (ns):      243881118
# bytes allocated:   13905496
# pool allocs:       314368
# non-pool GC allocs:412
# malloc() calls:    550
# realloc() calls:   493
# GC pauses:         2

# Making dispatch call: amplification dual channel
#  13.940308 seconds (690.34 k allocations: 29.719 MiB, 3.35% gc time)
# elapsed time (ns): 13940307816
# gc time (ns):      467057318
# bytes allocated:   31162632
# pool allocs:       687248
# non-pool GC allocs:997
# malloc() calls:    1114
# realloc() calls:   984
# GC pauses:         5

# Making dispatch call: meltcurve single channel
#   4.068145 seconds (700.27 k allocations: 88.775 MiB, 11.46% gc time)
# elapsed time (ns): 4068144765
# gc time (ns):      466131773
# bytes allocated:   93087712
# pool allocs:       699164
# non-pool GC allocs:644
# malloc() calls:    355
# realloc() calls:   105
# GC pauses:         14

# Making dispatch call: meltcurve dual channel
#  48.444924 seconds (8.27 M allocations: 393.947 MiB, 10.21% gc time)
# elapsed time (ns): 48444923803
# gc time (ns):      4946285937
# bytes allocated:   413083800
# pool allocs:       8271275
# non-pool GC allocs:1419
# malloc() calls:    702
# realloc() calls:   202
# GC pauses:         63
# full collections:  1

# Making dispatch call: standard curve single channel
#   0.084199 seconds (1.90 k allocations: 86.211 KiB)
# elapsed time (ns): 84198636
# bytes allocated:   88280
# pool allocs:       1900
# non-pool GC allocs:2

# Making dispatch call: optical cal single channel
#   0.008529 seconds (752 allocations: 28.359 KiB)
# elapsed time (ns): 8529347
# bytes allocated:   29040
# pool allocs:       752

# Making dispatch call: optical cal dual channel
#   0.023195 seconds (2.73 k allocations: 112.789 KiB)
# elapsed time (ns): 23195099
# bytes allocated:   115496
# pool allocs:       2734

# Making dispatch call: thermal performance diagnostic single channel
#   0.024872 seconds (6.20 k allocations: 206.281 KiB)
# elapsed time (ns): 24872069
# bytes allocated:   211232
# pool allocs:       6197
# non-pool GC allocs:3

# Making dispatch call: thermal performance diagnostic dual channel
#   0.020857 seconds (4.77 k allocations: 159.742 KiB)
# elapsed time (ns): 20857365
# bytes allocated:   163576
# pool allocs:       4767
# non-pool GC allocs:3

# Making dispatch call: thermal consistency single channel
#   1.408891 seconds (196.18 k allocations: 27.054 MiB, 16.47% gc time)
# elapsed time (ns): 1408891471
# gc time (ns):      232104860
# bytes allocated:   28368096
# pool allocs:       195557
# non-pool GC allocs:451
# malloc() calls:    90
# realloc() calls:   83
# GC pauses:         4

# Making dispatch call: thermal consistency dual channel
#  19.650871 seconds (3.35 M allocations: 168.469 MiB, 8.63% gc time)
# elapsed time (ns): 19650871467
# gc time (ns):      1696740250
# bytes allocated:   176652784
# pool allocs:       3345526
# non-pool GC allocs:986
# malloc() calls:    178
# realloc() calls:   163
# GC pauses:         27

# Making dispatch call: optical test single channel
#   0.002865 seconds (309 allocations: 12.695 KiB)
# elapsed time (ns): 2865254
# bytes allocated:   13000
# pool allocs:       309

# Making dispatch call: optical test dual channel
#   0.024774 seconds (4.14 k allocations: 175.016 KiB)
# elapsed time (ns): 24774074
# bytes allocated:   179216
# pool allocs:       4135
# non-pool GC allocs:1
# realloc() calls:   1

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
export AbstractSequence, AbstractMultiEchoSequence, echoAmplitudes, numEchoes, echoTimes, sequence, sequenceInfo

abstract type AbstractSequence end

# define interface for a sequence
@mustimplement echoAmplitudes(seq::AbstractSequence,R1::AbstractFloat,R2::AbstractFloat)
@mustimplement numContrasts(seq::AbstractSequence)
@mustimplement string(seq::AbstractSequence)

include("BlochSimulation.jl")
include("MultiEcho.jl")

function sequence(seqName::AbstractString, numProfiles, numSamplingPerProfile; kargs...)
  if seqName == "ME"
    return MESequence(; kargs...)
  else
    error("The sequence $seqName is not yet supported!")
  end

  return MESequence(; kargs...)
end

function sequenceInfo(seq::AbstractSequence)
  param = Dict{Symbol,Any}()
  param[:name] = string(seq)
  for field in fieldnames(typeof(seq))
    param[field] = getfield(seq,field)
  end
  return param
end


text:
using DataStructures

function tp_sort(G)
	n = length(G)

	# Create Reversed Graph G′
	G′ = [Vector{Int}() for _ in 1:n]
	for i = 1:n, e ∈ G[i]
		push!(G′[e], i)
	end

	num_ins = [length(g) for g ∈ G]

	que = Queue{Int}()
	for i = 1:n
		num_ins[i] == 0 && enqueue!(que, i)
	end

	ans = MutableLinkedList{Int}()
	while !isempty(que)
		now = dequeue!(que)
		pushfirst!(ans, now)
		for e ∈ G′[now]
			num_ins[e] -= 1
			num_ins[e] == 0 && enqueue!(que, e)
		end
	end

	return collect(ans)
end

function main()
    N, M = parse.(Int, split(readline()))

    G = [Vector{Int}() for _ in 1:N]
	for _ in 1:M
		src, dst = parse.(Int, split(readline()))
		push!(G[src], dst)
	end

    dp = zeros(Int, N) #dp[i] := the length of the longest path if i is the last node.
    tps = tp_sort(G)
    for v ∈ tps, v′ ∈ G[v]
        dp[v′] = max(dp[v′], dp[v] + 1)
    end

    println(maximum(dp))
end

main()


text:
# # [Ecosystem Models](@id ecosystem_example)
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/examples/Ecosystem.ipynb)

using AlgebraicDynamics
using AlgebraicDynamics.DWDDynam
using AlgebraicDynamics.UWDDynam

using Catlab.CategoricalAlgebra
using Catlab.WiringDiagrams
using Catlab.Graphics
using Catlab.Programs

using OrdinaryDiffEq
using Plots, Plots.PlotMeasures

# ## Land Ecosystem

# ### Rabbits and foxes

# A standard Lotka Volterra predator-prey model is the composition of three primitive resource sharers:

# 1. a model of rabbit growth: this resource sharer has dynamics $\dot r(t) = \alpha r(t)$ and one port which exposes the rabbit population.
# 2. a model of rabbit/fox predation: this resource sharer has dynamics $$\dot r(t) = -\beta r(t) f(t), \dot f(t) = \gamma r(t)f(t)$$ and two ports which expose the rabbit and fox populations respectively.
# 3. a model of fox population decline: this resource sharer has dynamics $\dot f(t) = -\delta f(t)$ and one port which exposes the fox population.

# However, there are not two independent rabbit populations -- one that grows and one that gets eaten by foxes. Likewise, there are not two independent fox populations -- one that declines and one that feasts on rabbits. To capture these interactions between the trio of resource sharers, we compose them by identifying the exposed rabbit populations and identifying the exposed fox populations.
# The syntax for this undirected composition is defined by an undirected wiring diagram.

## Define the primitive systems
dotr(u,p,t) = p[1]*u
dotrf(u,p,t) = [-p[2]*u[1]*u[2], p[3]*u[1]*u[2]]
dotf(u,p,t) = -p[4]*u

rabbit_growth = ContinuousResourceSharer{Float64}(1, dotr)
rabbitfox_predation = ContinuousResourceSharer{Float64}(2, dotrf)
fox_decline = ContinuousResourceSharer{Float64}(1, dotf)

## Define the composition pattern
rabbitfox_pattern = @relation (rabbits, foxes) begin
    rabbit_growth(rabbits)
    rabbitfox_predation(rabbits,foxes)
    fox_decline(foxes)
end

## Compose
rabbitfox_system = oapply(rabbitfox_pattern, [rabbit_growth, rabbitfox_predation, fox_decline])

# Previously, when we derived the Lotka-Volterra model via [undirected composition](https://algebraicjulia.github.io/AlgebraicDynamics.jl/dev/examples/Lotka-Volterra/#Undirected-composition), we by-hand defined the undirected wiring diagram that implements the composition pattern. In contrast, here we implement the same composition pattern as before but this time using the [`@relation` macro](https://algebraicjulia.github.io/Catlab.jl/stable/apis/programs/#Catlab.Programs.RelationalPrograms.@relation-Tuple). This strategy simplifies the definition and explicitly names the boxes and variables. We  visualize the composition pattern below.

to_graphviz(rabbitfox_pattern, box_labels = :name, junction_labels = :variable, edge_attrs=Dict(:len => ".75"))

# We can now construct an `ODEProblem` from the resource sharer `rabbitfox_system` and plot the solution.

α, β, γ, δ = 0.3, 0.015, 0.015, 0.7
params = [α, β, γ, δ]

u0 = [10.0, 100.0]
tspan = (0.0, 100.0)

prob = ODEProblem(rabbitfox_system, u0, tspan, params)
sol = solve(prob, Tsit5())

plot(sol, lw=2, title = "Lotka-Volterra Predator-Prey Model", bottom_margin=10mm, left_margin=10mm, label=["rabbits" "foxes"])
xlabel!("Time")
ylabel!("Population size")

# ### Rabbits, foxes, and hawks
# Suppose we now have a three species ecosystem containing rabbits, foxes, and hawks. Foxes and hawks both prey upon rabbits but do not interact with each other. This ecosystem consists of five primitive systems which share variables.
# 1. rabbit growth:  $\dot r(t) = \alpha r(t)$
# 2. rabbit/fox predation:  $\dot r(t) = -\beta r(t) f(t), \dot f(t) = \delta r(t)f(t)$
# 3. fox decline:  $\dot f(t) = -\gamma f(t)$
# 4. rabbit/hawk predation: $\dot r(t) = -\beta' r(t)h(t), \dot h(t) = \delta' r(t)h(t)$
# 5. hawk decline:  $\dot h(t) = -\gamma' h(t)$

# This means the desired composition pattern has five boxes and many ports and wires to keep track of. Instead of implementing this composition pattern by hand, we construct it as a pushout.

## Define the composition pattern for rabbit growth
rabbit_pattern = @relation (rabbits,) -> rabbit_growth(rabbits)

## Define the composition pattern for the rabbit/hawk Lotka Volterra model
rabbithawk_pattern = @relation (rabbits, hawks) begin
    rabbit_growth(rabbits)
    rabbit_hawk_predation(rabbits,hawks)
    hawk_decline(hawks)
end

## Define transformations between the composition patterns
rabbitfox_transform  = ACSetTransformation((Box=[1], Junction=[1], Port=[1], OuterPort=[1]), rabbit_pattern, rabbitfox_pattern)
rabbithawk_transform = ACSetTransformation((Box=[1], Junction=[1], Port=[1], OuterPort=[1]), rabbit_pattern, rabbithawk_pattern)

## Take the pushout to define the composition pattern for the rabbit, fox, hawk system
rabbitfoxhawk_pattern = ob(pushout(rabbitfox_transform, rabbithawk_transform))

## Visualize the compsition pattern
to_graphviz(rabbitfoxhawk_pattern, box_labels = :name, junction_labels = :variable, edge_attrs=Dict(:len => ".9"))

#-
## Define the additional primitive systems
dotrh(x, p, t) = [-p[5]*x[1]*x[2], p[6]*x[1]*x[2]]
doth(x, p, t)  = -p[7]*x

rabbithawk_predation = ContinuousResourceSharer{Float64}(2, dotrh)
hawk_decline         = ContinuousResourceSharer{Float64}(1, doth)

## Compose
land_system = oapply(rabbitfoxhawk_pattern,
                        [rabbit_growth, rabbitfox_predation, fox_decline,
                         rabbithawk_predation, hawk_decline])

## Solve and plot
β′, γ′, δ′ = .01, .01, .5
params = vcat(params, [β′, γ′, δ′])

u0 = [10.0, 100.0, 50.0]
tspan = (0.0, 100.0)

prob = ODEProblem(land_system, u0, tspan, params)
sol = solve(prob, Tsit5())

plot(sol, lw=2, title = "Land Ecosystem", bottom_margin=10mm, left_margin=10mm, label=["rabbits" "foxes" "hawks"])
xlabel!("Time")
ylabel!("Population size")

# Unfortunately, the hawks are going extinct in this model. We'll have to give hawks something else to eat!
#-
# ## Ocean Ecosystem

# Consider a ocean ecosystem containing three species —- little fish, big fish, and sharks -— with two predation interactions —- sharks eat big fish and big fish eat little fish.

# This ecosystem can be modeled as the composition of 3 machines:
# 1. Evolution of the little fish population:  this machine has one exogenous variable which represents a population of predators $h(t)$ that hunt little fish. This machine has one output which emits the little fish population. The dynamics of this machine is the driven ODE $$\dot f(t) = \alpha f(t) - \beta f(t)h(t)$$
# 2. Evolution of the big fish population:  this machine has two exogenous variables which represent a population of prey $e(t)$ that are eaten by big fish and a population of predators $h(t)$ which hunt big fish. This machine has one output which emits the big fish population. The dynamics of this machine is the drive ODE $$\dot F(t) = \gamma F(t)e(t) - \delta F(t) - \beta'F(t)h(t)$$
# 3. Evolution of the shark population:  this machine has one exogenous variable which represents a population of prey $e(t)$ that are eaten by sharks. This machine has one output which emits the shark population. The dynamics of this machine is the driven ODE $$\dot s(t) = \gamma's(t)e(t) - \delta's(t)$$

## Define the primitive systems
dotfish(f, x, p, t) = [p[1]*f[1] - p[2]*x[1]*f[1]]
dotFISH(F, x, p, t) = [p[3]*x[1]*F[1] - p[4]*F[1] - p[5]*x[2]*F[1]]
dotsharks(s, x, p, t) = [p[6]*s[1]*x[1]-p[7]*s[1]]

fish   = ContinuousMachine{Float64}(1,1,1, dotfish,   f->f)
FISH   = ContinuousMachine{Float64}(2,1,1, dotFISH,   F->F)
sharks = ContinuousMachine{Float64}(1,1,1, dotsharks, s->s)

# We compose these machines by (1) sending the output of the big fish machine as the input to both the little fish and shark machines and (2) sending the output of the little fish and shark machines as the inputs to the big fish machine.
# The syntax for this directed composition is given by a directed wiring diagram.

## Define the composition pattern
ocean_pattern = WiringDiagram([], [])
fish_box = add_box!(ocean_pattern, Box(:fish, [:pop], [:pop]))
Fish_box = add_box!(ocean_pattern, Box(:Fish, [:pop, :pop], [:pop]))
shark_box = add_box!(ocean_pattern, Box(:shark, [:pop], [:pop]))

add_wires!(ocean_pattern, Pair[
    (fish_box, 1)  => (Fish_box, 1),
    (shark_box, 1) => (Fish_box, 2),
    (Fish_box, 1)  => (fish_box, 1),
    (Fish_box, 1)  => (shark_box, 1)
])

## Visualize the composition pattern
to_graphviz(ocean_pattern, orientation=TopToBottom)

#-

## Compose
ocean_system = oapply(ocean_pattern, [fish, FISH, sharks])

## Solve and plot
α, β, γ, δ, β′, γ′, δ′ = 0.35, 0.015, 0.015, 0.7, 0.017, 0.017, 0.35
params = [α, β, γ, δ, β′, γ′, δ′]

u0 = [100.0, 10, 2.0]
tspan = (0.0, 100.0)

prob = ODEProblem(ocean_system, u0, tspan, params)
sol = solve(prob, FRK65(0))

plot(sol, lw=2, title = "Ocean Ecosystem", bottom_margin=10mm, left_margin=10mm, label=["little fish" "big fish" "sharks"])
xlabel!("Time")
ylabel!("Population size")

# ## Total ecosystem
# ### Another layer of composition

# We will introduce a final predation interaction -- hawks eat little fish --  which will combine the land and ocean ecosystems.

# There will be 16 parameters in to the total ecosystem.
# - parameters 1-7 will determine the land ecosystem
# - parameters 8 and 9 will determine the hawk/little fish predation. Parameter 8 gives the rate of hawk growth and parameter 8 gives the rate of little fish decline.
# - parameter 10-16 will determine the ocean ecosystem.

# The composition will be as resource shareres so the first thing we will do is use the dynamics of the machine `ocean_system` to define the dynamics of a resource sharer. We will also define a resource sharer that models hawk/little fish predation.

## Define the additional primitive systems
ocean_system_rs = ContinuousResourceSharer{Float64}(3, (u,p,t)->eval_dynamics(ocean_system, u, [], p[10:16]))

dothf(u,p,t) = [p[8]*u[1]*u[2], -p[9]*u[1]*u[2]]
fishhawk_predation = ContinuousResourceSharer{Float64}(2, dothf)

## Define the composition pattern
eco_pattern = @relation () where (rabbits, foxes, hawks, littlefish, BigFish, sharks)  begin
    turf(rabbits,foxes,hawks)
    air(hawks, littlefish)
    surf(littlefish, BigFish, sharks)
end

## Visualize the composition pattern
to_graphviz(eco_pattern, box_labels = :name, junction_labels = :variable, edge_attrs=Dict(:len => ".75"))

#-
## Compose
ecosystem=oapply(eco_pattern, [land_system, fishhawk_predation, ocean_system_rs])

# We can now plot the evolution of the total ecosystem.

## Solve and plot
u0 = [100.0, 50.0, 20.0, 100, 10, 2.0]
tspan = (0.0, 100.0)

params = [0.3, 0.015, 0.015, 0.7, .01, .01, .5,
          0.001, 0.003,
          0.35, 0.015, 0.015, 0.7, 0.017, 0.017, 0.35]

prob = ODEProblem(ecosystem, u0, tspan, params)
sol = solve(prob, Tsit5())
plot(sol, lw=2, bottom_margin=10mm, left_margin=10mm, label = ["rabbits" "foxes" "hawks" "little fish" "big fish" "sharks"])

# Let's zoom in on a narrower time-window.
tspan = (0.0, 30.0)

prob = ODEProblem(ecosystem, u0, tspan, params)
sol = solve(prob, Tsit5())
plot(sol, lw=2, bottom_margin=10mm, left_margin=10mm, label = ["rabbits" "foxes" "hawks" "little fish" "big fish" "sharks"])

# As a sanity check we can define the rates for the hawk/little fish predation to be 0. This decouples the land and ocean ecosystems. As expected, the plot shows the original evolution of the land ecosystem overlayed with the original evolution of the ocean ecosystem. This shows that they two ecosystems now evolve independently.

tspan = (0.0, 100.0)
params = [0.3, 0.015, 0.015, 0.7, .01, .01, .5,
          0, 0,
          0.35, 0.015, 0.015, 0.7, 0.017, 0.017, 0.35]
prob = ODEProblem(ecosystem, u0, tspan, params)
sol = solve(prob, Tsit5())
plot(sol, lw=2, bottom_margin=10mm, left_margin=10mm, label = ["rabbits" "foxes" "hawks" "little fish" "big fish" "sharks"])


text:
module sphsld_acous_coupl_examples
using FinEtools
using FinEtools.MeshExportModule
using UnicodePlots
import Arpack: eigs

function sphsld_acous_coupl_examples_H8()
    R = 0.5 * phun("m")
    E = 205000 * phun("MPa")
    nu = 0.3
    rho = 7850 * phun("kg/m^3");# mass density
    dummybulk, dummyrho = (1.0, 1.0)
    neigvs=6+18;
    OmegaShift=(2*pi*100.0)^2;
    nperradius = 2
    tolerance = R/10000; # geometrical tolerance

    # Construct the mesh
    r(c) = c[[1, 4, 3, 2, 5, 8, 7, 6]]
    origin = [0.0, 0.0, 0.0]
    fens,fes = H8spheren(R, nperradius); # Mesh
    fens1,fes1  =  mirrormesh(fens, fes, [-1.0, 0.0, 0.0], origin,  renumb = r);
    fens,newfes1,fes2 =  mergemeshes(fens1, fes1, fens, fes, tolerance)
    fes = cat(newfes1,fes2)
    fens1,fes1  =  mirrormesh(fens, fes, [0.0, -1.0, 0.0], origin,  renumb = r);
    fens,newfes1,fes2 =  mergemeshes(fens1, fes1, fens, fes, tolerance)
    fes = cat(newfes1,fes2)
    fens1,fes1  =  mirrormesh(fens, fes, [0.0, 0.0, -1.0], origin,  renumb = r);
    fens,newfes1,fes2 =  mergemeshes(fens1, fes1, fens, fes, tolerance)
    fes = cat(newfes1,fes2)

    # Extract the mesh of the boundary
    bfes  =  meshboundary(fes)

    # Debugging graphics
    # File  =   "Sphere.vtk"
    # vtkexportmesh(File, fes.conn, fens.xyz, MeshExportModule.H8)
    # @async run(`"paraview.exe" $File`)

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field

    numberdofs!(u)

    MR = DeforModelRed3D
    material = MatDeforElastIso(MR, rho, E, nu, 0.0)

    femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3, 2)), material)
    associategeometry!(femm,  geom)
    K  = stiffness(femm, geom, u)
    M = mass(femm, geom, u)

    d,v,nev,nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM)
    d = d .- OmegaShift;
    fs = real(sqrt.(complex(d)))/(2*pi)
    println("Eigenvalues: $fs [Hz]")

    femm = FEMMAcoustSurf(IntegDomain(bfes, TrapezoidalRule(2)), MatAcoustFluid(dummybulk, dummyrho))
    G = acousticcouplingpanels(femm, geom, u);
    @show G

    true

end # sphsld_acous_coupl_examples_H8

function allrun()
    println("#####################################################")
    println("# sphsld_acous_coupl_examples_H8 ")
    sphsld_acous_coupl_examples_H8()
end # function allrun

end # module sphsld_acous_coupl_examples


text:
#
# Implementation of MRG32k3a RNG Streams
#
########################################

using Random

import Base.rand

const m1 = Int64(2^32 - 209)
const m2 = Int64(2^32 - 22853)
const a12 = Int64(1403580)
const a13 = Int64(-810728)
const a21 = Int64(527612)
const a23 = Int64(-1370589)
const norm = Float64(1.0 / (1 + m1))
const two17 = Int64(131072)
const two53 = Int64(9007199254740992)

const DEFAULT_SEED = [ 12345, 12345, 12345, 12345, 12345, 12345 ]

const A1p0 =  [  0     1     0   ;
                 0     0     1   ;
                 a13   a12   0   ]

const A2p0 =  [  0     1     0   ;
                 0     0     1   ;
                 a23   0     a21 ]

const A1p76 =  [   82758667  1871391091  4127413238 ;
                 3672831523    69195019  1871391091 ;
                 3672091415  3528743235    69195019 ]

const A2p76 =  [ 1511326704  3759209742  1610795712 ;
                 4292754251  1511326704  3889917532 ;
                 3859662829  4292754251  3708466080 ]

const A1p127 = [ 2427906178  3580155704   949770784 ;
                  226153695  1230515664  3580155704 ;
                 1988835001   986791581  1230515664 ]

const A2p127 = [ 1464411153   277697599  1610723613 ;
                   32183930  1464411153  1022607788 ;
                 2824425944    32183930  2093834863 ]

const InvA1  = [  184888585           0  1945170933 ;
                          1           0           0 ;
                          0           1           0 ]

const InvA2 =  [          0   360363334  4225571728 ;
                          1           0           0 ;
                          0           1           0 ]

"""
Ensures a given seed is valid for MRG32k3a random number generator.
"""
function checkseed(x::Vector{Int})
    return length(x) == 6     &&
           all(x[1:6] .>= 0)  &&
           all(x[1:3] .< m1)  &&
           all(x[4:6] .< m2)  &&
          ~all(x[1:3] .== 0)   &&
          ~all(x[4:6] .== 0)
end

mutable struct MRG32k3a <: AbstractStreamableRNG
    Cg::Vector{Int64}  # the current state of the RNG
    Bg::Vector{Int64}  # the start point of the current substream
    Ig::Vector{Int64}  # the start point of the current stream

    function MRG32k3a(x::Vector{Int})
        @assert(checkseed(x))
        new(copy(x),copy(x),copy(x))
    end
    function MRG32k3a(x::Vector{Int}, y::Vector{Int}, z::Vector{Int})
        @assert(checkseed(x))
        @assert(checkseed(y))
        @assert(checkseed(z))
        return new(copy(x),copy(y),copy(z))
    end
end

"""
Produces a raw random number with 32 bits of precision.
"""

function rand(rng::MRG32k3a)

    p1::Int64 = (a12 * rng.Cg[2] + a13 * rng.Cg[1]) % m1
    p1 += p1 < 0 ? m1 : 0

    rng.Cg[1] = rng.Cg[2]
    rng.Cg[2] = rng.Cg[3]
    rng.Cg[3] = p1

    p2::Int64 = (a21 * rng.Cg[6] + a23 * rng.Cg[4]) % m2
    p2 += p2 < 0 ? m2 : 0

    rng.Cg[4] = rng.Cg[5]
    rng.Cg[5] = rng.Cg[6]
    rng.Cg[6] = p2

    u::Float64 = p1 > p2 ? (p1 - p2) * norm : (p1 + m1 - p2) * norm
end

"""
Seeds a given random number generator with seed x.
"""

function srand(rng::MRG32k3a,seed::Vector{Int64})
    @assert(checkseed(seed))
    for i = 1:6
        rng.Cg[i] = rng.Bg[i] = rng.Ig[i] = seed[i]
    end
end
"""
Resets a given random number generator to the beginning of the current stream.
"""
function reset_stream!(rng::MRG32k3a)
    for i = 1:6
        rng.Cg[i] = rng.Bg[i] = rng.Ig[i]
    end
end

"""
Resets a given random number generator to the beginning of the current substream.
"""
function reset_substream!(rng::MRG32k3a)
    for i = 1:6
        rng.Cg[i] = rng.Bg[i]
    end
end
"""
Takes a random number generator and shifts seed to next substream.
"""
function next_substream!(rng::MRG32k3a)
    rng.Bg[1:3] = MatVecModM(A1p76,rng.Bg[1:3],m1)
    rng.Bg[4:6] = MatVecModM(A2p76,rng.Bg[4:6],m2)
    for i = 1:6
        rng.Cg[i] = rng.Bg[i]
    end
end

# Required to extend randn(), randexp(), and other built-in random functions

#rand(rng::MRG32k3a, ::Type{Base.Random.Close1Open2}) = rand(rng, Base.Random.CloseOpen) + 1.0
#rand(rng::MRG32k3a, ::Type{Base.Random.CloseOpen}) = rand(rng::MRG32k3a, Float64)

############################################################
"""
An object that generates independent random number streams.
"""
mutable struct MRG32k3aGen <: AbstractRNGStream
    nextSeed::Vector{Int64}

    function MRG32k3aGen()
        new(copy(DEFAULT_SEED))
    end

    function MRG32k3aGen(x::Vector{Int})
        @assert(checkseed(x))
        new(copy(x))
    end

end

function show(io::IO,rng_gen::MRG32k3aGen)
    print(io,"Seed for next MRG32k3a generator:\n$(rng_gen.nextSeed)")
end

function show(io::IO,rng::MRG32k3a)
    print(io,"Full state of MRG32k3a generator:\nCg = $(rng.Cg)\nBg = $(rng.Bg)\nIg = $(rng.Ig)")
end

get_state(rng_gen::MRG32k3aGen) = copy(rng_gen.nextSeed)

get_state(rng::MRG32k3a) = copy(rng.Cg)

"""
Reseeds the RNG generator object.
"""
function srand(rng_gen::MRG32k3aGen,seed::Vector{Int64})
    @assert(checkseed(seed))
    for i = 1:6
        rng_gen.nextSeed[i] = seed[i]
    end
end

"""
Given an RNG generator object, returns the next RNG stream.
"""
function next_stream(rng_gen::MRG32k3aGen)
    rng = MRG32k3a(copy(rng_gen.nextSeed))

    rng_gen.nextSeed[1:3] = MatVecModM(A1p127,rng_gen.nextSeed[1:3],m1)
    rng_gen.nextSeed[4:6] = MatVecModM(A2p127,rng_gen.nextSeed[4:6],m2)

    return rng
end

"""
Computes (a*s + c) % m, all must be < 2^35. Overflow-safe.
"""
function MultModM(a::Int64, s::Int64, c::Int64, m::Int64)
    val = a * Float64(s) + c
    if abs(val) < two53
        v = Int64(val)
    else
        a1 = a ÷ two17
        a -= a1 * two17
        v  = (a1 * s) % m
        v  = v * two17 + a * s + c
    end
    v %= m
    v += (v < 0) ? m : 0
    return v
end

"""
Computes A*s % m, assuming abs(s[i]) < m. Overflow-safe.
"""
function MatVecModM(A::Array{Int64,2},s::Array{Int64,1}, m::Int64)

    v = [0,0,0]
    for i = 1:3
        for j = 1:3
            v[i] = MultModM(A[i,j], s[j], v[i], m)
        end
    end

    return v
end

"""
Computes matrix A*B % m, assuming abs(s[i]) < m. Overflow-safe.
"""
function MatMatModM(A::Array{Int64,2}, B::Array{Int64,2}, m::Int64)

    C = diagm([0,0,0])
    for i = 1:3
        C[:,i] = MatVecModM(A,B[:,i],m)
    end

    return C
end

"""
Computes the matrix A^(2^e) % m. Overflow-safe.
"""
function MatTwoPowModM(A::Array{Int64,2}, e::Int64, m::Int64)

    B = A
    for i = 1:e
        B = MatMatModM(B, B, m)
    end

    return B
end

"""
Computes the matrix  (A^n % m). Overflow-safe.
"""
function MatPowModM(A::Array{Int64,2}, n::Int64, m::Int64)
    W = A
    B = diagm([1,1,1])

    while n > 0
        if ( n % 2 == 1 )
            B = MatMatModM(W, B, m)
        end
        W = MatMatModM(W, W, m)
        n ÷= 2
    end
    return B
end

"""
Given a random number generator, jumps n steps forward if n > 0
(or -n steps backwards if n < 0), where

if e > 0, let n = 2^e + c;
if e < 0, let n = -2^(-e) + c;
if e = 0, let n = c.
"""
function advance_state!(rng::MRG32k3a, e::Int64, c::Int64)
    if c >= 0
        C1 = MatPowModM(A1p0,c,m1)
        C2 = MatPowModM(A2p0,c,m2)
    else
        C1 = MatPowModM(InvA1,-c,m1)
        C2 = MatPowModM(InvA2,-c,m2)
    end

    if e > 0
        B1 = MatTwoPowModM(A1p0,e,m1)
        B2 = MatTwoPowModM(A2p0,e,m2)
    elseif e < 0
        B1 = MatTwoPowModM(InvA1,-e,m1)
        B2 = MatTwoPowModM(InvA2,-e,m2)
    end

    if ~(e == 0)
        C1 = MatMatModM(B1,C1,m1)
        C2 = MatMatModM(B2,C2,m2)
    end

    rng.Cg[1:3] = MatVecModM(C1,rng.Cg[1:3],m1)
    rng.Cg[4:6] = MatVecModM(C2,rng.Cg[4:6],m2)
end


text:
using MeshTools
using Test

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


text:
function formatmatrix(inputmatrix)

  numrows = size(inputmatrix,1)
  numcols = size(inputmatrix,2)

  trans = inputmatrix'

  firstcol = trans[1,:]
  outputmatrix = map(x -> x/sum(firstcol),firstcol)

  for i=2:numcols
    row = trans[i,:]
    newrow = map(x -> x/sum(row),row)
      outputmatrix = vcat(outputmatrix,newrow)
  end
  outputmatrix'
end

mat = [1 2 3; 4 5 6]

println(formatmatrix(mat))